Computational Fluid Dynamics

Computational Fluid Dynamics

http://www.nd.edu/~gtryggva/CFD-Course/! The Euler Equations!

The Euler ⎛ρ ⎞ ⎛ ρu ⎞
equations for ∂⎜ ⎟ ∂ ⎜ 2 ⎟
1D flow:! ⎜ ρu ⎟ + ⎜ ρu + p ⎟=0
∂t ⎜ ⎟ ∂x ⎜ ⎟
⎝ ρE ⎠ ⎝ ρu ( E + p / ρ ) ⎠
Solving the
E = e + u 2 /2
Euler where!
Define! h = e + p / ρ; H = h + u 2 / 2 = E + p / ρ;
Equations! Ideal Gas:! p = ρRT; e = e(T ); c v = de /dT

Grétar Tryggvason ! h = h (T ); c p = dh /dT

Spring 2011!
R = c p − c v ; γ = c p /c v ; p = (γ −1) ρe

Computational Fluid Dynamics

Computational Fluid Dynamics

The Euler Equations! The Euler Equations!

⎛ ρu Predictor-corrector method: Limiting the variables!

The Euler ⎛ρ ⎞ ⎞
equations for ∂⎜ ⎟ ∂ ⎜ 2 ⎟ Δt n
1D flow:! ⎜ ρu ⎟ + ⎜ ρu + p
∂t ⎜ ⎟ ∂x ⎜
⎟=0
⎟
n +1/ 2
Predictor step! f j = f jn −
2h
(F j +1/ 2 − F j−1/
n
2)

⎝ ρE ⎠ ⎝ ρu ( E + p / ρ ) ⎠
Ψ (fj )
1 L n +1/ 2
Variables! f jL+1/ 2 = f jn +1/ 2 + n +1/ 2
− f j−1
2
or! ∂ ∂
f j +1/ 2 = f j +1 − Ψ ( f j +1 − f j )
f+ F=0 R n +1/ 2 1 R n +1/ 2 n +1/ 2
∂t ∂x 2

( )
where!
2 = F ( f ) j +1/ 2 , ( f ) j +1/ 2
+1/ 2 n +1/ 2 n +1/ 2
Find:! F jn+1/ L R

⎛ρ ⎞ ⎛ ρu ⎞
⎜ ⎟
f = ⎜ ρu ⎟ ; F = ⎜ ρu 2 + p Δt n +1/ 2
⎜ ⎟
⎜⎝ ρ E ⎟⎠
⎟
⎜⎝ ρu(E + p / ρ )⎟⎠
Final step! f jn +1 = f jn −
h
(F j +1/ 2 − F j−1/ 2 )
n +1/ 2

5! 6!
Computational Fluid Dynamics
Computational Fluid Dynamics

The Euler Equations! The Euler Equations!

The Euler equations can be solved using the flux Although the Riemann problem can be solved for the
limited high order methods described earlier by finding general case when the fluids are not stationary, the
the fluxes using solutions to the Riemann problem! solution is expensive (involving solving a nonlinear
equation for the pressure ratio across the shock) and
+1/ 2
F jn+1/ (
2 = F (f )
L n +1/ 2

j +1/ 2
,( f R )
n +1/ 2

j +1/ 2 ) the only information that we need from the solution is

the flux across the cell boundary.!

In principle we can solve this problem, the Riemann Therefore, usually we use approximate Riemann
problem, exactly by assuming constant states and then solvers. Different approximations are possible but
integrate the fluxes over the time step. In the last generally the result is a complex process. Here we will
lecture we did so for the case when the fluids are only outline it briefly!
initially at rest.!
 7! 8!
Computational Fluid Dynamics
Computational Fluid Dynamics

The Euler Equations! The Euler Equations!

The Roe approximate Riemann solver was one of the Roe introduced the following matrix A!
first method to compute the fluxes in a “simpler” way.!
⎡ 0 1 0 ⎤
It is based on approximating the Euler equation by a ⎢ ⎥
γ −3 2
linear equation ! ⎢
A= ⎢ 2
u (3 − γ ) u γ −1 ⎥
⎥
∂f ∂f ⎢ ⎥
+A =0 γ − 1
∂t ∂x ⎢
⎢⎣ 2
u − uH
3
H − γ −1 u2( ) γu ⎥
⎥⎦
whose fluxes can be found analytically. The
where the average variables are defined by!
linearization is done in such a way that the correct
wave speed is preserved. The jump in the solution
across a shock can be written in terms of the ρ L uL + ρ R uR ρL H L + ρR H R
u= ; H=
eigenvectors:! 3 ρL + ρR ρL + ρR
f R − f L = ∑ α p Rp
p =1

