kill(all)$

maxima : fvm : 2D convection-diffusion (start :discrete schema fvm), part 2, version 1.0
peter.vlasschaert@gmail.com.14/12/2022

rem : Model 2D :without source , S=0, steady → discrete 2D scheme.

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

rem : FVM = (F)inite (V)olume      (M)ethod  (use : conservation law ,use for CFD)
rem : FDM = (F)inite (D)ifference  (M)ethod
rem : FEM  = (F)inite (E)lement     (M)ethod
CFD = (C)omputational (F)luid (D)ynamics

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

explain →  conservation law , "see below"

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

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

load(lrats)$

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

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

ref:1
maxima : fvm :convection -diffusion relation peclet cell number (max :limit),version 1.0
ref:2, use p3b
maxima : fvm : 2D convection-diffusion (start :1d transporteq), part 1, version 1.0
ref : 3,use for 2D scheme
maxima :diffusion+source problem,a[p],a[E]a[W],b,version 1.0 

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

1e) 2D :diffision-convection equation,steady,S=0
-----------------------------------------------------------------
u : velocity  → x-direction
v : velocity   → y-direction
2D "dimension =" :Φ = Φ(x,y)
----------------------------------------------------------------

model : 2D
---------------
use :p3b

p1:'diff((v*Φ*ρ),y,1)+'diff((u*Φ*ρ),x,1)=Γ*('diff(Φ,y,2)+'diff(Φ,x,2));

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

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

scheme 2D :FVM
-----------------------
CV : (c)ontrol (v)olume
CS : (c)ontrol  (s)urface

legend :
======

--------------------------------------------------------------------------------
W :west
E  : east
S  : south 
N . north
--------------------------------------------------------------------------------
u: x-direction. WW→W → w.p,e→ E→  EE, pos x-direction
v: y-direction.  SS→  S →    s.p,n → N→NN, pos y-direction
--------------------------------------------------------------------------------
p = p(i,j) →  Φ:Φ(i,j)
--------------------------------------------------------------------------------
W :Φ:Φ(i-1,j)
E  :Φ:Φ(i+1,j)
S  :Φ:Φ(i,j-1)
N :Φ:Φ(i,j+1)
--------------------------------------------------------------------------------

||.|| = distance
------------------
||LM|| = Δx
||PQ|| = Δx(w)
||QR|| = |Δx(e)
||KJ|| = Δy
||TU|| = Δx(s)
||UV|| = Δx(n)
-------------------------------------------------------------------------------

CV ,CS
=====

CV : dimension = Δx*Δy
--------------------------------
CV : w ≤ x ≤ e &  s ≤ y ≤ n

CS : control surface
--------------------------
x-direction
face : || se,en || =  || ws,nw ||  = Δy
y-direction
face : || nw,en || =  || ws,se ||  = Δx

x-y plane : cartesian grid.
===================

========================================
pos : x-direction :
----------------------
in   :F(w) = u(w)*ρ(w) → Φ (w)
out :F(e)  = u(e)*ρ(e)   → Φ (e)
pos : y-direction :
----------------------
in   : F(s) = v(s)*ρ(s)  → Φ (s)
out : F(n)= v(n)*ρ(n) →Φ (n)
---------------------------------------------------------------------
interface :{Φ (w),Φ (e),Φ (s),Φ (n)},Φ(f), f : faces CV
nodal      :{Φ (W),Φ (p),Φ (E),Φ (S),Φ (E)}:  know values
---------------------------------------------------------------------
=======================================

use.ref1, 1D convection -diffusion
--------------------------------------------
start : p2

x-direction : u ( = positive direction)
-------------------

p2a:part(lhs(p1a),2) = part(rhs(p1a),2) ;

y-direction : v (= positive direction)
-------------------

p2b:part(lhs(p1a),1) = part(rhs(p1a),1) ;

x-direction :from : ref 1:({p4s4,p4s4a},{p4sa,p4sa1})
----------------------------------------------------------------
a [p]→ a1[p]

p3a:[D[e]=Γ(e)/Δx(e),D[w]=Γ(w)/Δx(w)];
p3b:subst([a[p]=a1[p]],[-F[w]/2+D[w]+F[e]/2+D[e]=a[p],F[w]/2+D[w]=a[W],D[e]-F[e]/2=a[E]]);

 rhs(p3b[1]);

p3c:makelist( rhs(p3b[i]) = lhs(p3b[i]) ,i,1,length(p3b));

y-direction :
----------------
a1[p]→ a2[p]

p3a1:subst([e=n,w=s,a1[p]=a2[p],Δx=Δy],p3a);
p3b1:subst([e=n,w=s,a1=a2,W=S,E=N],p3b);

p3c1:makelist( rhs(p3b1[i]) = lhs(p3b1[i]) ,i,1,length(p3b1));

→ combine :p3c & p3c1

p4:[a[p] =rhs(p3c[1])+rhs(p3c1[1])];

p4a:makelist([p3c[i],p3c1[i]],i,2,length(p3b));

p4b:flatten([p4,p4a]);

general eaquation: 2D problem .
-----------------------------------------
use : p4

rem : 1D ,ref 1: a[p]*Φ[p] = a[W]*Φ[W] + a[E]*Φ[E]

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

→ 2D

p5a:lhs(p5) = rhs(p5)+subst([W=S,E=N],rhs(p5));

rem :1D : D(.)

p5b1:flatten([p3a,p3a1]);

rem :1D : F(.) → F[.]

p5b2a:[F[e]=u(e)·ρ(e),F[w]=u(w)·ρ(w)];
p5b2b:lratsubst([e=n,w=s],p5b2a);

p5b2c:flatten([p5b2a,p5b2b]);

p5b3:p4b;

rem :

p6:makelist(p5b3[i],i,2,length(p5b3));

p6a:sum(p6[i],i,1,length(p6));

?  is

p7a:rhs(p5b3[1]);
p7b:rhs(p6a);

is(p7a=p7b);

because : take account ,continuity equation :G
-------------------------------------------------------------

p8a:p7a;
p8b: p7b+G;

p8c:p8a-p8b=0;

p8d :solve(p8c,G);

p9d1:a[p] =lhs(p6a)+G;
p9d2:subst(p8d,p9d1);

include : conservation law
-----------------------------------
in = out

X-DIRECTION (u) : ref 1,p3yy

p9:'diff((u·ρ),x,1)=0;
p9a:(F[w]-F[e]=0)*(-1);

Y-DIRECTION (v)

p10:subst([u=v,x=y],p9);
p10a:subst([w=s,e=n],p9a);

declare(ρ,constant)$
depends(u,x,v,y);

→ combine : vec(z) = u*e[x]+v*e[y]

s1:p9+p10;
s1a:factor(ev(s1,diff));
s2:p9a+p10a;

rem :out flow -in flow =0,see : Figure 1 " see :arrows ,around CV ,entering CS"
rem :ρ = constant ,∇.vec(z) = 0 ,for incompressible flow.

final :a[p] include "conservation law"
------------------------------------------------
use : p9d2

p11:lratsubst([p9a,p10a],p9d2);

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