kill(all)$

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

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

option 
====

ratprint:false$

load packages
===========

→ general:

load(lrats)$

→special

1e) numerical : nonlinear equations
---------------------------------------------
Function: mnewton (FuncList,VarList,GuessList)

               (fpprec : 25, newtonepsilon : bfloat(10^(-fpprec+5)))$
example :mnewton([x1+3*log(x1)-x2^2, 2*x1^2-x1*x2-5*x1+1],[x1, x2], [5, 5])           

 load(mnewton)$

definition
======

see :https://stackoverflow.com/questions/43870997/how-to-get-vector-from-matrix-in-maxima

This function returns the i-th row of a matrix similar to Matlabs M(i,:) 

This function returns the i-th column of a matrix similar to Matlabs M(:,i) 

col(M,i) := 
    block(
        return(transpose(transpose(M)[i])) 
    )$

row(M,i) := 
    block(
        return(transpose(transpose(M[i])))
    )$

Distillation example :(b)uble (p)oint method, bp-method
=======================================
n = number of trays
nf = number ot the Feed tray
noc = number of components

st:1;
n :4;
nf:3;
noc:2;

D = destillate (=flow rate)
---------------------------------

D:0.5;

F= Feed (= flow rate)
-----------------------------

F:1;

rr=reflux ratio

rr:1.5;

pa:rr=L[1]/D;

pb:solve(pa,L[1])[1];

1e) general : tridiagonal system
-----------------------------------------

p1:eq[st]:=A[st]*l[i,st-1]+B[st]*l[i,st]+C[st]*l[i,st+1]=DD[i,st]$

for k:2 thru n-1 do(
    p2:eq[k]:=A[k]*l[i,k-1]+B[k]*l[i,k]+C[k]*l[i,k+1]=DD[i,k]
)$

p3:eq[n]:=A[n]*l[i,n-1]+B[n]*l[i,n]+C[n]*l[i,n+1]=DD[i,n]$

a) list of equations:

leq:makelist(eq[ka],ka,1,n);
leq1:lratsubst([l[i,0]=0,l[i,n+1]=0],leq);

correction :'leq1' for Feed Tray ( j = nf)

→  rhs: from equations , 'above'

la:makelist(if j=nf then 
   DD[i,nf]=-f[i,j] else
   DD[i,j]=0,j,1,n);

→ new equations :after correction

leq2:lratsubst(la,leq1);

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

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

c) general A[i],B[i],C[i]

MA:coefmatrix(leq2,lvar);

MB:augcoefmatrix(leq2,lvar);

for j:1 thru n do(
B[j]:-S[i,j]-1)$

for j:2 thru n do(
A[j]:1    
)$

for j:1 thru n-1 do(
C[j]:S[i,j+1]    
)$

stripping factor:S[i,1]
---------------------------

S[i,1]:lratsubst([pb],D/L[1]);

MA1:ev(MA,nouns);

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

MA2:ev(MA1,nouns);

see :theoretical buble point method ,'distillation'

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

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

MA3[i]:=ev(MA2,nouns);

componenten: i =1,2,..,noc

first component: i=1 ,j =1,2,..,n,feed tray (j=nf)
=================================

MC1:ev(MA3[1],nouns);

physical :properties :

MC1a:[a[1]:14.3916,b[1]:2795.82,c[1]:230.00];

MC1b:ev(MC1,nouns);

1) vapor flow rate 
=============

a)  around first tray (j=1),top of column

pv1:V[2]=lratsubst([pb],L[1]+D);

b) other values : V[j], j = 3,4,..n

result : vapor flow rate
----------------------------

pv2:makelist(V[j]=rhs(pv1),j,2,n);

2) liquid flow rate 
=============

a) see 'above'

pl1:[pb];

b) other values : L[2], circle around :combine (first tray,second tray)

pl2:[L[2]=lratsubst(pv2,V[3]-D)];

c) last tray (j=n)

c1)  circle around:combine  'all trays'

pl3:F = D+L[n];

pl3a:solve(pl3,L[n]);

c2) circle around 'last tray'

p1ll:flatten(append(pv2,pl3a));

pl3b:[L[3]=lratsubst(p1ll,V[4]+L[4])];

pll:flatten([pl1,pl2,pl3a,pl3b]);