9! 10!
Computational Fluid Dynamics
Computational Fluid Dynamics

The Euler Equations! The Euler Equations!

The fluxes are given by! The alphas are given by!
⎡ ⎤⎡
⎡ α ⎤ ⎢ u ⎛2 + γ −1 u ⎞ 1 ⎛ u⎞ γ −1 1 ⎤
Fi+1/ 2 = (
1 L 1
)
F + F R − ∑ λ p α p Rp
3
⎢ 1 ⎥ ⎢ 4c ⎜⎝ c ⎟⎠
( ) −
2c ⎜⎝
(
1+ γ − 1 ⎟
c⎠
) 2 c2
⎥ ⎢ Δρ j +1/ 2
⎥⎢
⎥
2 2 p =1 ⎢ ⎥ ⎢ ⎥
⎥
⎢ ⎥ ⎢ γ −1u2 u ⎢ ⎥
where the eigenvalues are the same as for the original! ⎢ α2 ⎥ = ⎢ 1−
2 c2
(γ − 1) cu 2 ( c ⎥⎢
)
− γ − 1 2 ⎥ ⎢ Δρu j +1/ 2 ⎥
⎥
⎢ ⎥ ⎢ ⎥
⎢ ⎥ ⎢ u ⎛ u⎞ 1 ⎛ u⎞ γ −1 1 ⎥⎢ ⎥
λ1 = u − c ; λ2 = c ; λ3 = u + c ; ⎢⎣ α 3 ⎥⎦ ⎢ − ⎜ 2 − γ − 1 ⎟
4c ⎝ c⎠
( ) 2c ⎜⎝
(
1− γ − 1 ⎟
c⎠
) ⎢
2 c 2 ⎥ ⎢⎣ Δρ E j +1/ 2
⎥
⎥⎦
⎣ ⎦
and the scaled eigenvectors are given by!
Where the jumps across the cell boundary are:!
⎛ 1 ⎞ Δρ j +1/ 2 = ρ Rj +1/ 2 − ρ Lj +1/ 2 ; Δρu j +1/ 2 = ρu Rj +1/ 2 − ρu Lj +1/ 2 ; Δρ E j +1/ 2 = ρ E jR+1/ 2 − ρ E jL+1/ 2
⎛ 1 ⎞ ⎛ 1 ⎞
⎜ ⎟
R1 = ⎜ ⎟; R = ⎜ u
u −c ⎟ ⎟ ; R1 = ⎜ u + c ⎟; and (again)!
⎜ 1
⎜ 1 2 ⎟ ⎜ ⎟
⎜⎝ ⎟
H −uc ⎠ ⎜⎝ H + u c ⎟⎠
⎜⎝ 2 u ⎟⎠ ρ L uL + ρ R uR ρL H L + ρR H R ⎛ 1 ⎞
u= ; H= ( ⎝ 2 ⎠
)
; c 2 = γ −1 ⎜ H − u2⎟
ρL + ρR ρL + ρR

11! 12!
Computational Fluid Dynamics
Computational Fluid Dynamics

The Euler Equations! The Euler Equations!

The fluxes 2rd order Runge-Kutta!

are given by! Time Integration!
( )
f (1) = f n + ΔtL f n
F 1
i+1/ 2
=
1
2 (( R
)
L
( ) 1
2 )
1
2
1
ρu j +1/ 2 + ρu j +1/ 2 − λ1α1 − λ2α 2 − λ3α 3
2 f (2) = f + ΔtL ( f )
(1) (1)

