kill(all)$

maxima :theoretical :Bubble point method for Distillation,version 1.0
peter.vlasschaert@gmail.com,16/07/2020

load :packages
------------------

load(lrats)$

distillation column : with CMO
=====================
=====================

→ CMO = constant molar overflow (from tray to tray)

subscript
=======

i : component
j : tray( = plate)

parameter :
========
n = number of trays
nf = number of the feed tray

n:6;
nf:4;

Definitions
========

definition 1:equilibrium relation

y[i,j]:=ev(K[i,j]*x[i,j],nouns);

y[1,2];

definition 2:

p1:y[i,j]=v[i,j]/V[j];

definition 3:

p2a:x[i,j]=l[i,j]/L[j];

p2b:rhs(p2a)=lhs(p2a);

p3:lratsubst([p2a],p1);

p4:solve(p3,v[i,j])[1];

definition 4:

p5a:S[i,j]=K[i,j]*V[j]/L[j];

p5b:rhs(p5a)=lhs(p5a);

definition 5:

p6:lratsubst([p5b],p4);

condensor ( top column,total condensor ,'p7c' )
========

mass balance :total condensor
----------------------------------------

p7:V[2]=L[1]+D;
p7a:expand(p7/L[1]);

a)  stripping factor : j=1,use 'p5a'

p7b:lratsubst([j=1],lhs(p5a))=D/L[1];

p7b1:rhs(p7b)=lhs(p7b);

b) vapor(j=2),entering = liquid(j=1),leaving

p7c:v[i,2]/V[2]  = l[i,1]/L[1];

p7c1:solve(p7c,v[i,2])[1];

c)  combination : mass balance  &  stripping factor

p7d:lratsubst([p7b1],p7a);

d) combination : c)  & b)

→ first equation : tridiagonal system

p7e:lratsubst([p7d],p7c1);

p7e1:lhs(p7e)-rhs(p7e)=0;

e)  v[i,2]  as function l[i,2],S[i,2]

use :p6, tray(j=2)

p7f:lratsubst([j=2],p6);

f) combination :d) & e)

→ first equation : tridiagonal system

p7g:lratsubst([p7f],p7e1);

eq[i,1]:''p7g;

general tray=j , ( = not included , j≠ 1,nf,n)
============================
j=1 :eq[i,1], 'total condensor'

stage : j  (no reaction : in =out)
======
in    : molar  rate :L[i,j-1]  (down) ,V[i,j+1] (up)
out :  molar rate  : L[i,j],V[i,j]

with use of notation :l[i,j],v[i,j],mass balance around a  tray

p8a:l[i,j-1]+v[i,j+1]-l[i,j]-v[i,j]=0;

use :'p6'

p8a1:lratsubst([j=j+1],p6);

p8a2:p6;

p8a3:lratsubst([p8a1,p8a2],p8a);

→ equation : 'below' the feed tray  ( = nf) , j≠1,nf

for j:2 thru nf-1 do(
eq[i,j]:''p8a3
)$

eq[i,2];
eq[i,3];

→ equation : 'above' the feed tray  ( = nf),j≠n,nf

for j:nf+1 thru n-1 do(
eq[i,j]:''p8a3
)$

eq[i,5];

feed tray=j=nf
============================
j=nf :eq[i,nf], 'feed tray'

stage : j =nf  (no reaction : in =out)
======================================
in    : molar  rate :L[i,j-1]  (down) ,V[i,j+1] (up), F(i,nf)
out :  molar rate  : L[i,j],V[i,j]
rem : input saturated liquid, quality factor = 1
======================================
we use : p8a3 , rhs ≠ 0 => (-) sgn front (on the rhs)

p9:lhs(p8a3)= - f[i,j];

eq[i,j]:=''p9;

eq[i,nf];

last tray=j=n ,
============================
j=n :eq[i,n], 'reboilerr'

stage : j =n  (no reaction : in =out)
======================================
in    : molar  rate :L[i,j-1]  (down) 
out :  molar rate  : L[i,j],V[i,j]
rem : input saturated liquid, quality factor = 1
======================================
in -out = 0

p10:l[i,j-1]-l[i,j]-v[i,j]=0;

use :p6

p10a:lratsubst([p6,j=n],p10);

p10b:eq[i,n]:''p10a;

eq[i,n];

output : distillation equation
===================

print(eq[i,1])$

for j:2 thru n-1 do(
  print(eq[i,j])
)$

print(eq[i,n])$

makelist : equations eq[i,j], j=1,2..≤nf..n-1≤n

la:makelist(eq[i,j],j,1,n);

makelist : unknown l[i,j], j=1,2..≤nf..n-1≤n

lb:makelist(l[i,j],j,1,n);

tridiagonal :system ( = distillation)
==============

MA:coefmatrix(la,lb);

MB:augcoefmatrix(la,lb);

rem : augcoefmatrix, different sign (= last column,see equations 'la')
S[i,j]=(K[i,j]*V[j])/L[j]    i= component , j= tray(=plate)

noc : (n)umber (o)f (c)omponents

noc:2;

liquid phase composition  (stage =j),j=1..n
==============================

p11:x[i,j]:=l[i,j]/sum(l[i,j],i,1,noc);

x[1,1];

vapor phase composition  (stage =j),j=1..n
===============================

p12:v[i,j]/sum(v[i,j],i,1,noc);

p12a:y[i,j]:=(S[i,j]*l[i,j])/sum(S[i,j]*l[i,j],i,1,noc);

y[1,1];

vapor pressure, psat(i,j)
=================

p13:psat(T[j],a[i],b[i],c[i])=exp(a[i]-b[i]/(T[j]+c[i]));

p13a:rhs(p13);

equilibrium relation: K[i,j]
===================
P( = units KPa (example),other units ,)

K[i,j]:=''p13a/P;

Temperature
==========
P = pressure of all stages in the distillation column :
(P =101,325 kPa)
i= component
j= tray( =plate)

p14:pp[j]:=sum(x[i,j]*p13a,i,1,noc)-P = 0;

→T values can be found by : Newton- Rhaphson method ("nonlinear equations"),see :example.

pp[1];

pp[2];