kill(all)$

maxima : fvm : example 3:1d momentum equation, 1d continuity equation,pipeflow,version 1.0
peter.vlasschaert@gmail.com,09/02/2023

load packages
----------------------

load(lrats)$

see also : python ,see website

CV ={E[1],F[1],G[1],D[1]} = Control Volume , used mass balance .

Given : mathematical model : pipeflow "with two cells"
------------------------------------------------------------------------

1e) continuity equation : mass balance.
...........................................................................

q1:'diff(u(x)*ρ*A,x) = 0;

ρ = density
A = area
u = velocity (x-direction),1d

polygon ( = CV),in E[1] D[1] and out F[1] G[1]

q1a:integrate(q1, x,A,B);

set : F =u*ρ*A

q1b1:F[A] = A(A)*u(A)*ρ;
q1b2::F[B] = A(B)*u(B)*ρ;

q1b:F[A]-F[B] =0;

2e) momentum equation.
---------------------------------

q2:u(x)*ρ*A(x) =F(x);

q2a:'diff(ρ*u(x)*A(x)*u(x),x)=-A(x)*diff(p(x),x);

q2b:lratsubst([q2],q2a);

cell 1: integration x1 to x2
-------

q2c1:integrate(lhs(q2b),x,x1,x2)=-A(A)*integrate(part(rhs(q2b),1,2),x,x1,x2);

set : F(x1) =0, u(x2) =u(A),

assume F(2)≥ 0

q2c1a:subst([F(x1) =0, F(x2)=F2,u(x2) =u[A]],q2c1);

→ different : notation

q2c2a:subst([F(x2)=F2,p(x1)=p1,p(x2)=p2],q2c1a);

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

cell 2:integration x2 to x3
--------

q2c2:psubst([x1=x2,x2=x3,A(A)=A(B)],q2c1);

set :  u(x2) =u(A),u(x3) =u(B)

q2c1b:subst([u(x2) =u[A], u(x3) =u[B]],q2c2);

→ different : notation

assume F(3)≥ 0

q2c2b:subst([F(x3)=F3,F(x2)=F2,p(x3)=p3,p(x2)=p2],q2c1b);

rem : psubst ,parallel subst

→ first -step :

result : two momentum equations.
--------------------------------------------

→ general equation

q3a:q2c2a;
q3b:q2c2b;

change : q3b use with simple algorithm.

q3b1:q3b-part(lhs(q3b),2);

→ °equation : guess values ,° ( books = starred)

q4:[u[A]=u[A]^(°),u[B]=u[B]^(°),p2=p2^(°),p3=p3^(°),p1=p1^(°)];

q4a:subst(q4,q3a);
q4b:subst(q4,q3b1);

→∄ equation: see q5, ∄ (books = prime)

q5:X = X^(°)+X^(∄);
q5a:solve(q5,X^(∄))[1];

make : prime equations , fixed values 

1e) suppose :{p1,p3} is known         = 0
2e) when use simple algo : F2         = 0 in q4b

q5b:[u[A]^(°)=u[A]^(∄),u[B]^(°)=u[B]^(∄),p2^°=p2^(∄),p3^°=p3^(∄),p1^(°)=p1^(∄)];

q5b1:subst(q5b,q4a);
q5b2:subst(q5b,q4b);

use : 1e & 2e

q5c1:subst([p1^(∄)=0],q5b1)/F2;
q5c2:subst([p3^(∄)=0,F2=0],q5b2)/F3;

set : d[A],d[B]

q5d1:lratsubst([A(A)/F2=d[A]],q5c1);
q5d2:lratsubst([A(B)/F3=d[B]],q5c2);

→ second step ;

from continuity equation.

use : q1b, 

q6:q1b;

q6a:subst([F[A]=F[A]^(∄),F[B]=F[B]^(∄)],lhs(q6)) =subst([F[A]=F[A]^(°),F[B]=F[B]^(°)],lhs(q6));

assume :ρ,constant =1

q6b1:u[A]^(∄)*A(A)-u[B]^(∄)*A(B) =rhs(q6a);

q6b2:subst([F[A]^(°) =u[A]^(°)*A(A),F[B]^(°) =u[B]^(°)*A(B)],q6b1);

How find :p2^(∄)
----------------------

use :q6b2  and {q5b1,q5b2}

1e) subst :

q7a:subst([(q5d1),(q5d2)],q6b2);

2e) solve

q7b:solve(q7a,p2^(∄));

How find :u[A] & u[B]
------------------------------

use: q5

q8a:q5;

q8b1:subst([X=u[A]],q8a);
q8b2:subst([X=u[B]],q8a);

use :{q5d1,q5d2}

q8c1:subst([q5d1],q8b1);
q8c2:subst([q5d2],q8b2);

info :relaxation

use under relaxation :α[p]
---------------------------------

q9a:p2=p2^(°)+α[p]*p2^(∄);

use under relaxation :α[u]
---------------------------------

cell 1: use :q4a
-------------

q10a1:q4a;

q10b1:lhs(q10a1)/α[u] =rhs(q10a1) +(1-α[u])/α[u]*lhs(q10a1)^(°);

cell 2: use :q4b
-------------

q10a2:q4b;

q10b2:lhs(q10a2)/α[u] =rhs(q10a2) +(1-α[u])/α[u]*lhs(q10a2)^(°);

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