# 代写MCEN90008 Fluid Dynamics

- 首页 >> WebTHE UNIVERSITY OF MELBOURNE

DEPARTMENT OF MECHANICAL ENGINEERING

MCEN90008 Fluid Dynamics 2023

PART I OF ASSIGNMENT FOR POTENTIAL FLOW

Instructions:

• Assignment to be handed in by 23:59 on Sunday 17th September 2023.

• This assignment should be done in groups of 2 students. Both students

in the group will get the same mark for this assignment. If you choose

to do the assignment alone, no concession will be given. Your assignment

will be marked the same as an assignment done by two students.

• Please choose your group partner carefully.

• Only hand in one assignment per group.

• You can use any programming language to complete this assignment.

Should you decide to write your computer program in MATLAB, you

are not allowed to use the streamline, odexx or similar functions in MATLAB i.e. you need to write the program to solve the ordinary differential

equations yourself and not simply use the functions in MATLAB. Also do

not merely use the contour function of ψ to visualise streamlines. The

aim of this assignment is for you to produce a set of tools to enable you

compute the coordinates of pathlines. The same rules apply if you use

other languages.

• Please submit your assignment as aa zipped package containing

– Copies of all computer programs (appended as .m files if using Matlab).

– Written documentation with computer generated graphs and sketches.

This documentation must contain all discussions and all the exercises

that you were asked to do in the main part of this assignment.

• Marks will be deducted for incorrect or absent axis labels.

• Unless stated otherwise, please use equal scaling along each axis (this is

achieved by setting daspect([1 1 1]) In Matlab).

1

Part I

For the first part of this assignment, you will write codes to accurately plot pathlines from given starting

positions. For steady flows, pathlines are the same as streamlines. Make sure you get this question right.

For the later parts of this assignment you will build on the computer program which you have written

here to study much more complex flows.

1. (Total marks for question = 27) In order to write a computer code to sketch streamlines, it is

best to first start by writing a computer program to solve a simple problem.

(a) Write down the analytical solution to the following ordinary differential equation (ODE) in

the form x = f(t)

dx

dt = −

x

3

4

(1)

with the initial condition

x = 8 at t = 0. (2)

(1 Mark)

(b) Write a computer program to solve Eq. (1) using Euler’s (EU) and the 4th order Runge Kutta

(RK-4) method (see the appendix for more information on EU and RK-4). Solve the equation

from 0 ≤ t ≤ tf where tf = 0.2, 0.4, 0.6, 0.8 and 1.0. Perform your computations with various

step sizes, (h = 0.1, 0.01, 0.001), and tabulate your results in the table shown below (show 4

decimal places).

EU RK-4

tf Analytical answer h = 0.1 h = 0.01 h = 0.001 h = 0.1 h = 0.01 h = 0.001

0.2

0.4

0.6

0.8

1.0

The output from your computer program must approach the analytical answer as h gets

smaller. Plot the analytical solution (x vs t) along with the EU and the RK4 solution for

h = 0.01 and 0 < t < 1. Use the axis limits axis([0 1 0 8]) and daspect([1 10 1]).

Which is the more accurate numerical method? EU or RK-4? (8 Marks)

(c) Write down the analytical solution (write them as equations of the form x = f(t) and y = f(t))

to the following set of coupled ODE

dx

dt = x + y (3)

dy

dt = −2x − y

You are given that at time, t = 0

x

y

=

1

2

. (4)

(3 Marks)

2

(d) Extend your program from part (1b) so that it can solve the set of equations given by the set

of Equations (3) for 0 ≤ t ≤ tf . Check your program by completing the table below for both

x and y for various tf and time-steps h. To save time in future questions, you should try to

write your program such that you can easily change the set of equations f(x, y) and g(x, y) as

defined in the Appendix.

EU RK-4

tf Analytical answer h = 0.1 h = 0.01 h = 0.001 h = 0.1 h = 0.01 h = 0.001

