kill(all)$

laminar flow in a pipe ( N-S),version 1.0
peter.vlasschaert@gmail.com,24/04/2019.

load packages ?

load(lrats)$

see: 'website'

====================================================
fully developed flow , steady, incompressible, v=0,version 1.0
peter.vlasschaert@gmail.com,01/04/2019
====================================================
=>  we start from : final result : three equations.
================================================================
'diff(v,y,1)+'diff(u,x,1)=0     'incompressible flow',condition                             (1)
('diff(u,y,1))*v+u*('diff(u,x,1))=Γ*('diff(u,y,2))+Γ*('diff(u,x,2))-'diff(p,x,1)           (2) ('diff(v,y,1))*v+u*('diff(v,x,1))=Γ*('diff(v,y,2))+Γ*('diff(v,x,2))-'diff(p,y,1)            (3)
================================================================

2D version (N-S) 
==============

p1a :'diff(v,y,1)+'diff(u,x,1)=0;
p1b:('diff(u,y,1))*v+u*('diff(u,x,1))=Γ*('diff(u,y,2))+Γ*('diff(u,x,2))-'diff(p,x,1);
p1c:('diff(v,y,1))*v+u*('diff(v,x,1))=Γ*('diff(v,y,2))+Γ*('diff(v,x,2))-'diff(p,y,1);   

depending : u = u(x,y)
                        v = v(x,y)
                        p = p(x,y)

depends ([u,v,p],[x,y]);

fully developed flow  :  v=0  (because perpendicular flow direction , 2D)
rem : y - direction , 'perpendicular'.

p2:[v=0];
p2a:diff(v,y,1)+diff(u,x,1)=0;
p2b:(diff(u,y,1))*v+u*(diff(u,x,1))=Γ*(diff(u,y,2))+Γ*(diff(u,x,2))-diff(p,x,1);
p2c:(diff(v,y,1))*v+u*(diff(v,x,1))=Γ*(diff(v,y,2))+Γ*(diff(v,x,2))-diff(p,y,1); 

depends : take action then ' can be removed.

p3a:subst(p2,p2a);
p3b:subst(p2,p2b);
p3c:subst(p2,p2c);

differentation p3a , with respect  'x'

p4:diff(p3a,x);

use : 'p3a' and 'p4'  on  'p3b' and 'p3c'

p5a:[p3a,p4];
p5b:subst(p5a,p3b);
p5c:subst(p5a,p3c);

2D : only get , momentum equation 'x -direction ', direction flow .
=> 'p5b'

p6:p5b;

rem :  y - direction flow plane ( perpendicular to flow direction = x)
-----------------------------------------------------------------------------------------------------
=>  flow plane : line  '2D'  =>  line : y -drection
=>  flow plane : area '3D' =>  area : z -direction and y -direction
------------------------------------------------------------------------------------------------------
circular pipe : area  =>  circle  :  radius = r
                                            circle :  equation ,  y^2 + z^2 = r^2
                                            direction flow        ,  x -direction (u = velocity)
------------------------------------------------------------------------------------------------------

3D version (N-S) 
==============

what will be changed : second order differential u with respect to y
=>  second order differential u with respect to  y
=>  second order  differential u with respect to y
general : Lapace equation 2D  = Δu = Δu(y,z)

1e) depends must be changed   u(x,y)   to  u(y,z)
------------------------------------------------------------------------

depends([u],[y,z]);

dependencies;

p7: Δu = diff(u,y,2)+diff(u,z,2);

2e) circle equation  ( needed for area : circle)
-------------------------------------------------------------------
rem => simple circle to origin with radius = r

p8:r^2 = y^2+z^2; 

p8a:sqrt(p8);

p8b: r=lhs(p8a);

p8c: lhs(p8b) = rhs(p8a);

p8c1: rhs(p8c) = lhs(p8c);

depends([u],[y,z]);

depends([u],r);

dependencies;

3e)  modified  the equation  from  2D N-S
-------------------------------------------------------------
=>  use laplace equation  Δu , see above.

