kill(all)$

maxima : fvm :nonuniform mesh diffusion 2D with sourc S gauss divg theorema (J[q] flux,LHS  wall),version 1.0
peter.vlasschaert@gmail.com,21/22/2022

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

load(lrats)$

nonuniform mesh : diffusion 2D with sourc S,cartesian (LHS wall : C[2],I[1],C[1],D[2])
=========================================

general figure (LHS wall =C[2],I[1],C[1],D[2]),ΔV  =Δy*Δx/2

x-direction : unit vector = e[x] = e[1]
y-direction : unit vector = e[y] = e[2] → perpendicular to (x -direction)
---------------------------------------------------------------------------------------------
math:
====
dot product : → e[i]*e[j] =δ(i,j)  → (i=j → δ(i,j) =1 ,i ≠ j →  δ(i,j) =0)
---------------------------------------------------------------------------------------------
===============================
CV = Control Volume (around : p)
===============================
volume : 1D =Δx*1*1       =Δx 
volume : 2D =Δx*Δy*1    =Δx*Δy
volume : 3D =Δx*Δy*Δz = Δx*Δy*Δz 
CS  = Control Surface A[m] ,m={w,e,s,n}
===============================
faces             : w=west, e=east,s=south,n=north  ⇒ unknown
nodal values : W,E,S,N                                          ⇒     known
---------------------------------------------------------------------------------------------
x-direction : W,w,p,e,E
y-direction : S,s,p,n,N
---------------------------------------------------------------------------------------------
A[m] → m                                                           ⇒ vector
--------------------------------------------------------------------------------------------
distance : ||.||
Δx(e) =||Ep||
Δx(w)=||Wp||
Δy(s) =||Sp||
Δy(n)=||Np||
p surround with CV
....................................................................................................................

LHS (Wall J[q] flux, wall =||C[1],C[2]||)

LHS wall(C[2],n,s,C[1]), diffusion problem with source
=========================================

px:[w,e,s,n];

dot product rules ,unit(n) vectors.
-----------------------

⇒ unit vector :without index

p1:[e1*e1=1,e2*e2=1,e3*e3=1,e1*e2=0,e2*e3=0,e1*e3=0];

⇒ unit vector :with     index ,later used ,for  simplification,diffusion term "see:below".

p1a:[e[1]*e[1]=1,e[2]*e[2]=1,e[3]*e[3]=1,e[1]*e[2]=0,e[2]*e[3]=0,e[1]*e[3]=0];

simp :false$

p2a:'integrate(∇.'(Γ*∇(Φ)),V) ='integrate('(Γ*'∇(Φ)).'unit(n),A);

p2b:lhs(p2a)='∑*'(( B1)[f].(B2)[f]);

simp:true$

f ={w,e,s,n} ,faces(=CS) of CV
--------------------------------------
B1:Γ*∇(Φ)       ⇒ vector
B2:(vec(dA))   ⇒ vector
B1[f] . B2[f]    ⇒ dot product,see : definition with infix operator
--------------------------------------
vec(dA) = unit(n)*dA
-------------------------------------

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

computation : (p2b),2D f ,f1,f2
----------------------------------------

depends(Φ,[x,y]);

dependencies;

→ (Γ*∇(Φ))[f],2D,f ={w,e,s,n} 
---------------------------------------

p4a:∇(Φ) =diff(Φ,x)*e[1]+diff(Φ,y)*e[2];

p4a1:ratcoeff(rhs(p4a),e[1]) =at (diff(Φ,x), [x = f1]);
p4a2:ratcoeff(rhs(p4a),e[2]) =at (diff(Φ,y), [y= f2]);

p4a3:Γ=Γ[f];

A[f],f ={w,e,s,n},f1=[w,e],f2=[s,n]

p4bx:f=[w,e,s,n];
p4bx1:f1=makelist(rhs(p4bx)[i],i,1,2);
p4bxa:makelist(rhs(p4a1),f1,rhs(p4bx1));
p4bxa1:makelist(subst([f1=rhs(p4bx)[i]],rhs(p4a1)),i,1,2)*e[1];
p4bx2:f2=makelist(rhs(p4bx)[i],i,3,4);
p4bxb:makelist(rhs(p4a2),f2,rhs(p4bx2));
p4bxb1:makelist(subst([f2=rhs(p4bx)[i]],rhs(p4a2)),i,3,4)*e[2];

