kill(all)$

maxima : intro to residue curve ( = simple batch distil,euler,multiple),version 1.0a ,part 6b
peter.vlasschaert@gmail.com,04/10/2020

ref 1:conceptual design of distillation systems (derivation :simple batch)
       M.F.Doherty & M.F.Malone
       McGRAW-Hill International Editions.
ref 2:systematic methods of chemical process design
      Lorenz T.Biegler,....................................,Arthur Westerberg
      Prentice Hall International series
ref 3:maxima : intro to residue curve ( = distillation curve,fixed α),version 1.0 ,part  1 to 5,(a)

ref3:my website.
other :data, wankat

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

Euler method 'simple batch distillation'
==============================
ref:maxima : intro to residue curve ( = distillation curve,fixed α),version 1.0 ,part  2

r1:'diff(x[i],z,1)=x[i]-y[i];

a) theoretical: simple batch , 'ODE's'  = " Ternary system "
--------------------------------------------------------------------------------
Ternary plot :A = index 1,B= index 2,C = index 3

r1a:lratsubst([i=1],r1);
r1b:lratsubst([i=2],r1);
r1c:lratsubst([i=3],r1);

→ rhs for ODEs

r1d:[rhs(r1a),rhs(r1b),rhs(r1c)];

→ constraint : x,y

r1d1:'sum(x[i],i,1,3)=1;
r1d2:'sum(y[i],i,1,3)=1;

reason : only need solve two ODE' s ,constraint{r1d1,r1d2}

rem : z=warped time
rem : α = relative volatilities ( = constant)

b)  How to solve system ODEs , numerical with Euler.
------------------------------------------------------------------------

theory : Euler Method
-------------------------------

 scheme x[i,k+1] = x[i,k]+h*(x[i,k]-y[i,k])
-----------------------------------------------------
 h = stepsize  , ( h ≤ 0.01),see  ref :1
 y = vapour
 x = liquid
------------------------------------------------------
otherwise : process simulators , Chemcad,..................

→  step 1: y & α

ref :maxima : intro to residue curve ( = distillation curve,fixed α),version 1.0 ,part 5

r2a:y[i]=(x[i]*α[i,3])/(x[2]*(α[2,3]-1)+x[1]*(α[1,3]-1)+1);

use : r1 & r2a

→ step 2:pratical ODEs

r2b:subst([r2a],r1);

system to solve ODEs:Euler method

r2c1:subst([i=1],r2b);
r2c2:subst([i=2],r2b);
r2c3:solve(ev(r1d1,nouns),x[3])[1];

rem :index A,1,B,2,C,3  & z,h (theoretical)

→ Euler equations :ODEs

r3:xb*(α2-1)+xa*(α1-1)+1;

r3a: xa = xa+h*(xa-(xa*α1)/r3);
r3b: xb = xb+h*(xb-(xb*α2)/r3);
r3c: xc = 1-(xa+xb);

h =step size.

c) source code :to solve {r3a,r3b,r3c}
---------------------------------------------------

ref : values,'example', x10 = xa ,x20=xb
maxima : intro to residue curve ( = distillation curve,fixed α),version 1.0 ,part 5
------------------------------------------------------------------------------------------------------------
ref : ' VBA  files',pg 285 to 290., 'use not the constraint in the code'
seperation process engineering ( thrid ed)
Philip C.Wankat
Pearson
data :
-------------------------------------------------------------------------------------------------------------

initial values:j=0 ( reboiler,bottom )
xa,xb,xc
constant :α1,α2

for i=1 to N
( j =i+1,
xa = xa+h*(xa-(xa*α1)/r3) 
xb = xb+h*(xb-(xb*α2)/r3)
xc=-xb-xa+1
)

ref :https://themaximalist.org/tag/euler-method/
===================================

(    /* Euler Method for  initial value problem
          y'=xy ,  y(0)=0.5   */ 
x0:0,
y0:.5,
h:.25,
nsteps:10,
xeuler:makelist(0,i,1,nsteps+1,1),
yeuler:makelist(0,i,1,nsteps+1,1),          
xeuler[1]:x0,
yeuler[1]:y0,
for i:1 thru nsteps do (
     xeuler[i+1]:xeuler[i]+h,
     yeuler[i+1]:yeuler[i]+h*xeuler[i]*yeuler[i]  
     )                                        
);

print (" Euler Method, simple batch distillation (A=>xA,B=>xB,C=>xC)")$

same :example data ;'distillation curve'
------------------------------------------------------

fpprintprec : 4;

