kill(all)$

maxima : fvm : example 5, 2d simple algorithm ,version 1.0
peter.vlasschaert@gmail.com,17/02/2023

load package
-----------------

load(lrats)$

ref :
lgo :(maxima)
new:10/02/2023:maxima :  fvm:simple algo.1D steady (semi-implicit method for pressure linked equation), version 1.0
download  (maxima,code)       :maximafvmsimplealgo1dsteadypressurelinkedversion1

Given : Δx = Δy
----------

→ parameter :

p1a:u[w] = 50;
p2a:v[s] = 20;
p3a:P[N] = 0; /*need no correction*/
p4a:P[E] = 10;/*need no correction*/

p4:[p1a,p2a,p3a,p4a];

→ momentum equation :

p5a:u[e]=de*(P[p]-P[E]);
p5b:v[n]=dn*(P[p]-P[N]);

p5c:[de=1,dn=.6];

p5a1:subst(p5c,p5a);
p5b1:subst(p5c,p5b);

find : u[e],v[n],p[p]
==============

suppose : values

→ step 1 :

p6a:P[W] =20;
p6b:P[S] =20;
p6c:P[p] =15;

→ step 2:

p7a:u[e]^(°) =de*(P[p]^(°)-P[E]^(°));
p7b:subst([u=v,e=n,E=N,de=dn],p7a);

p7a1:subst([de=rhs(p5c[1]),P[p]^(°)=rhs(p6c),P[E]^(°)=rhs(p4a)],p7a);
p7b1:subst([dn=rhs(p5c[2]),P[p]^(°)=rhs(p6c),P[N]^(°)=rhs(p3a)],p7b);

→ step 3:

∫∫∂u/∂x*dx*dy +∫∫∂u/∂y*dy*dy 
1e) integration,x-direction  : w → e :
2e) integration,x-direction  : s  → n :
=> 

p8:(u[e]-u[w])*Δy+(v[n]-v[s])*Δx = 0;

p8a:u[e] =u[e]^(°)+de*(P[p]^(∃)-P[E]^(∃));
p8b:u[w] =u[w]^(°)+dw*(P[W]^(∃)-P[p]^(∃));
p8c:v[n] =v[n]^(°)+dn*(P[p]^(∃)-P[N]^(∃));
p8d:v[s] =v[s]^(°)+ds*(P[S]^(∃)-P[p]^(∃));

p8e:[p8a,p8b,p8c,p8d];

→ pressure corrected equation :

p9:expand(subst(p8e,p8));

p9a:(de+dw+ds+dn)*P[p]^ (∃) = de*P[E]^(∃)+dw*P[W]^(∃)+dn*P[N]^(∃) +ds*P[S]^(∃) +(u[w]^(°)-u[e]^(°)+v[s]^(°)-v[n]^(°));

computation : P[p]^∄
----------------------------

p10:[p5c,p7a1,p7b1,P[N]^(∄)=0,P[E]^(∄)=0,dw=0,ds=0];

p10a:flatten(p10);

p10b:subst(p10a,p9a);

use : p1a,p2a

p10c:lhs(p10b) =u[w]^°+v[s]^°-14.0;

p10d:subst([u[w]^°=rhs(p1a),v[s]^°=rhs(p2a)],p10c)/1.6;

find : u[e]
-------------

use :p8a,P[E]^(∄) =0

p11a:p8a;

p11b:u[e]=1*35+5;

find : v[n]
-------------

use :p8c,P[N]^(∄) =0

p12a:p8c;
p12b:v[n]=.6*35+9;

step 5: pressure correction.

p13a:P[p] =P[p]^(∄)+P[p]^(°);
p13b:lhs(p13a) =35+15;

check cvg : b =0

p14a: b=(u[w]^(°) -u[e]^(°))+(v[s]^(°) -v[n]^(°));

see :p4,p11b,p12b

p14b :lhs(p14a) =(50-40)+(20-30);

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