kill(all)$

maxima:fvm:simple algo.1D steady (semi-implicit method for pressure linked equation), version 1.0
peter.vlasschaert@gmail.com,10/02/2023

how to use :
new:10/02/2023:maxima   : fvm : example 3:1d momentum equation, 1d continuity equation,pipeflow,version 1.0
youtube : https://www.youtube.com/watch?v=N85HmIIgm4s&t=1333s
download  (maxima,code)       :maximafvmexample3moment1dcont1dversion12023
download  (maxima,png  )       :maximafvmexample3moment1dcont1dversion12023outp
new:10/02/2023:python     : fvm : example 3:1d momentum equation, 1d continuity equation,pipeflow,version 1.0
download  (python,code)        :pythonfvmexample3contmomentum1dpipeflowversion12023
download  (python,outp )       :pythonfvmexample3contmomentum1dpipeflowversion12023outp
new:10/02/2023:python     : fvm : example 3:1d momentum equation, 1d continuity equation,pipeflow,tabulate,version 1.0
download  (python,code)        :pythonfvmexample3contmomentum1dpipeflowversion12023tabulate
download  (python,outp  )      :pythonfvmexample3contmomentum1dpipeflowversion12023tabulateoutp

load packages,option.
---------------------------

load(lrats)$

ref : patankar notation used.(1980)
CV = Control Volume.

pressure : -∫∫∫∇p*dV 
pos : x-direction ,see figure
----------------------------------
FRAME :
========================================================
=>  1D : -∫∫∫∂p/∂x*dV = -∫∫∫∂p/∂x*dx*dy*dz =-Δy*Δz*∫∂p/∂x*dx = -Δy*Δz*p[x]
integration : ref u[e] , p → E   ,figure
result 1:-Δy*Δz*(p[E]-p[p]) =Δy*Δz*(p[p]-p[E])    ,explain : qi  = > force 
integration : ref u[w] , W → p ,figure
result 2:-Δy*Δz*(p[p]-p[W]) =Δy*Δz*(p[W]-p[p]), explain :qi  => force
=========================================================
rem : p[p] : p:pressure ,[p], p inside centroid of CV (patankar notation:1980)

qi:' integrate(f(x),x,a,b) =F(b)-F(a);

→ simple integration : cos(x), explain :qi

qi1:subst([f(x)=cos(x)],qi);

qi2:lhs(qi1)=lhs(ev(qi1,integrate));

general momentum 1D  , see above "FRAME"
-----------------------------

1e) exact equations :
------------------------

equation : velocity
-----------------------

ref :u[w] , see : FRAME,result :2

q1a:a[w]*u[w] =sum(a[nb]*u[nb],nb,1,n) + b +A[w]*(p[W]-p[p]);

q1a1:q1a/a[w];

set :d[w] =A[w]/a[w]

q1a2:expand(lratsubst([reverse(d[w] =A[w]/a[w])],q1a1));

q1a2:lhs(q1a2) =collectterms(rhs(q1a2),d[w],1/a[w]);

ref :u[e]   , see : FRAME,result :1

q1b:a[e]*u[e] =sum(a[nb]*u[nb],nb,1,n) + b +A[e]*(p[p]-p[E]);

q1b1:q1b/a[e];

set :d[e] =A[e]/a[e]

q1b2:expand(lratsubst([reverse(d[e] =A[e]/a[e])],q1b1));

q1b2:lhs(q1b2) =collectterms(rhs(q1b2),d[e],1/a[e]);

2e) starred equations : guess value contain,°
-----------------------------------------------------

q2a:[u[w]=u[w]^°,u[nb]=u[nb]^°,p[p]=p[p]^°,p[W]=p[W]^°];

q2aa:makelist(lhs(q2a[i]),i,1,length(q2a));

ref :u[w] , see : FRAME,result :2 → u^(°)[w]
use :q1a

q2a1:lratsubst(q2a,q1a);

q2a1:lhs(q2a1) =collectterms(rhs(q2a1),A[w]);

q2b:subst([w=e,W=E],q2a);

q2bb:makelist(lhs(q2b[i]),i,1,length(q2b));

ref :u[e] , see : FRAME,result :1 → u^(°)[e]
use :q1b

q2b1:lratsubst(q2b,q1b);

q2b1:lhs(q2b1) =collectterms(rhs(q2b1),A[e]);

3e) prime equations : corrected values,∄
-------------------------------------------------
u[]^∄ = corrected velocity
p[]^∄ = corrected pressure

q3x:(X^(∄) =X-X^°);

simple algo :(sum(a[nb]*u[nb]^°,nb,1,n)) =0,sum(a[nb]*u[nb],nb,1,n) =0

