kill(all);

fluid : Continuity equation  cylindrical  coordinates,part1.
peter.vlasschaert@gmail.com,30/04/2019

load packages.
=============

load(lrats)$

elementary volume element :cylindrical

=>    (x,y,z)   <->  (r,θ,z)

What is a continuity equation :  ?

Answer : conservation mass in the system.
======================================================
(rate mass "in"  -  rate mass  "out " =  rate mass "accumulation")
======================================================
--------------------------------------------------------------------------------------------------
rem :  [Kg/s = Kg/m^3*m^3/s = Kg/m^3*m^2*m/s]
             = ρ *Q = ρ *A*v
--------------------------------------------------------------------------------------------------
symbol :  ρ  = density     [Kg/m^3]
                   Q = flow rate   [m^3/s]
                   A = Area           [m^2]
                    v = velocity     [m/s]
---------------------------------------------------------------------------------------------------
rate mass " "  =  [Kg/s]  ,rate =  [1/s]
---------------------------------------------------------------------------------------------------

General equation : continuity
========================

p1:rate[mass,in] -rate[mass,out] =  rate[mass,accum];

rem :accum: accumulation

principle ,stationary volume element : ΔV = Δx*Δy*Δz
(design with geogebra 5)

green  wall : area  =>   Δy*Δz
purple wall : area =>    Δy*Δz
distance = green wall between purple wall = Δx  ,( red - direction)
---------------------------------------------------------------------------------------------------------------------------
Δx : red-direction
Δy : green-direction
Δz : blue-direction
----------------------------------------------------------------------------------------------------------------------------

differential equation of continuity     (x-direction ,see Figure 2)
====================================================

1e) Derivation 'continuity equation'
-----------------------------------------------------

Flow into the green wall  and accumulation between the green wall and
purple wall then flow goes out through purple wall.

How to get taylor series: ?

Menu : Calculus : Get Series.

a) with two terms you have linear approximation .
---------------------------------------------------------------------------

@ point :0

p2a:taylor(f(x), x, 0, 1);

@point : 0

p2b:taylor(f(x), x, a, 1);

b)  f(x+Δx) :  ?
---------------------

b1) taylor :serie , with three terms
---------------------------------------------------

p3a:f(x+Δx) = f(x) + 'diff(f(x),x)*Δx/1! + 'diff(f(x),x,2)*(Δx)^2/2!;

b2) substitution x => x0
------------------------------------

p3b1:[f(x) = f(x0)];
p3b2:subst([x=x0],lhs(p3a)) = subst(p3b1,rhs(p3a));

b3)  set : Δx+x0 = x
-----------------------------

p3c1:Δx+x0 = x;
p3c2:solve(p3c1,Δx);

b4)  substitution  Δx=x-x0 into p3b2
-----------------------------------------------------

p4d1:[p3c1];
p4d2:subst(p4d1,lhs(p3b2)) = subst(p3c2,rhs(p3b2));

rem : see above ' p2b' relation with taylor series  from  Maxima,a =x0

b5) find the coef .of derivatives :'diff(f(x0),x,2), ('diff(f(x0),x,1))
-------------------------------------------------------------------------------------------

rhs : 'f(x)'

coef : f(x0)
----------------

p4e1:ratcoeff(rhs(p4d2),('diff(f(x0),x,0)));

coef : f'(x0)
-----------------

p4e2:coeff(rhs(p4d2),('diff(f(x0),x,1)));

coef : f''(x0)
-----------------

p4e3:coeff(rhs(p4d2),('diff(f(x0),x,2)));

b6) taylor realation between f(x+Δx)
-----------------------------------------------------

Maxima : 'command',three terms
--------------------------------------------------
above : a =x0,taylor

p4f1:taylor(f(x), x, x0, 2);

Math : 'p4d2' ,above.
-------------------------------

p4f2:p4d2;

linear approximation : math. 'p5a'
============================

p5:p3a;

set : f''(x)  = 0
----------------------

p5a:lhs(p5) = subst([('diff(f(x),x,2)) =0],rhs(p5));

relation : between p5a and (ρ*u)
                                                                          (x+Δx)
======================================

see :Figure 2:
f : ρ*u,see out of the control volume : (ρ*u)
                                                                                (x+Δx)
           , see into the control  volume    :(ρ*u)
                                                                                (x)
---------------------------------------------------------------------

what about definition : 'p1'
======================

into wall green : f(x), 'Δy*Δz'
------------------------------------------

lhs: 'p1'
------------

p6a:rate[mass,in] = f(x) *Δy*Δz;
p6b:rate[mass,out] = f(x+Δx) *Δy*Δz;

p6c:subst([p6a,p6b],lhs(p1)),factor;

use : 'p5a'

p6d:subst([p5a],p6c);

use : Δx*Δy*Δz = ΔV

doesn't works : how to solve

p6e:subst([Δx*Δy*Δz = ΔV],p6d);

solve: 'p6e'

p6f:[Δx*Δy*Δz = ΔV];
p6f1:lratsubst(p6f,p6d);

rem : because part of packages. 'lrats'

rhs:'p1'
-----------

p7:rate[mass,accum] =Δρ/Δt*ΔV;

Δρ/Δt*ΔV =[Kg/m^3*m^3/s]

p7a:rhs(p7);

relation between 'lhs' +'rhs' = 0
-----------------------------------------------

p8:p7a-p6f1= 0;

p8a:p8,factor;

Δv≠ 0,easy way to remove : ΔV =1,unit cube.

p8b:subst([ΔV=1],p8a);

p8c:p8b,expand;

limit  Δt  ->  0

p8d:'diff(ρ,t) +part(p8c,1,2) = rhs(p8c);

depends ([u],[x]);

dependencies;

from :relation : between p5a and (ρ*u).
use : 'p8d'

p8e:[f(x) = ρ*u];
p8e1:subst(p8e,p8d);

continuity equation : x-direction
-------------------------------------------------
x-direction  : u
y-direction  : v
z-direction  : w
U=U(u,v,w)

1e) general
------------------

p9:vec(U) =[u,v,w];

2e) x-direction 
----------------------

p9a:subst([v=0,w=0],p9);

conitinuity : x,y,z  : direction
------------------------------------------
because symmetric

because :rhs,same for all:diff(ρ,t)

y-direction
----------------

depends([v],[y]);

p9b:subst([u=0,w=0],p9);
p9b1:subst([u=v,x=y],part(lhs(p8e1),1));

z-direction
----------------

depends([w],[z]);

p10b:subst([u=0,v=0],p9);
p10b1:subst([u=w,x=z],part(lhs(p8e1),1));

combination : x-direction, y-direction,z-direction
=========================================

p11:lhs(p9) =rhs(p9a+p9b+p10b);

p11a:lhs(p8e1)+p9b1+p10b1 = 0;

ρ =constant in time , incompressible flow : trick ρ=1.

p11b1:subst([ρ=1],p11a);

rem : div(U) = 0. delta.U =0