kill(all)$

fluid flow gravity (down ,inclined plane),incompressible flow (NS).'version1.0'
peter.vlasschaert@gmail.com, 23/04/2019

NS : Navier -Stokes.

load packages

load(lrats)$

rem:

∂(ρ*Φ)/∂t + div(ρ vec(U)*Φ) =div(Γ*grad(Φ)) + S(Φ) ' 1'
problem : 2D
============
Φ = 1  :  mass conservation
Φ = u  :  momentum conservation ' x -direction'
Φ = v  :  momentum conservation ' y -direction'
----------------------------------------------------------------------
info : incompressible flow , ρ =constant

rem : reference document (see website)

fully developed flow , steady, incompressible, v=0,version 1.0
use : p4a,p4b,p4c, correct with  'see '1',above'.

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

p1a :  Φ = 1
p1b :  Φ = u
p1c :  Φ = v 

because ρ =constant, multiply both sides. 
p1a,p1b,p1c by ρ

S(Φ)  =>  Sx,Sy from  p1b ,p1c
=========================

1e)  How to find , Γ = 0 , u=0, v=0 into p1b,p1c.
--------------------------------------------------------------------

p2:[Γ= 0, v= 0, u= 0];

p2a:Sx=subst(p2,rhs(p1b));
p2a1:[rhs(p2a) = lhs(p2a)];
p2b:Sy=subst(p2,rhs(p1c));
p2b1:[rhs(p2b) = lhs(p2b)];

2e)  how to put Sx,Sy into 'p1b,p1c'.
-----------------------------------------------------

p3a:subst(p2a1,p1b);
p3b:subst(p2b1,p1c);

3e) correction ,gravitation of inclined plane (angle = α).
-----------------------------------------------------------------
S(Φ) extra factor , to include inclined plane.
see : figure ' vec(g),components x,y direction'

p4:G=[g*ρ*sin(α),-g*ρ*cos(α)];

see : Figure 1 ,below .

p4a:lhs(p2a) = rhs(p2a)+rhs(p4)[1];
p4b:lhs(p2b) = rhs(p2b)+rhs(p4)[2];

Modified equations 'p1a,p1b,p1c'
===========================

1e)  use 'p3a,p3b' put in ' p4a,p4b'
----------------------------------------------------

p5:[p4a,p4b];

p5a:subst(p5,p3a);
p5b:subst(p5,p3b);

2e) Modified equations  
-----------------------------------

p6a:p1a;
p6b:p5a;
p6c:p5b;

ρ = constant ,  use incompressible flood.

3e) Modified equations  
-----------------------------------

a) functional relation  between p, u,v,x,y
-----------------------------------------------------------

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

depends(p,[x,y]);

declare(ρ,constant);

b) use  functional relation into 'p6a,p6b,pc'
----------------------------------------------------------------

p7a:ev(p6a,nouns);
p7b:ev(p6b,nouns);
p7c:ev(p6c,nouns);

c) final form after dividing by ρ
----------------------------------------------

p8a:p7a/ρ,ratsimp;
p8b:p7b/ρ,fullratsimp;
p8c:p7c/ρ,fullratsimp;

fully developed flow , v=0
=====================

1e)  v=0
-----------

p9a:subst([v=0],p8a);
p9b:subst([v=0],p8b),expand;
p9c:subst([v=0],p8c),expand;

2e) put 'p9a' into 'p9b' = condition see above.
--------------------------------------------------------------------
result : two equations.

a) elimination first order differential ;'diff(u,x,1)=0
-------------------------------------------------------------------------

p10:[p9a];
p10a:subst(p10,p9b);
p10b:subst(p10,p9c);

b) elimination second order differential ;'diff(u,x,2)=0
-------------------------------------------------------------------------------

take : differential from p9a to x
----------------------------------------------

p11:[diff(p9a,x)];

subst: p11 into 'p10a,p10b'
-----------------------------------------

p11a:subst(p11,p10a);
p11b:subst(p11,p10b);

subst:Γ/ρ =ν  into 'p11a,p11b'
--------------------------------------------

p12:[Γ/ρ =ν ];

