THE UNIVERSITY OF MELBOURNEDEPARTMENT OF MECHANICAL ENGINEERINGMCEN90008 Fluid Dynamics 2023PART I OF ASSIGNMENT FOR POTENTIAL FLOWInstructions:• Assignment to be handed in by 23:59 on Sunday 17th September 2023.• 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 use any programming language to complete this assignment.Should you decide to write your computer program in MATLAB, youare not allowed to use the streamline, odexx or similar functions in MATLAB i.e. you need to write the program to solve the ordinary differentialequations yourself and not simply use the functions in MATLAB. Also donot merely use the contour function of ψ to visualise streamlines. Theaim of this assignment is for you to produce a set of tools to enable youcompute the coordinates of pathlines. The same rules apply if you useother 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 exercisesthat 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 isachieved by setting daspect([1 1 1]) In Matlab).1Part IFor the first part of this assignment, you will write codes to accurately plot pathlines from given startingpositions. 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 writtenhere to study much more complex flows.1. (Total marks for question = 27) In order to write a computer code to sketch streamlines, it isbest 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) inthe form x = f(t)dxdt = −x34(1)with the initial conditionx = 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 equationfrom 0 ≤ t ≤ tf where tf = 0.2, 0.4, 0.6, 0.8 and 1.0. Perform your computations with variousstep sizes, (h = 0.1, 0.01, 0.001), and tabulate your results in the table shown below (show 4decimal places).EU RK-4tf Analytical answer h = 0.1 h = 0.01 h = 0.001 h = 0.1 h = 0.01 h = 0.0010.20.40.60.81.0The output from your computer program must approach the analytical answer as h getssmaller. Plot the analytical solution (x vs t) along with the EU and the RK4 solution forh = 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 ODEdxdt = x + y (3)dydt = −2x − yYou are given that at time, t = 0xy=12. (4)(3 Marks)2(d) Extend your program from part (1b) so that it can solve the set of equations given by the setof Equations (3) for 0 ≤ t ≤ tf . Check your program by completing the table below for bothx and y for various tf and time-steps h. To save time in future questions, you should try towrite your program such that you can easily change the set of equations f(x, y) and g(x, y) asdefined in the Appendix.EU RK-4tf Analytical answer h = 0.1 h = 0.01 h = 0.001 h = 0.1 h = 0.01 h = 0.001246810(8 Marks)Which is the more accurate numerical method? Use the more accurate numerical method inALL 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π. Plotthe EU and RK4 calculations on separate figures. Should the lines join up? Do the lines joinup? (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 modifyyour program such that you can pass it a number of starting positions xs and ys, the time steph and total time of integration tf . In matlab, this could be achieved by writing your programas 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 arethe 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 thenumber of streamlines, Nl, the time of integration tf and the time step h. The subsequentNllines 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.011.0 21.5 32.0 4(3 Marks)2. (Total marks for question = 19) This question is based on an exam question from the 2022potential flow exam. We wish to model a Pitot-static tube as shown in figure 1(a) using a sourceand a uniform flow as shown in figure 1(b). The source is offset in the positive x direction fromthe origin by distance s such that the stagnation point is located at the origin (x, y) = (0, 0). Thetotal 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 andobtain values for the source strength Q and the offset s. (4 marks)3U∞p∞ p0xtapps(a) D = 0.01mU∞p∞(b)syxFigure 1(b) Use the codes that you developed in question 2 to produce a plot showing the streamlinepattern associated with your source and sink. Please show the streamlines external of thePitot-static tube in black, and the internal streamlines (from the source) as blue, and lines ofseparatrix 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 timestep (h ≈ 1×10−5s) 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 stagnationstreamline shown in blue in figure 1(a). Show that,U∞ =s2(p0 − p∞)ρHowever, the Pitot-tube does not directly measure p∞, instead it probes the pressure ps froma static pressure tap some distance xtap downstream of the leading edge stagnation point.Ups =s2(p0 − ps)ρIf the velocity measured by the Pitot-static tube (Ups) is 2% higher than U∞, show that thepressure coefficient Cpps at the location of the static tap must be,Cpps =ps − p∞12ρU2∞≈ −0.04Hence plot filled colour contours of Cp around the Pitot-static tube body (the half Rankinebody) and determine the minimum distance xtap such that the Pitot-static tube measured4velocity Ups is only 2%. greater than the true oncoming freestream velocity U∞. Do not showCp 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 Cppsshow that if we assume xtap is large relative to D, we canapproximate that,xtap ≈Q2πU∞(p1 − Cp − 1)(3 marks)Part II3. (Total marks for question = 20) The streamfunction for a potential vortex (with positiveanti-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 coordinatesare given by (x0,y0).(a) Show that the time it takes for a fluid particle to circulate around this vortex, tc, is given bytc =4π2r2Γ(6)(2 Marks)The Cartesian components of velocity, u and v for this potential vortex are given byu(x, y) = −Γ2πy − y0(x − x0)2 + (y − y0)2(7)andv(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 anissue. 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 pr2 + δ2 (9)for the vortex instead of Eq. (5).Show that the Cartesian components of velocity, u and v can be written asuδ(x, y) = −Γ2πy − y0(x − x0)2 + (y − y0)2 + δ2(10)andvδ(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) shifttext (3 Marks)5(c) Four vortices of the form given in equation (9) with δ = 0.2, are arranged in an unboundedspace with the following properties to model a mixing flow due to four fixed stirrers.Vortex x0 y0 Γ δ1 -1.5 0 -1 0.22 1.5 0 -1 0.53 0 -1.5 1 0.54 0 1.5 1 0.5Use the quiver command in Matlab to plot the velocity vector field associated with thisarrangement. (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.matThis file contains the x and y locations of 10050 particles of fluid at time t = 0 within themixing 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 particlelocations and submit this in your report. The stirrers (vortices) are fixed in physical space anddo not mutually induce motion in each other. You only need to track the particle trajectoriesunder 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 useyour Euler codes to track these particles over 1001 steps of h = 0.1. Plot the final particlelocations 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 givenby the following expression,yt = 0.90.2969rxc− 0.1260 xc− 0.3516 xc2+ 0.2843 xc3− 0.1015 xc4(a) The 0018 designation tells us that the total maximum thickness to chord ratio of this profile is18%. 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 transformthe shifted circle to an airfoil with the same thickness to chord ratio (it is likely that you willneed to seek a numerical solution to this question). Hint: In the Joukowsky transformed zjplane, 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 shiftedcircle.zj = (zsc) + b2(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 overa 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 thecomplex potential function is,w = U∞z +a2z| {z }cylinder flow−a2ziΓ2πlog(z)| {z }potential vortex6Note 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 surfaceof the cylinder, stagnation points and any separatrix lines. It will be much easier for thisquestion if you work in the complex plane. (4 Marks)We are going to move towards a Joukowski transformation of this cylinder flow. Weeventually 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 + iy0Jowkowski transformation z3 = z2 +b2z2Clockwise rotation (to ensure horizontal freestream) z4 = z3e−iΘwhere Θ is a positive value. We also need to adjust the circulation strength Γ to ensure thatthe Kutta condition is satisfied (see parts d and e). We will perform these 4 transformationsover 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 + x0Assume an angle of attack α = 8◦. Use the bodily shift x0 such that you will obtain a sharptrailing 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 +b2z2and then rotate the flow to ensure that we have horizontal freestream velocity in the farupstream of the airfoil.z4 = z3e−iΘ(where Θ is the original angle applied in part c). Plot the results in the z4 plane, including thetransformed 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 aboutitself a circulation of sufficient strength to hold the rear stagnation point at the trailingedge.Derive an expression to alter the circulation Γ in the z plane such that this condition is met inthe z4 plane. Replot the results (including cylinder, streamlines, separatrix and the Jowkowskicircle) 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 −V2JU2∞where VJ is the magnitude of the velocity on the surface of the Jowkowski airfoil. The velocity onthe surface of the airfoil can be found from the following,VJ = Vcdzdz4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 thevelocity magnitude field around your airfoil (V4 =pu24 + v24). 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 fillfill 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 theJoukowski airfoil at α = 8◦. Plot also on the same figure (with a different colour line) the Cpvs x/c obtained experimentally from your airfoil laboratory at the same angle of attack (if youdo not have lab data at α = 8◦, you may compare your lab data and Jowkowski result at adifferent matched α, as close as possible to 8◦). Comment on differences and similarities. fillfill fill fill filll (6 Marks)(h) Finally, compute lift coefficients CL via integration of your obtained Cp distribution for theJoukowski airfoil. Perform this computation for a range of angles of attack α from 0◦to 20◦insteps of 2◦. Compare this result with the experimental data and with the analytical solutionL = ρΓU∞. Plot all three curves on the same figure using different colours. Compare anddiscuss the results from the three curves. (8 Marks)8Appendix: Numerical solution to ordinary differential equationsConsider the following ordinary differential equation,ddtx(t) = f(x). (12)If f(x) is a relatively simple function, then one can obtain exact solution to Eq. (12) 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.(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, 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) + 16k4hwherek1 = f(xn)k2 = fxn +h2k1k3 = fxn +h2k2k4 = 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)h94th order Runge-Kutta method:xn+1 = xn +16k11 +13(k21 + k31) + 16k41hyn+1 = yn +16k12 +13(k22 + k32) + 16k42hwhere 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 = fxn +h2k11, yn +h2k12k22 = gxn +h2k11, yn +h2k12k31 = fxn +h2k21, yn +h2k22k32 = gxn +h2k21, yn +h2k22andk41 = f (xn + hk31, yn + hk32)k42 = g (xn + hk31, yn + hk32)10References[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