kill(all)$

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

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

load(lrats)$

numumerical :'diff(u,y,2)='diff(p,x,1)/Γ'
================================
=>  numerical : finite differences

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

representation :second order derivative  (grid)
=============================================

website :see : SPECIAL TOPICS  :  PDE  :    numerical
---------------------------------------------------------------------------------
u(x,y)   =  u[i,j]

ii:1;
jj:2;

u[ii,jj];

def :   diff(u,y,2)= Uyy
==================
uniform grid : Δx = Δy  ( = cartesian grid)

p2: Uyy = (u[i,j+1] -2*u[i,j]+u[i,j-1])/Δy^2;

diff(p,x,1) = pressure gradient  constant
Γ                = constant
=================================

p2a: rhs(p2) = rhs(p1);

solve : u[i,j]
==========

p3:solve(p2a,u[i,j])[1];

u[i,j]  : define  =>  function
======================

p4:define(U(i,j),rhs(p3)-lhs(p3)=0),expand;

rem : 

How define domain    j = 1..n  (y - direction)
====================================
-----------------------------------------------------------------
interior : 2,3,4,...., n-1
outside : 1,n
------------------------------------------------------------------
result : n equations 
------------------------------------------------------------------

n: 7;
m : n -2;

1e)  interior nodes domain  ( m = nodes)
======================
m=2,3,..(n-1)

p5a:U(i,2);
p5b:U(i,3);
p5c:U(i,4);
p5d:U(i,5);
p5e:U(i,6);

p6a:[p5a,p5b,p5c,p5d,p5e]$

p6b:[u[i,1],u[i,2],u[i,3],u[i,4],u[i,5],u[i,6],u[i,7]];

coef : matrix
===========

p7:coefmatrix(p6a,p6b);

how to find coef :u[i,1],u[i,7] , 'given values'
====================================

coef : u[i,1]    (bottom wall)
======================
u[i,1] : given value on the bottom wall

p7a:diff(lhs(p5a),u[i,1]);

rem :  value first column of 'p7'

result : lhs 'final system' ( = first equation)
==================================

p7a1:p7a*u[i,1];

coef : u[i,7]    (top wall)
======================
u[i,7] : given value on the top wall

p7b:diff(lhs(p5e),u[i,n]);

rem :  value last column of 'p7'

result : lhs 'final system' ( = last equation)
==================================

p7b1:p7b*u[i,7];

how to find : constant values 
=========================
=>    ' p5a,p5b,p5c,p5d,p5e '

p8:[0,0,0,0,0,0,0];

p8a:map( "=",p6b,p8);

rem : how to use  'map' command : Maxima.

p9a:subst(p8a,lhs(p5a));
p9b:subst(p8a,lhs(p5b));
p9c:subst(p8a,lhs(p5c));
p9d:subst(p8a,lhs(p5d));
p9e:subst(p8a,lhs(p5e));

vector : p9a,p9b,p9c,p9d,p9e.

p9:[p9a,p9b,p9c,p9d,p9e];

how to build 'final system' to solve : problem above.
============================================

system :  A uu = L     'final sytem'
==========================

1e) How to find 'A' from 'p7'.
-----------------------------------------

solution : to removes first column (=1) and last column(=7)

A:submatrix(p7,1,7);

rem : tridiagonal matrix can be solved  TDMA algorithm.

rem : submatrix(rows,matrix,columns)
columns = 1,7

2e)  How to find uu : column vector
-----------------------------------------------------
solution : a) first ,last element  remove from 'p6b'
                   b) transpose   a)

p10:matrix(p6b);

a) first ,last element  remove from 'p6b'
================================

p10a:submatrix(p10,1,n);

b) transpose a)
=============

uu:transpose(p10a);

3e)  How to find L : column vector  (lhs,final system)
-------------------------------------------------------------------------------
solution :  use  'p9' and 'p7a1' and 'p7b1'

p11a:p9;
p11b:p7a1;
p11c:p7b1;

a)  combination of 'p11a' and 'p11b' and 'p11c'
=======================================

1a) put values empty vector 'length m' from 'p11b,p11c'
=========================================

p12:[0,0,0,0,0];

p12[1] : p11b;
p12[m]:p11c;

p13:p12;

2a) lhs vector : 'vector : p13' -'vector: p11a'
====================================

p14:p13-p11a;

p14a:matrix(p14);

3a) transpose 'p14a'
=================

L:transpose(p14a);

rem : all values  L must be given from : data

How to solve 'final' system.
=======================
u[1,j] for  j = 2,3,4,5,6

p15:'A.uu=L;

p15a:A.uu=L;

How to find  : solution 'uu'
======================
1e) find  inverse matrix A = A^(-1)
2e)  uu  = A^(-1) .L   , ' matrix multplication'

1e) find  inverse matrix A = A^(-1)
--------------------------------------------------

p16: invert(A);

2e)  uu  = A^(-1) .L   , ' matrix multplication'
---------------------------------------------------------------

p16a:uu=p16.L;

numerical values    (no moving plates)
===============
Γ  = 0.02 kg/m*s
'diff(p,x,1)=-500 Pa/m
Δy  = 1*10^(-3) m
u[i,1] = 0   m/s              
u[i,7] = 0   m/s
because : u[i,1] = u[i,7] =wall , ussed   channel flow.

p17:['diff(p,x,1) =-500,Δy =10^(-3),Γ = 0.02,u[i,1]=0,u[i,n]=0];

p17a:subst([p17],p16a);

flow velocity : maximum at center of the channel.
                           minimum at walls of the channel.
(result : symmetric parabolic profile.)
---------------------------------------------------------------------------------
newton viscous fluid trough two dimensional channel
----------------------------------------------------------------------------------
=>  channel formed between two infinite plates.

numerical values    ( moving plates)
===============
Γ  = 0.02 kg/m*s
'diff(p,x,1)=-500 Pa/m
Δy  = 1*10^(-3) m
u[i,1] = 0   m/s                  (steady plate  )
u[i,7] = 10* 10^(-3)  m/s  (moving plate)
because : u[i,1] = u[i,7] =wall , ussed   channel flow.

p18:['diff(p,x,1) =-500,Δy =10^(-3),Γ = 0.02,u[i,1]=0,u[i,n]=10*10^(-3)];

p18a:subst([p18],p16a);

(result: non symmetric parabolic profile.)
rem : one of plates is moving and other is not moving.

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