maxima:general  cv 2D steady heat conduction.(fourier law),finite difference & fvm,version 1
peter.vlasschaert@gmail.com,31/10/2022

ref :An introduction to computational fluid dynamics
       The Finite Volume Method (=  FVM ),second edition
       H.versteeg
       Pearson
ref : figures  are constructed with the use of "geogebra" classic.

load packages,option
==========

load(lrats)$ ratprint:false$

general interior node "laplace equation = ΔT",2D version
==========================================

pp1: 'diff(T,x,2)+'diff(T,y,2) =0;
pp2: ΔT = 0;

pp1 or pp2 or steady state  "independent of time" = dT/dt =0
T = Temperature
t  =  time

1D:  x - direction 
============
cv = control volume = Δx*1*1 = A*Δx

2D . x - direction & y-direction
======================
cv =control volume = Δx*Δy*1

3D : x - direction & y-direction & z-direction
================================
cv = control volume = Δx*Δy*Δz  

reason : volume general 3D use of the ones (= unity deep)
FVM = FINITE VOLUME METHOD (see: website)
rem : A = surface can be unit surface

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

rem : FEM : FINITE ELEMENT METHOD       (see:website)
rem : FDM : FINITE DIFFERENCE METHOD  (see:website)

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

general interior node "laplace equation = ΔT", control volume ={A1,B1,C1,D1}

dotqi, i = integer  , "FOURIER LAW"  = Transfer of Energy
==========================
cv = control volume around (i,j)

p1:'diff(q[i],t)=-k*AA*'diff(T,nn);

AA = surface cv
k     = heat conduction coef. ( medium )
nn   = outward normal on the surface in consideration
t      = time 
T     = Temperature

→ constraint : cv
m= number of heat flows  cv

p2:sum('diff(q[i],t),i,1,m) = 0;

rem : m = 4 for ΔT, internal node in the grid :E1

CONTROL VOLUME  METHOD
===========================

p1a:subst([m=4],p2);

simpsum :true;

p1b:ev(p1a,ratsimp);

rem :https://stackoverflow.com/questions/28340690/maxima-expanding-nested-summation

→ surface :B1 D1 (see : figure 1),differentation @ x=i+1/2

p2a:'diff(q[1],t,1)=-k*Δy*at('diff(T,x),x=i+1/2);

p2a1:lhs(p2a)=-k*Δy*(T[i,j]-T[i+1,j])/Δx;

→ surface :B1 A1 (see : figure 1),differentation @ y=j+1/2

p2b:'diff(q[2],t,1)=-k*Δx*at('diff(T,y),y=j+1/2);

p2b1:lhs(p2b)=-k*Δx*(T[i,j]-T[i,j+1])/Δy;

→ surface :A1 C1 (see : figure 1),differentation @ x=i-1/2

p2c:'diff(q[3],t,1)=-k*Δx*at('diff(T,x),x=i-1/2);

p2c1:lhs(p2c)=-k*Δy*(T[i,j]-T[i-1,j])/Δx;

→ surface :C1 D1 (see : figure 1),differentation @ y=j-1/2

p2d:'diff(q[4],t,1)=-k*Δx*at('diff(T,y),y=j-1/2);

p2d1:lhs(p2d)=-k*Δx*(T[i,j]-T[i,j-1])/Δy;

use  : constraint:cv

p3:[p2a1,p2b1,p2c1,p2d1];

p3a:subst(p3,p1b);

set β =Δx/Δy

p3b:lratsubst([Δx/Δy=β],p3a);

p3b1:p3b*β/k;

→ find coef,

 { T[i,j],T[i+1,j],T[i-1,j],T[i,j+1],T[i,j-1]}

p3b1a:coeff(lhs(p3b1),T[i,j]);
p3b1b:coeff(lhs(p3b1),T[i+1,j]);
p3b1c:coeff(lhs(p3b1),T[i-1,j]);
p3b1d:coeff(lhs(p3b1),T[i,j+1]);
p3b1e:coeff(lhs(p3b1),T[i,j-1]);

list: coef.

p3b2a:[p3b1a,p3b1b,p3b1c,p3b1d,p3b1e];

set .coef,

p3b2b:{p3b1a,p3b1b,p3b1c,p3b1d,p3b1e};

commands : listify,setify

fig 1 : set β =1

p3c:subst([β=1],p3b);

assume(k>0);

p3d:p3c/k;

USE : p3b2a

p3d1:subst([β=1],p3b2a);/* {T[i,j],T[i+1,j],T[i-1,j],T[i,j+1],T[i,j-1]*/

rem : see comments "p3d1"

ΔT =div(grad(T))=0, five point stencil (= molecule) for laplacian 

5-point laplacian (stencil or molecule)

How to find : T[i,j]
==============

p3e:solve(p3d,T[i,j])[1];

explain : We take four neighbours around the centralpoint T[i,j],compute the meanvalue.

================================================================================