Chapter 11: Methods for dealing with complex geometries Ib hi Sezai Ibrahim S i Department of Mechanical Engineering Eastern Mediterranean University

Fall 2011-2012

Introduction When the solution domain is not rectangular Cartesian grids cannot be used. When the boundary does not coincide with the co-ordinate lines we can use stepwise approximation near the boundary. p pp y

Such an approximation is tedious and time consuming. consuming Introduces errors in the computation of wall shear stresses and heat transfer. Cells in the solid cylinder do not take part in the calculations Æ waste of computer storage. ME555 : Computational Fluid Dynamics

2

I. Sezai – Eastern Mediterranean University

1

Body fitted grids CFD methods for complex geometries: 1. Structured curvilinear grids (body fitted grids) 2. Unstructured grids

St t d curvilinear ili id are based b d on mapping i off the th Structured grids flow domain onto a computational domain. It is very difficult to find viable mappings when the geometry is complex. In such cases the domain may be sub-divided into several sub-regions or blocks. b i bl k Æ Block Bl k structured t t d grids. id For the most complex geometries unstructured grids are used. Most commercial packages use unstructured grids. ME555 : Computational Fluid Dynamics

3

I. Sezai – Eastern Mediterranean University

Body fitted grids An example of an orthogonal curvilinear mesh for calculating g flow around an aerofoil

Use of a non-orthogonal body-fitted grid arrangement for the prediction of flow over a half-cylinder ME555 : Computational Fluid Dynamics

4

I. Sezai – Eastern Mediterranean University

2

Cartesian vs. curvilinear grids Flow over a heat exchanger tube bank (only a part shown) Æ

Body fitted grids give better results with no waste of grids

(a) Cartesian grid using an approximated profile to represent cylindrical surfaces; (b) predicted flow pattern using a 40 X 15 Cartesian grid ME555 : Computational Fluid Dynamics

(a) Non-orthogonal body-fitted grid for the same problem; (b) predicted flow pattern using a 40 X 15 structured body-fitted grid

5

I. Sezai – Eastern Mediterranean University

Curvilinear grids - Difficulties The governing equations in curvilinear co-ordinate system are much more complex. The equations in 2-D are: Continuity

∂ ∂ ( ρU ) + ( ρV ) = 0 ∂ξ ∂η

x-momentum

J

y-momentum

J

⎧ ⎪ Metric coefficients ⎨ ⎪ ⎩ Jacobian of transformation Contravariant velocities

(11.1)

∂(ρu) ∂ ∂ ⎡ ρUu − μ Jg 11uξ ⎤⎦ + ⎡ ρVu − μ Jg 22uη ⎤⎦ + (11.2) ∂t ∂ξ ⎣ ∂η ⎣ ∂p ∂p ∂ ∂ 12 12 ⎡ − μ Jg uη ⎤⎦ + ⎡ − μ Jg uξ ⎤⎦ + JSu = − yη + yξ + ∂ξ ∂η ∂ξ ⎣ ∂η ⎣ ∂ ( ρ υ) ∂ ∂ ⎡ ρUυ − μ Jg 11υξ ⎤⎦ + ⎡ ρVυ − μ Jg 22 υη ⎤⎦ + (11.3) ∂t ∂ξ ⎣ ∂η ⎣ ∂p ∂p ∂ ∂ ⎡ − μ Jg 12 υη ⎤⎦ + ⎡ − μ Jg 12 υξ ⎤⎦ + JS υ = + xη − xξ + ∂ξ ∂η ∂ξ ⎣ ∂η ⎣

g12 = −( xξ xη + yξ yη ) / J 2 g = ( xη + yη ) / J 11

2

2

η

2

g 22 = ( xξ2 + yξ2 ) / J 2 J = xξ yη − xη yξ

U = yη u − xη υ

ME555 : Computational Fluid Dynamics

6

V = xξ υ − yξ u

ξ

I. Sezai – Eastern Mediterranean University

3

Curvilinear grids - difficulties Body-fitted grids are structured, so grid refinement is generally not local.

The refinement needed to resolve the boundary layers persists in the interior mesh, which represents a waste of computer storage. ME555 : Computational Fluid Dynamics

7

I. Sezai – Eastern Mediterranean University

Curvilinear grids - difficulties To generate curvilinear meshes it is necessary to map the physical geometry into a computational geometry.

Mapping of physical geometry to computational geometry in structured meshes.Æ

(a) physical grid in x-y co-ordinates

Mapping may be a tedious task if the domain is more complex. ME555 : Computational Fluid Dynamics

(b) the mapped structure for (a) in the computational domain 8

I. Sezai – Eastern Mediterranean University

4

Curvilinear grids - difficulties For complex geometries structured grid generation may be very difficult resulting in skewed cells, which can lead to stability problems for CFD solvers.

A structured non-orthogonal mesh for a pent-roof i.c. engine geometry ME555 : Computational Fluid Dynamics

9

I. Sezai – Eastern Mediterranean University

Curvilinear grids - difficulties Problems encountered with structured grids: Still difficult and time consuming to generate If the geometry cannot be readily mapped into a rectangle (in 2D) or parallellepiped (3D) this can result in skewed grid lines. Unnecessary grid resolutions can result when mapping is difficult Mapping is sometimes impossible with complex 3D geometries with internal objects/parts.

ME555 : Computational Fluid Dynamics

10

I. Sezai – Eastern Mediterranean University

5

Block-structured grids In block-structured grids the domain is sub-divided into regions, each of which has a structured mesh. Such meshes are more flexible than (single block) structured meshes meshes. The block structured approach allows the use of fine grids in regions where greater resolution is required.

ME555 : Computational Fluid Dynamics

11

I. Sezai – Eastern Mediterranean University

Block-structured grids

Block-structured mesh arrangement for an engine geometry, including inlet and exhaust ports, used in engine simulations with KIVA-3V ME555 : Computational Fluid Dynamics

12

I. Sezai – Eastern Mediterranean University

6

