kill(all)$

maxima : fvm :convection -diffusion relation peclet cell number (max :limit),version 1.0
peter.vlasschaert@gmail.com,12/11/2022

load :packages , option
------------------------------

load(lrats)$ ratprint:false$

convection - diffusion problem ( integration),exact
-------------------------------------------------------------------

ref :updintegrate (option xy = 1,2)
---------------------------------------------
new:10/11/2022:maxima : updintegrate for differential inside the differential ,version 1.0a

grid spacing : different
=================
same:
new:31/10/2022:maxima :  appendix : fvm : Grid generation for CV in 1D,version 1.0

rem : method introduced by Patankar & Spalding

x-dirction : graph

legend :
------------------------------------------------------------------
||.|| = distance
||pe|| = ||pw|| = Δx
||a1a2|| = ||Wp|| =Δx(w)
||a2a3|| = ||Ep|| =  Δx(e)
=====================================
rem :||aiaj|| = ||ajai|| ,{i,j} = integer (≠ 0)
------------------------------------------------------------------
nodal points : {W,p,E}
interface : {w,e}
-------------------------------------------------------------------
orange : CV = Control Volume ,CS = Control Surface
→ in : w  (CS)
→ out : e (CS)
CV : w ≤ x ≤ e
--------------------------------------------------------------------
pos : direction :
in   :F(w) = u(w)*ρ(w) → Φ (w)
out :F(e)  = u(e)*ρ(e)   → Φ (e)
---------------------------------------------------------------------
interface :{Φ (w),Φ (w)}
nodal      :{Φ (W),Φ (p),Φ (E)}:know
----------------------------------------------------------------------

1e) convection -diffusion model.
------------------------------------------

p1:'diff(u*ρ*Φ,x) = 'diff(Γ*'diff(Φ,x),x);

ρ = density
u = velocity x-direction 
Φ = function ( see below ,Φ(x))
Γ = diffusion coef.

conservation law :
----------------------

p1a:'diff(u*ρ,x) = 0;

explain :conservation law
------------------------------------
====================================================
see figure : F(.) =u(.)*ρ(.) , . ={w,e}
---------------------------------------------
CV :  orange Cell
CS  :  in CS (= F(w)) = out CS(=F(e))
====================================================
---------------------------------------------------------------------------------------------
general : ρ = constant , div(ρ*v) = 0 → div(v) = 0 for 1D : v=u
---------------------------------------------------------------------------------------------
div = divergence operator , see : vector analysis

2e) integration of model
--------------------------------

a) definition : updintegrate
-----------------------------------

updintegrate(expr,il,ol,co,xy):= block([zz],  
if  xy =1  then 
 (
        t1aa:setify(co) ,
    t1b:args(expr),t1c:args(t1b[1]),
    t1d:t1b[2],t1e:setify(t1c),
    t1f:listify(setdifference(t1e,t1aa)),
    t1g:makelist(atom(t1f[i]),i,1,length(t1f)),
    t1g1:makelist( if t1g[i] = false then nx:i ,i,1,length(t1g)),
    t1h1:delete(t1f[nx],t1f),
     /* no differential part inside differential */   
    t1h2:makelist(subst([ttt=t1h1[i]],ttt(t1d)),i,1,length(t1h1)),
    t1h3:product(t1h2[i],i,1,length(t1h2)),
    t1h4:subst([x=zz],t1h3),    
     /* differential part inside differential      */
    t1i1:t1f[nx] ,
    t1i2:args(t1i1),
    t1i3:t1i2[1](t1i2[2]),   
    t1i4: at(subst([t1i2[1]=t1i3],t1i1),[x=zz]),
    /* combine : inside */
    t1j1:t1h4*t1i4,
    /* constant.list product */
    t1k0:product(co[i],i,1,length(co)),    
    /* integration interval : [il,ol] */    
    t1k1:at(t1j1,[zz=ol]) -at(t1j1,[zz=il]),
    t1k2:t1k0*t1k1    
    )
    else
    (xc1:args(expr),
     xc2:args(xc1[1]),
     xc3:product(xc2[i](xc1[2]),i,1,length(xc2)),
     xc4:at(xc3,[xc1[2]=zz]),
     xc5:at(xc4,[zz=ol]) - at(xc4,[zz=il])   
    )  
)$

