kill(all)$

fully developed flow , steady, incompressible, v=0,version 1.0
peter.vlasschaert@gmail.com,01/04/2019

rem :website : more :modelling +advanced 69-
----------------------------------------------------------------------
( equation : history, see below.)  => 
69e) The viscosity (mu) and distance (y)between two parallel plates(laminar flow,heated)

∂(ρ*Φ)/∂t + div(ρ vec(U)*Φ) =div(Γ*grad(Φ)) + S(Φ)
------------------------------------------------------------------------------------
terminology : fluid flow
===============================================
steady term + convection term = diffusion term + Source
===============================================
vec(U) = [u,v,w]  : 3D    = >  U(u,v,w)
vec(U) = [u,v]       :2D    = >  U(u,v)
vec(U) = u             :1D   = >   U(u) = u
D : dimensional 
ρ = density ( incompressible flow , ρ = constant).
vec [u] = vector u
S : pressure gradient , nuclear source, chemical reaction
================================================

problem : 2D
============
Φ = 1  :  mass conservation
Φ = u  :  momentum conservation ' x -direction'
Φ = v  :  momentum conservation ' y -direction'
----------------------------------------------------------------------
mass conservation : 
======================
∂(ρ)/∂t + div(ρ vec(u)) = 0  ... (1)
rem : grad(constant) = 0
momentum conservation :
======================
=>  'x -direction '
∂(ρ*u)/∂t + div(ρ vec(U)*u) =div(Γ*grad(u)) + S(Φ) ...(2)
=>  'y -direction '
∂(ρ*v)/∂t + div(ρ vec(U)*v) =div(Γ*grad(v)) + S(Φ) ...(3) 
------------------------------------------------------------------------------
S : general source explain above.(Φ)

depends(U,[u,v]);

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

mass conservation : see equation (1)
--------------------------------------------------------

1e) incompressible flow , ρ = constant
================================

declare(ρ,constant);

math : (1)

lhs equation (1)
-----------------------

p1: 'diff(ρ,t) + 'diff(ρ*u,x) + 'diff(ρ*v,y) = diff(ρ,t) + diff(ρ*u,x) + diff(ρ*v,y) ;

p1a: factor(p1);

result
=====

p1b: rhs(p1a) = 0;

incompressible flow
-------------------------------

p1c:part(lhs(p1b),2) = 0;

mathematics :  'p1c' = div(U)

momentum conservation : see equation (2,3)
-------------------------------------------------------------------
Γ  = constant , u = independent of time

declare (Γ,constant);

1e) momentum conservation ' x -direction' (2)

lhs : equation (2)
==============

p2:'diff(ρ*u,t) + u*'diff(u,x) + v*'diff(u,y) =diff(ρ*u,t) + u*diff(u,x) + v*diff(u,y);

rhs : equation (2)
==============

source : ' x direction ' : pressure gradient  : Sx

depends(p,[x,y]);

p2a1:'diff(Γ*diff(u,x),x) +'diff(Γ*diff(u,y),y) +'Sx = diff(Γ*diff(u,x),x) +diff(Γ*diff(u,y),y) + Sx;

p2b1:Sx = - diff(p,x);

p2c1:subst([p2b1],rhs(p2a1));

result : momentum equation ' x-direction'
===================================

p2d1:rhs(p2) = p2c1;

2e) momentum conservation ' y -direction' (3)

lhs : equation (3)
==============

p3:'diff(ρ*v,t) + u*'diff(v,x) + v*'diff(v,y) =diff(ρ*v,t) + u*diff(v,x) + v*diff(v,y);

rhs : equation (3)
===============

source : ' y direction ' : pressure gradient  : Sy

p3a1:'diff(Γ*diff(v,x),x) +'diff(Γ*diff(v,y),y) +'Sy = diff(Γ*diff(v,x),x) +diff(Γ*diff(v,y),y) + Sy;

p3b1:Sy = - diff(p,y);

p3c1:subst([p3b1],rhs(p3a1));

result : momentum equation ' y-direction'
===================================

p3d1:rhs(p3) = p3c1;

final result : three equations.
========================

p4a:p1c;
p4b:p2d1;
p4c:p3d1;

1e) fully developed flow. ( v=0)
=========================

p5a:subst([v=0],p4a);
p5b:subst([v=0],p4b);
p5c:subst([v=0],p4c);

2e) suppose : convection term = 0
---------------------------------------------------
rem : only 'p5b'  on 'lhs' = 0

p6a:p5a;
p6b:subst([lhs(p5b) = 0],p5b);

rem : p  only function x  : no 'p5c'

solution :problem
===============

1e)  second differentation of u , 'p6a'
===============================

p7:diff(p6a,x);

2e) subst 1e) into 'p6b'
===================

p7a:subst([p7],p6b);

3e) solve for 'diff(u,y,2)
===================

p7b:solve(p7a,'diff(u,y,2))[1];

find : u
======

p8:ode2(p7b,u,y);

u=0 , y=0  (wall)
============

p8a:subst([u=0,y=0],p8);

p8a1:[rhs(p8a) =lhs(p8a)];

p8b:subst([p8a1],p8);

moving : plate : u=u_plate (velocity) , y=height =h
=========================================
rem : height of moving plate above the bottom.

p8c:subst([y=h,u=u[plate]],p8b);

final result: moving plate @ top
===========================

p9:solve(p8c,%k2);

p9a:subst(p9,p8b);

p9b:expand(p9a);

velocity profile : u = u(y, others)
=>  parabolic profile.

value source: Sx
==============

p10:diff(p,x,1) = -500;

viscosity
========

p10a:Γ = 0.02;

height of channel
==============

p10b:h=6*10^(-3);

velocity of moving plate
====================

p10c:u[plate] = 10*10^(-3);

result : graphics (find properties of the function)
========================================

p11a:subst([p10,p10a,p10b,p10c],p9b);

p11b:wxplot2d([rhs(p11a)],[y,0,6*10^(-3)]);

max parabolic profile :'diff(u,y,1) = 0;
==============================

p11c:'diff(u,y) = diff(rhs(p11a),y);

wxplot2d([rhs(p11a),rhs(p11c)],[y,0,6*10^(-3)]);

intersection : blue line 'and' red line.  ( value for extrema:(min or max))

p11d:rhs(p11c) = 0;

p11e:solve(p11d,y)[1];

p11e:float(p11e);

min or max (second order differentation)
=========
max :  diff(u,y,2) < 0
min  :  diff(u,y,2) > 0

p12:'diff(u,y,2) = diff(rhs(p11c),y);

if rhs(p12) < 0  then disp(maximum)  else disp(
 minimum)$

when : u[plate] = 0

p13::subst([p10,p10a,p10b,u[plate] = 0],p9b);

graphics : plate doesn't move , ' channel flow,laminair '.
========

p13a:wxplot2d([rhs(p13)],[y,0,6*10^(-3)]);

max:maximum.
-----------------------

p13b:diff(p13,y);

p13c:solve(rhs(p13b) =0,y);

p13d:float(p13c);

wxplot2d([rhs(p13),rhs(p13b)],[y,0,6*10^(-3)]);

intersection : blue line 'and' red line.  ( value for extrema:(max)).

info :u[plate] = 0
=> parabolic profile =  'perfect symmetric'.

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