Department of Aeronautical Engineering, KCT

1 UNIT – V

2 INTRODUCTION TO SIMPLE, PISO SCHEMES & EXAMPLE PROBLEMS

4 Pressure Correction Techniques (SIMPLE ALGORITHM)

5
∂ρ ∂(ρu) ∂(ρv) ∂(ρw)
6 + + + =0 (A)
∂t ∂x ∂y ∂z
∂P ∂2 u ∂2 u ∂2 u ∂u ∂u ∂u ∂u
7 ρg x − + μ ( 2 + 2 + 2) = ρ ( + u + v + w ) (B)
∂x ∂x ∂y ∂z ∂t ∂x ∂y ∂z
8

9 Introduction:

10 Generally, the inviscid incompressible flow is governed by elliptic nature based

11 Partial Differnetial Equations, which has been solved through an iterative prcoess

12 such as relation technique, etc. But most of us executing the incompressible viscous

13 flow based problems rather than other kinds of flows. Mathematically, viscous

14 incompressible is exhibit a mixed elliptic-parabolic nature so advanced techniques

15 such as pressure correction technique, etc have been used to solve this kind of

16 problems. This pressure correction technique has been widely used in numerical

17 simulation of incompressible Navier-Strokes Equation. The pressure technique is

18 embodied is an algorithm called SIMPLE [Semi-Implicit Method for Pressure Linked

19 Equations]. SIMPLE algorithm is applicable for both compressible and incompressible

20 flow. The concept called ‘average pressure’ is used to construct the next iteration from

21 the results of the presvious iteration.

22

1
Computational Fluid Dynamics
 Department of Aeronautical Engineering, KCT

23
24 Figure 0. SIMPLE selection facility at ANSYS FLUENT

25 If we go through Navier-Strokes equations, which contains the term pressure

26 gradient. But the representation of pressure gradient has been included in a depended

27 manner. Thereofre the seperate care needs to be taken in the estimation of pressure.

28 For that reason, the pressure correction technique [SIMPLE Algorithm] was emerged.

29 The complete philosophy of Presssure Correction Technique is reavelaed in Figure 1.

30

2
Computational Fluid Dynamics
 Department of Aeronautical Engineering, KCT

31
32 Figure [Link] of Pressure Correction Technique

33 In this case, we can take 2-D problem for our simplicity also finalized to neglect the

34 body forces, and thus the 2-D momentum equations are as follows:

35

∂P ∂2 u ∂2 u ∂u ∂u ∂u
36 − + μ ( 2 + 2) = ρ ( + u + v ) (1)
∂x ∂x ∂y ∂t ∂x ∂y
37

∂P ∂2 v ∂2 v ∂v ∂v ∂v
38 − + μ ( 2 + 2) = ρ ( + u + v ) (2)
∂y ∂x ∂y ∂t ∂x ∂y
39

40 Above equations follow directly from the model on an infinitely small volume

41 fixed in space. So finite difference form of above equations will be somewhat akin to

42 discretized equations obtained from finite volume approach. but the original work

43 involved finite volume approach.

44

45 In this case, instead of Fintie Volume approach, the finite difference method is

46 implemented in this equation. In which, forward difference scheme is used for the

47 representation of time derivatives and central difference scheme is used for the

48 representation of spatial derivatives.

3
Computational Fluid Dynamics
 Department of Aeronautical Engineering, KCT

49

50 Collocated Grid: where all flow variables are stored in the same positions.

51

52 Staggered Grid: Scalar variables [pressure, density, total enthalphy, etc] are stored in

53 the cell centres of the control volumes, whereas the velocity or momentum variables

54 are located at the cell faces.

55
56 Figure 2. The staggered grid in the direction of X

57

58 The pressure values are stored at the RED points and the velocity values are stored at

59 the BLUE points [Figure 2].

60

61 Note: Equiavlent FVM will deal with shaded cell [Figure 2]

62
1
63 vj+1 = 2 (vi,j+1 + vi+1,j+1 )
At points ‘a’: ̅̅̅̅̅
2 2 2

64
1
65 At points ‘b’: v̅̅̅̅̅
j−
1 = (vi,j−1 + vi+1,j−1 )
2 2 2 2

66

4
Computational Fluid Dynamics
 Department of Aeronautical Engineering, KCT

∂[ρu] ∂[ρu2 ] ∂[ρuv] ∂P ∂2 u ∂2 u

67 (1) ⇒ + + = − + μ ( 2 + 2)
∂t ∂x ∂y ∂x ∂x ∂y
68

∂[ρu] ∂[ρu2 ] ∂[ρuv] ∂P ∂2 u ∂2 u

69 = −[ + ]− + μ ( 2 + 2) (3)
∂t ∂x ∂y ∂x ∂x ∂y
70

71 Forward Difference Scheme for time derivative,

72

[ρu]t+∆t
1 − [ρu]
t
1
∂[ρu] i+ ,j
2
i+ ,j
2
73 = (4)
∂t ∆t
74

75 Central Difference Scheme for spatial derivative,

[ρu2 ]t 3 − [ρu2 ]t 1
∂[ρu2 ] i+ ,j
2
i− ,j
2
76 = (5)
∂x 2∆x
77

[ρuv]t 1 − [ρuv]t 1
∂[ρuv] i+ ,j+1
2
i+ ,j−1
2
78 = (6)
∂y 2∆y
79

80 Forward Difference Scheme for pressure derivative,

81

∂P [p]ti+1,j − [p]ti,j
82 = (7)
∂x ∆x
83

84 Second order Finite Difference derivative for higher order spatial derivative,

85

ut 3 − 2ut 1 + ut 1
∂2 u i+ ,j
2
i+ ,j
2
i− ,j
2
86 = (8)
∂x 2 [∆x]2
87

2 ut 1 − 2ut 1 + ut 1
∂ u i+ ,j+1
2
i+ ,j
2
i+ ,j−1
2
88 = (9)
∂y 2 [∆y]2

5
Computational Fluid Dynamics
 Department of Aeronautical Engineering, KCT

89

90 Substitue the equations (4), (5), (6), (7), (8), and (9) in equation (3),

91

[ρu]t+∆t
1 − [ρu]
t
1
i+ ,j i+ ,j
2 2
92
∆t
[ρu2 ]t 3 − [ρu2 ]t 1 [ρuv]t 1 − [ρuv]t 1
i+ ,j i− ,j i+ ,j+1 i+ ,j−1
2 2 2 2
93 = −[ + ]
2∆x 2∆y