N:1000  $ /*number of steps              */
NN:2 $ /*number of distillation curves */
h: 0.01 $ /* step size ,see above values */
/*parameter */
α1:2.4$  /* A = benzene */
α2:1.0$  /* B = toluene */
α3:0.21$ /* C = cumene */
/* data structure : output */
/*changed mutiple : distillation curve,matrix xA,xB,xC =double list*/
xA:makelist(makelist(0,i,1,N),j,1,NN)$
xB:makelist(makelist(0,i,1,N),j,1,NN)$
xC:makelist(makelist(0,i,1,N),j,1,NN)$
/*initial  conditions for components A,B,C,matrix=double list */
/*----------------------------------------------------------- */
/*initial condition : first distillation curve */
xA[1][1]:0.9905$
xB[1][1]:0.0085$
xC[1][1]:1-xA[1][1]-xB[1][1]$
/*initial condition : second distillation curve */
xA[2][1]:0.9875$
xB[2][1]:0.002$
xC[2][1]:1-xA[2][1]-xB[2][1]$
print (" → initial condition : A,B,C")$
print ("")$
xA$
xB$
xC$ 
for j:1 thru NN do(
      for i:1 thru N-1 do(
sumx :xB[j][i]*α2+xA[j][i]*α1+xC[j][i]*α3, 
xA[j][i+1] : xA[j][i]+h*(xA[j][i]-(xA[j][i]*α1)/sumx),
xB[j][i+1] : xB[j][i]+h*(xB[j][i]-(xB[j][i]*α2)/sumx),
xC[j][i+1] : xC[j][i]+h*(xC[j][i]-(xC[j][i]*α3)/sumx)
       )
                  )$
/* output : triangle triple (A,B,C) */
print (" → output  'triangle triple' : A,B,C")$
print ("")$
xA$
xB$
xC$

output ternary plot
==============

a) conversion : data 'Euler Method'.
------------------------------------------------

conversion :triangle triple (A,B,C) → cartesian coordinates(x_t,y_t)
ref :maxima : intro to residue curve ( = ternary plot,load draw),version 1.0b ,part 4a
      peter.vlasschaert@gmail.com,01/10/2020

x_t(aa_,bb_,cc_) :=-(2*bb_+aa_-2)/2$
y_t(aa_,bb_,cc_) := sqrt(3)/2*aa_  $

eu:[]$

/*number of distillation curves:NN */

for j:1 thru NN do(
for i:1 thru N do(
eu:endcons([x_t(xA[j][i],xB[j][i],xC[j][i]),y_t(xA[j][i],xB[j][i],xC[j][i])],eu)           
))$

eu:float(eu)$

hh[j]:=makelist(eu[N*(j-1)+i],i,1,N);

→ first : distillation curve ' cartesian coordinates' (initial condition,"see above"),

eu1:hh[1]$

→second: distillation curve ' cartesian coordinates'(initial condition,"see above")

eu2:hh[2]$

→ NN :'distillation curve ' 

mm:[]$
for j:1 thru NN do(
mm:endcons(hh[j],mm)    
)$

mm$

b) output : ternary plot
------------------------------

lp(a_,b_):=a_*t+b_*(1-t)$
f:sqrt(3)/20 $/* scaling factor,*/
s1x:[(k/sqrt(3))*f,((20-k)/sqrt(3))*f]$

s1y:k*f$ /* y values are fixed values */
s1:[k/sqrt(3)*t+(20-k)/sqrt(3)*(1-t),k]$
p1:makelist(makelist
  ([ev(lp(s1x[1],s1x[2]),nouns),ev(s1y,nouns)],
   t,0,1,0.1),k,0,10),float$

s2x:[((k+10)/sqrt(3))*f,((2*k)/sqrt(3))*f]$
s2y:[(10-k)*f,0*f]$
p2:makelist(makelist
([ev(lp(s2x[1],s2x[2]),nouns),ev(lp(s2y[1],s2y[2]),nouns)],
t,0,1,0.1),k,0,10)$

s3x:[((k)/sqrt(3))*f,((2*k)/sqrt(3))*f]$
s3y:[k*f,0*f]$
p3:makelist(makelist
([ev(lp(s3x[1],s3x[2]),nouns),ev(lp(s3y[1],s3y[2]),nouns)],
t,0,1,0.1),k,0,10)$

wxdraw2d( 
    font = Arial,
    font_size = 20,
    line_width =4,
    color = blue,
    point_size = 0.01, 
    points_joined = true,     
    makelist(points(p1[i]),i,1,length(p1)),    
    makelist(points(p2[i]),i,1,length(p2)),    
    makelist(points(p3[i]),i,1,length(p3)),
    /* here the source code ,eu :data ,see above */
    color = red,     
    makelist(points(mm[k]),k,1,NN), 
    label(["distillation curve",1/2,1/2]),
    /* end source code */
    /* info on the plot */
    color = black,   
    title = "Ternary plot {A,B,C}",
    color = red,
    label(["A",1/2+0.025,sqrt(3)/2-0.01]),
    label(["B",0+0.025,0+0.025]),
    label(["C",1-0.025,0+0.025]),       
    xtics =false,  /* removes numbers */
    ytics =false,  /* removes numbers */ 
    /* remove the frame     */
    axis_top    = false, 
    axis_bottom = false,
    axis_left   = false,
    axis_right  = false, 
    /* end remove the frame */   
    proportional_axes = 'xy 
         ),wxplot_size=[600,600]$