kill(all)$

maxima:reaction kinetics A-> B,stirred tank , process,version 1.0
peter.vlasschaert@gmail.com,17/05/2023

load package
------------------

load(lrats)$

lumped system : C[A] =f(t) , t=time.
-----------------------------------------------
assumption : CSTR model for reactor.

balance : 
-----------

p1:Accum = In-Out+Generation;

model : startup "reactor"
--------------------------------

first order reaction :
===============
---------------------------------------------------
A → B  ,r[A] = -k*C[A]V, reaction rate 
input concentration :C[Ao] ,
---------------------------------------------------
set :C[A]=CA,C[Ao]=CA0,n[A]=nA

set :

pp1:[C[A]=CA,C[Ao]=CA0,n[A]=nA,F[in]=Fin];

pp2:[Out =0];

q1:C[A] = n[A]/V;/*mol/m^3*/

q1a:reverse(solve(q1,n[A])[1]);

model : startup "reactor"
--------------------------------

q2:flatten([Accum='diff(n[A],t) ,Generation = -k*C[A]*V,In =F[in]*C[Ain],pp2,q1a]);

q2a:lratsubst(q2,p1);

q2a1:subst(pp1,q2a);

@t=0 then C[Ain] =CA0

q2a11:part(args(q2a1)[1],1);

q2a22:subst([C[Ain]=CA0],q2a1);

integrate : 

q2b:'integrate(1/rhs(q2a22),q2a11,0,q2a11)='integrate(1,t,0,t);

assume(nA>0,Fin>0,k>0,CA0>0);

assume(CA0*Fin-k*nA>0);

q2b1:ev(q2b,nouns);

q2b2:expand(solve(q2b1,q2a11))[1];

q2b3:factor(q2b2);

@ startup : reactor volume = V is not constant 

model :unsteady after startup , V=VT = [m^3] → until steady.
--------------------------------------------------------------------------------
unsteady : depends on time 

definition : residence time 
-----------------------------------

declare(VT,constant)$

q3a:t[res] = VT/Fin;

before : t ≤ t[res], start up
after    : t[res]

use :q1a

q3a1:subst([V=VT],q1a);

use : pp1

q3a2:subst(pp1,q3a1);

q3b1a:[Accum = VT*'diff(CA,t)];

q3b1b:[In =Fin*CA0,Out=Fin*CA,Generation=-k*CA*VT];

q3b1c:flatten([q3b1a,q3b1b]);

model : unsteady 

use :p1

q3c:subst(q3b1c,p1);

assume(CA0>CA);

q3b:expand(ode2(q3c,CA,t));

→ need find %c

m1:sort(listofvars(rhs(q3b)));

condition allowed when : V=VT
-------------------------------------------

→ first step : integration , "q3b"

t =tbar "see definition : `res`idence time =t[res]"
CA=CAbar 

use : q3b

q3b1:subst([t=tbar,CA=CAbar],q3b);

→second step : solve for %C

q3b2:solve(q3b1,m1[1]);

→thrid step : subst "q3b2" into "q3b"

q3b3:subst(q3b2,q3b);

q3b33:ratsimp(q3b3);

→ fourth step: find "time constant = τ " system

assume(VT*k+Fin>0,VT>0);

q3b3a:limit(lhs(q3b3),t,inf) =limit(rhs(q3b3),t,inf);

q3b3a1:num(rhs(q3b3a))/Fin;
q3b3a2:expand(denom(rhs(q3b3a))/Fin);

q3bb:τ =args(q3b3a2)[1];

q3cc:reverse(q3bb);

→ fifth step: final solution 
t : tbar → t;  CA : CAbar  → CA

q3dd:lhs(q3b3) =expand(lratsubst([q3cc],rhs(q3b3)));

q33ee:collectterms(q3dd,CAbar);

→ sixth step : steady solution  , t→∞

use :q3b3a

q33ff:q3b3a;

use :q33ee

q33ff1:CA=limit(rhs(q33ee),t,inf);

What @ end startup period
======================
------------------------------------------------------------
condition : (time,concentration) =(tbar.CAbar)
------------------------------------------------------------

q44:tbar = VT/Fin;

use:q2b2

q44a:ratsimp(q2b2/VT);

use :nA/VT =CAbar

q44b:CAbar = expand(subst([t=tbar],rhs(q44a)));

→ use : definition tbar

q44c:subst([q44],q44b);

→ use : definition "time constant =τ"

use :q3bb

q44d:expand(lratsubst([q3cc],q44c));

--------------------------------------------------------------------------------------------------------------