[p]ti+1,j − [p]ti,j
94 −[ ]
∆x
ut 3 − 2ut 1 + ut 1 ut 1 − 2ut 1 + ut 1
i+ ,j i+ ,j i− ,j i+ ,j+1 i+ ,j i+ ,j−1
2 2 2 2 2 2
95 + μ[ + ] (10)
[∆x]2 [∆y]2

96

97 For our simplicity let us take the big component of equation (10) is equal to A,

98

[ρu2 ]t 3 − [ρu2 ]t 1 [ρuv]t 1 − [ρuv]t 1

i+ ,j i− ,j i+ ,j+1 i+ ,j−1
2 2 2 2
99 A = −[ + ]
2∆x 2∆y

ut 3 − 2ut 1 + ut 1 ut 1 − 2ut 1 + ut 1
i+ ,j i+ ,j i− ,j i+ ,j+1 i+ ,j i+ ,j−1
2 2 2 2 2 2
100 +μ[ + ]
[∆x]2 [∆y]2

101 Equation (10) becomes,

102

[ρu]t+∆t
1 − [ρu]
t
i+ ,j
2
1
i+ ,j
2
[p]ti+1,j − [p]ti,j
103 =A−[ ]
∆t ∆x
104
∆t
105 [ρu]t+∆t
1 = [ρu]
t
1 + A∆t − {( ) [[p]ti+1,j − [p]ti,j ]} (11)
i+ ,j
2
i+ ,j
2 ∆x
106

6
Computational Fluid Dynamics
 Department of Aeronautical Engineering, KCT

107
108 Figure 3. The staggered grid in the direction of Y

109

110 Note: Equiavlent FVM will deal with shaded cell [Figure 3]

111
1
112 At points ‘c’: u
̅̅̅̅̅
i−
1 = (ui−1,j + ui−1,j+1 )
22 2 2

113
1
114 ui+1 = 2 (ui+1,j + ui+1,j+1 )
At points ‘d’: ̅̅̅̅̅
2 2 2

115

116 Similar from equation (11), we can also form the equation for ‘y’ momentum,

117
∆t
118 [ρv]t+∆t1 = [ρv]t 1 + B∆t − {( ) [[p]ti,j+1 − [p]ti,j ]} (12)
i ,j+
2
i ,j+
2 ∆y
119

[ρv 2 ]t 3 − [ρv 2 ]t 1 [ρuv]t 1 − [ρuv]t 1

i ,j+ i ,j− i+1 ,j+ i−1 ,j+
2 2 2 2
120 B = −[ + ]
2∆y 2∆x

vt 3 − 2v t 1 + vt 1 vt 1 − 2v t 1 + vt 1
i ,j+ i,j+ i ,j− i+1 ,j+ i,j+ i−1,j+
2 2 2 2 2 2
121 +μ[ + ]
[∆x]2 [∆y]2

122

7
Computational Fluid Dynamics
 Department of Aeronautical Engineering, KCT

∆t
123 (11) ⇒ [ρu]t+∆t
1 = [ρu]
t
1 + A∆t − {( ) [[p]ti+1,j − [p]ti,j ]}
i+ ,j
2
i+ ,j
2 ∆x
124

125 We know that, the equation (11) has been estimated through the assumption of guess

126 pressure, so the equation (11) becomes,

127
∆t
128 [ρu∗ ]t+∆t ∗ t
1 = [ρu ] 1 + A∗ ∆t − {( ) [[p∗ ]ti+1,j − [p∗ ]ti,j ]} (13)
i+ ,j
2
i+ ,j
2 ∆x
129

130 From the philosopy of SIMPLE algorithm,

131

132 p = p∗ + p′ ⇒ p′ = p − p∗ (14)
133
∆t
134 [ρu′ ]t+∆t ′ t
1 = [ρu ] 1 + A′ ∆t − {( ) [[p′ ]ti+1,j − [p′ ]ti,j ]} (15)
i+ ,j
2
i+ ,j
2 ∆x
135

136 [ρu′ ]t+∆t t+∆t ∗ t+∆t

1 = [ρu] 1 − [ρu ] 1 (16)
i+ ,j i+ ,j i+ ,j
2 2 2

137

138 [ρu′ ]t 1 = [ρu]t 1 − [ρu∗ ]t 1 (17)

i+ ,j i+ ,j i+ ,j
2 2 2

139

140 From the equation (14), we can also derive author paramters,

141

142 A = A∗ + A′ ⇒ A′ = A − A∗ (18)
143

144 [p′ ]ti+1,j = [p]ti+1,j − [p∗ ]ti+1,j (19)

145

146 [p′ ]ti,j = [p]ti,j − [p∗ ]ti,j (20)

147

148 Similarly, we can expand the above equation for ‘y’ Momentum,

149
8
Computational Fluid Dynamics
 Department of Aeronautical Engineering, KCT

150 B = B∗ + B′ ⇒ B′ = B − B∗ (21)
151

152 [p′ ]ti,j+1 = [p]ti,j+1 − [p∗ ]ti,j+1 (22)

153
∆t
154 [ρv ∗ ]t+∆t1 = [ρv ∗ ]t 1 + B∗ ∆t − {( ) [[p∗ ]ti,j+1 − [p∗ ]ti,j ]} (23)
i,j+
2
i ,j+
2 ∆y
155
∆t
156 (12) ⇒ [ρv]t+∆t1 = [ρv]t 1 + B∆t − {( ) [[p]ti,j+1 − [p]ti,j ]}
i ,j+
2
i ,j+
2 ∆y
157
∆t
158 [ρv ′ ]t+∆t1 = [ρv ′ ]t 1 + B′ ∆t − {( ) [[p′ ]ti,j+1 − [p′ ]ti,j ]} (24)
i,j+
2
i,j+
2 ∆y
159

160 [ρv ′ ]t+∆t1 = [ρv]t+∆t1 − [ρv ∗ ]t+∆t1 (25)

i,j+ i,j+ i,j+
2 2 2

161

162 [ρv ′ ]t 1 = [ρv]t 1 − [ρv ∗ ]t 1 (26)

i,j+ i,j+ i,j+
2 2 2

163

164 The ultimate aim of this SIMPLE algorithm is to provide convergence of this CFD

165 computation. In order to attain convergence, the value of pressure correction factor

166 ‘p′ ’ sholud be tends to zero. [p′ → 0].

167

168 If [p′ → 0] means, arbitrarly, the values of A′ , B ′ , [ρv ′ ]t , and [ρu′ ]t are equal to zero.

169

170 Already We derived that,

