;

maxima :appendix:fvm:stencil for finite volume method 2D (heat),version 1.0
peter.vlasschaert@gmail.com,31/10/2022

load(stringproc)$
load(lrats)$

ordergreat(a,b,c,d)$

step 1:  general equation : stencil equation

p1:T(x,y):=a+b*x+c*y+d*x*y;

step 2: find {a,b,c,d}

→ cell 1: dimension (Δx,Δy)

1e) corners: 0,1,2,5

pa1:TT0=T(0,0);
pa2:TT1=T(Δx,0);
pa3:TT2=T(0,Δy);
pa4:TT5=T(Δx,Δy);

2e) find corners:back from the equations above,list

→ corner : numbers

pa5:makelist(ev(sconcat(pa,i),nouns),i,1,4);

stringp(pa5[1]);

pa6:makelist(eval_string(pa5[i]),i,1,length(pa5));

atom(lhs(pa6[1]));

rem :convert atom to string : ""

stringp("lhs(pa6[1])");

sa:eval_string("lhs(pa6[1])");

string(%);

stringp(%);

sremove("TT",string(sa));

→ string :list: corner numbers

pa7:makelist(eval_string("lhs(pa6[i])"),i,1,4);

pa71:map(string,pa7);

sremove("TT",pa71[1]);

pa72:makelist(sremove("TT",pa71[i]),i,1,length(pa5));

check: string

pa72[2]+pa72[4];

→ number :list: corner numbers

pa8:map(eval_string,pa72);

check : numbers

pa8[2]+pa8[4];

3e) find variables.

pa9:T(1,1);

pa91:makelist(part(pa9,i),i,1,length(pa5));

→ solve : eqs:pa6 for pa91

pa10:solve(pa6,pa91);
pa101:flatten(pa10);

step 3: find T(x,y)

pa11:subst(pa101,T(x,y));

cell 1

pa12:T1(x,y):=''pa11;

computation : q15(out),q25(out)
=======================
→ fourier law 

q(ii,jj,AA,Tf):=-k*AA*Tf;