Unstructured grids There is no implicit structure of co-ordinate lines in unstructured grids. The mesh can easily be concentrated where necessary storage without wasting computer storage. Control volumes can be of any shape.

ME555 : Computational Fluid Dynamics

13

I. Sezai – Eastern Mediterranean University

Unstructured grids A mixture of different cell shapes can be used.

Unstructured mesh generation is fairly g straightforward. Complex domains can be meshed. Mesh refinement and adaption (in regions with high gradients) are much easier. ME555 : Computational Fluid Dynamics

14

I. Sezai – Eastern Mediterranean University

7

Discretisation in unstructured grids There are two ways of defining CV’s in unstructured grids: 1. Cell-centered control volumes 2. Vertex-centered control volumes.

(a) cell-centred control volumes

(b) vertex-centred control volumes

We will use cell-centered method. ME555 : Computational Fluid Dynamics

15

I. Sezai – Eastern Mediterranean University

Discretisation in unstructured grids The general transport equation (2.39) is

∂ ( ρφ ) + div( ρφ v ) = div(Γ grad φ ) + sφ ∂t Integrating g g Eq. q (11.4) ( ) over a 3-D control volume ∂ ( ρφ ) ∫CV ∂t dV + CV∫ div( ρφ v)dV =CV∫ div(Γ grad φ )dV + CV∫ sφ dV For a vector a Gauss’ divergence theorem states ∫ div adV =∫ n ⋅ adA CV

(11.4) (11.5)

(11.6)

A

n.a Æ component p of vector a in the direction of the outward unit vector n normal to infinitesimal surface element dA.

ME555 : Computational Fluid Dynamics

16

I. Sezai – Eastern Mediterranean University

8

Discretisation in unstructured grids Application of Gauss’ divergence theorem to equation (11.1) gives ⎞ ∂⎛ (11.7) ⎜ ∫ ρφ dV ⎟ + ∫ n ⋅ ( ρφ v )dA = ∫ n ⋅ (Γ grad φ )dA + ∫ sφ dV ∂t ⎝ CV A CV ⎠ A U i A =An, Using A where h A = Axi+A i Ayj+A j Azk, k A = |A| = ( Ax2 + Ay2 + Az2 )1/ 2 n = nxi+nyj+nzk (unit vector normal to surface) A = surface area vector A = magnitude of the surface area

⎞ ∂⎛ φ dV ⎟ + ∫ ( ρφφ v ) ⋅ dA = ∫ (Γ gradd φ ) ⋅ dA + ∫ sφ dV ⎜ ∫ ρφ ∂t ⎝ CV A A CV 4244 3 1442443 1 424 3 144 244⎠3 14 Transient term

Convection term

Diffusion term

(11 8) (11.8)

Source term

Let us discretize each term separately ME555 : Computational Fluid Dynamics

17

I. Sezai – Eastern Mediterranean University

Transient term Using first order time derivative, ⎞ ( ρφV ) P − ( ρφV ) nP−1 ∂⎛ ⎜ ∫ ρφ dV ⎟ ≈ ∂t ⎝ CV Δt ⎠

f2 N3

N2 N1

f3

f1

P

Superscript S i n–11 refers f to previous i time i step f N Terms with no superscript refer to present time. Cell volume VP can be calculated by applying Gauss theorem on the identity 1 = div(xi): 4

4

V = ∫ dV = ∫ div( xi )dV = ∫ xi ⋅ ndA ≈ ∑ x f Afx V

V

A

((11.9))

f

f → refers to cell faces. For the cell of node P the summation is over the faces f1, f2, f3 and f4. Afx = x-component of the cell face surface vector.

xf = x-component of the position vector of the center of face f ME555 : Computational Fluid Dynamics

18

I. Sezai – Eastern Mediterranean University

9

Volumes and surface areas A f = Af n = Afx i + Afy j + Afz k

(11.10)

Instead of xi, one can also use yj or zk, that is: V ≈ ∑ y f A fy ≈∑ z f A fz f

(11.11)

f

ME555 : Computational Fluid Dynamics

19

I. Sezai – Eastern Mediterranean University

Surface areas Consider the cell of node P. Let the top surface be f1 with vertices v1,v2,v3,v4,v5.

A1 v4 v5

Calculation of area of surface f1: D Decompose the h surface f f1 into i triangles i l with one common vertex. For example choose common vertex to be v1.

v3 f1

v2

v1

The surface normal vector of the whole cell face f1 is the sum of the surface vectors of all triangles. A1 =

1 Nv ∑ ⎡( ri −1 − r1 ) × ( ri − r1 )⎤⎦ 2 i =3 ⎣

(11 12) (11.12)

Nv = number of vertices in the cell face (5 vertices for face f1) ri = position vector of vertex i Note that there are (Nv – 2) triangles. Also note that the numbering of the vertices is counter-clockwise when looking from outside the cell. ME555 : Computational Fluid Dynamics

20

I. Sezai – Eastern Mediterranean University

10

The cell face center can be found by averaging the coordinates of the vertices of the face. ( xc ) f =

Nv

1 Nv

∑x Nv

1 ( yc ) f = Nv

∑y

1 ( zc ) f = Nv

Nv

ME555 : Computational Fluid Dynamics

(11 13) (11.13)

i

i =1

∑z i =1

(11.14)

i

i =1

(11.15)

i

21

I. Sezai – Eastern Mediterranean University

Diffusion term The diffusion term is approximated as

∫ (Γ grad φ ) ⋅ dA ≈ ∑ A

f = nb ( P )

Γ f ∇φ f .A f

(11.17)

f2 N3

where the summation is over the faces of a cell. (over the faces f1, f2, f4, f4 for cell P in the figure). Define: PN = rN – rP (11.18) rP , rN: position vector of centroids of cells P and N ePN : unit vector along line PN n: unit normal vector of surface f Af: surface area vector

N2 N1 f1

P

f3

f4 N4

Af ePN P

et n

θ

f

PN

N

We wish to approximate the diffusive flux D f = Γ f ∇φ f ⋅ A f at the surface ME555 : Computational Fluid Dynamics

22

I. Sezai – Eastern Mediterranean University

