-
Notifications
You must be signed in to change notification settings - Fork 1
/
ChemicalEquilibrium.m
70 lines (70 loc) · 2.92 KB
/
ChemicalEquilibrium.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
%% 1:CO2 2:MEA 3:H2O 4:HCO3- 5:CO3-2 6:MEACOO- 7:MEAH+ 8:H3O+ 9:OH-
%% reaction1: 2H2O <===> H3O+ + OH-
%% reaction2: CO2 + 2H2O <===> H3O+ + HCO3-
%% reaction3: HCO3- + H2O <===> H3O+ + CO3-2
%% reaction4: MEACOO- + H2O <===> MEA + HCO3-
%% reaction5: MEAH+ + H2O <===> H3O+ + MEA
%% syms n1 n2 n3 n4 n5 n6 n7 n8 n9 nCO2 nMEA nH2O Keq1 Keq2 Keq3 Keq4 Keq5
% F=[n8*n9/n3^2-Keq1
% n8*n4/n1/n3^2*(n1+n2+n3+n4+n5+n6+n7+n8+n9)-Keq2
% n8*n5/n4/n3-Keq3
% n2*n4/n6/n3-Keq4
% n2*n8/n7/n3-Keq5
% n2+n6+n7-nMEA
% n1+n4+n5+n6-nCO2
% n3+n4+n5+n8+n9-nH2O
% -n4-2*n5-n6+n7+n8-n9];
% V=[n1 n2 n3 n4 n5 n6 n7 n8 n9];
% jacobian(F,V)
%% clear
% clc
% n0=[0.031606090373281 0.319253438113949 4.467258601553829];
% T=313;
%%
function [x n]=ChemicalEquilibrium(n0,T)
Keq1=exp(132.899-13445.9/T-22.4773*log(T));
Keq2=exp(231.456-12092.1/T-36.7816*log(T));
Keq3=exp(216.049-12431.7/T-35.4819*log(T));
Keq4=exp(2.8898-3635.09/T);
Keq5=exp(2.1211-8189.38/T-0.007484*T);
n=[n0(1) n0(2) n0(3) 0.0001 0.0001 0.0001 0.0001 0.0001 0.0001]';
nCO2=n0(1);nMEA=n0(2);nH2O=n0(3);
F=[n(8)*n(9)/n(3)^2-Keq1
n(8)*n(4)/n(1)/n(3)^2*(n(1)+n(2)+n(3)+n(4)+n(5)+n(6)+n(7)+n(8)+n(9))-Keq2
n(8)*n(5)/n(4)/n(3)-Keq3
n(2)*n(4)/n(6)/n(3)-Keq4
n(2)*n(8)/n(7)/n(3)-Keq5
n(2)+n(6)+n(7)-nMEA
n(1)+n(4)+n(5)+n(6)-nCO2
n(3)+n(4)+n(5)+n(8)+n(9)-nH2O
-n(4)-2*n(5)-n(6)+n(7)+n(8)-n(9)];
tol=10^(-5);
error=1;
while error > tol
J=[0 0 -(2*n(8)*n(9))/n(3)^3 0 0 0 0 n(9)/n(3)^2 n(8)/n(3)^2
(n(4)*n(8))/(n(1)*n(3)^2)-(n(4)*n(8)*(n(1)+n(2)+n(3)+n(4)+n(5)+n(6)+n(7)+n(8)+n(9)))/(n(1)^2*n(3)^2) (n(4)*n(8))/(n(1)*n(3)^2) (n(4)*n(8))/(n(1)*n(3)^2)-(2*n(4)*n(8)*(n(1)+n(2)+n(3)+n(4)+n(5)+n(6)+n(7)+n(8)+n(9)))/(n(1)*n(3)^3) (n(8)*(n(1)+n(2)+n(3)+n(4)+n(5)+n(6)+n(7)+n(8)+n(9)))/(n(1)*n(3)^2)+(n(4)*n(8))/(n(1)*n(3)^2) (n(4)*n(8))/(n(1)*n(3)^2) (n(4)*n(8))/(n(1)*n(3)^2) (n(4)*n(8))/(n(1)*n(3)^2) (n(4)*(n(1)+n(2)+n(3)+n(4)+n(5)+n(6)+n(7)+n(8)+n(9)))/(n(1)*n(3)^2)+(n(4)*n(8))/(n(1)*n(3)^2) (n(4)*n(8))/(n(1)*n(3)^2)
0 0 -(n(5)*n(8))/(n(3)^2*n(4)) -(n(5)*n(8))/(n(3)*n(4)^2) n(8)/(n(3)*n(4)) 0 0 n(5)/(n(3)*n(4)) 0
0 n(4)/(n(3)*n(6)) -(n(2)*n(4))/(n(3)^2*n(6)) n(2)/(n(3)*n(6)) 0 -(n(2)*n(4))/(n(3)*n(6)^2) 0 0 0
0 n(8)/(n(3)*n(7)) -(n(2)*n(8))/(n(3)^2*n(7)) 0 0 0 -(n(2)*n(8))/(n(3)*n(7)^2) n(2)/(n(3)*n(7)) 0
0 1 0 0 0 1 1 0 0
1 0 0 1 1 1 0 0 0
0 0 1 1 1 0 0 1 1
0 0 0 -1 -2 -1 1 1 -1];
%dn=J\(-F);
dn=pinv(J)*(-F);
n=n+dn;
F=[n(8)*n(9)/n(3)^2-Keq1
n(8)*n(4)/n(1)/n(3)^2*(n(1)+n(2)+n(3)+n(4)+n(5)+n(6)+n(7)+n(8)+n(9))-Keq2
n(8)*n(5)/n(4)/n(3)-Keq3
n(2)*n(4)/n(6)/n(3)-Keq4
n(2)*n(8)/n(7)/n(3)-Keq5
n(2)+n(6)+n(7)-nMEA
n(1)+n(4)+n(5)+n(6)-nCO2
n(3)+n(4)+n(5)+n(8)+n(9)-nH2O
-n(4)-2*n(5)-n(6)+n(7)+n(8)-n(9)];
error=max(abs(F));
end
x=zeros(9,1);
for i=1:9
x(i)=n(i)/sum(n);
end