kill(all)$

virial equation of state : binary mixtures.  'version 1.0'  (thermodynamics :part 2a)
peter.vlasschaert@gmail.com,13/07/2017

B_[i,i]: second virial component 'i'  pure component, i=1..2
B_[j,j]: second virial component 'i'  pure component, j=1..2
B_[i,j] =B_[j,i]: second virial component , i<>j 'interaction'
given : y_i   ,i=1..2           '       sum(y_i)=1       '

p1:B=sum(sum(y_[i]*y_[j]*B_[i,j], i, 1, 2), j, 1, 2);

p2:B_[2,1]=B_[1,2];

p3:ev(p1,p2);

The compressibility factor of the mixture

p4:Z=1+(B*P)/(R*T);

p5:ev(p4,p3);

l_r i= l/l_ci
index : r = reduced
            c = critical
            i = component
l = T,P            (temperature : Kelvin ,pressure)
omega           (acentric factor =function of the  reduced sat pressure by fixed  reduced temperature)
omega =0     (simple fluid)
rem : omega = -1-log10(P_r'sat')    '(fixed by  T_r)'
rem : P_sat                                       '( vapour pressure)'
ln (P_sat )  = A-B/(T+C)                   '(Antoine equation)'
rem : use table (T  ?  °C,K) 'check'

Z_c = (P_c*V_c)/(R*T_c)    (critical compressibility factor),index : c=critical
Z = Z(T_r,P_r,omega)         (Pitzer and co..)
Z = Z(T_r,P_r,Z_c)               (Lyndersen and co..)

given :  i = component ,T_ci,V_c,Z_ci,omega_i    , i=1..2
             mixture :  B_[i,i] at temperature = T

p6:T_c12 = sqrt(T_c1*T_c2);

p7:ω_12=(ω_1+ω_2)/2;

p8:Z_c12=(Z_c1+Z_c2)/2;

p9:V_c12=((V_c1^(1/3)+V_c2^(1/3))/2)^3;

p10:P_c12=(Z_c12*R*T_c12)/V_c12;

p11:ev(p10,p6,p8,p9);

p12:T_r12= T/T_c12;

p13:ev(p12,p6);

given  : B0 = 0.083-0.422/(T_r12)^(1.6)
             B1 = 0.139-0.172/(T_r12)^(4.2)

p14:B0 = 0.083-0.422/(T_r12)^(1.6);
p15: B1 = 0.139-0.172/(T_r12)^(4.2);

p141:ev(p14,p13);
p151:ev(p15,p13);

find :  B_[1,2]
=>  (B_[1,2]*P_c12)/(R.T_c12) = B0+ω_12*B1

p16: (B_[1,2]*P_c12)/(R.T_c12) = B0+ω_12*B1;
p17:part(solve(p16,B_[1,2]),1);

p18:ev(p17,p6,p7,p11,p141,p151);

p19:ev(p5,p18);

molar volume : V=B+RT/P ( =>  when   Z=(B*P)/(R*T)+1)

p20:V=B+RT/P;
p21: ev(p20,p3);

rem: V=B+R.T/P+(C-B^2)*P/RT ( => when  Z=1+B*P/RT+(C-B^2)*(P/RT)^2)
          C = given
rem: dimension volume = R.T/P (one mole)
=>  P.V/(RT)=1 => V= (R*T)/P
rem:  eq : 'p4'    =>   Z = B/V+1

p22: V=R.T/P;

 V replace by  (R.T)/P   in  'p4'

p23:part(p4,1)=ratsubst(part(p22,1),(R*T)/P,part(p4,2));

 Z= (P.V)/(R.T)

p24:Z= (P.V)/(R.T);

 'p23' = 'p24'

p25:part(p23,2)-part(p24,2) = 0;
p26:ratsimp(p25);
p27:(-1)*expand(num(part(p26,1)));
p28:denom(part(p26,1));

we want solve : quadratic equation
y = a*V^2+b*V+c=0

V=-(sqrt(b^2-4*a*c)+b)/(2*a),v=(sqrt(b^2-4*a*c)-b)/(2*a)

p29:p27=0;