Remark: foils with

1
Computer Fluid Dynamics E181107
2181106
CFD7
Solvers, schemes SIMPLEx, upwind,
Remark foils with black background could be
skipped, they are aimed to the more advanced
courses
Rudolf Žitný, Ústav procesní a zpracovatelské
techniky CVUT FS 2010
2
FINITE VOLUME METHOD
CFD7
Principle of the Finite Volume Method is in
application of Gauss theorem upon integral of the
solved transport equation. Thus the volumetric
integral inside the fluid element is replaced by
surface integral, evaluated as the sum of
integrals at faces of the fluid element.
3
FINITE VOLUME METHOD
CFD7
FINITE CONTROL VOLUME ?x Fluid element dx
Finite size
Infinitely small size
4
FVM diffusion problem 1D
CFD7
CV - control volume A surface of CV
Gauss theorem
Diffusive flux
1D case
5
FVM diffusion problem 1D
CFD7
Overall flux from the west side
Linearization of source term
aE
aW
aP

• Check properties
• aP aWaE for S0 (without sources)
• All coefficients are positive

6
Tutorial sphere drying
CFD7
Enthalpy of evaporation
Heat transfer from air at temperature To
W P E
Mass transfer coefficient
Diffusion coefficient
Internal control volumes
TR temperature at surface of sphere must be
extrapolated (using TW and TP)
?x
?x/2
W P R
7
FVM diffusion problem 2D
CFD7
N
W P E
A
S
Boundary conditions
Example insulation q0

• First kind (Dirichlet BC)
• Second kind (Neumann BC)
• Third kind (Newtons BC)

Example Finite thermal resistance (heat transfer
coefficient ?)
8
FVM convection diffusion 1D
CFD7
CV - control volume A surface of CV
Gauss theorem
1D case
9
FVM convection diffusion 1D
CFD7
F-mass flux (through faces A) D-diffusion
conductance
Continuity equation (mass flowrate conservation)
FeFw
Relative importance of convective transport is
characterised by Peclet number of cell (of
control volume, because Pe depends upon size of
cell)
Thermal conductivity
Prandtl
Binary diffusion coef
Schmidt
10
FVM convection diffusion 1D
CFD7

• Methods differ in the way how the unknown
  transported values at the control volume faces
  (?w , ?e) are calculated (different interpolation
  techniques)
• Result can be always expressed in the form
• Properties of resulting schemes should be
  evaluated
• Conservativeness
• Boundedness (positivity of coefficients)
• Transportivity (schemes should depend upon the
  Peclet number of cell)

11
FVM central scheme 1D
CFD7

• Conservativeness (yes)
• Boundedness (positivity of coefficients)
• Transportivness (no)

Very fine mesh is necessary
12
FVM upwind 1st order 1D
CFD7

• Conservativeness (yes)
• Boundedness (positivity of coefficients) YES for
  any Pe
• Transportivness (YES)

But at a prize of decreased order of accuracy
(only 1st order!!)
13
FVM Upwind FLUENT
CFD7
Upwind of the first order Upwind of the second
order
?1
?f
?0
14
FVM hybrid upwind Spalding 1972 1D
CFD7

• Central scheme for Pelt2
• Upwind with diffusion supressed for Pegt2

• Conservativeness (yes)
• Boundedness (positivity of coefficients) YES for
  any Pe
• Transportivness (YES)

But at a prize of decreased order of accuracy at
high Pe
15
FVM power law Patankar 1980 1D
CFD7

• Polynomial flow for Pelt10
• Upwind with zero diffusion for Pegt10

• Conservativeness (yes)
• Boundedness (positivity of coefficients) YES for
  any Pe
• Transportivness (YES)

16
FVM QUICK Leonard 1979 1D
CFD7
Quadratic Upwind Interpolation more nodal
points necessary
w e
WW W P E
EE

• Third order of accuracy very low numerical
  diffusion comparing with other schemes
• Scheme is not bounded (stability problems are
  frequently encountered)

