kill(all)$

maxima :fvm : generalized relation convection-diffusion (Pe=Peclet) ,version1.0a
peter.vlasschaert@gmail.com,17/12/2022

see website :

ref:maxima :  convection -diffusion problem ((constraint =continuity equation)),version 1.0

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

load(lrats)$

general mathematical model : convection-diffusion model.
==========================================

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

Jx = convection -diffusion

p2a:Jtot =lhs(p1*Δx(f)/Γ);
p2b:Jtot =expand(rhs(p1*Δx(f)/Γ));

def inition : peclet cell number ,x-direction
--------------------------------------------------------

p2c:reverse(Pe=ratcoeff(rhs(p2b),Φ));

p2d:lratsubst([p2c],p2b);

rem : Δx(f), f =e,w ,constant value outside the differential.
see below : convection →  with parameter α       ,see below.
                   diffusion     →  difference  in general ,see below. (differential )

general scheme :{α,β} as scheme. ( = discretization scheme)
===========================================

x-direction : u> 0 :flow : Φ[i]→Φ[i+1]
-------------------------------------------------

rem : Φ[i] →[f]→Φ[i+1], f = e,w

p3a:Jconv = Pe*(α*Φ[i]+(1-α)*Φ[i+1]);
p3b: Jdiff   =β*(Φ[i]-Φ[i+1]);

p3: Jtot =lhs(p3a-p3b);

p4:Jtot =rhs(p3a-p3b);

A=A(Pe), B=B(Pe)

p5: Jtot =B*Φ[i] -A*Φ[i+1];

→ conversion : p4 & p5

1e)  find B 
-------------

p6a:B(Pe)=ratcoeff(rhs(p4),Φ[i]);

2e) find A
----------------

p6b:A(Pe)=collectterms(ratcoeff(rhs(p4),-1*Φ[i+1]),Pe);

relation : B(Pe) & A(Pe)
-------------------------------

p6c1:ratsimp(p6a-p6b);

p6c:solve(p6c1,B(Pe))[1];

a) special case : Φ[i] =Φ[i+1]
------------------------------------

use :p5,relation A ? B

p7a:facsum(subst([Φ[i] =Φ[i+1]],p5));

use : p4 ?

p7b:ratsimp(subst([Φ[i] =Φ[i+1]],p4));

result :only convection , because no β in Jtot.

p7c:[A=A(Pe),B=B(Pe)];/*see : above*/

→ p7a =p7b

p7c1:subst(p7c,rhs(p7a)=rhs(p7b));

p7c2:solve(p7c1,B(Pe))[1];

check :p6c & p7c2

p7c3: is(rhs(p6c)=rhs(p7c2));

conclusion  : Φ[i] = Φ[i+1]
-----------------------------------
1e) only convection  "p7b"
2e) general relation : B(Pe) = function(A(Pe))
-----------------------------------------------------------------
rem : Pe depent flow direction , contains  u inside 

→ reverse : flow direction
===================

x-direction : u< 0 :flow : Φ[i+1]→Φ[i]
-----------------------------------------------------
1e) replace : Φ[i] → Φ[i+1] & Φ[i+1] → Φ[i] in "p5"
2e) Pe → -Pe in 2e)
3e) front sign = - " to get equation between two direction "
two direction ={u>0 & u<0}

[Φ[i] = Φ[i+1],Φ[i+1] = Φ[i] ]

use : rhs(p5)

orginal

p8a:rhs(p5);

p8aa:subst(p7c,p8a);

→ 1e) replace

p8a1:psubst([Φ[i] = Φ[i+1],Φ[i+1] = Φ[i] ],p8aa);

→ 2e) Pe → -Pe, "see above :rem,depent flow direction"

p8a2: subst([Pe =-Pe],p8a1);

→ 4e) equation between : p8aa ,we need front sign before : p8a2

p8a3 : p8aa=(-1)*p8a2;

p8a31:expand(p8a3);

final result : flow direction depent
---------------------------------------------

