kill(all)$

maxima :diffusion+source =0,example:1,BC,,a[p],a[E]a[W],b,version 1.0
peter.vlasschaert@gmail.com,19/12/2022

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

load(lrats)$ ratprint:false$

x-dirction  → T(x) ; T= Temperature

depends([T],x);

check :

 dependencies;

px:A=1; /*area A=Δy=1*/

steady case : heat conduction with source (S = linear equations)
-----------------------------------------------------------------------------------

p1: 'diff(k*'diff(T,x),x)+S=0;

Φ = function ( T = temperature)
Γ = k  (heat ), 1D transport equation, k when T ,fourier.
k = thermal conductivity

no source : S = 0

p1a:[S=0];

p2:subst(p1a,p1);

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

a) need integration : variable from differential.
--------------------------------------------------------------

p3a1:args(part(lhs(p2)));

p3a2:args(p3a1[1]);

p3a3:part(args(p3a1[1])[1],1);

→ orginal version : wrong ode2

p3b:ode2(p2,T,x);

→ new version :ode2

p3b1:diff(k*T,x,2)=0;
p3b2:ode2(p3b1,p3a3,p3a1[2]);

solution : second order differential "p3b2"-> linear profile for T

b) discritization : central  difference.
------------------------------------------------
m = CS ,{e,w}

p4a:p3a2[1];

p4b1:at(p4a,x=m)*k(m);

p4b2:makelist(subst([m=i],p4b1),i,[w,e]);

p4b3:makelist(part(p4b2[i],2),i,1,length(p4b2));

fixed : value over total domain : 

p4b4: makelist(part(p4b3[i],0),i,1,length(p4b3));

→ uniform : k

p4b5:map("=",p4b3,p4b4);

→ method :central difference "differential"

p4bb1:p4b2[1]=(T[p]-T[W])/Δx(w)*k(w);
p4bb2:p4b2[2]=(T[E]-T[p])/Δx(e)*k(e);

p4bb3:[p4bb1,p4bb2];

use :

p4bb4:subst(p4b5,p4bb3);

p5:p4bb4[2]-p4bb4[1];

2e) general equation 1D
---------------------------

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

p5a:T[p]*a[p] =T[E]*a[E] +T[W]*a[W];

p5b:rhs(p5) =0;

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

⇒  non-uniform

p5b1:ratcoeff(p5b,T[p])*(-1)$
p5b1a:a[p] =(expand(lhs(p5b1)));
p5b2:ratcoeff(p5b,T[W])$
p5b2a:a[W] =(expand(lhs(p5b2)));
p5b3:ratcoeff(p5b,T[E])$
p5b3a:a[E] =(expand(lhs(p5b3)));

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

p5c:[p5b1a,p5b2a,p5b3a];

final :equation :use p5a

p5d: subst(p5c,p5a);

⇒  uniform

here : Δx = Ax(e) =Δ(w),interval is uniform.

p6:[Δx(e)=Δx,Δx(w)=Δx];

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

example 1: 1D steady , diffusion. ( = uniform )
-----------------------------------------
see : Figure 1

N:5;/*number CV*/

⇒ nodal values :  points (regular CV),{2,3,4}

→ p6,p5d

p7a:subst(p6,p5d);

p7b:[p=i,W=i-1,E=i+1];

p7c:subst(p7b,p7a);

definition : equation , regular : nodal values.
---------------------------------------------------------

eq[i]:=''p7c;

principal : build system of general equations
----------------------------------------------------------
eq [i]:= lhs(..) = rhs(..)

p8:makelist(eq[i],i,[2,3,4]);

left BC , T[W]=Tb, Δx(w) = Δx/2 ,Δx(e)=Δx  ⇒ {1}
-----------------------------------------------------------------

p9a:[T[W]=Tb,Δx(w) = Δx/2,Δx(e)=Δx,p=1,E=2];

p9a1:subst(p9a,p5d);

right BC , T[E]=Te, Δx(e) = Δx/2,Δx(w) = Δx   ⇒ {5}
---------------------------------------------------------------------

p9b:[T[E]=Te, Δx(e) = Δx/2,Δx(w) = Δx,p=5,W=4];

p9b1:subst(p9b,p5d);

final : result : example 1:
--------------------------------

final list : equations
--------------------------

p10:[p9a1,p8,p9b1];

p10a:flatten(p10);

p10a1:makelist(lhs(p10a[i]) -rhs(p10a[i]) =0,i,1,N);

p11:makelist(T[i],i,1,N);

p11a:coefmatrix(p10a,p11);

p11b:transpose(p11);

p11c:p11a.p11b;

p12:coefmatrix(p10a,flatten([Tb,p11,Te]));

p12a:transpose(flatten([Tb,p11,Te]));

p12b:p12.p12a;

p12c:transpose(makelist(0,i,1,N));

p12d: p12b =p12c;

p12d1:p12b[1,1];

check :lhs(p10a1) =lhs(p12b)

relation : matrix form & list of equations
------------------------------------------------------

p12d2:makelist(is(lhs(p12b[i,1])=lhs(p10a1[i])),i,1,N);

find :[T[1],T[2],T[3],T[4],T[5]]

→ solve :

p13:expand(solve(p10a1,p11)[1]);

temperature : x-direction 
--------------------------------

p13a:flatten([Tb=Tb,p13,Te=Te]);

p13aa:makelist(rhs(p13a[i]),i,1,N+2);

values to make : to make plot
--------------------------------------

x-direction

p14:makelist(x=Δx/2+(i-1)*Δx,i,1,N);

p14a:flatten([x=0,p14,x=5*Δx]);

values : Tb=100,Te=500,Δx=0.1,k = 1000.

p15:[Tb=100,Te=500,Δx=0.1,k = 1000];

x-values :

xxp:float(makelist(subst(p15,rhs(p14a[i])),i,1,N+2));

xpa:length(xxp);

y-values :

yyp:float(makelist(subst(p15,p13aa[i]),i,1,N+2));

ypa:length(yyp);

yx:makelist([xxp[ii],yyp[ii]],ii,1,N+2);

exact :solution :
===========

p16:p3b2;

sa:listofvars(p16);

p16b1:delete(sa[1],%);
p16b2:delete(sa[4],p16b1);

→ left : BC

x1:subst(p15,[T=Tb,x=0]);

p16a:subst(x1,p16);

→ right :BC

x2:subst(p15,[T=Te,x=N*Δx]);

p16b:subst(x2,p16);

p17:solve([p16a,p16b],p16b2);

p17a:ev(p16,p17);

p17b:rhs(p17a);

plot in range  "p17b" : 0 → 5*(.1) =0.5

rem :Menu :plot →plot2d 

plot :exact solution, "p17a"

wxplot2d([rhs(p17a)], [x,0,.5],[y,-100,500])$

values computing : fvm
fvm : finite volume method.

p17c:yx;

p17c1:length(p17c);

wxplot2d([[discrete, yx]],
        [style, points], [color, red],
        [point_type, asterisk],
        [legend, "computation=fvm"],
        [xlabel, "x"],
        [ylabel, "T=800*x+100"])$;

wxplot2d([[discrete, yx], p17b], [x,0,0.50],
        [style, points, lines], [color, red, blue],
        [point_type, asterisk],
        [legend, "fvm", "exact"],
        [xlabel, "x"],
        [ylabel, "T=800*x+100"])$

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