11

Diffusion term The gradient expression should involve centroids P and N of the cells. e PN =

By definition

PN d PN

Af

ePN

(11.19)

d PN = PN

et n

θ N

f

P

The gradient in the direction of e can be written as ∂φ f φ − φP ∇φ f .e PN = = N ∂ePN

At

Ad

(11.20)

d PN

If the surface vector Af is written as the sum of two vectors Ad and At A f = A d + At (11.21) where Ad is in the direction of line PN,, then (11.22) ∇φ f .A f = ∇φ f .A d + ∇φ f .A t = A d ∇φ f .e PN + ∇φ f .A t φ − φP ∇φ f .A f = A d N + ∇φ f .A t (11.23) 1 424 3 d PN 14243 orthogonal term

non-orthogonal (cross diffusion) term

Cross diffusion term is treated as source term. ME555 : Computational Fluid Dynamics

23

I. Sezai – Eastern Mediterranean University

Diffusion term Several options can be used in defining the direction of At. Af ePN P

et n f

Af

At ePN

θ N Ad

(a) Minimum correction

P

et n f

Af

At

θ

ePN N

P

Ad

(b) Orthogonal correction

A d = Af cos θ e PN = ( A f .e PN )e PN A t = A f − A d = Af (n − cos θ e PN ) = A f − ( A f .e PN )e PN

(b) Set At such that |Ad| is equal to |Af| (orthogonal correction). A t = A f − A d = Af (n − e PN ) A d = Af e PN (c) Set At to be parallel to et (overrelaxed correction). Af ⋅Af Af Ad = e PN = e PN A f ⋅ e PN cos θ Af ⋅Af Af At = A f − A d = A f − e PN = A f − e PN A f ⋅ e PN cos θ 24

f

At

θ N

Ad

(c) overrelaxed correction

(a) Set At to be perpendicular to Ad (minimum correction).

ME555 : Computational Fluid Dynamics

et n

(11.24) (11.25) (11 26) (11.26) (11.27) (11.28)

I. Sezai – Eastern Mediterranean University

12

Diffusion term Substituting the terms in Eq.(11.28) for the overrelaxed approach into Eq. (11. 23)

φ − φP ∇φ f .A f = Ad N + ( ∇φ ) f .A t d PN 1424 3 14243 non-orthogonal orthogonal th l term t

where

Ad = A d =

(11.29) Linear interpolation p factor

( (cross diffusion) diff i ) term

Af ⋅Af

A t = A f − Ad e PN

A f ⋅ e PN

The overbar term in Eq. (11.29) is obtained via linear interpolation from the adjacent nodal values:

∇φ f = (1 − f P )∇φP + f P ∇φN fP =

Pf PN

or preferably: f P =

(11.30)

Pf ⋅ e PN Pf ⋅ e PN = 1/2 PN ( PN ⋅ PN )

A more accurate approximation is to use Eq. (11.98) Note that Ad is the magnitude of area vector Ad. Also note that Ad is in the direction of PN and At is perpendicular to Af. ME555 : Computational Fluid Dynamics

25

I. Sezai – Eastern Mediterranean University

Diffusion term Orthogonal term is treated implicitly, whereas the cross diffusion term is treated as a source term in a deferred correction approach. Overrelaxed approach (method (c)) is preferred, because the coefficient of the implicit term is greater than that of methods (a) and (b), and it increases with the skew angle θ. Thus, the diagonal dominance of the coefficient matrix is enhanced.

ME555 : Computational Fluid Dynamics

26

I. Sezai – Eastern Mediterranean University

13

Calculation of Gradients Calculation of Gradient at Nodal Points An expression is required to calculate gradient at points P and N in equation (11.29). One popular method is to use Gauss’ rule as: 1 1 φ dA = ∫v V ∫s VP where φf is calculated from ∇φ P =

1 V

∇φ dV =

nf

∑φ f =1

f

Af

f2 N3

φ f = φ f + ∇φ f ⋅fof o

(11.32)

o

where

φ f = (1 − f P )φP + f P ∇φN , o

(11 31) (11.31)

fP =

∇φ ffo = (1 − f P )∇φP + f P ∇φN

f3

P

N2 f1

rPf

N1

rNf

f4

Pf ⋅ e PN

N4

( PN ⋅ PN )

1/ 2

fof = Pf − Pfo = Pf − (Pf ⋅ e PN )e PN

f

P

N

ePN f o

An iterative process is required to calculate gradients: Step 1: Calculate gradient from Eq. (11.31) Step 2: Calculate φf from Eq. (11.32) Step 3: Repeat steps 1 and 2 until convergence. (4-5 repetitions required ) ME555 : Computational Fluid Dynamics

27

I. Sezai – Eastern Mediterranean University

Diffusion source term Note that in the general transport equation (11.4), only part of the diffusion term appears explicitly. The other part is buried in the source term: ∂ ( ρφ ) + div( ρφ v ) = div(Γ grad φ ) + sφ ∂t For example, in the x-momentum equation (2.32a) the diffusion term is ∂ ⎡ ∂u ⎤ ∂ ⎡ ⎛ ∂u ∂v ⎞ ⎤ ∂ ⎡ ⎛ ∂v ∂w ⎞ ⎤ + λ div u ⎥ + ⎢ μ ⎜ + ⎟ ⎥ + ⎢ μ ⎜ + 2μ ⎟⎥ ⎢ ∂x ⎣ ∂x ⎦ ∂y ⎣ ⎝ ∂y ∂x ⎠ ⎦ ∂z ⎣ ⎝ ∂z ∂x ⎠ ⎦ =

∂ ⎛ ∂u ⎞ ∂ ⎛ ∂u ⎞ ∂ ⎛ ∂u ⎞ ⎜μ ⎟+ ⎜μ ⎟+ ⎜μ ⎟ ∂x ⎝ ∂x ⎠ ∂y ⎝ ∂y ⎠ ∂z ⎝ ∂z ⎠

div( μ grad u )