combine list :pv2 & pll

pvl:flatten([pv2,pll]);

MC1c:lratsubst(pvl,MC1b);

P :pressure of 'all stages':distillation column
================================
P( = kPa units)

pk1:[P=101.325];

T[j], j=1..n,assumed values for temperature profile,units :°C

pk2:makelist(T[j]=50+(j-1)*10,j,1,n);

pk:flatten([pk1,pk2]);

MC1d:lratsubst(pk,MC1c);

MC1e:float(MC1d);

lvar1[i]:=''lvar;

lvar1[1];

f1:0.5;

ff:f[1,nf]=f1;

MB1:lratsubst([i=1],MB);

MB1a:lratsubst([ff],MB1);

to find :column ,you need to change sign:see' lqe2'

pc1:-col(MB1a,n+1);

--------------------------------------------------------------------
math      : A*x = b, A = matrix
maxima :  x = A^(-1)*b =invert(A)*b
--------------------------------------------------------------------
info :x,b :column vectors
======================================
see : appendix : 'Thomas Method',tridiagonal system.
======================================

res1:invert(MC1e).pc1;

res1a:flatten(args(res1));

res1b:maplist("=",lvar1[1],res1a);

second component: i=2 ,j =1,2,..,n,feed tray (j=nf)
====================================

physical properties:

MC2a:[a[2]:16.2620,b[2]:3799.89,c[2]:226.35];

MC2b:ev(MA3[2],nouns);

MC2c:lratsubst(pvl,MC2b);

MC2d:lratsubst(pk,MC2c);

MC2e:float(MC2d);

list unknown : i=2

h1:lvar1[2];

f2:0.5;
ff1:f[2,nf]=f2;

MB2:lratsubst([i=2],MB);

MB2a:lratsubst([ff1],MB2);

to find :column ,you need to change sign:see' lqe2'

pc2:-col(MB2a,n+1);

res2:invert(MC2e).pc2;

res2a:flatten(args(res2));

res2b:maplist("=",lvar1[2],res2a);

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

component: i=1
===========

pxa1:lratsubst([i=1],pxa);

use :res1b,res2b

resx:flatten([res1b,res2b]);

pxaa:lratsubst(resx,pxa1),float;

component :i=2
===========

pxa2:lratsubst([i=2],pxa);

pxbb:lratsubst(resx,pxa2),float;

pxcc:flatten([pxaa,pxbb]);

check :stage:= j, x[i,j]  ,i=1,2,..,noc
=========================

stage[j]:=lratsubst(pxcc,sum(x[i,j],i,1,noc));

ev(stage[1],nouns),float;

mt:makelist(ev(stage[1],nouns),i,1,n),float;

update : Temperature profile
=====================

orginal :see 'above'
-------------------------

pt1:pk2;

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

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

p13a:rhs(p13);

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

stage j=1     (top :distillation column)
=======

z1:ev(pp[1],nouns);

z1a:lratsubst(pxcc,z1),float;

z1b:lratsubst([P=101.325,j=1],z1a);

use : package :'see above',numerical solve , 'nonlinear equation'
==============================================

z1c:mnewton(z1b,T[1], 50)[1] ;

stage j=2
=======

z2:ev(pp[2],nouns)$
z2a:lratsubst(pxcc,z2),float$
z2b:lratsubst([P=101.325,j=2],z2a)$
z2c:mnewton(z2b,T[2],rhs(pt[2]) )[1] ;

stage j=3  (feed tray)
=======

z3:ev(pp[3],nouns)$
z3a:lratsubst(pxcc,z3),float$
z3b:lratsubst([P=101.325,j=3],z3a)$
z3c:mnewton(z3b,T[3],rhs(pt[3]) )[1] ;

stage j=4  (last tray), bottom distillation column
=======

z4:ev(pp[4],nouns)$
z4a:lratsubst(pxcc,z4),float$
z4b:lratsubst([P=101.325,j=4],z4a)$
z4c:mnewton(z4b,T[4],rhs(pt[4]) )[1] ;

update :temperature profile
------------------------------------

upt:flatten([z1c,z2c,z3c,z4c]);

repeat the equation above: for updated temperature profile.(see :check,'above' ≠ 1.0)