maxima :tridiagonal : intro 1 :distillation multicomponent,version 1.0
peter.vlasschaert@gmail.com, 10/7/2019

kill(all)$

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

load(lrats)$

option:
=====

ratprint:false;

K(T) : vapor-liquid equilibrium ratio  , y,x
x : mole fraction (liquid phase)
y : mole fraction (vapor phase)
*******************************************

database : K(T)
===========

K(T) = K[i,j,α,β,γ,T]
i : component ( = products)
j : stage          ( = plate)

component : i=1..5 (rows:M)
------------------------
first       column :α ,for each i=1..5
second column :β ,for each i=1..5
third      column γ , for each i=1..5

→ 1e) data

M: matrix(
 [0.70,0.3e-2,0.6e-4], 
 [2.21,1.95e-2,0.6e-4], 
 [1.50,-1.6e-2,0.8e-4], 
 [0.86,-0.97e-2,0.46e-4], 
 [0.71,-0.87e-2,0.42e-4]
);

see:https://themaximalist.org/2017/02/28/a-little-maxima-function-to-find-the-dimensions-of-a-matrix/

definition : matlab like function
--------------------------------------

matsize(A1):=[length(A1),length(transpose(A1))];

p1:matsize(M);

rem : first       element : p1 (= rows)     =>  (n)umber (o)f (c)omponents ("noc")
         second element : p1 (= column)  => parameter (α,β,γ)

noc: p1[1];

→ 2e) T ( = Temerature °F = 'Farad')
--------------------------------------------------
number of stages : j : 1...10
----------------------------------------------------------------------
result : initial temperature guess for the column 'given'
=======================================
T[1]  =stage 1   : condensor : 'total'
T[10]=stage 10 : partial          'reboiler'
T[6]   = stage 6 : Feed Stage
=======================================

T:[132,164,172,180,188,196,204,212,220,228];

n = number of stages

n:matsize(T)[1];

function : K(T)
===========
i : component ( = products)
j : stage          ( = plate)

K[M,T,i,j]:= M[i,1]+M[i,2]*T[j]+M[i,3]*(T[j])^2;

example : input

K[M,T,2,1];

data :distilation column 
=================
--------------------------------------------------------------
F = flow rate of the feed : value
K = vapor-liquid equilibrium ratio  , y,x
L = flow rate of the liquid stream  (list: maxima)
V = flow rate of the vapor stream (list : maxima)
n      = number of equilibrium stages  'see above'
noc  = number of components           'see above'
w1   = flow rate of the overhead product,50
w2   = flow rate of the bottom     product,50
x : mole fraction (liquid phase), matrix
y : mole fraction (vapor phase), martrix
z :mole fraction (feed), list
-------------------------------------------------------------
subscript
------------
i : component ( = products)
j : stage          ( = plate)

property : (c)onstant (m)olal (o)verflow ( CMO)
rr : reflux ratio ( = 2.5)
-------------------------------------------------------------
F = w1+w2
-------------------------------------------------------------
w1 :top of column
w2 : bottom of column

rr :2.5;
w1 :50;
w2 :50;

F:w1+w2;

1e) initial guessed : liquid flow
---------------------------------------

L:makelist(0,i,1,n);

=>  position : feed tray : nf =6

nf:6;

a)  liquid flow : below , feed tray 
-----------------------------------------

for i:1 thru nf-1 do(
    L[i]:rr*w1
)$

L;

b)  liquid flow :above , feed tray (included)  & last tray (excluded)
-----------------------------------------------------------------------------------
last tray : top of the column

for i:nf thru n-1 do(
    L[i]:rr*w1+F
)$

L;

c)  liquid flow :last tray
-----------------------------------------------------------------------------------

L[n]:w2;

result : liquid flow rate profile of the column
===============================

L;

2e) initial guessed : vapor flow
---------------------------------------

V:makelist(0,i,1,n);

a) V[1] =0
--------------

V[1]:0;

b) vapor flow rate : 2 to n stage
----------------------------------------

for i:2 thru n do(
    V[i]:(1+1/rr)*L[1]
)$

result :vapor flow rate profile of the column
===============================

V;

3e) mole fractions from the component of the feed
------------------------------------------------------------------

z:makelist(0,i,1,noc);
z[1]:.055;
z[2]:.165;
z[3]:.227;
z[4]:.218;
z[5]:.335;

z;

check :sum of the mole fractions = 1
==========================

zz:sum(z[i],i,1,noc);

build :tridiagonal system
==================
rem : equations, component balances, 'around a tray  : input & output'

→ stage : j=1

eq[1]:-(L[1]+w1)*x[i,1]+V[2]*y[i,2]=0;

→ stage : j=2 to j = nf-1

for k:2 thru nf-1 do(
eq[k]:=L[k-1]*x[i,k-1]-(L[k]*x[i,k]+V[k]*y[i,k])+V[k+1]*y[i,k+1]=0  
)$

eq[2];