⎡ ∂ ⎛ ∂u ⎞ ∂ ⎛ ∂v ⎞ ∂ ⎛ ∂w ⎞ ⎤ ∂ + ⎢ ⎜μ ⎟+ ⎜μ ⎟+ ⎜μ ⎟ ⎥ + (λ div u) ⎣ ∂x ⎝ ∂x ⎠ ∂y ⎝ ∂x ⎠ ∂z ⎝ ∂x ⎠ ⎦ ∂x = div( μ grad u ) + sMx

sMx

For incompressible fluids divu = 0. If further, μ = 0, then sMx = 0 (from discussion just after equation (2.34)) . ME555 : Computational Fluid Dynamics

28

I. Sezai – Eastern Mediterranean University

14

Diffusion source term However, sometimes μ may be a function of temperature. Morever, in the momentum equations of turbulence models, μ is replaced by (μ + μt), where μt is the turbulent viscosity, which is not constant at all. For those cases, sMx will not be zero. The source term sMx can be expressed in vector form as x sMx M = div ( μ gradu ⋅ i ) = div ( μ∇ u )

∂ ∂ ∂ where ∇ x is the x-component of the gradient operator ∇ = i + j + k ∂x ∂y ∂z and u = ui + vj + wk is the velocity vector. Integrating over the control volume gives:

Sudiff =

∫ div( μ∇ u)dV = ∫ (μ∇ u) ⋅ ndA ≈ ∑ x

x

CV

f = nb ( P )

S

( μ∇ x u ) f Afx + ( μ∇ x v) f Afy + ( μ∇ x w) f Afz

Similarly, the corresponding terms for the y- and z-momentum equations are

Svdiff =

∫ div(μ∇ u)dV = ∫ (μ∇ u) ⋅ ndA ≈ ∑

( μ∇ y u ) f Afx + ( μ∇ y v) f Afy + ( μ∇ y w) f Afz

∫ div( μ∇ u)dV = ∫ (μ∇ u) ⋅ ndA ≈ ∑

( μ∇ z u ) f A fx + ( μ∇ z v) f Afy + ( μ∇ z w) f Afz

y

y

CV

S

diff w

=

f = nb ( P )

S

z

z

f = nb ( P )

CV S (11.33) The gradients ∇u , ∇v and ∇w at the nodal points are found using Gauss rule (eqn.(11.31) and their face values are found via linear interpolation using eqn. (11.30).

ME555 : Computational Fluid Dynamics

29

I. Sezai – Eastern Mediterranean University

Convection term The convection term is approximated as n fP

n fP

∫ ( ρφ v) ⋅ dA ≈ ∑ ( ρ v ⋅ A) f φ f = ∑ m& f φ f A

f =1

f2

(11.33)

f =1

N3

nfP: number b off ffaces off cell ll P, P (faces (f f1, f2, f3, f4 for f cell ll P in the figure).

m& f = ( ρ v ⋅ A) f (Mass flow rate through face f)

N2 N1

f3

f1

P

f4 N4

(11.34)

To calculate mf , vf at the cell faces should be found. vf is calculated using the momentum interpolation method (MIM) of Rhie and Chow, to be discussed later.

Calculation of φf for Convection Terms a) First order methods 1) Upwind difference method (UD) φ f = φP for m& f > 0

(11.35)

φ f = φN for m& f < 0

(11.36)

ME555 : Computational Fluid Dynamics

30

I. Sezai – Eastern Mediterranean University

15

Convection term High Order Methods N f 1) Linear Upwind Difference Scheme (LUD), N nd f or 2 order Upwind Difference Scheme (SOU) N P f r r LUD formulation along a line in 1-D is: f ⎛ φ −φ ⎞ 1 φ f = φ Uf + ⎜ P U ⎟ Δx N (11.37) ⎝ Δx ⎠ 2 where subscript U refers to upwind node. This can be extended to 3D space using Taylor series expansion around point U: (11.38) φ f = φ Uf + ∇φU ⋅rU U where φ f = φP , ∇φU = ∇φP , rU = rPf = r f − rP for m& > 0 2

2

1

1

3

3

Pf

Nf

4

4

φ Uf = φN , ∇φU = ∇φN , rU = rNf = r f − rN

φ f = φP + ∇φP ⋅ rPf

or

for m& < 0

for m& > 0

(11.39)

φ f = φN + ∇φN ⋅ rNf for m& < 0 In Eq. (11.39) ∇φP and ∇φN are calculated using Eq. (11.31) ME555 : Computational Fluid Dynamics

31

I. Sezai – Eastern Mediterranean University

Method of deferred correction The convection term φf at the cell faces is implemented as the sum of the upwind value and other higher order terms which are evaluated at the previous iteration. HO

φ f = φ Uf + Δφ of = φ Uf + (φ of − φ oU f ) Δφ of is added to source term Su φ uf = face value of φ to be computed by 1st order upwind method φ of = face value of φ computed by high order scheme from previous old values φ of = face value φ computed by 1st order upwind method from previous old values HO u

For unstructured grids:

φ f = φP + ∇φP ⋅ rPf

for m& > 0

φ f = φN + ∇φN ⋅ rNf

for m& < 0 HO

U

= (φ of − φ of ) is added to source term

ME555 : Computational Fluid Dynamics

32

I. Sezai – Eastern Mediterranean University

16

Source term Source term in Eq. (11.8) is discretized as

∫ sφ dV = sφV

(11.40)

P

∂p CV x For x-momentum equation: sφ = − ∂x = −(∇p) P ∂p y For y-momentum equation: sφ = − = −(∇p) P ∂y ∂p ∂p ∇p = i + j = (∇p ) x i + (∇p ) y j where ∂x ∂y

(11.41)

