kill(all)$

maxima :fvm:appendix:double integration diffusion term (spatial,time),part 2,version 1.0a
peter.vlasschaert@gmail.com,28/11/2022

see website :
ref : 1
new:20/10/2022:maxima :  fvm: code d(c)integrate(expr,var,lo,hi),version 1.0

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

depends(T,x);

diffusion term
--------------------

p1: 'diff(k*'diff(T,x),x);

Φ = function ( T = temperature)
Γ = k  (heat )

p1a:args(p1);

→ integration variable:from "w→e"

p1b:p1a[2];

p1b:args(p1)[2];

see : ref 1, "above"

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

1e) integration : between interface ([w]&[e])
------------------------------------------------------------

result : integration over the CV

p2a:dintegrate(p1,x,w,e);

p2b:part(p2a,2,1);
p2c:part(p2a,1);

integration : diffusion model

How to integrate : use  "dcintegrate"

p2d:[p2b,p2c];

2e) discritization : central  difference. (=CDS)
------------------------------------------------

p3a:integrate(p1,x);

p3b:[((T[p]-T[W])*k(w))/Δx(w),(k(e)*(T[E]-T[p]))/Δx(e)];

3e) combine 
-----------------

p3a:map("=",p2d,p3b);

4e) integration p1: spatial

p3ax:subst(p3a,p2a);

p3a1a:listofvars(p3ax);

p3a1b:makelist(p3a1a[i],i,2,length(p3a1a)-1);

p3a1b1:part(p3a1b[1],0);
p3a1b2:part(p3a1b[1],1);

p3a1:makelist(part(p3a1b[i],1),i,1,length(p3a1b));

p3a1x:T[r]^°;

p3a1xx:makelist(subst([r=v],p3a1x),v,p3a1);

p3a1xxx:map("=",p3a1b,p3a1xx);

simp:false;

representation : spatial integration

p3a2:'integrate(p1,x,w,e) =p3ax;

CV = see Figure 1

rem : "old time level",see below :T[r]^°,rhs(p3a2)

p3a2y:expand(subst(p3a1xxx,rhs(p3a2)));

5e) final integration p1:time ,t→t+Δt

p4a:'integrate(lhs(p3a2),t,t,t+Δt)='integrate(rhs(p3a2),t,t,t+Δt);

princepal time integration.
-----------------------------------

p5:'integrate(T[r],t,t,t+Δt);

p5a:ratsimp(ev(p5,nouns));

p5b:p5=subst([T[r]=θ*T[r]+(1-θ)*T[r]^°],p5a);

p5c:makelist(subst([r=xx],p5b),xx,[W,p,E]);

p5c1:makelist(subst([r=sa],T[r] = rhs(p5b)),sa,[W,p,E]);

θ values
------------

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

θ =0 : old time level (=t)

p5c1a:ratsimp(subst([θ =0],p5c1));

θ =1/2 : equal weighted  (between(t & t+Δt))

p5c1b:ratsimp(subst([θ =1/2],p5c1));

θ =1 : new time level  (=t+Δt)

p5c1c:ratsimp(subst([θ =1],p5c1));

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

final time integration : use {p3a1,p5c1}

p5d:lhs(p4a) =subst(p5c1,rhs(p3a2));

→ simplification : not from "p5d"

p5dx:lhs(p5d) =(θ*rhs(p3a2)+(1-θ)*p3a2y)*Δt;

→ θ time level (= double integration)

θ =0 : old time level (=t), used Explicit scheme

p5d1a:lhs(p5d) =expand(ratsimp(subst([θ =0],rhs(p5d))));

θ =1/2 : equal weighted  (between(t & t+Δt)),used Crank -Nicolson scheme

p5d1b:lhs(p5d) =expand(ratsimp(subst([θ =1/2],rhs(p5d))));

θ =1 : new time level  (=t+Δt),used Implicit scheme

p5d1c:lhs(p5d) =expand(ratsimp(subst([θ =1],rhs(p5d))));

rem : sum, integrate
---------------------------

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

properties(sum);

properties(integrate);

declare(integrate,linear);

properties(integrate);

general : rem
-----------------

sum : symbol usefull (= Σ)
simp : false → can't evaluation integral on the lhs of expression.

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

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