rem : nn = "normal direction flux"
          Tf = ('diff(Tf,nn,1))

rem q15:q(1,5,Δy/2,..)

→ q15   x-coordinate = Δx/2,  y-coordinate = Δy/4
      nn = x direction

first step :('diff(Tf,nn,1))

pb1:'diff('T1(x,y),x)=diff(T1(x,y),x);

pb1a:ratsimp(subst([x=Δx/2,y=Δy/4],rhs(pb1)));

pb1aa:at(lhs(pb1),[x=Δx/2,y=Δy/4]);

pb1aaa:pb1aa=pb1a;

second step: q

pb2: 'q(1,5,Δy/2,pb1a)=(q(1,5,Δy/2,pb1a));

rem q25:q(2,5,x,Δx/2,..)

→ q25   x-coordinate = Δx/4,  y-coordinate = Δy/2
       nn = y direction

first step :('diff(Tf,nn,1))

pc1: 'diff('T1(x,y),y)=diff(T1(x,y),y);

pc1a:ratsimp(subst([x=Δx/4,y=Δy/2],rhs(pc1)));

pc1aa:at(lhs(pc1),[x=Δx/4,y=Δy/2]);

pc1aaa:pc1aa=pc1a;

second step: q

pc2: 'q(2,5,Δx/2,pc1a)=(q(1,5,Δx/2,pc1a));

cell 2.
q26: x-coordinate = -Δx/4,y-coordinate = Δy/2,AA=Δx/2,nn=y
q36: x-coordinate = -Δx/2,y-coordinate = Δy/4,AA=Δy/2,nn=x

cel22:[x= -Δx/4,y=Δy/2,AA=Δx/2];
cel23::[x= -Δx/2,y=Δy/4,AA=Δy/2];
xx:pa91;

→ equations : needed cel 2

pd1:TT3=T(-Δx,0);
pd2:TT6=T(-Δx,Δy);

→ TT[0,2,3,5]

pd3:[pa1,pa3,pd1,pd2];

→ solve equations ,{a,b,c,d}

pd4:flatten(solve(pd3,pa91));

→  input in T(x,y)

pd5:subst(pd4,T(x,y));

→ T2 (equation : cell 2)

pd6:T2(x,y):=''pd5;

→ cell26

→ normal direction:y  , differentation  T2 

pd7a1:'diff(T2(x,y),y)=diff(T2(x,y),y);

pd7a2:subst(cel22,rhs(pd7a1));

→ final result : q26

pd7a3:'q(2,6,Δy/2,pd7a2)=ratsimp(q(2,6,Δy/2,pd7a2));

→ cell23

→ normal direction:x  , differentation  T2 

pe7a1:'diff(T2(x,y),x)=diff(T2(x,y),x);

pe7a2:subst(cel23,rhs(pe7a1));

→ final result : q36

pe7a3:'q(3,6,Δy/2,pe7a2)=ratsimp(q(3,5,Δy/2,pe7a2));

cell 3.
q37: x-coordinate = -Δx/2,y-coordinate = -Δy/4,AA=Δy/2,nn=x
q47: x-coordinate = -Δx/4,y-coordinate =- Δy/2,AA=Δx/2,nn=y

cel33:[x= -Δx/2,y=-Δy/4,AA=Δy/2];
cel34::[x= -Δx/4,y=-Δy/2,AA=Δx/2];

→ equations : needed cel 3

pf1:TT4=T(0,-Δy);
pf2:TT7=T(-Δx,-Δy);

→ TT[0,3,4,7]

pf3:[pa1,pd1,pf1,pf2];

→ solve equations ,{a,b,c,d}

pf4:flatten(solve(pf3,pa91));

→  input in T(x,y)

pf5:subst(pf4,T(x,y));

→ T3 (equation : cell 3)

pf6:T3(x,y):=''pf5;

→ cell33

→ normal direction:x  , differentation  T3

pf7a1:'diff(T3(x,y),x)=diff(T3(x,y),x);

pf7a2:subst(cel33,rhs(pf7a1));

→ final result : q37

pf7a3:'q(3,7,Δy/2,pf7a2)=ratsimp(q(3,7,Δy/2,pf7a2));

→ cell34

→ normal direction:y  , differentation  T3

pf8a1:'diff(T3(x,y),y)=diff(T3(x,y),y);

pf8a2:subst(cel34,rhs(pf8a1));

→ final result : q47

pf8a3:'q(4,7,Δx/2,pf8a2)=ratsimp(q(3,7,Δx/2,pf8a2));

cell 4.
q18: x-coordinate = Δx/2,y-coordinate = -Δy/4,AA=Δy/2,nn=x
q48: x-coordinate = Δx/4,y-coordinate =- Δy/2,AA=Δx/2,nn=y

cel41:[x= Δx/2,y=-Δy/4,AA=Δy/2];
cel44::[x= Δx/4,y=--Δy/2,AA=Δx/2];

→ equations : needed cel 4

pg1:TT8=T(Δx,-Δy);

→ TT[0,1,4,8]

pg3:[pa1,pa2,pf1,pg1];

→ solve equations ,{a,b,c,d}

pg4:flatten(solve(pg3,pa91));

→  input in T(x,y)

pg5:subst(pg4,T(x,y));

→ T4 (equation : cell 4)

pg6:T4(x,y):=''pg5;

→ cell41

→ normal direction:x  , differentation  T4

pg7a1:'diff(T4(x,y),x)=diff(T4(x,y),x);

pg7a2:subst(cel41,rhs(pg7a1));

→ final result : q18

pg8a3:'q(1,8,Δy/2,pg7a2)=ratsimp(q(1,8,Δy/2,pg7a2));

→ cell44

→ normal direction:y  , differentation  T4

ph7a1:'diff(T4(x,y),y)=diff(T4(x,y),y);

ph7a2:subst(cel44,rhs(ph7a1));

→ final result : q48

ph8a3:'q(4,8,Δx/2,ph7a2)=ratsimp(q(4,8,Δx/2,ph7a2));

Total balances 
===========

q15+q18+q25+q26-q36-q37-q47-q48 =0
out (flux)-in(flux) =0

rem: heat flux : fourier laws of conduction

step 1 : set α

pl1a:α=Δx/Δy;
pl1b:rhs(pl1a)=lhs(pl1a);

step 2: balances

→ non unity grid aspect ratio "α"

pl2:pb2+pc2+pd7a3-pe7a3-pf7a3-pf8a3+pg8a3-ph8a3=0$

pl2a:map(rhs,pl2)$

pl3a:lratsubst([pl1b],pl2a)/k;

→ unity grid aspect ratio "α=1"

pl4a:subst([α=1],pl3a);

pl5a:pl4a*(-8);

pl5b:factor(pl5a);

pl5c:part(lhs(pl5b),1);

General stencil equation
---------------------------------

result : stencil .
===========

pl5d:pl5b/pl5c;

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