p12a:subst(p12,p11a);
p12b:subst(p12,p11b);

rem : p12a,p12b , ODEs,
Maxima : ode2 to solve.
===================

first : solve p12b 'ODE'
====================

1e) solve  ode2
-----------------------

p13a:ode2(p12b,p,y);

2e) how to find  %c. ( must be function of x )
------------------------------------------------------------------

a ) constant function of x .
======================

p13b: [ %c = f(x) ];

b) from BC on top of free surface.
============================
y=height = h above bottom of the inclined channel, p=p0
p = pressure.

p13c:subst(p13b,p13a);

p13c1:[y=h,p=p0];

p13c2:subst(p13c1,p13c);

c) find f(x)
========

p13d:solve(p13c2,f(x));

d) final result  'p12a'
=================
from p13d into p13c

declare(g,constant);

p13e:subst(p13d,p13c);

p13f:p13e-p0;

p13g:lhs(p13f) =factor( rhs(p13f));

What is tangential stress:τ
======================

definition :τ  = force /surface = newton /m^2
-----------------
https://en.wikipedia.org/wiki/Shear_stress

fully developed flow : definition  (v=0)
u :  x-direction
v :  y-direction

p14:τ=Γ*diff(u,y);

Γ= dynamic viscosity 

second: solve p12b 'ODE'
======================

from before 'p13g' pressure function of y only .

a) solve p12b to take pressure gradient  x - direction =0
-----------------------------------------------------------------------------------
so equation  'p12a' be changed .

p15:[diff(p,x)=0];
p15a:subst(p15,p12a);

b) solved u(y) for p15a
-----------------------------------

p15b:ode2(p15a,u,y);

c) how to find %k1 and %k2 from 'p15b'
------------------------------------------------------------
BC : Boundary condition
bottom : inclined channel  y=0
top         : inclined channel y=h ( = free surface)
info : tangential stress concept for the free surface

  BC:bottom  (u=0 when y = 0)
=========================
use 'p15b'

p16a:[u=0,y=0];
p16b:subst(p16a,p15b);
p16b1:rhs(p16b) = lhs(p16b);

p16c:subst(p16b1,p15b);

  BC:top  (τ=0 when y = h)
=========================
use 'p15b','p14'

rem : tangential stress = 0 at the free surface.(y=h)

p17:Γ*diff(u,y) =0;

a)  use 'p17'  on  p16c
---------------------------------

p17a:diff(p16c,y);

p17b:p17a*Γ;

p17c:[lhs(p17b) =0,y=h];
p17d:subst(p17c,p17b);

b) find  % k2 for p15b
--------------------------------

p18:solve(p17d,%k2)[1];

c) final result 'p12a'
-----------------------------
use 'p16b1','p18' into 'p15b'

p19:subst([p16b1,p18],p15b),factor;

p19a:p19,expand;

rem : 'p19a' parabolic profile  for u
----------------------------------------------------
u = a1*y^2+a2*y+a3
----------------------------------------------------

a1=coeff(rhs(p19a),y,2);
a2=coeff(rhs(p19a),y,1);
a3=ev(rhs(p19a),y=0);

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

formulation  problem, inclined plane = α
=================================

problem :  
========

0 = g*sin(α)+(('diff(u,y,2))*Γ)/ρ-'diff(p,x,1)/ρ  (1)
0 =-g*cos(α)-'diff(p,y,1)/ρ                                  (2)


===========================================
BC (1) , (y=0, u=0; y = h  , τ = 0)  ,second order diffeq.
=>  %k1, %k2 ,see above.
==========================================
BC (2) , ( y=h  ,p = p0)                   , first  order diffeq.
=>  %c =f(x) ,see above.
==========================================

info :
====
u = velocity (x - direction)
Γ =dynamic viscosity 
p = pressure
p0 = surface pressure (y= h)
g   = 9.81 m/s^2 
α   = angle of inclined plane
h  = height of the channel = distance (K,L)
ρ = density of the fluid.

inclined plane α,  

solution  for  (1),(2)
================

p20a:p19;
p20b:p13g;

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