maxima,kinetics : -('diff(C,x,1))*u-C*k+('diff(C,x,2))*D=0,symbolic,part 1,version 1
peter.vlasschaert@gmail.com,27/01/2023

Model :Tubular reactor and differential section
-----------------------------------------------------------
D = diffusion coef.            ,[m^2/s]
u  = velocity x-direction   ,[m/s]
C = concentration             ,[mol/m^3]
x  = dfistance                     , [m]

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

show this below : m1

m1:-('diff(C,x,1))*u-C*k+('diff(C,x,2))*D=0;

boundary conditions : steay problem (=independent of time)

@ x=0, inlet

m2a:u*Cin = u*C -D*'diff(C,x);

@ x=L,outlet   : L = length of reactor.

m2b:'diff(C,x) =0;

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

mass balance : Control Volume @x
--------------------------------------------
CV = differential section (A*dx)
A = Area , A*dx=dV

kill(all)$

→ without chemical reaction :
---------------------------------------

p1: IN = OUT +ACCUM;
p1a:solve(p1,ACCUM);

rem :entering differential section [A] or leaving differential [A] =[mol/s]

IN : convection & dispersion [mol/s]

→ convection : (u velocity in x-direction:[m/s])

IN1:u*A*C; /*[m/s*m^2*mol/m^3] =[mol/s]*/

→ dispersion : (D = diffusion coef. [m^2/s])

rem : IN2: D*C*dx =[m^2/s*mol/m^3*m] =[mol/s],wrong why , C=C(x) 

IN2:-D*'diff(C,x)*A;/*[m^2/s*mol/m^4*m^2=mol/s]*/

rem : C depends on x, fick law : minus in front : -D*diff(C,x)
---------------------------------------------------------------------------
suppose : A,D are constant

OUT : convection & dispersion [mol/s] , displacement x to x+Δx, leave : x+Δx

rem : Taylor expansion :  leaving CV
------------------------------

s1:args(taylor(f(x), x, x0, 3));

s2:subst([x-x0=Δx],s1);

s3:subst([x0=x],s2);

approximation f(x+Δx) : Δx displacement

s4:f(x+Δx) = sum(s3[i],i,1,2);

s4a: subst([ratcoeff(part(rhs(s4),1),Δx) =diff(f(x),x)],s4);

reason : f : f(x) → C:C(x)

s4b: subst([f=C],s4a);

use : maxima depends

depends(C,x);

s4c:lhs(s4b) =subst([C(x)=C],rhs(s4b));

→ convection : (u velocity in x-direction:[m/s]) , @ x+Δx

OUT1:subst([C=lhs(s4c)],IN1);

OUT1:subst([s4c],OUT1);

→ dispersion : (D = diffusion coef. [m^2/s]) ,      @ x+Δx

OUT2:subst([C=lhs(s4c)],IN2);

OUT2:subst([s4c],OUT2);

OUT2:ev(OUT2,nouns);

ACCUM : accumulation [mol/s], unsteady system 

ACCUM : 'diff(C,t)*A*Δx;

mass balance : reactor 
---------------------------

v1a:INx    = IN1+IN2;
v1b:OUTx =OUT1+OUT2;
v1c: ACCUM = rhs(v1a) -rhs(v1b);

v1c:expand(v1c);

unsteady : time depent , steady : no time depent

v1c1:ratcoeff(v1c,A*Δx);

unsteady : C =C(x,t)  -> v1c1 , C =concentration [mol/s]

v1c2:subst([lhs(v1c1)=0],v1c1);

v1c2:reverse(v1c2);

steady : C =C(x)  -> v1c2 , C =concentration [mol/s]

→ with reaction : first order

Simple reaction : AA -> BB
----------------------------------
AA react : lhs : v1c, -k*C = -k*CAA, minus sign 

second order differential :constant coef , see assume : below

v1c3:lhs(v1c2) -k*C = rhs(v1c2);

assume(k>0,D>0,u>0);

v1c4:ode2(v1c3,C,x);

method;

dispersion model changes the reactor analysis from an initial value problem to
a boundary value problem. The boundary conditions proposed by Danckwerts are :

dimension : reactor  , 0<x<L, inlet < x < outlet
----------------------------------------------------------

→  Boundary Condition :

 begin : reactor , @ x =0
------------------------------
→ inlet :

bc1: u*Cin = u*C -D*diff(C,x);

 end : reactor , @ x =L
------------------------------
→ outlet :

bc2: diff(C,x) =0;

what about : implement BC, v1c4
-----------------------------------------

@ x=0,bc1, steady

→ first step : differentation:v1c4 with respect to x

w1:diff(v1c4,x);

w1a:w1$

@  x=0

w1:at(w1,x=0);

→ second step : w1*(-D)

w1:-D*w1;

→ thrid step : both sides +C*u

w1:w1+C*u;

→ fourth step : rhs(w1)= u*Cin

w1:rhs(w1) = lhs(bc1);

@ x=L,bc2, steady

→ first step : use w1a

w2:w1a;

@  x=L

w2:at(w2,x=L);

→ second : use bc2

w2:rhs(w2) = rhs(bc2);

find : {%k1,%k2}

aa:flatten(solve([w1,w2],[%k1,%k2]))$

find → %k1

a1:aa[1];

find → %k2

a2:aa[2];

solution found by ode2 ,maxima :v1c4

sol:v1c4;

→ first step : input{%k1,%k2} found from {bc1,bc2},subst : aa, "see above"

sol:subst(aa,sol);

→ solve C :sol

solC:solve(sol,C)[1];

solCA:expand(ratsimp(rhs(solC)));

rem : u>0,D>0,k>0 ,then this solution

error : lisp , ratsimp en then expand.

make function : concentration(parameters)
------------------------------------------------------

Cf(Cin,k,L,u,D):=''solCA;

given : Reaction A -> B , k = reaction constant ,Cin = concentration A
------------------------------------------------------------------------------------------
u = 1           [m/s]          , velocity
D = 0.001   [m^2/s]      , diffusion constant
k = 1           [1/s]            , reaction constant
L = 10         [m]              , length of the reactor
Cin = 1       [mol/m^3] , begin concentration CA
------------------------------------------------------------------------------------------
use : Cf ,"see above"

example : fuction given parameters

sp1:[u=1.0,D=0.001,k=1.0,L=10,Cin=1];
sp1a:CA = subst(sp1,Cf(Cin,k,L,u,D));

sp1b:rhs(sp1a);

plot: error :arithmetic error FLOATING-POINT-OVERFLOW signalled
wxplot2d([sp1b], [x,0,10])$

other methods : to solve by ,finite difference,part 2, how to do
---------------------------------------------------------------------------------------------------------------------------------------------------