Introduction to Computational Fluid Dynamics

Fluid Mechanics and Heat Transfer. Basic equations.

Continuity (Mass Conservation)
  
divv   0  div  v   0
t
incompressible

Momentum Equation (Navier-Stokes equations)

Incompressible viscous fluid (   constant ,   constant )

v    
    v   v    F  grad p   v
 t 

Lecture 2 –Basic equations. Dimensionless parameters 1

 Introduction to Computational Fluid Dynamics

Energy equation
dE d 1
   p      divk grad T 
dt dt   
where E is the internal energy of a fluid particle (viscous compressible fluid)

- gas heating due to the compression

d 1 
 p     p div v
dt    continuity equation
- mechanical work of the viscous forces transformation in heat

Lecture 2 –Basic equations. Dimensionless parameters 2

 Introduction to Computational Fluid Dynamics
p
By taking into account that the enthalpy of a unit of mass is i  E  we

obtain
di dp
     divk grad T 
dt dt
But for a perfect gas we have
di  c p dT
and energy equation becomes


d
c pT   dp    divk grad T 
dt dt
where
df Df  f 
    v   f
dt Dt t
is the substantial derivative.

Lecture 2 –Basic equations. Dimensionless parameters 3

 Introduction to Computational Fluid Dynamics
Momentum Equation (for fluid saturated porous media)
By a porous medium we mean a material consisting of a solid matrix
with an interconnected void.

Porous media are present almost always in the surrounding medium, very
few materials excepting fluids being non-porous.

heat exchanger processor cooler porous rock

Lecture 2 –Basic equations. Dimensionless parameters 4

 Introduction to Computational Fluid Dynamics

Darcy’s (law) equation (momentum equation)

 K 
v   grad p   g

where
 v is the average velocity (filtration velocity, superficial velocity,
seepage velocity or Darcian velocity)
 K  k g is the permeability of the porous medium

Darcy’ law expresses a linear dependence between the pressure gradient

and the filtration velocity. It was reported that this linear law is not valid
for large values of the pressure gradient.

Lecture 2 –Basic equations. Dimensionless parameters 5

 Introduction to Computational Fluid Dynamics
Darcy’s law extensions
Forchheimer’s extension
cF  
 
p    v v v 
K K
where cF is a dimensionless parameter depending on the porous medium

Brinkman’s extension
  
p   v  ~  v
K
where ~ is an effective viscosity

Brinkman-Forchheimer’s extension
  cF   
p   v   v v  ~  v
K K

Lecture 2 –Basic equations. Dimensionless parameters 6

 Introduction to Computational Fluid Dynamics

Energy equation for porous media

T
 c p m   c p  fluid v  T  k e  T


t
where

 c p m    c p  fluid  1    c p solid

k e   k fluid  1   k solid

Lecture 2 –Basic equations. Dimensionless parameters 7

 Introduction to Computational Fluid Dynamics
Bouyancy driven flow1

 There exists a large class of fluid flows in which the motion is caused by
buoyancy in the fluid.
 Buoyancy is the force experienced in a fluid due to a variation of density in
the presence of a gravitational field. According to the definition of an
incompressible fluid, variations in the density normally mean that the fluid is
compressible, rather than incompressible.
 For many of the fluid flows of the type mentioned above, the density variation
is important only in the body-force term of the momentum equations. In all
other places in which the density appears in the governing equations, the
variation of density leads to an insignificant effect.
 Buoyancy results in a force acting on the fluid, and the fluid would accelerate
continuously if it were not for the existence of the viscous forces.
 The situation depicted above occurs in natural convection.

1
I.G. Currie, Fundamental Mechanics of Fluids, 3rd ed., Marcel Dekker, New York, 2003

Lecture 2 –Basic equations. Dimensionless parameters 8

 Introduction to Computational Fluid Dynamics

Generally, density is a function of temperature and concentration

   T , c 

Black Sea salinity, temperature

and density
(Ecological Modelling, 221,
2010, p. 2287-2301)

Lecture 2 –Basic equations. Dimensionless parameters 9

 Introduction to Computational Fluid Dynamics

Boussinesq’s approximation

,
Momentum

Lecture 2 –Basic equations. Dimensionless parameters 10

 Introduction to Computational Fluid Dynamics
Dimensionless equations
u v w
  0
x y z

 u u u u  p   2u  2u  2u 
   u  v  w   Fx     2  2  2 
 t x y z  x  x y z 

 v v v v  p   2v  2v  2v 
   u  v  w   Fy     2  2  2 
 t x y z  y  x y z 

 w w w w  p  2w 2w 2w 

   u v  w   Fz     2  2  2 
 t x y z  z  x y z 

 T T T T    2T  2T  2T 
