MONASH UNIVERSITY

Department of Mechanical and Aerospace Engineering

MEC4447 — Computers in Fluids and Energy —

Assignment 3

Due Monday 8th June 5 pm through Moodle. (Do it much earlier if possible...)

It contributes 10% to the overall subject mark.

Submit modified Matlab code and writeup.

Solving the Navier-Stokes Equations in Streamfunction-Vorticity Form

to Predict Vortex Breakdown in a Cylindrical Container

The unsteady incompressible Navier-Stokes and continuity equations are

∂u

∂t

+ 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 some problems, it is possible to circumvent this difficulty by writing the equations in

terms of the streamfunction and the vorticity. Because using the streamfunction automati-

cally satisfies the continuity equation by design, this reduces the equation set. The pressure

does not occur in the vorticity equation.

An interesting case is for flows in axisymmetric geometries using cylindrical polar coordinates

when the velocity components do not depend on the azimuthal angle. In this case the

governing equations reduce to

∂ωθ

∂t

+ uz

∂ωθ

∂z

+ ur

∂ωθ

∂r

= ν(

∂2ωθ

∂z2

+

1

r

∂

∂r

(r

∂ωθ

∂r

)−

ωθ

r2

) +

2uθ

r

∂uθ

∂z

+

urωθ

r

, (1)

∂uθ

∂t

+ uz

∂uθ

∂z

+ ur

∂uθ

∂r

+

uruθ

r

= ν(

∂2uθ

∂z2

+

1

r

∂

∂r

(r

∂uθ

∂r

)−

uθ

r2

), (2)

and

−

∂

∂z

(

1

r

∂Ψ

∂z

)−

∂

∂r

(

1

r

∂Ψ

∂r

) = −

1

r

∂2Ψ

∂z2

−

1

r

∂2Ψ

∂r2

+

1

r2

∂Ψ

∂r

= ωθ. (3)

Here, u = (uz, ur, uθ) is the three-dimensional velocity field, with each component only a

function of the axial and radial coordinates z and r, respectively. Also, ωθ is the θ component

of the vorticity and Ψ is the axisymmetric streamfunction, defined below.

1

The vorticity components in cylindrical polar coordinates, and noting there is no variation

with θ, are given by

ωz =

1

r

∂(ruθ)

∂r

−

1

r

∂ur

∂θ

=

1

r

∂(ruθ)

∂r

, ωr =

1

r

∂uz

∂θ

−

∂uθ

∂z

= −

∂uθ

∂z

, ωθ =

∂ur

∂z

−

∂uz

∂r

,

and the relationships between the axial and radial velocity components and streamfunction

are

uz =

1

r

∂Ψ

∂r

ur = −

1

r

∂Ψ

∂z

. (4)

Note that the continuity equation in cylindrical polar coordinates is

∂uz

∂z

+

1

r

∂(rur)

∂r

+

1

r

∂uθ

∂θ

= 0.

Substituting the definitions of ur and uz in terms of the streamfunction into this expression

shows that it is satisfied identically.

So, the first two equations, (1) and( 2) are (Parabolic) time-stepping equations, whilst the

last equation, (3), is an elliptic equation. Each equation can be discretised using central

difference approximations for both first and second spatial derivatives.

For this project we are going to predict the flow inside a cylindrical container driven by a

spinning lid rotating with a constant angular velocity, as shown in the figure below. The

spinning lid centrifuges the fluid near the lid towards the outer radial wall. It then moves

left along the wall to the other end, then back towards the axis. At this stage it travels back

up towards the spinning lid as a strongly vortical (spinning) flow. Under certain conditions,

this vortex can burst or break down, meaning that it expands radially very quickly. This

situation can occur in aeronautical flows with the trailing vortices that occur from wing and

winglet tips. The original F/A-18 jets owned by the Australian Defence Force had a problem

that at high angles of attack the trailing vortex coming off the leading edge extension would

burst over the wings, leading to huge unsteady pressure forces.

rotating lid

H

R

Ω

solution plane

