我用立方型状态方程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

%     break;

%     end

end

plot(k12_record,F,'.')

%%