c p   u v  w   k  2  2  2     q( 0)
 t x y z   x y z 

  u  2 
 v
2

 w
2
 v u 
2
 w v 
2
 u w  
2
   2   2   2             
  x   x   x   x y   y z   z x  

Lecture 2 –Basic equations. Dimensionless parameters 11

 Introduction to Computational Fluid Dynamics
We introduce further the dimensionless variables

t x y z u v w  F p T  T0
  , X  ,Y  , Z  ,U  ,V  ,W  , F ' , P  ,  
t0 L L L U0 U0 U0 g p0 T

into the governing equations

UU 0    U 0 U ,
u 
u 
  UU 0  X  U 0 U ,
t  t t0  x X x L X
p 
 P p0  X  p0 P ,
x X x L X

 2u   u    U 0 U  U 0  2U X U 0  2U
      2
x 2 x  x  x  L X  L X x L X 2
2

T


  T  T0    T  ,
T


  T  T0  X  T 
t  t t0  x X x L X

 2T   T    T   T  2 X T  2U
      2
x 2 x 
 
x x  L X  L X 2 x
L X 2

Lecture 2 –Basic equations. Dimensionless parameters 12

 Introduction to Computational Fluid Dynamics
and we obtain the dimensionless equations

U V W
  0
X Y Z

L U U V W gL p0 P    2U  2U  2U 
U V W  2 F ' x     2  2
t0U 0  X Y Z U0  U0 2 X U 0  X
L 2
Y Z 

 k 
  
L    
Q   c p    2  2  2  U 0
 
t0U 0 
U
X
V
Y
W
Z

U0L   X 2 Z 2 Z 2   c LT 
 
  p

Lecture 2 –Basic equations. Dimensionless parameters 13

 Introduction to Computational Fluid Dynamics
Dimensionless numbers
A. Strouhal number
The Strouhal number (St) is a dimensionless number describing oscillating flow
mechanisms:
L L oscillation
St   
t0U 0 U 0 average velocity
B. Froude number
The Froude number (Fr) is a dimensionless number defined as the ratio of the flow inertia
to the external field (the latter in many applications simply due to gravity).
In naval architecture the Froude number is a very significant figure used to determine the
resistance of a partially submerged object moving through water, and permits the
comparison of similar objects of different sizes, because the wave pattern generated is
similar at the same Froude number only.
U0
Fr 
gL

Lecture 2 –Basic equations. Dimensionless parameters 14

 Introduction to Computational Fluid Dynamics
C. Euler number
It expresses the relationship between a local pressure drop e.g. over a restriction and the
kinetic energy per volume, and is used to characterize losses in the flow, where a perfect
frictionless flow corresponds to an Euler number of 1.
p0
Eu 
U 0 2
Usually, we take p0  U 0 2 and we obtain Eu = 1.

D. Reynolds number
U 0 L
U0L
Re  
 
The Reynolds number is defined as the ratio of inertial forces to viscous forces and
consequently quantifies the relative importance of these two types of forces for given
flow conditions.
 laminar flow occurs at low Reynolds numbers, where viscous forces are dominant,
and is characterized by smooth, constant fluid motion;
 turbulent flow occurs at high Reynolds numbers and is dominated by inertial forces,
which tend to produce chaotic eddies, vortices and other flow instabilities.

Lecture 2 –Basic equations. Dimensionless parameters 15

 Introduction to Computational Fluid Dynamics
E. Prandtl number
cp / 
Pr   
k k / c p 

It is the ratio of momentum diffusivity to thermal diffusivity. The Prandtl number

contains no such length scale in its definition and is dependent only on the fluid and the
fluid state. As such, the Prandtl number is often found in property tables alongside other
properties such as viscosity and thermal conductivity.
•gases - Pr ranges 0.7 - 1.0 •water - Pr ranges 1 - 10
•liquid metals - Pr ranges 0.001 - 0.03 •oils - Pr ranges 50 - 2000

E. Eckert number
It expresses the relationship between a flow's kinetic energy and the boundary layer
enthalpy difference, and is used to characterize heat dissipation.

U 02
Εc 
c p T

Lecture 2 –Basic equations. Dimensionless parameters 16

 Introduction to Computational Fluid Dynamics
Using the above dimensionless numbers we obtain:
U V W
  0
X Y Z
U U V W 1 P 1   2U  2U  2U 
St U V W  F ' x  Eu    2  2 
 X Y Z Fr 
X Re  X 2
Y Z 
….
   Q 1   2  2  2  Ec
 