rotation axis

r

z

Figure 1: Left: Spinning lid rig. Right: Vortex breakdown over an F18 simulated in a water

channel.

2

Boundary Conditions

At solid walls a streamline has to follow the wall, unless there is flow through the wall. In

addition, the same streamline must also pass along the axis due to cylindrical symmetry.

Thus, at all boundaries the streamfunction can be taken as Ψ = 0. The azimuthal vorticity

is not zero - except on the axis. In general, the 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 w using

Taylor-series gives

Ψp = Ψw +∆z

∂Ψ

∂z

+

(∆z)2

2

∂2Ψ

∂z2

+

(∆z)3

6

∂3Ψ

∂z3

+O((∆z)4).

Here ∂Ψ/∂z = −rur,w, where ur,w is equal to tangential velocity at the wall. In addition,

at the wall, ωθ = ∂ur/∂z (the other component is zero at the wall). Note that terms up to

third-order are required to maintain second-order spatial accuracy, since the vorticity comes

from the second-order term in this expansion, as is seen below. Hence,

(Ψp −Ψw) + ∆z ur,w ≃ −

(∆z)2

2

∂(r ur,w)

∂z

−

(∆z)3

6

∂

∂z

∂(r uz,w)

∂z

,

(Ψp −Ψw) + ∆z ur,w = −

(∆z)2r

2

∂ur,w

∂z

−

(∆r)3r

6

∂2ur,w

∂z2

,

(Ψp −Ψw) + ∆z ur,w = −

(∆z)2r ωθ,w

2

−

(∆z)3r

6

∂ωθ,w

∂z

,

(Ψp −Ψw) + ∆z ur,w = −

(∆z)2r ωθ,w

2

−

(∆z)3r

6

(ωθ,p − ωθ,w)

∆z

,

(Ψp −Ψw) + ∆z ur,w = −

(∆z)2r ωθ,w

3

−

(∆z)2r ωθ,p

6

.

Here, the last term, involving the time derivative of the vorticity, has been replaced by a

first-order finite-difference approximation. This maintains the second-order overall accuracy

because the error in this last term is then O(∆z3 × ∆z) = O(∆z4), consistent with terms

that we are omitting. (To get the expression for the vorticity, we divide by ∆z2, hence the

accuracy is O(∆z4/∆z2) = O(∆z2)). Thus, we can write

ωθ,w = −3((Ψp −Ψw) + ∆r ur,w)/(∆z

2r)−

1

2

ωθ,p +O(∆z

2).

Since ur,w = 0 on the wall, and Ψw = 0, this reduces to

ωθ,w = −3

Ψp

∆z2r

−

1

2

ωθ,p +O(∆z

2).

A similar expression applies at the radial wall (with a few added terms), and at the righthand

wall.

3

wz∆

computational

boundary

p

Figure 2: Evaluation of boundary conditions for vorticity at lefthand wall.

On the central axis, the radial derivative of the axial velocity must vanish by symmetry.

Hence ∂uz/∂r = 0. We need a second-order approximation at the boundary to give overall

second-order accuracy of the solution. An appropriate approximation is

∂uz

∂r

≃

3uz,r=0 − 4uz,r=∆r + uz,r=2∆r

2∆r

= 0,

⇒ uz,r=0 = (4uz,r=∆r − uz,r=2∆r)/3.

The azimuthal velocity is zero on all boundaries except the spinning lid, where uθ = Ωr.

In addition, the azimuthal vorticity is given by ωθ =

∂ur

∂z

−

∂uz

∂r

. Given the constraint on

the axial velocity and that ur = 0 on the boundary because of symmetry, this implies the

azimuthal vorticity is zero on the axis.

Implementation

The 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 streamfunction and internal vorticity, as above.

Taking equation (1) and using a forward-time centred-space approximation gives:

ω

(n+1)

θ,i,j − ω

(n)

θ,i,j

∆t

+ u

(n)

