翻译文章中的化学方程式,结果不同(Matlab)

Translating chemical equations from article, results differ (Matlab)

我一直在尝试将一组化学方程式转换为 MATLAB 代码,以便能够求解不同的化学物质。我有近似的解决方案(因为它来自图表)但是在输入所有数据并多次检查后我仍然无法找到问题所在。我想知道出了什么问题,如果有人可以帮助我。 graph/equation 的来源是这篇link 的文章:The chemistry of co-injected BOE. 我稍后要重现的图表是论文中的图2,见下图:

现在我得到的10cc、40cc和90cc的结果分别是:
HF 43%, H2F2 48%, F- 3%, HF2- 6% 对比~28%, 63%, 2%, 7% (10cc).
HF 35%, H2F2 33%, F- 14%, HF2- 18% 对比~24%, 44%, 6%, 26% (40cc).
HF 21%, H2F2 12%, F- 37%, HF2- 30% 对比~18%, 23%, 20%, 45% (90cc).

脚本如下:

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-

需要一个额外的函数 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


文章中使用的方程式如下图所示:

(HF)2 在我的脚本中是 H2F2

看来问题不在于 Matlab,在这方面也有一些帮助。

可以在此处找到最终解决方案和更新的 Matlab 代码: https://chemistry.stackexchange.com/questions/98306/why-do-my-equilibrium-calculations-on-this-hf-nh4oh-buffer-system-not-match-thos