St

U
X
V
Y
W
Z
 
Re Pr  X 2

Z 2

Z 2  Re 


or using differential operators


   V 0 (23)
V   1  1 
St   V   V  F ' Eu P  V (24)
 Fr Re
 
  V    
1 Ec
St    (25)
 Re Pr Re

Lecture 2 –Basic equations. Dimensionless parameters 17

 Introduction to Computational Fluid Dynamics
Comments
 Dimensionless numbers allow for comparisons between very different systems and
tell you how the system will behave
 Many useful relationships exist between dimensionless numbers that tell you how
specific things influence the system
 When you need to solve a problem numerically, dimensionless groups help you to
scale your problem.
 Analytical studies can be performed for limiting values of the dimensionless
parameters

Express the governing equations in dimensionless form to:

(1) identify the governing parameters
(2) plan experiments
(3) guide in the presentation of experimental results and theoretical solutions

If the flow is not oscillatory, we take usually t0  L U 0 and thus St =1. If t0   we have
St = 0 and the flow is steady.

Lecture 2 –Basic equations. Dimensionless parameters 18

 Introduction to Computational Fluid Dynamics
Numerical methods for initial values problems (IVP or Cauchy problem)
Pendulum
Mathematical model:
d 2 g d
L  sin   0 ,  (t 0 )   0 , (t 0 )   ' 0
dt 2 L dt

Radioactive elements decay:

dm
  m , m(t0 )  m0
dt

where m is the mass of the radioactive element and  is the decay constant.

Lecture 2 –Basic equations. Dimensionless parameters 19

 Introduction to Computational Fluid Dynamics
Theoretical aspects

The general Cauchy problem:

 y ' (t )  f (t , y ), at b

 y (a)  y a
Using numerical methods we obtain a discrete approximation of y in
some points, called nodes, which form a grid (mesh).

A grid for the interval [a, b] having a constant step h is:

a, a  h, a  2h, ..., a  ih,..., a  ( N  1)h, b
or
t i  a  (i  1)h, i  1, ..., N
where N is the number of subintervals, and the step is constant:
h  (b  a) /( N  1) .

Lecture 2 –Basic equations. Dimensionless parameters 20

 Introduction to Computational Fluid Dynamics

By using a numerical method we obtain the approximated discrete values

not
y (t i )  yi , i  1, ..., N

y(t)

yi+2
yi+1
yi-1 yi
yi-2
t
i-2 i-1 i i+1 i+2

Lecture 2 –Basic equations. Dimensionless parameters 21

 Introduction to Computational Fluid Dynamics
One step methods
An example: The Euler method (Leonhard Euler-(1707 – 1783))
 y ' (t )  f (t , y ), at b

 y (a)  y a
By using a Taylor expansion (y should be of
class C 2 ) we have:
h2
y (ti 1 )  y (ti  h)  y (ti )  h y ' (ti )  y" (i ), i  (ti , ti 1 )
2
Thus, one get:
h2
y (ti 1 )  y (ti )  h f (ti , y (ti ))  y" ( i )
2
By dropping the last term the Euler formula is obtained:
y0  y a
yi 1  yi  h f (ti , yi ), i  0,1, ..., N  1
Remark: yi 1 depends on yi , t i and h .

Lecture 2 –Basic equations. Dimensionless parameters 22

 Introduction to Computational Fluid Dynamics
(Richard L. Burden and J. Douglas Faires , Numerical Analysis,Ninth Edition, 2011).

Lecture 2 –Basic equations. Dimensionless parameters 23

 Introduction to Computational Fluid Dynamics

One may notice that the

Euler method follow the
tangent for the current
node to approximate the
value in the next node.

Lecture 2 –Basic equations. Dimensionless parameters 24

 Introduction to Computational Fluid Dynamics
Same problem for h = 0.2 (11 nodes). Comparison with the

Lecture 2 –Basic equations. Dimensionless parameters 25

 Introduction to Computational Fluid Dynamics
An one step method for the Cauchy problem is given by
yi 1  yi  h  (ti , yi , h), [t , y ]  [a, b]  R n , h  0
 : [a, b]  R n 1  R n
Considering the exact solution (in the grid points) ~ y (t i ) we define the
local truncation error as follow:
yi 1  yi ~y (ti 1 )  ~
y (ti ) ~
y (ti 1 )  ~
y (ti )
T (ti , yi , h)     (ti , yi , h) 
h h h
(the difference between the approximation increment and the exact increment for one step)

A method  is consistent if T (t , y, h)  0 uniformly when h  0 for