→ stage : j=nf '6'

eq[6]:L[5]*x[i,5]-(L[6]*x[i,6]+V[6]*y[i,6])+V[7]*y[i,7]=-F*mm;

→ stage : j=nf+1 to j = n-1

for k:nf+1 thru n-1 do(
eq[k]:=L[k-1]*x[i,k-1]-(L[k]*x[i,k]+V[k]*y[i,k])+V[k+1]*y[i,k+1]=0  
)$

→ stage : j=n

eq[10]:L[9]*x[i,9]-L[10]*x[i,10]-V[10]*y[i,10]=0;

equilibrium relation : y[i,j] =KK[i,j]*x[i,j]
i : 1..noc
j : 1..n
-----------------------------------------------------------------------
K[M,T,i,j]:= M[i,1]+M[i,2]*T[j]+M[i,3]*(T[j])^2  same KK[i,j]
-----------------------------------------------------------------------
KK[i,j]=α[i]+β[i]*T[j]+γ[i]*(T[j])^2

Q:makelist(eq[kk],kk,1,n);

length(Q);

p4:makelist(y[i,k1]=KK[i,k1]*x[i,k1],k1,1,n);

QQ:lratsubst(p4,Q);

QQ[6];

when i=1..noc → mm=z[1]..z[noc]

solution vector
==========

Q1:makelist(x[i,k2],k2,1,n);

lhs : equation
==========

Q2:coefmatrix(QQ,Q1);

KK[i,j]: vapor-liquid equilibrium ratio, y/x
----------------------------------------------------

Q3:makelist(KK[i,k2]=α[i]+β[i]*T[k2]+γ[i]*(T[k2])^2,k2,1,n);

Q4:lratsubst(Q3,Q2);

component : 1,  α[1] = M[1,1], β[1]=M[1,2],γ[1] =M[1,3] 
component : 2,  α[2] = M[2,1], β[2]=M[2,2],γ[2] =M[2,3] 
---------------------------------------------------------------------
→ component : k,  α[k] = M[k,1], β[k]=M[k,2],γ[k] =M[k,3] 

Q4a:lratsubst([i=k],Q4);

example : k=1
==========

Q4b:lratsubst([k=1],Q4a);

Q4b1:[α[1] = M[1,1], β[1]=M[1,2],γ[1]=M[1,3]]; 

Q4c:lratsubst(Q4b1,Q4b),float;

rhs : equation
==========
correction : component balance : feed tray (nf)
→ mm : i=1..noc, mm=z[1]..z[i]..z[noc]

Q5:makelist(d[i,k2],k2,1,n);

Q5:makelist(d[i,k3]=rhs(eq[k3]),k3,1,n);

Q5a:makelist(rhs(Q5[k4]),k4,1,n);

How to convert: list to matrix
-------------------------------------

Q5b:apply(matrix,[Q5a]);

Q5c:transpose(%);

→ component : k, mm= z[k]

example : k=1
==========

Q5d1:float(lratsubst([mm=z[1]],Q5c));

method : use inverse matrix: A
=======================
------------------------------------------------------------------------------
math         :  A*x = b  => A^(-1)*A*x=A^(-1)*b = I*x=A^(-1)*b
=>  x=A^(-1)*b 
maxima    
=>  x = invert(A)*b
------------------------------------------------------------------------------
b :column vector

Q6a:invert(Q4c).Q5d1;

Q6b1:lratsubst([i=1],Q1);
Q6b2:makelist(Q6b1[k5]=Q6a[k5][1],k5,1,n);

meaning :
=======

Q6a1:Q6a[1][1];

→ k=component=1 , stage =1

Q6b1:Q6a[n][1];

→ k=component=1 , stage =n=10

check :cvg 
========
---------------------------------------------------------------------------
1e) rhs: solve  Q4a, k=1,2,3..noc
      lhs : mm:= z[1],Z[2],..z[noc]
2e) method of inverse k =1 => Q6a
     repeat : k=2,k=3,...k = noc ,see 'above:example k=1'
3e)  sum (of all components : k=1...noc , j =stage )  = 1
--------------------------------------------------------------------------
meaning
=======
sum of all mole fraction of the components on stage = 1
------------------------------------------------------------------------
x[1,1]+x[2,1]+x[3,1]+x[4,1]+....+x[noc,1] = 1
------------------------------------------------------------------------
see above : x[1,1]    ,see below  
sum of all mole fraction of the components on stage =j
------------------------------------------------------------------------
x[1,j]+x[2,j]+x[3,j]+x[4,j]+....+x[noc,j] = 1
------------------------------------------------------------------------
see above : x[1,n]    ,see below  
============================================
when :sum ≠ 1 , need update 'Temperature profile'
→ bubblepoint subroutine ,appendix,later(→iterative process)
============================================

→ stage = 1,component = 1

Q6bc1:Q6b2[1];

→ stage = n,component = 1

Q6bc2:Q6b2[n];