z,i,j(

ω

(n)

θ,i+1,j − ω

(n)

θ,i−1,j

2∆z

) + ur,i,j(

ω

(n)

θ,i,j+1 − ω

(n)

θ,i,j−1

2∆r

) =

ν(

ω

(n)

θ,i+1,j − 2ω

(n)

θ,i,j + ω

(n)

θ,i−1,j

∆z2

+

ω

(n)

θ,i,j+1 − 2ω

(n)

θ,i,j + ω

(n)

θ,i,j−1

∆r2

+

1

ri,j

ω

(n)

θ,i,j+1 − ω

(n)

θ,i,j−1

2∆r

−

ω

(n)

θ,i,j

r2i,j

)+

2u

(n)

θ,i,j

ri,j

u

(n)

θ,i+1,j − u

(n)

θ,i−1,j

2∆z

+

u

(n)

r,i,jω

(n)

θ,i,j

ri,j

. (5)

4

This can be rearranged slightly to give an update equation for ωθ, to calculate the values at

the new timestep:

ω

(n+1)

θ,i,j = ω

(n)

θ,i,j +∆t

−u(n)z,i,j(ω

(n)

θ,i+1,j − ω

(n)

θ,i−1,j

2∆z

)− ur,i,j(

ω

(n)

θ,i,j+1 − ω

(n)

θ,i,j−1

2∆r

)+

ν(

ω

(n)

θ,i+1,j − 2ω

(n)

θ,i,j + ω

(n)

θ,i−1,j

∆z2

+

ω

(n)

θ,i,j+1 − 2ω

(n)

θ,i,j + ω

(n)

θ,i,j−1

∆r2

+

1

ri,j

ω

(n)

θ,i,j+1 − ω

(n)

θ,i,j−1

2∆r

−

ω

(n)

θ,i,j

r2i,j

)+

2u

(n)

θ,i,j

ri,j

u

(n)

θ,i+1,j − u

(n)

θ,i−1,j

2∆z

+

u

(n)

r,i,jω

(n)

θ,i,j

ri,j

. (6)

A similar approximation can be written for the uθ equation (2). You need to do this for this

assignment and insert your expression into the code template.

Equation (3) for the streamfunction Ψ is an elliptic equation that can be solved by relaxation.

Using second-order finite-difference approximations for the spatial derivatives gives:

−

1

ri,j

Ψi+1,j − 2Ψi,j +Ψi−1,j

∆z2

−

1

ri,j

Ψi,j+1 − 2Ψi,j +Ψi,j−1

∆r2

+

1

r2i,j

Ψi,j+1 −Ψi,j−1

2∆r

= ωθ,i,j, (7)

⇒ −

Ψi+1,j − 2Ψi,j +Ψi−1,j

∆z2

−

Ψi,j+1 − 2Ψi,j +Ψi,j−1

∆r2

+

1

ri,j

Ψi,j+1 −Ψi,j−1

2∆r

= ri,jωθ,i,j, (8)

To use successive over-relaxation (SOR), we need to obtain an expression for Ψi,j. Rear-

ranging

⇒ (

2

∆z2i,j

+

2

∆r2i,j

)Ψi,j =

Ψi+1,j +Ψi−1,j

∆z2

+

Ψi,j+1 +Ψi,j−1

∆r2

−

1

ri,j

Ψi,j+1 −Ψi,j−1

2∆r

+ri,jωθ,i,j, (9)

⇒ Ψ∗i,j =

(

Ψi+1,j +Ψi−1,j

∆z2

+

Ψi,j+1 +Ψi,j−1

∆r2

−

1

ri,j

Ψi,j+1 −Ψi,j−1

2∆r

+ ri,jωθ,i,j

)

/(

2

∆z2i,j

+

2

∆r2i,j

),

(10)

This gives a new estimate for the streamfunction at each point on the mesh given current

values. It is basically the Gauss-Seidel update equation. SOR just involves using this new

estimate (Ψ∗i,j) to give ∆Ψi,j = Ψ

∗

i,j−Ψ

old

i,j , then the new value of the streamfunction is given

by Ψnewi,j = Ψ

