kill(all)$

maxima :fvm:double integration one dim.heat conduction (unsteady) with source =S,version 1.0a
peter.vlasschaert@gmail.com,28/11/2022

load packages
-------------------

load(lrats)$

rem : unsteady : time dependence of temperature (=T)

see website :
ref : 1
new:27/11/2022:maxima : appendix : fvm : ,double integration (time,cv) for time differential (part 1),version 1.0
ref : 2
new:28/11/2022:maxima : appendix : fvm : ,double integration (time,cv) for diffusion term (part 2),version 1.0a
ref : 3 (a[p]=a[W]+a[E] =sum(a[nb])
new:12/11/2022:maxima : :convection -diffusion relation peclet cell number (max :limit),version 1.0
ref : 4
new:16/11/2022:maxima : :simple formulation of 4 basic rules , 1D  fvm (W,w,p,e,E) ,CV =[w,p,e],version1.0

rem : neighbourhood cells =nb, see: 1D ,nb =W & E

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

simp:false$

model : unsteady 1dim. heat conduction
-----------------------------------------------------
c = specefic heat of material
ρ = density
t = time
k =specefic heat conduction ( = fourier law)
x = 1 dim. 

rem : T=T(x,t)

p1:ρ*c*'diff(T,t) = 'diff(k*'diff(T,x),x)+S;

'integrate(lhs(p3a2),t,t,t+Δt)='integrate(rhs(p3a2),t,t,t+Δt)
:'integrate(p1,x,w,e) =p3ax;

→ double integration : model

1e) integration ,space : x  (Control Volume =ΔV = A*Δx),see figure.
---------------------------------

p1a:'integrate(lhs(p1),x,w,e) ='integrate(part(rhs(p1),1),x,w,e)+barS*ΔV ;

p1a1:part(rhs(p1a),1);

p1a2:part(rhs(p1a),2);

2e) integration,time : t  (between(t& t+Δt))

p1b:'integrate(lhs(p1a),t,t,t+Δt)='integrate(rhs(p1a),t,t,t+Δt);

p1c:lhs(p1b)='integrate(p1a1,t,t,t+Δt)+'integrate(p1a2,t,t,t+Δt);

see :ref 2: "above"
-----------------------
use :p5dx

p1c1:part(rhs(p1c),1);

→ simplification

p1c22:part(rhs(p1c),2) =barS*A*Δx*Δt;

p1c2a:rhs(p1c2);

see : ref 1, "above"
-----------------------
use : p4d

p1c23:lhs(p1b);

model : after integration 

A:1;

reason : A not implimented in diffusion term

→ lhs : model

see :ref 2

p2a:p1c1=(θ·(k(e)·(T[E]-T[p])/Δx(e)-(T[p]-T[W])·k(w)/Δx(w))+(1-θ)·(k(e)·1/Δx(e)·(T[E]^°+(-1)·T[p]^°)+(-1)·((-1)·T[W]^°+T[p]^°)·k(w)·1/Δx(w)))*Δt;

p2b:lhs(p1c22)=barS*Δt*Δx;

→ rhs : model

see:ref 1

p2c:p1c23 =c·ρ·A·(T[p]-T[p]^°)·Δx;

simp:true;

model : computation
----------------------------

p3a:rhs(p2c) =rhs(p2a)+rhs(p2b);

→ simplification : collectterms , /Δt

p3a1:lhs(p3a) =collectterms(rhs(p3a),Δt);

p3a2:p3a1/Δt;

model:general equation:
--------------------------------

a[p]*T[p] =a[W]*T[W]+a[E]*T[E]+a[p]^°*T[p]^°+b;

find : coef. T[p]
---------------------

p4a:ratcoeff(p3a2,T[p]);

p4a1:lhs(p4a) = collectterms(rhs(p4a),θ);

p4af:a[p] =lhs(p4a1)-rhs(p4a1);

find : coef. T[W]
------------------------

find :coef :T[W]^°

p4b:ratcoeff(p3a2,T[W]^°);

p4b1:facsum(rhs(p4b),k(w))*T[W]^°;

find :coef :T[W]

p4b2:ratcoeff(p3a2,T[W]);

p4b3:rhs(p4b2)*T[W];

→ combine

p4b4:a[W] = combine(p4b3+p4b1);

find : coef. T[E]
------------------------

find :coef :T[E]^°

p5b:ratcoeff(p3a2,T[E]^°);

p5b1:facsum(rhs(p5b),k(e))*T[E]^°;

find :coef :T[E]

p5b2:ratcoeff(p3a2,T[E]);

p5b3:rhs(p5b2)*T[E];

→ combine

p5b4:a[E] =combine( p5b3+p5b1);

find : coef. T[p]^°
---------------------

p5a:ratcoeff(p3a2,T[p]^°);

p5a1:a[p]^°=expand( rhs(p5a))-(lhs(p5a));

find : coef:

p6:b=last(rhs(p3a2));

rem : equation not good one a[p]: because "seen before"

a[p] =f(a[W],a[E],a[p]^°)

see : ref 3,why ?

use :p4af →

p7a:p4af;

because : Δt, first part (depent : on time)

p7a1:a[p]^° = part(rhs(p7a),1);

find(a[W],a[E])

p7a2:expand(part(rhs(p7a),2));

p7a3:factor(a[W]*θ+a[E]*θ);

→ final result :a[p]

p7af:a[p] =lhs(p7a1)+p7a3;

pa7af1:a[W]=(ratcoeff(part(p7a2,1),θ));
pa7af2:a[E]=(ratcoeff(part(p7a2,2),θ));

back to : p5a1

→ find :coef : T[p]^°
----------------------------

p8a:rhs(p5a1);

p8a1:[pa7af1,pa7af2,p7a1];

→ reverse :direction

p8a2:makelist(rhs(p8a1[i])=lhs(p8a1[i]),i,1,length(p8a1));

p8af:expand(lratsubst(p8a2,p8a));

p8af1:p8af*T[p]^°;

find coef : T[W]
---------------------
use : p8a2

p9a1:rhs(p4b4);

p9af:lratsubst(p8a2,p9a1);

p9af1:facsum(p9af,a[W]);

find coef : T[E]
---------------------
use :p8a2

p9b1:rhs(p5b4);

p9bf:lratsubst(p8a2,p9b1);

p9bf1:facsum(p9bf,a[W]);

final equation : f(θ)
--------------------------

p10a:a[p]*T[p] = p9af1+p9bf1+p8af1+b;

p10b:p8a1;

p10c:p7af;

p10d:p6;

rem :barS "mean value over CV",A=1 ( = area CV)

θ   values :
--------------

θ  = 0   , Explicit Scheme
θ  = 1/2, Crank-Nicolson Scheme
θ  = 1   , Implicit Scheme

→ θ  = 0   , Explicit Scheme

p10af:subst([θ=0],p10a);
p10af1:subst([θ=0],p10c);

→ θ  = 1/2, Crank-Nicolson Scheme

p10bf:subst([θ=1/2],p10a);
p10bf1:subst([θ=1/2],p10c);

→ θ  = 1   , Implicit Scheme

p10cf:subst([θ=1],p10a);
p10bf1:subst([θ=1],p10c);

linear source term : S[Φ]+S[p]*T[p]
==========================

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

p11a: b =S[Φ];
p11b:lhs(p10c) =rhs(p10c) -S[p];

→ see : ref 3,ref 4

negative sign (S[p])

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

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