| 查看: 425 | 回复: 0 | ||
beyondyzs铜虫 (小有名气)
|
[求助]
二元作用参数
|
|
我用立方型状态方程P-R方程以及vdW-1混合规则进行气液相平衡计算，但是通过计算得到，无法求出二元相互作用参数，请各位高手及前辈指教。 clc; clear; close all; %% %êäèëêy¾Y N = 2; R = 8.314; Tc1 = 304.12; Tc2 = 774.30; omga1 = 0.228; omga2 = 0.786; T1 = 185.10; T2 = 600.8; Tr1 = T1 / Tc1; Tr2 = T2 / Tc2; Pc1 = 7.38*10^6; Pc2 = 1.56*10^6; P=[7.38 8.30 9.12 10.28 11.80 13.20 14.38 15.95].*10^6; T=328.2; y1=[0.9952 0.9913 0.9905 0.9887 0.9867 0.9851 0.9834 0.9816]; y2 = 1 - y1; x1=[0.4365 0.4788 0.5136 0.5507 0.6132 0.6521 0.7006 0.7998]; x2 = 1 - x1; %% %¼ÆËã a1 a2 b1 b2 (a1 a2¼ÆËã1«ê½óëÔ-1«ê½ËÆoõ2»ò»ÖÂ￡© k1 = 0.37464 + 1.54226 * omga1 - 0.26992 * omga1^2; k2 = 0.37464 + 1.54226 * omga2 - 0.26992 * omga2^2; a1 = 0.45724 * R^2 * Tc1^2 / Pc1 * (1 + k1 * (1 - sqrt(Tr1)))^2; a2 = 0.45724 * R^2 * Tc2^2 / Pc2 * (1 + k2 * (1 - sqrt(Tr2)))^2; b1 = 0.0778 * R * Tc1 / Pc1; b2 = 0.0778 * R * Tc2 / Pc2; %% k12=0; %¼ùéèk12 error=1; %¼ùéèÎó2î StepLength = 2/0.01; %% %½¨á￠·½3ì syms A B Z Vs; GasEq = Z.^ 3 - (1 - B).* Z.^2 + (A - 3.*B.^2 - 2*B).*Z - (A * B - B.^2 - B.^3); eq_mid = sym('eq_mid'); %3õê¼»ˉ±äá¿ Vs1_out=zeros(length(P),1); Vs2_out=zeros(length(P),1); Z1_out=zeros(length(P),1); Z2_out=zeros(length(P),1); ln_fai_l=zeros(length(P),1); ln_fai_g=zeros(length(P),1); fai_l=zeros(length(P),1); fai_g=zeros(length(P),1); F=zeros(length(1:1:StepLength),1); fg=zeros(length(1:1:StepLength),1); fl=zeros(length(1:1:StepLength),1); k12_record=zeros(length(1:1:StepLength),1); %% %Ñ-»·¿aê¼ for I=150:1:150 I %é趨k12μıä»ˉ k12 = k12 + 0.01; k12_record(I) = k12; %Çóêy×éam bm am1 = a1 .* y1.^2 + 2 .*sqrt( a1 .* a2 ).* (1 - k12) + a2 .* y2.^2; bm1 = b1 .* y1 + b2 .* y2; am2 = a1 .* x1.^2 + 2 .*sqrt( a1 .* a2 ).* (1 - k12) + a2 .* x2.^2; bm2 = b1 .* x1 + b2 .* x2; %′úèë·½3ì×é2￠Çó½a for j=1:1:length(P) %Çó½aÆøìå·½3ì eq_mid = subs(GasEq, A, am1(j) * P(j) / (R^2 * T^2)); eq_mid = subs(eq_mid, B, bm1(j) * P(j) / (R * T)); eq_mid = subs(eq_mid, Z, Vs * P(j) / (R * T)); t_mid = subs(vpa(solve(eq_mid, Vs),6)); for t0=1:1:length(t_mid) % Vs_out(j) = max(subs(vpa(solve(eq_mid, Vs),6))); %Vs¸ÃèçoÎè¡Öμ￡¿￡¿3ìDòÔYê±Ä¬èÏè¡×î′óÖμ % Vs_out(j) = min(subs(vpa(solve(eq_mid, Vs),6))); %Vs¸ÃèçoÎè¡Öμ￡¿￡¿3ìDòÔYê±Ä¬èÏè¡×î′óÖμ if isreal(t_mid(t0)) Vs1_out(j) = t_mid(t0); end end Z1_out(j) = Vs1_out(j) * P(j) / (R * T); A1_out = am1(j) * P(j) / (R^2 * T^2); B1_out = bm1(j) * P(j) / (R * T); %òoÏà eq_mid = subs(GasEq, A, am2(j) * P(j) / (R^2 * T^2)); eq_mid = subs(eq_mid, B, bm2(j) * P(j) / (R * T)); eq_mid = subs(eq_mid, Z, Vs * P(j) / (R * T)); t_mid = subs(vpa(solve(eq_mid, Vs),6)); for t0=1:1:length(t_mid) % Vs_out(j) = max(subs(vpa(solve(eq_mid, Vs),6))); %Vs¸ÃèçoÎè¡Öμ￡¿￡¿3ìDòÔYê±Ä¬èÏè¡×î′óÖμ % Vs_out(j) = min(subs(vpa(solve(eq_mid, Vs),6))); %Vs¸ÃèçoÎè¡Öμ￡¿￡¿3ìDòÔYê±Ä¬èÏè¡×î′óÖμ if isreal(t_mid(t0)) Vs2_out(j) = t_mid(t0); end end Z2_out(j) = Vs2_out(j) * P(j) / (R * T); A2_out = am2(j) * P(j) / (R^2 * T^2); B2_out = bm2(j) * P(j) / (R * T); %Çó½afai ln_fai_g(j) = b1 / bm1(j) * (Z1_out(j)-1)-... log(P(j) *(Vs1_out(j) -bm1(j))/R*T) +... am1(j)/(2*sqrt(2)*bm1(j)*R*T) *... ( b1/bm1(j)-2/am1(j)*(a1*y1+sqrt(a1*a2)*(1 - k12 )* y2))*... log((Vs1_out(j) + (1 + sqrt(2)) * bm1(j))/(Vs1_out(j)-(sqrt(2)+1)*bm1(j))); fai_g(j) = exp(ln_fai_g(j)); ln_fai_l(j) =b1 / bm1(j) * (Z2_out(j)-1)-... log(P(j) *(Vs2_out(j) -bm1(j))/R*T) +... am1(j)/(2*sqrt(2)*bm1(j)*R*T) *... ( b1/bm1(j)-2/am1(j)*(a1*x1+sqrt(a1*a2)*(1 - k12 )* x2))*... log((Vs2_out(j) + (1 + sqrt(2)) * bm1(j))/(Vs2_out(j)-(sqrt(2)+1)*bm1(j)));; fai_l(j) = exp(ln_fai_l(j)); end fg = P .* y1 .* fai_g'; fl = P .* x1 .* fai_l'; F(I) = sum(x1 .* fai_l' - y1 .* fai_g'); %èç1ûD¡óúÎó2î¾ííË3öÑ-»· % if F<error % break; % end end plot(k12_record,F,'.') %% |
回复此楼» 猜你喜欢
拟解决的关键科学问题还要不要写
已经有7人回复
-
请教限项目规定
已经有3人回复
-
存款400万可以在学校里躺平吗
已经有15人回复
-
Materials Today Chemistry审稿周期
已经有6人回复
-
基金委咋了？2026年的指南还没有出来？
已经有10人回复
-
基金申报
已经有6人回复
-
推荐一本书
已经有13人回复
-
国自然申请面上模板最新2026版出了吗？
已经有17人回复
-
纳米粒子粒径的测量
已经有8人回复
-
疑惑？
已经有5人回复