2

I've been trying to translate a set of chemical equations to MATLAB code, to be able to solve for different chemical species. I have the approximate solution (as it's from a graph) but after entering all the data and checking multiple times I still haven't been able to find what is wrong. I'm wondering what is going wrong and if anyone could please help me out. The source for the graph/equation is the article at this link: The chemistry of co-injected BOE. The graph I want to reproduce later on is figure 2 in the paper, see the image below:

Graph of fraction of species from article

Now the results I get for 10cc, 40cc and 90cc are respectively:
HF 43%, H2F2 48%, F- 3%, HF2- 6% in comparison ~28%, 63%, 2%, 7% (10cc).
HF 35%, H2F2 33%, F- 14%, HF2- 18% in comparison ~24%, 44%, 6%, 26% (40cc).
HF 21%, H2F2 12%, F- 37%, HF2- 30% in comparison ~18%, 23%, 20%, 45% (90cc).

The script is the following:

clc;
clear all;
%Units to be used 
%Volume is in CC also cm^3, 1 litre is 1000 CC, 1 cc = 1 ml
%density is in g/cm^3
%weigth percentages are in fractions of 0 to 1
%Molecular weight is in g/mol
% pts=10; %number of points for linear spacing

%weight percentages of NH4OH and HF
xhf=0.49;
xnh3=0.28;

%H2O
Vh2o=1800;
dh2o=1.00; %0.997 at 25C when rounded 1
mh2o=18.02;

%HF values
Vhf=100;
dhf49=1.15;
dhf=dh2o+(dhf49-dh2o)*xhf/0.49; %@ 25C
Mhf=20.01;
nhf=mols(Vhf,dhf,xhf,Mhf);

%NH4OH (NH3) values
% Vnh3=linspace(0.1*Vhf,1.9*Vhf,pts);
Vnh3=10;
dnh3=0.9; %for ~20-31% @~20-25C
Mnh3=17.03; %The wt% of NH4OH actually refers to the wt% of NH3 dissolved in H2O
nnh3=mols(Vnh3,dnh3,xnh3,Mnh3);

if max(nnh3)>=nhf
    error(['There are more mols NH4OH,',num2str(max(nnh3)),', than mols HF,',num2str(nhf),'.'])
end

%% Calculations for species 
Vt=(Vhf+Vh2o+Vnh3)/1000; %litre
A=nhf/Vt; %mol/l
B=nnh3/Vt; %mol/l
syms HF F H2F2 HF2 NH3 NH4 H OH
eq2= H*F/HF==6.85*10^(-4);
eq3= NH3*H/NH4==6.31*10^(-10);
eq4= H*OH==10^(-14);
eq5= HF2/(HF*F)==3.963;
eq6= H2F2/(HF^2)==2.7;
eq7= H+NH4==OH+F+HF2;
eq8= HF+F+2*H2F2+2*HF2==A;
eq9= NH3+NH4==B;
eqns=[eq2,eq5,eq6,eq8,eq4,eq3,eq9,eq7];
varias=[HF, F, H2F2, HF2, NH3, NH4, H, OH];
assume(HF> 0 & F>= 0 & H2F2>= 0 & HF2>= 0& NH3>= 0 & NH4>= 0 & H>= 0 & OH>= 0)
[HF, F, H2F2, HF2, NH3, NH4, H, OH]=vpasolve(eqns,varias);% [0 max([A,B])])

totalHF=double(HF)+double(F)+double(H2F2)+double(HF2);
HFf=double(HF)/totalHF %fraction of species for HF
H2F2f=double(H2F2)/totalHF %fraction of species for H2F2
Ff=double(F)/totalHF %fraction of species for F-
HF2f=double(HF2)/totalHF %fraction of species for HF2-

an extra function needed is called mols.m

%%%% amount of mol, Vol=volume, d=density, pwt=%weight, M=molecularweight
function mol=mols(Vol, d, pwt, M)
mol=(Vol*d*pwt)/M;
end


The equations being used from the article are in the image below:
Equations from article

(HF)2 is H2F2 in my script

Bob van de Voort
  • 211
  • 1
  • 11
  • Welcome to Stack Overflow! Please show the equations you are trying to translate, it will help people spot any differences. Also, I've improved (or tried to improve) your title, I hope you don't mind. – Cris Luengo May 23 '18 at 18:48
  • Thanks, I'll add the equations as an image since that way I can't make any typos. – Bob van de Voort May 23 '18 at 18:54
  • This has been solved here ;) https://chemistry.stackexchange.com/questions/98306/why-do-my-equilibrium-calculations-on-this-hf-nh4oh-buffer-system-not-match-thos – Bob van de Voort Jun 27 '18 at 23:13

1 Answers1

0

So appears the issue wasn't so much with Matlab, had some help in that area as well.

Final solution and updated Matlab code can be found here: https://chemistry.stackexchange.com/questions/98306/why-do-my-equilibrium-calculations-on-this-hf-nh4oh-buffer-system-not-match-thos

Bob van de Voort
  • 211
  • 1
  • 11