p8a5a:ratcoeff(p8a31,Φ[i]);
p8a5b:ratcoeff(p8a31,Φ[i+1])*(-1);
p8a5c:[p8a5a,p8a5b];

analysis : A(Pe) & B(Pe)
-------------------------------

p9a:p7c2;
p9b:p8a5a;
p9c:p8a5b;

rem : Δx(e) →   p=Φ[i] & E=Φ[i+1]
rem : Δx(w) → W=Φ[i] & p=Φ[i+1]

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

genaral equation: for a scheme.
-----------------------

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

F() = F[]
------------
rem : relation convection term , F(w),F(e)
D() =D[]
------------
rem : relation diffusion term ,   D(w),D(e)

relation between : Jw = Je
----------------------------------
conservation law → Jw = Je

 Δx(e) →   p=>Φ[i] & E=>Φ[i+1]
 Δx(w) →  W=>Φ[i] & p=>Φ[i+1]

use :p8aa

p10a:p8aa;

→ Je, e = out for CV

p10b1:[Pe=Pe(Δx(e)),Φ[i] =Φ[p],Φ[i+1] =Φ[E]];
p10b:subst(p10b1,p10a);

→ Jw,w = in for CV

p11b1:[Pe=Pe(Δx(w)),Φ[i] =Φ[W],Φ[i+1] =Φ[p]];
p11b:subst(p11b1,p10a);

dimensioless:: Jw=Je
---------------------------

p12:p10b =p11b;

use : general equation. (x-direction),p9d,cartesian

p12a1:ratcoeff(p12,Φ[p]);
p12a11:a[p]=lhs(p12a1)-rhs(p12a1);
p12a2:ratcoeff(p12,Φ[W]);
p12a22:a[W]=rhs(p12a2)-lhs(p12a2);
p12a3:ratcoeff(p12,Φ[E]);
p12a33:a[E]=rhs(p12a3)-lhs(p12a3);

p12b:[p12a11,p12a22,p12a33];

use :p9d

p12c:subst(p12b,p9d);

non-dimension : general equation.

p12d:p9d;

rem : Pe(.) =F(.)/D(.), cell peclet number

p13:[A(Pe(Δx(e))) = a[E]/D[e],B(Pe(Δx(w)))=a[W]/D[w]];

p13a:makelist(rhs(p13[i])=lhs(p13[i]),i,1,length(p13));

see why : ref 1 (exponential : scheme)

p14a:a[p]=F[w]/(%e^Pe(Δx(w))-1)+F[e]·%e^Pe(Δx(e))/(%e^Pe(Δx(e))-1);
p14b:a[W]=F[w]·%e^Pe(Δx(w))/(%e^Pe(Δx(w))-1);
p14c:a[E]=F[e]/(%e^Pe(Δx(e))-1);

definition : Peclet Cell number
----------------------------------------

→ see above ;

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

→ face : m=w,m1=Δx(w)

p15a1:subst([m=w,m1=Δx(w)],p15a);

p15a1x:p14b/D[w];
p15a1xx:lratsubst([reverse(p15a1)],(p15a1x));/*result number*/
p15aa:lratsubst(p13a,%);
p15aa1:subst([Pe(Δx(w))=x],rhs(p15aa));

f1(x):=''p15aa1;/*blue : color in plot*/

→ face : m=e,m1=Δx(e)

p15a2:subst([m=e,m1=Δx(e)],p15a);

p15a2x:p14c/D[e];
p15a2xx:lratsubst([reverse(p15a2)],(p15a2x));/*result number*/
p15bb:lratsubst(p13a,%);
p15bb1:subst([Pe(Δx(e))=x],rhs(p15bb));

f2(x):=''p15bb1;/*red : color in plot*/

reason:all numbers = {Pe(Δx(e)),Pe(Δx(w)),%e^Pe(Δx(e)),%e^Pe(Δx(w))}

make plot : relation between : A(Pe) & B(Pe)
----------------------------------------------------------

wxplot2d([f1(x),f2(x)], [x,-5,5],
 [gnuplot_postamble, "set size ratio 1; set zeroaxis;"])$