2

4

6

8

10

(8 Marks)

Which is the more accurate numerical method? Use the more accurate numerical method in

ALL subsequent parts of this assignment.

(e) Plot y vs x for h = 0.1, h = 0.01 and h = 0.001 for both EU and RK-4 with tf = 2π. Plot

the EU and RK4 calculations on separate figures. Should the lines join up? Do the lines join

up? (4 Marks)

(f) Use your computer program to plot the streamline pattern for the flow described by Eq. (3)

i.e. plot y vs x for a set of initial conditions, where tf = 10 and h = 0.01. You should modify

your program such that you can pass it a number of starting positions xs and ys, the time step

h and total time of integration tf . In matlab, this could be achieved by writing your program

as a function where you pass the variables xs and ys (which could be vectors) and tf and h.

This function would then be used as follows,

[xr yr] = my_streamline_function(xs, ys, tf, h)

where xr and yr are the returned coordinates (matrices) of the streamline and xs and ys are

the start locations,

xs = [1 1.5 2];

ys = [2 3 4];

For non-Matlabers, your program could read an input file where the first line contains the

number of streamlines, Nl

, the time of integration tf and the time step h. The subsequent

Nl

lines in the input file should contain the location of the initial points of the streamlines,

(x0i y0i). For example, the input file to draw three solution trajectories that start from (1, 2),

(1.5, 3) and (2, 4) with h = 0.01 and 0 < t < 10 will look something like,

3 10 0.01

1.0 2

1.5 3

2.0 4

(3 Marks)

2. (Total marks for question = 19) This question is based on an exam question from the 2022

potential flow exam. We wish to model a Pitot-static tube as shown in figure 1(a) using a source

and a uniform flow as shown in figure 1(b). The source is offset in the positive x direction from

the origin by distance s such that the stagnation point is located at the origin (x, y) = (0, 0). The

total thickness of the Pitot-static tube in the limit far downstream of the leading edge (as x → ∞)

is D = 0.01 m. The Pitot-static tube is inserted into a flow with U∞ = 50 ms−1

.

(a) Starting from the arrangement of flow features shown in figure 1(b), derive expressions and

obtain values for the source strength Q and the offset s. (4 marks)

3

U∞

p∞ p0

xtap

ps

(a) D = 0.01m

U∞

p∞

(b)

s

y

x

Figure 1

(b) Use the codes that you developed in question 2 to produce a plot showing the streamline

pattern associated with your source and sink. Please show the streamlines external of the

Pitot-static tube in black, and the internal streamlines (from the source) as blue, and lines of

separatrix in red. Use the following axis limits for your plot axis([-0.08 0.1 -0.02 0.02]),

and ensure that you use a 1:1 aspect ratio (daspect([1 1 1])). Hint: because your Pitotstatic tube is small and the freestream velocity is high, you will need to use a very small time

step (h ≈ 1×10−5

s) to accurately obtain streamlines (even with the RK4 method).(6 marks)

(c) The Pitot static tube measures velocity using Bernoulli’s equation applied along the stagnation

streamline shown in blue in figure 1(a). Show that,

U∞ =

s

2(p0 − p∞)

ρ

However, the Pitot-tube does not directly measure p∞, instead it probes the pressure ps from

a static pressure tap some distance xtap downstream of the leading edge stagnation point.

Ups =

s

2(p0 − ps)

ρ

If the velocity measured by the Pitot-static tube (Ups) is 2% higher than U∞, show that the

pressure coefficient Cpps at the location of the static tap must be,

Cpps =

ps − p∞

1

2

ρU2∞

≈ −0.04

Hence plot filled colour contours of Cp around the Pitot-static tube body (the half Rankine

body) and determine the minimum distance xtap such that the Pitot-static tube measured

4

velocity Ups is only 2%. greater than the true oncoming freestream velocity U∞. Do not show