(t , y )  [a, b]  R n .
(it is necessary that  (t , y ,0)  f (t , y ) )

Lecture 2 –Basic equations. Dimensionless parameters 26

 Introduction to Computational Fluid Dynamics
A method  is of order p if
T (t , y, h)  Kh p

uniformly on [a, b]  R n , where is a vector norm, and K is constant.

We may use the following notation:
T (t , y, h)  O(h p ), h  0
(for p  1 the method is consistent).
It can be shown that Euler method is a method of order one, O(h) .

Lecture 2 –Basic equations. Dimensionless parameters 27

 Introduction to Computational Fluid Dynamics
Improved Euler methods
evaluation of the tangent is made in an intermediary point of the interval ti , ti 1 

y(t)

slope = f(ti+1, yi+hf(ti, yi))

yexact
ymodified
yEuler
yi average slope
ti ti+1/2 ti+1 t slope = f(ti, yi)

Modified Euler

ti ti+1
Heun
Lecture 2 –Basic equations. Dimensionless parameters 28
 Introduction to Computational Fluid Dynamics
not
-evaluation in the middle of t i , t i 1 , ti 1 / 2  ti  h
1
2
 
 1 1 
y (t i 1 )  y (t i )  h f t i 1 / 2 , yi 1 / 2   y (t i )  h f  t i  h, yi  hf (t i , yi ) 
 2  2  

 Euler 
 
 1 1 
y (t i 1 )  y (t i )  h f  t i  h, yi  h f (t i , yi ) 
 2 2   

K1 (ti , yi ) 
  
K 2 (ti , yi ,h )

K1 (ti , yi )  f (ti , yi ) K1 (t i , yi )  f (t i , yi )
 1 1  K 2 (t i , yi , h)  f t i  h, yi  hK1 
K 2 (ti , yi , h)  f  ti  h, yi  hK1 
 2 2  1
yi 1  y i  h K 1  K 2 
yi 1  yi  hK 2 2
modified Euler method Heun’s method (trapezoidal rule)

Lecture 2 –Basic equations. Dimensionless parameters 29

 Introduction to Computational Fluid Dynamics
Runge-Kutta type methods
evaluation in intermediary points of t i , t i 1 

Standard Runge-Kutta method

h
y i 1  y i  K1  2 K 2  2 K 3  K 4 
6
K 1  f (t i , y i )
 1 1 
K 2  f  t i  h, y i  hK1
 2 2 
 1 1 
K 3  f  t i  h, y i  hK 2 
 2 2 
K 4  f t i  h, y i  hK 3 
Carl David Tolmé Runge Martin Wilhelm Kutta
(1856 – 1927) (1867 – 1944)

Lecture 2 –Basic equations. Dimensionless parameters 30

 Introduction to Computational Fluid Dynamics
Matlab solvers for Cauchy problems (ODE)
Syntax

[t,Y] = solver(odefun,tspan,y0,options, p1, p2, ...)

or
sol = solver(odefun,[t0 tf],y0...)
where solver can be
ode45, ode23, ode113, ode15s, ode23s, ode23t or ode23tb.

Input parameters (selection)

odefun – right hand member of the Cauchy problem
tspan –integration interval.
y0– initial value
options – solver options.

Lecture 2 –Basic equations. Dimensionless parameters 31

 Introduction to Computational Fluid Dynamics
Output parameters:
t – column vector of time points;
y – solution array
sol – solution structure

Lecture 2 –Basic equations. Dimensionless parameters 32

 Introduction to Computational Fluid Dynamics

Lecture 2 –Basic equations. Dimensionless parameters 33

 Introduction to Computational Fluid Dynamics
Example: Solve pendulum equation using ode45:

d 2 g d
2
 sin   0 ,  ( 0)  0, (0)  0.1
dt L dt
We note y1   and rewrite the system in the form
dy1
 y2 ;
dt
dy2 g
  sin y1
dt L
dy1
y1 0   0; 0  0.1
dt

function dy=fPendul(t,y,flag,g,L)
dy=[y(2);-g/L*sin(y(1))];

Lecture 2 –Basic equations. Dimensionless parameters 34

 Introduction to Computational Fluid Dynamics
%pendulum equation
a=0;b=pi/2;%integration interval
g=9.8;%accelaration due to the gravity
L=0.1;%length of the pendulum
y0=[0,0.1];%initial conditions
options=odeset('RelTol',1e-8);%modify the options
%use the solver
[t,y]=ode45('fPendul',[a,b],y0,options,g,L)
plot(t,y(:,1));

Lecture 2 –Basic equations. Dimensionless parameters 35