?e1 for Fegt0 else ?e0.
?w1 for Fwgt0 else ?w0.
17
FVM exponential Patankar 1980 1D
CFD7

• Exact solution for constant velocity and
  diffusion coefficient ?

• Conservativeness (yes)
• Boundedness (positivity of coefficients) YES for
  any Pe
• Transportivness (YES)

18
FVM exponential FLUENT
CFD7
Called power law model in Fluent

• Exact solution of

Projection of velocity to the direction r
?1
r
?0

• For Peltlt1 the exponential profile is reduced to
  linear profile

19
FVM exponential proof
CFD7

• Exact solution for constant velocity and
  diffusion coefficient ?
• Assume constant mass flux F and constant PeF/D.

Analytical solution
Boundary conditions
Identified constants
Exact solution at central point
20
FVM central scheme with modified viscosity
CFD7
Almost any scheme (with the exception of QUICK)
can be converted to central scheme introducing
modified transport coefficient ? Modified
diffusion conductances are and .
Central scheme
Upwind scheme
Exponential scheme
21
FVM dispersion / diffusion
CFD7

• There are two basic errors caused by
  approximation of convective terms
• False diffusion (or numerical diffusion)
  artificial smearing of jumps, discontinuities.
  Neglected terms in the Taylor series expansion of
  first derivatives (convective terms) represent in
  fact the terms with second derivatives (diffusion
  terms). So when using low order formula for
  convection terms (for example upwind of the first
  order), their inaccuracy is manifested by the
  artificial increase of diffusive terms (e.g. by
  increase of viscosity). The false diffusion
  effect is important first of all in the case when
  the flow direction is not aligned with mesh, see
  example
• Dispersion (or aliasing) causing overshoots,
  artificial oscillations. Dispersion means that
  different components of Fourier expansion (of
  numerical solution) move with different
  velocities, for example shorter wavelengths move
  slower than the velocity of flow.