p4c:makelist(A[f],f,rhs(p4bx));

special : rem (= LHS wall,blue part)
----------------------------------------------
CV → ΔV =Δy*Δx/2
CS  → west side (Figure:2) = "C1,C2"

ΔV :Δx/2*Δy;

p4c1:[-Δy*e[1] , Δy*e[1],-Δx/2*e[2],Δx/2*e[2]];

p4c1aa:map("=",p4c,p4c1);

p4c1a:makelist(p4bxa1[i]*p4c1[i],i,[1,2]);
p4c1b:makelist(p4bxb1[i]*p4c1[i+2],i,[1,2]);

p4c1:flatten([p4c1a,p4c1b]);

p4c2:map("=",p4c,p4c1);

assume(Γ,constant)

→ part 1:diffusion

p4c3:expand(Γ*sum(p4c1[i],i,1,length(p4c1)));

p4c3a:subst(p1a,p4c3);

→ part 2:source

p4c4:barS*ΔV;

2D problem : after integration : Γ*∇.(∇(Φ))+S = 0, S = source.
---------------------------------------------------------------------------------------

→ step 1: result integration over CV ,see : Figure 1

p4c5:p4c3a+p4c4 =0;

barS = mean S value  over the CV

→ step 2: fvm ,CDS for differential with respect to x or y 

p5a1:(at('diff(Φ,y,1),y=s)) =(Φ[p]-Φ[S])/Δy(s);
p5a2:(at('diff(Φ,y,1),y=n)) =(Φ[N]-Φ[p])/Δy(n);
p5a3:(at('diff(Φ,x,1),x=w)) =-J[q]/Γ;/*lhs wall :J[q] flux */
p5a4:(at('diff(Φ,x,1),x=e)) =(Φ[E]-Φ[p])/Δx(e);
p5a5: [p5a1,p5a2,p5a3,p5a4];

p5:subst(p5a5,p4c5);

→ step 3: general equations ( with S = source.)

→ step 3a) : general equations , barS = fixed 

p6:a[p]*Φ[p] =a[W]*Φ[W] +a[E]*Φ[E] +a[S]*Φ[S] +a[N]*Φ[N] +b;

p6a1:a[p]=expand(ratcoeff(lhs(p5),Φ[p]))*(-1);
p6a2:a[W]=expand(ratcoeff(lhs(p5),Φ[W]));/*lhs wall :J[q] flux */
p6a3:a[E]=expand(ratcoeff(lhs(p5),Φ[E]));
p6a4:a[S]=expand(ratcoeff(lhs(p5),Φ[S]));
p6a5:a[N]=expand(ratcoeff(lhs(p5),Φ[N]));
p6a6:b = barS*ΔV+J[q]*Δy;

→ step 3b) : general equations , barS =S[c]+S[p]*Φ[p], linear source.

use :p5,CV =Δx*Δy ,see : Figure 1

p6b1:subst([barS=S[c]+S[p]*Φ[p]],p5);

p6aa1:a[p]=expand(ratcoeff(lhs(p6b1),Φ[p]))*(-1);
p6aa2:a[W]=expand(ratcoeff(lhs(p6b1),Φ[W]));/*lhs wall :J[q] flux */
p6aa3:a[E]=expand(ratcoeff(lhs(p6b1),Φ[E]));
p6aa4:a[S]=expand(ratcoeff(lhs(p6b1),Φ[S]));
p6aa5:a[N]=expand(ratcoeff(lhs(p6b1),Φ[N]));
p6aa6 : b =S[c] *Δx/2*Δy+J[q]*Δy;

representation : J[q], f=w, LHS Wall, Dirichlet condition
--------------------------------------------------------------------------
BC :Dirichlet condition,BC = Boundary Conidition (LHS wall)

p7:J[q] =-Γ(Φ[p]-Φ[f])/Δx(f);
p7a:subst([f=w],p7);

"→ from : step 3a"

→ change b : S=Sbar*ΔV, new :b

bnew1 :facsum(subst(p7a,p6a6));

"→ from : step 3b"

→ change b : S=Sbar*ΔV , barS =S[c]+S[p]*Φ[p], new :b

bnew2 :facsum(subst(p7a,p6aa6));

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