kill(all)$

maxima:introduction : θ method ' distillation , part 1,version 1.0
peter.vlasschaert@gmail.com,14/10/2020

load packages ,option
-------------------------------

load(lrats)$ load(draw)$ ratprint:false$

1e ) step, notation
==============

a)  symbol 
---------------

pa: X[i,j];
print( " stage          =",j)$/*nl :trap,en:plate,tray */
print( " component  =",i)$    /*product*/

rem : example  X = molefraction , x[i,j] = liquid mole fraction , X → x
,                                                         y[i,j] = vapor mole fraction,  X→ y 
rem : j = j th stage of distillation column

pb:X[j]; /*Flow rate,stage j */
print("Flow rate @ stage ",j)$

pc:X; /*Flow rate */
print("Flow rate ")$

b)  notation    i th component  ,one feed entrance = X = F
------------------------------------------------------------------------------------
                        bottom column :xb = x[i,N+1]  →  B  ( =bottom)
                        top       column  xt = x[ i,1]        →  D  (=destillate)
                        feed     column  xf  = x[ i,f ]       →   F  (=feed)
rem  j:  N+1  <  f  < 1< 0 , bottom →  top ' distillation column'
--------------------------------------------------------------------------------------------
j = 0 ,accumulator , outside column  in combination with condensor 
f-1 = below feed  stage
f+1 = above feed  stage
-------------------------------------------------------------------------------------------
use : convention  :CHARLES D.Holland, ji → ij
-------------------------------------------------------------------------------------------
nc : number of components

p1a:[xb = x[i,N+1],B];
p1b:[xt = x[i,1],D];
p1c:[xf = x[i,f],F];  

→ more than one entrance : F[j],notation ( notation : in & out)
→ p1a,p1b,p1c : list notation in maxima , list =[], set = {},list sart : index 1

c)  mass balance ,component balance

→ Mass balance

p2a: p1c[2] = p1b[2]+p1a[2];

→ Componet balance

p2b:lhs(p2a)*rhs(p1c[1])=part(rhs(p2a),1)*rhs(p1b[1])+part(rhs(p2a),2)*rhs(p1a[1]);

2e ) step  theta method  'notation'
=========================
co = corrected
ca = calculated

a) set values
------------------

p3a: part(rhs(p2b),1)=b[i,co];
p3b: part(rhs(p2b),2)=d[i,co];

b) new mass balance
------------------------------
use :p2b

p3c:lratsubst([p3a,p3b],p2b);

c) constraint : mass balance
---------------------------------------
specified : p3b, summation of components.

p3d:sum(p3b,i,1,nc);
p3d1:reverse(p3d);

def :general equation : for every j -stage  → equality
-----------------------------------------------------------------------

p3:sum(x[i,j],i,1,nc)= 1;

p33a:print(" every stage = j")$

p3d2:subst([j=1],p3);

→ def :final constraint

p3d3:subst([p3d2],p3d1);

3e) step  theta equality ' definition'
===========================
co =corrected
ca = calculated , 'tridiagonale system : equations'

→ def :theta aquality

p4:b[i,co]/d[i,co] = θ*b[i,ca]/d[i,ca];

4e) step  use : component balance , θ,find d[i,co]
==================================:===

a)  mass balance :solve for d[i,co]
-----------------------------------------------
use :p3c

p5a:expand(p3c/d[i,co]);

b) θ definition
--------------------
relation between : co & ca

p5b:expand(lratsubst([p4],p5a));

c) find : d[i,co]
--------------------

p5c:solve(p5b,d[i,co])[1];

p5c1:part(rhs(p5c),1)/d[i,ca];

p5c2:expand(part(rhs(p5c),2)/d[i,ca]);

→ final d[i,co]

p5c3:lhs(p5c)= p5c1/p5c2;

5e) step use :final d[i,co], apply constraint p3d3
====================================

p6a:sum(p5c3,i,1,nc);

p6b:subst([p3d3],p6a);

6e) step use : p6b ,construct function  G(θ)
=================================
→ find θ

p7a:G(θ) = rhs(p6b)-D;

→ final def : ratio m[i,ca]

p7a1:m[i,ca]=lratsubst([θ=1],rhs(p4));

p7a2:reverse(p7a1);

p7a3:lratsubst([p7a2],p7a);

p7a31:part(rhs(p7a3),1,2,1);

p7a32:part(rhs(p7a3),1,1);

p7a33:part(rhs(p7a3),2);

7e) step use : G(θ) → G'(θ): needed: newton method
========================================
→ find θ

p7b:diff(p7a,θ);

→ ratio 'm'

p7b1:lratsubst([p7a2],p7b);

?  How  use this  : inside  G'(θ)
-------------------------------------------

factorsum(m[i,ca]^2*θ^2+2*m[i,ca]*θ+1);

→ analysis : rhs' p7b1'

p8a:part(rhs(p7b1),1,2);

p8aa1:part(p8a,0);
p8aa2:part(p8a,1);

p8a1:args(p8a);

p8a2:factorsum(p8a1);

p7b2:lratsubst([p8aa2=p8a2[1]],p7b1);

8e ) step :newton method 
------------------------------------

→ old values :θ , k: index, rhs

p9a:[θ = θ[k]];

→ new values :θ , k+1: index,lhs

p9b:subst([k=k+1],p9a);

→  symbolic, ' newton'  , find from taylor expansion
---------------------------------------------------------------------

p9c:θ = θ - lhs(p7a3)/lhs(p7b2);

p9c1:subst([p9b],lhs(p9c)) = subst([p9a],rhs(p9c));

→ computation , 'newton',values
--------------------------------------------
use :p7a3,p7b2 & p9c1

p9c2:subst([p9a],p7a3);
p9c3:subst([p9a],p7b2);

p9c4:θ[k+1] = θ[k] -1* rhs(p9c2)/rhs(p9c3);

analysis : G(θ)
--------------------
use :p7a3

assume(F>0,m[i,ca]>0,D>0,nC>0,x[i,f]>0);

→ final value , G(θ) ,θ → ∞
-------------------------------------

p10a:limit(rhs(p7a3),θ,inf) ;

use :p7a31,remove summation.

p10b:p7a32*p7a31+p7a33;

p10c:'limit(p10b,θ,inf)=limit(p10b,θ,inf);

def :final value : θ →∞

p10c1:'limit(G(θ),θ,inf) = rhs(p10c);

→ initial value, G(θ),θ → 0
-------------------------------------

p10d:'sum('limit(p10b,θ,0),i,1,nc)=sum(limit(p10b,θ,0),i,1,nc);

p10d1:lhs(p10d) = expand(ev(rhs(p10d),nouns));

→ general constraint : mole fraction, p3

p10d2:subst([j=f],p3);

→ use :p10d2,rhs(p10d1)

p10d3:lhs(p10d)= F*lhs(p10d2)-D;

 ?  problem : p10d1 to p10d3, F & x[i,f] , 'change properties of sum ?'

p10d4:lratsubst([p10d2],p10d3);

→ use  c)  mass balance ,component balance,see above,p2a

p10d5:p2a-D;/* general mass balance */

p10d6:lhs(p10d3)=lratsubst([p10d5],rhs(p10d4));

def :initial value : θ →0

p10d7:limit(G(θ),θ,0) = rhs(p10d6);

→   B >  θ >  -D   ,assume B >0,D > 0

=>  Newton's method :cvg positive root of G(θ),θ =0,start value ( k=0)

p11:subst([k=0],p9c4);

p11a:subst([θ[0]=0],p11);