kill(all)$

maxima : fvm: polair coordinates (cilindrical),grad,diffusion+S(=source) "version 1.0"
peter.vlasschaert@gmail.com,26/12/2022

ref :
--------------------------------------------------------------------------------------------------------------------------------------------------
maxima : fvm :nonuniform mesh diffusion 2D with sourc S gauss divg theorema (J[q] flux,LHS  wall),version 1.0
---------------------------------------------------------------------------------------------------------------------------------------------------

load(lrats)$

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

⇒ unit vector :without index

pd:[e[r]*e[θ]=0,e[θ]*e[θ]=1,e[r]*e[r]=1];

gauss divergence theorema
----------------------------------

simp :false$

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

p1bb:lhs(p1aa)='∑*'(( B1)[f].(B2)[f]);

f ={e,w,n,s}

distance : s = r*θ
---------------------

qa: s=r*Δθ; /*s=distance */

CV = control volume = HIJKH
CS = control surface 
w → HI → A[w]
e  → KJ → A[e]
s  → IJ  → A[s]
n → HK→ A[n]

paw : A[w] = -Δr*e[θ];
pae : A[e]  =  Δr*e[θ];
pas : A[s] =  - r[s]*Δθ*e[r];
pan : A[n] =   r[n]*Δθ*e[r];

paa:[paw,pae,pas,pan];

nabla :∇ operator "cilindrical coordinates"
----------------------------------------------------

p1:∇(q) = 'diff(q,r)*e[r] + 1/r*'diff(q,θ)*e[θ];

p1a:subst([q=Φ],p1);

∇(Φ)= ∇Φ

simp:true;

depends(Φ,[r,θ]);

declare(Γ,constant)$

diffusion part :

p1c:expand(p1a*Γ);

→ different faces :k,k={e,w,n,s}

p1ce: ∇[k](Φ) :=(Γ*('diff(Φ,θ,1))*e[θ])/r[k]+Γ*e[r]*('diff(Φ,r,1));

use :p1bb

→ B1[f]

p1d1:lhs(p1c);
p1d1a:Γ*∇[f](Φ);

→ B2[f] =A[f]

final : result dot product
-----------------------------------

p1e:p1d1a*A[f];

p1e1a:makelist(subst([f=m],p1e),m,[w,e,n,s]);

pbw : subst(paa,p1e1a[1]);
pbe : subst(paa,p1e1a[2]);
pbs : subst(paa,p1e1a[4]);
pbn : subst(paa,p1e1a[3]);

→ dot rule : pd

pcw : lratsubst(pd,pbw)/Γ;
pce : lratsubst(pd,pbe)/Γ;
pcs : lratsubst(pd,pbs)/Γ;
pcn : lratsubst(pd,pbn)/Γ;

pcc:[pcw,pce,pcs,pcn];

→ discritezation,see :figure 2, nomenclatuur

pdd:[('diff(Φ,θ,1))=(Φ[p]-Φ[W])/Δθ(w),
    ('diff(Φ,θ,1))=(Φ[E]-Φ[p])/Δθ(e),('diff(Φ,r,1))=(Φ[p]-Φ[S])/Δr(s)
    ,('diff(Φ,r,1))=(Φ[N]-Φ[p])/Δr(n)];

pee:makelist(lratsubst([pdd[kk]],pcc[kk]),kk,1,length(pdd));

→ sum*(B1[f] . B2[f]),see above :

diffusion part :

ppd:apply("+",pee); 

source    part :

pps:Smean*ΔV;

final result :model : diffusion+source =0, result : after the integration.
----------------------------------------------------------------------------------------

pdsi:ppd+pps=0;

How to find : ΔV , Smean = S[c]+S[p]*Φ[p]
-------------------------------------------------------

find ΔV = ?
---------------

1e) step :area triangle 
---------------------------

s1a: area[tri] = 1/2*h*b;

2e) find area : ADEA
--------------------------

s2a:As[tri]=1/2*r*(r*Δθ);

3e) find area :ABCA
--------------------------

s2b:Ab[tri]=subst([r=r+Δr],1/2*r*(r*Δθ));

s2b1:expand(s2b);

approx :
----------

sb1a:subst([part(rhs(s2b1),1)=0],s2b1);

4e) find : ΔV
-----------------

s3:ΔV = rhs(sb1a) - rhs(s2a);

s3a:subst([r=r[p]],s3);

find :Smean*ΔV = ?
-------------------------

s4a:Smean = S[c]+S[p]*Φ[p]; 

s4b:s3a;

s4c:s4a*s4b;

final result : after discretization
---------------------------------------

use :pdsi

fn :subst([s4c],pdsi);

p5:fn;

general equation :
=============

→ equation 

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

→  {a[p],b,a[W],a[S],a[N],a[E]}

p6a1:a[p]=expand(ratcoeff(lhs(p5),Φ[p]))*(-1);
p6a2:a[W]=expand(ratcoeff(lhs(p5),Φ[W]));
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=r[p]*(S[c])*Δr*Δθ;

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