Cp within the probe body. Use the same axis limits suggested for part (b) axis([-0.08 0.1

-0.02 0.02]), and use colour axis limits caxis([-0.5 0.5]) for Cp. (6 marks).

(d) For the assumed value of Cpps

show that if we assume xtap is large relative to D, we can

approximate that,

xtap ≈

Q

2πU∞(

p

1 − Cp − 1)

(3 marks)

Part II

3. (Total marks for question = 20) The streamfunction for a potential vortex (with positive

anti-clockwise circulation) is given by

ψ(r) = −

Γ

2π

ln r = −

Γ

2π

ln p

(x − x0)

2 + (y − y0)

2 (5)

where Γ is the circulation and r is the distance from the centre of the vortex. The center coordinates

are given by (x0,y0).

(a) Show that the time it takes for a fluid particle to circulate around this vortex, tc, is given by

tc =

4π

2

r

2

Γ

(6)

(2 Marks)

The Cartesian components of velocity, u and v for this potential vortex are given by

u(x, y) = −

Γ

2π

y − y0

(x − x0)

2 + (y − y0)

2

(7)

and

v(x, y) = Γ

2π

x − x0

(x − x0)

2 + (y − y0)

2

. (8)

Note that u(x, y) and v(x, y) are singular at (x, y) = (x0, y0). Analytically, this is not usually an

issue. However, when you are writing a computer program to solve problems or track particles,

functions that have singularities will usually cause your program to “blow up”.

(b) In order to regularise u and v, it is proposed to use the following streamfunction,

ψδ(r) = −

Γ

2π

ln p

r

2 + δ

2 (9)

for the vortex instead of Eq. (5).

Show that the Cartesian components of velocity, u and v can be written as

uδ(x, y) = −

Γ

2π

y − y0

(x − x0)

2 + (y − y0)

2 + δ

2

(10)

and

vδ(x, y) = Γ

2π

x − x0

(x − x0)

2 + (y − y0)

2 + δ

2

. (11)

Comment on the assumption of irrotational flow for a vortex as defined in equation (9) shift

text (3 Marks)

5

(c) Four vortices of the form given in equation (9) with δ = 0.2, are arranged in an unbounded

space with the following properties to model a mixing flow due to four fixed stirrers.

Vortex x0 y0 Γ δ

1 -1.5 0 -1 0.2

2 1.5 0 -1 0.5

3 0 -1.5 1 0.5

4 0 1.5 1 0.5

Use the quiver command in Matlab to plot the velocity vector field associated with this

arrangement. (use a vector spacing of 0.2 in the x and y directions). (3 Marks)

(d) On Canvas you will find a Matlab data file called:

particle_positions_2023.mat

This file contains the x and y locations of 10050 particles of fluid at time t = 0 within the

mixing system described above. You may read in the data using the command:

load(‘particle_positions_2023.mat’);

You can view the particle locations using:

plot(xp,yp, ‘k.’)

Use your RK4 codes to track these particles over 1001 steps of h = 0.1. Plot the final particle

locations and submit this in your report. The stirrers (vortices) are fixed in physical space and

do not mutually induce motion in each other. You only need to track the particle trajectories

under the influence of the 4 stirrers. (6 Marks)

(e) Comment on how realistic such a scenario might be? (2 Marks)

(f) Again using the particle starting locations given in ‘particle_positions.mat’, this time use

your Euler codes to track these particles over 1001 steps of h = 0.1. Plot the final particle

locations and submit this in your report. Comment on the comparison to the RK4 result.

shift text (4 Marks)

4. (Total marks for question = 36).

The half thickness yt as a function of x/c (where 0 ≤ x/c ≤ 1) for the NACA0018 profile is given

by the following expression,

yt = 0.9

0.2969r

x

c

− 0.1260 x

c

− 0.3516 x

c

2

+ 0.2843 x

c

3

− 0.1015 x

c

4

(a) The 0018 designation tells us that the total maximum thickness to chord ratio of this profile is