b) integration of model 
------------------------------

px1:[F(e)=u(e)*ρ(e),F(w)=u(w)*ρ(w)];
px1a:makelist(rhs(px1[i]) = lhs(px1[i]),i,1,length(px1));

→ convection + diffusion

p2:updintegrate(lhs(p1),e,w,[1],2) =updintegrate(rhs(p1),e,w,[1],1);

p2x:lratsubst(px1a,p2);

rem :[1] no constant in the differential.

→ conservation law.

p2a:updintegrate(lhs(p1a),e,w,[1],2) =rhs(p1a);

p2ax:subst(px1a,p2a);

3e) need profile assuption (discritization scheme)
---------------------------------
use : cds = central difference scheme

→ advection : part

p3:[Φ[w]=(Φ[p]+Φ[W])/2,Φ[e]=(Φ[p]+Φ[E])/2];

p3h:[Φ(w)=Φ[w],Φ(p)=Φ[p],Φ(e)=Φ[e],F(w)=F[w],F(e)=F[e]];

p3h1:lratsubst(p3h,lhs(p2x));

p3h2:expand(lratsubst(p3,p3h1));

use :p2x ,lhs

→ diffusion : part

p3a:[part(rhs(p2x),1,2) =1/Δx(w)*(Φ[p]-Φ[W]),part(rhs(p2x),2,1,2) =1/Δx(e)*(Φ[E]-Φ[p])];

p3ax1:expand(lratsubst(p3a,rhs(p2x)));

use:p2x ,rhs

→ combine

p3xx:p3h2 =p3ax1;

p3xx1:p3h2-p3ax1;

rem : no source (S=0)

→ conservation law (see : above)

use :p3h

p3yy:lratsubst(p3h,p2ax);

4e) general rule
--------------------
see :new:10/11/2022:maxima :  diffusion+source problem,a[p],a[E],a[W],b,version 1.0

find :{a[p],a[W],a[E],}

p4:a[p]*Φ[p] = a[W]*Φ[W] + a[E]*Φ[E] ;

use :p3xx1

find : coef. by differentation 
------------------------------------
neg sign ,a[p],lhs  "see:p4"

p4s1:a[p]=-1*diff(p3xx1,Φ[p]);
p4s2:a[W]=diff(p3xx1,Φ[W]);
p4s3:a[E]=diff(p3xx1,Φ[E]);
p4s :[p4s1,p4s2,p4s3];

set : D[.] =Γ(.)/Δx(.) , "Diffusion coef.",{e,w}

p4s4:[D[e] = Γ(e)/Δx(e),D[w] = Γ(w)/Δx(w)];
p4s4a:makelist(rhs(p4s4[i])=lhs(p4s4[i]),i,1,length(p4s4));

p4sa:expand(lratsubst(p4s4a,p4s));

p4sa1:makelist(rhs(p4sa[i])=lhs(p4sa[i]),i,1,length(p4sa));

find : real solution for problem : "same sign" for all :  {a[p],a[W],a[E],}
--------------------------------------------------------------------------------------------

1e) step : check a[W]
----------------------------

p4b1:p4sa[2];

2e) step : F[w] > 0,D[w] > 0
-------------------------------------

assume(F[w]>0,D[w]>0);
p4b2:lhs(p4b1) = sign(rhs(p4b1));

3e) step because a[W] > 0 → follows : a[E] > 0
--------------------------------------------------------------

p4b3: lhs(p4sa[3]) = sign(1);

4e) step : P(e) = F[e]/D[e], P(w) =F(w)/D[w]
------------------------------------------------------------
Cell Peclet number = advection/diffusion

p4b41:[F[e]/D[e] =P(e)];
p4b42:lratsubst(p4b41,p4sa[3]);

