Departments of Chemical Engineering and Mathematics
Mathematica Session 2

Mathematica Graphics:
Visualising Fluid Flows

In this session we are going to use Mathematica's graphical functions to visualise some simple fluid flows using streamlines. In later sessions, you'll be introduced to some of Mathematica's numerical solving functions, and apply them to the calculation of particle paths. Note that we will not explain any of the details of the fluid mechanics here -- you have covered, or will cover, some of these in your Fluid Mechanics course.

A new feature of Mathematica's graphics that you will meet in this session is its capability to generate and display moving, animated graphics. These are called "movies", and indeed they work on the same principle as the cinema: the user creates a series of still pictures (frames) and uses Mathematica to display them in rapid succession, creating the illusion of motion.

The input statements used in this session are rather substantial. However, they are largely repetitive. We suggest that, if you have not already done so,  you get the hang of the "Copy" (Ctrl-C) and "Paste"  (Ctrl-V)  features (in the Edit menu), and "Copy Input from Above" (Ctrl-L) in the Input menu.

1. The Streamlines for Flow Past a Cylinder

The streamlines for steady, laminar flow past a cylinder are described approximately by the contours of the following stream function, expressed in polar coordinates (r, θ):

ψ(r,θ)=v (r-((a^2)/r))sin θ  .

a is the radius of the cylinder and v represents the velocity of the flow. (Note that neither a nor v will affect the qualitative structure of the streamlines.) Set up a function definition for ψ in Mathematica with the command:

ψ[r_, θ_] := v (r - ((a^2)/r)) Sin[θ]

Notice the use of the combined "colon-equals" sign , ":=", and (on the left-hand side only) the underscore characters, "_". A function definition does not produce any output, so the only visual indication that you'll get is that a blue "In[...]" number will appear at the left. You can test this new function by typing things like

ψ[a, π]ψ[2, π/5]ψ[2.3231, 0.3452]

We are interested in making a contour plot of ψ, and that requires that we use Cartesian coordinates, (x, y). Fortunately, the transformation between the two coordinate systems is straightforward:

r=(Sqrt[(x^2)+(y^2)]) , θ = (tan^(-1))(y/x).

So ψ can be evaluated with Cartesian coordinates in the form

ψ[(Sqrt[(x^2) + (y^2)]), ArcTan[x, y]]

We have used a different syntax for ArcTan, instead of the usual ArcTan[y/x]. Why? Compare:

ArcTan[1/0]ArcTan[0, 1]

ArcTan[1/-1]ArcTan[-1, 1]

Now for the contour plot. We'll use the "/." substitution operation introduced last week to give values to a and v, and we'll remove the contour shading and increase the plot resolution:

ContourPlot[
      ψ[(Sqrt[(x^2) + (y^2)]), ArcTan[x, y]] /. {v -> 10, a -> 0.5}, {x, 
        -2, 2}, {y, -2, 2}, PlotPoints -> 30, ContourShading -> False]

The contours around (0, 0) seem rather strange. But, in fact, they are not physically relevant because there is no fluid present for r < a. A good way to screen out this irrelevant information is to insert a coloured disk of radius a onto the plot, thus:

ContourPlot[
      ψ[(Sqrt[(x^2) + (y^2)]), ArcTan[x, y]] /. {v -> 10, a -> 0.5}, {x, 
        -2, 2}, {y, -2, 2}, PlotPoints -> 30, 
      ContourShading -> False, Epilog -> {GrayLevel[0.5], Disk[{0, 0}, 0.5]}]

Experiment with using the Contours option in this plot. There are two different modes of use: if you specify a number, say

Contours -> 20

you will get 20 contours, equally spaced over the whole range of ψ values in the plot (as chosen by Mathematica, based on the x and y plot ranges). Alternatively, if you specify a list of values, say

Contours -> {-15, -10, -5, 0, 5, 10, 15}

you will get contours plotted for those specific values of ψ.

Exercises and explorations 1

1.	Investigate what the Table command does -- for example,

