kill(all)$

maxima:fvm:example 4:1d,more cells,use simple algo,version 1.0
peter.vlasschaert@gmail.com,16/02/2023

load(packages)
-------------------

load(lrats)$

ref : algo :(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

example 4:
-------------

F : definition (A=constant) :
---------------------------------

FW[jj]:=(ρ[jj]*u[jj]*A);
FE[jj]:=(ρ[jj+1]*u[jj+1]*A);

a:definition :
----------------

aW[jj]:=(ρ[jj]*d[jj]*A);
aE[jj]:= (ρ[jj+1]*d[jj+1]*A);

 principal : simple algo
------------------------------

q1a:u^(∃) =d*(p[i]^(∃)-p[i+1]^(∃));

q1b:u =u^(°)+u^(∃);

→ subst :

q:subst([q1a],q1b);

→ pressure correction equation ,1d

qq:ap*P[p]^(∃)=aw*P[W]^(∃) +ae*P[E]^(∃) +b^(∃);

pp:subst([P=pA],qq);/*general : correction eq*/
pp:subst([p=jj+1,W=jj,E=jj+2],pp);

values :

qa:subst([u=X,i=1],q);

Given : pipe flow
---------------------

parameter: general
==============

→ d:1

r1:makelist(d[i]=1,i,1,5);/*assume = constant*/

→ ρ:1

r2:makelist(ρ[i]=1,i,1,5);

rr:flatten([r1,r2]);

→ boundary condition :

r2a:u[1] = 10; /*[m/s]*/
r2b:pA[4]^∃= 0;      /*pressure =[N/m^2]*/

Solve :
--------

→ guess :

r2c:[u[2]=8.0,u[3]=11.0,u[4]=7.0]$

internal nodes    :Ai, i= 2,3
-----------------------------------

→ A[2]

s2a:aw=subst(rr,aW[1]);
s2b:ae=subst(rr,aE[1]);
s2:ap=facsum(rhs(s2a)+rhs(s2b));

→ use : guess value."r2c"

s2c:b^(∃) =subst(rr,FW[2]-FE[2]);
s2c:subst(r2c,s2c);

→ use:pressure correction equation: pp

vA2:subst([jj=1],pp);
vA2:subst([s2a,s2b,s2c,s2],vA2);
vA2:subst([r2a],vA2);

→ A[3]

s3a:aw=subst(rr,aW[2]);
s3b:ae=subst(rr,aE[2]);
s3:ap=facsum(rhs(s3a)+rhs(s3b));

→ use : guess value."r2c"

s3c:b^(∃) =subst(rr,FW[3]-FE[3]);
s3c:subst(r2c,s3c);

→ use:pressure correction equation: pp

vA3:subst([jj=2],pp);
vA3:subst([s3a,s3b,s3c,s3],vA3);

boundary nodes :Ai,i=1,4
-------------------------------

→ A[1]

s1a:aw=subst(rr,aW[1]);
s1b:ae=subst(rr,aE[1]);
s1:ap=facsum(rhs(s1a)+rhs(s1b));

corrected : bc

→ use : guess value."r2c"

s1a:aw =0;
s1:ap=facsum(rhs(s1a)+rhs(s1b));

reason :u[1],given @begin pipe.

corrected : FW[1] =0 ,why outside pipe. (use : bc condition)

s1c:b^(∃) =subst(rr,-FE[1]+u[1]*ρ[1]*A);/*u[1] = from boundary, not corrected*/
s1c:subst(r2c,s1c);
s1c:subst([r2a],s1c);

→ use:pressure correction equation: pp

vA1:subst([jj=0],pp);
vA1:subst([s1a,s1b,s1c,s1],vA1);

→ A[4]

→ use:pressure correction equation: pp

vA4:r2b;

use : bc,r2b

Solve : vAi,for corrected pressure, i=1,2,3,4

→ equations :

w1: makelist(sconcat(vA,i),i,1,4);

stringp(w1[1]);

eval_string(w1[1]);

vv1:ratsimp(makelist(eval_string(w1[i]),i,1,length(w1))/A);

vv1:subst([pA[4]^(∃) =0],vv1);

vv1a:makelist(vv1[i],i,1,length(w1)-1);

→ variables :

vv2:sort(listofvars(vv1a));

vv2a:makelist(vv2[i]^(∈)=v[i],i,1,length(w1)-1);

→ subst

vv3:subst(vv2a,vv1);

use : vv3 can't use solve

eqs :[p1=p2+2,2*p2=p3+p1-3,2*p3=p2+4];

var1:listofvars(eqs);

sol:solve(eqs,var1)[1];

info:

var2:makelist( p[i]^(∃) =rhs(sol[i]),i,1,length(sol));

use : q

vv4:q;

vv4a:subst([u=u[i+1]],vv4);

final solution :

v5:makelist(subst([i=j],vv4a),j,1,3);

v5a:subst(var2,v5);

→ use : bc,d = given

v5b:subst([p[4]^(∃)=0,d=1],v5a);

input : guess value

v5c:r2c;

v5c1:u[2]^(°) =rhs(v5c[1]);
v5c2:u[3]^(°) =rhs(v5c[2]);
v5c3:u[4]^(°) =rhs(v5c[3]);

v5d :[v5c1,v5c2,v5c3];

v5dx:makelist(u[i]^(°) =rhs(v5c[i-1]),i,[2,3,4]);

is(v5d=v5dx);

final solution : velocities
-------------------------------

v6:subst(v5d,v5b);

v6a:flatten([r2a,v6]);

→ velocities : are same @ every intersection.

=> result :single iteration ,the exact velocities.

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