Normalized maximum values for the * u* and

*velocity components.*w

Open access peer-reviewed chapter

Submitted: October 23rd, 2015 Reviewed: March 14th, 2016 Published: August 24th, 2016

DOI: 10.5772/63077

The present work addresses the numerical simulation of fluid flow for 2D problems. The physical principles and numerical models implemented in the software package EasyCFD are presented in a synthetic and clear way. The 2D form of the Navier-Stokes equations is considered, using the eddy-viscosity concept to take into account turbulence effects upon the mean flow field. The k-ε and the k-ω Shear Stress Transport (SST) turbulence models allow for the calculation of the turbulent viscosity. The numerical model is based on a control volume approach, using the SIMPLEC algorithm on an unstructured quadrilateral mesh. The mesh arrangement is a non-staggered type. The coordinate transformation, integration discretization and solution method for the governing equations are fully described. As an example of application, the airflow around a NACA 0012 airfoil is calculated and the results for the aerodynamic coefficients are compared with available experimental data.

- CFD
- SIMPLEC
- unstructured mesh
- fluid flow
- 2D simulation

The software EasyCFD [1] is a 2D simulation tool aimed at an initiation in the field of computational fluid dynamics. The main guiding lines on its development were the simplicity and the intuitiveness of utilization, in a self-contained package. The physical domain is discretized with quadrilateral unstructured meshes, allowing the simulation to deal with complex geometries of any configuration virtually. The present work synthesizes the main physical principles and numerical models implemented for the solution of 2D fluid flow problems, including heat transfer, in arbitrarily shaped geometries.

Advertisement
## 2. Basic transport equations

∂ ( ρ u ) ∂ t + ∂ ∂ x ( ρ u 2 ) + ∂ ∂ z ( ρ u w ) = ∂ ∂ x [ Γ ( 2 ∂ u ∂ x − 2 3 d i v V → ) ] + ∂ ∂ z [ Γ ( ∂ u ∂ z + ∂ w ∂ x ) ] − ∂ p ∂ x ![]()

E1
∂ ( ρ u ) ∂ t + ∂ ∂ x ( ρ u 2 ) + ∂ ∂ z ( ρ u w ) = ∂ ∂ x [ Γ ( 2 ∂ u ∂ x − 2 3 d i ν V → ) ] + ∂ ∂ z [ Γ ( ∂ u ∂ z + ∂ w ∂ x ) ] − ∂ p ∂ x ![]()

E1a
∂ ρ ∂ t + ∂ ∂ x ( ρ u ) + ∂ ∂ z ( ρ w ) = 0 ![]()

E2
∂ ( ρ c p T ) ∂ t + ∂ ∂ x ( ρ c p u T ) + ∂ ∂ z ( ρ c p w T ) = ∂ ∂ x [ Γ ∂ T ∂ x ] + ∂ ∂ z [ Γ ∂ T ∂ z ] + S T ![]()

E3
∂ ( ρ c p T ) ∂ t = ∂ ∂ x ( λ ∂ T ∂ x ) + ∂ ∂ z ( λ ∂ T ∂ z ) + S T ![]()

E4
∂ ( ρ ϕ c 2 ) ∂ t + ∂ ∂ x ( ρ u ϕ c 2 ) + ∂ ∂ z ( ρ w ϕ c 2 ) = ∂ ∂ x ( Γ ∂ ϕ 2 ∂ x ) + ∂ ∂ z ( Γ ∂ ϕ c 2 ∂ z ) ![]()

E5
Γ = ρ D ϕ c + μ t S c t ![]()

E6
X = ϕ c 1 X 1 + ϕ c 2 X 2 ![]()

E7
ρ = ( ϕ c 1 ρ 1 + ϕ c 2 ρ 2 ) − 1 ![]()

E8

For the description of the transport equations, the * x* and

The Navier-Stokes equations describe momentum conservation and, for a 2D situation, may be stated as follows:

where * p* [N/m

In turn, the conservation of mass law, or continuity equation, is

The energy conservation equation is obtained considering the transport equation for the enthalpy * T* [K] is its temperature. For a fluid medium, the corresponding equation is

with the source term _{T} representing the heat generation rate per unit volume [W/m^{3}] and _{t} as the laminar and the turbulent Prandtl number. For solid regions, heat transfer is governed by conduction:

where

where * ϕ* is the mass fraction, or concentration, of the secondary fluid. The diffusion coefficient is given by

where * D*[m

where the constraint * ρ* is given by

Advertisement
## 3. Turbulence modelling

### 3.1. The k-ε turbulence model

μ t = C μ ρ k 2 ε ![]()

E9
∂ ( ρ k ) ∂ t + ∂ ∂ x ( ρ u k ) + ∂ ∂ z ( ρ w k ) = ∂ ∂ x [ ( μ + μ t σ k ) ∂ k ∂ x ] + ∂ ∂ z [ ( μ + μ t σ k ) ∂ k ∂ z ] + P k − ρ ε + G T ![]()

E10
∂ ( ρ ε ) ∂ t + ∂ ∂ x ( ρ u ε ) + ∂ ∂ z ( ρ w ε ) = ∂ ∂ x [ ( μ + μ t σ ε ) ∂ ε ∂ x ] + ∂ ∂ z [ ( μ + μ t σ ε ) ∂ ε ∂ z ] + ε k ( C 1 P k − C 2 ρ ε + C 3 G T ) ![]()

E11
P k = μ t [ 2 ( ∂ u ∂ x ) 2 + 2 ( ∂ w ∂ z ) 2 + ( ∂ u ∂ z + ∂ w ∂ x ) 2 ] ![]()

E12
G T = − β g μ t pr t ∂ T ∂ z ![]()

E13
C μ = 0 . 0 9 ; σ k = 1 . 0 ; σ ε = 1 . 3 ; C 1 = 1 . 4 4 ; C 2 = 1 . 9 2 ; C 3 = 1 . 4 4 ![]()

E14
y + = ρ u τ y μ = ρ C μ 0 . 2 5 k y μ ![]()

E15
y + ≤ 1 1 . 6 3 ⇒ τ 0 = μ v 0 − v y ![]()

E16
y + > 1 1 . 6 3 ⇒ τ 0 = μ ρ C μ 0 . 2 5 χ k ln ( E * y + ) ( v 0 − v ) ![]()

E16a
P k = μ t ( ∂ v ∂ y ) 2 ![]()

E17
y + ≤ 1 1 . 6 3 ⇒ ρ ε = ρ C μ 0 . 7 5 k 1 . 5 y + y ![]()

E18
y + > 1 1 . 6 3 ⇒ ρ ε = ρ C μ 0 . 7 5 k 1 . 5 ln ( E * y + ) χ y ![]()

E18a
ε 0 = C μ 0 . 7 5 k 1 . 5 x y ![]()

E19
y + ≤ 1 2 ⇒ q ˙ = μ c p pr y ( T 0 − T ) ![]()

E20
y + > 1 2 ⇒ q ˙ = ρ c p C l μ 0 . 2 5 k pr t [ 1 χ ln ( E * y + ) + p * ] ( T 0 − T ) ![]()

E20a
P * = 9 . 2 4 [ ( pr pr t ) 0 . 7 5 − 1 ] [ 1 + 0.28 exp ( − 0 . 0 0 7 pr pr t ) ] ![]()

E21
### 3.2. The k-ω SST turbulence model

∂ ( ρ k ) ∂ t + ∂ ( ρ u k ) ∂ x + ∂ ( ρ w k ) ∂ z = P k ¯ − β * ρ ω k + ∂ ∂ x ( ( μ + σ k μ t ∂ k ∂ x ) ) + ∂ ∂ z ( ( μ + σ k μ t ∂ k ∂ z ) ) ![]()

E22
∂ ( ρ ω ) ∂ t + ∂ ( ρ u ω ) ∂ x + ∂ ( ρ w ω ) ∂ z = ∂ ∂ x ( ( μ + σ ω μ t ∂ ω ∂ x ) ) + ∂ ∂ z ( ( μ + σ ω μ t ∂ ω ∂ z ) ) + ![]()

E23
+ α P k ¯ v t − β ρ g ω 2 + 2 ( 1 − F 1 ) ρ σ ω 2 1 ω ( ∂ k ∂ x ∂ ω ∂ x + ∂ k ∂ z ∂ ω ∂ z ) ![]()

E23
P k ¯ = min ( P k , 1 0 β * ρ k ω ) ![]()

E24
F 1 = tanh { { min [ max ( k β * ω y ; 5 0 0 v y 2 ω ) ; 4 ρ σ ω 2 k C D k ω y 2 ] } } ![]()

E25
C D k ω = max ( 2 ρ σ ω 2 1 ω ∂ k ∂ x j ∂ ω ∂ x j , 1 0 − 1 0 ) ![]()

E26
v t = a 1 k max ( a 1 ω ; S F 2 ) ![]()

E27
S = S i j S i j ; S i j = 1 2 ( ∂ u i ∂ x j + ∂ u j ∂ x i ) ![]()

E28
F 2 = tanh { [ max ( 2 k β * ω y ; 5 0 0 v y 2 ω ) ] 2 } ![]()

E29
α = F 1 α 1 + ( 1 − F 1 ) α 2 ![]()

E30
ω v i s = 6 v 0 . 0 7 5 y 2 ; ω log = μ τ 0 . 3 χ y ![]()

E31
ω 1 = ω v i s 2 + ω log 2 ![]()

E32
u τ v i s = U 1 y + ; u τ log = x U 1 ln ( E y + ) ![]()

E33
u τ = ( u τ v i s ) 4 + ( u τ log ) 4 4 ![]()

E34

Two of the most popular turbulence models are presented next.

The standard formulation of this turbulence model is described in, e.g., [2]. The turbulent viscosity is given by

The turbulence kinetic energy, * k* [m

The term _{k} represents the production rate of * k* as the results of the velocity gradients:

while the term _{T} accounts for the production or destruction of * k* and

with * g* [m/s

In the proximity of a wall, the previous equations should be modified to account for the viscous effects that become predominant. Wall functions ensure the connection between the viscous sub-layer and the inertia layer, at a location established by the ^{+} value:

where * u* represents the friction velocity and

where _{0} is the wall velocity, ^{*}=* 9.793* for smooth walls and

where, once again, * v* is the generic velocity component parallel to the wall and

For the turbulence kinetic energy, a zero flux along the direction perpendicular to the wall is assigned. For the dissipation rate, Eq. (11) is not employed in the node adjacent to the wall. Instead, the dissipation rate is given by

As for momentum, the energy flux is computed differently depending on the ^{+} values. In this case, we have

where _{0} is the wall temperature and ^{*} is the so-called Jayatillaka function:

The * k-ω* SST model [3] represents a combination of the

where * ω* is the frequency of dissipation of turbulent kinetic energy [s

The weighting function _{1} is given by

and

where * y* represents the distance to the neighbour wall and

where * S* is the invariant measure of the strain rate given by

and

The constants are computed as a blend of the * k-ε* and

The constants are _{1}= 5/9; _{1}= 3/40; _{k1} = 0.85; _{ω1} = 0.5; _{2}= 0.44; _{2} = 0.0828; _{k2} = 1; _{ω2} = 0.856; * β** = 0.09.

The near wall treatment for momentum and turbulence equations follows the proposal described in [5]. The basic principle behind the automatic wall functions is to switch from a low-Reynolds number formulation to a wall function based on the grid nodes proximity to the wall. According to these authors, the automatic wall treatment avoids the deterioration of results typical of low-Reynolds models when applied on too coarse meshes.

The known solutions for * ω* in the viscous (linear) and in the logarithmic near wall region are

The imposed value for * ω* at the first node close to a wall is

For the turbulence kinetic energy, a zero flux along the direction perpendicular to the wall is assigned. In turn, for the momentum equations, a similar reasoning applies, with expressions for the shear velocity in the viscous and in the logarithmic region:

with _{1} being the fluid velocity relative to the wall velocity. The wall shear stress is computed as follows:

Advertisement
## 4. Numerical method

### 4.1. Transformation of coordinates

∂ ϕ ∂ x = ∂ ϕ ∂ ξ ∂ ξ ∂ x + ∂ ϕ ∂ ζ ∂ ζ ∂ x = ξ x ∂ ϕ ∂ ξ + ζ x ∂ ϕ ∂ ζ ![]()

E35
J = | x ξ x ζ z ξ z ζ | = x ξ z ζ − x ζ z ξ ![]()

E36
ξ x = z ζ / J ; ξ z = − x ζ / J ; ζ x = − z ξ / J ; ζ z = x ξ / J ![]()

E37
∂ ∂ ξ ( J ξ x ) + ∂ ∂ ζ ( J ζ x ) = 0 ; ∂ ∂ ξ ( J ξ z ) + ∂ ∂ ζ ( J ζ z ) = 0 ![]()

E38
J ∂ ( ρ ϕ ) ∂ t + ∂ ∂ ξ ( J ρ U ϕ ) + ∂ ∂ ζ ( J ρ W ϕ ) = ∂ ∂ ξ [ Γ J g 1 1 ∂ ϕ ∂ ξ ] + ∂ ∂ ζ [ Γ J g 3 3 ∂ ϕ ∂ ζ ] + ∂ ∂ ξ [ Γ J g 1 3 ∂ ϕ ∂ ξ ] + ∂ ∂ ζ [ Γ J g 1 3 ∂ ϕ ∂ ζ ] + J S ϕ ![]()

E39
g 1 1 = ξ x 2 + ξ z 2 ; g 3 3 = ζ x 2 + ζ z 2 ; g 1 3 = ξ x ζ x + ξ z ζ z ![]()

E40
F ξ = J ρ U = ρ ( z ζ u − x ζ w ) ![]()

E41
F ζ = J ρ W = ρ ( x ξ w − z ξ u ) ![]()

E42
J ∂ ( ρ ϕ ) ∂ t + ∂ ∂ ξ ( F ξ ϕ ) + ∂ ∂ ζ ( F ζ ϕ ) = ∂ ∂ ξ [ Γ J g 1 1 ∂ ϕ ∂ ξ ] + ∂ ∂ ζ [ Γ J g 3 3 ∂ ϕ ∂ ζ ] + J ( S c r o s s + S ϕ ) ![]()

E43
J ∂ ( ρ u ) ∂ t + ∂ ∂ ξ ( F ξ u ) + ∂ ∂ ζ ( F ζ u ) = ∂ ∂ ξ [ Γ J g 1 1 ∂ u ∂ ξ ] + ∂ ∂ ζ [ Γ J g 3 3 ∂ u ∂ ζ ] − z ζ ∂ p ∂ ξ + z ξ ∂ p ∂ ζ + S c r o s s u ![]()

E44
J ∂ ( ρ w ) ∂ t + ∂ ∂ ξ ( F ξ w ) + ∂ ∂ ζ ( F ζ w ) = ∂ ∂ ξ [ Γ J g 1 1 ∂ w ∂ ξ ] + ∂ ∂ ζ [ Γ J g 3 3 ∂ w ∂ ζ ] + x ζ ∂ p ∂ ξ − x ξ ∂ p ∂ ζ + S c r o s s w + I ![]()

E45
J ∂ ρ ∂ t + ∂ ∂ ξ [ ρ ( z ζ u − x ζ w ) ] + ∂ ∂ ζ [ ρ ( x ξ w − z ξ u ) ] = J ∂ p ∂ t + ∂ ∂ ξ ( J ρ U ) + ∂ ∂ ζ ( J ρ W ) = 0 ![]()

E46
### 4.2. Integration

J ∂ ( ρ ϕ ) ∂ t + ∂ ∂ ξ ( F ξ ϕ − Γ J g 1 1 ∂ ϕ ∂ ξ ) + ∂ ∂ ζ ( F ζ ϕ − Γ J g 1 1 ∂ ϕ ∂ ζ ) = J ( S c r o s s + S ϕ ) ![]()

E47
J ∂ ( ρ ϕ ) ∂ t + [ F e ϕ e − D e ( ϕ E − ϕ P ) ] − [ F o ϕ o − D o ( ϕ P − ϕ O ) ] + [ F t ϕ t − D t ( ϕ T − ϕ P ) ] − [ F b ϕ b − D b ( ϕ P − ϕ B ) ] = J ( S 1 P + S 2 P ϕ P ) ![]()

E48
F e = F ξ e = [ ρ ( z ζ u − x ζ w ) ] e F o = F ξ o = [ ρ ( z ζ u − x ζ w ) ] o ![]()

E49
F t = F ζ t = [ ρ ( x ξ w − z ζ u ) ] t F b = F ζ b = [ ρ ( x ξ w − z ζ u ) ] b ![]()

E49
D e = ( Γ J g 1 1 ) e ; D o = ( Γ J g 1 1 ) o ; D t = ( Γ J g 3 3 ) t ; D b = ( Γ J g 3 3 ) b ![]()

E50
a P ϕ P = a E ϕ E + a O ϕ O + a T ϕ T + a B ϕ B + b ![]()

E51
a P ϕ P = ∑ n b a n b u n b + b ![]()

E51a
a
E
=
D
e
+ 〚
−
F
e
, 0 〛 ;
a
O
=
D
o
+ 〚
F
o
, 0 〛
a
T
=
D
t
+ 〚
−
F
t
, 0 〛 ;
a
B
=
D
b
+ 〚
F
b
, 0 〛
![]()

E52
b q u i c k = 1 8 〚 F e , 0 〛 ( − ϕ O − 2 ϕ P + 3 ϕ E ) − 1 8 〚 − F e , 0 〛 ( 3 ϕ P − 2 ϕ E − 3 ϕ E E ) + 1 8 〚 F o , 0 〛 ( 3 ϕ O − 2 ϕ P − ϕ E ) − 1 8 〚 − F o , 0 〛 ( − ϕ O O − 2 ϕ O + 3 ϕ P ) + 1 8 〚 F t , 0 〛 ( − ϕ B − 2 ϕ P + 3 ϕ T ) − 1 8 〚 − F t , 0 〛 ( 3 ϕ P − 2 ϕ T − 3 ϕ T T ) + 1 8 〚 F b , 0 〛 ( 3 ϕ B − 2 ϕ P − ϕ T ) − 1 8 〚 − F b , 0 〛 ( − ϕ B B − 2 ϕ B + 3 ϕ P ) + ![]()

E53
### 4.3. Pressure-velocity coupling

a P u P = ∑ n b a n b u n b − z ζ ∂ p ∂ ξ + z ξ ∂ p ∂ ζ + S c r o s s u ![]()

E54
a P ( 1 + 1 E ) u P ∗ = ∑ n b a n b u n b ∗ + a P E u P m − z ζ ∂ p ∗ ∂ ξ + z ξ ∂ p ∗ ∂ ζ + S c r o s s u ![]()

E55
a P ( 1 + 1 E ) u P = ∑ n b a n b u n b + a P E u P m − z ζ ∂ p ∂ ξ + z ξ ∂ p ∂ ζ + b ![]()

E56
p = p ∗ + p ' ; u P = u P ∗ + u ′ P ![]()

E57
a P ( 1 + 1 E ) u P ' = ∑ n b a n b u n b ' − z ζ ∂ p ' ∂ ξ + z ξ ∂ p ' ∂ ζ ![]()

E58
[ a P ( 1 + 1 E ) − ∑ n b a n b ] u P ' = ∑ n b a n b ( u n b ' − u P ' ) − z ζ ∂ p ' ∂ ξ + z ξ ∂ p ' ∂ ζ ![]()

E59
a ˜ P u P ' = − z ζ ∂ p ' ∂ ξ + z ξ ∂ p ' ∂ ζ ![]()

E60
a ˜ P = a P ( 1 + 1 E ) − ∑ n b a n b ![]()

E61
u P ' = − z ζ a ˜ P ∂ p ' ∂ ξ + z ξ a ˜ P ∂ p ' ∂ ζ ; w P ' = − x ξ a ˜ P ∂ p ' ∂ ζ + x ζ a ˜ P ∂ p ' ∂ ξ ![]()

E62
u P = u P ∗ − z ζ a ˜ P ∂ p ' ∂ ξ + z ξ a ˜ P ∂ p ' ∂ ζ ; w P = w P ∗ + x ζ a ˜ P ∂ p ' ∂ ξ − x ξ a ˜ P ∂ p ' ∂ ζ ![]()

E63
### 4.4. The pressure correction equation

J ρ P − ρ P 0 Δ t + [ ρ ( z ζ u − x ζ w ) ] e − [ ρ ( z ζ u − x ζ w ) ] o + [ ρ ( x ξ w − z ξ u ) ] t − [ ρ ( x ξ w − z ξ u ) ] b = 0 ![]()

E64
J ρ P − ρ P 0 Δ t + F e * − F o * + F t * − F b * − ( ρ g 3 3 a ˜ P ∂ p ' ∂ ξ ) e + ( ρ g 3 3 a ˜ P ∂ p ' ∂ ξ ) o − ( ρ g 1 1 a ˜ P ∂ p ' ∂ ζ ) t + ( ρ g 1 1 a ˜ P ∂ p ' ∂ ζ ) b + ( ρ g 1 3 a ˜ P ∂ p ' ∂ ζ ) e − ( ρ g 1 3 a ˜ P ∂ p ' ∂ ζ ) o + ( ρ g 1 3 a ˜ P ∂ p ' ∂ ξ ) t − ( ρ g 1 3 a ˜ P ∂ p ' ∂ ξ ) b = 0 ![]()

E65
g 1 1 = x ξ 2 + z ξ 2 = g 3 3 J 2 ; g 3 3 = x ζ 2 + z ζ 2 = g 1 1 J 2 ; g 1 3 = x ξ x ζ + z ξ z ζ = − g 1 3 J 2 ![]()

E66
( ∂ p ' ∂ ξ ) e = p E ' − p P ' ; ( ∂ p ' ∂ ξ ) o = p P ' − p O ' ; ( ∂ p ' ∂ ζ ) t = p T ' − p P ' ; ( ∂ p ' ∂ ζ ) b = p P ' − p B ' ![]()

E67
( ∂ p ' ∂ ξ ) t = 0 . 2 5 ( p E ' + p T E ' − p O ' − p T O ' ) ; ( ∂ p ' ∂ ξ ) b = 0 . 2 5 ( p E ' + p B E ' − p O ' − p B O ' ) ( ∂ p ' ∂ ζ ) e = 0 . 2 5 ( p T ' + p T E ' − p B ' − p B E ' ) ; ( ∂ p ' ∂ ζ ) o = 0 . 2 5 ( p T ' + p T O ' − p B ' − p B O ' ) ![]()

E68
a P p P ' = a E p E ' + a O p O ' + a T p T ' + a B p B ' + b ![]()

E69
a E = ( g 3 3 ⟨ a ˜ P / ρ ⟩ ) e ; a O = ( g 3 3 ⟨ a ˜ P / ρ ⟩ ) o ; a T = ( g 1 1 ⟨ a ˜ P / ρ ⟩ ) t ; a B = ( g 1 1 ⟨ a ˜ P / ρ ⟩ ) b ![]()

E70
a P = a E + a O + a T + a B ![]()

E71
b = J ρ P 0 − ρ P Δ t − F e * + F o * − F t * + F b * − ( g 1 3 ⟨ a ˜ P / ρ ⟩ ∂ p ' ∂ ζ ) e + ( g 1 3 ⟨ a ˜ P / ρ ⟩ ∂ p ' ∂ ζ ) o − ( g 1 3 ⟨ a ˜ P / ρ ⟩ ∂ p ' ∂ ξ ) t + ( g 1 3 ⟨ a ˜ P / ρ ⟩ ∂ p ' ∂ ξ ) b ![]()

E72
( ∂ p ' ∂ ξ ) P = 0 . 5 ( p E ' − p O ' ) ; ( ∂ p ' ∂ ζ ) P = 0 . 5 ( p T ' − p B ' ) ![]()

E73
F e , o = F e , o ∗ − [ ⟨ ρ a ˜ P ⟩ ( g 3 3 ∂ p ' ∂ ξ − g 1 3 ∂ p ' ∂ ζ ) ] e , o ![]()

E74
F t , b = F t , b ∗ − [ ⟨ ρ a ˜ P ⟩ ( g 1 1 ∂ p ' ∂ ζ − g 1 3 ∂ p ' ∂ ξ ) ] t , b ![]()

E75
F e , o ∗ = F e , o m 1 + E + F ^ e , o + ⟨ ρ a ˜ P ⟩ e , o ( g 3 3 ∂ p ∗ ∂ ξ − g 1 3 ∂ p ∗ ∂ ζ ) e , o ![]()

E76
F t , b ∗ = F t , b m 1 + E + F ^ t , b + ⟨ ρ a ˜ P ⟩ t , b ( g 1 1 ∂ p ∂ ξ − g 1 3 ∂ p ∂ ζ ) t , b ![]()

E77
F ^ e , o = z ζ ⟨ ρ u ^ ⟩ e , o − x ζ ⟨ ρ w ^ ⟩ e , o ![]()

E78
F ^ t , b = x ξ ⟨ ρ w ^ ⟩ t , b − z ζ ⟨ ρ u ^ ⟩ t , b ![]()

E79
a ˜ P u ^ P = ∑ n b a n b u ^ n b + a P E u P m − z ζ ∂ p ∗ ∂ ξ + z ξ ∂ p ∗ ∂ ζ + b ⇒ u ^ P = ∑ n b a n b u ^ n b + b a ˜ P ![]()

E80
a ˜ P w ^ P = ∑ n b a n b w ^ n b + a P E w ^ P m + x ζ ∂ p ∗ ∂ ξ − x ξ ∂ p ∗ ∂ ζ + b ⇒ w ^ P = ∑ n b a n b w ^ n b + b a ˜ P ![]()

E81
### 4.5. Solution of the equations

ϕ = f ϕ c o m p u t e d + ( 1 − f ) ϕ p r e v i o u s ![]()

E82
〚 R u , R w , R m , R T , R k , R ε 〛 < R c o n v ![]()

E83
R ϕ = Σ a l l ( | a P ϕ P − Σ n b a n b u n b + b | J a p ( ϕ max − ϕ min ) ) / Σ a l l J ![]()

E84

Discretizationand integration of the transport equations described previously are performed using a non-orthogonal generalized mesh as shown in Figure 1. The independent Cartesian coordinates (* x, z*), describing the physical domain, are thus replaced by a boundary-fitted coordinate system (

Transformation of the original equations is accomplished by replacing the independent variables, using the chain rule, which states that, generically:

The Jacobian of the transformation

represents the ratio between the physical size of the control volume and its computational size (unitary). The derivatives

To obtain the strong conservative form in the boundary-fitted coordinate system (* ξ, ζ*), the transport equations are transformed through the application of the chain rule (35) and multiplied by the Jacobian of the transformation. The metric identity

is then used to recast some terms. The result, for a generic variable * ϕ*, is

The terms ^{11}, ^{13} and ^{33} are the contravariant metric relations given by

The non-orthogonal term ^{13} is null if the mesh is locally orthogonal. * U* and

Note that, in this case, the sub-index for the fluxes (such as in

where the cross-derivatives were incorporated into the source term _{cross}. A similar procedure is applied to obtain the generalized form of remaining equations, leading to the following form for the Navier-Stokes and continuity equations:

The integration and solution method for the transport equations are entirely based on the methodology in [6], with some of the suggestions described in [7] and incorporating the necessary modifications for the generalized mesh approach. The general Eq. (43) may be written as

The integration of the previous equation in its control volume leads to

where the source term has been written as a linear combination involving * ϕ* and

In the previous expressions, the subscripts indicate the location relative to the * CV* centre, in the computational domain, with the uppercase denoting neighbour nodes and the lowercase denoting neighbour faces:

For the solution of the equations, it is necessary to evaluate the values of * ϕ* in the

or, in a more compact manner,

with “* nb*” indicating that the sum is to be performed for all the neighbouring locations.

It is necessary to compute the face values * ϕ*, being advected, tends to assume a value closer to the upwind value. Several advection schemes may be adopted, being the simplest one the

Due to its first-order character, the upwind scheme is often not used due to associated numerical diffusion. The Quick scheme is third-order accurate. The deferred correction version of Hayase [8] combines the first-order upwind scheme with a third-order correction, _{quick}, added to the source term as follows:

The Quick scheme, although third-order accurate, presents oscillations that may lead to some unrealistic behaviour. Total variation diminishing (TVD) schemes, also implemented in the present code, were developed to provide second-order accurate solutions that are free or nearlyfree from oscillations. For further information on this, please refer to, e.g., [9].

EasyCFD adopts the SIMPLEC algorithm (Semi-Implicit Method for Pressure-Linked Equations-Consistent) [7], which is based on the original formulation * SIMPLE* [6]. Due to the non-staggered mesh arrangement (collocated mesh), the Rie-Chow interpolation procedure [10], with the modifications proposed in [11] and [12], is implemented. Let us consider the

During the iterative process, velocities ^{*} are computed from the available velocity field ^{m} and pressure field ^{*} obtained at the previous step as follows:

where * E* is the under-relaxation factor [7]. During the iterative process, unless convergence is reached, the velocity field

Thus, the pressure and velocity fields ^{*} and ^{*} should be corrected by a certain amount * p*' and

Subtracting Eq. (56) from Eq. (55) and taking into account Eq. (57), one obtains

The keystone of the SIMPLECalgorithm consists on the subtraction of the term

or, for simplicity:

where

The equations for the velocity correction are obtained through the pressure correction field, recurring to the previous equation:

leading to

As previously stated, the objective of the pressure correction is to produce a pressure field such that the solution of the momentum equations is a mass-conservative velocity field. Consequently, the equations for solving the * p’* field must be obtained from the continuity equation. The discretized form of this equation is obtained directly from the integration of Eq. (2):

As one may see, velocities are, now, needed at the control volume faces. Taking Eq. (63), for the * u* velocity component, at the “

The terms

The derivatives are evaluated as

and, for the cross-derivatives,

Introducing the discretization expressed by Eqs. (68) and 67 into Eq. (65) allows us to obtain the pressure correction equation:

where

The * p’* obtained from the solution of Eq. (69) is employed for correcting pressure and velocity corrections are obtained from Eq. (63). Note that, since these equations are defined at the control volume centre,

One may note that pressure values at the control volume centre are not included in the evaluation of these derivatives (the same applies for the pressure derivatives in the momentum equations). This may lead to the well-known checkerboard pattern for the pressure field. To avoid this effect, the Chie-Row interpolation method proposes that the mass fluxes, to be evaluated at the control volume interfaces for all the transport equations (_{e}, _{o}, _{t} and _{b}, in Eq. (48)), be corrected using the pressure correction * p’* field (instead of being computed from the corrected control volume centre velocities). The correction equations for fluxes are obtained from the correction equations for velocities:

The “starred” fluxes ^{*} at the control volume interfaces are obtained by interpolating the momentum equation. The keystone of the method is that the pressure gradient is not interpolated, but, instead, is obtained directly from the pressure at contiguous control volume centres. Considering Eq. (55), evaluated in terms of fluxes, one may write

The terms

where the revised velocities * û* and

Note that the source term * b* includes all the contributions (e.g., transient term, cross-derivatives and buoyancy for the

The solution of the equations previously described is carried out with an iterative procedure. For accelerating the convergence rate, two relaxation factors (described next) are applied.

The solution of the equation is sub- or over-relaxed in the following manner:

Values of * f*lower than unity lead to sub-relaxation, while values greater than unity over-relax the solution process. In the present code, the value

The whole flow field calculation is considered to be converged when all the normalized residuals are lower that a predefined value

The total normalized residual for the transport equations of

where

Advertisement
## 5. Examples of application

### 5.1. Flow past an airfoil

### 5.2. Natural convection inside a cavity

Prandtl number , P r = v α ![]()

E85
Rayleigh number , R a = g β Δ T L 3 Pr v 2 ![]()

E86
u ∗ = u L α ; w ∗ = w L α ![]()

E87

To exemplify the numerical method described previously, two test cases are presented next along with a comparison of results with published data.

A calculation was performed to compute the aerodynamic coefficients of a NACA 0012 airfoil operating at a Reynolds number of 6×10^{6}. The obtained values for the drag and lift coefficients are compared with existing experimental data. The first step is to define the calculation domain, which should be large enough in order to avoid numerical blockage effects. Figure 2 represents the domain, for a 1 m airfoil cord length. Lateral boundaries are assigned a free slip condition and a uniform velocity profile with 5% turbulence intensity is imposed at the inlet. A mass conservative condition is applied at the outlet boundary. After a mesh independency study, a total of approximately 250,000 mesh nodes were employed, with three mesh refinement regions. Particular care was taken near the airfoil surface, were ^{+} values ranging from, typically, 0.1 to 6, with an average value of 1.7 all around, were obtained. Figure 3 depicts the mesh employed.

For the present simulations, both the first-order hybrid [6] and the Quick advection schemes were employed, along with the * k-ε* and

Results for the dependence of the lift coefficient with the airfoil angle of attack * α* are shown in Figure 4. As expected, the lift coefficient presents a linear dependence with the angle of attack

Figure 5 depicts the relation between lift and drag coefficients. In this case, the * k-ω* SST turbulence model performs substantially better than the

The natural convection flow in a cavity is a classical test for numerical methods in fluid dynamics. The cavity is a square shape (cf. Figure 6) with adiabatic horizontal walls. A constant temperature is imposed in each the vertical wall.

The problem is governed by the following non-dimensional parameters:

where the thermal diffusivity is * ν* is the kinematic viscosity,

where * β* is the thermal expansion coefficient,

Figure 7(a) and (b) displays isothermal lines generated using a constant value spacing between the minimum and the maximum verified within the domain. Figure 8(a) and (b) shows the flow streamlines. The flow, in the steady-state situation, is characterized by a large vortex filling the cavity, rotating in the clockwise direction. Two small vortices rotating in the same direction are located near the cavity centre. For this case, the minimum and the maximum streamline values used in the visualization do not correspond to the total amplitude of the stream function within the domain. These values were, instead, adjusted in EasyCFD to correspond to those employed in [19]. The agreement between the calculations and those reported in the literature is very good. Vahl Davis and Jones [18] present normalized maximum values for the * u* and for

occurring in the vertical and horizontal symmetry lines, respectively. Table 1 shows the results obtained with EasyCFD, the reference values in [17] and the range of variation for the * 37* contributions reported in [18]. This range does not include the minimum and maximum reported values since these clearly fall outside the general trend of the remaining contributions.

Advertisement
## 6. Concluding remarks

The numerical simulation of fluid flow for 2D problems was addressed. The physical principles and numerical models here presented correspond to the implementation in the software package EasyCFD. Transformation of the original equations to cope with a non-orthogonal generalized mesh is described in detail, along with the coupling of momentum and continuity with an adapted SIMPLEC algorithm for non-staggered meshes. Although not addressed in the present chapter, this software package was developed entirely based on a graphical interface, aiming at an easy and intuitive utilization. With a fast learning curve, this package is very suitable for learning the principles and application methods in computational fluid dynamics and has a great value both as a didactic and an applied tool. Although, at first, the restriction to 2D situations may seem very limitative, a great number of practical situations may be addressed with this approach.

- 1.
Lopes AMG. A 2D software system for expedite analysis of CFD problems in complex geometries. Computer Applications in Engineering Education. 2015; 3 (2):27–38. DOI: 10.1002/cae.21668. - 2.
Launder BE, Spalding DB. The numerical computation of turbulent flows. Computer Methods in Applied Mechanics and Engineering. 1974; 3 (2):269–289. DOI: 10.1016/0045-7825(74)90029-2. - 3.
Menter FR. Zonal two equation kappa-omega turbulence models for aerodynamic flows. AIAA Paper 1993–2906, 24th Fluid Dynamics Conference, 6–9 July 1993; Orlando, Florida. - 4.
Menter, Kuntz M, Langtry R. Ten years of industrial experience with the SST turbulence model. Turbulence, Heat and Mass Transfer. 2003;4:625–632. - 5.
Menter FR, Ferreira JC, Esch T. The SST turbulence model with improved wall treatment for heat transfer predictions in gas turbines. International Gas Turbine Congress. 2003;(1992):1–7. - 6.
Patankar S. Numerical Heat Transfer and Fluid Flow. Series in computational methods in mechanics and thermal sciences; 1980. New York, McGraw-Hill. - 7.
Van Doormaal JP, Raithby GD. Enhancements of the SIMPLE method for predicting incompressible fluid flows. Numerical Heat Transfer. 2007; 7 (2):147–163. DOI: 10.1080/01495728408961817. - 8.
Hayase T, Humphrey JA, Greif R. A consistently formulated QUICK scheme for fast and stable convergence using finite-volume iterative calculation procedures. Journal of Computational Physics. 1992; 98 (1):108–118. DOI: 10.1016/0021-9991(92)90177-Z. - 9.
Vertseeg HK, Malalasekera W. An Introduction to Computational Fluid Dynamics. The Finite Volume Method. Pearson Education Limited; Essex, England, 2007. - 10.
Rhie CM, Chow WL. Numerical study of the turbulent flow past an airfoil with trailing edge separation. AIAA Journal. 1983; 21 (11):1525–1532. DOI: 10.2514/3.8284. - 11.
Majumdar S. Role of underrelaxation in momentum interpolation for calculation of flow with nonstaggered grids. Numerical Heat Transfer. 1998; 13 (1):125–132. DOI: 10.1080/10407788808913607. - 12.
Shen WZ, Michelsen JA, Sørensen NN, Sørensen JN. An improved SIMPLEC method on collocated grids for steady and unsteady flow computations. Numerical Heat Transfer: Part B. 2003; 43 :221–239. DOI: 10.1080/713836202. - 13.
Ehrlich LW. An ad hoc SOR method. Journal of Computational Physics. 1981; 44 (1):31–45. DOI: 10.1016/0021-9991(81)90036-X. - 14.
Hutchinson BR, Raithby GD. A multigrid method based on the additive correction strategy. Numerical Heat Transfer. 1986; 9 (5):511–537. DOI: 10.1080/10407788608913491. - 15.
Abbott IH, Von Doenhoff AE. Theory of Wing Sections, Including a Summary of Airfoil Data. Dover Publications; New York, 1959. - 16.
Ladson CL. Effects of independent variation of Mach and Reynolds numbers on the low-speed aerodynamic characteristics of the NACA 0012 airfoil section. NASA TM 4074, October 1988. - 17.
Vahl Davis GD. Natural convection in a square cavity: a benchmark numerical solution. International Journal for Numerical Methods in Fluids. 1983; 3 :249–264. - 18.
Vahl Davis GD, Jones IP. Natural convection in a square cavity: a comparison exercise. International Journal for Numerical Methods in Fluids. 1983; 3 :227–248. - 19.
Dixit HN, Babu V. Simulation of high Rayleigh number natural convection in a square cavity using the lattice Boltzmann method. International Journal of Heat and Mass Transfer. 2006; 49 :727–739.

Submitted: October 23rd, 2015 Reviewed: March 14th, 2016 Published: August 24th, 2016

© 2016 The Author(s). Licensee IntechOpen. This chapter is distributed under the terms of the Creative Commons Attribution 3.0 License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.