old

i,j + r∆Ψi,j, with r the relaxation parameter.

This defines a Gauss-Seidel update equation for Ψ, which can be used with SOR to iterate Ψ

given the current estimate of the vorticity field. Typically, this equation should be solved to

some convergence criterion at each timestep. In the code, the minimum number of iterations

is set to 5. That is, it completes at least 5 iterations to update Ψ before stepping forward

ωθ and uθ to the next timestep.

So, in summary a timestep consists of

5

• updating the Ψ field given the current vorticity field using SOR (only internal points

are updated since Ψ = 0 on the boundaries of the computational domain);

• then stepping forward the vorticity and azimuthal velocity to the next timestep using

the updated value of Ψ. Boundary conditions on ωθ (called vort in the code) are also

implemented during this update.

These two steps are repeated until the solution reaches its asymptotic state.

For the assignment you only need to fill in the update equation for uθ, which is similar in form

to equation 6, and then set other parameters such as the Reynolds number: Re = ΩR2/ν,

the physical dimensions: zlen = H and rlen = R, the relaxation parameter: relax , the

timestep: dt = ∆t, etc.

Program

You can download the template code from Moodle. You need to add some lines

defining the update equation for uθ and one of the ωθ boundary conditions. (These bits

are marked by ???? in the template). The rest should be OK. Please have a look at the

structure of the program to understand broadly how it works.

Note that the routine is supplied as a Matlab function rather than a script file (because

often this makes the code run much faster). However, make sure you are running the latest

version of MATLAB. Older versions can be very much slower...

You can call the routine by name (typing cyl cavity time in the command window - as long

as the file is in your matlab path) or just load it and run it.

You should supply the arguments required.

These arguments include:

• the number of cells in z and r: nz, nr;

• the domain size in z and r: zlen, rlen;

• the Reynolds number: re;

• the relaxation parameter for the Ψ equation: relax;

• the convergence criteria: eps psi, eps omega, eps ut;

• the timestep: dt;

To call the program as a function, supply (name, value) pairs as the arguments.

6

For example, typing

cyl cavity time(’re’, 1000, ’nz’, 60, ’nr’, 60, ’relax’, 1.5, ’dt’, 0.05);

calls the function with the values: re=100, nz=60, nr=60, relax=1.5, dt = 0.05.

Alternatively, just change the values directly in the program and rerun it.

Note that there is a primitive restart facility. After running through a case, the fields are

saved and can be used to restart the next run by changing the restart parameter to 1. If you

alter the grid size, Matlab kindly interpolates the old solution onto the new grid for you.

Note that the maximum timestep is controlled by both a Courant condition and a diffusion

condition. For fine grids you may need to use a smaller timestep to get the timestepping to

work.

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 ut=0.000001, nz=60, nr=60,

re=1000, zlen=2.5, rlen = 1.0, Ω = 1.

The Project

1. Use the streamfunction vorticity formulation to find the maximum and minimum

streamfunction for Re = 1994. Do a grid resolution study. Explain how you have

ensured that the predictions are accurate to within 2%. Supply images of the con-

verged streamfunction, azimuthal vorticity and azimuthal velocity fields. [2 marks].

2. Run the model at Re = 2765. In this case, the solution is periodic. Determine the

period of oscillation to within 2%. Again explain how you did this. (Note you should

be able to get the period to within sufficient accuracy from the convergence history

plot). [2 marks].

Project Writeup

I want a short report containing the following information.

• Derivation of the finite-difference equation for uθ, finally expressed as an update equa-

tion for uθ at the new timestep. (You can write this by hand and photograph them to

include in your writeup) [2 marks].

• Final Matlab program (submitted through Moodle) [2 marks].

• Brief answers to the questions above.

• Discretionary [2 marks] for the writeup.

7

A standard reference for this problem is

(1) Lopez, J. M., Axisymmetric vortex breakdown. Part I: Confined swirling flow. Journal

of Fluid Mechanics, 221, 533-552, 1990. (Copy on Moodle page).

8