kill(all)$

maxima : appendix:kinetics: convection-diffusion ,-u*dC/dx+D*d/dx(dC/dx) = 0,version 1.0
peter.vlasschaert@gmail.com,29/02/2023

load packages
------------------

load(lrats)$

Mathematical model : without chemical reaction.
--------------------------------------------------------------

unsteady problem : time depent

p1:-u*'diff(C,x) + D*'diff(C,x,2) ='diff(C,t);

 D = diffusion coef.           ,[m^2/s]
u  = velocity x-direction   ,[m/s]
C = concentration              ,[mol/m^3]
x  = distance                       , [m]

rem : time depent term 

p1a:rhs(p1);

steady problem : time indepent

p1b:lhs(p1) = 0;

steady problem : diffusion + convection = 0

p2:args(lhs(p1b));

→ find : diffusion term , contain D

p2a:p2[1];

→ find : convection term , contain u

p2b:part(p2[2],1);

rem : to find minus sign :

p2b1:part(p2[2],0);

rem : How to analyse parts of ode (ordinary differential equations)

upwind scheme 
--------------------

→ explain : relation between u & x

print("upwind : depent whats happening before")$

→ differential :

dC/dx :

first order scheme : finite difference scheme : "ref : C[i],x=0"

u,x same direction : suppose (u>0,x>0)  , C[i] > C[i-1]

p3a:'diff(C,x) = (C[i] -C[i-1])/Δx;/*backward difference*/

u,x different direction : suppose (u<0,x>0)  , C[i+1] > C[i]

p3b:'diff(C,x) = (C[i+1] -C[i])/Δx;/*forward difference*/

rem : x>0, positive direction cartesian system.

second order scheme : finite difference scheme : "ref : C[i] , x=0"

→ first step : need quadratic equation
C depent on x

p4:C = a+b*x + c*x^2;

p4a:sort(listofvars(rhs(p4)));

p4a1:makelist(p4a[i],i,1,3);
p4a2:last(p4a);

p4b:'diff(C,x)=diff(rhs(p4),p4a2);

→ u>0  , point upwind : i-1,i-2 , before i

point:i,i-1,i-2

p4b1:C[i] =subst([x=0],rhs(p4));
p4b2:C[i-1] =subst([x=-Δx],rhs(p4));
p4b3:C[i-2] =subst([x=-2*Δx],rhs(p4));

find: a,b,c ; for p4a :

sol1:flatten(solve([p4b1,p4b2,p4b3],p4a1));

final result : u>0 ,use p4b

df1:subst(sol1,p4b);

ref :C[i],at x=0

df1:lhs(df1) = limit(rhs(df1),p4a2,0);

→ u<0  , point upwind : i+1,i+2 , before i

point:i,i+1,i+2

p5b1:C[i] =subst([x=0],rhs(p4));
p5b2:C[i+1] =subst([x=Δx],rhs(p4));
p5b3:C[i+2] =subst([x=2*Δx],rhs(p4));

find: a,b,c ; for p4a :

sol2:flatten(solve([p5b1,p5b2,p5b3],p4a1));

final result : u<0 ,use p4b

df2:subst(sol2,p4b);

ref :C[i],at x=0

df2:lhs(df2) = limit(rhs(df2),p4a2,0);

upwind scheme :fvm see finite volume method

How to solve : p1b
-----------------------

p6:p1b;

assume(u>0,D>0);

use : ode solve symbolic ->

p6a:ode2(p6,C,x);

p6b:subst([C=C(x)],p6a);

what about : Peclet number : Pe ( no dimension)

suppose : boundary condition 
-------------------------------------
begin : x=0 → bc1
end    : x=L → bc2

@ x=0 ,bc1

p6b1:subst([x=0],p6b);

@ x=L ,bc2

p6b2:subst([x=L],p6b);

solve :{p6b1,p6b2} → {%k1,%k2}

p6b3:listofvars(rhs(p6b)); 

p6b3a:makelist(p6b3[i],i,1,2);

p6b4:flatten(solve([p6b1,p6b2],p6b3a));

use :p6b

sol3:factorsum(subst(p6b4,p6b));

definition : Pe number
----------------------------

explain  : D [m^2/s] & L*u [m*m/s]

p7:Pe = convterm/diffterm;

p7a:[convterm=L*u,diffterm=D];

p7b:subst(p7a,p7);

p7b1:reverse(p7b);

exact : solution for p1b when u >0,D>0
============================
u from assume above.

solution : sol3 use Pe

sol3:subst([p7b1],sol3);

rem : finite difference scheme both : central difference
-------------------------------------------------------------------
for first order ,second order differential

p8:-u*(C[i+1]-C[i-1])/(2*Δx)+D*(C[i+1]+C[i-1]-2*C[i])/Δx^2 =0;

first step : *Δx^2/D

p8a:expand(p8*Δx^2/D);

second step : cell peclet number = u*Δx/D =Pe(n)

p8b:lratsubst([(u*Δx)/D=Pe(n)],p8a);

p8c:ratsimp(p8b);

this system : result in tridiagonal system : solved with TDMA,
(also know Thomas Algorithm)

coef all  pos : Pe(n) = 2 , limit number ,explain ,detail ,see : fvm ,website.

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