MONASH UNIVERSITY
Department of Mechanical and Aerospace Engineering
MEC4447 — Computers in Fluids and Energy —
Assignment 3
Due Monday 5pm 24th May submitted through Moodle. It contributes 10% to
the overall subject mark.
Solving the Navier-Stokes Equations in Streamfunction-Vorticity Form
The steady incompressible Navier-Stokes and continuity equations are
u · ∇u = −∇P + ν∇2u,
and
∇ · u = 0.
Note, that there is no explicit equation for the pressure. This creates some difficulty in
solving the system numerically. This will be discussed further in lectures.
For two-dimensional problems, it is possible to circumvent this difficulty by writing the
equations in terms of the streamfunction and the vorticity. Because using the streamfunction
automatically satisfies the continuity equation, this reduces the equation set down to two.
The pressure does not occur in the vorticity equation.
The governing equations are
u · ∇ω = ν∇2ω,
and
∇2Ψ = −ω.
Here, u = (u, v), the two-dimensional velocity field, ω is the z component of the vorticity (the
only non-zero component, since the flow is two-dimensional) and Ψ is the streamfunction.
The vorticity is given by
ω =
∂v
∂x
− ∂u
∂y
,
and the relationship between the velocity and streamfunction is
u =
∂Ψ
∂y
v = −∂Ψ
∂x
.
1
Hence, eliminating the velocity components, allows the governing equations to be written as
∂Ψ
∂y
∂ω
∂x
− ∂Ψ
∂x
∂ω
∂y
− ν(∂
2ω
∂x2
+
∂2ω
∂y2
) = 0,
∂2Ψ
∂x2
+
∂2Ψ
∂y2
+ ω = 0.
Note, that the first equation is non-linear and it becomes difficult to solve as the Reynolds
number increases. The second equation is linear — just Poisson’s equation — and is relatively
easy to solve numerically.
Both equations can be discretised using central difference approximations for both first and
second derivatives.
Boundary Conditions
At solid walls a streamline has to follow the wall, unless there is flow through the wall. The
vorticity is not zero. A relationship between the streamfunction and the vorticity can be
obtained by expanding out the streamfunction at a point next to a wall.
Expanding the streamfunction at the point p about the wall (boundary) point b using Taylor-
Series gives
Ψp = Ψb −∆y∂Ψ
∂y
+
(∆y)2
2
∂2Ψ
∂y2
− (∆y)
3
6
∂3Ψ
∂y3
+O((∆y)4).
∂Ψ/∂y is the tangential velocity of the wall, i.e., ub. In addition, at the wall, ω = −∂u/∂y
(the other component is zero). Hence,
(Ψp −Ψb) + ∆y ub = (∆y)
2
2
∂u
∂y
− (∆y)
3
6
∂
∂y
(
∂u
∂y
) +O((∆y)4),
(Ψp −Ψb) + ∆y ub = −(∆y)
2
2
ωb +
(∆y)3
6
∂ωb
∂y
+O((∆y)4),
(Ψp −Ψb) + ∆y ub = −(∆y)
2
2
ωb +
(∆y)3
6
(
ωb − ωp
∆y
+O(∆y)) +O((∆y)4),
(Ψp −Ψb) + ∆y ub = −(∆y)
2
3
(ωb +
1
2
ωp) +O((∆y)
4),
ωb = − 3
∆y2
((Ψp −Ψb) + ∆y ub)− 1
2
ωp +O(∆y
2).
and by symmetry, for the right hand boundary
ωb = − 3
∆x2
((Ψp −Ψb)−∆x vb)− 1
2
ωp +O(∆x
2).
The expressions for the bottom and left boundaries can be obtained by reversing the sign of
the increment (i.e., ∆x→ −∆x or ∆y → −∆y).
2
bstreamfunction zero on wall
velocity given here too
py∆
Figure 1: Evaluation of boundary conditions for vorticity.
Generally, the value of the streamfunction at the wall is known, so this equation relates the
vorticity at the wall to the value of the streamfunction and vorticity just off the wall.
Implementation
The two governing equations can be expressed in finite-difference form for each internal point.
Generally, the streamfunction is known at all points on the boundary, while the vorticity
depends on the boundary velocity and streamfunction as above.
The equations can be solved by relaxation, i.e., successive over-relaxation for the stream-
function equation and successive under-relaxation for the vorticity equation.
Program
You can download the skeleton source from Moodle. You need to add some lines as indicated
in the program to set the equations, boundary conditions, etc.
Note that the routine is supplied as a matlab function rather than a script file (because
this makes the code run much faster). You can call the routine by name (typing cavity in
the command window - as long as the file is in your matlab path). You should supply the
arguments required.
These arguments are:
• the number of cells in x and y: nx, ny;
• the Reynolds number: re;
3
• the Pecle´t number: pe;
• the relaxation parameters for the ω, Ψ and temperature equations: relax1, relax2,
relax3;
• the convergence criteria: eps psi, eps omega, eps temp;
To call the function, supply (name, value) pairs as the arguments.
For example, typing
cavity(’re’, 100, ’pe’, 100, ’nx’, 32, ’ny’, 32, ’relax1’, 1.0, ’relax2’, 1.5, ’relax3’, 1.0);
calls the function with the values: re=100, pe=100, nx=32, ny=32, relax1=1.0, relax2=1.5,relax3=1.0.
Other parameters not specified are taken as default values defined at the top of the function.
The defaults are: eps psi=0.000001, eps omega=0.000001, eps temp=0.000001, relax1=1.0,
relax2=1.5, relax3=1.0 nx=32, ny=32, re=100, pe=100.
The Project
You are required to obtain results to an accuracy of 1% or better.
Part 1 (6 marks):
Use the streamfunction vorticity formulation to find the solution for the driven cavity flow for
Reynolds numbers: Re = 100, 1000, 3200. The Reynolds number is Re = Utop×(cavity width)
ν
.
You can use a series of grids to estimate the grid-independent value of the minimum stream-
function for each case.
Part 2 (4 marks):
Extend the code to also calculate the temperature field, which satisfies the equation
u · ∇T = κ∇2T,
where T is the temperature, κ is the thermal diffusivity (conduction coefficient) and u the
velocity vector. Take the (scaled) temperature over the top of the cavity as Tt = 0 and
the temperature at the bottom as Tb = 1. Assume the side walls are adiabatic, i.e., there
is no heat exchange to/from the walls, so that dT/dn ≡ n · ∇T = (dT/dx) = 0. Take
the ratio of viscosity to thermal diffusivity to be unity, i.e. Re = Pe (the Reynolds and
Pecle´t numbers are the same). Note that the flow field does not depend on the temperature
- so that the latter can be evaluated after the velocity field is calculated. (It is also OK to
calculate the temperature in the same loop as the streamfunction, of course, as is done in
the skeleton code.) The Pecle´t number is Pe = Utop×(cavity width)
κ
. Note that to implement
normal derivative boundary conditions you can use a second-order one-sided finite-difference
approximation:
dT
dx
|b = (−3Tb + 4Tp − Tpp)
2∆x
+O(∆x2),
4
with b, p, pp, indicating the boundary, and first and second internal points respectively.
Project Writeup
I want a short report containing only the following information.
• Finite-difference equations used + boundary conditions (You can hand write these).
• Final program.
• Two tables listing grid sizes, relaxation parameters, iterations taken, minimum stream-
function + temperature at the same point for each Reynolds number. Include extrapo-
lated values of the minimum streamfunction using Richardson extrapolation. Describe
how you obtained the grid-independent values.
• Streamfunction, vorticity and temperature contour plots corresponding Re = 100 and
Re = 3200.
1 unit
u = 1
Streamfunction = 0 on entire boundary
Velocity zero on boundary, except at top lid.
Sliding lid sets up internal recirculation
Driven Cavity Flow
1 unit
Figure 2: Problem definitions.
You can find lots of references on this in the library.
The standard reference for the driven cavity is
5
• Ghia U, Ghia KN, Shin CT, High-Re Solutions For Incompressible-Flow Using The
Navier Stokes Equations And A Multigrid Method, Journal Of Computational Physics
48 (3): 387-411 1982.
Also, you can find a discussion of boundary conditions in Computational Fluid Dynamics,
by P Roache.
Obtain the skeleton code from Moodle.
6

欢迎咨询51作业君