1⎛
) (( ) )
2
1
( )( ) ( ( ) 4rd order Runge-Kutta!
R R R
2
= ⎜ γ − 1 ρ E j +1/ 2 − γ − 3 ρu j +1/ 2 ρ j +1/ 2 +
= (f + f )
Fi+1/ 1
2
2⎝ 2 f n+1 n (2)

2 f (0) = f n
)(( ρu ) )
1 2
⎞
(γ − 1)( ρ E ) ( (ρ)
L L L
− γ −3

1 1
j +1/ 2

1
2 j +1/ 2 j +1/ 2 ⎟
⎠
3rd order Runge-Kutta!
f (1) = f (0) +
1
4
ΔtL f n ( )
( )
− λ1α1 u − c − λ2α 2u − λ3α 3 u + c ( ) ( )
2 2 2 f (1) = f n + ΔtL f n 1
f (2) = f (0) + ΔtL f (1)( )
3
Fi+1/ =
1
(( R
)
L
( )
ρu j +1/ 2 H R + ρu j +1/ 2 H L ) f (2) =
3 n 1 (1) 1
f + f + ΔtL f (1) ( )
3

( )
2
2 1
4 4 4 f = f + ΔtL f (2)
(3) (0)
1
( 1
) 1
− λ1α1 h − u c − λ2α 2u 2 − λ3α 3 h + u c ( ) 2
2 2 2 1 2 2
f n+1 = f n + f (2) + ΔtL f (2)
3 3 3
( ) f n+1 = f (0) + ΔtL f (2)( )
 13! 14!
Computational Fluid Dynamics
Computational Fluid Dynamics

The Euler Equations! The Euler Equations!

Results from:!
48 gridpoints!
PRINCIPLES OF COMPUTATIONAL !
FLUID DYNAMICS by P. Wesseling, !
Delft University of Technology, !
The Netherlands !

Matlab codes available at the website!

http://ta.twi.tudelft.nl/users/wesselin/cfdbook.html!

15! 16!
Computational Fluid Dynamics
Computational Fluid Dynamics

The Euler Equations! The Euler Equations!

The Roe approximate Riemann solver generally gives

148 gridpoints! well behaved results but it does allow for expansion
shocks in some cases. This can be corrected by the
so-called “entropy fix” (see Wesseling)!

17! 18!
Computational Fluid Dynamics
Computational Fluid Dynamics

The Euler Equations! The Euler Equations!

Several variants of similar methods have been proposed.!

The TVD property gives monotonic behavior for 1D
scalar conservation laws but generally not for One of the better known is the AUSM (Advection
nonlinear systems or multidimensional flows. However, Upstream Splitting Method) family of methods, originally
generally it is found that schemes that satisfy the TVD due to Liou and Steffen (1993), where the fluxes are split
property when applied to 1D scalar equation do well into convective and a pressure part. The method and its
for systems. When used for multidimensional flows by derivatives have been used widely for many aerospace
splitting, they are also generally found to be well and other applications, including multiphase flows. !
behaved.!
 Computational Fluid Dynamics
Computational Fluid Dynamics

The Euler Equations!

For the Euler equations:!

⎛ ρ⎞ ⎛ ρu ⎞ ⎛ ρv ⎞
⎜ ⎟ ⎜ ⎟ ⎜ ⎟
∂ ⎜ ρu ⎟ ∂ ⎜ ρu + p ⎟ ∂ ⎜ ρuv
2

+ + ⎟=0
The Euler ∂t ⎜ ρv ⎟ ∂x ⎜ ρuv ⎟ ∂y ⎜ ρv 2 + p ⎟
⎜ ⎟ ⎜ ⎟ ⎜ ⎟
⎝ ρE ⎠ ⎝ ρu(E + p / ρ)⎠ ⎝ ρv(E + p / ρ)⎠
Equations in two-
E = e + (u 2 + v 2 ) /2; p = (γ −1) ρe
dimensions! where!

In vector form!
∂f ∂F ∂G
+ + =0
∂t ∂x ∂y

Computational Fluid Dynamics

Computational Fluid Dynamics

The Euler Equations! The Euler Equations!
∂f ∂F + ∂F − ∂G + ∂G−
Flux vector splitting! + + + + =0
∂t ∂x ∂x ∂y ∂y
Two-dimensional problems are almost ⎛ 1 ⎞
⎜ ⎟
always solved by splitting where we ⎜ (γ −1)u ± 2c ⎟
apply one-dimensional methods to each ⎜ ⎟
ρ γ
F ± = ± (u ± c) ⎜ ⎟
2
direction in sequence!
4c ⎜ v ⎟
⎜ 2⎟
⎜ v 2 [(γ −1) u ± 2c ] ⎟ ⎛ 1 ⎞
⎜ ⎟
⎜2 + 2(γ 2 −1) ⎟⎠ ⎜ u ⎟
⎝
ρ ⎜ (γ −1)v ± 2c ⎟
±
(v ± c ) ⎜ ⎟
2
G =±
4c ⎜ γ ⎟
⎜ u 2 [(γ −1)v ± 2c ] 2 ⎟
⎜ ⎟
⎜2 + 2(γ 2 −1) ⎟⎠
⎝

Computational Fluid Dynamics

Computational Fluid Dynamics

The Euler Equations! The Euler Equations!

Results from a
Computational Shock! h=1/30

review paper:! h

domain!
P. Woodward and P.
Colella:!
The Numerical h=1/60

Simulation of two-
Dimensional Fluid
Flow with Strong 0

Shocks. ! Given inflow! Outflow!
1!
J. Comput. Phys. 3! h=1/120

54 (1984), 115-173.! wall!
Godunovʼs method! MacCormacʼs method!
 26!
Computational Fluid Dynamics
Computational Fluid Dynamics

The Euler Equations!

h=1/30

h

Higher Order and more
h=1/60

recent methods!
0

h=1/120

MUSCL! PPM!

27! 28!
Computational Fluid Dynamics
Computational Fluid Dynamics

ENO/WENO! ENO/WENO!

Beyond linear: Reconstruction of higher order Constructing an interpolation polynomial from the cell
approximations for the function in each cell averages: For anything higher than second order (linear)
(ENO and WENO). ! the problem is that the average value in the cell is not
equal to the value at the center. !
The critical step in the methods discussed so far
is the construction of a linear slope in each cell To get around this we look at the primitive function:!

( ) ∫ v (ξ ) d ξ
and the limitation of this slope to prevent x
V x ≡
oscillations. For higher order methods, a higher −∞
order profile needs to be constructed!
The lower bound is arbitrary and can be replaced!

j-2 j-1 j j+1 j+2

j-3/2 j-1/2 j+1/2 j+3/2

29! 30!
Computational Fluid Dynamics
Computational Fluid Dynamics

ENO/WENO! ENO/WENO!

Since this is the integral over the cells, the discrete version A polynomial interpolating the edge values is given by P(x)
is exact at the cell boundaries! and we denote its derivative by p(x)

i i

( ) ∑∫ () ∑
x j+1/2
V xi+1/ 2 = v ξ dξ = f i Δx p(x) = P'(x)
x j−1/2
j = −∞ j = −∞

Then it can be shown that!

( ) ( )
V xi−3/ 2 ; V xi−1/ 2 ; V xi+1/ 2 ; V xi+3/ 2 ; ( ) ( )
() ()
x j+1/2 x j+1/2
∫ x j−1/2
p ξ dξ = ∫ x j−1/2
P' ξ dξ = P(x j +1/ 2 ) − P(x j −1/ 2 )

()
x j+1/2
= V (x j +1/ 2 ) − V (x j −1/ 2 ) = ∫ x j−1/2
v ξ dξ = f i Δx

That is, the integral of p(x) over the cell is equal to the cell
j-2 j-1 j j+1 j+2
average fi

j-3/2 j-1/2 j+1/2 j+3/2

 31! 32!
Computational Fluid Dynamics
Computational Fluid Dynamics

ENO/WENO! ENO/WENO!

Thus, p(x) gives the correct average value in

each cell and the integrated value gives the
exact values of the primitive function at the cell
boundaries. ! j-2 j-1 j j+1 j+2

j-3/2 j-1/2 j+1/2 j+3/2

We need to write down a polynomial P(x) that
interpolates the values of the primitive function The interpolation polynomial is often taken to be the
of the cell boundaries and then differentiate Lagrangian Polynomial !
this polynomial to get p(x), which lets us
compute the variables at the cell boundary! k k
x − xi− r + l −1/ 2
P(x) = ∑ V (xi− r + m−1/ 2 )∏
m=0 l =0 xi− r −1/ 2 − xi− r + l −1/ 2
l≠m

Where r determines where we start and k is the order!

33!
Computational Fluid Dynamics
Computational Fluid Dynamics

ENO/WENO! ENO/WENO!

Example: Second order ENO!

f j −1/2 f j +1/2
j-2 j-1 j j+1 j+2

1. Construct left and right slopes f j −1 fj
by connecting the average
j-3/2 j-1/2 j+1/2 j+3/2
values in adjacent cells!

The question is now which point we select. We start by 2. Select the downstream flux by U>0! f j +1
interpolating over one cell (linear). To add one point we using the smaller slope!
can add either the point to the left or the right. In ENO we
select the points based on the minimum absolute value of
the divided differences of the function values! Δf j+ = f j +1 − f j Δf j− = f j − f j −1 j-1 ! j! j+1 !

⎪ j 2 j (
⎧ f + 1 amin Δf + , Δf − , if
j ) 1
2 (u + u ) > 0
j j +1
f j +1/2 = ⎨
⎩
1 +
( −
⎪ f j − 2 amin Δf j +1 , Δf j +1 , if ) 1
2(u + u ) < 0
j j +1

Computational Fluid Dynamics

Computational Fluid Dynamics

ENO/WENO! ENO/WENO!
Second order ENO scheme for the linear advection equation! Second order ENO scheme for the linear advection equation!
∂f ∂f ⎧ a,
+u =0 ⎪ a < b ∂f ∂f
∂t ∂x amin ( a,b ) = ⎨ +u =0
b, b ≤ a ∂t ∂x
Δt n n ⎩⎪
* n

h
(
f j = f j − u j f j +1/2 − f j −1/2
n
)
Blue: 2nd order ENO!
f jn +1 = f jn −
Δt 1 n n
h 2
( ( )
u j f j +1/2 − f jn−1/2 + u *j f j*+1/2 − f j*−1/2 ( ))
Red: 1st Upwind!
⎧ f + 1 amin Δf + , Δf − , if
⎪ j 2 j j ( ) 1
2 (u + u ) > 0
j j +1
f j +1/2 =⎨
⎩
1 +
(−
⎪ f j − 2 amin Δf j +1 , Δf j +1 , if ) 1
2(u + u ) < 0
j j +1

+ −
Δf = f j +1 − f j
j Δf = f j − f j −1
j
 37! 38!
Computational Fluid Dynamics
Computational Fluid Dynamics

ENO/WENO!

Increasingly we see methods developed for the

inviscid Euler equation with shocks being used for the
j-2 j-1 j j+1 j+2

j-3/2 j-1/2 j+1/2 j+3/2
advection part of the Navier-Stokes solvers. !
In the WENO (weighted essentially non-oscillating) scheme
we use all the points but weigh the contribution of each
according to a smoothness criteria. High-order WENO
represents the current state-of-the-art in computing of flows
with sharp interfaces!

For a recent reference on WENO see, for example:!

C-W Shu. High Order Weighted Essential Nonoscillatory
Schemes for Convection Dominated Problems. SIAM
Review, Vol. 51 (2009), 82-126.!

39! 40!
Computational Fluid Dynamics
Computational Fluid Dynamics

CIP-gradient augmentation!

The CIP (Constrained Interpolation Polynomial) !

Method (Yabe)!
In addition to advecting the marker function f, its derivative is
advected by fitting a third order polynomial through the
Other function and its derivatives.!
f and g given!
New f and g!

Start with!
Approaches! Introduce!
In 1D, the advection of
the derivative is given by!

Therefore, the derivative is translated with velocity u,

just as the function. In 2D splitting is used to separate
translation and deformation!

41!
Computational Fluid Dynamics
Computational Fluid Dynamics

CIP-gradient augmentation!

Enormous progress has been made in solution

techniques for hyperbolic systems with shocks in the last
The CIP method results
twenty years. Advanced methods are now able to
in very accurate
resolve complex shocks within a grid space or two, even
advection and for a
in multidimensional situations for a large range of
sharp interface it greatly
governing parameters and physical complexity.!
reduces overshoots, but
does not eliminate them Here, we have only examined relatively elementary
completely! aspects of methods for hyperbolic systems, but this short
introduction should have taught you methods to solve
such systems and introduced you to literature.!
 43!
Computational Fluid Dynamics

A few words about the project:!

I asked you to select your own parameters. In doing so

keep in mind that while Nature usually presents us with
problems that have multiphysics problems with a large
range of scales and extreme values of the physical
properties, in CFD we prefer single physics, small range of
scales and O(1) values for the physical properties. In the
project you are simply showing that the code works and
looking at trends. Using Matlab you are limited to modest
resolution and run-time so you should select you
parameters commensurate with the resolution you are
using! !