q1:wxplot2d([f1(x),f2(x)], [x,-5,5],
 [gnuplot_postamble, "set size ratio 1; set zeroaxis;"])$

blue : B(Pe(Δx(w))), face w of CV
red   : A(Pe(Δx(e))), face e   of CV 

plot : reverse direction ,relation from A & B
----------------------------------------------------------

1e) A = f2(x), x = Peclet number =Pe(.)
--------------------------------------------------

result :A(Pe) =f2(Pe)

q2:wxplot2d([f2(x)], [x,-5,5],
 [gnuplot_postamble, "set size ratio 1; set zeroaxis;"])$

from properties : A & B
=================
B(Pe)=A(Pe)+Pe , B(Pe)=A(-Pe)  ,  A(Pe)=B(-Pe)

q3:wxplot2d([f2(x),f2(x)+x], [x,-5,5],
 [gnuplot_postamble, "set size ratio 1; set zeroaxis;"])$

blue : A(Pe) red : B(Pe)

p16: ratsimp(f2(x)+x);

rem : see relation A = function(B) or  visaversa
rem : plot q1 same as q2,reason to see direct same : intersect for x=0

what when we reverse the flow Pe → -Pe
-------------------------------------------------------

q3:wxplot2d([f2(x),f2(-x)], [x,-5,5],
 [gnuplot_postamble, "set size ratio 1; set zeroaxis;"])$

blue : A(Pe) ,red : A(-Pe) 
------------------------------
==========================================
relations :B(Pe)=A(Pe)+Pe , B(Pe)=A(-Pe)  ,  A(Pe)=B(-Pe)
==========================================
follows : A(-Pe) =B(Pe), show mathematical computation above.

info :three relation

p17a:B(Pe)=A(Pe)+Pe;
p17b:B(Pe)=A(-Pe);
p17c:A(Pe)=B(-Pe);

relation : A(Pe) & A(|Pe|) =A(abs(Pe))
===========================

p18: A(abs(Pe));

step 1:assume : Pe>0

p18a:A(Pe)  = p18;

step 2:assume : Pe < 0

use :p17a

p18b:reverse(solve(p17a,Pe)[1]);

use :p17b

p18c:subst([p17b],p18b);

p18c1:solve(p18c,A(Pe))[1];

when : Pe < 0 

p18d:lhs(p18c1) =rhs(p18a)-Pe;

p18d1:part(rhs(p18d),2);

result :
Pe > 0  → A(Pe)=A(abs(Pe))
Pe < 0 →  A(Pe)=A(abs(Pe))-Pe

combine :use max function: for any Pe
---------------------------------------------------

p18e: A(Pe)=A(abs(Pe))+max(p18d1,0);

assume(Pe>0);/*check pos.values*/

p18f:ev(p18e,nouns);

forget(Pe>0);

check : neg.values
------------------------

assume(Pe<0);
p18f1:ev(p18e,nouns);
forget(Pe<0);

relation : B(Pe) & A(|Pe|) =A(abs(Pe))
===========================

p19:B(abs(Pe));

use :p17b

p19a:p17b;

step 1:assume : Pe<0

p19b:lhs(p19a) =rhs( p18a);

step 2:assume : Pe>0

use:p17a

p19c:p17a;

→ Pe>0

p19d:lhs(p19c) =p18+part(rhs(p19c),2);

p19d1:part(rhs(p19d),2);

result :
Pe > 0  → B(Pe)=A(abs(Pe))+Pe
Pe < 0 →  B(Pe)=A(abs(Pe))

combine :use max function: for any Pe
---------------------------------------------------

p19e:lhs(p19d) =part(rhs(p19d),1) +max(p19d1,0);

check :

→ pos value : Pe  →p19e

assume(Pe>0);
ev(p19e,nouns);
forget(Pe>0);

→ neg value : Pe  →p19e

assume(Pe<0);
ev(p19e,nouns);
forget(Pe<0);

info : Menu → Maxima → variables

values;

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