kill(al)$

maxima:fvm:how to use simple algorithme (example: 2). version 1.0
peter.vlasschaert@gmail.com,08/02/2023

example 2:

Given : parameters ,equations.
---------------------------------------

u[A] = 15 , velocity @ A , face of the cell
p3     = 10 ,pressure
momentum : equations 
u[B] = 5+2.5*(p1-p2)
u[C] = 5+7.5*(p2-p3)

pa:u[A] = 15 ;                     /* velocity @ A , face of the cell*/
pb:p3     = 10;                     /*pressure given : don't make correction  =0 */

pc1:u[B] = 5+2.5*(p1-p2)$
pc2:u[C] = 5+7.5*(p2-p3)$

find : p1 and p2
---------------------

→ first step : guess , pressure 

s1:p1^(°) = 0;
s2:p2^(°) = 0;
s3:[s1,s2,pb];

→ second step: use s3,pb into momentum 

s1a:subst(s3,u[B] = 5+2.5*(p1^(°)-p2^(°)));
s1b:subst(s3,u[C] = 5+7.5*(p2^(°)-p3));

→ third step: mass balance ( = continuity equation ),1D

Patankar (1980)

pk1:u[e] - u[w] = 0;

momentum equation:

pk1a:u[e] = u[e]^(°)+d[e]*(p[p]^(∄)-p[E]^(∄));
pk1b:u[w] = u[w]^(°)+d[w]*(p[W]^(∄)-p[p]^(∄));

pk1c:rhs(pk1a)=rhs(pk1b);

corrected pressure : with Pantakar notation.

pkk:solve([pk1c],p[p]^(∄));

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

fourth step : cell i = 1,2

ps1:d[A] =0;
ps2:d[B] =2.5;
ps3:d[C] =7.5;

ps:[ps1,ps2,ps3];

cell 1: around p1^(.)
-------------
see : figure 1, use pkk

psa1:flatten([w=A,e=B,u[A]^(°)=15,[ps]]);/* include fixed values*/
psa2:subst(psa1,pkk );

first ref : figure 1 ,E = p2^(.) ,u[B]^° = see :s1a

psa3:[p[p]^(∄)=p1^(∄),u[B]^(°)=rhs(s1a), p[E]^(∄)=p2^(∄)];

psa:expand(subst(psa3,psa2));

psf1:psa;

cell 2: around p2^(.)
-------------
see : figure 1, use pkk

psb1:flatten([w=B,e=C,[ps]]);/* include fixed values*/
psb2:subst(psb1,pkk );

second ref : figure 1 ,E = p2^(.) ,u[B]^° = see :s1a,s1b

psb3:[p[p]^(∄)=p2^(∄),u[B]^(°)=rhs(s1a),u[C]^(°)=rhs(s1b), p[W]^(∄)=p1^(∄)];

psb:expand(subst(psb3,psb2));

rem : given value :p[E]^∄ = 0 → p3

psf2:subst([p[E]^(∄)=0],psb);

fifth step : solve psaf1,psbf2

pr1:sort(listofvars([psf1,psf2]));

pr2a:flatten([psf1,psf2]);

pr2b:solve(pr2a,[p1^(∄),p2^(∄)])[1];

when you don't want : ratprint:false;

sixth step : check , cvg continuity equation

use : second step ,use :ps, p3 need no correct =0,pr2b

pz1:u[B]=u[B]^(°)+d[B]*(p1^(∄)-p2^(∄));
pz2:u[C]=u[C]^(°)+d[C]*(p2^(∄)-p3^(∄));

pzz:flatten([u[B]^(°)=rhs(s1a),u[C]^(°)=rhs(s1b),pr2b,ps,p3^(∄)=0]);

pzf :subst(pzz,[pz1,pz2]);

use :pk1

pzc : pzf[1]-pzf[2];

seventh step : update pressure.
use : s3

pw1:p1=p1^(°)+p1^(∄);
pw2:p2=p2^(°)+p2^(∄);
pw:[pw1,pw2];
pwf1:flatten([s3,pr2b]);
pwf2:subst(pwf1,pw);

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