Section 5: Finite Volume Methods for the Navier Stokes Equations
In this initial lecture we introduce the two-dimensional Navier-Stokes Equations
continuity equation:
𝜕𝜌 𝜕(𝜌𝑢) 𝜕(𝜌𝑣)
+ + =0
𝜕𝑡 𝜕𝑥 𝜕𝑦
x-momentum:
𝜕𝑢 𝜕𝑢 𝜕𝑢 𝜕𝑝 𝜕 𝜕𝑢 𝜕 𝜕𝑢 𝜕𝑣
𝜌( +𝑢 +𝑣 )=− + (2𝜇 + 𝜆∇ ∙ 𝑢̅) + [𝜇 ( + )]
𝜕𝑡 𝜕𝑥 𝜕𝑦 𝜕𝑥 𝜕𝑥 𝜕𝑥 𝜕𝑦 𝜕𝑦 𝜕𝑥
y-momentum:
𝜕𝑣 𝜕𝑣 𝜕𝑣 𝜕𝑝 𝜕 𝜕𝑢 𝜕𝑣 𝜕 𝜕𝑣
𝜌( +𝑢 +𝑣 )=− + [𝜇 ( + )] + (2𝜇 + 𝜆∇ ∙ 𝑢̅)
𝜕𝑡 𝜕𝑥 𝜕𝑦 𝜕𝑦 𝜕𝑥 𝜕𝑦 𝜕𝑥 𝜕𝑦 𝜕𝑦
The above equations are expressed in what is known as the non-conservative
form. For implementation of the finite-volume method, the conservation-law
form is appropriate. The equivalent equations in conservation law form appear as
follows:
x-momentum
𝜕(𝜌𝑢) 𝜕(𝜌𝑢𝑢) 𝜕(𝜌𝑣𝑢) 𝜕𝑝 𝜕 𝜕𝑢 𝜕 𝜕𝑢 𝜕𝑣
+ + =− + (2𝜇 + 𝜆∇ ∙ 𝑢̅) + [𝜇 ( + )]
𝜕𝑡 𝜕𝑥 𝜕𝑦 𝜕𝑥 𝜕𝑥 𝜕𝑥 𝜕𝑦 𝜕𝑦 𝜕𝑥
y-momentum
𝜕(𝜌𝑣) 𝜕(𝜌𝑢𝑣) 𝜕(𝜌𝑣𝑣) 𝜕𝑝 𝜕 𝜕𝑢 𝜕𝑣 𝜕 𝜕𝑣
+ + =− + [𝜇 ( + )] + (2𝜇 + 𝜆∇ ∙ 𝑢̅)
𝜕𝑡 𝜕𝑥 𝜕𝑦 𝜕𝑦 𝜕𝑥 𝜕𝑦 𝜕𝑥 𝜕𝑦 𝜕𝑦
How? Let’s look at the unsteady and convection terms of the x-momentum
equation:
𝜕(𝜌𝑢) 𝜕(𝜌𝑢𝑢) 𝜕(𝜌𝑣𝑢)
+ + ≡
𝜕𝑡 𝜕𝑥 𝜕𝑦
𝜕𝑢 𝜕𝜌 𝜕(𝜌𝑢) 𝜕𝑢 𝜕𝑢 𝜕(𝜌𝑣)
𝜌 +𝑢 +𝑢 + 𝜌𝑢 + 𝜌𝑣 +𝑢
𝜕𝑡 𝜕𝑡 𝜕𝑥 𝜕𝑥 𝜕𝑦 𝜕𝑦
1
Note the terms highlighted in red represent the continuity equation multiplied by
u, hence sum to zero. Hence the two forms of expressing the left side of the
momentum equations are equivalent.
The momentum equations may also be rearranged on the right side as:
𝜕𝑢 𝜕𝑢 𝜕𝑢
𝜌( +𝑢 +𝑣 )=
𝜕𝑡 𝜕𝑥 𝜕𝑦
𝜕𝑝 𝜕 𝜕𝑢 𝜕 𝜕𝑢 𝜕 𝜕𝑢 𝜕 𝜕𝑣 𝜕
− + (𝜇 ) + (𝜇 ) + [ (𝜇 ) + (𝜇 )] + (𝜆∇ ∙ 𝑢̅)
𝜕𝑥 𝜕𝑥 𝜕𝑥 𝜕𝑦 𝜕𝑦 𝜕𝑥 𝜕𝑥 𝜕𝑦 𝜕𝑥 𝜕𝑥
𝜕𝑣 𝜕𝑣 𝜕𝑣
𝜌( +𝑢 +𝑣 )=
𝜕𝑡 𝜕𝑥 𝜕𝑦
𝜕𝑝 𝜕 𝜕𝑣 𝜕 𝜕𝑣 𝜕 𝜕𝑣 𝜕 𝜕𝑢 𝜕
− + (𝜇 ) + (𝜇 ) + [ (𝜇 ) + (𝜇 )] + (𝜆∇ ∙ 𝑢̅)
𝜕𝑦 𝜕𝑥 𝜕𝑥 𝜕𝑦 𝜕𝑦 𝜕𝑦 𝜕𝑦 𝜕𝑥 𝜕𝑦 𝜕𝑦
For each equation, the viscous terms (in the square brackets) sum to zero for a
constant property, incompressible flow, as does the last term in the equations.
How? Let’s look at the terms from the x-momentum equation:
𝜕 𝜕𝑢 𝜕 𝜕𝑣 𝜕 𝜕𝑢 𝜕 𝜕𝑣 𝜕 𝜕𝑢 𝜕𝑣
[ (𝜇 ) + (𝜇 )] ≡ (𝜇 ) + (𝜇 ) ≡ 𝜇 ( + )
𝜕𝑥 𝜕𝑥 𝜕𝑦 𝜕𝑥 𝜕𝑥 𝜕𝑥 𝜕𝑥 𝜕𝑦 𝜕𝑥 𝜕𝑥 𝜕𝑦
Finally, denoting the bracketed viscous terms and the 2nd viscosity coefficient
terms as Sx and Sy respectively, the momentum equations may be written in
conservation-law form as:
x-momentum
𝜕(𝜌𝑢) 𝜕𝑝
+ ∇ ∙ (𝜌𝑢̅𝑢) = − + ∇ ∙ (𝜇∇𝑢) + 𝑆𝑥
𝜕𝑡 𝜕𝑥
y-momentum
𝜕(𝜌𝑣) 𝜕𝑝
+ ∇ ∙ (𝜌𝑢̅𝑣) = − + ∇ ∙ (𝜇∇𝑣) + 𝑆𝑦
𝜕𝑡 𝜕𝑦
For a constant property, incompressible flow the source terms are zero.
2
Finite Volume Method
To implement the finite volume method, we integrate this form of the governing
equations over a control volume. By applying the Gauss divergence theorem we
convert the convective and viscous volume integral terms to surface integrals.
For the x-momentum equation:
∫[∇ ∙ (𝜌𝑢̅𝑢)] 𝑑𝑉 = ∫ 𝜌𝑢(𝑢̅ ∙ 𝑛̂) 𝑑𝐴
and
∫[∇ ∙ (𝜇∇𝑢)]𝑑𝑉 = ∫(𝜇∇𝑢) ∙ 𝑛̂𝑑𝐴
The pressure gradient term is evaluated as a volume integral.
The x-momentum equation becomes:
𝜕𝑝
∫ 𝜌𝑢(𝑢̅ ∙ 𝑛̂) 𝑑𝐴 = − ∫ 𝑑𝑉 + ∫(𝜇∇𝑢) ∙ 𝑛̂𝑑𝐴
𝜕𝑥
Similarly for the y-momentum equation:
𝜕𝑝
∫ 𝜌𝑣(𝑢̅ ∙ 𝑛̂) 𝑑𝐴 = − ∫ 𝑑𝑉 + ∫(𝜇∇𝑣) ∙ 𝑛̂𝑑𝐴
𝜕𝑦
In the following lectures we will look at the evaluation of these integrals over a
control volume. However, our first step, discussed in the next lecture, is to look
at the placement of the unknowns on the mesh.
3
4
5
The Semi-Implicit Method for Pressure-Linked Equation (SIMPLE)
Procedure
1) Guess values for u, v, p.
2) Solve the momentum equations for u and v. The resulting velocities
will satisfy momentum, but not mass conservation.
3) Solve a pressure correction equation for a pressure correction p’.
4) Update the guessed pressure using the just-computed pressure
correction, p’.
5) Update the velocities resulting from the solution to the momentum
equations with velocity corrections, u’ and v’ (which are directly related
to the computed pressure correction, p’). These new velocities will
now satisfy mass conservation, but no longer satisfy momentum.
6) Go back to step 1 where the guessed values are now the results from
steps 4 (for pressure) and 5 (for velocities).
The process continues until convergence to a pressure distribution that,
when used in the momentum equations, results in both mass
conserving and momentum conserving velocities out of the momentum
solvers.
In the following lectures we will derive:
1) the discretized momentum equations
2) expressions for velocity corrections in terms of pressure corrections
3) the discretized continuity equation, which leads to the pressure
correction equation
6
Discretization of x-momentum equation:
𝜕𝑝
∫ 𝜌𝑢(𝑢̅ ∙ 𝑛̂) 𝑑𝐴 = − ∫ 𝑑𝑉 + ∫(𝜇∇𝑢) ∙ 𝑛̂𝑑𝐴
𝜕𝑥
𝜌𝑢𝑢|𝑒𝑤 Δ𝑦 + 𝜌𝑣𝑢|𝑛𝑠 Δ𝑥 =
𝜕𝑢 𝑒 𝜕𝑢 𝑛
(𝑃𝑤 − 𝑃𝑒 )Δ𝑦 + 𝜇 | ∆𝑦 + 𝜇 | ∆𝑥
𝜕𝑥 𝑤 𝜕𝑦 𝑠
Define:
𝑚̇𝑒 = (𝜌𝑢)𝑒 Δ𝑦 𝑚̇𝑤 = (𝜌𝑢)𝑤 Δ𝑦
𝑚̇𝑛 = (𝜌𝑣 )𝑛 Δ𝑥 𝑚̇𝑠 = (𝜌𝑣 )𝑠 Δ𝑥
Rewrite the discretized equation as:
𝑚̇𝑒 𝑢𝑒 − 𝑚̇𝑤 𝑢𝑤 + 𝑚̇𝑛 𝑢𝑛 − 𝑚̇𝑠 𝑢𝑠 =
𝑈𝐸 − 𝑈𝑃 𝑈𝑃 − 𝑈𝑊
(𝑃𝑤 − 𝑃𝑒 )Δ𝑦 + 𝜇𝑒 Δ𝑦 − 𝜇𝑤 Δ𝑦 +
(𝛿𝑥 )𝑒 (𝛿𝑥 )𝑤
𝑈𝑁 − 𝑈𝑃 𝑈𝑃 − 𝑈𝑆
𝜇𝑛 Δ𝑥 − 𝜇𝑠 Δ𝑥
(𝛿𝑦)𝑛 (𝛿𝑦)𝑠
Group the coefficients of 𝑈𝑃 :
𝑚̇𝑒 𝑢𝑒 − 𝑚̇𝑤 𝑢𝑤 + 𝑚̇𝑛 𝑢𝑛 − 𝑚̇𝑠 𝑢𝑠 =
𝜇𝑒 Δ𝑦 𝜇𝑤 Δ𝑦 𝜇𝑛 Δ𝑥 𝜇𝑠 Δ𝑥
(𝑃𝑤 − 𝑃𝑒 )Δ𝑦 − [ + + + ]𝑈 +
(𝛿𝑥 )𝑒 (𝛿𝑥 )𝑤 (𝛿𝑦)𝑛 (𝛿𝑦)𝑠 𝑃
𝜇𝑒 Δ𝑦 𝜇𝑤 Δ𝑦 𝜇𝑛 Δ𝑥 𝜇𝑠 Δ𝑥
𝑈 + 𝑈 + 𝑈 + 𝑈
(𝛿𝑥 )𝑒 𝐸 (𝛿𝑥 )𝑤 𝑊 (𝛿𝑦)𝑛 𝑁 (𝛿𝑦)𝑠 𝑆
7
Need to interpolate to find face values for 𝑈𝑒 , 𝑈𝑤 , 𝑈𝑛 , 𝑈𝑠
We will use 1st order upwinding hence flow direction must be
taken into account. We use the max function for that:
𝑚𝑎𝑥[𝑚̇𝑒 , 0]𝑈𝑃 − 𝑚𝑎𝑥 [−𝑚 ̇ 𝑒 , 0]𝑈𝐸 +
𝑚𝑎𝑥[−𝑚̇𝑤 , 0]𝑈𝑃 − 𝑚𝑎𝑥[𝑚̇𝑤 , 0]𝑈𝑊 +
𝑚𝑎𝑥[𝑚̇𝑛 , 0]𝑈𝑃 − 𝑚𝑎𝑥[−𝑚 ̇ 𝑛 , 0]𝑈𝑁 +
𝑚𝑎𝑥[−𝑚̇𝑠 , 0]𝑈𝑃 − 𝑚𝑎𝑥[𝑚̇𝑠 , 0]𝑈𝑆 =
𝜇𝑒 Δ𝑦 𝜇𝑤 Δ𝑦 𝜇𝑛 Δ𝑥 𝜇𝑠 Δ𝑥
(𝑃𝑤 − 𝑃𝑒 )Δ𝑦 − [ + + + ]𝑈 +
(𝛿𝑥 )𝑒 (𝛿𝑥 )𝑤 (𝛿𝑦)𝑛 (𝛿𝑦)𝑠 𝑃
𝜇𝑒 Δ𝑦 𝜇𝑤 Δ𝑦 𝜇𝑛 Δ𝑥 𝜇𝑠 Δ𝑥
𝑈 + 𝑈 + 𝑈 + 𝑈
(𝛿𝑥 )𝑒 𝐸 (𝛿𝑥 )𝑤 𝑊 (𝛿𝑦)𝑛 𝑁 (𝛿𝑦)𝑠 𝑆
Recall we need to ensure that
𝐴𝑃 = 𝐴𝐸 + 𝐴𝑊 + 𝐴𝑁 + 𝐴𝑆 + (𝑚̇𝑒 − 𝑚̇𝑤 + 𝑚̇𝑛 − 𝑚̇𝑠 )
It currently doesn’t. Note the difference in the signs within the
max functions.
8
We can rewrite the terms multiplying 𝑈𝑃 as:
𝑚𝑎𝑥[𝑚̇𝑒 , 0] + 𝑚𝑎𝑥 [−𝑚̇𝑤 , 0] +
{ } 𝑈𝑃 ≡
𝑚𝑎𝑥[𝑚̇𝑛 , 0] + 𝑚𝑎𝑥[−𝑚̇𝑠 , 0]
𝑚𝑎𝑥[−𝑚̇𝑒 , 0] + 𝑚𝑎𝑥[𝑚̇𝑤 , 0] +
{ 𝑚𝑎𝑥[−𝑚̇𝑛 , 0] + 𝑚𝑎𝑥[𝑚̇𝑠 , 0] + } 𝑈𝑃
(𝑚̇𝑒 − 𝑚̇𝑤 + 𝑚̇𝑛 − 𝑚̇𝑠 )
To check this, let’s assume a 1D flow west-to-east so that 𝑚̇𝑒
and 𝑚̇𝑤 are both positive:
𝑚̇𝑒 𝑈𝑃 = 𝑚̇𝑤 𝑈𝑃 + (𝑚̇𝑒 − 𝑚̇𝑤 )𝑈𝑃
or
𝑚̇𝑒 𝑈𝑃 = 𝑚̇𝑒 𝑈𝑃
If we assume flow from east-to-west so that 𝑚̇𝑒 and 𝑚̇𝑤 are
both negative:
𝑚̇𝑤 𝑈𝑃 = 𝑚̇𝑒 𝑈𝑃 + (−𝑚̇𝑒 + 𝑚̇𝑤 )𝑈𝑃
or
𝑚̇𝑤 𝑈𝑃 = 𝑚̇𝑤 𝑈𝑃
9
Now, we can define our coefficients as:
𝜇𝑒 Δ𝑦
𝐴𝑢𝐸 = 𝑚𝑎𝑥[−𝑚̇𝑒 , 0] +
(𝛿𝑥 )𝑒
𝜇𝑤 Δ𝑦
𝐴𝑢𝑊 = 𝑚𝑎𝑥[𝑚̇𝑤 , 0] +
(𝛿𝑥 )𝑤
𝜇𝑛 Δ𝑥
𝐴𝑢𝑁 = 𝑚𝑎𝑥[−𝑚̇𝑛 , 0] +
(𝛿𝑦)𝑛
𝜇𝑠 Δ𝑥
𝐴𝑆𝑢 = 𝑚𝑎𝑥[𝑚̇𝑠 , 0] +
(𝛿𝑦)𝑠
So that:
𝐴𝑢𝑃 = 𝐴𝑢𝐸 + 𝐴𝑢𝑊 + 𝐴𝑢𝑁 + 𝐴𝑆𝑢 + (𝑚̇𝑒 − 𝑚̇𝑤 + 𝑚̇𝑛 − 𝑚̇𝑠 )
The discretized x-momentum equation is then written as:
𝐴𝑢𝑃 𝑈𝑃 = 𝐴𝑢𝐸 𝑈𝐸 + 𝐴𝑢𝑊 𝑈𝑊 + 𝐴𝑢𝑁 𝑈𝑁 + 𝐴𝑆𝑢 𝑈𝑆 + (𝑃𝑤 − 𝑃𝑒 )Δ𝑦
(Note that if you sufficiently converge your pressure corrections
(see pg. 13-15) such that the corrected velocities conserve mass
over each cell, then the (𝑚̇𝑒 − 𝑚̇𝑤 + 𝑚̇𝑛 − 𝑚̇𝑠 ) term in the
definition of 𝐴𝑢𝑃 will be very close to zero. In that case, I have
found that the code will converge even if those terms are
neglected.)
10
Solution Procedure
Next, we will build an under-relaxation procedure into the
iterative solver, since the equation is nonlinear.
1 𝑢 𝑢 𝑢 𝑢
𝑈𝑃 = 𝑢 (𝐴𝐸 𝑈𝐸 + 𝐴𝑊 𝑈𝑊 + 𝐴𝑁 𝑈𝑁 + 𝐴𝑆 𝑈𝑆 + (𝑃𝑤 − 𝑃𝑒 )Δ𝑦)
𝐴𝑃
Add and subtract 𝑈𝑃𝑜𝑙𝑑 from the right side:
𝑈𝑃 = 𝑈𝑃𝑜𝑙𝑑 +
Ω 𝑢 𝑢 𝑢 𝑢
𝑢 [(𝐴𝐸 𝑈𝐸 + 𝐴𝑊 𝑈𝑊 + 𝐴𝑁 𝑈𝑁 + 𝐴𝑆 𝑈𝑆 + (𝑃𝑤 − 𝑃𝑒 )Δ𝑦) −
𝐴𝑃
𝐴𝑢𝑃 𝑈𝑃𝑜𝑙𝑑 ]
The above is an SOR type scheme. Let’s keep working on it.
Multiply through by 𝐴𝑢𝑃 ⁄Ω.
𝐴𝑢𝑃 𝐴𝑢𝑃 𝑜𝑙𝑑
𝑈 = 𝑈 +
Ω 𝑃 Ω 𝑃
[𝐴𝑢𝐸 𝑈𝐸 + 𝐴𝑢𝑊 𝑈𝑊 + 𝐴𝑢𝑁 𝑈𝑁 + 𝐴𝑆𝑢 𝑈𝑆 + (𝑃𝑤 − 𝑃𝑒 )Δ𝑦] − 𝐴𝑢𝑃 𝑈𝑃𝑜𝑙𝑑
or, rearranging:
𝐴𝑢𝑃 𝐴𝑢𝑃 𝑜𝑙𝑑
𝑈 = (1 − Ω) 𝑈𝑃 +
Ω 𝑃 Ω
[𝐴𝐸 𝑈𝐸 + 𝐴𝑢𝑊 𝑈𝑊 + 𝐴𝑢𝑁 𝑈𝑁 + 𝐴𝑆𝑢 𝑈𝑆 + (𝑃𝑤 − 𝑃𝑒 )Δ𝑦]
𝑢
Now, redefine ̃𝑢 = 𝐴𝑢𝑃 ⁄Ω so that:
𝐴 𝑃
𝑈𝑃 = (1 − Ω)𝑈𝑃𝑜𝑙𝑑 +
1 𝑢
[𝐴𝐸 𝑈𝐸 + 𝐴𝑢𝑊 𝑈𝑊 + 𝐴𝑢𝑁 𝑈𝑁 + 𝐴𝑆𝑢 𝑈𝑆 + (𝑃𝑤 − 𝑃𝑒 )Δ𝑦]
̃
𝐴𝑢
𝑃
11
where, to under relax the solution, Ω < 1.
Following the same procedure, the discretized y-momentum
equation becomes (where 𝐴̃𝑣 = 𝐴𝑣𝑃 ⁄Ω):
𝑃
𝑉𝑃 = (1 − Ω)𝑉𝑃𝑜𝑙𝑑 +
1 𝑣 𝑣 𝑣 𝑣
[ 𝐴 𝐸 𝑉𝐸 + 𝐴 𝑊 𝑉𝑊 + 𝐴 𝑁 𝑉𝑁 + 𝐴𝑆 𝑉𝑆 + (𝑃𝑠 − 𝑃𝑛 )Δ𝑥]
𝐴̃𝑣
𝑃
where:
𝜇𝑒 Δ𝑦
𝐴𝑣𝐸 = 𝑚𝑎𝑥[−𝑚̇𝑒 , 0] +
(𝛿𝑥 )𝑒
𝜇𝑤 Δ𝑦
𝐴𝑣𝑊 = 𝑚𝑎𝑥[𝑚̇𝑤 , 0] +
(𝛿𝑥 )𝑤
𝜇𝑛 Δ𝑥
𝐴𝑣𝑁 = 𝑚𝑎𝑥[−𝑚̇𝑛 , 0] +
(𝛿𝑦)𝑛
𝜇𝑠 Δ𝑥
𝐴𝑆𝑣 = 𝑚𝑎𝑥[𝑚̇𝑠 , 0] +
(𝛿𝑦)𝑠
So that:
𝐴𝑣𝑃 = 𝐴𝑣𝐸 + 𝐴𝑣𝑊 + 𝐴𝑣𝑁 + 𝐴𝑆𝑣 + (𝑚̇𝑒 − 𝑚̇𝑤 + 𝑚̇𝑛 − 𝑚̇𝑠 )
12
Pressure Correction Procedure
𝐴𝑢𝑃𝑤 𝑈𝑤 = ∑ 𝐴𝑢𝑛𝑏 𝑈𝑛𝑏 + (𝑃𝑊 − 𝑃𝑃 )∆𝑦
𝑛𝑏
𝐴𝑢𝑃𝑒 𝑈𝑒 = ∑ 𝐴𝑢𝑛𝑏 𝑈𝑛𝑏 + (𝑃𝑃 − 𝑃𝐸 )∆𝑦
𝑛𝑏
𝐴𝑣𝑃𝑠 𝑉𝑠 = ∑ 𝐴𝑣𝑛𝑏 𝑉𝑛𝑏 + (𝑃𝑆 − 𝑃𝑃 )∆𝑥
𝑛𝑏
𝐴𝑣𝑃𝑛 𝑉𝑛 = ∑ 𝐴𝑣𝑛𝑏 𝑉𝑛𝑏 + (𝑃𝑃 − 𝑃𝑁 )∆𝑥
𝑛𝑏
13
Let * denote a guessed value, then write:
𝐴𝑢𝑃𝑤 𝑈𝑤∗ = ∑ 𝐴𝑢𝑛𝑏 𝑈𝑛𝑏
∗ ∗
+ (𝑃𝑊 − 𝑃𝑃∗ )∆𝑦
𝑛𝑏
𝐴𝑢𝑃𝑒 𝑈𝑒∗ = ∑ 𝐴𝑢𝑛𝑏 𝑈𝑛𝑏
∗
+ (𝑃𝑃∗ − 𝑃𝐸∗ )∆𝑦
𝑛𝑏
𝐴𝑣𝑃𝑛 𝑉𝑛∗ = ∑ 𝐴𝑣𝑛𝑏 𝑉𝑛𝑏
∗
+ (𝑃𝑃∗ − 𝑃𝑁∗ )∆𝑥
𝑛𝑏
𝐴𝑣𝑃𝑠 𝑉𝑠∗ = ∑ 𝐴𝑣𝑛𝑏 𝑉𝑛𝑏
∗
+ (𝑃𝑆∗ − 𝑃𝑃∗ )∆𝑥
𝑛𝑏
Now define:
𝑈 = 𝑈∗ + 𝑈′
𝑉 = 𝑉∗ + 𝑉′
𝑃 = 𝑃∗ + 𝑃′
(i.e., U= guessed value + correction)
Subtract the second set of equations from the first:
𝐴𝑢𝑃𝑤 𝑈𝑤′ = ∑ 𝐴𝑢𝑛𝑏 𝑈𝑛𝑏
′ ′
+ (𝑃𝑊 − 𝑃𝑃′ )∆𝑦
𝑛𝑏
etc.
This is an equation for the corrections to velocity and pressure.
14
If we neglect the summation term:
∆𝑦 ′
𝑈𝑤′ = 𝑢 (𝑃𝑊 − 𝑃𝑃′ )
𝐴̃ 𝑃𝑤
∆𝑦 ′
𝑈𝑒′ = 𝑢 (𝑃𝑃 − 𝑃𝐸′ )
𝐴̃𝑃𝑒
∆𝑥 ′
𝑉𝑠′ = (𝑃𝑆 − 𝑃𝑃′ )
𝐴̃𝑣
𝑃𝑠
∆𝑥 ′
𝑉𝑛′ = 𝑣 (𝑃𝑃 − 𝑃𝑁′ )
𝐴̃𝑃𝑛
These are known as the velocity correction equations.
Note that all the 𝐴𝑢𝑃𝑤 , etc. for the corrections need to be the
tilde values in the computer code, i.e., 𝐴̃ 𝑢
𝑃𝑤 , our iteration
equation use the tilde value of the Ap’s.
15
Next, let’s look at the (steady) continuity equation:
𝜕(𝜌𝑢) 𝜕(𝜌𝑣)
+ =0
𝜕𝑥 𝜕𝑦
Discretize using the finite volume approach:
[(𝜌𝑈)𝑒 − (𝜌𝑈)𝑤 ]Δ𝑦 + [(𝜌𝑉 )𝑛 − (𝜌𝑉 )𝑠 ]Δ𝑥 = 0
Substitute:
𝑈 = 𝑈∗ + 𝑈′
𝑉 = 𝑉∗ + 𝑉′
𝑃 = 𝑃∗ + 𝑃′
(𝜌𝑈𝑒∗ + 𝜌𝑈𝑒′ − 𝜌𝑈𝑤∗ − 𝜌𝑈𝑤′ )Δ𝑦 +
(𝜌𝑉𝑛∗ + 𝜌𝑉𝑛′ − 𝜌𝑉𝑠∗ − 𝜌𝑉𝑠′ )Δ𝑥 = 0
Gather the “starred” and “prime” terms:
(𝜌𝑈𝑒′ − 𝜌𝑈𝑤′ )Δ𝑦 + (𝜌𝑉𝑛′ − 𝜌𝑉𝑠′ )Δ𝑥 +
[(𝜌𝑈𝑒∗ − 𝜌𝑈𝑤∗ )Δ𝑦 + (𝜌𝑉𝑛∗ − 𝜌𝑉𝑠∗ )Δ𝑥 ] = 0
The portion in yellow represents a “mass source term,” S.
Now substitute velocity corrections:
𝜌Δ𝑦 2 ′ ′
𝜌Δ𝑦 2 ′ ′
𝜌Δ𝑥 2 ′ ′)
(𝑃𝑃 − 𝑃𝐸 ) − ( 𝑃𝑊 − 𝑃𝑃 ) + ( 𝑃𝑃 − 𝑃𝑁 −
̃
𝐴𝑢
𝐴̃
𝑢
𝐴̃
𝑣
𝑃𝑒 𝑃𝑤 𝑃𝑛
𝜌Δ𝑥 2 ′
(𝑃𝑆 − 𝑃𝑃′ ) + 𝑆 = 0
𝐴̃
𝑣
𝑃𝑠
Gather the 𝑃𝑃′ terms on the left side, the remainder on the
right:
16
𝜌Δ𝑦 2 𝜌Δ𝑦 2 𝜌Δ𝑥 2 𝜌Δ𝑥 2 ′
( 𝑢 + 𝑢 + 𝑣 + 𝑣 ) 𝑃𝑃 =
𝐴̃ 𝐴̃
𝑃𝑒 𝐴̃
𝑃𝑤 𝐴̃ 𝑃𝑛 𝑃𝑠
2 2
𝜌Δ𝑦 ′ 𝜌Δ𝑦 ′ 𝜌Δ𝑥 2 ′ 𝜌Δ𝑥 2 ′
𝑃𝐸 + 𝑢 𝑃𝑊 + 𝑣 𝑃𝑁 + 𝑣 𝑃𝑆 − 𝑆 = 0
̃
𝐴𝑢
𝐴̃ 𝐴̃ 𝐴̃
𝑃𝑒 𝑃𝑤 𝑃𝑛 𝑃𝑠
or:
𝐴𝑃 𝑃𝑃′ = 𝐴𝐸 𝑃𝐸′ + 𝐴𝑊 𝑃𝑊
′
+ 𝐴𝑁 𝑃𝑁′ + 𝐴𝑆 𝑃𝑆′ − 𝑆
Now, putting into the form of an SOR scheme:
𝜔
𝑃𝑃′ = 𝑃𝑃′ + (𝐴𝐸 𝑃𝐸′ + 𝐴𝑊 𝑃𝑊
′
+ 𝐴𝑁 𝑃𝑁′ + 𝐴𝑆 𝑃𝑆′ − 𝑆 − 𝐴𝑃 𝑃𝑃′ )
𝐴𝑃
This is our equation over which we will iterate to find the
pressure corrections.
17
What about boundary conditions?
We do not change velocities on boundaries using velocity
correction terms. Hence, velocities normal to the boundaries
must be set so that mass is conserved globally.
That is, on the east/west boundaries:
𝑈 ∗ = 𝑈 and 𝑈 ′ = 0
This means that, for instance, on the east boundary:
Δ𝑦 ′ ′
𝑈𝑒′ = 𝑒 (𝑃𝑃 − 𝑃𝐸 ) ≡ 0
𝐴𝑃𝑢
Consequently, 𝑃𝑃′ = 𝑃𝐸′
We can force this in our code by setting the coefficient 𝐴𝐸 to
zero along the entire east boundary, or equivalently by setting
𝐴𝑒𝑃𝑢 to a very large number (𝑖. 𝑒. , 1030 ). Similarly for the other
boundaries.
Note that the pressure update must be under relaxed:
𝑃𝑛𝑒𝑤 = 𝑃𝑜𝑙𝑑 + 𝛼𝑃′ where 0 < 𝛼 < 1
Velocities are also corrected using the velocity correction
equation. These corrections are not under relaxed as they
conserve mass. The corrections look like:
Δ𝑦 ′
𝑈𝑒𝑛𝑒𝑤 = 𝑈𝑒𝑜𝑙𝑑 + 𝑈𝑒′ = 𝑈𝑒𝑜𝑙𝑑 + (𝑃𝑃 − 𝑃𝐸′ )
𝐴̃𝑢
𝑃𝑒
Note the pressure correction equation itself is linear, so it can
be over relaxed, that is, 1 < 𝜔 < 2
18
Suggested Index Notation
(and what is used in Fortran code)
19
The Driven Cavity Problem
0<=x<=1
0<=y<=1
Let density = 1 kg/m^3; dynamic viscosity = 0.01 kg/m s.
North side moves left-to-right at 1 m/s.
Results in a Reynolds number of 100.
20
21
Developing Flow in a Channel
0<=x<=10
0<=y<=1
Let density = 1 kg/m^3; dynamic viscosity = 0.01 kg/m s.
Inlet velocity =1 m/s.
Zero derivative velocity specification at outlet.
Reynolds number = 100.
22
Other Boundary Conditions
Pressure Boundary Condition
1) Set Ap in pressure correction equation to 10^30 on pressure
boundary condition cells. This will force 𝑃′ to remain at 0,
hence pressure will not be altered from initial specification.
2) Then, say on east end of domain, compute 𝑈𝑒 by enforcing
continuity:
−(𝜌𝑉𝐴)𝑛 + (𝜌𝑉𝐴)𝑠 + (𝜌𝑈𝐴)𝑤
𝑈𝑒 =
(𝜌𝐴)𝑒
23
Symmetry Boundary Condition
i.e., Symmetry Plane on North Boundary
1) Set V=0 along north face of top wall cells (as a boundary
condition).
2) Set du/dy = 0 along north face of top wall cells.
3) In pressure correction routine, set the 𝐴𝑁 coefficient along
top wall cells to 0.
24
Procedure to Generate Fully Developed Flow
Velocity at inlet plane updated to velocity at outlet plane as the
solver proceeds through the outer iterations.
Useful for creating a fully developed flow profile for input into
another simulation as inlet boundary condition.
1) Set inlet and outlet velocity profiles from initial guess so that
mass is conserved overall.
2) Set a zero streamwise derivative boundary condition at the
outlet
3) Complete an outer iteration, then impose the outlet velocity
profile onto the inlet profile.
4) Continue process until flow is fully developed.
See example channel flow next page:
25
26
Place Solid Object within Flow
Block out middle cell.
1) In the x-momentum solver set 𝐴𝑢𝑃 to a very large number
(10^30) so that 𝑈𝑃 does not change from initial guess (which
should be zero for a blocked out cell).
2) Same procedure for y-momentum to fix 𝑉𝑃 at zero.
Can construct fairly elaborate geometries by this basic
approach.
Example of Flow Blockage within Channel Flow (left to right).
Note distance downstream outflow boundary is from blockage.
27