q3i:[(sum(a[nb]*u[nb]^°,nb,1,n)) =0,sum(a[nb]*u[nb],nb,1,n) =0];

→ ref :u[w] →u[w]^∄
use : q3 ,q1a-q2a1

q3a:q1a-q2a1;

q3aa:subst(q3i,q3a);

q3aa:factorsum(lratsubst([u[w]-u[w]^°=u[w]^∄,p[W]-p[W]^°=p[W]^∄,p[p]-p[p]^°=p[p]^∄],q3aa));

q3aa:collectterms(q3aa/a[w],A[w]);

set : A[w]/a[w] = d[w]

q3aa:factorsum(lratsubst([A[w]/a[w] = d[w]],q3aa));

→ ref :u[e] →u[e]^∄
use : q3 ,q1b-q2b1

q3b:q1b-q2b1;

q3bb:(subst(q3i,q3b));

q3bb:factorsum(lratsubst([u[e]-u[e]^°=u[e]^∄,p[E]-p[E]^°=p[E]^∄,p[p]-p[p]^°=p[p]^∄],q3bb));

q3bb:q3bb/a[e];

set : A[e]/a[e] = d[e]

q3bb:factorsum(lratsubst([A[e]/a[e] = d[e]],q3bb));

→ first step :q3x

q4:q3x;

q4a:reverse(q4+X^°);

→ ref :u[w]

q4a1:subst([X=u[w]],q4a);

use : q3aa

q4a2:subst([q3aa],q4a1);

→ ref :u[e]

q4b1:subst([X=u[e]],q4a);

use : q3bb

q4b2:subst([q3bb],q4b1);

→ second step :continuity equation

FRAME : 1
3D :∂(ρ*u)/∂x+∂(ρ*v)/∂y+∂(ρ*w)/∂z
u,v,z : velocity x-direction,y-direction,z-direction
integration : CV,dV=dx*dy*dz
--------------------------------------
∫∫∫(∂(ρ*u)/∂x*dV+∂(ρ*v)/∂y*dV+∂(ρ*w)/∂z*dV)
Patankar notation:
x-integration : w→e, (w=west,e=east )
y-integration : s→n , (s=south,n=north )
z-integration : b→t  , (b=bottom,t=top )
-------------------------------------------------------------------------------
1D :∂(ρ*u)/∂x*dV →∫∫∫∂(ρ*u)/∂x*dx*dy*dz =ρ[x]*u[x]*Δy*Δz
x-integration : w→e, (w=west,e=east )
∫∫∫∂(ρ*u)/∂x*dx*dy*dz = ρ[e]*u[e]*Δy*Δz - ρ[w]*u[w]*Δy*Δz =0
------------------------------------------------------------------------------

q5:ρ[e]*u[e]*Δy*Δz - ρ[w]*u[w]*Δy*Δz =0;

→ third step :check continuity equation 

a) updated : velocity
b) updated : pressure
use : q5 & {q4a2,q4b2}

q5a:q4a2;
q5b:q4b2;

q55:[q5a,q5b];

input : {q55} in q5

q5c:subst(q55,q5);

q5c1:expand(q5c);

→ fourth step :general equation, "prime equation : "

q6:a[p]*p[p]^∄ =a[E]*p[E]^∄+a[W]*p[W]^∄+b;

moment : b=0 ,then we find : pressure ,velocity

q6a:collectterms(lhs(q5c1),p[p]^∄,p[E]^∄,p[W]^∄) =0;

q6a1:part(lhs(q6a),1);

q6a2:reverse(q6a1-q6a);

find :p[p]^∄
---------------

q6a3:q6a2/part(lhs(q6a2),2);

find :{a[p],a[W],a[E]}
--------------------------

v1:a[p] =factorsum( part(lhs(q6a2),2));
v2:b = factorsum(subst([p[W]^∄ = 0,p[E]^∄ =0],rhs(q6a2)));
v3:a[W] = ratcoeff(rhs(q6a2),p[W]^∄);
v4:a[E] = ratcoeff(rhs(q6a2),p[E]^∄);

v5:is((rhs(a[p])) =(rhs(a[W]+a[E])));

v5a:a[p]=a[W]+a[E];

when : b=0, follows what
-------------------------------

w1:subst([b=0],v2);

rem :Δy*Δz ≠0

w1b:part(rhs(w1),1) =0;/*continuity,b=0,find solution. */

when : b=0

use :vi ,i=1,2,3,4 into q6,

w2b:lratsubst([b=0],q6);

w3b:solve(w2b,p[p]^∄)[1];

use : v5a

w4b:subst([v5a],w3b);

what about : boundary 
-----------------------------
velocity define  :d[.] =0
pressure define :p[]^∄ =0

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