Numerical dispersion
Numerical diffusion
1
0 - boundary condition, all values are zero in
the right triangle (without diffusion)
22
FVM unstructural mesh generation
CFD7
Finite Volume Method Finite Element Method
Voronoi polygons Delaunay triangulation each side
of Voronoi polygon has a node associated with the
node i and these nodes are vertices of Delaunay
triangles
Voronoi control volume associated with node i.
Set of points x,y having distance to i less than
the distance to any other node j.
Property of Delaunay triangles Circumcenter
remains inside the Delaunay triangle
i
j
23
FVM structural mesh generation
CFD7
Laplace equation solvers
24
FVM Navier Stokes equations
CFD7
Different methods are used for numerical solution
of Navier Stokes equations at Compressible flows
(finite speed of sound, existence of velocity
discontinuities at shocks) Incompressible
flows (infinite speed of sound)
Explicit methods (e.g. MacCormax) but also
implicit methods (Beam-Warming scheme) are used.
These methods will not be discussed here.
Vorticity and stream function methods (pressure
is eliminated from NS equations). These methods
are suitable especially for 2D flows. Primitive
variable approach (formulation in terms of
velocities and pressure). Pressure represents a
continuity equation constraint. Only these
methods (pressure corrections methods SIMPLE,
SIMPLER, SIMPLEC and PISO) using staggered and
collocated grid will be discussed in this lecture.
Beckman
25
FVM Navier Stokes equations
CFD7
Example Steady state and 2D (velocities u,v and
pressure p are unknown functions)
This is equation for u
This is equation for v
This is equation for p?
But pressure p is not in the continuity equation!
Solution of this problem is in SIMPLE methods,
described later
26
FVM checkerboard pattern
CFD7
Other problem is called checkerboard pattern This
nonuniform distribution of pressure has no effect
upon NS equations
The problem appears as soon as the u,v,p values
are concentrated in the same control volume
W P E
27
FVM staggered grid
CFD7
Different control volumes for different equations
Control volume for momentum balance in y-direction
Control volume for continuity equation,
temperature, concentrations
Control volume for momentum balance in x-direction
28
FVM staggered grid
CFD7
Location of velocities and pressure in staggered
grid
Control volume for momentum balance in y-direction
j2 j1 j j-1
Control volume for continuity equation,
temperature, concentrations
J1 J J-1 J-2
Pressure I,J U i,J V
I,j
u
p
v
Control volume for momentum balance in x-direction
I-2 I-1 I I1 I2
i-1 i i1 i2
29
FVM continuity equation
CFD7
j2 j1 j j-1
Control volume for continuity equation,
temperature, concentrations
?x
J1 J J-1 J-2
Pressure I,J U i,J V
I,j
p
u
v
I-2 I-1 I I1 I2
i-1 i i1 i2
30
FVM momentum x
CFD7
j2 j1 j j-1
J1 J J-1 J-2
Pressure I,J U i,J V
I,j
u
p
p
Control volume for momentum balance in x-direction
I-2 I-1 I I1 I2
i-1 i i1 i2
31
FVM momentum y
CFD7
Control volume for momentum balance in y-direction
j2 j1 j j-1
J1 J J-1 J-2
Pressure I,J U i,J V
I,j
p
v
p
I-2 I-1 I I1 I2
i-1 i i1 i2
32
FVM SIMPLE step 1 (velocities)
CFD7
Input Approximation of pressure p
33
FVM SIMPLE step 2 (pressure correction)
CFD7
true solution
Neglect!!!
subtract
Substitute uu and vv into continuity
equation
Solve equations for pressure corrections
34
FVM SIMPLE step 3 (update p,u,v)
CFD7
Only these new pressures are necessary for next
iteration
Continue with the improved pressures to the step
1
35
FVM SIMPLEC (SIMPLE corrected)
CFD7
Neglect!!! Good estimate as soon as neighbours
are not far from center
Improved diJ
36
FVM SIMPLER (SIMPLE Revised)
CFD7
Idea calculate pressure distribution directly
from the Poissons equation
37
FVM Rhie Chow (Fluent)
CFD7
The way how to avoid staggering (difficult
implementation in unstructured meshes) was
suggested in the paper
Rhie C.M., Chow W.L. Numerical study of the
turbulent flow past an airfoil with trailing edge
separation. AIAA Journal, Vol.21, No.11, 1983
This technique is called Rhie-Chow interpolation.
38
FVM Rhie Chow (Fluent)
CFD7
Rhie C.M., Chow W.L. Numerical study of the
turbulent flow past an airfoil with trailing edge
separation. AIAA Journal, Vol.21, No.11, 1983
How to interpolate ue velocity at a control
volume face
39
Fluent - solvers
CFD7
Segregated solver
Coupled solver
?, ?, cp,
?, ?, cp,
u?
v?
w?
p? (SIMPLE)
CFL Courant Friedrichs Lewy criterion
k,?,?
T,k,?,?
Nonstationary problem always solved by implicit
Euler
Nonstationary problem solved by implicit Euler
(CFL5), explicit Euler (CFL1) or Runge Kutta
Algebraic eq. by multigrid AMG, or FAS (Full
Aproximation Storage)
Algebraic eq. by multigrid AMG
40
FVM solvers
CFD7
Dalí
41
FVM solvers
CFD7
Solution of linear algebraic equation system Axb

• TDMA Gauss elimination tridiagonal matrix
  (obvious choice for 1D problems, but suitable for
  2D and 3D problems too, iteratively along the
  x,y,z grid lines method of alternating
  directions ADI)
• PDMA pentadiagonal matrix (suitable for QUICK),
  Fortran version
• CGM Conjugated Gradient Method (iterative
  method each iteration calculates increment of
  vector ? in the direction of gradient of
  minimised function square of residual vector)
• Multigrid method (Fluent)