18%. Assuming that we are transforming a circle with radius a = 1m into a Joukowski airfoil,

determine the required horizontal offset x0 and the Joukowski circle radius b that will transform

the shifted circle to an airfoil with the same thickness to chord ratio (it is likely that you will

need to seek a numerical solution to this question). Hint: In the Joukowsky transformed zj

plane, the chord is along the real axis, and the thickness is defined in the imaginary direction.

The Joukowski transform is given below where zsc are the complex coordinates of the shifted

circle.

zj = (zsc) + b

2

(zsc)

Plot the Joukowski and NACA0018 profile on the same figure. Comment on similarities /

differences. (2 Marks)

Next we will use this value of b to step our way through the process of transforming the flow over

a shifted cylinder with rotation to the flow over an airfoil.

First consider a negative doublet in uniform flow from left to right with circulation, for which the

complex potential function is,

w = U∞

z +

a

2

z

| {z }

cylinder flow

−

a

2

z

iΓ

2π

log(z)

| {z }

potential vortex

6

Note that a positive value for Γ will result in counter clockwise rotation.

(b) If a = 1, U∞ = 100 and Γ = −3πaU∞, plot the streamlines for this flow. Also plot the surface

of the cylinder, stagnation points and any separatrix lines. It will be much easier for this

question if you work in the complex plane. (4 Marks)

We are going to move towards a Joukowski transformation of this cylinder flow. We

eventually wish to get the flow for a symmetrical airfoil at an angle-of-attack of α = 8◦

.

This involves 4 transformations.

Counter-clockwise rotation (to get angle of attack over cylinder) z1 = zeiΘ

Bodily shift z2 = z1 − x0 + iy0

Jowkowski transformation z3 = z2 +

b

2

z2

Clockwise rotation (to ensure horizontal freestream) z4 = z3e

−iΘ

where Θ is a positive value. We also need to adjust the circulation strength Γ to ensure that

the Kutta condition is satisfied (see parts d and e). We will perform these 4 transformations

over parts (c) and (e).

(c) Transformations 1 & 2. Rotate and then bodily shift the original flow plotted in part (a)

using the following transformations (to obtain a horizontally shifted cylinder at angle of attack

α)

z1 = zeiΘ z2 = z1 + x0

Assume an angle of attack α = 8◦

. Use the bodily shift x0 such that you will obtain a sharp

trailing edge after application of the Joukowski transformation. Plot the results in the z2 plane

(including cylinder, streamlines, separatrix and the Joukowski circle). (4 Marks)

(d) Transformations 3 & 4. Now apply the Joukowski transform,

z3 = z2 +

b

2

z2

and then rotate the flow to ensure that we have horizontal freestream velocity in the far

upstream of the airfoil.

z4 = z3e

−iΘ

(where Θ is the original angle applied in part c). Plot the results in the z4 plane, including the

transformed cylinder (which should now be a Jowkowski airfoil), streamlines and separatrix.

Comment on the flow at the trailing edge. What is wrong with it? dummy text (4 Marks)

(e) Adjust the circulation Γ. The Kutta-Joukowski condition states that,

A body with a sharp trailing edge which is moving through a fluid will create about

itself a circulation of sufficient strength to hold the rear stagnation point at the trailing

edge.

Derive an expression to alter the circulation Γ in the z plane such that this condition is met in

the z4 plane. Replot the results (including cylinder, streamlines, separatrix and the Jowkowski

circle) in the z4 plane with the new determined value of Γ. (2 Marks)

The pressure distribution around the Jowkowski airfoil can be written as,

Cp = 1 −

V

2

J

U2∞

where VJ is the magnitude of the velocity on the surface of the Jowkowski airfoil. The velocity on

the surface of the airfoil can be found from the following,

VJ = Vc

dz

dz4

Where Vc is the magnitude of the velocity on the surface of the cylinder in the z plane.

7

(f) Plot a new figure that is the same as the previous, but now includes colour contours of the

