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

ref1:maxima,kinetics :
-('diff(C,x,1))*u-C*k+('diff(C,x,2))*D=0,symbolic,part 1,version 1

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;

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

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;

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

error : lisp , ratsimp en then expand.

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

general equation :finite differenc method
use :m1

→ first step : find differentials :

p1:m1;

use :args

p1a:args(lhs(p1));

find :'diff(C,x,1)

p1a1:args(args(p1a[1])[1])[1];

find :'diff(C,x,2)

p1b1:args(args(p1a[3]))[1];

→ second step: finite differenc approx. for differential:

finite difference : forward scheme, 'diff(C,x,1):
---------------------------------------------------------
differential section : along x-axes : Δx

p1c1:p1a1 = (C[i] -C[i-1])/Δx;

finite difference : central difference scheme, 'diff(C,x,2):
--------------------------------------------------------------------
differential section : along x-axes : Δx

p1c2:p1b1 = (C[i+1]-2*C[i] +C[i-1])/Δx^2;

p1c:[p1c1,p1c2];

→ third step: finite differenc of differential equation ( = ode):

use : p1

p1d:subst(p1c,p1);

p1da:lhs(p1d);

p1db1:part(p1da,3);

p1db2:subst([C=C[i]],p1db1);

p1db:[p1db1=p1db2];

p1e:subst(p1db,p1da)=0;

→ build : tridiagonal matrix
                  =============

all nodes : inlet (node),general node ,outlet (node)

n:7$

rem : real simulation take n =50

general node :

find : equation internal nodes : i, " node inside reactor: measure "
--------------------------------------------------------------------------------

i: 2 to i:n-1

eq(i):=''p1e;

example : general node (=2)

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

eq(2);

find coef. before:C[i], i=1,2,3

p1a:factor(expand(makelist(ratcoef(lhs(eq(2)),C[i])*C[i],i,1,3)));

p1b:sum(p1a[i],i,1,3) =0;

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

p2:makelist(eq(i),i,2,n);

begin (=node):@ x=0
------------------
use : bc1

p2a:bc1;

finite difference : central difference scheme, 'diff(C,x,1):
---------------------------------------------------------------------
result : p2a1

p2a1:(C[i+1]-C[i-1])/(2*Δx)$
p2a2:lhs(p2a) = C[i]*u-D*p2a1;

first node : i=1, begin node.

→ i:1

p2b:subst([i=1],p2a2);

→ solve :C[0]

p2b1:solve(p2b,C[0]);

use : eq[i], general 

p2b2:eq(1);

→ final equation : node(1) ,i=1

use : eliminate : C[0] from p2b1 & p2b2

p2c:expand(subst(p2b1,p2b2));

end (=node):@ x=L
------------------
use : bc2,central difference.

p3a:reverse(p2a1 =0);

→ i:n

p3a1:subst([i=n],p3a);

→ solve C[n+1]

p3a2x:reverse(solve(p3a1,C[n+1]));

→ use eq[i] , i=n

p3a3:eq(n);

p3a4:subst([C[n-1]=C[n+1]],p3a3);

find variable list: i = 1..n
-------------------------------

p4a:makelist(C[i],i,1,n+1);

all equation : tridiagonal system.
------------------------------------------

p4b:flatten([p2c,p2,p3a4]);

H1*C = H2 , (A*x=b)
---------------------------

p5a:coefmatrix(p4b,p4a);

given : Reaction A -> B , k = reaction constant ,Cin = concentration A
------------------------------------------------------------------------------------------
u = 1           [m/s]          , velocity
D = 0.0001  [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"

use ref1:sp1

→ example :

p5b:[u=1.0,D=0.0001,k=1.0,L=10,Cin=1];

p5c:float(subst(p5b,p5a));

use : Δx =L/n , differential section.

p5c1:Δx =float(subst(p5b,L/n));

find:

→ H1, lhs

p5c2:subst([p5c1],p5c);

→ H2, rhs , Cin coef, first equation : p4b[1]

p6:p4b[1];

→ coef . Cin

lhs :

p6a:ratcoeff(lhs(p6),Cin)*Cin;

→ lhs to rhs :

p6b:expand((p6-p6a));

→ H2, rhs , Cin coef, first equation : p4b[1]

p6c:flatten([rhs(p6b),makelist(0,i,2,n+1)]);

use :p5b

p6c1:subst(p5b,transpose(p6c));

use :p5c1

p6c2:subst(p5c1,p6c1);

→ H2, rhs , Cin coef, from first equation : p4b[1]

find:C[i] ,i=1..n

solf:solve(p4b,p4a)$

example : see ref1,see : "above"

solfn:flatten(subst(p5b,%))$
solfn:subst(p5c1,solfn);

table : output : Concentration [A] ,C[i=1..n+1]
----------------------------------------------------------

sd :table_form(%)$

make plot : displacement begin to end reactor
---------------------------------------------------------

→ step 1 :x-values

x-values : distance : begin to end : a0=0→ L
-------------------------------------------------------
use :p5c1

xv:subst([p5c1,a0=0],makelist(a0+(i-1)*Δx,i,1,n+1));

nq:length(xv);

→ step 2 :y-values

y-values : concentration A : CA ,here C[i], i=1,2,....,n+1
--------------------------------------------------------------------
use: solfn 

yv: makelist(rhs(solfn[i]),i,1,n+1);

nq1:length(yv);

→ step 3 :discrete plot of point

How to this : Menu : Plot →Special →Discrete plot ,correct [xv]⟶xv,[yv]
rem : 10 is here length reactor

wxplot2d([['discrete, xv, yv]], [x,0,10])$

When : CA -> 0 , reaction : A → B , B  will reach value one by full conversion:
result : CB = 1

info : how to do this xv ,yv
----------------------------------
(%i1) xy: [[10, .6], [20, .9], [30, 1.1], [40, 1.3], [50, 1.4]]$
(%i2) plot2d([[discrete, xy], 2*%pi*sqrt(l/980)], [l,0,50],
        [style, points, lines], [color, red, blue],
        [point_type, asterisk],
        [legend, "experiment", "theory"],
        [xlabel, "pendulum's length (cm)"],
        [ylabel, "period (s)"])$

rem : wx to makeplot inline document

wxplot2d([[discrete, xv,yv], [discrete, xv,yv]],
        [style, points, lines], [color, red, blue],
        [point_type, asterisk],
        [title, "reaction: A→B"],
        [legend, "experiment", "theory"],
        [xlabel, "distance along the reactor"],
        [ylabel, "CA ,concentration of A"])$