- May 15, 2020

THE UNIVERSITY OF MELBOURNEDEPARTMENT OF MECHANICAL ENGINEERINGMCEN90008 FLUID DYNAMICSASSIGNMENT FOR POTENTIAL FLOWInstructions:• Assignment to be handed in by 23:59 on Friday 13th September 2018.• This assignment should be done in groups of 2 students. Both studentsin the group will get the same mark for this assignment. If you chooseto do the assignment alone, no concession will be given. Your assignmentwill 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 thisassignment. If you should decide to write your computer program inMATLAB, you are not allowed to use the streamline, odexx or similarfunctions in MATLAB i.e. you need to write the program to solve theordinary differential equations yourself and not simply use the functionsin MATLAB. Also do not merely use the contour function of ψ to visualisestreamlines. The aim of this assignment is for you to produce a set of toolsto 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 exercisesthat 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 achievedby setting daspect([1 1 1]) In Matlab).11. (20% of total mark) In order to write a computer code to sketch flow patterns and predictunsteady 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 thecomputer 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)dxdt= −x3(1)with the initial conditionx(t = 0) = 10.0. (2)(1 Marks)(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 equationfrom t = 0 to t = 5. Perform your computations with various step sizes, (h = 1, 0.1, 0.01), andtabulate your results in the table shown belowEU RK-4t Analytical answer h = 1 h = 0.1 h = 0.01 h = 1 h = 0.1 h = 0.011.02.03.04.05.0The output from your computer program must approach the analytical answer as h getssmaller. 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 ODEdxdt= −2x− 5ydydt= 4x+ 2y (3)You are given that at time, t = 0{xy}={−16}. (4)(3 Marks)(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 andy. To save time in future questions, you should try to write your program such that you caneasily change the set of equations f(x, y) and g(x, y) as defined in the Appendix.EU RK-4t Analytical answer h = 0.1 h = 0.01 h = 0.001 h = 0.1 h = 0.01 h = 0.0010.40.81.21.62.02Which is the more accurate numerical method? Use the more accurate numerical method inALL 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. Shouldthe 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, thestreamline 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 thatyou can pass it a number of starting positions xs and ys, the time step h and total time ofintegration tf . In matlab, this could be achieved by writing your program as a function whereyou pass the variables xs and ys (which could be vectors) and tf and h. This function wouldthen 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 fileshould contain the location of the initial points of the streamlines, (x0i y0i). For example, theinput file to draw two solution trajectories that start from (0, 1) and (2, 0) with h = 0.02 and0 < t < 2 will look something like,2 0.02 20.0 1.02.0 0.0(3 Marks)2. (20% of total mark) In lectures, it was demonstrated that using the method of images, the flowof 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ψ(x, y) =Qs2piarctan(y − lsx)+Qs2piarctan(y + lsx)−Qe2piarctan(y − lex)−Qe2piarctan(y + lex).(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 hasstrength Qe and the strength of the source is Qs.(a) Show that the velocity components associated with the above stream function shown in Eq.(5) are,u(x, y) =dxdt=Qs2pi(xx2 + (y − ls)2+xx2 + (y + ls)2)−Qe2pi(xx2 + (y − le)2+xx2 + (y + le)2)(6)v(x, y) =dydt=Qs2pi((y − ls)x2 + (y − ls)2+(y + ls)x2 + (y + ls)2)−Qe2pi((y − le)x2 + (y − le)2+(y + le)x2 + (y + le)2)(7)(3 Marks)3(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 programshould accept an input that gives the start positions in x and y for the streamlines, the timestep and the time of integration. Assume h = 0.01Please start your streamlines in sensible places to really show the important flow features. Tobegin with I suggest 20 starting positions around the smokers at radius = 0.1, equally spacedin the azimuthal direction. Also add other streamlines coming from the right, left and top ofyour figure. Refine this with streamlines placed close to any separatrix lines and stagnationpointsFor 10 radially spaced positions about the smokersx y0.5000 1.50000.4045 1.79390.1545 1.9755-0.1545 1.9755-0.4045 1.7939-0.5000 1.5000-0.4045 1.2061-0.1545 1.02450.1545 1.02450.4045 1.20610.5000 1.5000There are singularities at the center of the source and sink. You will need to add some errortrapping to catch these. One rather unsophisticated way is to prevent the velocities going toohigh,% if V > threshold, then set it to zeroV(V > threshold) = 0;It would be better to prevent your streamlines from getting too close to the singularity. ghostghost (7 Marks)(c) Derive analytical expressions for the location of all of the stagnation points in the flow. Plotall separatrix lines onto your diagram. What is the safe distance for the non-smoker? ghostghost (4 Marks)(d) Sketch the flow pattern given by Eq. 5 for Qs = 2 and Qe = 1. Again derive expressions forthe location of the stagnation points. Plot all separatrix lines onto your diagram (this mightrequire a modification to your program). (6 Marks)3. (10% of total mark) After carefully designating a safe distance for the non-smoker, someonecomes along and opens the window by a tiny amount, believing that this will help. A very gentle0.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, markswill be given for starting streamlines in places that allow us to best visualise the flow. ghostghost (4 Marks)(c) Write out expressions for the stagnation points at ±x on the floor. Solve this equation by anymeans 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-clockwisecirculation) is given by4ψ(r) = −Γ2piln r (8)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 asu(x, y) = −Γ2piyx2 + y2(9)andv(x, y) =Γ2pixx2 + y2. (10)(1 Mark)(b) Show that the time it takes for a fluid particle to circulate around this vortex, tc, is given bytc =4pi2r2Γ(11)(1 Mark)(c) Note that u(x, y) and v(x, y) are singular at (x, y) = (0, 0). Analyticaly, this is not usuallyan issue. However, when you are writing computer program to solve problems, functions thathave singularities will usually cause your program to “blow up”. In order to regularize u andv, it is proposed to use following streamfunctionψδ(r) = −Γ2piln√r2 + δ2as the streamfunction for the vortex instead of Eq. (8).i. Show that the Cartesian components of velocity, u and v can be written asuδ(x, y) = −Γ2piyx2 + y2 + δ2(12)andvδ(x, y) =Γ2pixx2 + y2 + δ2. (13)(1 Mark)ii. Show that the distribution of vorticity, ωδ = −∇2ψδ, is given byωδ(r) =Γ δ2pi (r2 + δ2)2(14)(1 Mark)iii. Plot uδ(0, y) for −1 ≤ y ≤ 1 for δ = 0.2, 0.1 and 0.05 all on the same plot. Let Γ = 2pi.ghost ghost (2 Marks)iv. Based on your analysis in this section, comment on the effects of δ on your solution. ghostghost (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 twovortices of opposite signs (i.e. Γ = 2pi for one vortex and Γ = −2pi for another vortex) locatedat (0, 1) and (0,−1). Your program should read an input file where the first line contains thenumber of vortices, N . The subsequent N lines in the program should contain the location ofvortices given by x, y and the circulation of each vortex, given by Γ. For example, the inputfile for the situation given above will look something like520.0 1.0 6.2831853070.0 -1.0 -6.283185307(8 Marks)5. (13% of total mark) The streamfunction for a potential vortex placed at a position (b, 0) can bewritten asψ(x, y) = −Γ2piln√(x− b)2 + y2 (15)In subsequent parts of this assignment, it is desired to have a streamfunction which produces astreamline that corresponds to a circle of radius a located at the origin. The streamfunction givenby Eq. (15) does not have this property.(a) Show that in order to obtain a streamfunction that produces a streamline which correspondto the circle of radius a located at the origin, an image vortex of opposite strength must beplaced at (a2/b, 0), i.e.ψ(x, y) = −Γ2piln√(x − b)2 + y2 +Γ2piln√(x−a2b)2+ y2 (16)Hint: Use equation 16 as a starting point to prove that the streamline is circular. Work inpolar coordinates. (3 Marks)(b) Obtain the x and y velocity components for the streamfunction given by Eq. (16). ghost ghostghost ghost (1 Mark)(c) Use Γ = 2pi, 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.x2 + y2 = a2 = 1. (2 Marks)(d) Because the function is singular at (x, y) = (b, 0) and (x, y) = (a2/b, 0), it is proposed thatwe regularize the functions by adding δ2 in the argument for the logarithmic function, i.e. wewould express the streamfunction asψδ(x, y) = −Γ2piln√(x− b)2 + y2 + δ2 +Γ2piln√(x−a2b)2+ y2 + δ2. (17)i. Show that the streamfunction given by Eq. (17) has a constant value on the circle x2+y2 =a2 provided that b2δ2 is small. (2 Marks)ii. Obtain the x and y velocity components for the streamfunction given by Eq. (17). ghostghost ghost (1 Mark)iii. Use Γ = 2pi, 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 tothe perfect circle of radius 1, i.e. x2 + y2 = 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 placedat (x, y) = (±b, 0) in an attempt to mix up the flow. A schematic of this practical situation isshown in Fig. 1.6Cylinder wallFigure 1: Schematic of the practical problem investigated in this assignmentMathematically, 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(x, y) (18)whereψ1(x, y) = −Γ2piln√(x − b)2 + y2 + δ2 (19)ψ′1(x, y) =Γ2piln√(x−a2b)2+ y2 + δ2 (20)ψ2(x, y) = −Γ2piln√(x + b)2 + y2 + δ2 (21)ψ′2(x, y) =Γ2piln√(x+a2b)2+ y2 + δ2 (22)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,that isdxpdt= udypdt= vwhere 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). ghostghost (1 Mark)(b) Using your program written in part 3(e) to draw streamlines corresponding to the streamfunc-tion given by Eq. (18). Use Γ = 2pi, a = 1, b = 0.5 and δ = 0.1. ghost ghost ghost ghost ghostghost (2 Marks)(c) Write a computer program to simulate the evolution of Np fluid particles representing the“blob” of dye. Again, use Γ = 2pi, a = 1, b = 0.5 and δ = 0.1. As a guide, your solution atat time t = 0.5 with Np = 1, 000 should look similar to Fig. 2. Note that the computationalresources that you need to perform the calculations increases dramatically with Np. If youhave access to a good computer, you might like to use Np = 1, 000. For those of you who haveaccess 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 ghostghost (1 Marks)−1 −0.5 0 0.5 1−1−0.500.51xyFigure 2: Solution at time t = 0.587. (10% of total mark) You have recently graduated with an engineering degree from the Universityof Melbourne and have been hired by an engineering consulting company called Chan Inc. Thiscompany has been given a lucrative contract to look at different ways of mixing fluids in a tanksimilar to that shown in Fig. 1. Being a young and dynamic engineer, you would like to proposethat a better method of mixing the “blob” of dye might be to alternately run only one of theagitators for a given time interval. In other words, you think that the company should run onlyagitator A for (1/2)T and then turn it off and run only agitator B for another (1/2)T and keepalternating them. However, being a new employee in this company, you do not want to make a foolof yourself so you want to make sure that what you are proposing is feasible before you go to yoursupervisor with this proposal. Thus, you would like to run some computer simulations to test outyour hypothesis.Mathematically, the flow situation mentioned above can be modelled by the following streamfunc-tionψ(x, y) ={ψ1(x, y) + ψ′1(x, y) nT ≤ t <(n+ 12)Tψ2(x, y) + ψ′2(x, y)(n+ 12)T ≤ t < (n+ 1)Twhere 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 streamfunctionabove. (2 Marks)(b) Write a computer program to simulate the evolution of fluid particles. Carry out simulationswith Γ = 2pi, a = 1 and the following parametersCase T bA 0.1 0.5B 0.5 0.5C 1.0 0.5D 2.0 0.5E 4.0 0.1Use δ = 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 ? ghostghost ghost (1 Mark)(TOTAL MARKS AVAILABLE = 100)9Appendix: Numerical solution to ordinary differential equationsConsider the following ordinary differential equation,ddtx(t) = f(x). (23)If f(x) is a relatively simple function, then one can obtain exact solution to Eq. (23) using analyticaltechniques which you should have learnt in your maths course. However, when f(x) is a complicatedfunction, 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 ordi-nary 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, seereferences [1], [2] and [3]). The formula for both methods areEuler’s method:xn+1 = xn + hf(xn)4th order Runge-Kutta:xn+1 = xn +(16k1 +13(k2 + k3) +16k4)hwherek1 = f(xn)k2 = f(xn +h2k1)k3 = f(xn +h2k2)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 thesimulation. Generally, the smaller the value of h, the numerical solution that you obtain is more accurateand stable.If you are asked to solve a system consisting of a set of coupled ordinary differential equationsddtx(t) = f(x, y)ddty(t) = g(x, y),the approximate numerical solution can be obtained using the Euler and 4th order Runge-Kutta methods.The relevant formulae areEuler’s method:xn+1 = xn + f(xn, yn)hyn+1 = yn + g(xn, yn)h104th order Runge-Kutta method:xn+1 = xn +(16k11 +13(k21 + k31) +16k41)hyn+1 = yn +(16k12 +13(k22 + k32) +16k42)hwhere the approximated solution of x(t = tn) and y(t = tn) are denoted as xn and yn respectively. kij iscalculated as followsk11 = f(xn, yn)k12 = g(xn, yn)k21 = f(xn +h2k11, yn +h2k12)k22 = g(xn +h2k11, yn +h2k12)k31 = f(xn +h2k21, yn +h2k22)k32 = g(xn +h2k21, yn +h2k22)andk41 = f (xn + hk31, yn + hk32)k42 = g (xn + hk31, yn + hk32)11References[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.12