代写MCEN90008 Fluid Dynamics

THE 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