∆t
171 (15) ⇒ [ρu′ ]t+∆t ′ t
1 = [ρu ] 1 + A′ ∆t − {( ) [[p′ ]ti+1,j − [p′ ]ti,j ]}
i+ ,j
2
i+ ,j
2 ∆x
∆t
172 (24) ⇒ [ρv ′ ]t+∆t1 = [ρv ′ ]t 1 + B′ ∆t − {( ) [[p′ ]ti,j+1 − [p′ ]ti,j ]}
i,j+
2
i,j+
2 ∆y
173 Apply the boundary conditions and neglection concept in the equations (15), and (24),

174 the new mosidified equations as,

9
Computational Fluid Dynamics
 Department of Aeronautical Engineering, KCT

175
∆t
176 [ρu′ ]t+∆t
1 = − {( ) [[p′ ]ti+1,j − [p′ ]ti,j ]} (27)
i+ ,j
2 ∆x
177
∆t
178 [ρv ′ ]t+∆t1 = − {( ) [[p′ ]ti,j+1 − [p′ ]ti,j ]} (28)
i,j+
2 ∆y
179 Subsittute the Equation (27) in Equation (16) and Equation (28) in Equation (25),

180

181 Equation (16) becomes,

∆t
182 − {( ) [[p′ ]ti+1,j − [p′ ]ti,j ]} = [ρu]t+∆t ∗ t+∆t
1 − [ρu ] 1
∆x i+ ,j
2
i+ ,j
2
∆t
183 [ρu]t+∆t ∗ t+∆t
1 = [ρu ] 1 − {( ) [[p′ ]ti+1,j − [p′ ]ti,j ]} (29)
i+ ,j
2
i+ ,j
2 ∆x
184

185 Equation (25) becomes,

∆t
186 − {( ) [[p′ ]ti,j+1 − [p′ ]ti,j ]} = [ρv]t+∆t1 − [ρv ∗ ]t+∆t1
∆y i,j+
2
i,j+
2
∆t
187 [ρv]t+∆t1 = [ρv ∗ ]t+∆t1 − {( ) [[p′ ]ti,j+1 − [p′ ]ti,j ]} (30)
i,j+
2
i,j+
2 ∆y
188 From the NS equation, we know that,

189

190 The continuity equation is,

∂(ρu) ∂(ρv)
191 + =0 (31)
∂x ∂y
192 The differential equation of (31) can be easily converted into difference equation with

193 the help of central finite difference scheme at the point of (i, j),

194

[ρu]t+∆t t−∆t
1 − [ρu] 1 [ρv]t+∆t1 − [ρv]t−∆t1
i+ ,j i− ,j i,j+ i,j−
2 2 2 2
195 + =0 (32)
2∆x 2∆y
196

197 With the help of Equations (29) and (30) are we can modify the equation (32),

10
Computational Fluid Dynamics
 Department of Aeronautical Engineering, KCT

