kill(all)$

maxima :concentration respons for distilation column ,version 1.0
peter.vlasschaert@gmail.com,26/08/2020

ref : chemical and energy process engineering 
        sigurd skogestad
        crc  press
see page : p314-315

load packages & option
===================

load(lrats)$
ratprint :false $

y[i],relative volatility =α
=================
y[i]=f(α,x[i])

definition : vapour -liquid 
===================

m(α,i):=(α*x[i])/(1+(α-1)*x[i]);

y[i], i = 1,2,..4

m1:makelist(y[i]:m(α,i),i,1,3);

y[4]:x[4]; /*total condensor */

differential equations : distillation column (four stage :distillation : column)
========================================================

q1[1]:dxdt(1) = L2*x[2]-V1*y[1]-B*x[1];         /*bottom of column */
q1[2]:dxdt(2) = L3*x[3]+V1*y[1]-L2*x[2]-V2*y[2]+F*zf;  /*feed tray*/
q1[3]:dxdt(3) = L4*x[4]+V2*y[2]-L3*x[3]-V3*y[3];/*above feed tray*/
q1[4]:dxdt(4) = V3*y[3]-L4*x[4]-D*y[3];        /*top of column     */

m2:makelist(rhs(q1[i]),i,1,4);

m2a:listofvars(m2);

input : differential for distillation
========================

→ constant

α:4.78 $     /*constant : relative volatility */
F:1$         /* Feed Rate                     */
zf:0.5 $     /*Feed Composition               */
;

→  flows in column

=>  Vapour flow

V:3.55 $ V1:V $ V2:V $ V3:V $

 => Liquid flow

L:3.05 $ L4:L $ L3:L $ L2:L+F $

→  equation  Bottom,Distillate

D : V3-L4$
B : L2-V1$

differential equations with (input)
=========================

m3:m2;

m3a:ev(m3,nouns);

m3b:listofvars(m3a);

start condition : distillation ( ≠ IC intial condition)
====================================

m4:[0.5,0.5,0.5,0.5];

rem : length list  (m3b) = length list (m4)

solve : system of differential equations ( 'm3a')
====================================

how to use :
=========
Function: rk
    rk (ODE, var, initial, domain)
    rk ([ODE1, …, ODEm], [v1, …, vm], [init1, …, initm], domain)

→ domain

tbegin : 0;
tend   : 25;
stepvalue : 0.1;

domain:[t,tbegin,tend,stepvalue];

rk(m3a,m3b,m4,domain); error 
→  x[1] ,x[2],x[3],x[4] ,not allowed.

change : x[i] => xi , i=1..length(m3b)

new differential : equations (use :m3a)
=============================

m5a:makelist(x[i]=concat(x,i),i,1,length(m3b));

m5b:lratsubst(m5a,m3a),ratsimp;

m5c:listofvars(m5b);

solve : new differential equations (m5b)
==============================

sol:rk(m5b,m5c,m4,domain)$

1e) output : plot2d
==============

IC : [0.5,0.5,0.5,0.5];

How to plot 2d   : 
==========================================================
  → plot result ( xaxis = t ,yaxis =y)
given     : dx/dt = 4-x^2-4*y^2     dy/dt = y^2-x^2+1
solution : sol: rk([4-x^2-4*y^2,y^2-x^2+1],[x,y],[-1.25,0.75],[t,0,4,0.02])$
plot        : plot2d ([discrete,makelist([p[1],p[3]],p,sol)], [xlabel,"t"],[ylabel,"y"])$
==========================================================
info : p[1] →t,p[2]  → x,p[3] → y

s1:[discrete,makelist([p[1],p[2]],p,sol)]$
s2:[discrete,makelist([p[1],p[3]],p,sol)]$
s3:[discrete,makelist([p[1],p[4]],p,sol)]$
s4:[discrete,makelist([p[1],p[5]],p,sol)]$

wxplot2d ([s1,s2,s3,s4], [xlabel,"t"],[ylabel,"xi,i=1..4"])$

2e) output :plot2d
==============

IC : [0.0998,0.3160,0.6536,0.9002];

m6a:[0.0998,0.3160,0.6536,0.9002]$

sol1:rk(m5b,m5c,m6a,domain)$

ss1:[discrete,makelist([pp[1],pp[2]],pp,sol1)]$
ss2:[discrete,makelist([pp[1],pp[3]],pp,sol1)]$
ss3:[discrete,makelist([pp[1],pp[4]],pp,sol1)]$
ss4:[discrete,makelist([pp[1],pp[5]],pp,sol1)]$

m6b:[ss1,ss2,ss3,ss4]$

wxplot2d (m6b, [xlabel,"time"],[ylabel,"xi,i=1..4"])$