THE UNIVERSITY OF MELBOURNE
DEPARTMENT OF MECHANICAL ENGINEERING
MCEN90008 FLUID DYNAMICS
ASSIGNMENT FOR POTENTIAL FLOW
Instructions:
• Assignment to be handed in by 23:59 on Friday 13th September 2018.
• 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 either use the C programming language or MATLAB to do this
assignment. If you should 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.
• When you submit your assignment, please include
– Copies of all computer programs.
– 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.
• For your plots please have equal scaling along each axis (this is achieved
by setting daspect([1 1 1]) In Matlab).
1
1. (20% of total mark) In order to write a computer code to sketch flow patterns and predict
代写MCEN90008作业、代写MATLAB实验作业、代写MATLAB编程语言作业
unsteady flows, it is best to first start by writing a computer program to solve a simple problem.
Make sure you get this question right. In the later part of this assignment you will build on the
computer program which you have written here to study much more complex flows.
(a) Write down the analytical solution to the following ordinary differential equation (ODE)
(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 t = 0 to t = 5. Perform your computations with various step sizes, (h = 1, 0.1, 0.01), and
tabulate your results in the table shown below
EU RK-4
t Analytical answer h = 1 h = 0.1 h = 0.01 h = 1 h = 0.1 h = 0.01
The output from your computer program must approach the analytical answer as h gets
smaller. Which is the more accurate numerical method? EU or RK-4? (5 Marks)
(c) Write down the analytical solution to the following set of coupled ODE
You are given that at time, t = 0
(d) Extend your program from part (1b) so that it can solve the set of equations given by Eq. (3)
for 0 ≤ t ≤ tf . Use tf = 2. Check your program by completing the table below for both x and
y. 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
t Analytical answer h = 0.1 h = 0.01 h = 0.001 h = 0.1 h = 0.01 h = 0.001
Which is the more accurate numerical method? Use the more accurate numerical method in
ALL subsequent parts of this assignment. (5 Marks)
(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. Should
the lines join up? Do the lines join up? (2 Marks)
(f) From your analytical solution, what is the value of tf such that when you plot y vs x, the
streamline exactly returns to the initial conditions (i.e. there is no gap or “doubling up”)?
ghost ghost (1 Marks)
(g) 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. 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, f)
where xr and yr are the returned coordinates (matrices) of the streamline. 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 two solution trajectories that start from (0, 1) and (2, 0) with h = 0.02 and
0 < t < 2 will look something like,
2. (20% of total mark) In lectures, it was demonstrated that using the method of images, the flow
of smoke from smokers (modelled as a source) in a room with a suction fan (modelled as a sink),
can be modelled by the following stream function
(5)
Here the subscript s refers to the smokers (source), and the subscript e to the extraction fan (sink).
The fan is placed le from the floor directly above the smokers (who have height ls). The fan has
strength Qe and the strength of the source is Qs.
(a) Show that the velocity components associated with the above stream function shown in Eq.
(b) Use your computer program written in Question 1 and sketch the flow pattern given by Eqs.
(6) and (7). You can assume that le = 3, ls = 1.5, Qs = 2 and Qe = 4. Again, your program
should accept an input that gives the start positions in x and y for the streamlines, the time
step and the time of integration. Assume h = 0.01
Please start your streamlines in sensible places to really show the important flow features. To
begin with I suggest 20 starting positions around the smokers at radius = 0.1, equally spaced
in the azimuthal direction. Also add other streamlines coming from the right, left and top of
your figure. Refine this with streamlines placed close to any separatrix lines and stagnation
points
For 10 radially spaced positions about the smokers
There are singularities at the center of the source and sink. You will need to add some error
trapping to catch these. One rather unsophisticated way is to prevent the velocities going too
high,
% if V > threshold, then set it to zero
V(V > threshold) = 0;
It would be better to prevent your streamlines from getting too close to the singularity. ghost
ghost (7 Marks)
(c) Derive analytical expressions for the location of all of the stagnation points in the flow. Plot
all separatrix lines onto your diagram. What is the safe distance for the non-smoker? ghost
ghost (4 Marks)
(d) Sketch the flow pattern given by Eq. 5 for Qs = 2 and Qe = 1. Again derive expressions for
the location of the stagnation points. Plot all separatrix lines onto your diagram (this might
require a modification to your program). (6 Marks)
3. (10% of total mark) After carefully designating a safe distance for the non-smoker, someone
comes along and opens the window by a tiny amount, believing that this will help. A very gentle
0.07 ms
−1 uniform breeze is imposed in the negative x direction.
(a) Write out the new expression for the stream function and the u and v velocity components.
ghost ghost (2 Marks)
(b) Assuming Qs = 2 and Qe = 4, use your code to plot the new streamline pattern. Again, marks
will be given for starting streamlines in places that allow us to best visualise the flow. ghost
ghost (4 Marks)
(c) Write out expressions for the stagnation points at ±x on the floor. Solve this equation by any
means to find the stagnation point location (these expressions may not be solvable analytically).
Compare the new safe distance for a non-smoker downwind of the smokers. (4 Marks)
4. (15% of total mark) The streamfunction for a potential vortex (with positive anti-clockwise
circulation) is given by
where Γ is the circulation and r is the distance from the centre of the vortex.
(a) Show that the Cartesian components of velocity, u and v can be written as
(b) Show that the time it takes for a fluid particle to circulate around this vortex, tc, is given by
(c) Note that u(x, y) and v(x, y) are singular at (x, y) = (0, 0). Analyticaly, this is not usually
an issue. However, when you are writing computer program to solve problems, functions that
have singularities will usually cause your program to “blow up”. In order to regularize u and
v, it is proposed to use following streamfunction
as the streamfunction for the vortex instead of Eq. (8).
i. Show that the Cartesian components of velocity, u and v can be written as
ii. Show that the distribution of vorticity, ωδ = −∇2ψδ, is given by
iii. Plot uδ(0, y) for −1 ≤ y ≤ 1 for δ = 0.2, 0.1 and 0.05 all on the same plot. Let Γ = 2π.
ghost ghost (2 Marks)
iv. Based on your analysis in this section, comment on the effects of δ on your solution. ghost
ghost (1 Mark)
(d) Write a computer program to draw instantaneous streamline patterns from point vortices.
Show that it works by using your program to plot the instantaneous streamline pattern of two
vortices of opposite signs (i.e. Γ = 2π for one vortex and Γ = −2π for another vortex) located
at (0, 1) and (0, −1). Your program should read an input file where the first line contains the
number of vortices, N. The subsequent N lines in the program should contain the location of
vortices given by x, y and the circulation of each vortex, given by Γ. For example, the input
file for the situation given above will look something like
5. (13% of total mark) The streamfunction for a potential vortex placed at a position (b, 0) can be
In subsequent parts of this assignment, it is desired to have a streamfunction which produces a
streamline that corresponds to a circle of radius a located at the origin. The streamfunction given
by Eq. (15) does not have this property.
(a) Show that in order to obtain a streamfunction that produces a streamline which correspond
to the circle of radius a located at the origin, an image vortex of opposite strength must be
Hint: Use equation 16 as a starting point to prove that the streamline is circular. Work in
polar coordinates. (3 Marks)
(b) Obtain the x and y velocity components for the streamfunction given by Eq. (16). ghost ghost
ghost ghost (1 Mark)
(c) Use Γ = 2π, b = 0.5 and a = 1 and draw the streamlines for the flow as described by Eq.
(16) and show that there is a streamline that correspond to the perfect circle of radius 1, i.e.
(d) Because the function is singular at (x, y) = (b, 0) and (x, y) = (a
2/b, 0), it is proposed that
we regularize the functions by adding δ
2
in the argument for the logarithmic function, i.e. we
would express the streamfunction as
i. Show that the streamfunction given by Eq. (17) has a constant value on the
ii. Obtain the x and y velocity components for the streamfunction given by Eq. (17). ghost
ghost ghost (1 Mark)
iii. Use Γ = 2π, b = 0.5 and a = 1 and draw the streamlines for the flow as described by Eq.
(17) for δ = 1.0, 0.2, 0.1 and 0.05 and show that there is a streamline that is very close to
the perfect circle of radius 1, i.e. x
2 + y
2 = 1 when δ is small. (4 Mark)
6. (12% of total mark) It is desired to mix a “blob” of dye placed in a square −0.1 ≤ x ≤ 0.1,
0.0 ≤ y ≤ 0.2 in a stirring (or mixing) tank of diameter 2m. Two electrical agitators will be placed
at (x, y) = (±b, 0) in an attempt to mix up the flow. A schematic of this practical situation is
shown in Fig. 1.
6
Cylinder wall
Figure 1: Schematic of the practical problem investigated in this assignment
Mathematically, the flow shown in Fig. 1 can be modelled by the following streamfunction
ψ(x, y) = ψ1(x, y) + ψ′1
(x, y) + ψ2(x, y) + ψ′2
The mixing of the “blob” of dye can be simulated by randomly placing fluid particles in the region
−0.1 ≤ x ≤ 0.1, 0.0 ≤ y ≤ 0.2. Each of the fluid particle will be convected by the local fluid velocity,
where xp and yp are the x and y coordinate of the corresponding fluid particle respectively.
7
(a) Obtain expressions for the x and y component of the velocity field given by Eq. (18). ghost
ghost (1 Mark)
(b) Using your program written in part 3(e) to draw streamlines corresponding to the streamfunction
given by Eq. (18). Use Γ = 2π, a = 1, b = 0.5 and δ = 0.1. ghost ghost ghost ghost ghost
ghost (2 Marks)
(c) Write a computer program to simulate the evolution of Np fluid particles representing the
“blob” of dye. Again, use Γ = 2π, a = 1, b = 0.5 and δ = 0.1. As a guide, your solution at
at time t = 0.5 with Np = 1, 000 should look similar to Fig. 2. Note that the computational
resources that you need to perform the calculations increases dramatically with Np. If you
have access to a good computer, you might like to use Np = 1, 000. For those of you who have
access to a very good computer, and have lots of time to do the calculations, use Np = 10, 000.
ghost ghost (5 Marks)
(d) Plot the location of all fluid particles at times t = 0, 1, 2, 3, 4 and 12. (3 Marks)
(e) Is this a good configuration for mixing the “blob” of dye at the specified location? ghost ghost
Figure 2: Solution at time t = 0.5
8
7. (10% of total mark) You have recently graduated with an engineering degree from the University
of Melbourne and have been hired by an engineering consulting company called Chan Inc. This
company has been given a lucrative contract to look at different ways of mixing fluids in a tank
similar to that shown in Fig. 1. Being a young and dynamic engineer, you would like to propose
that a better method of mixing the “blob” of dye might be to alternately run only one of the
agitators for a given time interval. In other words, you think that the company should run only
agitator A for (1/2)T and then turn it off and run only agitator B for another (1/2)T and keep
alternating them. However, being a new employee in this company, you do not want to make a fool
of yourself so you want to make sure that what you are proposing is feasible before you go to your
supervisor with this proposal. Thus, you would like to run some computer simulations to test out
your hypothesis.
Mathematically, the flow situation mentioned above can be modelled by the following streamfunction
where n = 0, 1, 2, 3......, and T is the period.
(a) Write down the x and y velocity components for the unsteady flow given by the streamfunction
above. (2 Marks)
(b) Write a computer program to simulate the evolution of fluid particles. Carry out simulations
with Γ = 2π, a = 1 and the following parameters
Use δ = 0.1 and a = 1.0 for all cases. Plot the position of all fluid particles at time t =
0, 1, 2, 3, 4, 5, 6, 9, 12 for all cases mentioned above. (7 Marks)
(c) Which of the cases above do you think results in the best mixing of the “blob” of dye ? ghost
ghost ghost (1 Mark)
(TOTAL MARKS AVAILABLE = 100)
9
Appendix: Numerical solution to ordinary differential equations
Consider the following ordinary differential equation,
d
dtx(t) = f(x). (23)
If f(x) is a relatively simple function, then one can obtain exact solution to Eq. (23) 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.
(23). 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:
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
dtx(t) = f(x, y)
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
10
4th order Runge-Kutta method:
where the approximated solution of x(t = tn) and y(t = tn) are denoted as xn and yn respectively. kij is
calculated as follows
k41 = f (xn + hk31, yn + hk32)
k42 = g (xn + hk31, yn + hk32)
11
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.
因为专业，所以值得信赖。如有需要，请加QQ：99515681 或邮箱：99515681@qq.com
微信：codehelp