∆t ∆t
[ρu∗ ]t+∆t ′ t ′ t ∗ t−∆t ′ t ′ t
1 − {(∆x) [[p ]i+1,j − [p ]i,j ]} − [ρu ] 1 + {(∆x) [[p ]i,j − [p ]i−1,j ]}
i+ ,j i− ,j
2 2
198
2∆x
∆t ∆t
[ρv ∗ ]t+∆t1 − {( ′ ]t [p′ ]ti,j ]} − [ρv ∗ ]t−∆t1 + {( ) [[p′ ]ti,j − [p′ ]ti,j−1 ]}
i,j+
2
∆y) [[p i,j+1 − i,j−
2
∆y
199 + =0 (33)
2∆y
200

201 The equation (33) is very complicated one so, for simplicity the equation (33) is

202 becomes,
203 ap′i,j + bp′i+1,j + bp′i−1,j + cp′i,j+1 + cp′i,j−1 + d = 0 (34)

204

205 Where,
∆t ∆t ∆t ∆t 1 1
206 𝑎 = 2 [∆x2 + ∆y2 ] ; 𝑏 = − ∆x2 ; 𝑐 = − ∆y2 ; 𝑑 = ∆x [[ρu∗ ]t+∆t ∗ t−∆t
1 − [ρu ] 1 ] +
i+ ,j i− ,j
[[ρv ∗ ]t+∆t
i,j+
1 −
∆y
2 2 2

207 [ρv ∗ ]t−∆t1 ] ;

i,j−
2

208

209 The equation (34) is commonly known as Pressure Correction Equation, in which all

210 the coeeficents are known to compuer so with the help of boundary conditions we can

211 easily calculate the pressure correction [p′i,j , p′i+1,j , p′i−1,j , p′i,j+1 , p′i,j−1 ] at various locations

212 have been estimated. If we know the value of pressure correction then we can easily

213 estimate the ‘ corrected pressure’ value, which supports the calculation of velocities

214 through continuity and momentum equations.

215

216 References

217

218 1. [Link]

219

220

221

222

11
Computational Fluid Dynamics
 Department of Aeronautical Engineering, KCT

223 PISO [Pressure Implicit with Split Operator] – Algorithm

224 The PISO stands for Pressure Implicit with Splitting Operators.

225 PISO is a pressure – velocity calculation procedure developed originally for non-

226 iterative computation of unsteady compressible flows.

227 The PISO also has been adapted successfully for the iterative solution of steady state

228 problems.

229 PISO involves one predictor step and two corrector steps and may be seen as an

230 extension of SIMPLE, with a further corrector step to enhance it.

231

232
233 Figure 1. Flowchart of PISO Algorithm [Source: Wikipedia]

234

235 Procedures Involved PISO Scheme

12
Computational Fluid Dynamics
 Department of Aeronautical Engineering, KCT

236

237 [1]. Assume the Guess Pressure [p∗ ]

238 [2]. Substitute the guess pressure [p∗ ] in the momentum equation and thereby

239 estimate the ‘u∗ ’, ‘v ∗ ’and ‘w ∗ ’.

240 [3]. Substitute the values of ‘u∗ ’, ‘v ∗ ’and ‘w ∗ ’ in the continuity equation and thereby

241 estimate the first pressure correction [p′ ].

242 [4]. Update the corrected pressure ‘p∗∗ = p∗ + p′ ’

243 [5]. Substitute the corrected pressure [p∗∗ ] in the momentum equation and thereby

244 estimate the velocities of ‘u∗∗ ’, ‘v ∗∗ ’and ‘w ∗∗ ’.

245 [6]. Substitute the values of ‘u∗∗ ’, ‘v ∗∗ ’and ‘w ∗∗ ’in the continuity equation and

246 thereby estimate the second pressure correction [p′′ ].

247 [7]. Update the corrected pressure ‘p∗∗∗ = p∗∗ + p′′ ⇒ p∗ + p′′ + p′ ’

248 [8]. Substitute the corrected pressure [p∗∗∗ ] in the momentum equation and thereby

249 estimate the velocities of ‘u∗∗∗ ’, ‘v ∗∗∗ ’and ‘w ∗∗∗ ’.

250 [9]. The consolidated outputs of ‘ p∗∗∗ ’, ‘u∗∗∗ ’, ‘v ∗∗∗ ’and ‘w ∗∗∗ ’are taken as final

251 outputs in the case of unsteady compressible flows whereas in the case of steady

252 state problems, the iteration are needs to be executed until correction factor tends

253 to be zero.

254

255 Predictor Step

256 Discretized momentum equations (11) – (12) are solved with a guessed or intermediate

257 pressure field ‘p∗ ’ to give velocity components ‘u∗ ’ and ‘v ∗ ’ using the same method as

258 the SIMPLE algorithm.

259

260 Corrector Step 1

261 The ‘u∗ ’ and ‘v ∗ ’ fields will not satisfy continuity unless the pressure field ‘p∗ ’ is

262 correct. The first corrector step of SIMPLE is introduced to give a velocity field (‘u∗∗ ’

263 and ‘v ∗∗ ’), which satisfies the discretized continuity equation. The resulting equations

13
Computational Fluid Dynamics
 Department of Aeronautical Engineering, KCT

264 are the same as the velocity correction equations (29)-(30) of SIMPLE but, since there

265 is a further correction step in the PISO algorithm.

266 We use a slightly different notation,

267

268 p∗∗ = p∗ + p′ (35)

269

270 u∗∗ = u∗ + u′ (36)

271

272 v ∗∗ = v ∗ + v ′ (37)
273

274 These equations (35), (36), and (37) are used to define corrected velocities ‘u∗∗ ’ and

275 ‘v ∗∗ ’.

276

277 From SIMPLE algorithm derivation,

∆t
278 (29) ⇒ [ρu]t+∆t ∗ t+∆t
1 = [ρu ] 1 − {( ) [[p′ ]ti+1,j − [p′ ]ti,j ]}
i+ ,j
2
i+ ,j
2 ∆x
∆t
279 (30) ⇒ [ρv]t+∆t1 = [ρv ∗ ]t+∆t1 − {( ) [[p′ ]ti,j+1 − [p′ ]ti,j ]}
i,j+
2
i,j+
2 ∆y
280 The above two equations can be modified as,
∆t
281 [ρu∗∗ ]t+∆t ∗ t+∆t
1 = [ρu ] 1 − {( ) [[p′ ]ti+1,j − [p′ ]ti,j ]} (38)
i+ ,j
2
i+ ,j
2 ∆x
282
∆t
283 [ρv ∗∗ ]t+∆t1 = [ρv ∗ ]t+∆t1 − {( ) [[p′ ]ti,j+1 − [p′ ]ti,j ]} (39)
i,j+
2
i,j+
2 ∆y
284

285 As in the SIMPLE algorithm, the above equations (38) and (39) are substituted into the

286 discretized equation (32) to yield pressure correction equation (34) with its coefficients

287 and source terms.

[ρu]t+∆t t−∆t
1 − [ρu] 1 [ρv]t+∆t1 − [ρv]t−∆t1
i+ ,j i− ,j i,j+ i,j−
2 2 2 2
288 (32) ⇒ + =0
2∆x 2∆y
289 Comparative term analysis:
14
Computational Fluid Dynamics
 Department of Aeronautical Engineering, KCT

Terms in the SIMPLE Terms in the PISO algorithm

algorithm

Corrected Pressures are Corrected Pressures are

[ρu]t+∆t
1 , [ρu∗∗ ]t+∆t
1
i+ ,j i+ ,j
2 2

[ρu]t−∆t
1 , [ρu∗∗ ]t−∆t
1
i− ,j i− ,j
2 2

[ρv]t+∆t1 , [ρv ∗∗ ]t+∆t1

i,j+ i,j+
2 2

[ρv]t−∆t1 [ρv ∗∗ ]t−∆t1

i,j− i,j−
2 2

Both the terms are analytically same magnitude but the only

difference is naming. In PISO scheme, the simply represent the

corrected pressure through ‘double star’ on superscripts.

290
∆t ∆t
[ρu∗ ]t+∆t ′ t ′ t ∗ t−∆t ′ t ′ t
1 − {(∆x) [[p ]i+1,j − [p ]i,j ]} − [ρu ] 1 + {(∆x) [[p ]i,j − [p ]i−1,j ]}
i+ ,j i− ,j
2 2
291
2∆x
∆t ∆t
[ρv ∗ ]t+∆t1 − {( ′ ]t [p′ ]ti,j ]} − [ρv ∗ ]t−∆t1 + {( ) [[p′ ]ti,j − [p′ ]ti,j−1 ]}
i,j+
2
∆y) [[p i,j+1 − i,j−
2
∆y
292 + =0 (40)
2∆y
293

294 (34) ⇒ ap′i,j + bp′i+1,j + bp′i−1,j + cp′i,j+1 + cp′i,j−1 + d = 0 [

295 ∵ 𝐀𝐬 𝐩𝐞𝐫 𝐏𝐈𝐒𝐎 𝐦𝐞𝐭𝐡𝐨𝐝, 𝐭𝐡𝐢𝐬 𝐞𝐪𝐮𝐚𝐭𝐢𝐨𝐧 𝐢𝐬 𝐜𝐚𝐥𝐥𝐞𝐝 𝐚𝐬 𝐟𝐢𝐫𝐬𝐭 𝐩𝐫𝐞𝐬𝐬𝐮𝐫𝐞 𝐜𝐨𝐫𝐫𝐞𝐜𝐭𝐢𝐨𝐧 𝐞𝐪𝐮𝐚𝐭𝐢𝐨𝐧]
296

297 ap′i,j + bp′i+1,j + bp′i−1,j + cp′i,j+1 + cp′i,j−1 + d = 0 (41)

298

299 The equations (41) and (34) are analytically same that is why the same equations are

300 used for the both the schemes. The equation (41) is solved to yield the first pressure

301 correction field ‘p′ ’. Once the ‘pressure corrections’ are known, the velocity

302 components ‘u∗∗ ’ and ‘v ∗∗ ’ can be obtained through the equations (38) and (39).
15
Computational Fluid Dynamics
 Department of Aeronautical Engineering, KCT

303

304 Corrector Step 2

305

306 To enhance the SIMPLE procedure, the PISO performs a second corrector step.

307 After the estimation of the velocities ‘u∗∗ ’, ‘v ∗∗ ’and ‘w ∗∗ ’, this needs to be submitted in

308 the continuity equation. And thus we can easily estimate the value of second pressure

309 correction [p′′ ].

310

311 p∗∗∗ = p∗∗ + p′′ ⇒ p∗ + p′′ + p′ (42)

312

313 At now, we know the values of first pressure correction and second pressure

314 correction, so with the help of equation (42), we can easily update the pressure ‘p∗∗∗ ’,

315 wherein the initial guess pressure ‘p∗ ’still maintained without any change.

316

317 As per the algorithm, the updated pressure ‘p∗∗∗ ’ needs to be submitted in momentum

318 equation and thus the updated velocities are as follows:

319
∆t
320 [ρu∗∗∗ ]t+∆t ∗∗ t+∆t
1 = [ρu ] 1 − {( ) [[p′′ ]ti+1,j − [p′′ ]ti,j ]} (43)
i+ ,j
2
i+ ,j
2 ∆x
321
∆t
322 [ρv ∗∗∗ ]t+∆t1 = [ρv ∗∗ ]t+∆t1 − {( ) [[p′′ ]ti,j+1 − [p′′ ]ti,j ]} (44)
i,j+
2
i,j+
2 ∆y
323

324 Substitute the equation (38) in (43) and (39) in (44),

325
∆t ∆t
326 [ρu∗∗∗ ]t+∆t ∗ t+∆t
1 = [ρu ] 1 − {( ) [[p′ ]ti+1,j − [p′ ]ti,j ]} − {( ) [[p′′ ]ti+1,j − [p′′ ]ti,j ]} (45)
i+ ,j
2
i+ ,j
2 ∆x ∆x
327
∆t ∆t
328 [ρv ∗∗∗ ]t+∆t1 = [ρv ∗ ]t+∆t1 − {( ) [[p′ ]ti,j+1 − [p′ ]ti,j ]} − {( ) [[p′′ ]ti,j+1 − [p′′ ]ti,j ]} (46)
i,j+
2
i,j+
2 ∆y ∆y

16
Computational Fluid Dynamics
 Department of Aeronautical Engineering, KCT

329

330 Substitute the above values in the discretized continuity equation (32), if the given

331 problem is steady state problem. Whereas the given problem is unsteady problem, the

332 equation (45), (46) and (42) are considered as outcomes [according to non-iterative

333 principle].

334 Introduction to Turbulence Modelling

335

336 For low fluid velocities, regular orderly motion of the fluid elements are fluid

337 dynamically called as Laminar Flow.

338 Random disordered motion in the fluid due to the increment of velocity is

339 qualitatively called as Turbulent Flow.

340 Basic entities in the turbulent flow are Eddies, which are commonly called as Lams of

341 rotating fluid masses.

342

343 Some important features of turbulent flow

344

345 ➢ Randomness of transport variables such as velocity, pressure, etc have been

346 generated with respect to time and also space.

347 ➢ Strong mixing is there in turbulent flow because of this interaction of the

348 fluctuation components of velocities.

349 ➢ Wind range of length scale and time scale, which makes the modelling of the

350 turbulence very complicated

351

352 Step by Step Process involved in the complete cycle of turbulence

353

354 Step – 1: Extracts Kinetic Energy from mean flow

355 Step – 2: Formation of Large Eddy

356 Step – 3: Formation of smaller Eddy

17
Computational Fluid Dynamics
 Department of Aeronautical Engineering, KCT

357 Step – 4: Formation of further smaller Eddy

358 Step – 5: Formation of further and further smaller Eddy

359 Step – 6: Formation of smallest Eddy

360 Step – 7: Energy mopped up due to Viscous Dissipation

361

362 For the large eddies, the inertial force is more dominant than viscous force that is why,

363 whatever energy extracted from the mean flow that cannot be dissipated by the large

364 eddy in the form viscous dissipation. These aforementioned seven steps are simply

365 known as energy cascading.

366

367 To understand the complication involved in the turbulence modelling, we will go for

368 small derivation, which will help us lot to appreciate the various turbulences and its

369 capturing complexity.

370

371 Possible of formation of Turbulence sizes are from largest to smallest so the order of

372 analysis based derivation is important in this section. [Note: Already we did order of

373 analysis in boundary layer equation]

374

Assumptions [Denotations] for order of analysis

Terms involved in the Largest Eddy Smallest

derivation Eddy

Length-Scale L 𝑙

Velocity-Scale U u

L 𝑙
Time-Scale
U 𝑢

375

18
Computational Fluid Dynamics
 Department of Aeronautical Engineering, KCT

376 Information – 1: Generally, the length scale of largest eddy is depends upon the system

377 length scale [simply we can dimensions of the overall control volume]

378

379 Let us take Rate of extraction of turbulent Kinetic Energy from mean flow per unit

380 mass is ‘π’.

1
𝑚U 2
381 𝜋~ [ 2 ]
𝑚
𝑝𝑒𝑟 𝑢𝑛𝑖𝑡 𝑡𝑖𝑚𝑒

382
1 2 1 2
U
383 𝜋~ 2 ⇒𝜋~ 2U
𝑡𝑖𝑚𝑒 𝑠𝑐𝑎𝑙𝑒 𝑡𝑖𝑚𝑒 𝑠𝑐𝑎𝑙𝑒
384

385 Here the time scale is based on the length and velocity scales of the eddy.
L
386 Time scale for Largest Eddy =
U
387 And this time scale is typically called as turnover time scale, which is a time scale

388 characteristics to the large eddy turn over.

1 2
U
389 𝜋~2
L
U
U3
390 𝜋~ (1)
L
391

392 [Note: in the order of time scale we have been not considered half and full values,

393 positive and negative values, the magnitude is the major concern]

394

395 Rate of Dissipation of Turbulent Kinetic Energy (𝜀),

396 𝜀 ~ 𝜐 𝑒𝑖𝑗 𝑒𝑖𝑗
397 In which, the ‘𝑒𝑖𝑗 ’ is rate of deformation tensor associated with Smallest Eddy, so this

398 rate of deformation is given by the velocity gradient in terms of the scale. And the

399 term ‘𝜐’ is represented the kinematic viscosity.

19
Computational Fluid Dynamics
 Department of Aeronautical Engineering, KCT

u
400 𝑒𝑖𝑗 ~
𝑙
401 Therefore,
u2
402 𝜀~𝜐 2 (2)
𝑙
403

404 Information – 2: Generally, the inertial force is more dominant in large eddy so the

405 kinetic energy extraction process is executed through Large Eddy. Whereas, in the

406 smallest eddy the viscous force more dominant therefore the viscous dissipation

407 process is executed through smallest eddy.

408

409 The energy extract from mean flow for the formation of large eddy is must be

410 transferred into small eddies and thereby it must be dissipated through the physics of

411 viscous dissipation.

412

413 Therefore,

414

415 𝜋~𝜀 (3)

416

417 If the viscous force of the eddy dominate high means, the eddies have been dissipated

418 surely so we will implement the order of analysis like this,

419
𝐼𝑛𝑒𝑟𝑡𝑖𝑎𝑙 𝐹𝑜𝑟𝑐𝑒
420 ~1
𝑉𝑖𝑠𝑐𝑜𝑢𝑠 𝐹𝑜𝑟𝑐𝑒
421

422 [𝑅𝑒𝑦𝑛𝑜𝑙𝑑𝑠 𝑁𝑢𝑚𝑏𝑒𝑟]𝑅𝑒𝑆𝑚𝑎𝑙𝑙𝑒𝑠𝑡 𝐸𝑑𝑑𝑦 ~ 1

423
𝑢∗𝑙
424 ~1 (4)
𝜐
425

20
Computational Fluid Dynamics
 Department of Aeronautical Engineering, KCT

426 Relationship between Smallest Eddy Length Scale, Rate of Dissipation and Kinematic

427 Viscosity:

428

429 From (2), we can write,

430

u2 𝜐
431 (2) ⇒ 𝜀 ~ 𝜐 2
⇒ 𝑙 2 ~ u2
𝑙 𝜀
𝜐 2
432 𝑙2~ u (5)
𝜀
433 From (4), we can write,
𝑢∗𝑙 𝜐
434 (4) ⇒ ~1⇒𝑢~
𝜐 𝑙
𝜐
435 𝑢~ (6)
𝑙
436 Substitute the equation (6) in equation (5),

437

2
𝜐 𝜐2
438 𝑙 ~
𝜀 𝑙2
439

4
𝜐3
440 𝑙 ~
𝜀
441
1
𝜐3 4
442 𝑙~ { } (7)
𝜀
443 From this comparative analysis, it is understood that the length scale of Smallest Eddy

444 is fully depends on fluid’s kinematic viscosity and Rate of Dissipation of Turbulent

445 Kinetic Energy. Also this length scale is known as Kolmogorov Length Scale.

446

447 Relationship between Largest Eddy Length Scale and Smallest Eddy Scale

448 From (3),

449

450 𝜋~𝜀

21
Computational Fluid Dynamics
 Department of Aeronautical Engineering, KCT

451 From (1) and (2),

U3 u2
452 ~𝜐 2
L 𝑙
453 From (6),
𝜐2
U3 [ ]
454 ~ 𝜐 𝑙2
L 𝑙
𝜐3L
455 𝑙4~ (8)
U3
456

457 Divide the equation (8) by ‘L4 ’ on both of its sides,

458

𝑙4 𝜐3L
459 ~
L4 U 3 L4
𝑙4 𝜐3
460 ~
L4 U 3 L3
𝑙4 1
461 4
~ 3 3
L U L
[ 3 ]
𝜐
462 From the definition of Reynolds number, we can further modify the above equation,
𝑙4 1
463 ~
L4 [[𝑅𝑒𝑦𝑛𝑜𝑙𝑑𝑠 𝑁𝑢𝑚𝑏𝑒𝑟]𝑅𝑒 3 ]
𝐿𝑎𝑟𝑔𝑒 𝐸𝑑𝑑𝑦

𝑙4
464 ~ 𝑅𝑒 −3 𝐿𝑎𝑟𝑔𝑒 𝐸𝑑𝑑𝑦 (9)
L4
465

466 From equation (9), we observed that, the Relationship between Largest Eddy Length

467 Scale and Smallest Eddy Scale is completely depends upon the Reynolds Number of

468 Large Eddy with respect to its system length scale.

469

470 References

471 1. [Link]

472

473

22
Computational Fluid Dynamics
 Department of Aeronautical Engineering, KCT

474 Reynolds Average Navier Strokes (RANS) Equation

475

476 Continuity Equation [1-D]

∂ui
477 =0 (1)
∂xi
478 Substituting ui = u̅i + u′ i in the Equation (1) and taking the average of the Equation,
∂[u̅i + u′ i ]
479 =0
∂xi
480 Average of above equation is,
̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅
∂[u̅i ] ∂[u′ i ]
481 + =0
∂xi ∂xi
̅̅̅̅̅̅
∂[u̅i ] ̅̅̅̅̅̅̅
∂[u′ i ]
482 + =0 (2)
∂xi ∂xi
483 From the basic average rule,
̅̅̅̅̅̅̅
∂[u′ i ]
484 =0 (3)
∂xi
485 [Note – 1: Any turbulent fluctuation quantity will have a zero average ]

486 Substitute the Equation (3) in Equation (2),

̅̅̅̅̅̅
∂[u̅i ]
487 =0
∂xi
∂[u̅i ]
488 =0 (4)
∂xi
489

490 Equation (4) is Continuity Equation for the mean flow.

491

492 Subtract Equation (4) from Equation (1),

∂ui ̅̅̅̅̅̅
∂[u̅i ]
493 − =0
∂xi ∂xi
494 From the basis of RANS principle, we can say,

∂ui ̅̅̅̅̅̅
∂[u̅i ] ∂[u′ i ]
495 − = =0
∂xi ∂xi ∂xi
496 Thus,

23
Computational Fluid Dynamics
 Department of Aeronautical Engineering, KCT

∂[u′ i ]
497 =0 (5)
∂xi
498 Equation (5) is Continuity Equation for the Turbulent Fluctuation Field.

499

500 Similarly we can expand the Equations (4) and (5) for 3-D,

501

502 Continuity Equation [3-D] for the mean flow

∂[u̅] ∂[v̅] ∂[w
̅]
503 + + =0
∂x ∂y ∂z
504 Continuity Equation [3-D] for the Turbulent Fluctuation Field
∂[u′ ] ∂[v ′ ] ∂[w ′ ]
505 + + =0
∂x ∂y ∂z
506

507 Momentum Equation in ‘i’ direction

∂ ∂ui ∂P ∂ui ∂ui
508 [μ ( )] − = ρ( + ui )
∂xi ∂xj ∂xi ∂t ∂xj
509

∂ ∂ui ∂P ∂ui ∂[ui uj ]

510 [μ ( )] − = ρ( + ) (6)
∂xj ∂xj ∂xi ∂t ∂xj
511

512 Substitute uj = u̅j + u′ j ; ui = u̅i + u′ i and P = ̅

P + P′

513

∂ ∂u̅i + u′ i ̅ + P′]
∂[P ∂[u̅i + u′ i ] ∂[{u̅i + u′ i }{u̅j + u′ j }]
514 [μ ( )] − = ρ( + )
∂xj ∂xj ∂xi ∂t ∂xj
515

516 Take time average on the both side of the above equation,

517

̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅
∂ ∂u̅i + u′ i ̅̅̅̅̅̅̅̅̅̅̅̅
∂[P̅ + P′] ̅̅̅̅̅̅̅̅̅̅̅̅̅ ̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅
∂[u̅i + u′ i ] ∂[{u ̅i + u′ i }{u̅j + u′ j }]
518 [μ ( )] − = ρ( + )
∂xj ∂xj ∂xi ∂t ∂xj

519

24
Computational Fluid Dynamics
 Department of Aeronautical Engineering, KCT

∂ ̅̅̅̅̅̅̅̅̅̅̅
∂u̅i + u′ i ̅̅̅̅̅̅̅̅
̅ + P′]
∂[P ̅̅̅̅̅̅̅̅̅̅ ̅̅̅̅̅̅̅̅̅̅
̅i + u′ i }{u
̅i + u′ i ] ∂[{u
∂[u ̅̅̅̅̅̅̅̅̅̅
̅j + u′ j }]
520 [μ ( )] − = ρ( + )
∂xj ∂xj ∂xi ∂t ∂xj
521

∂ ∂u̅i ∂ ∂u̅̅̅̅
′
i
̅] ∂[P̅′ ]
∂[P
522 [μ [ ]] + [μ [ ]] − −
∂xj ∂xj ∂xj ∂xj ∂xi ∂xi

∂u̅i ̅̅̅′̅i
∂u ∂
523 = ρ {[( ) + ( )] + [ (u u + ̅̅̅̅̅̅̅
̅̅̅̅̅ u′ i u′ j + ̅̅̅̅̅̅
ui u′ j + ̅̅̅̅̅̅
u′ i uj )]}
∂t ∂t ∂xj i j
524 From the Note-1, the following terms are zero,
̅̅̅′̅i
∂u ∂[P̅′ ] ∂u̅̅̅′̅i
525 = ̅̅̅̅̅̅
ui u′ j = ̅̅̅̅̅̅
u′ i uj = = μ[ ]=0
∂xj ∂xi ∂xj

526 Therefore,

∂ ∂u̅i ̅]
∂[P ∂u̅i ∂
527 [μ [ ]] − = ρ {[( )] + [ (u u + ̅̅̅̅̅̅̅
̅̅̅̅̅ u′ i u′ j )]}
∂xj ∂xj ∂xi ∂t ∂xj i j

528 Just rearrange the above equation into a same manner of original ‘x’ momentum

529 equation,

530

∂u̅i ∂ ∂ ∂u̅i ̅]
∂[P ∂
531 ρ {[( )] + [ (u ̅̅̅̅̅)]}
i uj = [μ [ ]] − − ̅̅̅̅̅̅̅
(ρu ′ u′ )
i j (7)
∂t ∂xj ∂xj ∂xj ∂xi ∂xj

532 The equation (7) is called as Reynolds Averaged Navier Strokes (RANS) Equation. In
∂
533 ̅̅̅̅̅̅̅
which, the process of time averaging has introduced a new term [∂x (ρu ′ u′ )]. The
i j
j

534 ̅̅̅̅̅̅̅
term (−ρu ′ u′ ) is known as Reynolds Stress or Turbulent Stress.
i j

535

536 For 3-D,

2
u′1 u′1 u′ 2 u′1 u′ 3
̅̅̅̅̅̅̅ 2
537 −ρu ′ u′ = −ρ [ ′ ′
i j u 2u 1 u′ 2 u′ 2 u′ 3 ] (8)
2
u′ 3 u′1 u′ 3 u′ 2 u′ 3
538

539 The equation (8) is contained six unknowns, which needs to be solved with the help

540 of sufficient number of equations.

541

25
Computational Fluid Dynamics
 Department of Aeronautical Engineering, KCT

542 ➢ In the RANS equations, the Reynolds Stress give additional unknowns

543 ̅̅̅̅̅̅̅
[(−ρu ′ u′ ], but there are no explicit governing differential equations for the
i j

544 additional unknowns.

545 ➢ As, Number of unknowns > Number of Equations, the problem is

546 indeterminate. One needs to close the problem to obtain a solution. This is

547 known as Closure Problem in Turbulence.

548 ➢ The Turbulence Modelling tries to represent the Reynolds stresses in terms of

549 the time averaged velocity components.

550 ➢ The common turbulence models are classified on the basis of the number of

551 additional transport equations that need to be solved along with the RANS

552 equations.

553 ➢ Several approaches have evolved to model Reynolds stress tensor. The

554 commonly followed methodologies include

555 o Eddy Viscosity Models, and

556 o Reynolds Stress Transport Models

557

Number of Extra Name of the Model

Equations

Zero Mixing Length

One Spalart-Allmaras

Two k-epsilon

Two k-omega

Three Transition-k-kl-omega

Four Transition SST

Five Reynolds Stress

26
Computational Fluid Dynamics
 Department of Aeronautical Engineering, KCT

558

559 ̅̅̅̅̅̅̅
The term Reynolds Stress or Turbulent Stress (−ρu ′ u′ ) are further expanded for
i j

560 simplification,

2ρ ∂[ui ] ∂[uj ]
561 ̅̅̅̅̅̅̅
−ρu ′ u′ = −
i j kδij + μt ( + ) (9)
3 ∂xj ∂xi
∂[u ] ∂[uj ]
562 Boussinesq Eddy Viscosity Approximation [μt ( ∂xi + )] is incorporated in the
j ∂xi

563 above equation (9),

564 Where, k is Kinetic Energy of Turbulent Fluctuations, and μt is Eddy Viscosity

2 2 2
u′1 + u′ 2 + u′ 3
565 k= (10)
2
566 Substitute the Equation (9) in Equation (7)
∂u̅i ∂
567 ρ {[( )] + [ (u ̅̅̅̅̅)]}
u
∂t ∂xj i j

∂ ∂u̅i ̅]
∂[P ∂ 2ρ ∂[ui ] ∂[uj ]
568 = [μ [ ]] − + (− kδij + μt ( + )) (11)
∂xj ∂xj ∂xi ∂xj 3 ∂xj ∂xi

569

∂u̅i ∂ ∂ ∂u̅i ̅effective ] ∗ δij

∂[P
570 ρ {[( )] + [ (u ̅̅̅̅̅)]}
i uj = [μeffective [ ]] − (12)
∂t ∂xj ∂xj ∂xj ∂xi

571

572 Where,
2ρ
573 ̅
Peffective = P + k (13)
3
574

575 μeffective = μ + μt (14)

576

577 All the other terms can be estimated through average velocities except Eddy Viscosity

578 [μt ]. So we have to further model the term Eddy Viscosity [μt ].

579

580 Under k-Epsilon turbulence model

581
27
Computational Fluid Dynamics
 Department of Aeronautical Engineering, KCT

582 The standard values of all the model constants as fitted with benchmark experiments

583 [Launder and Sharma, Letters in Heat and Mass Transfer, 1 (1974), (131 - 138)],

584

k2
585 μt = ρCμ (15)
ε
586

587 Where, Cμ = 0.09, ′k′ is kinetic energy of turbulence, ′ε′ is turbulence dissipation rate

588

589 The Model equation for the Turbulent Kinetic Energy [k] is as follows:

590

𝐷𝑘 𝜕𝑘 𝜕𝑘 𝜕 𝑣𝑡 𝜕𝑘
591 = + u̅j = [ ]+𝑃−𝜀 (16)
𝐷𝑡 𝜕𝑡 ∂xj ∂xj 𝜎𝑘 ∂xj
592 Where,
𝜕𝑘 𝜕𝑘
593 = 𝑅𝑎𝑡𝑒 𝑜𝑓 𝐼𝑛𝑐𝑟𝑒𝑎𝑠𝑒 ′𝑘 ′ ; u̅j = 𝐶𝑜𝑛𝑣𝑒𝑐𝑡𝑖𝑣𝑒 𝑇𝑟𝑎𝑛𝑠𝑝𝑜𝑟𝑡;
𝜕𝑡 ∂xj
𝜕 𝑣𝑡 𝜕𝑘
594 [ ] = 𝐷𝑖𝑓𝑓𝑢𝑠𝑖𝑣𝑒 𝑇𝑟𝑎𝑛𝑠𝑝𝑜𝑟𝑡; 𝑃 = 𝑅𝑎𝑡𝑒 𝑜𝑓 𝑃𝑟𝑜𝑑𝑢𝑐𝑡𝑖𝑜𝑛;
∂xj 𝜎𝑘 ∂xj
595 𝜀 = 𝑅𝑎𝑡𝑒 𝑜𝑓𝑑𝑒𝑠𝑡𝑟𝑢𝑐𝑡𝑖𝑜𝑛
596

597 The Model equation for the Turbulent Dissipation [𝜀] is as follows:
𝐷𝜀 𝜕𝜀 𝜕𝜀 𝜕 𝑣𝑡 𝜕𝜀 𝑃𝜀 𝜀2
598 = + u̅j = [ ] + 𝐶𝑔1 − 𝐶𝑔2 (17)
𝐷𝑡 𝜕𝑡 ∂xj ∂xj 𝜎𝑔 ∂xj 𝑘 𝑘
599 Where,
𝜕𝜀 𝜕𝜀
600 = 𝑅𝑎𝑡𝑒 𝑜𝑓 𝐼𝑛𝑐𝑟𝑒𝑎𝑠𝑒 ′𝜀 ′ ; u̅j = 𝐶𝑜𝑛𝑣𝑒𝑐𝑡𝑖𝑣𝑒 𝑇𝑟𝑎𝑛𝑠𝑝𝑜𝑟𝑡;
𝜕𝑡 ∂xj
𝜕 𝑣𝑡 𝜕𝜀 𝑃𝜀
601 [ ] = 𝐷𝑖𝑓𝑓𝑢𝑠𝑖𝑣𝑒 𝑇𝑟𝑎𝑛𝑠𝑝𝑜𝑟𝑡; 𝐶𝑔1 = 𝑅𝑎𝑡𝑒 𝑜𝑓 𝑃𝑟𝑜𝑑𝑢𝑐𝑡𝑖𝑜𝑛;
∂xj 𝜎𝑔 ∂xj 𝑘
𝜀2
602 𝐶𝑔2 = 𝑅𝑎𝑡𝑒 𝑜𝑓𝑑𝑒𝑠𝑡𝑟𝑢𝑐𝑡𝑖𝑜𝑛
𝑘
603 Where, 𝐶𝑔1 = 1.44; 𝐶𝑔2 = 1.92; 𝜎𝑘 = 1; 𝜎𝑔 = 1.3 ;

604

605 Under k-Omega turbulence model

28
Computational Fluid Dynamics
 Department of Aeronautical Engineering, KCT

606

607 This model comes under two equation category, in which a modified version of the

608 ‘k’ equation used in the ‘k-ω’ model. A transport equation for ‘ω’ is dissipation per

609 unit kinetic energy.

k
610 μt = ρ (18)
ω
611

612 The equation for ‘k’ and ‘ω’ are as follows:

613

𝐷𝑘 𝜕𝑘 𝜕𝑘 𝜕 𝑣𝑡 𝜕𝑘 ∂[ui ] ∂[uj ] ∂[ui ]

614 = + u̅j = [𝑣 + ] + 𝑣t ( + ) − 𝛽 ∗ 𝑘𝜔 (19)
𝐷𝑡 𝜕𝑡 ∂xj ∂xj 𝜎𝑘 ∂xj ∂xj ∂xi ∂xj
615

𝐷𝜔 𝜕𝜔 𝜕𝜔 𝜕 𝑣𝑡 𝜕𝜔 ∂[ui ] ∂[uj ] ∂[ui ]

616 = + u̅j = [𝑣 + ] + 𝛼𝑣t ( + ) − 𝛽1 𝜔2 (20)
𝐷𝑡 𝜕𝑡 ∂xj ∂xj 𝜎∞ ∂xj ∂xj ∂xi ∂xj
617

618 The standard values of all the model constants are,

619

620 𝛽1 = 0.075, 𝛽 ∗ = 0.09, 𝑣t = 0.553, 𝜎𝑘 = 𝜎∞ = 2.0:

621

622 References

623

624 1. [Link]
625

29
Computational Fluid Dynamics