kill(all)$

maxima : appendix : fvm :Grid generation for CV in 1D,version 1.0
peter.vlasschaert@gmail.com,31/10/2022

rem : i notation (....i-1,i,i+1,.....) , e,w notation (E,W,P)

ref : see website
----------------------
-------------------------------------------------------------------------------------------------------------------------------------------------
→ integration
new:20/10/2022:maxima :  fvm: code d(c)integrate(expr,var,lo,hi),version 1.0
→ discretization, first order differential
new:20/10/2022:maxima :  appendix : fvm : first order ,finite difference (central,forward,backward),version 1.0
→ coef. diffusion (taylor expansion)
new:20/10/2022:maxima :  appendix : fvm : first order ,finite difference (central,forward,backward),version 1.0
--------------------------------------------------------------------------------------------------------------------------------------------------

legend:
======
-------------------------------------
E  = east
W = west
N  = north (j notation)
S  = south  (j notation)
1D : x-direction  (i -notation) , see : Figure 1
2D : x-direction  (i -notation)
        y-dirction   (j -notation)
in   = flux through  CV
out = flux out the    CV

process : discritezation 1D domain (single : CV).
===================================

control :volume:CV
-------------------------

p1:CV = A*Δx;

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

p1a: u;

integration : from w: (x(i-1/2)) to e:(x(i-1/2))
--------------------------------------------------------

p1b1: Δx =x(i+1/2) -x(i-1/2);

p1b2:reverse(p1b1);

ref : integration definition: diffusion part (= of the transport equation)

dintegrate(pp1,xx1,aa1,bb1):=block([s1],
s1:integrate(pp1,xx1),
s2:at(s1,xx1=bb1)-at(s1,xx1=aa1),
return(s2)    
)$

p1c1:'diff((Γ(x)·('diff(u(x),x,1))),x,1);

p1c2a:'integrate(p1c1,x,w,e) =dintegrate(p1c1,x,w,e);

p1c2b:'integrate(p1c1,x,x(i-1/2),x(i+1/2))=subst([w=x(i-1/2),e=x(i+1/2)],dintegrate(p1c1,x,w,e));

discritezation : rhs "differential", use : central difference, diffusion coef.
--------------------------------------------------------------------------------------------
known values : nodal points,see :ref

→ difference :

w,e : notation

p1d1a:part(rhs(p1c2a),1,2)  = (u(E)-u(P))/Δx;
p1d1b1:part(rhs(p1c2a),2,0);
p1d1b2:part(rhs(p1c2a),2,1,2) = (u(P)-u(W))/Δx; ;

→ p1c2a

p1e1:lhs(p1c2a) = subst([p1d1a,p1d1b2],rhs(p1c2a));

i : notation (see :Figure1)

p1f:[E=i+1,W=i-1,P=i];

p1fa:subst(p1f,[p1d1a,p1d1b2]);

use : p1c2a

p1fb:lhs(p1c2b) =subst(p1fa,rhs(p1c2a));

→ diffusion :

last step : What about →Γ(e),Γ(w) 
--------------------------------------------
we need relations: nodal values (=E;W,P) and not interface values (e,w)

→ 1e) first order approximation : taylor series.

p1e:f(Δx+x)=(at('diff(f(x),x,1),x=m))·Δx+f(m1);

distance : Δx in the taylor series,m = value @ interface.

Γ(e) :

→ 2e)  m = e , Δx =Δx/2, f=Γ

p1ea:[m=e,Δx=Δx/2,f=Γ];

p1ea1:subst(p1ea,p1e);

set :Γ(Δx/2+x) =Γ(e)

p1ea2:subst([Γ(Δx/2+x) =Γ(e),m1=P],p1ea1);

p1ea3:[p1d1a,p1d1b2];

p1ea3a:subst([u=Γ],p1ea3);

final :Γ(e)

p1ea3b1:ratsimp(subst(p1ea3a,p1ea2));

final :Γ(w)

m=w,Δx=-Δx/2,f=Γ

p1ea3b2a:subst([e=w,E=W],p1ea3b1);

final : result : Diffusion part (from transport equation)

p1g1a:[p1ea3b1,p1ea3b2a];

use : p1e1

→ e,w

p1g1b:subst(p1g1a,p1e1);

→ i:notation

p1g1ca:flatten([p1f]);

use :p1g1a

p1g1cb:subst(p1g1ca,p1g1a);

p1g1cc:subst(p1g1cb,p1fb);

use : p1b1

p1g1cd:subst([p1b1],p1g1cc);

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