Using Gauss’ rule (Eq. (11.31 with (11.32)): n

1 1 1 f ∇ pdV = pd A = ∑ pf A f V ∫v V ∫s VP f =1 nf nf ⎞ ⎞ 1 ⎛ 1 ⎛ ∇pP = ⎜ ∑ p f A xf ⎟ i + ⎜ ∑ p f A yf ⎟ j = (∇p) x i + (∇p) y j ⎜ ⎟ ⎜ ⎟ VP ⎝ f =1 ⎠ VP ⎝ f =1 ⎠n

∇p P =

f ∂p VP = −(∇p ) Px VP = −∑ p f A xf ∂x f =1 CV nf ∂p y p f A yf ∫CV sφ dV = sφVP = ∂y VP = −(∇p) P VP = −∑ f =1

∫ sφ dV = sφV

P

=−

ME555 : Computational Fluid Dynamics

33

(11.42) (11 43) (11.43)

for x-momentum equation (11.44)

for y-momentum equation I. Sezai – Eastern Mediterranean University

Substituting the discretized forms of the transient, convection and diffusion terms into the general equation (11.8) aP φ P = ∑ a N φ N + S c (11.45) N = nb ( P )

Sc = scbodyVP + Scbc + Sctrans + Scdc + Scpres + Scdiff , aN =

Γ f Ad

Sp = s

d PN

− min ⎡⎣ m& f , 0⎤⎦ ,

V +S +S

body p P

bc p

∑

f = nb ( P )

+

aP =

+S +S dc p

Sctrans = aPoφPo ,

S trans = −aPo , p Scdc =

trans p

pres p

,

aPo =

(

∑

( )

Γ f ∇φ

f

.At ,

aN − S p ,

(S pdc = S ppres = 0)

ρVP Δt

) ∑

min ⎣⎡ −m& f , 0⎦⎤ φ fHO − φP −

f = nb ( P )

∑

N = nb ( P )

(Scdiff is given by Eqn. (11.33)

f = nb ( P )

(

min ⎣⎡m& f , 0⎦⎤ φ fHO − φN

At = A f − Ad e PN ,

Ad =

)

(11.46)

Af ⋅Af A f ⋅ e PN

nf ⎧ x x for x-momentum equation ⎪= −(∇p) P VP = −∑ p f A f f =1 ⎪ Scpres = ⎨ nf ⎪= −(∇p) y V = − p A y for y-momentum equation ∑ P P f f ⎪ f =1 ⎩ , scbody = source terms per unit volume in the differential equation (body forces) s body p bc S bc p , Sc = source contributions of near-boundary cells

Generally, the source terms are linearized as: S = Sc + S pφ

ME555 : Computational Fluid Dynamics

34

and S = sV

I. Sezai – Eastern Mediterranean University

17

In solving steady flows, Eq. (11.45) is relaxed as aP

α

φP =

∑

N = nb ( P )

aN φN + S + (1 − α )

aP

α

φP*

(11.47)

0 ≤ α ≤ 1 underrelaxation factor factor. In that case, the terms in Eq. (11.46) are redefined as aP * ⎧ ⎪⎪ S ← S + (1 − α ) α φP ⎨ ⎪ a ← aP ⎪⎩ P α where superscript * refers to the previous iteration value.

ME555 : Computational Fluid Dynamics

35

(11.48)

I. Sezai – Eastern Mediterranean University

MIM method To find mf in Eq.(11.34), velocity vector vf at cell faces are calculated using the momentum interpolation method (MIM) originally proposed by Rhie and Chow. Rewriting the momentum equations for a 2D flow in the form of Eq. (11.45) u aP

α

a

v P

α

uP =

∑

aNu u N + bPu − V (∇p ) P .i

∑

aNv vN + bPv − V (∇p ) P .j

N = nb ( P )

vP =

(11.49) (11.50)

nb ( P ) Where, b is the source termsN =excluding pressure source. Above Eqns can be written as,

where

⎧u P ⎫ ⎧ H [u ]P ⎫ ⎧ D[u ]P ⎨ ⎬−⎨ ⎬ = −⎨ ⎩ vP ⎭ ⎩ H [ v ]P ⎭ ⎩ 0

H [φ ]P =

∑

aN φN + bP φ

N = nb ( P )

φ

0 ⎫ ⎧ (∇p ) P ⋅ i ⎫ ⎬⎨ ⎬ D[v]P ⎭ ⎩(∇p ) P ⋅ j⎭

D[φ ]P =

aPφ / α

V aφP / α

(11.51) (11.52)

Defining the vector forms of the above operators as ⎡ H [u ] P ⎤ H[ v ]P = ⎢ ⎥ ⎣ H [ v ]P ⎦

⎡ D[u ]P DP = ⎢ ⎣ 0

0 ⎤ D[v]P ⎥⎦

⎡ (∇P) P ⋅ i ⎤ ⎡(∇P ) Px ⎤ (∇P) P = ⎢ ⎥=⎢ y⎥ ⎣(∇P) P ⋅ j⎦ ⎣(∇P ) P ⎦

(11.53)

The momentum equation in vector form becomes

v P = H[ v ]P − D P (∇p ) P ME555 : Computational Fluid Dynamics

36

(11.54) I. Sezai – Eastern Mediterranean University

18

MIM method Similarly, for point N

v N = H[ v]N − D N (∇p ) N Eq. (11.54) can be written at the cell face f as v f = H f − D f (∇p) f

( (11.55) ) Assuming that the term Hf can be approximated by using an interpolation between points P and N, that is Hf = Hf where the overbar indicates linear average between point P and N. Then, Eq. (11.55) becomes (11.56) v f = H f − D f (∇p) f The linearly averaged form of Eq.(11.55) is (11.57)

H f = v f + D f ∇p f D f ∇p f = D f ∇p f

Using,

(11.58)

Eq (11.57) becomes, H f = v f + D f ∇p f ME555 : Computational Fluid Dynamics

37

I. Sezai – Eastern Mediterranean University

Mass flow rate where

v f = (1 − f P ) v P + f P v N D f = (1 − f P )D P + f P D N ∇p f = (1 − f P )(∇p) P + f P (∇p) N Df = Df

where, fP is the geometric interpolation factor defined in Eq. (11.30). A more accurate interpolated value can be found using Eqn. (11.32). Substituting eq (11.58) into eq (11.56)

(

v f = v f − D f (∇p ) f − D f ∇p f

Mass flow rate:

(11.59)

)

m& f = ( ρ v ⋅ A ) f

Substituting Eq. (11.59) for vf :

(

m& f = ρ f v f ⋅ A f − ρ f D f ∇p f ⋅ A f − ρ f D f (∇p ) f ⋅ A f ME555 : Computational Fluid Dynamics

38

)

(11.60)

I. Sezai – Eastern Mediterranean University

19

Mass flow rate The term ∇p f .A f is discretized using Eq. (11.29)

∇p f .A f = Ad

p N − pP + ∇p d PN

( )

f

.( A t ) f

(11.61)

where the overbar denotes linear interpolation between the two cells straddling the surface f. Substituting into Eq. (11.60)

⎛ p − pP m& f = ρ f v f ⋅ A f + ρ f D f (∇p ) f ⋅ A f − ρ f D f ⎜ N Ad + ∇p ⎝ d PN Af ⋅Af Ad or using = d PN A f ⋅ PN

( )

m& f = ρ f v f ⋅ A f + ρ f D f (∇p ) f ⋅ A f − ρ f

(D f A f ) ⋅ A f A f ⋅ PN

f

⎞ .( A t ) f ⎟ ⎠

(11.62)

( p N − pP ) − ρ f ( D f ∇p f ) ⋅ ( A t ) f (11.63)

= aNp′ from Eq. (11.77) ME555 : Computational Fluid Dynamics

39

I. Sezai – Eastern Mediterranean University

SIMPLE Algorithm Consider the momentum equations in Eq. (11.56) rewritten at the face f as

v f − H[ v ] f = − D f (∇p ) f

(11.64)

Let p', u', v' be the correction needed to correct the guessed pressure and velocity fields, fi ld i.e. i

p = p * + p′

and

v = v * + v′

(11.65)

For a guessed pressure field p* the corresponding velocity can be written as

v*f − H[ v* ] f = −D f (∇p* ) f

(11.66)

Subtracting Eq. (11.66) from Eq. (11.64) yields

v′f − H[ v′] f = −D f (∇p′) f ME555 : Computational Fluid Dynamics

40

(11.67) I. Sezai – Eastern Mediterranean University

20

SIMPLE Algorithm In SIMPLE method the second term on the left hand side of Eq. (11.67) is neglected. Then, Eq. (11.67) becomes

v′f = − D f ∇p′f

(11.68)

The discretized continuity equation is ( ρ P − ρ Po ) VP + ∑ m& f = 0 Δt f = nb ( P )

where

m& f = ρ v f ⋅ A f

(11.69)

( ρ P − ρ Po ) VP + ∑ m& *f + ∑ m& ′f = 0 Δt f = nb ( P ) f = nb ( P )

or

(11.70)

m& ′f = ( ρ v′) f ⋅ A f

where

Substituting v'f from Eq. (11.68) (11.71)

m& ′f = ( ρ v′) f ⋅ A f = − ρ f (D f ∇p′f ) ⋅ A f

ME555 : Computational Fluid Dynamics

41

I. Sezai – Eastern Mediterranean University

SIMPLE Algorithm The term ∇p′f .A f is discretized using Eq. (11.29) ∇p′f .A f =

A f ⋅ A f p′N − p′P + ∇p′ .( A t ) f f A f ⋅ e PN d PN

( )

(11.72)

Substituting Eq. (11.72) into Eq. (11.71) gives m& ′f = − ρ f

(D f A f ) ⋅ A f A f ⋅ d PN

( )

(11.73)

( p′N − p′P ) − ρ f D f ∇p′ .( A t ) f f

Substituting Eq. (11.73) into Eq. (11.70) , the continuity equation becomes

∑

f = nb ( P )

ρf

(D f A f ) ⋅ A f A f ⋅ d PN

( p′N − pP′ ) =

∑

f = nb ( P )

m& *f +

( ρ P − ρ Po ) VP − ∑ ρ f D f ∇p′ .( A t ) f f Δt f = nb ( P )

( )

(11.74)

Equation (11.74) is also called the pressure correction equation.

ME555 : Computational Fluid Dynamics

42

I. Sezai – Eastern Mediterranean University

21

SIMPLE Algorithm The pressure correction equation (11.74) can be written as: aPp′ pP′ = aNp′ = ρ f p′ P

a =

∑

N = nb ( P )

(11.77)

(D f A f ) ⋅ A f

∑

A f ⋅ PN

N = nb ( P )

bPp′ = −

aNp′ p′N + bPp′

aNp′

( ρ P − ρ Po ) VP − ∑ m& *f + ∑ ρ f D f ∇p′f ⋅ ( A t ) f Δt f = nb ( P ) f = nb ( P )

The last term on the right hand side of bP term can be neglected. ⎛ Du Df = ⎜ f ⎜ 0 ⎝

0 ⎞ ⎟ D vf ⎟⎠

ME555 : Computational Fluid Dynamics

(11.78)

43

I. Sezai – Eastern Mediterranean University

SIMPLE Algorithm From Eq. (11.73), the mass correction is

m& ′f = −aNp′ ( p′N − pP′ ) − ρ f D f ∇p′f ⋅ ( A t ) f

(11.79)

After solving the p' field from Eq. (11.77) the mass flow rate at the surfaces are corrected by using Eq. (11.79)

m& f = m& *f − aNp′ ( p′N − p′P ) − ρ f D f ∇p′f ⋅ ( A t ) f

(11.80)

The nodal velocities are corrected as v P = v*P + v′P = v*P − D P ∇p′P

(11.81)

and the pressure field is corrected by using

pP = pP* + α p pP′

(11.82)

αp = pressure under-relaxation factor ME555 : Computational Fluid Dynamics

44

I. Sezai – Eastern Mediterranean University

22

SIMPLE Algorithm Step 1: Solve the discretized transport equation (11.47) for φ = u and φ = v aP

a

φ = ∑ a φ + S + (1 − α ) P φP* α P N =nb ( P ) N N α Step 2: Calculate mass flow rate at cell faces (Eqn. (11.63))

m& f = ρ f v f ⋅ A f + ρ f D f (∇p ) f ⋅ A f − aNp ' ( pN − pP ) − ρ f ( D f ∇p f ) ⋅ ( A t ) f Step3: Solve pressure correction equation (11.77) aPp′ pP′ =

∑

N = nb ( P )

aNp′ p′N + bPp′

Step 4: Correct velocities and pressure at points P using Eqn’s (11.81), (11.82) pP = pP* + α p pP′ v P = v*P + v′P = v*P − D P ∇pP′ Step 5: Correct mass flow rate using equation (11.80) : m& f = m& *f − aNpp′ ( p′N − p′P ) − ρ f D f ∇p′f ⋅ ( A t ) f

Step 6: Solve the discretized transport equation (11.47) for other unknowns (i.e. for φ = Temperature) a a P φ = ∑ a φ + S + (1 − α ) P φP* α P N =nb ( P ) N N α Step 7: Repeat steps 1 to 6 until convergence. ME555 : Computational Fluid Dynamics

45

I. Sezai – Eastern Mediterranean University

Modifications to original MIM method In the original MIM method proposed by Rhie and Chow, using equations (11.5611.57) for the face values Hf and Df , makes the solution to depend on the relaxation parameter α as well as the time step size Δt. To overcome this, Eqns (11.49-11.50) are written as (Bo Yu et al. , Numer. Heat Transf., Part B, v42, pp. 141-166, 2002) :

vP = HP − where

α aP

VP (∇p ) P + (1 − α ) v nP−1 +

α aP

a oP v oP +

α aP

b d VP (sbody )P c

⎡⎛ ⎤ ⎞ dc H P = α a ⎢⎜ ∑ a N v N ⎟ + (Sbc c ) P + (S c ) P ⎥ ⎢⎣⎝ N = nb ( P ) ⎥⎦ ⎠

(a)

−1 P

Note that the source terms resulting from relaxation, unsteady term and body forces has been separated from H term. The correspondig equation at the face is α α α (b) v f = H f − V f (∇p ) f + (1 − α ) v nf −1 + a of v of + V f (sbody )f c af af af or

v f = H f − D f (∇p ) f + (1 − α ) v nf −1 + ME555 : Computational Fluid Dynamics

46

α af

a of v of + D f (sbody )f c

(c)

I. Sezai – Eastern Mediterranean University

23

Modifications to original MIM method Following the procedure as in Rhie and Chow method, the face velocity can be expressed as v f = v f − ⎡⎣ D f (∇p ) f − D f (∇p ) f ⎤⎦ + (1 − α ) ⎡ v nf −1 − v nf −1 ⎤ ⎣ ⎦

(d)

+ α (a f ) −1 ⎡a of v of − a of v of ⎤ ⎣ ⎦ ) f − D f (sbody )f ⎤ + ⎡ D f (sbody c c ⎣ ⎦

where af is calculated using

1 1 1 = (1 − f P ) + f P af aP aN It should be noted that Eqn. (d) is equivalent to Eqn. (c). However, Eqn (d) may be preferred since extra storage is required for Hf in Eqn (c). ME555 : Computational Fluid Dynamics

47

I. Sezai – Eastern Mediterranean University

Modifications to original MIM method The source term (sc)f is estimated using

(sbody ) f = (1 − f P )(sbody ) P + f P (s body )N c c c

(f)

However, if the source term is high, the following equation can be used to estimate ((sc)f: 1 1 (sbody ) f = ((sbody ) P + (s body ) N ) + (∇ (s body ) P ⋅ rPf + ∇(sbody ) N ⋅ rNf ) (g) c c c c c 2 2 where the gradient at the nodal point P can be calculated using the Gauss' rule. The use of the neighboring cell source values for the estimation of (sc)f, is expected to circumvent the convergence difficulties encountered in flows with large source terms.

ME555 : Computational Fluid Dynamics

48

I. Sezai – Eastern Mediterranean University

24

TVD Schemes in Unstructured Grids For Cartesian grids, for ue > 0, using a TVD scheme, φe may be written as

1 2 where φP − φW r= φE − φP

φe = φP + ψ (r )(φE − φP )

φe

φP

(11.83)

φE

φW e

w WW

W

P

E

EE

ve

ψ = flux limiter function ψ = 0 for UD scheme ψ = 1 for CD scheme ψ = r for LUD scheme ψ = 0.25r + 0.75 for QUICK scheme However, in unstructured grids, for face f1 of cell P, the upwind point U is not available. There are methods to construct U, but these are costly. In the absence of φU, the method proposed by Darwish and Moukalled (2003) will be used ME555 : Computational Fluid Dynamics

49

f2

N2 N1

N3

f3

P

f1

U f4 N4

I. Sezai – Eastern Mediterranean University

TVD Schemes in Unstructured Grids Since the index based notation used above is not suitable for unstructured grids, a more appropriate notation proposed by Darwish and Moukalled (2003) will be adopted.

N2

f2 N3

f1

P

D

U

Locate U such that rUD = 2rCD

f4

f

(P) N4

Using this notation:

C: upwind node around face f D: downwind node around face f U: (virtual) upwind node of C

(N) v

C N1

f3

f

(P)

(N)

U

C

v

D

1 2

φ f = φC + ψ (rf )(φD − φC ) where

rf =

φC − φU φD − φC

(11.84) (11.85)

The main difficulty in using TVD schemes in unstructured grids lies in the need for defining a virtual node U. ME555 : Computational Fluid Dynamics

50

I. Sezai – Eastern Mediterranean University

25

Exact r Formulation Returning to the definition of r:

rf =

φC − φU φD + (φC − φU ) − φD (φD − φU ) − (φD − φC ) = = φD − φC φD − φC φD − φC

(11.86)

Note that D and C are the nodes straddlingg the interface and thus are readilyy available. For node U we can write:

(φD − φU ) = ∇φC ⋅ rUD = (2∇φC ⋅ rCD )

(11.87)

since rUD = 2rCD (C is at the center of the UD segment). Other positions of U could also be chosen but with a loss of accuracy as the nodal gradient yields a second order accuracy only when the difference is centered at C. Th fformulation The l ti off r becomes b

rf =

(2∇φC ⋅ rCD ) − (φD − φC ) (2∇φC ⋅ rCD ) = −1 φD − φC φD − φC

(11.88)

After calculating rf from Eq. (11.88) find ψ(rf) from a TVD scheme (chapter 5) and then calculate φf from Eq. (11.84). ME555 : Computational Fluid Dynamics

51

I. Sezai – Eastern Mediterranean University

Exact r Formulation Note that this r-formulation can also be used for any high order scheme without TVD. For example for QUICK scheme: ψ(rf) = 0.25rf + 0.75 Then, the application of high order convection schemes with TVD in a deferred correction manner, is added as a source term to Sdc in accordance with Eq. (11.46) as S dc =

∑

f = nb ( P )

(

) ∑

min ⎡⎣ − m& f , 0 ⎤⎦ φ fHO − φP −

f = nb ( P )

(

min ⎡⎣ m& f , 0 ⎤⎦ φ fHO − φN

)

(11.89)

where, φfHO is the φf value computed from Eq. (11.84).

ME555 : Computational Fluid Dynamics

52

I. Sezai – Eastern Mediterranean University

26

High order Convection Schemes in Unstructured Grids For structured grids the following notation is suitable for the one dimensional upwind type convection schemes:

f UU

U

C

D

DD

vf

f DD

D

C

UU

vf

C: upwind node around face f D: downwind node around face f U: upwind node of C ME555 : Computational Fluid Dynamics

U

53

I. Sezai – Eastern Mediterranean University

Convection Schemes With the above notation, the basic convection schemes are given in the table below. Table 11.1 Convection Schemes

Convection Scheme

Formula

φ f = φC φC + φD

UPWIND

φf =

CENTRAL

2

3 2

1 2

SOU (LUD)

φ f = φC − φU

QUICK

φ f = φC + φU − φD

3 4

DOWNWIND ME555 : Computational Fluid Dynamics

54

3 8 φ f = φD

1 8

I. Sezai – Eastern Mediterranean University

27

Virtual Node Formulation for High Order Schemes For unstructured grids, there is no line structure. For convection across a face point U does not exist. If a virtual node U is defined, then the above convection schemes can be directly used for unstructured grids.

f2 N3

N2 N1

f3

f1

P

f

(P)

v

C

C: upwind node around face f D: downwind node around face f U: (virtual) upwind node of C

((N)) D

U

f4

f

(P) N4

(N)

U

C

v

D

If point U is located such that rUD = 2rCD, then from Eq. (11.87)

φU = φD − 2∇φC ⋅rCD

(11.90)

Since φD, φC and φU are available and ∇φC is calculated from Eq. (11.31), then φf can be calculated from the convection formulas given in the Table 11.1. ME555 : Computational Fluid Dynamics

55

I. Sezai – Eastern Mediterranean University

Gradient Formulation for High Order Schemes Another method of applying high order convection schemes in unstructured grids is by means of formulating the convection schemes in terms of gradients. Consider QUICK scheme

1 8 1 = φC + ((φD − φU ) + 2(φD − φC )) 8 φ −φ ⎞ ⎛ φ −φ = φC + ⎜ D U + D C ⎟ 4 ⎠ ⎝ 8

φ f = φC + (3φD − 2φC − φU )

φ −φ 1 ⎛ φ −φ = φC + ⎜ D U + D C 2 ⎜ 4 rCf 2 rCf ⎝ 1 = φC + ( ∇φC + ∇φ f ) ⋅ rCf 2 ME555 : Computational Fluid Dynamics

56

⎞ ⎟ ⋅ rCf ⎟ ⎠ (11.91) I. Sezai – Eastern Mediterranean University

28

Gradient Formulation for High Order Schemes Consider SOU (LUD) scheme 3 1 φ f = φC − φU 2 2 φC − φU = φC + 2 φ −φ φ −φ = φC + D U − D C 2 2 = φC + ( 2∇φC − ∇φ f ) ⋅rCf Consider CD scheme φ +φ φf = C D 2 φ −φ = φC + D C 2 = φC + ∇φ f ⋅rCf ME555 : Computational Fluid Dynamics

(11.92)

(11.93)

57

I. Sezai – Eastern Mediterranean University

Gradient Formulation for High Order Schemes Table 11.2 Gradient Formulation of High Order Schemes

Convection Scheme

Formula

UPWIND

φ f = φC φ f = φC + ∇φ f ⋅ rCf

CENTRAL SOU (LUD) QUICK

φ f = φC + ( 2∇φC − ∇φ f ) ⋅ rCf φ f = φC +

DOWNWIND

1 ( ∇φC + ∇φ f ) ⋅ rCf 2 φ f = φD

Values of φC and ∇φC are available. We need an expression for ∇φ f . ME555 : Computational Fluid Dynamics

58

I. Sezai – Eastern Mediterranean University

29

Gradient at a Cell Face Gradient vector of a variable φ at a cell face can be decomposed into two components along directions d and t, where d is along PN and t is tangential (perpendicular) to d.

∇φ f = ∇ d φ f + ∇ t φ f ∇dφ f =

φ N − φP rPN

(11.94) ∇φ f

(11 95) (11.95)

e PN

∇ φ f = (∇φ f ⋅e PN )e PN d

ePN P

(11.96)

f

∇ dφ f

∇tφ f = ∇φ f − ∇ d φ f = ∇φ f − (∇φ f ⋅e PN )e PN

∇ tφ f

N

(11.97)

Substituting Eq. (11.95) and (11.97) into Eq. (11.94),

∇φ f =

φ N − φP rPN

e PN + ∇φ f − (∇φ f ⋅ e PN )e PN

(11.98)

∇φ f is obtained via linear interpolation from adjacent nodal values using Eq. (11.30)

Values of φC , ∇φC and ∇φ f are available. Calculate φ f using the expressions listed in Table (11.2). ME555 : Computational Fluid Dynamics

59

I. Sezai – Eastern Mediterranean University

30