kill(all)$

maxima:fvm:appendix:exponential scheme "1D -convection-diffusion,u > 0 ",version 1.0
peter.vlasschaert@gmail.com,15/02/2022

load : packages
--------------------

load(lrats)$

→ see website :
---------------------

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

ref 1:
maxima : convection -diffusion problem (exact : solution)integration,version 1.0
ref 2:
maxima : fvm: 1D, (Pe)clet number and graph (from exact equation) (Γ=cts),version 1.0

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

convection - diffusion problem (exponential scheme).
=======================================

1e) grid spacing : non-uniform.
------------------------------------------

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
----------------------------------------------------------------------

2e) model : 1D convection-diffusion
-----------------------------------------------

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

→ same result : when use integration ,with respect to x

p1a:args(lhs(p1)) =args(rhs(p1));

f :face  CV , CS =[w],[e]. →f :{w,e}
----------------------------------------------

p1b1:J[f] =(lhs(p1a) - rhs(p1a))[1];

→ without : sort  ( preference)

p1b2:flatten([args(lhs(p1a)[1]),args(rhs(p1a)[1])]);

p1b22:delete(Φ,p1b2);

→ with    :sort

p1b3:sort(p1b2);/*change order in the list */

use : p1b2

p1b2a:[u(f),ρ(f),Γ(f),at(diff(Φ(x),x),x=f)];

p1b2b:map("=",p1b22,p1b2a);

 J[f] =total fux = convection flux - diffusion flux
---------------------------------------------------------------

p1c:lhs(p1b1) =psubst(p1b2b,rhs(p1b1));

see : figure 1,part of convection term.

p1c1:F[f] = ratcoeff(part(rhs(p1c),1),Φ);

integration over CV : CS = [w] →CS = [e]

p1c2:lratsubst([reverse(p1c1)],p1c);

p1c2a:lhs(p1c2)=subst([Φ=Φ(f)],part(rhs(p1c2),1))+part(rhs(p1c2),2);

→ final result : integration

p1c2b1:subst([f=e],lhs(p1c2a))-subst([f=w],lhs(p1c2a)) =0;
p1c2b1:subst([f=e],rhs(p1c2a))-subst([f=w],rhs(p1c2a)) =0;

exact : solution 1D -convection-diffusion (x-direction)
-----------------------------------------------------------------------

see : ref 2.use "p1"

p2:(Φ(x)-Φ(0))/(Φ(L)-Φ(0))=(%e^((Pe(L)*x)/L)-1)/(%e^Pe(L)-1);

domain : cell peclet number , 0≤ x≤ L in "p2"

a) domain : between  p≤e≤E ,see: Figure 1
------------------------------------------------------
p : reference → Φ(0) in p2 :Φ(p)
E :                 → Φ(L) in p2: Φ(E)
L :Δx(e)  → rhs

lhs → 

p2a1:[Φ(0)=Φ(p),Φ(L)=Φ(E)];

rhs →

p2a2:[L=Δx(e)];

result : profile assumption,interval

p2e:subst(p2a1,lhs(p2)) =subst(p2a2,rhs(p2));

p2e1:denom(lhs(p2e));

p2e2:p2e*(p2e1)-part(p2e1,2);

b) domain : between  W≤w≤p ,see: Figure 1
------------------------------------------------------
W : reference → Φ(0) in p2 :Φ(W)
p :                   → Φ(L) in p2: Φ(p)
L :Δx(w)  → rhs

lhs → 

p2b1:[Φ(0)=Φ(W),Φ(L)=Φ(p)];

rhs →

p2b2:[L=Δx(w)];

result : profile assumption,interval

p2w:subst(p2b1,lhs(p2)) =subst(p2b2,rhs(p2));

p2w1:denom(lhs(p2w));

p2w2:p2w*(p2w1)-part(p2w1,2);

first computation : J[e] ,f =e , into p1c2a, "see above"
--------------------------------------------------------------------
x = xe

p3a:p1c2a;

p3a1:subst([f=e],p3a);

step 1: differentation ,x = xe
-------------------------------------

use :p2e2,p1b2a

p4a1:subst([f=e],last(p1b2a)) =subst(x=xe,diff(rhs(p2e2),x));

step 2: Φ(e)
---------------

use :p1c2a

p4a2:ratcoeff(part(rhs(p1c2a),1),F[f]) $

p4a3:subst([f=e],p4a2) =subst([x=xe],rhs(p2e2));

step 3: convection @  interface [e]
-----------------------------------------

p5a: Jcon[e] = subst([f=e],part(rhs(p1c2a),1));

p5a1:subst([p4a3],p5a);

step 4: (-1)*diffusion @  interface [e]
-----------------------------------------------
use :p3a1,p4a1

p5b:Jdiff[e] = subst([f=e],part(rhs(p1c2a),2));

p5b1:subst([p4a1],p5b);

→ need : some simplification  "p5b1"

definition : Pe(.)
----------------------
f : interface of CV  : CS =[e]

definition : Diffusion (relation to peclet Cell Number)
=======================================

remember : 1D CV , f =e , here

p6a:D[f] =Γ(f)/Δx(f);

p6a1:reverse(p6a);

definition : Pe(.), Cell peclet number.
===========================

p6b:Pe(m1) = F[m]/D[m];

→ solve : p6b

p6b1a:solve(p6b,F[m])[1];
p6b1b:solve(p6b,D[m])[1];

→ simplification :p5b1

a) first use : f =e , p6a1
------------------------------

p7a:subst([f=e],p6a1);

into : p5b1

p7a1:lratsubst([p7a],p5b1);

b) need : F[e] inside  p6b1a,m1=Δx(e),m =e
---------------------------------------------------------

p7a3:subst([m=e,m1=Δx(e)],p6b1a);

p7a4:[reverse(p7a3)];

→ into : p7a1

p7a5:lratsubst(p7a4,p7a1);

final result: J[e]

use : {p5a1,p7a5}

p8a:J[e] =Jcon[e]+Jdiff[e];

p8b:subst([p5a1,p7a5],p8a);

p8b1:ratsimp((p8b));

p8b2:facsum(p8b1);

analog : J[w]
------------------
see : e =w, 

b) domain : between  W≤w≤p ,see: Figure 1
------------------------------------------------------
W : reference → Φ(0) in p2 :Φ(W)
p :                   → Φ(L) in p2: Φ(p)
L :Δx(w)  → rhs
( come from : above)

p8c1:subst([e=w],p8b2)$

p8c2:psubst([p=W,E=p],p8c1);

rem :psubst ,list p on lhs & rhs

general equation :
------------------------
a[p]*Φ(p) = a[W]*Φ(W) +a[E]*Φ(E)

inflow = outföw , from CV
------------------------------------
w → in CV
e → out CV

p9: rhs(p8b2) =rhs(p8c2);

p9a:expand(p9);

p9a1:ratcoeff(p9a,Φ(p));
p9b1:ratcoeff(p9a,Φ(W));
p9c1:ratcoeff(p9a,Φ(E));

result : general equation
---------------------------------

p10:a[p]*Φ(p) = a[W]*Φ(W) +a[E]*Φ(E);

p10a1:a[p] =lhs(p9a1)-rhs(p9a1);

p10b1:a[W] =rhs(p9b1)-lhs(p9b1);

p10b1:a[E] =rhs(p9c1)-lhs(p9c1);

rem : very good scheme ,computational use much more memory of the computer.
          for CFD.

CFD = (C)omputational (F)luid (D)ynamics.

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