CFD for the design of energy systems

Lecture 1: Overview of fluid dynamic equations and their numerical resolution

Introduction

Fluid dynamic involves in various energy systems as:

• Combustion
• Wind turbine, Hydraulic turbine, hydrokinetic turbine
• Nuclear power plant
• …

Design of energy systems Design of energy systems Design of energy systems

• Better understanding of flow physical behavior is needed for design of • Better understanding of flow physical behavior is needed for design of • Better understanding of flow physical behavior is needed for design of
energy systems energy systems energy systems

Temperature flux between wall

Energy dissipation due to
and flow and temperature mixing:
generation of turbulent flows in
safety of nuclear power plan
Atomisation and mixing when hydrokinetic turbine
chemistry reactions are used

When numerical tools (CFD for Computational Fluid Dynamic) are used the limitation of this tool has to be
known.

Real flow behavior ßà Numerical model prediction ??

The limitation can be due to:
Governing
• the modeling step:
modeling equations the assumptions done to derive the set of equations of
the problem
Real flow
behavior Numerical
resolution • the numerical resolution step:
the numerical technics used to solve this set of equations
An engineer performing CFD studies has to know:
Simulation
- the impact of the numerical resolution

- the limit of the model used

 I. Overview of fluid dynamic equations and their numerical resolution

I. a. Governing equations

Reference (books):
Introduction to Fluid Mechanics, by Y. NAKAYAMA
Fundamentals of Fluid Mechanics, by B. MUNSON, D. YOUNG and T. OKIISHI
Advanced Fluid Mechanics, by W. GRAEBEL

The fluid dynamic governing equations are obtained from the 3 basic conservation principles:
- conservation of mass
nd
- conservation of momentum (2 Newton principle)
- conservation of energy

Applying these conservation principles with constant density assumption lead to the system for incompressible
flows (with index notation and Einstein’s convention, i.e. sum on index appearing two times in the same term):

𝜕𝑢#
=0 (mass conservation)
𝜕𝑥#
𝜕𝑢# 𝜕𝑢# 𝑢7 1 𝜕𝑝 𝜕 = 𝑢#
+ =− +𝜈 (momentum conservation)
𝜕𝑡 𝜕𝑥7 𝜌 𝜕𝑥# 𝜕𝑥7 𝜕𝑥7

with 𝑢# the velocity component in the direction 𝑥# , 𝑝 the pressure, 𝜌 the density and 𝜈 = 𝜇/𝜌 the kinematic
viscosity.

This system is closed with 4 unknown terms 𝑝, 𝑢 for 4 equations, in 3D. This means that the energy conservation
principles is not needed for incompressible flows.

For compressible flows, the system writes (for perfect gas assumption),
𝜕𝜌 𝜕𝜌𝑢#
+ =0 (mass conservation)
𝜕𝑡 𝜕𝑥#
𝜕𝜌𝑢# 𝜕𝜌𝑢# 𝑢7 𝜕 1 𝜕𝑢7 𝜕 𝜕𝑢#
+ =− 𝑝+ 𝜇 + 𝜇 (momentum conservation)
𝜕𝑡 𝜕𝑥7 𝜕𝑥# 3 𝜕𝑥7 𝜕𝑥7 𝜕𝑥7
𝜕𝜌𝑇 𝜕𝜌𝑢# 𝑇 𝜕 𝜕𝑇
+ = 𝛼 (energy conservation)
𝜕𝑡 𝜕𝑥# 𝜕𝑥7 𝜕𝑥7
𝑝 = 𝜌𝑟𝑇 (thermodynamic relation)
with 𝑇 the temperature, 𝛼 the thermal diffusivity and 𝑟 the constant for perfect gas.

This system is closed with 6 unknown terms 𝑝, 𝑢, 𝑇, 𝜌 for 6 equations, in 3D. Thermodynamic relation has been
used to close the system.

Final remarks about incompressible and compressible flows:

• Why free divergence velocity (𝜕𝑢# /𝜕𝑥# = 0) means “incompressible flows”

 It can be shown that the divergence of the
velocity field is linked with the variation of
the volume of fluid

In incompressible case, the free

divergence velocity means that there is no
variation of fluid volume. The flow can
then not be “compressed” …

• Difference on the pressure between “incompressible” and “compressible” flows

Difference between incompressible and compressible flows appear when the propagation of pressure wave are
considered.
In a 1D, inviscid case (without viscosity), without fluid motion (without fluid velocity), the propagation of a
pressure disturbance, 𝑝K , is solution of
𝜕 = 𝑝K =
𝜕 = 𝑝K
− 𝑎 =0
𝜕𝑡 = 𝜕𝑥 =
=
with 𝑎 the speed of sounds (equation obtained from a linearization procedure, using definition of speed of
sound).

This means that the pressure waves are propagated at the speed of sound. In incompressible case, 𝑎 → ∞ (since
𝑀𝑎𝑐ℎ 𝑛𝑢𝑚𝑏𝑒𝑟 → 0), then the pressure waves are propagated instantaneously.
 I. b. Basics on numerical resolution

Reference (books):
Computational methods for fluid dynamics, by J. FERZIGER and M. PERIC
Numerical heat transfers and fluid flow, by S. PATANKAR
Computational fluid dynamics for engineers, by T. CEBECI, J. SHAO, F. KAFYEKE and E. LAURENDEAU

The non-linear equations governing the fluid motion have no analytical solution. Numerical methods are then
used to approximate solutions. Various methods can be used: Finite Differences, Finite Volumes, Finite
Elements, …

Let consider a typical partial differential equation (PDE) to solve (in 1D):
𝜕𝜙 𝜕𝜙
+𝑐 =0
𝜕𝑡 𝜕𝑥
with 𝑐 a constant and 𝜙 the unknown term.

The numerical method first defined a mesh (a finite number of points in the continuous domain) where 𝜙 will
be evaluated.

Discretization step

The numerical resolution consists in an “approximation” of the derivative terms. This is the discretization.

For example, using finite differences method (based on Taylor series expansion), it can be written for
𝜙 𝑥# , 𝑡 W = 𝜙#W :
𝜕𝜙 W 𝜕𝜙 W 𝜙#WXK − 𝜙#W W
𝜙#XK W
− 𝜙#ZK
+𝑐 = +𝑐 + 𝒪(Δ𝑡, Δ𝑥 = ) = 0
𝜕𝑡 # 𝜕𝑥 # Δ𝑡 2Δ𝑥
with Δ𝑡 and Δ𝑥 the time step and the space step, respectively.

The PDE is the discretized using a first order scheme in time and a second order scheme in space.

Accuracy of scheme

The accuracy of the scheme is given by its

order. The “modified wave number analysis”
allows to better understand the scheme
capabilities.

Starting from the function 𝑓 𝑥 = cos 𝑘𝑥

(exact derivative: 𝑓 _ 𝑥 = −𝑘 sin 𝑘𝑥 ), the
idea is to compare the exact derivative with
the “approximate” derivative written as,
𝑓 _mod 𝑥 = −𝑘`aa sin 𝑘𝑥 ).

This shows that the schemes are accurate for

low wave number but not for high wave
number (stiff gradient)
Resolution of linear system

The discretized equation writes:

𝜙#WXK − 𝜙#W W
𝜙#XK W
− 𝜙#ZK
+𝑐 =0 (1)
Δ𝑡 2Δ𝑥
leading to the matrix form:
1 −𝛼
0
−𝛼 ⋱ ⋱
⋱ ⋱ ⋱
𝜙 WXK = 𝜙W + 𝑆 (system 1)
⋱ ⋱ ⋱
⋱ −𝛼 ⋱
0
−𝛼 1
de
with 𝜙 W the vector composed the value of the unknown quantity in all points of the mesh and 𝛼 = −𝑐 .
=df

Another scheme could be:

𝜙#WXK − 𝜙#W WXK
𝜙#XK WXK
− 𝜙#ZK
+𝑐 =0
Δ𝑡 2Δ𝑥
leading to the matrix form:
1 𝛼
0
𝛼 ⋱ ⋱
⋱ ⋱ ⋱
𝜙 WXK = 𝜙 W + 𝑆 (system 2)
⋱ ⋱ ⋱
⋱ ⋱ 𝛼
0
𝛼 1

• (system 1) is an explicit scheme: knowing 𝜙 W , 𝜙 WXK is directly available

• (system 2) is an implicit scheme: knowing 𝜙 W , a linear system has to be solved to find 𝜙 WXK .

Remark: the algorithms to solve linear system are often iterative algorithms (Gauss-Seidel, Jacobi, Conjugate
gradient, …). A confusion between physical iteration (from the temporal advancement) and non-physical
iteration (from the iterative method) should be avoided.

Why implicit scheme? Notion of stability

Let consider 𝜀#W the difference between the numerical

approximation and the exact solution of the discrete
equations.
Exact solution of PDE
consistency
The goal is to know if 𝜀#WXK > 𝜀#W , meaning that the
numerical resolution leads to numerical instability.
Exact solution of discret (linear)
equations
Assuming 𝜀#W = 𝐴W 𝑒 jk#∆f .

stability The explicit formulation,

Numerical approximation
𝜙#WXK − 𝜙#W W
𝜙#XK W
− 𝜙#ZK
+𝑐 = 0,
LAX theorem: Consistency + Stability = Convergence
Δ𝑡 2Δ𝑥
leads to
=
𝜀#WXK Δ𝑡 = =
= 1 + 𝑐 sin 𝑘∆𝑥 > 1.
𝜀#W 2Δ𝑥
nd
The 1st order explicit formulation in time is always unstable (for advection equation, with 2 order central
scheme for spatial derivative)

The implicit formulation,

𝜙#WXK − 𝜙#W WXK
𝜙#XK WXK
− 𝜙#ZK
+𝑐 = 0,
Δ𝑡 2Δ𝑥
leads to
 =
𝜀#WXK 1
= < 1.
𝜀#W Δ𝑡 = =
1+ 𝑐 sin 𝑘∆𝑥
2Δ𝑥
nd
The 1st order implicit formulation in time is always stable (for advection equation, with 2 order central
scheme for spatial derivative). The major interest of implicit schemes is the larger range of stability.

Stability for compressible and incompressible flows resolution

For compressible flows, the stability condition is given by the CFL condition:
Δ𝑥
Δ𝑡 ≤ 𝐾
𝑢 +𝑎
with 𝐾 a coefficient depending of the numerical scheme used and 𝑎 the speed of sound (appearing because
the pressure waves have to be solved in compressible formulation).

If a compressible solver is used to solve incompressible case, Δ𝑡 → 0 because 𝑎 → ∞ in this case (impossible to
solve pressure waves because these waves are propagated instantaneously).

Specific treatment of pressure in incompressible solver: the fractional step method

The momentum conservation (Navier-Stokes equation) in incompressible case writes:

𝜕𝑢# 𝑢#WXK − 𝑢#W 𝜕𝑢# 𝑢7 𝜕 = 𝑢# 1 𝜕𝑝
= =− +𝜈 −
𝜕𝑡 Δ𝑡 𝜕𝑥7 𝜕𝑥7 𝜕𝑥7 𝜌 𝜕𝑥#
Remark: 𝑢# is the velocity component in the direction 𝑥# (the index 𝑖 does not means the spatial discretization
in this case)

step 1: Find the sub-step velocity field 𝑢#∗ , such as:

𝑢#∗ − 𝑢#W 𝜕𝑢# 𝑢7 𝜕 = 𝑢#
=− +𝜈
Δ𝑡 𝜕𝑥7 𝜕𝑥7 𝜕𝑥7

step 2: Find the pressure 𝑝, as a solution of the Poisson equation

𝜕𝑢#∗ Δ𝑡 𝜕 = 𝑝
=
𝜕𝑥# 𝜌 𝜕𝑥# 𝜕𝑥#

step 3: Finally, find 𝑢#WXK such as,

Δ𝑡 𝜕𝑝
𝑢#WXK = 𝑢#∗ −
𝜌 𝜕𝑥#

Remark: other methods exist.