42
Solver TDMA tridiagonal system MATLAB
CFD7
function x TDMAsolver(a,b,c,d) a, b, c, and d
are the column vectors for the compressed
tridiagonal matrix n length(b) n is the
number of rows Modify the first-row
coefficients c(1) c(1) / b(1) Division by
zero risk. d(1) d(1) / b(1) Division by zero
would imply a singular matrix. for i 2n
id 1 / (b(i) - c(i-1) a(i)) Division by
zero risk. c(i) c(i) id
Last value calculated is redundant. d(i)
(d(i) - d(i-1) a(i)) id end Now back
substitute. x(n) d(n) for i n-1-11
x(i) d(i) - c(i) x(i 1) end end
Forward step (elimination entries bellow diagonal)
43
Solver TDMA applied to 2D problems (ADI)
CFD7
Y-implicit step (repeated application of TDMA for
10 equations, proceeds from left to right)
X-implicit step (proceeds bottom up)
y
x
44
Solver PDMA pentagonal system MATLAB
CFD7
Backsubstitution L'xc x(N)c(N)
x(N-1)c(N-1)-gamma(N-1)x(N) for
kN-2-11 x(k)c(k)-gamma(k)x(k1)-delta
(k)x(k2) end else
Non-symmetric Matrix Scheme ddiag(A)
ediag(A,1) fdiag(A,2)
h0diag(A,-1) g00diag(A,-2)
alphazeros(N,1) gamzeros(N-1,1)
deltazeros(N-2,1) betzeros(N,1)
czeros(N,1) zzeros(N,1)
alpha(1)d(1) gam(1)e(1)/alpha(1)
delta(1)f(1)/alpha(1) bet(2)h(2)
alpha(2)d(2)-bet(2)gam(1) gam(2)(
e(2)-bet(2)delta(1) )/alpha(2)
delta(2)f(2)/alpha(2) for k3N-2
bet(k)h(k)-g(k)gam(k-2)
alpha(k)d(k)-g(k)delta(k-2)-bet(k)gam(k-1)
gam(k)( e(k)-bet(k)delta(k-1) )/alpha(k)
delta(k)f(k)/alpha(k) end
bet(N-1)h(N-1)-g(N-1)gam(N-3)
alpha(N-1)d(N-1)-g(N-1)delta(N-3)-bet(N-1)gam(N
-2) gam(N-1)( e(N-1)-bet(N-1)delta(N-2)
)/alpha(N-1) bet(N)h(N)-g(N)gam(N-2)
alpha(N)d(N)-g(N)delta(N-2)-bet(N)gam(N-1)
Update bLc c(1)b(1)/alpha(1)
c(2)(b(2)-bet(2)c(1))/alpha(2) for
k3N c(k)( b(k)-g(k)c(k-2)-bet(k)c(k-1
) )/alpha(k) end
function xpentsolve(A,b) Solve a
pentadiagonal system Axb where A is a strongly
nonsingular matrix Reference G.
Engeln-Muellges, F. Uhlig, "Numerical Algorithms
with C" Chapter 4.
Springer-Verlag Berlin (1996) Written by Greg
von Winckel 3/15/04 Contact gregvw_at_chtm.unm.edu
M,Nsize(A) xzeros(N,1) Check for
symmetry if AA' Symmetric Matrix Scheme
Extract bands ddiag(A) fdiag(A,1)
ediag(A,2) alphazeros(N,1)
gammazeros(N-1,1) deltazeros(N-2,1)
czeros(N,1) zzeros(N,1) Factor
ALDL' alpha(1)d(1) gamma(1)f(1)/alpha(
1) delta(1)e(1)/alpha(1)
alpha(2)d(2)-f(1)gamma(1)
gamma(2)(f(2)-e(1)gamma(1))/alpha(2)
delta(2)e(2)/alpha(2) for k3N-2
alpha(k)d(k)-e(k-2)delta(k-2)-alpha(k-1)gamma
(k-1)2 gamma(k)(f(k)-e(k-1)gamma(k-1))
/alpha(k) delta(k)e(k)/alpha(k)
end alpha(N-1)d(N-1)-e(N-3)delta(N-3)-al
pha(N-2)gamma(N-2)2 gamma(N-1)(f(N-1)-e(N-
2)gamma(N-2))/alpha(N-1) alpha(N)d(N)-e(N-2
)delta(N-2)-alpha(N-1)gamma(N-1)2
Update Lxb, Dcz z(1)b(1)
z(2)b(2)-gamma(1)z(1) for k3N
z(k)b(k)-gamma(k-1)z(k-1)-delta(k-2)z(k-2)
end cz./alpha
Back substitution Rxc x(N)c(N)
x(N-1)c(N-1)-gam(N-1)x(N) for kN-2-11
x(k)c(k)-gam(k)x(k1)-delta(k)x(k2)
end end
45
Solver CGM conjugated gradient GNU
CFD7
Solution of linear algebraic equation system Axb
Minimisation of quadratic function
function x conjgrad(A,b,x) rb-Ax
pr rsoldr'r for i1size(A)(1)
ApAp alpharsold/(p'Ap)
xxalphap rr-alphaAp
rsnewr'r if sqrt(rsnew)lt1e-10
break end
prrsnew/rsoldp rsoldrsnew
end end
Residual vector in k-th iteration
Vector of increment conjugated with previous
increments
46
Solver MULTIGRID
CFD7
Solution on a rough grid takes into account very
quickly long waves (distant boundaries etc), that
is refined on a finer grid.
Rough grid (small wave numbers)
Fine grid (large wave numbers details)
Interpolation from rough to fine grid
Averaging from fine to rough grid
47
Unsteady flows
CFD7
Discretisation like in the finite difference
methods discussed previously
Example Temperature field
EXPLICIT scheme
IMPLICIT scheme
CRANK NICHOLSON scheme
First order accuracy in time First
order accuracy in time Second order
accuracy in time
Stable only if
Unconditionally stable and bounded
Unconditionally stable, but bounded
only if
48
FVM Boundary conditions
CFD7

• Boundary conditions are specified at faces of
  cells (not at grid points)
• INLET specified u,v,w,T,k,? (but not
  pressure!)
• OUTLET nothing is specified (p is
  calculated from continuity eq.)
• PRESSURE only pressure is specified (not
  velocities)
• WALL zero velocities, kP, ?P
  calculated from the law of wall

P
yP
Recommended combinations for several outlets
Forbidden combinations for several outlets
there is no unique solution because no flowrate
is specified
49
Stream function-vorticity
CFD7
Introducing stream function and vorticity instead
of primitive variables (velocities and pressure)
eliminates the problem with the checkerboard
pattern. Continuity equation is exactly
satisfied. However the technique is used mostly
for 2D problems (it is sufficient to solve only 2
equations), because in 3D problems 3 components
of vorticity and 3 components of vector potential
must be solved (comparing with only 4 equations
when using primitive variable formulation of NS
equations).
Ghirlandaio
50
Stream function-vorticity
CFD7
The best method for solution of 2D problems of
incompressible flows
Pressure elimination
Introduction vorticity (z-component of vorticity
vector)
Velocities expressed in terms of scalar stream
function
Vorticity expressed in terms of stream function
51
Stream function-vorticity
CFD7
Three equations for u,v,p are replaced by two
equations for ?, ?. Advantage continuity
equation is exactly satisfied. Vorticity
transport is parabolic equation, Poisson equaion
for stream function is eliptic.
Lines of constant ? are streamlines (at wall
?const). The no-slip condition at wall
(prescribed velocities) must be reflected by
boundary condition for vorticity
How to do it, is shown in the next slide on
example of flow in a channel
52
Stream function-vorticity
CFD7
Stream function ?
? Increases from 0 to ?w for example ?uy
?N-1?N
M
y
H
N-1 N
2
1 2
x
1
?0
Vorticity
Zero vorticity at axis (follows from
definition) Prescribed velocity profile at inlet
u1(y) and v1(y)0
Vorticity at wall (prescribed velocity U, in this
case zero)
53
Example cavity flow 1/4
CFD7
Cavity flow. Lid of cavity is moving with a
constant velocity U. There is no in-out flow,
therefore stream function (proportional to
flowrate) is zero at boundary, at all walls of
cavity
Vorticity at moving lid
Vorticity at fixed walls
54
Example cavity flow 2/4
CFD7
FVM applied to vorticity transport equation
(upwind of the first order). Equidistant grid is
assumed ??x?y
Eliptic equation for stream function
This systém of algebraic equation (including
boundary conditions) is to be solved iteratively
(for example by alternating direction method ADI
described previously)
55
Example cavity flow 3/4
CFD7
Matlab tok v dutine. metoda vírivosti a
pokutové funkice. Upwind. Metoda strídavých
smeru. h0.01 um0.0001
visc1e-6 parametry site, a casovy krok
relax0.5 n21 dh/(n-1) d2d2
niter50 psizeros(n,n)
omegazeros(n,n) uzeros(n,n)
vzeros(n,n) azeros(n,1) bones(n,1)
czeros(n,1) rzeros(n,1) for
iter1niter for i1n
u(i,n)um end for i2n-1
for j2n-1
u(i,j)(psi(i,j1)-psi(i,j-1))/(2d)
v(i,j)(psi(i-1,j)-psi(i1,j))/(2d)
end end
stream function x-implicit for j2n-1 for
i2n-1 a(i)-1/d2b(i)4/d2c(i)-1/d2
r(i)(psi(i,j-1)psi(i,j1))/d2 omega(i,j)
end r(1)0r(n)0 pstridag(a,b,c,r,n)
for i2n-1 psi(i,j)(1-relax)psi(i,j)
relaxps(i) end end stream function
y-implicit for i2n-1 for j2n-1
a(j)-1/d2b(j)4/d2c(j)-1/d2r(j)(psi(i-1,j)p
si(i1,j))/d2 omega(i,j) end
r(1)0r(n)0 pstridag(a,b,c,r,n) for
j2n-1 psi(i,j)(1-relax)psi(i,j)relax
ps(j) end end vorticity boundary
conditions for i1n omega(i,n)relaxomega(i
,n)(1-relax)2/d2 (psi(i,n)-psi(i,n-1)-umd)
omega(i,1)relaxomega(i,1)(1-relax)2
(psi(i,1)-psi(i,2))/d2 omega(1,i)relaxomega
(1,i)(1-relax)2(psi(1,i)-psi(2,i))/d2
omega(n,i)relaxomega(n,i)(1-relax)2(psi(n,i)-
psi(n-1,i))/d2 end
vorticity x-implicit for j2n-1 for
i2n-1 upu(i,j)vpv(i,j)
a(i)-visc/d2-max(up,0)/d
b(i)4visc/d2(abs(up)abs(vp))/d
c(i)-visc/d2-max(-up,0)/d
r(i)omega(i,j-1)(visc/d2max(vp,0)/d)
omega(i,j1)(visc/d2max(-vp,0)/d) end
r(1)omega(1,j)r(n)omega(n,j)
pstridag(a,b,c,r,n) for i2n-1
omega(i,j)(1-relax)omega(i,j)relaxps(i)
end end vorticity y-implicit for i2n-1 for
j2n-1 upu(i,j)vpv(i,j)
a(j)-visc/d2-max(vp,0)/d
b(i)4visc/d2(abs(up)abs(vp))/d
c(j)-visc/d2-max(-vp,0)/d
r(j)omega(i-,j)(visc/d2max(up,0)/d)
omega(i1,j) (visc/d2max(-up,0)/d)
end r(1)omega(i,1)r(n)omega(i,n)
pstridag(a,b,c,r,n) for j2n-1
omega(i,j)(1-relax)omega(i,j)relaxps(j)
end end end
56
Example cavity flow 4/4
CFD7
Re1000
Re1
?
57
Example cavity flow FLUENT
CFD7
Re1
u
v