kill(all)$

maxima : convection -diffusion problem (exact : solution)integration,version 1.0
peter.vlasschaert@gmail.com,09/11/2022

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

load(lrats)$ ratprint:false$

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

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

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

bc : Φ(0) = 1  ... (1)
       Φ(L) = 0 ....(2)
bc = (b)oundary (c)onditions 

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

→lhs : p1

p1a1:part(lhs(p1),0);
p1a2:Jcond =part(lhs(p1),1);

cond =conduction

→rhs : p1

p1b1:part(lhs(p1),0);
p1b2:Jdif =-part(rhs(p1),1);

dif = diffusion

→definition " total flux : "
-----------------------------------

p2a:J = Jcond +Jdif;

assume : Γ constant

p2a1:[p1a2,p1b2];
p2b:subst(p2a1,p2a);

rem : d(J)/dx = 0 , J = c1, c1 = constant

→ definition a[.]*Φ(.) , S =SOURCE 
-----------------------------------------------
S = 0 , in p1

p3a:a[p]*Φ[p] = a[W]*Φ[W] + a[E]*Φ[E] +S;
p3b:a[p]*Φ[p] = a[nb]*Φ[nb] + S;

p3c:factor(rhs(-p3a+p3b));

p3c1:expand(solve(p3c,Φ[nb]))[1];

nb :neighbour hood of (p)

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

use : p2b

p4:rhs(p2b) = c1; /* see :rem above */;

p4a:Jcond-c1 = Jdif;

p4b:subst(p2a1,p4a);

→ separation : integration.

p4c:integrate(1/lhs(p4b) ,Φ)=integrate(1/part(rhs(p4b),1,1),x)+c2;

→ multiply : u*ρ ( = use ratcoeff).

p4c1:ratcoeff(lhs(p4b),Φ)*p4c;

→ operator : exp

p4c2:exp(lhs(p4c1)) =exp(rhs(p4c1));

→ solve :Φ

p4c3:solve(p4c2,Φ)[1];

p4c4:rhs(p4c3);

1e) Φ is function of x
----------------------------
use : p4c4

p5:Φ(x)=''p4c4 ;

2e) find :{c1,c2}
--------------------

→ find :c2

p6a:subst([x=0],p5);
p6b:subst([x=L],p5);

p6c1:expand(p6b-p6a);
p6c2:factor(expand(ratcoeff(lhs(p4b),Φ)*p6c1));

p6c3:log(p6c2);

p6c4:log((Φ(L)-Φ(0))*u*ρ)=log(%e^(c2*u*ρ)*(%e^((L*u*ρ)/Γ)-1)), logexpand=super;

p6c5:solve(p6c4,c2)[1];

p6c6:[lhs(p6c5) =logcontract(rhs(p6c5))];

rem : math : ln -> maxima : log

→  find : c1

p7:expand(subst([p6c6],p5));

solve : c1

p7:factor(solve(p7,c1));

subst :{c1,c2} into p5
---------------------------

p8:flatten([p6c6,p7]);

check :

p8a:factor(subst(p8,p5));

3e) find : (Φ(x) -Φ(0))/(Φ(L) -Φ(0))
---------------------------------------

p9:p5-p6a;

p9a:factor(subst(p8,p9));

4e) exact solution (p1) : ode (Γ =constant), BC (" see:above ")
--------------------------------------------------------------------------

p9b:p9a/lhs(p6c1);

p9c:solve(p9a,Φ(x))[1];

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

interface :[e]

→ Φ(e)

p10a:[Φ(x)=Φ(e),Φ(0)=Φ(p),Φ(L)=Φ(E),L=Δx(e),x=Δx];

p10b:subst(p10a,p9b);

interface :[w]

p11a:[Φ(x)=Φ(w),Φ(0)=Φ(W),Φ(L)=Φ(p),L=Δx(e),x=Δx];
p11b:subst(p11a,p9b);

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)

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