kill(al1);

Maxima :prove :result Azeotrope  x[i]=y[i], VLE (binary),version 1.0
peter.vlasschaert@gmail.com,25/08/2020

use : MODIFIED  RAOULT'S LAW & Gibbs-Duhem Equation

 ==================================================
symbol :γ_1,γ_2 = activity coefficients ( = binary system)
symbol :Pvap[i]   = vapour pressure ( use :Antione Equation),i=1,2
symbol : P           = total pressure
===================================================
x[i] ,mole fraction 'liquid'
y[i] ,mole fraction 'vapour'
===================================================

load packages : option
==================

load(lrats)$
ratprint:false$

to use : lratsubst = > load packages

TH = Thermodynamics

Definition : binary system : Gibbs-Duhem (=GD) equation   (TH)
================================================

depends(γ_1,[x[1]]);
depends(γ_2,[x[1]]);

rem : error when use γ[i]

p1:x[1]*diff(log(γ_1),x[1])+x[2]*diff(log(γ_2),x[1])=0;

Definition : Raoult's Law ( =RL) (TH)
============================

p2:y[i]*P = x[i]*Pvap[i];

p2a:makelist(lratsubst([i=j],p2),j,1,2);

Definition : Modified  Raoult's Law (=MRL) (TH)
====================================
use : activity coefficients :γ_1,γ_2   (needed to use : GD)

p3a:y[1]*P =x[1]*γ_1*Pvap[1];
p3b:y[2]*P =x[2]*γ_2*Pvap[2];

General equation : mole fractions (=mf)
==============================

p4a:sum(x[i],i,1,2)=1;
p4b:sum(y[i],i,1,2)=1;

To prove : Azeotrope x[i] = y[i], binary VLE mixture
---------------------------------------------------------------------
VLE = VAPOUR LIQUID EQUILIBRIUM

1e)  step : sum (p3a & p3b)
--------------------------------------

p5a:factorsum(p3a+p3b);

2e) use  mf 
-----------------

p5b:lratsubst([p4a,p4b],p5a);

p5b1:lratsubst([p4b],p5a);

3e)  differentation of both sides with respect x[1]
-------------------------------------------------------------------

depends(P,x[1]);
p5c:diff(p5b,x[1]);

4e)  1-x[1] =x[2]  from mf
----------------------------------

p5d:solve(p4a,x[2])[1];
p5d1:[rhs(p5d)=lhs(p5d)];

p5d2:subst(p5d1,p5c);

identical replace :use   'Maxima :subst  '

5e) both sides :divide by  'γ_1*γ_2'
-------------------------------------------------

p5e:p5d2/(γ_1*γ_2);

p5e1:p5e,expand;

6e)  input : GD equation.
---------------------------------

p5f:args(p1);

p5f1:args(p5f[1]);

p5f2:(p5f1[2]= -p5f1[1]);

p5f3:solve(p1,diff(γ_1,x[1]));

p5f4:lhs(p5e1)=lratsubst(p5f3,rhs(p5e1));

p5f5:lhs(p5f4)=expand(rhs(p5f4));

p5f6:lhs(p5f5) = factorsum(rhs(p5f5));

7e) Azeotrope point  dp/dx[1] = 0 ,bcause γ_i≠0, i=1,2
--------------------------------------------------------------------------

p5f7:part(rhs(p5f6),1);

a) try to find  solution
================

=>   ≠ 0

p5f7a:part(rhs(p5f6),1,1);

=>   = 0

m1:sol 1

p5f7b:part(rhs(p5f6),1,2)=0;

find solution   → x[2] =0  then  γ_1=1 :  suppose symmetric
                       → x[1] =0  then  γ_2=1 :   

p6a:ode2(p5f7b,γ_2,x[1]);

solve for  %c

p6a1:lratsubst([γ_2=1,x[1]=0],p6a);

p6a2:rhs(p6a1)=lhs(p6a1);

p6a3:subst([p6a2],p6a);

p6a3:log(p6a3);

p6a4:subst([p5d],p6a3);

m2:sol 2
use :GD  into p5f7b

p6a5:solve(p1,diff(γ_2,x[1]));

p6a6:subst(p6a5,p5f7b);

p6a7:ode2(p6a6,γ_1,x[1]);

p6a8:limit(p6a7,x[1],0);

assume(%c >0);

physical impossible:

p6a9:lhs(p6a8) = ev(rhs(p6a8),nouns);

b) try to find the solution
==================

=>  =  0

p6b1:p5f7a =0;

 => ≠ 0

p6b2:p5f7b;

use : p6b1

p6b3:part(lhs(p6b1),2);

p6b4:p6b1+(-1)*p6b3;

Definition : equilibrium constant = K[i], i=1..2
==================================

p6b4a:(Pvap[1]*γ_1)/P=K[1];
p6b4b:(Pvap[2]*γ_2)/P=K[2];

p6b4c:makelist(K[i]= y[i]/x[i],i,1,2);

 ==================================

q1) divide :total pressure = P

p6b5:p6b4/P;

q2) use {p6b4a,p6b4b}

p6b6:lratsubst([p6b4a,p6b4b],p6b5);

q3) use (p6b4c)

p6b7:lratsubst(p6b4c,p6b6);

q4)  combination : q3) & mf 

p6b8:[y[2]=1-y[1],x[2]= 1-x[1]];

p6b8a:lratsubst(p6b8,p6b7);

p6b8a1:solve(p6b8a,y[1]);

p6b8b1:[y[1]=1-y[2],x[1]= 1-x[2]];
p6b8b2:lratsubst(p6b8b1,p6b7);
p6b8b3:solve(p6b8b2,y[2]);

We Have The Proved : Azeotrope: x[i] = y[i]

rem : relation between (p6b4:allowed ,p5b1)

p7a:p6b4;
p7b:p5b1;

p7c:lratsubst([p7a],p7b);

p7c1:lhs(p7c)=factorsum(rhs(p7c));

use : mf

p7c2:lratsubst([p4a],p7c1);

solve:γ_1

p7c3:solve(p7c2,γ_1);