kill(all)$

multicomponent flash distillation ,version 1.0
peter.vlasschaert@gmail.com,11/08/2020

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

load(lrats)$
ratprint:false$

ref : website ( = software ,'drawing'  = Geogebra , Paint from Microsoft)
===================================================
online use :https://www.geogebra.org/m/MvdP8dju

F = Feed Flowrate ,z[i] = mole fraction
V = Vapor Flowrate,y[i] = mole fraction
L = Liquid Flowrate,x[i] = mole fraction
============================
NC = Number of Components

p1:Input1 = Output1+Output2;

info : sum of all input = sum of all output, see CutSet (in Figure 1)
Cut Set = general for generation of equations, 'see below'

1e) Mass Balance ( use CutSet)
=======================

pa1:lratsubst([Input1=F,Output1=L,Output2=V],p1);

2e) Component Mass Balance (use CutSet)
================================

pa2:lratsubst([Input1=F*z[i],Output1=L*x[i],Output2=V*y[i]],p1);

3e) constraint : mole fraction.
=====================

pa31:sum(x[i],i,1,NC) = 1;
pa32:sum(y[i],i,1,NC) = 1;
pa33:sum(z[i],i,1,NC) = 1;

pa3:[pa31,pa32,pa33];

4e) equilibrium equation = K[i]
======================
relation : between ,x[i] & y[i]

pa4:y[i]=K[i]*x[i];

pa41:sum(lhs(pa4),i,1,NC) = sum(rhs(pa4),i,1,NC);

use :3e)

pa42:lratsubst(pa3,pa41);

pa43:[rhs(pa42) = lhs(pa42)];

rem :[ ], list

5e) find x[i]  = f ( ψ =V/F,z[i],K[i])
=======================

step a: remove L ,pa1 & pa2
-------------------------------------

pa51:solve(pa1,L);

pa52:lratsubst(pa51,pa2);

step b:use 4e)
-------------------

pa53:lratsubst([pa4],pa52);

step c:divide (step b)/F ,expand
-----------------------------------------

pa54:pa53/F,expand;

step d:replace V/F by ψ in step c
-------------------------------------------

pa55:lratsubst([V/F=ψ],pa54);

step e: solve step d for x[i]
----------------------------------

pa56:solve(pa55,x[i])[1];

6e) take sum both site 'pa56' use :constraint (3e)
===================================

pa6:sum(pa56,i,1,NC);

pa61:lratsubst(pa3,pa6);

pa62:rhs(pa61)=lhs(pa61);

7e) find output  y[i]
=============

a) step : use pa56 *k[i]
-----------------------------

pa71:pa56*K[i];

b) step : use pa4
-----------------------------

pa72:lratsubst([rhs(pa4)= lhs(pa4)],pa71);

8e)  find Rachford–Rice equation
========================
use :pa72

ref :https://chempedia.info/info/rachford_rice_equation/

a) step : use sum operator on pa72,use : constraint(=pa3)
---------------------------------------------------------------------------

pa8:sum(pa72,i,1,NC);

pa81:lratsubst(pa3,pa8);

pa81a:args(rhs(pa81));

b) step :use pa6 & pa61
--------------------------------

pa82:pa61;

pa82a:args(rhs(pa82));

DEF :Rachford–Rice equation
======================

pa83:lhs(pa81-pa82)=sum(pa81a[1]-pa82a[1],i,1,last(pa81a));

reason complicated : process :simple gives ' error from  Maxima'
: p81a-pa82 ,see 'above' ,latest vesion of Maxima ,see Website 'Maxima'

pa84:factorsum(pa83);

ψ = V/F ( fraction of the feed vaporized)
=============================
rem : Flash Calculation :  0 ≤ ψ ≤ 1
1e) step : guess V/F 
2e) step  : use check relation pa61, when satisfied : you find solution

Rem : Calculation of K[i]  , given P (see Figure 1:total pressure inside the Flash  Drum)
                                            given Psat ( from : Antoine Equation (pa91) : find Tables for pure components)
                                              , only unknown T ( = given)

pa91:psat(i,A,B,C,T):=exp(A+B/(T-C));

rem :{A,B,C} triple : unic for chemical species.

rem : see , 'ref : book chemcal engineering ' ,for units , T(°C,R,T,F,..),psat ,( Pa,kPa,mmHg,..) ≠
----------------------------------------------------------------------------------------------------------------------------
(Antoine Equation:see ref book)

ideal solutions
==========

pa92:K[i]=psat(i,A,B,C,T)/P;

noni deal solutions
==============

pa92:K[i]=(psat(i,A,B,C,T)*γ[i])/P;

γ[i] = Acitivity Coefficient 
==================
 (anolog use: acid -base problem solving :see my website ,'analytical chemistry')
=====================================================================
Modelling  : -The Margulus Model
                   - The Vanlaar   Model
                   - The Wilson    Model
                   - The  NRTL     Model
                   - The  UNiQUAC  Model
                   - The Hildebrand  Model 
reason use  : Chemical plant simulators : Aspenplus ,KBC,Pro||,Hysys,Chemcad,
                       Simcentral ,coco ,....( used Chemical Factories,most of them are very expensive).   
=====================================================================                                                 

Given : F ,z[i],T,P   =>  find(V;L,y[i],x[i]), i=1..NC
==================================
suppose :  ψ , see above.

implimentation :see my website.(simple : calculation)

pa93:lratsubst([pa92],pa84);

use :pa51

pa94:V=ψ*F;
pa95:lratsubst([pa94],pa51[1]);

find : x[i],y[i] see pa56 & pa72

pa96:lratsubst([pa92],pa56);
pa97:lratsubst([pa92],pa72);

very important : How to solve Nonlinear equation (numerical), see pa84
use some way to find solution :cvg (instruction in your code to do optimalisation).