page1
NUMERICAL SOLUTION OF COUTTEE FLOW
USING FINITE DIFFERENCE METHOD
THIS PROJECT THESIS IS SUBMITTED IN PARTIAL FULLFILMENT OF
THE REQUIREMENT FOR THE DEGREE OF BACHELOR OF
TECHNOLOGY IN MECHANICAL ENGINEERING
Project submitted by
UMAKANTA BEHERA
BIPIN KUMAR
NIKESH HEMBRAM
ARIJIT AS
ANIMESH SARDAR
Under the guidance of
SUBHAS CHANDRA RANA
Assistant Professor
Department of Mechanical Engineering
NATIONAL INSTITUTE OF TECHNOLOGY DURGAPUR
DEPARTMENT OF MECHANICAL ENGINEERING
NATIONAL INSTITUTE OF TECHNOLOGY DURGAPUR - 713209
page2
CERTIFICATE
This is to certify that the thesis titled Numerical Solution of Couttee Flow
using finite difference method has been prepared by Umakanta Behera
( Roll No-09/ME/19), Nikesh Hembram (Roll No-09/ME/20), Animesh Sardar (Roll
No-07/ME/17), Bipin Kumar (Roll No-08/ME/18), Arijit As (Roll No-08/ME/10)
under my guidance in partial fulfillment of the requirement for the award of the
degree Bachelor Of Technology in Mechanical Engineering under the
Department Of Mechanical Engineering , National Institute Of Technology
Durgapur,India.
S.C. Rana
B. Halder
Assistant Professor
Head Of The department
Mechanical Engineering Department
Department
NIT Durgapur
Mechanical Engineering
NIT Durgapur
ACKNOWLEDGEMENT
page3
We take this opportunity to express our profound gratitude and deep regards to our guide.
Prof. Subhas Chandra Rana for his esteemed guidance, monitoring and constant
encouragement throughout the duration of this project. The blessings, help and guidance given
by him time to time shall carry us a long way in the journey of life on which we are about to
embark.
Again we would like to thank the Head of the Department for providing a good
environment and facilities to complete this project. Last but not the least we would like to
thank all those who were directly or indirectly associated with this project and helped us in
various ways.
UMAKANTA BEHERA
NIKESH HEMBRAM
ANIMESH SARDAR
BIPIN KUMAR
ARIJIT AS
CONTENTS
TOPICS
PAGE
Why is this project
5
Introduction to finite difference method
Introduction to the Computational Fluid Dynamics
page4
Derivation of basic conservation equations
7
Momentum equation
9
Derivation of Couttee flow equation from basic
13
Numerical methods and solution using
33
o Crank-Nicolson implicit scheme
Thomas algorithm for TDMA
36
Back substitution method
38
Code in matlab
39
Couette Flow Graph
41
Results and discussions
43
Why is this project
page5
Works demonstrate the ability of student to
apply engineering technique in specified area or in
any functional area in specific manner.
It provides the student to develop the idea of
designing and complete the work in stipulated
period.
It enables the student to know themselves
getting knowledge of doing in a proper way which is
easier way and economic.
It also develops the spirit of student for team
work.
It develops the ability to solve the problem while
doing the job with patience .
It develops the skills and enhances the
knowledge in a specified area.
Introduction to Computational fluid dynamics and Finite
Difference Method
Fluid mechanics is the branch of physics that studies fluid and forces
on them.
page6
Fluid mechanics is further subdivide into fluid statics (fluid at rest),
fluid dynamics (fluid at motion) and fluid dynamics (the study of
effect of forces on fluid motion).
Especially fluid dynamics is now the active field of research with
many unsolved problems and some solved problems. Fluid dynamics
can be very complex mathematically. Sometimes it can be best
solved by Numerical Methods with the aid of computers. A modern
approach called Computational Fluid Dynamics (CFD) is devoted to
solve the fluid mechanical problems.
In mathematics, Finite Difference Methods are numerical methods for
approximation the solutions to differential equations using finite
difference equations to approximate derivatives. Explicit method,
Implicit methods and Crank-Nicolson method are the common
methods in it.
Derivation of basic conservation equations
Derivation of Continuity Equation
Reynolds Transport Equation
ur uur
dN
d V .dA
dt t cv
cs
page7
Here N = extensive property and,
N
m
i.e. extensive property(mass based) per unit mass.
= Intensive property (independent of mass)
For continuity equation, N = m then
= 1.
So continuity equation is
dm
0
dt
(at control mass)
Or,
Or,
Or,
uruur
d V dA
cs
t cv
uruur
d
d V dA
cv dt
cs
ur
d
d Vd
cv dt
cv
ur
dP
V d
cv
dt
Or,
As the control volume is arbitrary control volume, then the
integrant must be zero.
(There is no guarantee that the sum of the integration of one
part of the positive side and the other part on the negative
side of the control volume is zero)
That is why
P ur ur
.V 0
t
page8
ur
i j k
x
y
z
And
ur
ui vj wk
And
ur ur
.V 0
t
Or
.(a)
ur ur
ur ur
.V V .
t
Or,
Or,
ur ur
ur ur
V . .V
ur ur
D
.V
Dt
, where
D
Dt
is material derivative.
When the density is constant, ( =constant),then it is said to
be compressible fluid. And , when 2% change in the density
or mach number 0.3 then the fluid is said to be
incompressible fluid.
urur
0 0 V
uu
r uu
r
or,
.V 0
j k . ui vj wk 0
y
z
x
i
Or,
Or,
u v w 0
x y z
From equation (a)
page9
ur
ur
V
F m
t
or,
u v w
0
x
y
z
This is the differential form of an continuity equation of
incompressible fluid.
Momentum Equation
From Newtons 2nd law of motion,
uuu
r
uur
F ma
Or,
Or,
So, N=m
ur
V
ur
ur
V
F m
t
ur d
ur
F
mV
dt
uu
r
and
Now from Reynolds Transport equation,
uuur
d mV
r
uuur
uuur uuuruuuuu
V d V V .dA
cs
t cv
dt
Along x-direction,
r
uuruuuu
d mu
ud uV dA
cs
t cv
dt
uruur
d mu
u )d uV dA
cv
cv t
dt
Or,
page10
d mu
(u ) ur ur
cv
uV d
dt
t
Or,
ur ur ur
d mu
u
u.V V . u d
cv
dt
t
Or,
uruu
r ur ur
d mu
u
u u .V V . u d
cv
dt
t
t
Or,
uruu
r
ur ur
ur ur
d mu
u
u
.V u V . V . u d
cv
dt
t
t
Or,
uruu
r
ur ur
D
u
ur ur d mu
D
u
.V V . u d
.V
cv
dt
t
Dt
Dt
Or,
uruu
r
D
.V 0
Dt
And ,
d mu
cv
dt
u ur ur
V . u d
t
Or,
ur ur
d mu
V . u d
cv
dt
Or,
page11
d mu
Du
d
cv
dt
Dt
Or,
direction
.. in x-
Similarly in y- direction:d mv
Dv
d
cv
dt
Dt
Then in z-direction:d mw
Dw
d
cv
dt
Dt
ur d
ur
F
mV
dt
And,
uuuu
r uuuu
r
uuuu
r
F F surface Fbody
uuuu
r uuuu
r
uuuu
r
uuuu
r
F F viscous F pressure Fbody
uur
uur
uur
uur
F .d A p.d A gd
cs
cs
cv
uur
uur
.d A
Here
uur
cs
p.d A
=gauss divergence rule and
ur
ur
uur
F .d p.d gd
cv
uur
cv
cv
ur
ur
ur
F . p g d
cv
cs
=gradient rule
page12
xx xy xz
yx yy yz
zx zy zz
For x-direction only
cv
i j k xxi xyj xxk P g x d Du
Dt
x y
z
x
cv
P gx xx xy xx
Du
Dt
x
x y z
P gx u
Du
Dt
x
. for incompressible fluid
P gy v
Dv
Dt
y
in x- direction
P gz w
Dw
Dt
z
in y- direction
in z- direction
u v w 0
x y z
..................continuity equation
page13
COUETTE FLOW DERIVATION
Consider the viscous flow between two parallel plates separated by
the vertical distance D .The upper plate is moving at a velocity u e ,
and the lower plate is stationary; i.e., its velocity is u=0. The flow in
the xy plane is sketched below. The Flow Field between the two plates
is driven exclusively by the shear stress exerted on the fluid by
moving the upper plate, resulting in a velocity profile across the flow,
u=u(y).
The governing equation for
this flow is the x-momentum equation, given by equation (1),
p xx yx zx
Du
f x
Dt
x x y z
.(1)
page14
When applied to the coquette flow, this equation is greatly simplified,
as follows. Examining the figure, we note that the model for coquette
flow stretches to plus and minus infinity in the x- direction. Since
there is no beginning and end of this flow, the flow field variables
( any quantity)
0
x
must be independent of x; that is,
for all quantities.
Moreover from the continuity equation, written for the steady flow,
u v
0
x
y
Since
u x 0
(a)
for couette flow,equation (a) becomes,
v
v
v 0
y
y
y
..(b)
Evaluated at the lower wall, where v=0 at y=0, equation (b) yields
v
0
y y 0
v
0
y 0
Or,
.(c)
If we extend v in a taylor series about the point y=0, we have
page15
2v
v
y 2 / 2 .....
y
y y 0
y 0
v y v 0
.(d)
Evaluated at the upper wall,equation(d) becomes
2v
v
D
D 2 / 2 .....
y 2
y 0
y 0
v D v 0
..(e)
Since both v(D)=0 and v(o)=0, as well as
v
y
0
y 0
from equation (c),
the only result that makes sense in equation (e) Is that
n, and hence
v=0
n v
y n
0
y 0
for all
..(f)
everywhere .This is a physical quantity of couette flow , namely, that
there is no vertical component of velocity anywhere. This states that
the streamlines for Couette flow are straight, parallel streamlines a
result which is almost intuitively obvious simply by inspecting the
figure. Finally the y- momentum equation,
Dv
P
xy yy zy fy
Dt
y x
y
z
We have, for couette flow with no body forces,
0
P yy
y y
(g)
In aerodynamics problems, fluids can be assumed to be newtonian,
for such fluids,
page16
yy
Hence, with
yy 0
u v
v
2 0
x y
y
..(h)
,equation(g) yields,
p
0
y
.(i)
Conclusion: For coquette flow, there are no pressure gradients in
either the x or y direction. With all the above information, we return
to the x-momentum equation. From this equation, for unsteady, twodimensional flow with no body forces, we have
u u u v v P xx yx
t
x
y
x x
y
..(j)
As from the equation of newtonian fluids,
xx u v 2 u 0
x
x y
yx v u u
x y
y
(k)
(l)
Substituting (k) and (l) into (j),we have, for couette flow,
u u
t y y
.(m)
At this stage, we will now assume an incompressible, steady,
constant- tempreture flow for which =constant. With this, (m)
becomes
page17
2u 0
y 2
(n)
The above equation is the governing equation for incompressible,
constant-tempreture, steady couette flow.
The analytical solution,
u c1 y c2
u 0, y 0
c2 0
u ue, y D
ue
c1
D
Therefore, velocity profile,
u
y
ue D
Numerical solution:
Let us assume that velocity profile is not linear in unsteady flow and it
may be defined as
u 0 for 0 y D at t 0
ue for y D
Therefore it is evident that as time proceeds velocity profile turned to
be linear and flow becomes steady.
At
t t3
,velocity profile linear & flow is steady.
Now from x-momentum equation (m) without omitting the unsteady
term we get,
page18
u u
t
y 2
, which is parabolic and therefore it must be time
marching problem with given initial conditions & boundary conditions.
Let us first non-dimensionalise the qualitiesu ' u , y' y and t ' t
ue
D
D / ue
Therefore,
2
u u
t
y 2
(u / ue ) ue2
2 (u / ue ) ue
2 D 2
t / ( D / ue ) D
(
y
/
D
)
u * ue D
2u *
t *
y *2
u
1 2u
t Re y 2
D
For simplicity dropping primesu
1 2u
.................... A
t Re y 2
D
Nature of PDE and its stability criteria
if the partial differential equation be as given below :
page19
Axx Bxy C yy f (x,y,x , y )
Then For B2 4ac Hyperbolic PDE
For B2 4ac Elliptic PDE
For B2 4ac Parabolic PDE
For this case, C= 1 ,A=0,B=0, so B2 4ac
ReD
so,
u
1 2u
t
1
is parabolic and its stability criteria is
2
2
t R eD y
ReD y 2
Finite Difference Method(General Approach):
Finite difference method is useful to discretise the partial differential equations
of all the three kinds to get an approximate algebraic equations.The details of
the method is discussed below:
(A)Forward difference formula :The Taylor series expansion ,
u ( x) x 2 2u ( x) x 3 u ( x)
u ( x x) u ( x) x
.
.......
x
2! x 2
3! x3
May be written in terms of displacement operator, E and derivative operator, D
x 2 2 x 3 3
Eu ( x ) [1 xD
.D
.D ...]u ( x )
2!
3!
(1)
Where
x-
;
X+
E=displacement operator i,e it will displace the function by next nodal value i.e. at node x the
function u has value of u(x) and Eu(x) has value of u(x+
From equation(1)
) which is the value of u at the node (x+
page20
Eu ( x) exD u ( x)
therefore E e xD
D
1
ln E
x
With these definition,the derivative of u(x) at i may be expressed as under
ui u xi
ui 1 u xi 1 u ( xi x )
Eui u ( xi 1 ) u ( xi x )
therefore Eui ui 1
Eui 1 ui 2 etc from the defintion of E
therefore E 2ui ui 2
E nui ui n
ui
n=no. of times by which the function
xi
has been displaced from the node
ui ui 1 ui Eui ui ( E 1)ui
E 1
E=1+
therefore D
1
ln E
x
1
ln E (ui )
x
1
u
or ,
ln(1 )(ui )
x i x
or , D(u i )=
1 2 3 4
u
... ui
.... A
2
3
4
x i x
E 1 E 1 E 1 ... u
1
u
( B)
2
3
4
x i x
Or,
page21
So, by definition
Eui ui 1 , Eui 1 E 2ui ui 2 ;.......E nui ui n
For 1st order accuracy we will consider only one term on RHS
1ST order accuracy :- Only 1
st
term of equation (A) or equation (B) is considered.
Eu u u u
1
u
[ E 1]ui i i i 1 i
x
x
x i x
2nd order accuracy:- upto 2
nd
term of equation (A) or equation (B) is considered.
E 1 u
1
u
2
x i x
2
1
1 2
1 2
u
i
Eui ui E 2 E 1 ui ui 1 ui E ui 2 Eui u
2
2
x i x
1
1
1
= ui 1 ui ui 2 ui 1 ui
x
2
2
3
1
3
2ui 1 ui ui 2 3u 4u u
E 1
i
i 1
i2
2
2
=
ui
x
2x
3x
xd
3ui 4ui 1 ui 2 e 1
2x
3x
ui
Putting E=e xD
xD
1 xD
2!
3ui 4ui 1 ui 2
2x
3ui 4ui 1 ui 2 x 3 3
.D (ui ).
2x
3x
xD
3!
3x
....... 1
(ui )
3ui 4ui 1 ui 2 x 2 3u
u
2x
3 x 3
x i
Similarly 3rd order,4th order .nth order accuracy is possible.
(b)
Backward Difference Formula :
page22
u ( x x) u ( x) x
u ( x) x 2 u ( x ) x3 u ( x)
...
x
2! x 2
3! x 2
E 1 and D
The expression may be written in terms of
E u ( x)
x.D
1 x.D
2!
xD
1 x.D
2!
xD
3!
xD
3!
xD
4!
xD
4!
.... u ( x)
....
=e xD
-xD ln E 1
1
D=- ln( E 1 )
x
... C
Again
E 1ui ui 1 ; E 1ui 1 E 2ui ui 2 ...........E nui ui n
ui ui ui 1 ui E 1ui 1 E 1 ui
1 E 1 E 1 1
Now from equation (C) we have,
1
ln(1 )
x
1
Dui ln(1 )(ui )
x
D
page23
1 ( ) 2 ( )3 ( ) 4
u
ui
x
2
3
4
x i
1
u
1 E 1
x i x
1 E 1
2
1 E 1
3
1ST order accuracy :- Only 1
st
(C
)
u
i
( D)
term of equation (A) or equation (B) is considered.
ui E 1ui ui ui 1
1
u
1
[1 E ]ui
x
x
x i x
nd
order accuracy:- upto 2
nd
term of equation (A) or equation (B) is considered.
2
2
1 E 1
1
u
1
1 E 2 ui
x i x
u
E 2 u i
ui E 1ui i E 1ui
u
2
2
x
i
3ui E 2ui
2 E 1ui
u
2
2
x i
x
3ui 4ui 1 ui 2 x 2 3u
u
x i
2x
3 x3
Central difference Formula
1
E 2 ui
ui u
1
2
1
E2
1
E 2 ui
1 u
2
1
E 2
1
2
1
E2
1
2
1
E ui
2
page24
Using the definition of mean,
1
E2
1
E 2
from forward difference we get, E ex.D
1
E 2
From backward difference ,
1
E 2
1
E2
ex.D /2
e x.D
ex.D/2
E1/2 E 1/2 exD/2 exD/2 2sinh xD / 2 using sinh x
xD / 2 sinh 1
D
2sinh
x
1
2
/2
1
3 3 5 5 7
u
........ ui
24 640 7168
x i x
2 2
1 .3 .
12.32.52
3
2
1 / 2
2
u
2
2
/ 2
3!
5!
7!
x i x
1/ 2
1/2
E
E
1
u
u
1/2
1/2
E E
i
24
x i x
Thus,from equation(E) or (F) we have:-
......
( E)
e x e x
2
page25
2nd order accuracy with 1st term only:
3
1 1/2
1 1
u
E ui E 1/2u i
. E1/2 E 1/2 ui
x 24
x i x
e xD /2 exD /2
u
u
i 1/2 i 1/2
ui
x
24.x
u
u
x 2 D 3ui
i 1/2 i 1/2
x
24 24
u
u
x 2 3u
i 1/2 i 1/2
x
24 x3
4th order accuracy with two terms:
1 1/2
E 3/2 3E.E1/2 3E1/2 .E 1 E 3/2
3 5ui
u
1/2
E
u
E
u
i
i
24
640x
x i x
24ui 1/2 24 Ei 1/2 ui 3/2 3ui 1/2 3ui 1/2 ui 3/2 3. x.d
24x
640.x
5
u
27ui 1/2 27ui 1/2 ui 3/2
3
u
i 3/2
x 4 5
24x
640
x
ui
To avoid the half integer mesh point, the definition of mean may be employed.where =
1
2
E1/2 E 1/2
E1/2 E 1/2
2
1
1/2 E 1/2
E
4
2
1 2 4
2
1 2 /4
Now from equation F we have
3
5
7
u
3 5 ......
x
24 640 7168
i x 1 2 /4
2 2
3
5
7
1 3 5 ......
x
4
24 640 7168
2 3 4 15 6
1
...
x 8 816 4864
3 3 5 5 7
ui
......
24 640 7168
15. 7
3 3. 5 5. 7
3 5
3. 7
3. 5
3. 7
...
...
...
... u
24 640 7168
8 824 8640 816 81624 4864 i
page26
1
3
5
3
3
15
1 1
3
7
5
ui
x
24 8
640 8 24 8 16
7168 8 640 8 16 24 48 64
3 12 22 5 50 42 70 350 7
ui
x
3!
5!
71680
3 12 22 5 512 9 7
..... ui
x
3!
5!
71680 9
3 12 22 5
1 9
7 ui
x
3!
5!
2 7 10 9
3 12 22 5
12.22.32
7 ui
x
3!
5!
2 4 7 9 10
3 12 22 5 12.22.32 7
u
.... ui ......(G )
3!
5!
7!
x i x
From equation (G) 2nd order accuracy may be obtained as :
2nd order accuracy(with one term from R.H.S.):u
ui . ui
x i x
x 3!
1
1
1
12
2
2
2
E
E
E
1
12
ui
E E ui
x
2 3! x
1
1
D
D
x
x
12
12
x D2
x D2
2
2
2
2
E
E
E
E
e
e
e
ui
2x
2 3!x
1 x.D ... 1 x.D ... x.D
E E 1
ui
2x
2 3 2 x
3
3
u u
2x .D .ui
i 1 i 1
2x
2 3 2 x
u u
x 2 3u
i 1 i 1
2x
6 x 3
page27
Fourth order accuracy (with two terms on R.H.S )
u
3
.12.22. 5
ui
x i x
3!
x.5!
2 x
5
E E
1
1
1
1
12.22.2 x.D
2
2
2
2
E E E E
ui 2 x 5!
3!
1
2
1
2
1
1
3
1
1
23
1
2
2
2
2
2
E
E
E
3.
E
.
E
3.
E
.
E
1
1
u 4 x
EE
i 120
2 x
3!
1 6 E 6 E 1 E 2 3E 3 E 1 E 3 3.E 1 E 2
x 4 5ui
ui 30 x 5
2x
6
8E 8E 1 E 2 E 2
x 4 5ui
.ui
12x
30 x 5
ui 2 8ui 1 8ui 1 ui 2 x 4 5ui
12x
30 x 5
Higher order derivatives : Forward difference
page28
u
1
ln 1 ui ...... First order forward difference
x i x
nu
1
n ln 1
n
x i x
ui
taking nth power on both the sides
n
2 3 4
.... ui
2
3
4
n n 1 n 2
nu
1 n
2 3 4
2
n 1
n
.
...
....
x n i x n
2
3
4
2!
2
n n 1 n 2 4 6 8
2
1
3
4
n n n. n 1 n. n 1 .
n. n 1 .
.....
...
... ui
x
2
3
4
2!
4
9
16
1
n n 1 n 3n 5 n 2 n n 2 n 3 n 3
n n
.... ui 1
x
2
24
48
nu
1
n
n
x
x
Backward difference:Similarly u
1
ln 1 ui
x i
x
n
nu
1
n ln 1 ui
n
x i
x
2 3 4
.... ui
2
3
4
n
n 3n 5 n 2 n n 2
u
1
n
n n n n 1
......(2)
x
x
2
24
48
1
n
x
Central difference:u 1
1
2sin h ui
x x
2
3 3 5 5 7
...... ui
24 640 7168
n
n
n 22 5n
n 5 n 1
4
n 1 2
x
24
64
90
4 7
5
nu
1
n n
x
x
.....(3)
If n is even , the difference formulae are obtained at the integer mesh points.
No need for any manipulation.
page29
But if n is odd, then for getting integer mesh points following manipulation is
required.
2 sin 1 u
i
2 x
2
1
4
n
n 3 2 5n 2 52n 135 4
n 1
ui
x
24
5760
nu
x n
Second order Derivative (n=2)
From equation (1): Forward:
2u
1 2
11 4 5 5
ui
2
2
x
x
12
6
(5)
2u
1 2
11 4 5 5
ui
2
2
x
x
12
6
(6)
2u
1 2 4 6 8
2
ui
2
x
x
12 90 560
(7)
2u
2 5 4 259 6
8
2
0 x
ui
2
x
x
24
90
Forward difference. from equation 5
1.1st order accuracy,
(8)
page30
2u 1 2 3
x 2
x 2
x 2
ex D 1
2
1
ui
E 1 u1
x 2
x 2
2
E 1 u x 3 D 3 u
E 22
i
i
x
x 2
u
2ui 1 ui
3u
i 2
x
x 2
x3
2.2nd order accuracy:-
2u 1 E 1 2 E 1 3 11 E 1 4
12.x 2
x2 x2
2 2 E 1 E 3 3 E 2 3E 1
E
ui 11 2 x 4 .D 4ui
2
x
12.x
4E 2 5E 2 u 11 x2 4u
i
12
x 2
x4
u 4ui 2 5ui 1 2ui 11 2 4u
i 3
x
12
x 2
x4
Backward Difference method from equation(6)
(1).1st order accuracy:-
2u 1 2 1 3 1 E u x3 D3 u
x2 x2 i x2 i
x2 x2
3u
ui 2ui1 ui2
x.
x2
x3
Similarly, 2nd order accuracy,
2u 2ui 5ui 1 4ui2 ui 3 11 x2 4u
12
x 2
x 2
x4
Central Difference (integer mesh points) from equation
(7):
2nd order accuracy:-
page31
2u 1 2 4
12.x2
x2 x2
2
1
x 2.D 2
1/2
1/2
E
E
u
u
i
x 2
12.x 2 i
2 4
1 2 E E 1 2 ui x 4u
x
12.x
2
u 2ui ui1 x 4u
i1
x 2
12.x 4
(2).Fourth order accuracy:Similarly,
2u ui2 16ui1 30ui ui2 x4 6u
90 x6
12x2
x2
Central difference (half integer points) from equation (8)
Second order accuracy:
2u 2 5 4
2 ui
x 2 i x 2
x 24
2
4
E1/2 E 1/2 5 1/2
1/2
1/2
1/2
2 E E
ui
.
E
E
ui
24
x
2x 2
E1/2 E 1/2
5 x 4 4
1
2
u
i 24 x2 D ui
2x 2
4
E 3/2 E 3/2 E1/2 2 E1/2 2E 1/2 E 1/2 5
2 u
x
24
2x 2
x 4
u
u
u
u
5
4u
i 3/2 i 1/2 2i 1/2 i 3/2 x 2 4
24
2x
x
4th order accuracy: 2u
1
259 4 6u
x
5ui 5/2 39ui 3/2 34ui 1/2 34ui 1/2 39ui 3/2 5u
5760
x 2 48x 2
x6
Third order derivative:1)forward difference from equation (1)
1st order accuracy:-
page32
3u 1 3 u 3 3u
i
x3 i x2 i x2
4
3
1
ui 3 .x 4 . u
E
x3
2x3
x4
1
3x 4u
3
2
3 E 3E 3E 1 ui
2 x4
x
u 3ui 2 3ui1 ui 3x 4u
i 3
2 x 4
x3
2) backward difference from equation (2) -same as above
3) Central Difference (integer mesh points) from equation(4)
nth order accuracy
1
2
u E E
x3
2
3
1
2
1
3 3 5
3
ui
3
x
24
E
5
6
1
1
12
E E 12
2
2
. E E ui
E E ui
2x 3
2 24 x 3
1
3
1
1
12
32
2
2
2
2 1
E E E 3EE 3E E E ui
1
.x 5 .D 5 ui
3
3
2 x
4 x
E 2 3E 3 E 1 E1 3 3E 1 E 2 ui 1 2 5u
x
2x 3
4
x 5
1
1 2 5u
ui 2 2ui 1 2ui 1 ui 2 x 5
2x 3
4
x
1
2
1
2
1
2
1
2
Multidimensional Finite Difference
page33
xi x0 i.x;
yi yo j.y
There forward operators in x-&y- directions are
central difference = x & y
therefore,
Similarly backward,
& y
&
u
u
u
1
= xuij 0 x i1, j i, j 0 x
x ij x
x
u u
u
1
= xuij 0 x i, j i1, j 0 x
x ij x
x
Now using the crank Nicolson implicit method from discretisation
which is based on the finite difference method which we have
dissussed earlier, we get
u 1 n1 n
u
ui
t Vt i
2u
1 1 n1 n
. ui1 ui1 2 uin1 u in uin11 u in1
y 2 Vy 2 2
From equation(a) by discretisation with Nicolson-crank method
uin1 uiu
1 . 1 . 1 uin11 uin1 2 uin1 uin 1 uin11 u in1
Vt
Re 2 Vy 2
D
Vt
n1
uin1 uin
ui1 uin1 2uin1 2uin uin11 uin1
2 Vy Re
D
Grouping all the terms at time level n+1 on L.H.S the following
discretised equation is obtained:
-
Vt
2
2 Vy Re
.uin11 1
Vt
u n1
u n1
i
2
2
2 Vy R i 1
Vy Re D
eD
Vt
page34
Vt n
Vt
n 1
u
uin1
ui
i 1
2
2
Vy ReD
2 Vy ReD
2 Vy ReD
Vt
or, Auin11 Buin1 Auin11 ki
...............( B)
Where, A= -E/2
B 1 E and E=
Vt
2 Vy ReD
where, Ki uin 1 E E uin1 uin1
2
From equation (B), we get for i=2 to N
For i=2, Au1n1 Bu2n1 Au3n1 K 2
Au2n1 Bu3n1 Au4n1 K3
Au 3n1 Bu4n1 Au5n1 K 4
For i 3,
For i= 4,
..........(c)
n 1
n1
n 1
For i N, AuN 1 Bu N AuN 1 KN
So the equations in (c) can be written in the matrix form as follows,
1
0
0
0
0
.
.
0
0
0
0
B
A
0
0
.
.
.
.
0
0
A
B
A
0
0
0
0
0
0
0
0
A
B
A
0
0
0
0
0
.
.
0
A
B
A
0
0
0
0
.
.
.
0
A
.
.
0
.
.
.
.
.
.
0
A
.
.
0
.
.
.
.
.
.
0
A
.
A
0
0
0
0
.
.
.
0
A
B
0
0
0
.
.
.
0
0
A
1
( N 1) X ( N 1)
u1
u2
u3
u4
.
.
.
uN 1
uN
uN 1
K2
K
3
K4
.
.
K N 1
KN
1
( N 1,1)
..........( D)
N 1,1
page35
Now applying TDMA, we can solve the above simultaneous linear
equations numerically:
For stability,
1 t 1
Re D y 2 2
therefore,
2
t 1 Re D y
2
lets us take ,t E Re D y where E varied from 1 to 4000.
2
Discretised equations (B)(whose matrix form is given in equation (C) for i= 2 to N)
is solved on a grid as shown in the y-direction.
The vertical distance across the duct is divided into N equal increments of length
y by distributing N+1 grid points over the height D.
D
Therefore, y=
N
From the boundary condition, u1 0 & u N 1 1.
Therefore,equation (B) represents the N-1 equations with N-1 unknowns name ly, u 2 , u 3.....u N .
Where 1st and (N+1)th values of the velocity u (i.e. u1 =0 and u N 1 =1)
are the boundary values as per the problem statement.
Now
applying TDMA, we can solve the above simultaneous linear
equations numerically:
Thomas Algorithm for tri-diagonal matrix(D)
Thomas Algorithm for TDMA(Tridiagonal matrix Algorithm)
page36
1
0
0
0
0
.
.
0
0
0
0
B
A
0
0
.
.
.
.
0
0
A
B
A
0
0
0
0
0
0
0
0
A
B
A
0
0
0
0
0
.
.
0
A
B
A
0
0
0
0
.
.
.
0
A
.
.
0
.
.
.
.
.
.
0
A
.
.
0
.
.
.
.
.
.
0
A
.
A
0
0
0
0
.
.
.
0
A
B
0
0
0
.
.
.
0
0
A
1
( N 1) X ( N 1)
u1
u2
u3
u4
.
.
.
uN 1
uN
uN 1
K2
K
3
K4
.
.
K N 1
KN
1
( N 1,1)
..........( D)
N 1,1
Equation (A.1,A.2,..A.N) represents a system for N linear simultaneous equations with N unknowns
x , x ...........x
1
d x a x
1
b x d x a x
bx d x a x
................. A.1
.............. A.2
................. A.3
b x d x a x
c
........................................................................
bN -1 x d x a x c
.............. A.4
4
N -2
N -1
b x d x
N
N 1
N -1
N -1
N -1
............... A.5
Standard method for solving a system of linear algebraic equation is
Gaussian elimination.
page37
The Thomas algorithm is essentially the result of applying Gaussian
A.2 d A.1 b d b x d d x
1
d a x c d
1
- b d x ab x cb
................................................................ ...
d d ab x d a x c d cb
2
b a / d 1 x a x c b c / d
2
........................ A.6
d ' x a x c'
d ' d -b a / d , c' c -b c / d
2
elimination as follows:
.............. A.6
d '2 x2 a2 x3 c '2
By
A.3 d ' - b A.6 we get
2
d '2 b3 x2 d '2 d3 x3 d '2 a3 x4
d '2 c3
b3c '
d '2 b3 x2 b3a2 x3
...............................................................................
d3 b3a2 / d '2 x3 a3 x4 c3 b3c '2 / d '2
.......................... A.7
d '3 x3 a3 x4 c '3
The pattern which is developing is-
d -b a / d d '
d2
(A.6) can be obtained from (A.2) by dropping first term, replacing
2
by
, keeping
by
d
Similarly,(A.7) is obtained from(A.3) by dropping first term, replacing
c by c - b c ' / d '
3
Keeping third term unchanged and by replacing
So if we generalize the pattern-
c -b c / d c' .
c
third term unchanged and replacing
d -b a / d d '
3
by
c'
3
page38
d 'i di - ba
i i-1 / d 'i-1 where i 2,3..............N
& c 'i ci - bc
i 'i-1 / d 'i-1 where i 2,3..............N
for i=3 to N
c(i 1,1)b(i, i 1)
d '(i 1, i 1)
a (i 1,1)b(i, i 1)
d '(i, i ) d (i, i )
d '(i 1, i 1)
c '(1,1) c(1,1)
0
0
0
0
.
.
0
0
0
0
0
. . . .
0
d 22 a23 0
. . . .
.
'
0 d33 a34 0 . . .
.
'
0
0 d 44 a45 0 . .
.
.
0
0
. . 0 .
.
.
.
0
0 . . 0
.
.
.
.
0 0 . .
0
.
.
.
. 0 0 . aN 1, N
'
.
.
.
. . 0 0 d NN
0
0
.
. . . .
0
0
0
0
.
.
.
.
.
.
.
aN , N 1
1
n 1
x1
x2
x3
x4
x5
N 1 XN 1
xN 1
xN
xN 1
N 1 X 1
0
C21
C '31
C '41
.
.
.
.
CN ,1
1
Back substitution method: Back-Substitution
The process of solving a linear system of equation that has been transformed row
echelon form. The is solved equation first, then the next-to-last, etc.
Example:
Consider a system with the given row-echelon form for its
augmented matrix.
The equations for this system are
page39
The last equation says z = 2. Substitute this into the second
equation to get
Now substitute z = 2 and y = 13 into the first equation to
get
Thus the solution is x = 24, y = 13, and z = 2.
x(i,1)
c(i,1) a(i, i 1).x(i 1,1)
for i=N to 2.
d (i, i)
But in our case,
u (i,1)
c(i,1) a(i, i 1).u (i 1,1)
for i=N to 2.
d (i, i)
MATLAB PROGRAM FOR COUETTE FLOW
n=50;
dy=1/n;
%pre alocation of the matrices
AA=zeros(n+1);
u=zeros(n+1,1);
c=zeros(n+1,1);
y=zeros(n+1,1);
%boundary conditions
y(1)=0;y(n+1)=1;
u(1)=0;u(n+1)=1;
%grid
page40
for i=2:n
y(i)=(i-1)*dy;
end
dt=0.01;Re=25;
E=dt/(Re*(dy)^2);
A=-E/2;B=1+E;
it=0;
dum=1;
eps=1.e-10;
while(dum==1)
dum=0;it=it+1;
for i=2:n
uold=u(i,1);
end
for i=1:n+1
d(i,i)=B;%main diagonal
end
d(1,1)=1;d(n+1,n+1)=1;
for i=2:n
a(i,i+1)=A;%super diagonal
end
for i=3:n
b(i,i-1)=A;% sub diagonal
end
for i=2:n
c(i,1)=(1-E)*u(i,1)+(E/2)*(u(i+1,1)+u(i-1,1));
end
c(n+1,1)=1;
% back substitution
for j=3:n
d(j,j)=d(j,j)-a(j-1,j)*b(j,j-1)/d(j-1,j-1);
c(j,1)=c(j,1)-c(j-1,1)*b(j,j-1)/d(j-1,j-1);
end
for i=n:-1:2
u(i,1)=(c(i,1)-a(i,i+1)*u(i+1,1))/d(i,i);
end
if it==1;plot(u,y);hold on;end
if it==20;plot(u,y);hold on;end
if it==40;plot(u,y);hold on;end
if it==60;plot(u,y);hold on;end
if it==80;plot(u,y);hold on;end
if it==100;plot(u,y);hold on;end
if it==200;plot(u,y);hold on;end
if it==300;plot(u,y);hold on;end
if it==400;plot(u,y);hold on;end
if it==500;plot(u,y);hold on;end
if it==600;plot(u,y);hold on;end
if it==700;plot(u,y);hold on;end
if it==800;plot(u,y);hold on;end
if abs(u(j)-uold)>eps; dum=1; end
page41
end
plot(u,y,'ro-');
COUETTE FLOW
1
0.9
0.8
0.7
Distance,y
0.6
0.5
0.4
0.3
0.2
0.1
0
0.1
0.2
0.3
0.4
0.5
Velocity,u
0.6
0.7
0.8
Result and disscussion
0.9
page42
From the governing equation of Couette Flow it is found that the PDE
obtained is parabolic in nature which has been discussed in the
nature of PDE article.Therefore it is expected that the solution of the
equation of parabolic PDE should be time marched.From the graph it
is found that at different time different velocity profile was
obtained.The first velocity profile is at time,t =0.01 (non-dimensional)
and as time marched forward the velocity, within the two plates as
discussed in the Couette Flow definition,is gradually developing and
finally at sufficient time (t=36.17) the velocity profile becomes fully
developed and merged with the steady state velocity profile as it is
u
y
ue D
given by
or by putting ue 1 and D=1 the non-dimensional form of the solution
comes to be u=y which is a straight line with slope=1.
page43
