For an inviscid luid (i.e.  luid low without energy losses - no frictional forces), the low around a solid cylinder can be described by a potential How solution, given by

φ(r,θ) = U (r + a2r) cos(θ).                                                   (1)

Here φ is the potential function, r = “x2 + y2  is the distance from the cylinder centre, and θ = tan — 1 y北  the angle with the x-direction (see igure 1a).  The parameters U and a represent the free stream velocity and the cylinder radius, respectively.

e

2

1.5

y

1

0.5

r

θ

−0.5

x

(a)

−1.5

−2

−2       −1.5       −1       −0.5         0         0.5         1         1.5         2

x

(b)

Figure 1:  (a) deinition sketch; (b) contour plot of the low kinetic enery e and velocity vectors.

The advantage of a potential function is that the information of the two velocity components is embedded in φ, namely as:

u = ∂φ

v = ∂φ

or more compactly

一u = ∇φ

(2)

The objective of this tutorial is to determine the velocity ield. Take U = 1 m/s and a = 0.5 m, and consider a domain −4a < x < 4a and −4a < y < 4a.

a. Download the master ile ex02a .m from Blackboard. The irst step is to determine 2-D arrays for r(x,y) and θ(x,y). Use the following strategy:

1.  Deine two one-dimensional arrays for xi  and yj   with  N北   =  100 and Ny   =  120 elements, respectively;

2.  Pre-allocate the two-dimensional arrays X , Y , r and θ using zeros. Here X and Y are coordi- nate arrays that are required for data visualisation later on.

3.  Loop over all xi  and yj . For each cell, calculate

(a) the 2-D x-coordinate Xi,j  = xi  and y-coordinate Yi,j  = yj ;

(b) the radius ri,j ;

(c)  θi,j  (hint: use the function atan2 to calculate θ, search the help for details);

Using subplot create a igure with two plots next to each other.  In the left plot, display r(x,y) with 10 contour lines, and in the right plot, display θ(x,y) with 40 contour lines.  Use the function contourf (use the 2-D arrays X and Y for the irst two input arguments), and use axis equal tight to set the correct igure aspect-ratio. Are these igures as you would expect?

b. Create another for-loop, and calculate φi,j . Check for each i, j whether ri,j  < a. If this is the case, set φi,j  = NaN (this tells MATLAB that φ is not deined at these points, and consequently MATLAB will not plot in this area).

Use contourf to plot φ as a function of x and y .

c. Create another for-loop. Calculate u and v by taking central diferences:

ui,j  =

vi,j  =

φi+1,j  − φi— 1,j 2x

φi,j+1  − φi,j — 1 2y

(3)

(4)