p9: ['diff(u,y,2)= lhs(p7)];

p9a:subst(p9,p6);

p9b:subst([p7],p9a);

4e)  u= u(r) rewrite both second order differntials
-------------------------------------------------------------------------
=> change rule : calculus

depends([r],[y,z]);

dependencies;

a) start : first order differential.

p10:'diff(u,y) = subst([p8c],diff(r,y))*diff(u,r);

p10a:lhs(p10) = ev(rhs(p10),nouns);

p10b:'diff(u,y,2) =' diff(rhs(p10a),y);

p10c:lhs(p10b) = ev(rhs(p10b),nouns);

p10d:lhs(p10c) = lratsubst([p8c1],rhs(p10c));

p10e:p10d,expand;

analog : 'diff(u,z,2)

p11:[y=z];

p11a:subst(p11,p10e);

combine : Δu

p12: p10e+p11a;

p12a:lhs(p12) = lratsubst([p8c1,p8],rhs(p12));

p12b:lhs(p12) = expand(rhs(p12a));

a) coef: 'diff(u,r,2)  from  rhs(p12b)

p13m1:'diff(u,r,2);

p13a: ratcoeff(rhs(p12b),'diff(u,r,2));

simplify 'p13a'
---------------------

p13a1:subst([p8c],p13a);

p13a2:ev(p13a1,nouns);

put : 'diff(u,r,2)

p13a3:ratsimp(p13a2)*p13m1;

b) coef:'diff(u,r)  from  rhs(p12b)

p13m2:'diff(u,r,1);

p13b: ratcoeff(rhs(p12b),'diff(u,r,1));

put : 'diff(u,r,1)

p13b1:p13b*p13m2;

result: Δu
=========

p14a:lhs(p12) = p13a3+p13b1;

p14b:subst([rhs(p7)=lhs(p7)],p14a);

solve : differential equation , pipeflow.
=======================================
BC (1)   r =0   =>   du/dr = 0, symmetry concept.
BC(2)    r=R   =>   u=0 ,velocity (@wall)
=======================================
use : 'p9a' and 'p9b'  subst : p14b

p15a:p9a;
p15b:p9b;

1e) solve 'p15a' for Δu
----------------------------------

p16a:solve(p15a,Δu);

2e) subst 'p16a' into 'p14b'
---------------------------------------
=>  ODE to solve

p16b:subst(p16a,p14b);

3e) solve ODE
----------------------

assume(R>0);

p16c:ode2(p16b,u,r);

find : %k1,%k2

rem : make test : BC

use BC (2)  ,  r=R  => u=0
=====================

p17:[u=0,r=R];
p17a:subst(p17,p16c);

so : better use BC(1) first  ?
=======================
How to : use 'p16c'

1e)  differentation of u with respect to r
----------------------------------------------------------
=>  r=0 =>  du/dr = 0

p18a: diff(rhs(p16c),r,1)=0;

p18b: solve(p18a,%k1)[1];

p18c:subst([r=0,diff(u,r)=0],p18b);

2e)  result 'p16c'
-------------------------
use : 'p18c'

p18d:subst([p18c],p16c);

3e)  find  %k2
--------------------

p18e:solve(p18d,%k2);

use BC(2)  
=========
=>  r=R  => u = 0

1e) use BC(2)
---------------------

p19:[r=R,u=0];
p19a:subst(p19,p18d);

2e) solve  %k2
-----------------------

p19b:solve(p19a,%k2);

final solution : 'p15a' or 'p15b'
=========================
use : p16c

p20:subst([p18c,p19b[1]],p16c);

p20a:factor(p20);

p20b:expand(p20a);

result : u is a parabolic profile.
=====

rem : profile parabolic 'r'
====================
u = b1*r^2+b2*r+b3

b1=ratcoeff(rhs(p20b),r,2);
b2=ratcoeff(rhs(p20b),r,1);
b3=ev(rhs(p20b),r=0);

see : 'website'

fluid :simple pipeflow (use shear force),version 1.0
peter.vlasschaert,23/04/2019

check solution    r=R  result  u = 0
=============

p21:[r=R,u=0];
p21a:subst(p21,p20a);

p21b:is(p21a);

'------------------------------------------------------------------------------------------------------'