p4b43:num(rhs(p4b42)) =0;

rem :D[e] ≠ 0 ,can never be zero.

p4b44:part(lhs(p4b43),1,2) =0;

max : limit (constraint Peclet Cell Number)

p4b45:solve(p4b44,listofvars(p4b44))[1];

rem : P(e) ≤ 2;Peclet Cell Number : P(e).

find : relation between a[p] = function(a[W] ,a[E]) and use : conservation law
------------------------------------------------------------------------------------------------------
see above :p3yy

1e) step : conservation law + a[p], "see:p4sa"
----------------------------------------------------------

→ conservation law

p4c1a:[p3yy];

→ a[p]

p4c1b:p4sa[1];

2e) step : a[p] = a[W]+a[E] +G
----------------------------------------
find :G

p4c2: a[p] =expand(lratsubst(p4sa,a[W]+a[E]+G));

p4c2a:[reverse(p4c2-p4c1b)];

use : conservation law "p4c1a"

p4c2b:lratsubst(p4c1a,p4c2a);

p4c3:makelist(p4sa1[i],i,2,length(p4sa));

p4c3a:lratsubst(p4c3,p4c2);

use : p4c2b

p4c3b:subst(p4c2b,p4c3a);

3e) sstep: sign of a[p] =pos
------------------------------------
use :p4b2,p4b3, "see:above"

p4d:[p4b2,p4b3];
p4d1:subst(p4d,p4c3b);

rem : a[p] > 0
result : "same sign" for all :  {a[p],a[W],a[E],}
---------------------------------------------------------------------------------
rem : coef different we get : physical solution is bad.
(math :diagonal dominant matrix ,domain solution,let to cvg.)
---------------------------------------------------------------------------------

explain : flow direction is taken into account.
------------------------------------------------------------

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

REM :{a[p],a[W],a[E]},"same sign" inconsitent values in interval [w] & [e]
-------------------------------------------------------------------------------------------------
reason : Peclet Cell Number → see above. (constraint)
-------------------------------------------------------------------------------------------------
REM : not general peclet number for the overall domain.
-------------------------------------------------------------------------------------------------

REM :values :know @ nodal points
REM : with this scheme : direction of flow 
See : other scheme  
1e) - upwind scheme
----------------------------
new:04/11/2022:maxima :  appendix : fvm : upwind-scheme Φ[w],version 1.0 when(u>0,u<0),version 1.0
new:04/11/2022:maxima :  appendix : fvm : upwind-scheme Φ[e],version 1.0 when(u>0,u<0),version 1.0
2e) - quick     scheme
----------------------------
new:04/11/2022:maxima :  appendix : fvm : quickscheme Φ[e],version 1.0 when(u>0,u<0),version 1.0
new:04/11/2022:maxima :  appendix : fvm : quickscheme Φ[w],version 1.0 when(u>0,u<0),version 1.0

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

5e) interface values :CV with " input = [w] & output = [e]"
----------------------------------------------------------------------------

shown : exact solution
------------------------------

interface :[e]
-----------------

→ Φ(e)

(Φ(e)-Φ(p))/(Φ(E)-Φ(p))=(%e^(u·Δx·ρ/Γ)-1)/(%e^(Δx(e)·u·ρ/Γ)-1);

interface :[w]
------------------

(Φ(w)-Φ(W))/(Φ(p)-Φ(W))=(%e^(u·Δx·ρ/Γ)-1)/(%e^(Δx(e)·u·ρ/Γ)-1);

rem : when use :
new:31/10/2022:maxima :  appendix : fvm : Grid generation for CV in 1D,version 1.0
new:18/10/2022:maxima :  convection -diffusion problem (exact : solution),version 1.0
---------------------------------------------------------------------------------------------------------------
Φ(e),Φ(w)
=======
when grid uniform         : into p9b
x                                       → Δx
L                                       → Δx/2

here :

Φ(e),Φ(w)
=======
when (non) uniform grid  : into p9b
x                                       → Δx
L                                       → Δx(e),Δx(w)

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