Table[n, {n, 1, 10}]Table[c, {c, -10, 10, 2}]

Use Table commands with Contours to specify some sets of contours.

2.	Show using contour plots that the streamlines become closer and closer to the lines y = constant when either a->0, or x and y become large.

2. The Streamlines for Flow Past a Rotating Cylinder

The following modification to the cylinder stream function serves to model (approximately) the case of steady flow past a cylinder when it is rotating with an angular speed ω:

ψ2[r_, θ_] := 
      v (r - ((a^2)/r)) Sin[θ] + ω (a^2) Log[r]

Exercises and explorations 2

1.	By making contour plots of ψ2, produce your own versions of the following streamline diagrams for the different flow configurations as ω is varied. We leave you to work out suitable plot ranges and Contours specifications for each case.

(Actually, it is difficult to produce the "spike" and "loop" contours in the bottom diagrams using ContourPlot alone. We used another command, ImplicitPlot, to get those particular curves and then combined the graphics with the Show command.)

2.	Indicate by hand on the diagrams above which way the fluid is going along each streamline -- assume that the fluid enters from the left. (This can be done with Mathematica as well, but it's a somewhat advanced exercise.)

3. The Streamlines for Flow Past a Flat Plate

Next, we will consider the case of steady flow past a flat plate, which is inclined at some angle α to the flow. The idea is the same as before: a stream function is defined in terms of coordinates which are not Cartesian, but by use of an appropriate transformation we can make contour plots in x and y. The difference here is that the stream function and the coordinate transformation are quite difficult to set up, making use of complex numbers and hyperbolic functions, so we are going to leave out the details and give you a "black box" result. This is the stream function expressed in terms of x and y, and defined as a function of the inclination angle, α:

ψ3[α_] := 
      Im[Cosh[ArcCosh[(E^(I α)) (x + I y)] - I α]]

N.B. the spaces between I and α, and I and y. I represents the unit imaginary number, (Sqrt[-1]), usually written as i -- but in Mathematica it must be input as a capital letter.

The length of the plate is going to be fixed at 2 units, and we have omitted a velocity factor. Here is the contour plot for α = π/4, where we have used Epilog to add in a picture of the plate:

ContourPlot[ψ3[π/4], {x, -2, 2}, {y, -2, 2}, 
      ContourShading -> False, Contours -> 15, PlotPoints -> 30, 
      Epilog -> {Thickness[0.015], RGBColor[1, 0, 0], 
          Line[{{-Cos[α], Sin[α]}, {Cos[α], 
                  -Sin[α]}}] /. α -> π/4}]

It's useful to wrap everything up into a single function of the angle α. Copy down the previous input and modify it to this definition:

plateFlow[α_] := 
      ContourPlot[ψ3[α], {x, -2, 2}, {y, -2, 2}, 
        ContourShading -> False, Contours -> 15, PlotPoints -> 30, 
        Epilog -> {Thickness[0.015], RGBColor[1, 0, 0], 
            Line[{{-Cos[α], Sin[α]}, {Cos[α], 
                  -Sin[α]}}]}]

The wiggling plate: A Mathematica movie

Now for the movie that you were promised at the beginning. We'll make use of Table to build up a sequence of frames as the inclination angle is varied from π/4 to 3π/4:

Table[plateFlow[angle], {angle, π/4, 3  π/4, π/8}]

To animate these frames, the first thing you have to do is double-click on the cell bracket that contains them all (on the right-hand edge of the Mathematica window). They should all disappear apart from the first frame in the sequence (consult a tutor if this isn't working for you). Then, either select Animate Selected Graphics from the Cell menu or (quicker) hold down the "Ctrl" key and press "Y".

You can change the speed or direction of the movie using the buttons at the bottom left of the Mathematica window. In this case, you want a wiggling (back-forward) motion: click on the "loop" button between the "back" and "forward" buttons. To stop, click the mouse button anywhere in the Mathematica window.

For a smoother animation (with more frames), change the angle step in the Table command from π/8 to a smaller value.