kill(all)$

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

load(lrats)$

laminar flow : between two stationary plates (num = galerkin method)
peter.vlasschaert@gmail.com, 10/04/2019

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

p1:'diff(u,y,2)='diff(p,x,1)/Γ;

p1a:Γ*p1;

p1b:p1a/ρ;

subst :  'ν =Γ/ρ'

p1b1:[Γ/ρ =ν];

p1b2:solve(p1b1,Γ)[1];

p1b3:[rhs(p1b2)=lhs(p1b2)];

p1c:lratsubst(p1b1,p1b);

p1d:lhs(p1c) - rhs(p1c)  = 0;

set :Cte ='diff(p,x,1)/ρ

p1e: 'diff(p,x,1)/ρ =Cte;

p1e1:subst(p1e,p1d);

fig 1:laminar flow between two parallel stationary plates.

u= 0,wall

how to find  u[c]   : center - line velocity
=================================
u = velocity  x-direction
v = velocity   y-direction ( v = 0,fully developed flow)

Galerkin weighted residual method (numerical method)
==============================================

problem : laminar flow ,steady state. (two parallel plates)
===============================================
--------------------------------------------
BC : (B)oundary (C)ondition
--------------------------------------------
-----------------------------------------------------------
wall :u(x,-h/2) = 0 ,below  the center-line .
wall :u(x, h/2) = 0 ,above  the center-line .
------------------------------------------------------------
v= 0 :  -h/2 < y ≤ h/2
------------------------------------------------------------

1e)  trial function  (find α[i])
==============

p2:u = sum(N[i]*α[i],i,1,n);

info :from before , parabolic profile .  'see above,website'
-------------------------------------------------------------------------------------
N[i],{1,y,y^2}  = > parabolic profile.

p2a:n=2;
p2b:N[1]=1;
p2c:N[2]=y^2;

p2d:[p2a,p2b,p2c];

p2e:subst(p2d,p2),simpsum;

general equation: Galerkin weighted residual
======================================
=>  Weighting function : WW(y) =W[i]  = > N[i] values

1e)  BC condition (bottom ,top).
-----------------------------------------------

p3a: BC1 = -h/2;
p3b: BC2 = h/2;

2e) integral representation.
----------------------------------------

depends(WW,y);

p4:WW(y) = [rhs(p2b),rhs(p2c)];

p5:integrate( WW(y)*lhs(p1e1),y,rhs(p3a),rhs(p3b)) = 0;

for   i=1
========
WW(y) = 1

p4a1:lhs(p4) = part(rhs(p4),1);

p5a1:subst([p4a1,p2e],p5);

p5a2:ev(p5a1,nouns);

h : real value ≠ 0

p5a3: part(lhs(p5a2),2) = rhs(p5a2); 

for   i=2
==========
WW(y) = y^2

p4b1:lhs(p4) = part(rhs(p4),2);

p5b1:subst([p4b1,p2e],p5);

p5b2:ev(p5b1,nouns);

p5b3:num(lhs(p5b2)) = rhs(p5b2);

p5b4: part(lhs(p5b3),2) = rhs(p5b3); 

rem : both solution are the same.

how to check : 'p5a3' ? 'p5b4'

if is(lhs(p5a3) = lhs(p5b4)) = true then
    q: lhs(p5b4) = rhs(p5b4);

3e) find α[i],i=1,2. for 'u'
-----------------------------------

find :  α[2] , use 'q:'
================

p6:solve(q,α[2]);

find :  α[1] , use 'p2e'
=================

p6a:subst(p6,p2e);

α[1] find from BC, u = 0 => y = ± h/2

p6b:subst([u=0,y=h/2],p6a);

p6c:solve(p6b,α[1]);

4e) find u  use 'p1e'
-----------------------------

p7a:[rhs(p1e) = lhs(p1e)];

p7b:subst(p7a,p6a);

p7c:subst(p6c,p7b);

p7d:subst(p7a,p7c);

p7e:factor(p7d);

use packages  :  'lratsubst'

p7f:lratsubst(p1b3,p7e);

check the solution : p7f into 'p1'
====================

p8:['diff(p,x,1)=c1];

p8a:subst(p8,p7f);

p8b:subst([p8a],lhs(p1)) =subst(p8,rhs(p1));

p8c:ev(p8b,nouns);

p8d:is(p8c);

?  ' How to solve :  Max (u) = u[c] ', 
============================
see : Figure 1

profile : parabolic
===============

p9:p7e;

Max(u) = u[c]  : 'center-line velocity'
===========

p9a:[y=0,u=u[c]];
p9b:subst(p9a,p9);

p9c:lratsubst(p1b3,p9b);

We  can also use arbitrary function :
==============================
symmetric interval : cos(-y) = cos(y),along the center-line.

1e) function must be constrained on BC
------------------------------------------------------------
uu = 0

p10:uu=λ*cos(π*y/h);

check :  BC conditions

p10a:rhs(p10)=0;
p10b:solve(p10a,y);

p10c:subst([y=-h/2],p10);

2e) use'p5' ,approximation: 'galerkin method'
------------------------------------------------------------------
WW(y) = rhs(p10)  , uu = u

p11:p5;

p11a:[WW(y) = rhs(p10),u=rhs(p10)];

p11b:subst(p11a,p11);

p11c:ev(p11b,nouns);

p11d:num(lhs(p11c)) = rhs(p11c);

p11e:solve(p11d,λ);

we have two solutions 
----------------------------------

p11e1:p11e[1];
p11e2:p11e[2];

rem : second solution trivial, only the first.

3e) find u, distribution function.
-----------------------------------------------

p12:[p11e1,p11a[2]];

p12a:subst([p11e1],p12[2]);

p12b:subst(p7a,p12a);

p12c:lratsubst(p1b3,p12b);

4e) find  u[c] = max(u)
--------------------------------

p12d:[y=0,u=u[c]];

p12d1:subst(p12d,p12c);

result : two solutions , problem ( galerkin method).
----------------------------------------------------------------------------

values : h =0.1,dp/dx=-5.0*10^(-3) N/m^2,Γ =10^(-3) Ns/m^2

p13b:[h=0.1,'diff(p,x)=- 5.0*10^(-3),Γ=10^(-3)];

1e)  exact solution .
================

p13:p7f;
p13a:p9c;

numerical
=========

une:subst(p13b,p13);
unec:subst(p13b,p13a);

2e) approximation solution
=======================

p14:p12c;
p14a:p12d1;

numerical
=========

una:subst(p13b,p14);
unac:subst(p13b,p14a),numer;

error = max(u[c])= abs(unec-unac)
=============================

error_uc=abs(rhs(unec)-rhs(unac));

graph (exact,approx)
==================

wxplot2d([rhs(une),rhs(una)],[y,-0.1/2,0.1/2]);

blue graph : exact 
red               :  approximation
h =0.1 m

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