velocity magnitude field around your airfoil (V4 =

p

u

2

4 + v

2

4

). You can use the Matlab commands pcolor or contourf to plot these contours.fill fill fill fill fill fill fill fill fill fill fill fill fill

fill fill fill fill fill (4 Marks)

(g) Plot Cp vs x/c (where x is distance along the chord line from leading to trailing edge) for the

Joukowski airfoil at α = 8◦

. Plot also on the same figure (with a different colour line) the Cp

vs x/c obtained experimentally from your airfoil laboratory at the same angle of attack (if you

do not have lab data at α = 8◦

, you may compare your lab data and Jowkowski result at a

different matched α, as close as possible to 8◦

). Comment on differences and similarities. fill

fill fill fill filll (6 Marks)

(h) Finally, compute lift coefficients CL via integration of your obtained Cp distribution for the

Joukowski airfoil. Perform this computation for a range of angles of attack α from 0◦

to 20◦

in

steps of 2◦

. Compare this result with the experimental data and with the analytical solution

L = ρΓU∞. Plot all three curves on the same figure using different colours. Compare and

discuss the results from the three curves. (8 Marks)

8

Appendix: Numerical solution to ordinary differential equations

Consider the following ordinary differential equation,

d

dtx(t) = f(x). (12)

If f(x) is a relatively simple function, then one can obtain exact solution to Eq. (12) using analytical

techniques which you should have learnt in your maths course. However, when f(x) is a complicated

function, one needs to use a computer and numerical algorithms to approximate the solution to Eq.

(12). There are many numerical algorithms that can be used to approximate the solution to an ordinary differential equation. In this assignment, you are asked to explore the use of two popular methods,

Euler’s method and the 4th order Runge-Kutta method (for more information about these methods, see

references [1], [2] and [3]). The formula for both methods are

Euler’s method:

xn+1 = xn + hf(xn)

4th order Runge-Kutta:

xn+1 = xn +

1

6

k1 +

1

3

(k2 + k3) + 1

6

k4

h

where

k1 = f(xn)

k2 = f

xn +

h

2

k1

k3 = f

xn +

h

2

k2

k4 = f (xn + hk3)

where xn is the approximate value of x(t = tn) where tn = nh. h is the time step that you choose for the

simulation. Generally, the smaller the value of h, the numerical solution that you obtain is more accurate

and stable.

If you are asked to solve a system consisting of a set of coupled ordinary differential equations

d

dtx(t) = f(x, y)

d

dty(t) = g(x, y),

the approximate numerical solution can be obtained using the Euler and 4th order Runge-Kutta methods.

The relevant formulae are

Euler’s method:

xn+1 = xn + f(xn, yn)h

yn+1 = yn + g(xn, yn)h

9

4th order Runge-Kutta method:

xn+1 = xn +

1

6

k11 +

1

3

(k21 + k31) + 1

6

k41

h

yn+1 = yn +

1

6

k12 +

1

3

(k22 + k32) + 1

6

k42

h

where the approximated solution of x(t = tn) and y(t = tn) are denoted as xn and yn respectively. kij is

calculated as follows

k11 = f(xn, yn)

k12 = g(xn, yn)

k21 = f

xn +

h

2

k11, yn +

h

2

k12

k22 = g

xn +

h

2

k11, yn +

h

2

k12

k31 = f

xn +

h

2

k21, yn +

h

2

k22

k32 = g

xn +

h

2

k21, yn +

h

2

k22

and

k41 = f (xn + hk31, yn + hk32)

k42 = g (xn + hk31, yn + hk32)

10

References

[1] Atkinson, K., Elementary Numerical Analysis, John Wylie & Sons, 1985.

[2] Chapra, S. C. and Canale, R. P., Numerical Methods for Engineers, McGraw-Hill, 2002.

[3] Kreyszig, E., Advanced Engineering Mathematics, John Wylie & Sons, 1993.

11