This preview shows page 1. Sign up to view the full content.
Unformatted text preview: ME 563  Intermediate Fluid Dynamics  Su Lecture 0  Visual fluids examples
Fluid dynamics is a unique subject because it's very visual. The book "An Album of Fluid Motion" by Milton van Dyke (on reserve at Wendt) is a collection of fascinating images from fluids experiments. You really can't say you understand fluids if you only think of it in terms of equations. Figure 1 is a top view of a triangular wing immersed in a flow of water (moving from left to right in the image). Colored fluid, appearing as white, is introduced near the leading edge. The wing Figure 1: Turbulent transition in flow over a wing. is inclined at a 20 angle of attack. Initially the fluid pattern is very smooth, then the filaments of colored fluid make a very abrupt transition to turbulence. The abruptness of the turbulent transition is interesting for many reasons, not least of which is that it's not really predicted by the equations of fluid flow. Another interesting property of fluid flows is that the organization of it is very persistent. The upper left photo in Fig. 2 is of the wake of a circular cylinder in a water flow. The cylinder is at the left edge of the photo. The mean flow is very slow about one cylinder diameter per second. The Reynolds number is 105. The alternating pattern of eddies (vortices) is called the `K`rm`n a a vortex street.' Intuitively, it kind of makes sense that a slow, laminar flow would be very organized. The upper right photo in Fig. 2, showing the wake of a plate at a 45 angle of attack, is taken at a flow that is, relatively, about 40 times faster (Reynolds number 4300). The flow in this case is turbulent, but even so, the alternating pattern of eddies is visible above the randomness. To drive home this point further, the lower photo in the figure is of the wake of a tanker inclined at roughly 45 to the mean current. The pattern assumed by the oil slick is amazingly similar to that in the upper right photo, even though the Reynolds number of the ship wake is on the order of 107 . The sheer range of length scales that are interesting in fluids problems is pretty remarkable too. Active research ranges from flow in microchannels (blood flow in capillaries, for example) all the way up to cosmic systems like nebulae. The upper left photo in Fig. 3 shows the KelvinHelmholtz instability, which arises in the interface between flow streams at different velocities. The photo shows a rectangular tube (the 18 inch ruler in the image shows the scale), with pure water moving lefttoright on top and colored salt water moving righttoleft below. The upper image on the right is of clouds near Denver. On the leeward (sheltered) side of mountain ranges, you will often have layers of high winds above relatively slow air. When clouds sit near the interfaces, these KelvinHelmholtz structures result. Needless to say, the cloud structures are much bigger than those seen in the lab, but the math is the same. (The cloud image is taken from "Hydrodynamic Stability" by Drazin and Reid, also on reserve.) 1 Figure 2: Persistence of flow organization, for Reynolds numbers spanning five orders of magnitude. Figure 3: Fluid problems span a vast range of length scales. 2 ME 563  Intermediate Fluid Dynamics  Su Lecture 1  Math fundamentals
The subject of fluid dynamics brings to mind pipes, pumps, wind tunnels, fans and any number of other engineering devices. However, to study fluid dynamics effectively requires some basic mathematical tools. (Interestingly, while fluid dynamics is thought to be very much the domain of engineers, much work in fluids, particularly in Europe, is performed in applied math, math, applied physics and physics departments.) In this class it will be assumed that you are familiar with the basic techniques of differentiation and integration, partial derivatives, differential equations and vector calculus. You may not remember everything from your math courses but you should at least know where to find things (tables of integrals, for example). What is important is your mathematical intuition. The following expression should look vaguely familiar: f (x0 + x)  f (x0  x) =? x0 2x lim If a function f (x) is differentiable at a point x0 , then this gives the value of the derivative, f (x0 ), at that point, that is f (x0 + x)  f (x0  x) = f (x0 ) if f is differentiable at x0 . x0 2x lim (1) This expression gives the value of f (x0 ) if f is differentiable, but just because the limit in (1) exists doesn't mean that f is differentiable at x0 (what's a counterexample?). It turns out, if the onesided limits are equal f (x0 + x)  f (x0 ) f (x0 )  f (x0  x) = lim x0 x0 x x lim (with x > 0) (2) then f is differentiable at x0 , and f (x0 ) is equal to the value of the onesided limits of (2), which is then equal to the twosided limit of (1). Thinking of derivatives in terms of differences (hopefully) seems completely trivial, but understanding (1) and (2) implicitly is key to subjects like fluids that rely heavily on differential equations (and if it were easy, it wouldn't have taken Isaac Newton to come up with calculus). This is especially true since data (whether experimental or computational) is digital, and discrete math is used extensively. 1 Vector analysis Vector analysis is vital to fluid mechanics because the most fundamental descriptive quantity in fluid flow, the velocity, is a vector. We'll start with the threedimensional Cartesian coordinate system with unit vectors ex , ey and ez . Consider the vectors a = ax ex + ay ey + az ez and b = ^ ^ ^ ^ ^ ^ bx ex + by ey + bz ez . The dot product (or scalar product) of a and b is ^ ^ ^ a b = a b cos = ax bx + ay by + az bz (3) where   denotes the length, or magnitude, of a vector, and is the angle between the vectors. The cross product (or vector product) of a and b, a b, is more complicated. Its magnitude is given by a b = a b sin 1 Figure 1: The direction of the cross product a b. where is the (smallest) angle between a and b (i.e. 180 ). The direction of a b is perpendicular to the plane of a and b, in the righthand sense of rotation from a to b through the angle (Fig. 1). Because of this direction definition, the cross product anticommutes, i.e. a b = b a. It turns out that the most convenient way to express the cross product is as a determinant ex ey ez ^ ^ ^ a b = ax ay az = (ay bz  az by )^x + (az bx  ax bz )^y + (ax by  ay bx )^z e e e bx by bz (4) Easily seen through (4) are the identities ^x ey = ez , ey ez = ex , and ez ex = ey . e ^ ^ ^ ^ ^ ^ ^ ^ Of particular interest is the operator (variously called the gradient operator, del operator, or grad operator). We can write as = ex ^ + ey ^ + ez ^ x y z (5) has no meaning unless it's operating on something. The familiar operations are the gradient, the divergence, the curl, and the Laplacian. 1.1 The gradient and directional derivative
operates on a scalar quantity, (x, y, z), the result is called the gradient of , and is = ex ^ + ey ^ + ez ^ x y z (6) When written To interpret this, consider a unit vector ^ = sx ex + sy ey + sz ez , and let s be the distance variable s ^ ^ ^ along this vector. Without loss of generality, we will assume that s = 0 at the origin, x = 0. Then, the coordinates of any point on the saxis are x = sx s y = sy s z = sz s (7) Suppose we want to know the rate of change of in the s direction, at the origin. This d/ds is also known as the directional derivative. By the chain rule, this can be written (note the selective use of partial derivatives) d ds = = dx dy dz + + x ds y ds z ds sx + sy + sz x y z 2 (8) Figure 2: Sample volume for illustrating the divergence. where we have used (7). By inspection, we have d = ds ^ s (9) that is, the rate of change of a scalar function in an arbitrary direction is equal to the scalar product of the gradient with the unit vector, ^, in that direction. Using (3), we can also write s (since ^=1) s d =   cos ds (10) where is the angle between and es . So we can also say that d/ds is the projection of ^ onto the direction s. Finally, observe from (10) that is perpendicular to lines of constant . 1.2 The divergence and divergence theorem Now let's consider what happens when we apply to vector quantities. Consider a vector field f = fx ex + fy ey + fz ez . (We call this a vector `field' because it's defined over a volume in space, ^ ^ ^ not just at a single point.) The dot product of with f is called the divergence, and is given by f = fx fy fz + + x y z (11) To understand the name `divergence', let the threedimensional volume V be a cube with infinitesimal side lengths dx = dy = dz, as in Fig. 2. We can write the volume of V as dx dy dz = dV . Assume that on each face of the cube, f is constant. Now consider face 1 in the figure. The component of f that points out of the cube on face 1 is fx . Because f is constant on the face, we can write this as fx (x = dx). On face 4, the component of f that points out of the cube is fx (x = 0). The area of each of these two faces is dy dz, so the net volume flux out of the cube through faces 1 and 4 can be written Flux out of faces 1 and 4 = [fx (x = dx)  fx (x = 0)] dy dz We can go through similar arguments for the remaining faces, and we get Flux out of faces 2 and 5 = [fy (y = dy)  fy (y = 0)] dx dz Flux out of faces 3 and 6 = [fz (z = dz)  fz (z = 0)] dx dy (13) (12) 3 Now notice that we can rewrite (12) as [fx (x = dx)  fx (x = 0)] dy dz = = fx (x = dx)  fx (x = 0) dx dy dz dx fx (x = dx)  fx (x = 0) dV dx fx dV x (14) where the last part is true in the limit of dx approaching zero. Similarly, the flux out of faces 2 and 5 can be written (fy /y) dV , and out of faces 3 and 6 can be written (fz /z) dV . Referring back to (11), given a vector f , the net volume flux out of the volume dV equals the product of f and dV . Thus f is a measure of the spread, or divergence, of the vector field f . This somewhat nonrigorous analysis can be generalized as the divergence theorem. Given an arbitrary volume V , enclosed by the surface S, with outward unit normal vector n at all points on S, the divergence theorem states f dV = f n dS. (15) V S This will be very useful when we get to the equations for fluid flow. 1.3 The curl and Stokes' theorem
with the vector f is called the curl, and is given by ey ^ y The cross product of f = ex ^ x ez ^ z = fx fy fz fz fy  y z ex + ^ fx fz  z x ey + ^ fy fx  x y ez ^ (16) To see where the term `curl' comes from, look at Fig. 2 again, but this time only consider face 6. Call this the surface S, with area dx dy = dS. We'll define a direction of travel around the border of S in the counterclockwise direction. Assume also that the vector field f is constant on each edge of S. Starting at the the origin, we first travel along the edge defined by y = 0. The component of f parallel to the direction of travel on this edge is fx (y = 0). We can then define a contour integral (a sort of net travel) as fx (y = 0) dx. The next edge is the one defined by x = dx, along which the net travel fy (x = dx) dy. Going all the way back around to the origin, we end up with (taking care with the signs) Net travel = fx (y = 0) dx + fy (x = dx) dy  fx (y = dy)) dx  fy (x = 0)) dy = [fy (x = dx)  fy (x = 0)] dy  [fx (y = dy)  fx (y = 0)] dx fy (x = dx)  fy (x = 0) fx (y = dy)  fx (y = 0) = dS  dS dx dy fy fx  dS x y (17) Again, this last expression is true for dx and dy approaching zero. Comparing with (16), (17) is just the zcomponent of f . By the direction convention for integration around a closed contour, ez is the normal vector for face 6 (our surface S) in Fig. 2 for the counterclockwise integration ^ direction. Thus, the integral of f around the contour enclosing S equals the component of f in the direction normal to S multiplied by dS. In the context of the integration around the closed contour, the use of the term `curl' is obvious. 4 ^ Figure 3: Direction convention for integration around a closed contour, C. The vector n is the unit normal to the surface. Stokes' theorem generalizes this. Let S be a twodimensional surface enclosed by the curve C. Then (
S f ) n dS = C f dx (18) where the direction of the line integral relates to the direction of the normal vector n in the righthand sense (Fig. 3). 1.4 The Laplacian
2 2 2 + 2 + 2 x2 y z If we take the divergence of a gradient, we get the Laplacian, which is defined ( is a scalar)
2 = ( ) = (19) applying (6) and (11). The Laplacian is interesting in relation to quantities that undergo gradient diffusion. An example is heat, which diffuses proportionally to the gradient in temperature. If we let be the temperature and be the thermal diffusivity, then the temperature flux vector can be written (note the minus sign) Temperature flux vector =  (20) The Laplacian can thus be interpreted by substituting  for f in the discussion of the divergence in Sec. 1.2. That is, the Laplacian describes the net flux of the scalar quantity into a volume. Consider a uniform temperature field with a sharp positive spike somewhere in it. The spike will have a strongly negative Laplacian value, which means that heat will flow out from the region of the spike, which makes sense, since diffusivity tends to smooth out sharp gradients. Another scalar quantity that undergoes gradient diffusion is species concentration a vector quantity that does is momentum, but we will need more math tools to consider the diffusion of vector quantities. 5 ME 563  Intermediate Fluid Dynamics  Su Lecture 2  More math, plus some basic physics
In the first lecture we went over some basic math concepts, in particular the operator . We'll finish up with basic math by going over tensors and index notation, then talk about some basic physical concepts. 1 Tensors and index notation In the last lecture we considered the dot product, where two vectors result in a scalar, and the cross product, where the two vectors yield a third vector. There is another way to multiply vectors together that gives rise not to a scalar or a vector, but to a tensor. To deal with those it's most convenient to use index notation (also called tensor notation). 1.1 Index notation Scalars and vectors are actually specific cases of tensors. A scalar is a tensor of order (or rank) zero, and a vector is a tensor of order one. (Generally, however, if we say that something is a tensor without specifying its order, we will mean that it is a tensor of order two.) The order of a tensor tells you the number of indices necessary to describe it. In ndimensional space, a tensor of order m has nm components. We can repeat some of the results of the last lecture using index notation. Consider a vector f defined in threedimensional space. Let the three orthogonal unit vectors be e1 , e2 , and e3 . (We're ^ ^ ^ not using x, y and z because the coordinate system is not necessarily Cartesian, or it could be rotated from x, y and z, etc.) Then we can write f as f = f1 e1 + f2 e2 + f3 e3 = fi ei ^ ^ ^ ^ This expression illustrates two key aspects of index notation: An index (in this case i) takes on values corresponding to the number of dimensions in the space being considered. If an index is repeated in a term, then that term is summed over that index (this is the summation convention, which physicists call the Einstein summation convention). The dot product and cross product of two vectors can be expressed conveniently in index notation, with the aid of two new operators. The dot product of two vectors a and b is a b = ai bj ij , where ij = 1 0 if i = j if i = j. (2) (1) This ij is called the Kronecker delta. (Of course, we could also have written a b = ai bi using the summation convention.) The cross product of a and b is 1 if ijk is cyclic, i.e. 123, 231, or 312 a b = ijk aj bk ei , where ijk = 1 if ijk is anticyclic, i.e. 321, 213, 132 ^ (3) 0 if any of the two indices are identical. The
ijk is called the permutation tensor or permutation operator. 1 1.2 Tensors Consider two vectors a = a1 e1 + a2 e2 + a3 e3 = ai ei and b = b1 e1 + b2 e2 + b3 e3 = bi ei . If we multiply ^ ^ ^ ^ ^ ^ ^ ^ them together, not as a dot product or cross product, but just ab, we get ab = T, or a1 b1 a1 b2 a1 b3 ab = ai bj = a2 b1 a2 b2 a2 b3 = Tij = T (4) a3 b1 a3 b2 a3 b3 where, by convention, the first index, here i, represents the row of the tensor and the second index, j, represents the column. Also, a quick notational point. The order of a tensor is sometimes identified by underlining. So a vector a is often written a, and a tensor T is often written T . The underlining convention is common with written work; in printed work, vectors and tensors are usually just represented in boldface (some authors use lowercase bold for vectors and uppercase bold for tensors, but this isn't universal). Because the first index in a tensor represents the row, and a vector is a tensor of order one, the vector a can be written a1 a2 . a = ai ei = (5) ^ a3 The matrix representation of a is not written with the unit vectors ei because in tensor form, the ^ unit vectors are assumed to go with the components i (this can sometimes be confusing). The transpose of a tensor is a tensor with its rows and columns switched. For a secondorder tensor T, this means the transpose TT is given by T11 T12 T13 T11 T21 T31 for T = Tij = T21 T22 T23 , the transpose is TT = Tji = T12 T22 T32 . (6) T31 T32 T33 T13 T23 T33 In the expression for the transpose, when we write TT = Tji , the first index still corresponds to the row even though we use j instead of i. A tensor T is called symmetric if T = TT , i.e. if Tij = Tji . A tensor T is antisymmetric if T = TT , or Tij = Tji . For a vector a, taking the transpose just means a1 for a = a2 , the transpose is aT = a1 a2 a3 . (7) a3 The standard rules of matrix multiplication can be applied to tensors. Matrix multiplication corresponds to the dot product. It is possible to take the dot product of a tensor T and a vector a a1 T11 a1 + T12 a2 + T13 a3 T11 T12 T13 T a = Tij aj = T21 T22 T23 a2 = T21 a1 + T22 a2 + T23 a3 (8) T31 T32 T33 a3 T31 a1 + T32 a2 + T33 a3 (Notice the application of the summation convention.) The dot product in this case doesn't yield a scalar, but instead a vector. The dot product essentially gives a result one order reduced from the highest order tensor in the product. This leads to the other form of multiplication we'll be interested in, the doubledot product (also known as the scalar product for tensors). This takes two secondorder tensors and yields a scalar, as T : T = Tij Tij . 2 (9) 1.3 Velocity gradient tensor
u. For Cartesian coordi The tensor we will be most interested in is the velocity gradient tensor, nates, with u = u^x + v^y + w^z , u is defined as e e e u x uj u= = u y xi u
z v x v y v z w x w y w z (10) The Laplacian of the velocity,
2 2 u, will also appear often in our discussions. This is found by u v w x y z u= ( u) = u = = 2 2u + u x2 y 2 2v 2v x2 + y2 2 2w + w x2 y 2 2 u + + + x v y y u v z z 2u z 2 2v z 2 2w z 2 x x w y w z (11) 2w 2v so taking the Laplacian of the velocity vector is equivalent to taking the Laplacian of each of the components. (We'll talk about the physical meaning of u in detail later.) 2 Physical concepts Let's now discuss some basic physical ideas. One of the goals of this class is to develop the ability to think through problems intuitively before having to bring math tools to bear. Suppose I have two airplanes. They are identical except one of them has winglets on the wingtips, like you see on some airliners (747's, A340's etc.). Figure 1: Two wings: wing 1, no winglet, wing 2, winglet. The first wing shows a tip vortex, while the second doesn't. Figure 2: Tip vortex on wing 1. Why does this tip vortex arise? First, we know that wings are designed to provide elevated pressures on the lower surface, and reduced pressures on the upper surface. However, if we pick a point just outside of the wingtip, the pressure has to be the same from above and below, because 3 the pressure can't have a discontinuity. As a result, the pressure on the underside of the wing is gradually reduced as we approach the tip, and the pressure above the wing is gradually increased as we approach the tip, so the pressures become equal outside the tip. Because fluids move in the Figure 3: Pressure distribution around a wing. direction of decreasing pressure, the airflow on the lower surface of the wing moves toward the tip, and the airflow on the lower surface moves toward the plane's fuselage. Put together, the result is a tip vortex with the sense of rotation shown in Fig. 2. (The winglet in wing 2 is designed to prevent the formation of this vortex.) Basically, two principles (and no math) were necessary to show this  the idea that pressure fields are continuous, and that flows are driven in the direction of decreasing pressure. Why, though, would a designer want to get rid of the tip vortex? Let's consider what effect the tip vortex has by imagining that we're stationary in space, and can somehow see the velocity of the air. A plane equipped with winglets passes, leaving the air essentially undisturbed. A plane without winglets passes, and in contrast we see the air spinning off the wingtips. The subsequent reasoning goes like this: The flow with the vortical structures has more kinetic energy than the undisturbed flow. The kinetic energy has to be provided by the plane's engines. The energy provided by the engines to the vortices is energy not provided to the plane's forward motion. The vortices thus correlate with an effective drag force on the plane. We've thus shown that the tip vortices cause drag, and we needed to apply only the conservation of energy, again with no math. Unfortunately...not all problems can be solved with just physical intuition, but physical intuition is still very handy. 4 ME 563  Intermediate Fluid Dynamics  Su Lecture 3  Basic concepts of ideal fluids / Euler's equations
Reading: Acheson, chapter 1. In describing fluid flow it is common to speak of the velocity field, u = u(x, t), which gives the value of the velocity vector, u, at all points x and times t. (In Cartesian coordinates, u = u ex + v ey + w ez = (u, v, w) and x = x ex + y ey + z ez = (x, y, z).) The following basic notions will ^ ^ ^ ^ ^ ^ be used extensively. A steady flow is one where the velocity is a function of spatial position only, i.e. u = u(x), and u/t = 0. A twodimensional flow is one in which the velocity is independent of one spatial direction, and there is no velocity component in that direction. In Cartesian coordinates, u = u(x, y, t) ex + ^ v(x, y, t) ey is an example of a twodimensional flow. ^ A streamline is a curve that is parallel to the local velocity vector at all points. Suppose s is a variable that measures the distance traveled along a particular streamline. In Cartesian coordinates, a streamline would be expressed x = x(s), y = y(s), z = z(s). At any point on the streamline, dx/dy = u/v, dy/dz = v/w, and dx/dz = u/w (compare Eq. 1.5 in the book). The material derivative Df /Dt describes the rate of change of a quantity f in a reference frame following the fluid. It is defined Df f f f f = +u +v +w Dt t x y z f f = + (u )f = + u ( f ). t t (1) The velocity components u, v and w appear courtesy of the chain rule. This Df /Dt is also known as the substantial derivative or particle derivative. The process described by Df /Dt is known as advection, and involves the change in f due only to the motion of the fluid. Other factors can affect f and must be dealt with separately (examples to follow). 1 Ideal fluids
1. It is incompressible. In practice, this generally means that the density, , is a constant throughout the fluid. 2. It is inviscid. Thus, the only force acting across any surface element dS in the fluid is the pressure. An ideal fluid is one for which the following conditions hold Strictly speaking, the property of incompressibility is better ascribed to the flow than the fluid. For example, gases are obviously compressible, but under certain conditions (specifically, low Mach number) their flows can be treated as incompressible. 1 Figure 1: An infinitesimal cubic volume V . 1.1 Equations of motion for an ideal fluid To get the equations of motion for an ideal fluid, we might reasonably expect to use the familiar conservation properties from physics: conservation of mass, momentum, and energy. Sure enough, the first equation of motion describes the conservation of mass. Let V be an arbitrary volume in the fluid, enclosed by the surface S. Because the density of the fluid is constant, the net volume flux across the surface S has to be zero. By the divergence theorem, this means u dV = 0 for any arbitrary volume V . (2) V We can make this statement stronger. Suppose that u is positive (or negative) at some point x0 . Assuming that the u is continuous, then we can define a volume V around x0 such that u is positive (or negative) at all points in V . But then this violates (2). Thus, the conservation of mass leads to the first equation of motion for an ideal fluid, namely u = 0 everywhere in an ideal fluid. (3) Now for the conservation of momentum. Again, we let V be some arbitrary volume enclosed by S. By the second property of ideal fluids listed above, the only forces acting across S are pressure forces. We are interested in the forces acting on the fluid in V ; we are not interested here in the forces applied by the fluid contained in V . The net force exerted on the fluid in V by the surrounding fluid is  pn dS (4) S where the minus sign indicates that the force vector points into V in the direction of the normal vector to S. To cast this in a more convenient form, we'll again do an analysis similar to what we did to illustrate the divergence in Lecture 1. Figure 1 is the same infinitesimal cubic volume V , with volume dV = dx dy dz, that we used before. The pressure force applied by the surrounding fluid on face 1 is p(x = dx) dy dz ex , and on face 4 is p(x = 0) dy dz ex . (As before, we have ^ ^ assumed that p is constant on each face.) The other faces can be treated analogously. We get Pressure force on faces 1 and 4 = [p(x = dx)  p(x = 0)] dy dz ex ^ Pressure force on faces 2 and 5 = [p(y = dy)  p(y = 0)] dx dz ey ^ Pressure force on faces 3 and 6 = [p(z = dz)  p(z = 0)] dx dy ez . ^ (5) 2 Look at the force on faces 1 and 4. We can rewrite that component as [p(x = dx)  p(x = 0)] dy dz ex =  ^ p(x = dx)  p(x = 0) dV ex ^ dx p  dV ex . ^ x (6) Compiling the x, y and zcomponents, and recalling the definition of the gradient, we get Net pressure force on V =  p p p dV ex  ^ dV ey  ^ dV ez = ( p) dV, ^ x y z (7) which tells us the net pressure force on the infinitesimal volume element. (Notational point: the text calls the volume V .) Besides the pressure force, we will also allow for gravitational forces on our volume element Gravitational force on V = mV g = g dV where mV = dV is the mass of fluid contained in V . Adding (8) to (7), we get Total force on V = ( p + g) dV (9) (8) This applied force has to be equal to the rate of change of the linear momentum of the fluid in V (the product of its mass and acceleration), that is, Mass acceleration of the fluid in V = dV Du Dt (10) where we have used the material derivative to express the acceleration of the fluid. We do this because the the forces described by (9) act on the fluid contained in V at a particular instant in time. Any changes in momentum effected by those forces are manifested in that parcel of fluid (what the book calls a `blob') as it moves. The material derivative is thus relevant because it is the derivative that follows the fluid. Equating (10) and (9), we get the remaining equations of motion for an ideal fluid, Du 1 = p + g. Dt (11) In general, (11) represents three equations, one for each dimension. They are nonlinear partial differential equations. Together, (3) and (11) are called Euler's equations of motion for an ideal fluid. (We didn't use the conservation of energy in deriving Euler's equations. Energy will become interesting when we talk about viscous fluids.) 3 Ҹ ֳ ݸ ֡ ֵ и ֵ ٸ غ ۸ ֵ ٵ ִ ٵ ٵ ١ٷ ֺ ݸ ص ظ ٺ Ѻ ݺ ٵ ظ Һ ٵ ٵ ظ Ӻ Ӹ ۺ ݸ ۸ ۸ Ӹ ٵ Ӹ ҳ غ ۸ к ݸ ֢ Ӹ ֢ ۺ г ظ ֢ ٵ и غ Ӹ ص ֢ ٵ ٵ ٵ ֵ ֵ ٵ ٴ ݸ ҳ ݺ к ظ ҳ ݺ ӹ ۸ ٴ ص ڴ ص ֵ ӹ ֵ ݸ غ ֵ ָ ݸ ֵ ٺ ӹ ҳ ֹ ݸ ӹ ӹ ӹ и Ҹ ׳ Ѹ ٸ ١ ٵ ١ ӹ и Ѻ ҳ ָ ׳ Ѹ Ҹ к Ӹ и Һ ׳ ME 563  Intermediate Fluid Dynamics  Su Lecture 5  Limitations of ideal fluid theory / role of viscosity
Reading: Acheson, chapter 1, chapter 2 (through 2.3). In the last lecture we introduced the concept of circulation, which we defined as =
C u dx where C is a closed curve. The circulation around an airfoil explains the lift generated. Consider the airfoil in Fig. 1, surrounded by a curve C, along which we will compute the circulation in the direction indicated. Assume that the mean airflow is horizontal, from left to right, in the reference Figure 1: Computation of the circulation around a wing. frame where the wing is stationary. To compute the circulation, first assume that the sides labeled 2 and 4 are far enough away from the wing that the flow has no vertical component there, so u dx = 0 on those sides. Also, because C is a rectangle, sides 1 and 3 have the same length, but on side 1 we integrate in the positive xdirection and on side 3 we integrate in the negative xdirection. The circulation around C then becomes
side 1 side 3 = = u dx +
side 1 u dx
side 3 ( u  u ) dx. (1) What do we know about the values of u above and below the wing? We can assume that the flow is steady and irrotational, and that gravity is negligible. By the (second) Bernoulli theorem, p + 1 u2 is constant everywhere. For the wing to have lift, the pressure on the lower surface has 2 to be greater than that on the upper surface, which then means that the air speed below the wing (along side 3 of C) has to be less than the air speed above the wing (side 1 of C). From (1), this means that the circulation, , is negative (for the direction of integration defined in Fig. 1). Bernoulli's theorems are formulated for ideal fluids, and while real fluids aren't ideal, ideal fluid theory does a good job of predicting lift for airfoils at small angles of attack. (At high angles of attack, ideal fluid theory fails rapidly, as in Fig. 1.11 in the text.) We have shown that the lift of an airfoil depends on the circulation around it. Let's think a little bit more about the circulation being nonzero. Specifically: We have assumed that the mean flow in which we've immersed our wing is irrotational (has zero vorticity), that the flow around the wing can be treated as twodimensional. These are perfectly valid assumptions. 1 By Stokes's theorem, we showed that the circulation around a curve C is zero if the vorticity is zero over the entire surface S bounded by C. By considering the vorticity equation in the last lecture, we showed that for twodimensional flows of ideal fluid subject to conservative body forces, the vorticity of each individual fluid element is conserved. Given all this, we might reasonably think: It appears that vorticity correlates with circulation, and the vorticity in the flow surrounding the wing is zero. So how does the circulation around the wing arise? 1 Viscosity To explain this properly, we have to consider what happens when a wing, initially at rest, is set into motion suddenly. (Described in 1.1 in the text.) It turns out that as the wing is set into motion, a vortex (with, of course, nonzero vorticity) forms at the trailing edge of the wing. As we will (hopefully) cover later, the circulation contained in the vortex is equal and opposite to the circulation around the wing. The vortex trails further and further behind the wing as time evolves. When we look at a wing in steady flow, we can assume that the vortex is, effectively, infinitely far behind the wing, so the wing can be considered to be surrounded by irrotational flow. However, the circulation around the wing, which balances the circulation in the vortex, remains. Here is a subtle point: while the existence of circulation around an airfoil is permissible by ideal fluid theory, we can't explain the origin of the circulation unless the viscosity of the fluid is considered. Viscosity can be a tricky concept because there are major qualitative differences in fluid flows with very small viscosity and those with no viscosity. That is, the behavior of a fluid with very small viscosity can be vastly different from the behavior of a fluid with zero viscosity. 1.1 Viscous forces One of the properties that define an ideal fluid is that the only force acting across any surface element in the fluid is the pressure, which acts in a direction normal to the plane. The simplest explanation of viscous fluids is that they have an additional force that acts tangential to a given surface element. This force depends on the velocity gradient across the surface element. A fast layer of fluid wants to drag a slow layer along with it, and the slow layer wants to slow down the fast layer. (Fig. 2.) As shown in the figure, while pressure is a normal force (meaning its direction), Figure 2: Pressure forces and viscous shear forces across a surface in a fluid. 2 viscosity imposes shear forces1 . Specifically, viscous fluids resist shear, i.e. they oppose velocity gradients. Many of the qualitative differences between viscous and inviscid fluids arise because viscous fluid flows have boundary layers near solid surfaces. Boundary layers arise because of the noslip condition for viscous fluids, which holds that a fluid element in contact with a solid boundary must have the same velocity as the boundary. (Essentially, this is because momentum is transferred across a fluidsolid interface similarly to the way momentum is transferred across a fluidfluid interface.) The result of the noslip condition is that while a flow can appear to be happily streaming past a body, say a wing or a car or something, at the actual solid surface the fluid and the body have the same velocity. One clear manifestation of this is the way you can't clean off a dusty car just by driving fast. In this class we will only be concerned with Newtonian fluids, where the viscous shear stress, , is proportional to the velocity gradient. In Fig. 1, this means that = du dy where is called the coefficient of dynamic viscosity. 1 There are also normal viscous forces, which we'll describe later, but the shear force is the one that's more intuitive. 3 ME 563  Intermediate Fluid Dynamics  Su Lecture 6  Basic viscous flow ideas
Reading: Acheson, chapter 2 (through 2.3). In the last lecture we introduced the concept of viscosity. The most intuitively understandable property of viscous fluids is that they resist shear. For example, given a parallel flow of viscous fluid with a velocity gradient perpendicular to the flow direction, viscous shear forces will be present that act to make the flow speed uniform, by simultaneously accelerating slower fluid and decelerating fast fluid (Fig. 1). Figure 1: Viscous shear forces in a parallel fluid flow. Viscosity is a diffusive property. In particular, owing to viscosity, momentum in a fluid diffuses from regions of high momentum to regions of low momentum. (The rate of momentum loss is equal to a negative force, and the rate of momentum increase is equal to a positive force.) As such, viscous diffusion of momentum is analogous to the diffusion of heat or species concentration. We will also see that viscosity is responsible for the diffusion of vorticity. In the last lecture, viscosity arose in our discussion through the notion of the boundary layer. Specifically, it was pointed out that while the flow around an airfoil could be treated for the most part using ideal fluid theory, the generation of circulation around an airfoil by which the wing generates lift can only be explained by considering the viscosity of the fluid. In the case of the wing, as with many flow problems, the influence of viscosity is confined to a thin boundary layer near the solid surface of the wing. The boundary layer arises because of the noslip condition for viscous flows, which states that at a solid surface, the fluid must have the same velocity as the surface. The result is a boundary layer region in which the flow makes the adjustment from the essentially inviscid main flow to the velocity values at the solid body. If the viscosity is small, this boundary layer can be very thin. Figure 1.11 in the book shows that ideal fluid theory is quite valid for a wing until its angle of attack exceeds a critical value. The failure of the ideal fluid theory arises when the boundary layer ceases to be thin and instead affects a substantial portion of the flow. This phenomenon is known as boundary layer separation, and occurs when the boundary layer experiences a strong increase in pressure along the solid boundary. Boundary layer separation is responsible for significant qualitative differences in flow properties between ideal fluids and those with even very small viscosity. The classic example of drag on a sphere (e.g. why a dimpled golf ball has less drag than a smooth ball) is fundamentally tied up with boundary layer separation. (More on this later.) 1 Equations of motion for viscous fluids We'll have a look at the equations of motion for viscous fluids here, go over some sample problems, and defer the detailed derivation of the equations for later. We will restrict our attention to 1 Newtonian fluids, where the viscous shear stress is directly proportional to the velocity gradient, that is (in the coordinate system of Fig. 1) = du , where is called the coefficient of dynamic viscosity. dy This shear stress has the units of pressure (force/area), so the viscous shear force, F , acting across a surface element in a flow with area A is F = A = du A. dy The equations of motion for a Newtonian fluid with constant density, , and constant viscosity, , are u =0 1 )u =  p+ Du u = + (u Dt t 2 u + g. (1) These are known as the NavierStokes equations, and are the same as the Euler equations we've been discussing except for the viscous term, (/) 2 u. We showed earlier that the Laplacian operator, 2 , could be interpreted in light of our understanding of the divergence, because the Laplacian is just the divergence of a gradient. That is, given a scalar quantity, ,
2 = ( ). To interpret the viscous term in the NavierStokes equations, first observe that the second equation in (1) is really three equations, one for each of the velocity components. Let's consider Cartesian coordinates, so u = (u, v, w), and write the equation for the ucomponent: u + (u t )u =  1 p + x 2 u, (2) where particular care has to be taken with the second term on the left hand side. The gravitational term has been assumed to lie in the zdirection, as is conventional in Cartesian coordinates. Let's now take a (somewhat crude) look at the Laplacian in (2), 2 u = 2 u/x2 + 2 u/y 2 + 2 u/z 2 . As with our discussions of the divergence, start with the familiar infinitesimal cubic volume element V (Fig. 2). We want to calculate the viscous forces on each face. Consider first the Figure 2: Our usual cubic volume element. shear forces. To evaluate these, we need the velocity derivatives u/y and u/z, since the fluid 2 is Newtonian. Specifically, on faces 2 and 5, where the normal vector is parallel to the yaxis, we need u/y, and on faces 3 and 6, where the normal is parallel to the zaxis, we need u/z. As usual, on each of these faces the derivative components of interest will be assumed to be constant. Figure 3 shows us how to interpret the shear terms in the Laplacian, 2 u/y 2 and 2 u/z 2 . Fig. 3a illustrates the variation of the ucomponent of the velocity in the ydirection, which is Figure 3: Interpretation of the viscous shear terms. relevant to faces 2 and 5. Consider face 2, i.e. the y = dy surface of V . The faster fluid at y > dy imposes a force in the positive xdirection on the fluid inside V , represented by Viscous shear force on face 2 = u (y = dy) dx dz. y Meanwhile, on face 5 (y = 0), the slower fluid at y < 0 imposes a force in the negative x direction on the fluid inside V : Viscous shear force on face 5 =  u (y = 0) dx dz. y The net viscous shear force on faces 2 and 5 is thus, with some manipulation, Net viscous shear force on faces 2 and 5 = = = u u (y = dy)  (y = 0) dx dz y y
u y (y = dy)  dy u y (y = 0) dx dy dz (3) 2u dV y 2 Since (2) describes the rate of change of the velocity component u, we need to divide the mass out of (3). The mass contained in V is just mV = dV , so, dividing this out of (3), we get the term 2u y 2 consistent with (2). The variation of the ucomponent of the velocity in the zdirection (Fig. 3b) 2 can be described similarly, and gives rise to the term u . z 2
2 The variation of the ucomponent in the xdirection leads to a similar term, u , but in this x2 case the viscous effect is not a shear effect, but instead a normal effect. This arises because viscous 3 Figure 4: Interpretation of the normal shear term. fluids resist not only resist shearing, but also resist expansion and compression (think of pulling apart or squeezing, say, a blob of caramel). An expansion is represented by (u/x) > 0, while a compression would have u/x < 0 (Fig. 4). It is somewhat unintuitive that compressions and expansions should be described mathematically in a similar fashion to shearing, but this turns out to be the case. That is, u describes a normal viscous force on a face, just as u and u x y z described shear viscous forces. Consider face 1 of our volume V , corresponding to the x = dx surface in Fig. 4. In the expansion case, u/x > 0 across x = dx, and the faster fluid at x > dx will act to pull on, or accelerate, the slower fluid in V . On face 4, in the expansion case, the slow fluid at x < 0 imposes a drag on the fluid in V . The net result is (skipping some steps that are exactly analogous to the development of (3)): Net normal force on faces 1 and 4 = So, the viscous term in (2) breaks down as
Normal force on faces 1,4 Shear force on faces 2,5 Shear force on faces 3,6 2u dV. x2 2 u= 2u x2 + 2u y 2 + 2u z 2 . (4) 4 ME 563  Intermediate Fluid Dynamics  Su Lecture 7  The Reynolds number/some viscous flow examples
Reading: Acheson, 2.22.4. In the last lecture we wrote down the NavierStokes equations of motion for incompressible, viscous fluids: u =0 p )u =  + u + (u t 2 u + g. (1) The first of these equations describes the incompressibility; the remaining equations (one for each spatial dimension) describe the evolution of the velocity field (and are sometimes referred to as the momentum equations). These equations differ from the the Euler equations discussed earlier by the presence of the viscous term 2 u. As we discussed in the last lecture, this term covers both viscous shear forces and viscous normal forces. 1 The Reynolds number and kinematic viscosity UL , The Reynolds number quantifies the importance of viscosity in a given fluid flow. It is defined as Re (2) where U is a characteristic velocity in the flow, L is a characteristic length scale, and and are the density and dynamic viscosity of the fluid, respectively. To see how to interpret the Reynolds number, let's look at where it comes from. We will want to rewrite the momentum equation in (1) in nondimensional form. Define x = i xi , L u = i ui , U t = t , L/U (3) which are the nondimensional forms of these variables. It is important to use the nondimensional form of the operators as well, i.e. =L , L = . t U t 1 p U + 2 L L
(4) Now use (3) and (4) in (1), and we get U 2 u U 2 + (u L t L )u =  2 u + g. Now multiply (L/U 2 ) through on both sides: u + (u t
`inertial term' )u =  L 2 u + 2g UL U 1 2 L =  p + u + 2 g, Re U p + (5) `viscous term' where we have written p = p/(U 2 ), which is the common nondimensionalization for the pressure. 1 If we have chosen the parameters U and L correctly for the problem, then the various nondimensional terms involving u and its derivatives in (5) should be of order 1. In that case, the value of the Reynolds number tells us the relative importance of the inertial and viscous terms O(Re) = inertial term . viscous term Thus, a high Reynolds number means viscous effects are less important, and a low Reynolds number means viscous effects are very important. For ideal fluids, then, it is required that the Reynolds number be very large. However, just because the Reynolds number in a flow is large doesn't mean that we can automatically assume that ideal fluid theory applies. This is because turbulence develops at high Re, and turbulence complicates matters immensely. You will have noticed that in the NavierStokes equations, viscosity appears divided by the density, so we always see /. This is called the kinematic viscosity, /, and is the relevant viscosity in most fluid problems. The kinematic viscosity is antiintuitive in some ways. For example, the dynamic viscosity, , of water is about 55 times that of air, which is probably consistent with your intuition. However, the kinematic viscosity of air is about 15 times that of water, so in that sense air is more viscous than water. 2 Plane Poiseuille flow Figure 1: Plane Poiseuille flow. Consider the system of Fig. 1, consisting of the flow of a viscous, incompressible fluid between two fixed, parallel plates, where the flow is driven by a constant pressure gradient P/x. (This is known as plane Poiseuille flow.) We want to know the velocity field, which in general is u = (u(x, y, z, t), v(x, y, z, t), w(x, y, z, t)). The first thing we assume is that the flow is steady, because the pressure gradient driving the flow is constant. We will also assume that the velocity doesn't vary in the x or zdirections. To justify this, consider that to an observer who's free to move around in the xz plane (for fixed values of y), all points in the plane look the same the pressure gradient is constant, and the boundaries don't vary with x or z. We can also assume that the wcomponent of the velocity is zero, because there is no pressure gradient or gravitational force that could potentially drive the flow in the zdirection. This leads to u = (u(y), v(y), 0) and we haven't needed the NavierStokes equations yet. Now we'll apply the NavierStokes equations. Because the flow is incompressible, u = u v w + + = 0. x y z 2 The first term u/x is zero because u is not a function of x, and the third term w/z is zero because w = 0, so this means that v/y is also zero. We already have v = v(y), so v/y = 0 means that v is constant everywhere. But v = 0 at y = h because the fixed walls are there, so v = 0 everywhere. (Observe that this also means we are ignoring gravity.) Thus, u = (u(y), 0, 0) We are left with the following terms in the NavierStokes equation for the ucomponent of velocity: 0= 1 p 2u 2u + 2 = C + 2 x x x where C is a constant (because the pressure gradient and density are constant). Integrating once, Cy + C1 = u . y To find C1 , note that the u profile needs to be symmetric around y = 0 (why?), so u/y = 0 at y = 0, and C1 = 0. Integrating again, we get 1 2 Cy + C2 = u. 2 The boundary condition on u is that u = 0 at y h. This allows us to evaluate C2 . The final result is u= C 2 (y  h2 ). 2 3 ME 563  Intermediate Fluid Dynamics  Su Lecture 8  More viscous flow examples
Reading: Acheson, 2.32.5. In the last class we illustrated some basic viscous flow ideas by finding the velocity field for the steady flow of viscous, incompressible fluid driven by a constant pressure gradient between two parallel plates (plane Poiseuille flow). In this lecture we will look at a more sophisticated example. The added sophistication doesn't arise particularly in the physics of the problem, but rather in the mathematical tools necessary for solving the differential equations. 1 Flow due to an impulsively moved plane boundary (This example is from the book. We'll work it out in full detail here.) Consider the system in Fig. 1. A viscous, incompressible fluid occupies the region y > 0, which is bounded by a solid wall at y = 0. The wall is initially at rest, but at t > 0 the wall has velocity U in the xdirection. The flow is not pressuredriven. We want to compute the velocity field u(x, t) = (u, v, w) for t > 0. Figure 1: An impulsively moved plane boundary. As usual we want to simplify the problem as much as possible before looking at equations. First of all, the geometry of the problem doesn't vary in the x and zdirections, so we assume that the velocity components are only functions of y. Also, there is nothing driving the flow in the zdirection, so we assume that w = 0. This gives us u = (u(y, t), v(y, t), 0). Because the flow is incompressible, u = u/x + v/y = 0. We have u/x = 0 since u = u(y, t), so v/y = 0. We know v = v(y, t), so v/y = 0 means v is constant everywhere but we have as a boundary condition that v = 0 at y = 0, so finally, we get v = 0 everywhere. Thus u(x, t) = (u(y, t), 0, 0). (This is a plane parallel shear flow.) The NavierStokes momentum equation is: u + (u t )u =  p + 2 (1) u + g. We are only interested in the ucomponent of this equation, which, for the velocity field of (1), works out as u 2u = 2. t y (2) 1 We have eliminated the pressure gradient term because the flow is not pressuredriven. This equation is called the onedimensional diffusion equation, and arises with diffusive properties such as heat. (The thermodynamic analogue of this problem would be if our fluid had initially constant temperature, and suddenly the temperature of the boundary was changed.) The conditions on this problem are: Initial condition: u(y, 0) = 0 for y > 0. Boundary conditions: u(0, t) = U and u(, t) = 0 for t > 0. (3) We are going to look for something called a similarity solution. There is a lot of lore about similarity solutions in the mathematical literature, but unfortunately there seem to be no general rules about when you can and can't look for similarity solutions and even if there were, the math would be way too complicated for us. Instead, we'll rely on reasoning. Let's consider what's happening in this system. Initially the fluid and the wall are at rest (Fig. 1). Suddenly, the wall is set in motion at velocity u = U . Because the fluid is viscous, the fluid in contact with the wall also moves at speed U by the noslip condition. This establishes a velocity difference, or gradient, with the layer of slower fluid immediately above. The resulting viscous shear stress acts to accelerate that layer of slower fluid. As its velocity increases, this layer, in turn, acts to accelerate the slower fluid above it, and on and on. At any time in the process, we get a velocity distribution that looks something like that on the right in Fig. 1. Basically, the moving wall introduces momentum to the system, which gets diffused upward by the fluid viscosity as time proceeds. What can we say about the shape of the velocity profile as time proceeds? One thing we can assume is that U is a scale factor for the velocity field, so a change in U would have a direct, linear effect on the velocity values u(y). For example, if we double U , then we would just double all of the velocity values in the fluid for a given t. Another thing we notice is that the process by which the different `layers' of fluid continually accelerate slower fluid above them has no stopping point, because the fluid domain is infinite in the positive ydirection. So there is no fixed length scale set by the boundary conditions in the ydirection. This suggests that if we look at the velocity profile at successively later times, the shape will be the same but the profile will just be continually stretched in the ydirection. What is the stretching factor for the profile height? The text describes one way to find this `stretching' (or normalizing) factor, . Another way is to use dimensional analysis. The height of the profile can only depend on the time, t, and the viscosity, . The only other parameter in the problem is U , and we've already reasoned that it just scales the velocity values. (Looking at the u profile in Fig. 1, this means that U scales the profile horizontally, and doesn't affect the height in y.) The viscosity, , has units of (length)2 /(time), so the only way to form a variable, , with units of length from and t is to set = t. So we are looking for a velocity field of the form u = U f (), where = y y = . t By the chain rule, we have (using the notation df /d = f ()) u y = f () = f () 1/2 3/2 , t t 2 t u 1 = f () = f () 1/2 1/2 , and y y t 2u 1 = f () . y 2 t 2 Inserting into (2), we get f () y 2 1/2 t3/2 1 = f () , or t (4) 1  f = f . 2 Observe that this is a firstorder differential equation, because one term involves the first derivative and the other term involves the second derivative. So we can define g() = f (), and inserting into (4), we get 1 dg  g = . 2 d We can use the separation of variables to write this as 1 dg  d = . 2 g Integrating both sides, we get 1  2 + C = ln(g), 4 where C is an arbitrary constant. Taking the exponential of both sides gives us g = eC e
2 /4 = Be 2 /4 = df , d (5) where B = eC and we have resubstituted g = df /d. Integrating (5) once, we get f () = A + B 0 es 2 /4 ds. (6) We use the dummy variable of integration s, because we want to evaluate the integral at specific limits so we can impose boundary conditions, and there is no analytical expression for the integral 2 of es /4 . Our boundary conditions were given in (3). Since the similarity function f was defined by u = U f (), our conditions on f are f (0) = 1 and f () = 0. The condition f (0) = 1 is easy, because the integral in (6) goes to zero, so A = 1. To evaluate f at = , we exploit the fact that the definite integral 0 es 2 /4 ds = is found in integral tables. So for f to satisfy the boundary condition at = requires that 1 B =  . So finally, we get as our solution for u u(y, t) = U 1 1 0 es 2 /4 y ds , where = . t (7) The velocity field defined by (7) basically looks like the one depicted in Fig. 1. 3 2 Diffusion of vorticity/flow with circular streamlines We said before that viscosity was responsible for the diffusion of both momentum and vorticity. This can be seen explicitly by looking at the equation for the evolution of the vorticity field, (x, t): + (u t ) = ( )u + 2 , (8) in which the term 2 describes the viscous diffusion of vorticity, analogously to the viscous term in the NavierStokes momentum equation. To illustrate the viscous diffusion of vorticity, it is convenient to consider flows with circular streamlines. These are most naturally described in cylindrical coordinates. In particular, we are interested in velocity fields of the form u = u (r, t) e . ^ (9) The book gives the NavierStokes equations in cylindrical coordinates in 2.4. The one we're interested in is the equation for the time evolution of the u velocity component, which is u + (u t u r u 2 ur 1 p u 2 u + 2 = +  2 r r r r u where u = ur + + uz , r r z 1 2 1 2 and 2 = r + 2 2 + 2. r r r r z ) u + ) u = 0 and
2u , For the velocity field of (9), we then have (u u 1 p = + t r = 2 u r 2 + 1 u r r , which gives us (10) 1 u 2 u u +  2 r 2 r r r . Consider the pressure term. Every other term in (10) is a function of r and t only, so we can write p = P (r, t). Integrating, we get p = P (r, t) + C(r, t) where C(r, t) is a constant of integration with respect to the variable , which means that it can vary with r and t. We can say that P (r, t) has to be zero, because otherwise p() = p( + 2), i.e. p would have multiple values at a given position. Thus p is not a function of , so (10) becomes u = t 1 u 2 u u +  2 2 r r r r . (11) We will use this to explore particular flows with circular streamlines in the next lecture. 4 ME 563  Intermediate Fluid Dynamics  Su Lecture 9  Even more viscous flow examples
Reading: Acheson, 2.42.5. In the last lecture, we began to consider flows with circular streamlines, for which the velocity fields are given by u = u (r, t) e . ^ For such a flow, the NavierStokes equation describing incompressibility, satisfied. The NavierStokes equation for the time evolution of u is u = t 1 u 2 u u +  2 r 2 r r r . u = 0, is obviously (1) Considering flows with circular streamlines is a convenient way to illustrate the viscous diffusion of vorticity. We'll work through one of the book examples, that of the viscous decay of a line vortex. 1 Viscous decay of a line vortex
0 e , ^ 2r Define a line vortex as u= (2) with 0 constant. For a flow with circular streamlines, the vorticity is computed as = u= 1 (ru ) . r r For the line vortex flow of (2), then, the vorticity works out as = 1 (0 /2) = r r undefined if r = 0 0 if r > 0. That is, there is a vorticity singularity at r = 0. If we set this line vortex as the initial condition, we want to know what happens to the flow as time proceeds. Instead of proceeding by looking at the velocity, u , the book chooses to proceed using the circulation: (r, t) = 2ru (r, t), where, obviously, u (r, t) = (r, t) . 2r (3) Expressed as an evolution equation for the circulation, the NavierStokes equation (1) is thus = t 2 1  r 2 r r , (4) 1 using u 1 = , t 2r t u 1 , and =  r 2r r 2r 2 2 u 1 2 1 =  2 + 3. 2 2 r 2r r r r r The initial condition is clear from the definition of the line vortex (2) and the expression for the circulation (3) (r, t = 0) = 0 . We also have the boundary condition that u be finite at r = 0 and t > 0, namely (r = 0, t) = 0 for t > 0. (6) (5) This problem is very similar in form to the one about the impulsively moved plane boundary we talked about in the last lecture. As we did there, we will look for a similarity solution, of the form = 0 f (), where = r/ t. (7) Plugging (7) into (4), we get f = f t r
2 +f 2 1  f , r 2 r r (8) using the notation df /d = f . For our definition of , we have r =  1/2 3/2 , t 2 t 1 = , and r t 2 = 0. r 2 Inserted in (8): f Multiplying through by
r r 2 1/2 t3/2 = f 1 1 1  f . t r t t on both sides: f r2 r = f  f , giving 2t t 2 f = f  f , or 2 1 f  =f . 2 (9) 2 As in the last lecture, we approach (9) first as a firstorder differential equation, by using g() = f (). This gives g By separation of variables, dg = g Integrating both sides, we get ln g = ln  2 + C, 4 1  2 d. 1  2 =g = dg . d where C is an arbitrary constant. Taking the exponential of both sides, g = eC e Integrating again to get f , we have f = 2eC e
2 /4 2 /4 = df . d + A. (10) Now, we need to apply the initial and boundary conditions. By the initial condition (5), at t = 0, r > 0 = , f = 1. Plugged into (10), this gives us f () = 2eC e()
2 /4 + A = 1, so A = 1. The boundary condition (6) requires f = 0 for r = 0 and t > 0, i.e. f ( = 0) = 0. This gives 1 f (0) = 2eC + 1 = 0, i.e. eC = . 2 With these conditions, our solution for the circulation, , is = 0 1  er giving the velocity field as u = 0 2 1  er /4t . 2r
2 /4t , The velocity, u , at different times (different values of t), plotted against r, is shown in Fig. 1. 3 Figure 1: Plot of u vs. r for t = 0, 0.5, 1, 2 and 4. From Batchelor, An Introduction to Fluid Dynamics. 4 ME 563  Intermediate Fluid Dynamics  Su
Lecture 10  D'Alembert's paradox / the velocity potential and stream function
Reading: Acheson, 4.12, 4.13, 4.1, 4.2. We will set aside the book's chapter on waves (some of which we'll touch on later) and spend a small amount of time on the chapter about classical airfoil theory. `Classical airfoil theory' amounts to the detailed use of ideal fluid ideas. We're interested only in taking a couple of simple results from this chapter.
1 D'Alembert's paradox Earlier we pointed out that ideal fluid theory predicts lift well for wings until the angle of attack exceeds some threshold. There is another instructive example where ideal fluid theory fails to match with physical observations  expressed as D'Alembert's paradox, which states that the steady, uniform flow of an ideal fluid past a fixed body imposes no drag on the body. This can be viewed as a statement that ideal fluids are frictionless. This has obvious implications for viscous fluids.... We will need to express the notion of the momentum balance for control volumes. The book derives this in an interesting way (which gives a good lesson in the use of index notation). For an ideal fluid, we can write Euler's equation for momentum as: (u )u =  P (1) where we can write (labeling our Cartesian directions with the subscripts 1, 2 and 3) u = (u1 , u2 , u3 ) = ui ei ^ = e1 ^ + e2 ^ + e3 ^ = ei ^ x1 x2 x3 xi u = u1 + u2 + u3 = ui x1 x2 x3 xi making use of the implied summation convention for repeated indices. We can then write (1) as uj p ui =  xj xi (2) which is strictly speaking the equation for the xi component of (1). Now we want to integrate this over a fixed volume V enclosed by a surface S with local outward normal vector n. For the left hand side, we will use the following: ^ uj ui (uj ui ) = uj + ui . xj xj xj But, by the implied summation convention, the second term on the right hand side is ui because uj = ui xj u1 u2 u3 + + x1 x2 x3 = ui ( u) = 0 u for ideal fluids. So we end up with ui (uj ui ) = uj . xj xj 1 (3) Now we want to integrate (3) over V , i.e. we want to find (using =const.) But observe that we can write (uj ui ) = xj (think about this), and then we have for (4) (uj ui ) = xj (ui u) dV = (ui u) n dS = ^ ui uj nj dS (5) (ui u) (uj ui ). xj (4) V V V S S where we have applied the implied summation to the u n dot product. ^ For the right hand side of (2), our integral is  p dV =  xi p ni dS =  p^ dS n V S S (using an identity in the text derived from the divergence theorem). The result of integrating (2) is thus 
S p^ dS = n S u(u n) dS ^ which is exactly what you'd find from a control volume analysis. 2 The velocity potential and streamfunction Section 4.2 in the text goes over the justification for the use of the velocity potential for the case of irrotational flow (i.e. vorticity, = u, = 0) u= , where is a scalar function in more detail than we did in the first homework. Note especially that in going from Eq. 4.2 to 4.3, Acheson makes use of Appendix A.2. The streamfunction is another useful formulation for the velocity field, and is valid for incompressible, twodimensional flow, where u = 0. The fundamental idea is to use the vector identity that ( f) = 0 for any vector f i.e. the divergence of the curl of any vector is zero. This justifies Eq. 4.8 in the text. Important note: the velocity potential and the streamfunction are not limited to ideal fluids. 2 ME 563  Intermediate Fluid Dynamics  Su
Lecture 11  Vortex motion: induced velocity, Kelvin's circulation theorem
Reading: Acheson, 5.15.3. It is sometimes convenient to formulate flow problems in terms of the vorticity rather than the velocity. Not only are a lot of interesting flow systems characterized by rotation, but the equations for the time evolution of the vorticity can be simpler because they don't involve pressure. If you're doing computer simulations, for example, not having to compute the pressure might simplify your task considerably. 1 Velocity field induced by a vortex Of course, computing a flow field in terms of the vorticity makes little sense if we can't get back the velocity from the vorticity field. The way we do this is through the BiotSavart law, which tells us the velocity field induced at a point x by a given vorticity field, as u(x) = 1 4 (x ) s dV. s3 (1) In this equation, x is the point at which we're computing the velocity field, the integrand tells us the velocity induced at x by the vorticity in a volume element dV located at x , s is the vector (with magnitude s) connecting the point x to x, and the integral is taken over all space (specifically, all space in which the vorticity is nonzero). The process by which a vorticity field induces a velocity field is exactly analogous to that by which an electric current induces a magnetic field. In the form of (1), the BiotSavart law is written for vorticity fields distributed over volumes, since the integrand includes the volume element dV . We are also interested in describing vorticity fields that are concentrated into lines (in the electromagnetic analogy, this corresponds to line currents). For such vorticity fields, the BiotSavart law becomes u(x) = 1 4 dl(x ) s s3 (2) where dl(x ) is a length element along a vortex line located at x , and the integral is taken over all vortex lines that might be present. The vortex line strength is denoted by , not by the vorticity, , as appeared in (1). There are two immediate explanations for this. The first is that the vorticity is undefined on vortex lines, as we recall from Chapter 2. The second is that the vorticity has the wrong units for (2). We will illustrate the use of (2) by computing the velocity field induced by a straight, infinite vortex line of uniform strength aligned with the positive zaxis (Fig. 1). Consider a point located a distance r from the vortex line. Without loss of generality, we can assume that the point lies at z = 0. The line element dl points in the positive zdirection, so dl = dz ez . ^ The crossproduct dl s therefore works out as dl s = dz s sin e , ^ where e points along the dotted loop in the figure in the sense indicated, by the righthand rule. ^ The magnitude s is given by s= z2 + r2, 1 Figure 1: Velocity induced by a vortex line. since the line element dl in question is located at z, so sin can be written sin = The relation (2) thus becomes e ^ u(r) = 4 = r r . = 2 + r2 s z r [r 2 + z 2 ]3/2 1 dz z= e ^ 4r 2 z= [1 + (z/r)2 ]3/2 dz (3) where we can take the r out of the integral since the integral is over z. To solve the integral, we make the substitution z/r = tan and use the following relations: dz = r sec2 d 1 + tan2 = sec2 z = = . 2 Equation (3) thus becomes
/2 /2 e ^ u(r) = 4r =/2 1 e ^ d = sec 4r cos d =
=/2 e ^ /2 [sin ]/2 = e . ^ 4r 2r (4) In Chapter 2 of the text (p.46) we were given the velocity field of a line vortex as u= 0 e . ^ 2r (5) By comparison of (4) and (5), then, we see that the measure of the strength of a line vortex, which we called in (2), is in fact the circulation 0 induced by the vortex. 2 2 Theorems of vortex motion: Kelvin's circulation theorem Kelvin's circulation theorem is expressed as follows: given an ideal fluid subject to conservative body forces (e.g. gravity), and with C(t) being a closed circuit consisting of the same fluid particles as time proceeds, the circulation =
C(t) u dx is not a function of time. The proof of this is outlined in the text (in 5.1, and also in Exercise 5.2, which is explained in the back of the book). A key point is that there is no condition in the theorem requiring the surface enclosed by C(t) to lie entirely in the fluid. In particular, it is permissible for C(t) to enclose a wing, for example. This completes the explanation of the generation of lift on a wing started in Chapter 1 (explained in 5.1). 3 ME 563  Intermediate Fluid Dynamics  Su
Lecture 11a  Proof of Kelvin's circulation theorem
Kelvin's circulation theorem. Consider an ideal fluid moving in the presence of a conservative body force (e.g. gravity). Let C(t) be a closed circuit consisting of the same fluid particles as time proceeds. Then the circulation =
C(t) u dx is independent of time. Proof. We will need first to prove the following relation: d dt u dx = Du dx. Dt C(t) C(t) We start by representing the points on the curve C(t) by x(s, t), where s [0, 1]. That is, we are parameterizing C with the variable s, so that every fluid particle on the curve C gets assigned a unique value of s between 0 and 1. An example of how this parameterization works is given on the left in Fig. 1. To express dx, the vector line element along the curve C, in terms of Figure 1: Sample parameterization of a curve C in terms of s (left); determination of a line element dx along C (right). this parameterization, look at the other cartoon in Fig. 1. Consider two points on C, x(s, t) and x(s + ds, t). Let dx be the vector that goes from x(s, t) to x(s + ds, t), i.e. x(s, t) + dx = x(s + ds, t). Thus, dx = x(s + ds, t)  x(s, t). Clearly, as ds becomes small, then dx becomes tangent to the curve C, and the above relation becomes dx = x(s + ds, t)  x(s, t) x ds = ds ds s as ds 0 (1) where the quotient is the partial derivative x because t is held constant. s This expression for dx allows us to write d dt d u dx = dt C(t)
1 u
s=0 x ds. s (2) The limits s = 0 and 1 on the second integral reflect the parameterization, whereby going from s = 0 to 1 takes us once around the curve C(t). Since s and t are independent variables in the 1 parametric representation of C, we can put the time derivative d/dt inside the s integral as long as we write it as the partial derivative /t, with s held constant. Thus d dt
1 s=0 x u ds = s 1 0 t u x s ds. (3) By the product rule, we have for the term in the second integrand t u x s = u t x s +u t x s . s t s t The term (u/t)s is the time derivative of the velocity for a particular fluid element (fixed value of s), which is by definition the same as the material derivative Du/Dt. For the second term on the right hand side, we can interchange the order of partial differentiation, so u t x s =u s x t =u u s =
t s t t s s 1 uu . 2 This allows us to rewrite (3) as
1 0 t x u s 1 ds =
0 Du x + Dt s s 1 uu 2 ds. (4) Note that the second term in the integrand just evaluates as
1 0 s 1 1 u u ds = uu 2 2 1 = 0,
0 because s = 0 and s = 1 represent the same point in space, and u is singlevalued in space (meaning it has only one value at any point). So, combining this result with (4), (3) and (2), we have d dt
1 C(t) u dx =
0 Du x ds = Dt s C(t) Du dx, Dt (5) where we have used (1) in the last step. Now we go back to Euler's equation for the velocity field (since we are dealing with an ideal fluid) Du 1 = p+g = Dt p + where we have written g =  since gravity is conservative, and we have moved inside the derivative because it's a constant for ideal fluids. Then, (5) becomes d = dt p + dx. (6) C(t) To intepret this integral, recall some of the properties of gradients. Let be a scalar function. If ^ is a unit vector, then we had r ^ = r 2 d , dr which means that the gradient of dotted into the unit vector ^ is equal to the derivative of in r the ^ direction. Now consider the vector line element dx, which we will write as dx = dx^ . So, dx r x ^ is the magnitude of dx, and x is its direction. We then have ^ dx = ( x)dx = d dx = d, dx ^ where this d is the change in on moving a distance dx in the x direction. Now we look again at (6). Using the results for the gradient above, we can write
at end of curve C(t) at start of curve C(t)  C(t) p + dx =  C(t) d p + = p +  p + = 0, where the result is zero because C(t) is a closed curve (so its start and end are the same point), and p, and are singlevalued functions. Thus d = 0, dt so the circulation is not a function of time. Q.E.D. 3 ME 563  Intermediate Fluid Dynamics  Su
Lecture 12  Vortex motion: vortex lines, Helmholtz's vortex theorems
Reading: Acheson, 5.2, 5.3. A vortex line is a curve which is parallel to the instantaneous vorticity vector = u at all of its points. A vortex line thus relates to the vorticity in the same way that a streamline relates to the velocity. So, if we define a vortex line as x = (x(s), y(s), z(s)), and the vorticity as = (x , y , z ), then the following relations hold: dx/ds dy/ds dz/ds = = x y z at all points on a vortex line. The slope of a vortex line, as projected into the x  y plane, is then dy/dx = y /x , with similar results for the y  z and x  z planes. We introduced the concept of vortex tubes in the second homework assignment. The book defines a vortex tube as the set of vortex lines that pass through a simple closed curve in space. (By simple closed curve we basically mean that the curve is connected i.e. one piece and doesn't cross itself.) One Helmholtz vortex theorem deals with vortex lines; the other, with vortex tubes. Both are formulated for ideal fluids in the presence of conservative body forces (like gravity), which are the same conditions for which we formulated Kelvin's circulation theorem. 1 The Helmholtz vortex theorems Helmholtz's first vortex theorem. The fluid elements that lie on a vortex line at some instant continue to lie on a vortex line, that is, vortex lines `move with the fluid.' To visualize this, suppose we mark a vortex line at some time by, say, coloring the specific fluid particles that make up the line. Then at all later times, the marked particles will continue to constitute a vortex line. By this theorem, we can also say that vortex tubes move with the fluid. The proof in the book proceeds by defining a vortex surface as a surface that is everywhere tangent to the vorticity vector . By applying Kelvin's circulation theorem and Stokes's theorem, you can show that the fluid particles that make up the vortex surface at any point in time will continue to make up a vortex surface for all times, that is, the vorticity will always be tangent to that set of fluid particles. Now, pick one vortex line out of that vortex surface, and construct another, separate surface that's tangent to the vorticity everywhere and includes that vortex line. Then the vorticity will always be tangent to those fluid particles also. The intersection of the two different surfaces thus remains a vortex line. Because the two surfaces are defined in terms of fluid particles and thus move with the flow, this means that the vortex line defined by the intersection of the surfaces moves with the flow. Helmholtz's second vortex theorem. The quantity =
S ^ n dS is the same for all crosssectional surfaces S cut through a vortex tube. Also, is not a function of time. 1 The proof of this theorem is very similar to the argument we applied in the virst homework assignment to show that vortex tubes can't end in a fluid. The added step is to equate the vorticity ^ `flux' through a crosssection of the tube ( S n dS) to the circulation around the crosssection using Stokes's theorem. Because vortex lines move with the fluid (by the above theorem), any curve enclosing a vortex tube also moves with the fluid. Then by Kelvin's circulation theorem, the circulation around such a curve remains constant. The interesting interpretation of this theorem is in terms of thin vortex tubes where the vorticity is essentially constant across any crosssection. In that case, the theorem tells us that if the crosssectional area decreases, the vorticity has to increase. Since the fluid is incompressible, the crosssectional area decreases if we stretch a particular set of fluid elements lengthwise (recall that vortex tubes consist of a fixed set of fluid elements, from the first theorem). Thus, stretching a vortex tube lengthwise results in increased vorticity. This has interesting interpretations in analogies with other physical examples  think especially of the conservation of angular momentum. 2 ME 563  Intermediate Fluid Dynamics  Su
Lecture 13  Vortex motion: interpretation of Helmholtz's vortex theorems
Reading: Acheson, 5.35.5. In the last lecture we covered Helmholtz's (two) vortex theorems, which apply to ideal fluids subject to conservative body forces. The theorems state: The fluid elements that lie on a vortex line at some point in time continue to define a vortex line. This means that vortex lines move with the fluid. Given a vortex tube, the circulation, =
S ^ n dS, (using Stokes's theorem) is the same for all crosssections S cut through the vortex tube, and is the same for all times. The second of these two theorems in particular lends itself to ready interpretation in terms of physical intuition. 1 Thin vortex tubes If a vortex tube is thin enough, then the vorticity will be virtually constant across the tube's crosssection. We will consider the special case of a thin, cylindrical vortex tube, defined as the set of vortex lines that pass through a small circle drawn in the flow, with the additional assumption that the vortex lines are parallel (Fig. 1) Consider a finite segment cut from such a vortex tube, with Figure 1: A thin, cylindrical vortex tube. length l. Label the circular crosssections that form the ends of this finite segment as S1 and S2 , and denote the area of these faces as S. We will focus our attention on face S1 . Because its area, S, is small, we will assume that the vorticity on any point on that face has the constant value 1 . Then the circulation around the circle C1 that encloses the surface S1 is 1 =
C1 u dx = (Stokes's theorem) S1 ^ n dS = S. (1) Now consider what happens if we change the length of the finite segment of the vortex tube. Specifically, we will consider the set of fluid elements initially contained in the segment of length l cut from the vortex tube. Because ideal fluids are incompressible, the volume occupied by that set of fluid elements has to be constant. Assume that the segment we're considering remains cylindrical even as its length changes. The volume of the fluid element is lS, so as depicted in Fig. 2, if we 1 Figure 2: Stretching/contracting a vortex tube. change its length from l to l1 > l, then the crosssectional area has to decrease (proportionally), giving S1 < S. Similarly, if l2 < l, then S2 > S (proportionally). But by (1), S is constant, so: If we stretch a thin segment of vortex tube (l1 > l), its vorticity increases (proportionally), giving 1 > . If we contract a thin segment of vortex tube (l2 < l), its vorticity decreases (proportionally), 2 < . 1.1 Tornadoes Figure 3: Tornado formation. This vortex stretching process is what happens in a type of tornado (the nonsupercell variety). Initially, suppose there's some small vorticity near the ground, distributed over a wide area, perhaps in the form of a large, slow circulation caused by the presence of hills or mountains or something. If that slow rotation meets with strong updrafts, maybe due to thermal convection, then the effect is to stretch the vortex tube defined by that circulation. As the tube gets stretched to the characteristic thin tornado profile, the vorticity increases, accounting for the strong rotation of the tornado. As Acheson points out, Helmholtz's first theorem is then illustrated as the tornado subsequently moves to follow the pattern of clouds overhead. 1.2 Relation to the vorticity equation We can also get at Helmholtz's second theorem by looking at the Euler vorticity equation: D = ( Dt 2 )u. Suppose that the vortex lines are almost perfectly aligned with the zaxis, so (0, 0, z ). In that case, z z (2) and we can thus write for the zcomponent of the equation Dz w z Dt z (3) where w is the velocity component in the zdirection. Now consider what it means for w/z to be positive or negative. If w/z is positive, then fluid elements located at higher z also move faster in the zdirection, i.e. there is stretching in the zdirection. If w/z is negative, then fluid elements located at higher z are moving more slowly in the zdirection, i.e. there is contraction in the zdirection. Thus by the above expression, stretching in the zdirection increases the zcomponent of vorticity, since Dz /Dt > 0, while contracting in the zdirection decreases z , since Dz /Dt < 0 in that case. (Question: do you see a possible inconsistency between (2) and (3)?) 1.3 Angular momentum analogy From basic physics you already know some principles that should give you some physical intuition for Helmholtz's second theorem. In particular, we know that angular momentum of a body is conserved in the absence of torques, where angular momentum is the moment of inertia multiplied by the angular velocity (which, like the vorticity, has units of (time)1 ). You will remember that the moment of inertia of an object is related to its radial spread from the axis of roation. That is, for a given mass, an object has a higher moment of inertia if it extends further from the axis of rotation. The angular velocity, by virtue of its units, describes a frequency of rotation (and thus does not relate directly to a linear velocity). Figure 4: Effect of changing the moment of inertia. As an illustration of this, consider Fig. 4. This is supposed to represent a table with a small hole cut in the center. A small weight attached to a string sits on the table, and the other end of the string reaches through the hole. Suppose we set the weight into motion so it swings around at the end of the string. If we pull on the string from below the table, we decrease the moment of inertia of the weight by decreasing its radial distance from the axis of rotation. Then, by conservation of angular momentum, the frequency of its rotation around the hole in the table must increase. This essentially describes what happens when you stretch or contract a segment vortex tube. The mass of a set of fluid elements that make up the segment is constant, but by stretching you decrease the radius and thus decrease the moment of inertia of the segment. It compensates by increasing its rotation rate (i.e. its vorticity). This is the gist of Helmholtz's second theorem. 3 1.4 Something to think about How does this relate to the vortex formed by draining a tub? What's the approximate angular velocity of the tub vortex? If we treat the water in the tub as ideal, where would the vorticity in the tub vortex come from (look at 5.2 in the book)? What about the venerable explanation that the Coriolis force causes tub vortices to rotate in different directions in the Northern and Southern Hemispheres? The Coriolis force is linearly related to the Earth's rotation rate. How does this compare with the tub vortex's rotation rate? Is it reasonable for the Coriolis force to affect the rotation in the tub vortex? 1.5 A cautionary note In discussing the Helmholtz theorems we have made a number of simplifying assumptions. The theorems themselves assume fluids are ideal; and we have also made things easier on ourselves by talking about `thin' vortex tubes, that are aligned with the zaxis, etc. These conditions won't necessarily apply with viscous fluids, but the qualitative conclusions are pretty general. Just remember what assumptions we've made and why they might fail. An example: what fluid property might make it problematic to apply conservation of angular momentum, and why? 4 ME 563  Intermediate Fluid Dynamics  Su
Lecture 14  Vortex motion: axisymmetric flow, spherical coordinates
Reading: Acheson, 5.35.5. We now want to look at some problems that represent axisymmetric flow. This is analogous to twodimensional flow in Cartesian coordinates, in that the flow has nonzero velocity components in two directions, and these components are only allowed to vary in those two directions. The difference is that with axisymmetric flow, the coordinates of interest are referred to a single axis or pole for example, cylindrical coordinates and the coordinates that are allowed to vary are the coordinate that measures distance along the polar axis and the coordinate that measures distance away from the polar axis. Mathematically, the velocity field for axisymmetric flow is expressed in cylindrical coordinates as: u = (ur (r, z, t), 0, uz (r, z, t)). (1) Because u = 0, the streamlines of an axisymmetric flow will be confined to any plane that includes the zaxis, and because u is not a function of , the streamlines are identical for all such planes (Fig. 1). This explains the nomenclature `axisymmetric.' Figure 1: Symmetry planes for axisymmetric flow. The vorticity field for an axisymmetric flow is given by u= 1 r r ur er ^ = ru r ur r^ e ez ^ z rz 1 = r = er ^ r^ e 0 0 ez ^ z rz ur uz  e . ^ z r So, for axisymmetric flow, the vorticity vectors are aligned only in the direction, which means that vortex lines are circular and centered on the zaxis. Vortex rings, commonly seen in the form of smoke rings, represent an axisymmetric flow where the vorticity is concentrated in a thin tube. (This is treated in detail in the text.) 1 1 Spherical coordinates Axisymmetric flow can also, and sometimes more conveniently, be described in spherical coordinates. As shown in Fig. 2, for an arbitrary point in space, the spherical coordinate r measures the distance of the point from the origin, the coordinate measures the angle between the Cartesian zaxis and the r vector, and the coordinate is the angle between the xaxis and the projection of r into the x  y plane. Remember that the cylindrical coordinate system (r,,z) is so named because if we Figure 2: The spherical coordinate system. let the coordinate r vary from 0 to R, let z vary from 0 to h, and let vary over the full range 0 to 2, we describe a cylinder of height h and crosssectional radius R. Similarly, for the spherical coordinate system (r,,), if we let r vary from 0 to R, let vary over its full range 0 to 2, and let vary over its full range 0 to (why not 2?), we describe a sphere of radius R. It can be seen that the spherical coordinate is identical to the cylindrical coordinate , which means that the cylindrical coordinates r and z can be transformed into the spherical coordinates r and . This means, then, that the cylindrical coordinate description of axisymmetric flow (1) can be written in spherical coordinates as: u = (ur (r, , t), u (r, , t), u = 0). (2) 2 The Stokes stream function For twodimensional, incompressible flow in Cartesian coordinates, we learned in Chapter 4 how to write the velocity in terms of a stream function that is constant along streamlines. Since axisymmetric flow is just twodimensional flow in cylindrical or spherical coordinates, we can reasonably expect to be able to write a stream function for axisymmetric flow too. Recall that the stream function representation arises because for incompressible flow because for such flows, u = 0. One of our vector identities tells us that the divergence of any curl is zero, or ( f ) = 0. This suggests that we can write u = f . Such a representation is especially convenient for twodimensional flows, because then the vector f only has a component in the third dimension, the one that has no velocity component, and also the direction in which the velocity terms don't 2 change. For axisymmetric flow expressed in cylindrical coordinates, u = ur er + uz ez , so we can ^ ^ write f = ^ , where . Then e u= 1 (r) er + ^ ez =  er + ^ ^ z r r z + r r ez . ^ (3) In Cartesian coordinates, this was the stream function, so called because was constant along streamlines. It turns out that is not constant along streamlines for the axisymmetric representation (3). In particular (u ) = ur + uz = + r z z r + r r = ur . z r Notice, however, that if we write = /r, so that the velocity field is written u= which gives (u ) = ur 1 1 + uz = + = 0, r z r z r r r z (/r) 1 1 1 er + ^ ez =  ^ er + ^ ez , ^ z r r r z r r then is constant along streamlines. This is called the Stokes stream function. To express these ideas in spherical coordinates, note that r in cylindrical coordinates is equal to r sin in spherical coordinates, and the direction in cylindrical coordinates is the same as the direction in sphericals. Thus, u= (^ ) = e e ^ r sin = 1 r 2 sin 1 er  ^ e . ^ r sin r (4) 3 Irrotational flow past a sphere Consider a sphere of radius a placed in a uniform flow (of velocity U , in the zdirection, as in Fig. 3) of ideal fluid, and realize that the flow is irrotational infinitely far upstream of the sphere. (The Figure 3: A sphere in uniform flow. flow is uniform there, and a uniform flow has no gradients, so the vorticity is zero.) From 5.2, all of the flow will thus be irrotational. For axisymmetric flow, the vorticity in spherical coordinates will only have a component in the direction, = e , which is written ^ = 1 (ru ) 1 ur  , r r r 3 which we can write using (4) in terms of the Stokes stream function, , as =  1 2 sin + 2 r sin r 2 r 1 sin . Specifying the (irrotational) flow for this problem thus comes down to solving 2 sin + 2 r 2 r 1 sin =0 (5) for the region outside the sphere, i.e. r a, subject to boundary conditions. The first condition applies at the surface: = 0 on r = a. (6) In particular, setting constant on the surface of the sphere ensures that the surface represents streamlines, i.e. that the velocity field at the surface is parallel to the surface. The choice of = 0 as opposed to any other constant is arbitrary. The other condition pertains for r , or far from the sphere, where the velocity has to be uniform in the zdirection with value U . In spherical coordinates, that means (referring to Fig. 3) ur = U cos u = U sin . To express these as conditions on the stream function, , defined by (4), we have from the rcomponent: ur = U cos = which gives us = U r 2 sin cos and, integrating over , d = U r 2 1 sin cos d = U r 2 sin2 + c1 (r). 2 (7) 1 r 2 sin Similarly, we have from the component: u = U sin =  which gives us (skipping to the integration) dr = U sin2 r 1 r dr = U sin2 r 2 + c2 (). 2 (8) 1 r sin r Comparing (7) and (8), we see that c1 and c2 can only be a constant that isn't a function of r or . We can, arbitrarily, set that constant to zero since the stream function is mainly of interest in terms of its derivatives (as in (4)). Thus, we are left with 1 = U r 2 sin2 . as r . 2 4 (9) This suggests a solution of the form = f (r) sin2 . (Never mind how we can assume this focus on the result.) Plugging this into (5), we get f 2 f = 0, r2 which suggests that f is a sum of terms of the form r n . We can show that these terms satisfy the differential equation only if n = 2 or n = 1. The general solution of the equation is a linear combination of these terms (again, don't sweat the details), or f = Ar 2 + B r which, applying our boundary conditions (6) and (9), gives = U 2 r2  a3 r sin2 . Applying (4), we get for the velocity components 1 a3 = U 1  3 cos r 2 sin r 1 a3 u =  = U sin 1 + 3 r sin r 2r ur = . (10) Let's look at what we have going on on the surface of the sphere. Figure 4 shows how the ur and u components look there. From (10), we have at r = a Figure 4: The velocity components on the surface of the sphere. ur (r = a) = 0 3 u (r = a) =  U sin , 2 so there is no normal component of velocity on the surface, as desired, but there is a tangential component, which is allowed because the fluid is inviscid. Consider a streamline that strikes the sphere at r = a, , or just off the zaxis. Initially, sin 0, so the tangential component u is small (and points in the upward direction along the surface, as shown in Fig. 5). At the equator of the sphere, = /2, u is maximized, and as we move along the surface to the top of the sphere, u drops sharply as sin 0. By Bernoulli's theorem (either one), this sharp drop in velocity as we move along the surface from the equator of the sphere to the top means that the pressure rises sharply along that streamline. This will become significant in terms of boundary layers, as we will see. 5 Figure 5: The velocity on a streamline at the surface of the sphere. 6 ME 563  Intermediate Fluid Dynamics  Su
Lecture 15  Vortex motion: the von Krmn vortex street, and some complex a a analysis
Reading: Acheson, 4.3, 5.65.8. The phenomenon of the von Krmn vortex street is one of the most interesting and bizarre in a a fluid mechanics. And as a bonus, it has great practical significance. The von Krmn vortex street is a particular case of the wake flow behind a circular cylinder a a whose long axis is perpendicular to a uniform main flow. The fascinating thing is that even if the main flow is steady (i.e. doesn't vary with time), and the cylinder is held fixed, an unsteady pattern of vortices will be shed behind the cylinder. The practical importance of this phenomenon is that the unsteady velocity field implies an unsteady pressure field, so that the cylinder experiences fluctuating forces. In a real flow, the cylinder in question might be a bridge pier in water, or a support cable for a television antenna, or a flagpole, where the fluctuating forces can have implications for structural stability, for example. 1 Some basic complex analysis To treat this problem analytically we will need some basic tools from complex analysis. Define a general complex number as z = x + iy, where i = 1 and x and y are real numbers. Then we refer to x as the real part of z, x = Re z, and to y as the imaginary part of z, y = Im z. We choose to call the real and imaginary parts x and y because the geometric interpretation of the complex number z is as a point in a plane, where the real part of z is the horizontal axis (hence x) and the imaginary part is the the vertical axis (hence y). There is also a convenient polar y z=x0+iy0=rei y0 r x0 x Figure 1: The complex plane. representation of z, as seen in Fig. 1. If we denote the length of z as r, and the angle it makes with the xaxis as , then x = r cos and y = r sin , so z = r(cos + i sin ) = rei , where we have used the expression ei = cos + i sin . Addition and multiplication for complex numbers looks like: (x1 + iy1 ) + (x2 + iy2 ) = (x1 + x2 ) + i(y1 + y2 ) (x1 + iy1 )(x2 + iy2 ) = (x1 x2  y1 y2 ) + i(x1 y2 + x2 y1 ) as you might expect. Note that multiplication (and division) are a little cleaner if you use the polar representation z = rei . 1 We can, of course, have functions of the complex variable z, which are written w = f (z) = g + ih, with g and h being realvalued functions. We call g the real part of the function w and h the imaginary part. We can take derivatives of w exactly the same way we do for functions of real variables (well, almost exactly, but the distinctions don't really matter for us). So if w = z n , for example, dw d n = z = nz n1 . dz dz There is one theorem that will be particularly useful to us. It states that if a function w = f (z) = g + ih is differentiable at some point z0 , then g h = x y and g h = , y x dw g h = +i . dz x x and (1) This can be verified easily for simple functions, say w = z 2 . 2 The complex potential
u= Why is this interesting? Remember that we earlier said that for irrotational flow, we could write where was called the velocity potential. Thus, in the twodimensional case, u= ,v = x y (for 2D, irrotational flow). (2) Also, we said that for twodimensional, incompressible flow, we could write u= ,v =  y x (for 2D, incompressible flow), (3) where was called the stream function. Now observe that (2) and (3) together look like (1), if we let the velocity potential be the function g and the stream function be h. For irrotational, incompressible twodimensional flow, we can then write the complex potential, w, as w(z) = + i, (4) where the complex variable z = x + iy, using the real spatial variables x and y. By (1), and using (2) and (3), we then have dw = +i = u  iv, dz x x (5) which tells us how to get the velocity components back from the complex potential in Cartesian coordinates. If we're using polar coordinates, then 1 = r r 1 u = = . r r ur = 2 (6) The text works through a couple of simple examples of the use of the complex potential. The one we'll use is the complex potential representation for the velocity field of a line vortex. For a line vortex located at the origin, u= e . ^ 2r The complex potential for this is i w =  ln(z). 2 If the line vortex is located at some point z = z0 = x0 + iy0 , then the complex potential is i w =  ln(z  z0 ). (7) 2 3 The von Krmn vortex street a a The text has a good discussion of the development of the von Krmn vortex street which I won't a a bother to duplicate. Then, Acheson models the vortex street (for a lefttoright mean flow) mathematically as an infinite set of line vortices, arranged in two rows. In the lower row, the vortices are a b Figure 2: A model of the von Krmn vortex street. The mean flow is from left to right. a a at z = na, where n is any integer, and have strength , corresponding to a counterclockwise rotation. In the upper row, the vortices are at z = (n + 1/2)a + ib and have strength , corresponding to a clockwise rotation. There are a couple of useful geometric observations. First, we point out that the velocity of any single line vortex will be that induced by every vortex other than itself. Suppose we consider a vortex on either row. By symmetry, the induced velocities of the vortices in the same row cancel in pairs, i.e. the induced velocity of the vortex located a distance ma away cancels with that located a distance ma away. Similarly, the ycomponents of the induced velocities of the vortices in the other row cancel symmetrically, but the xcomponents of the induced velocities add up to give a leftward velocity. Thus, the vortex street moves to the left under the action of its own vorticity. The text works through the math of determining the velocity of one of the vortices by adding the respective complex potentials together for the vortices in the other row. For a vortex on the top row, the complex potential of all of the vortices on the bottom row is i z w= ln sin + const. 2 z The resulting speed of the vortex street is V = tanh 2a b a . In the frame of reference where the cylinder is moving in a fixed stream, the vortex street follows the cylinder but gradually falls behind. 3 ME 563  Intermediate Fluid Dynamics  Su
Lecture 16  Vortex motion: effects of viscosity
Reading: Acheson, 5.9, 5.10. In Chapter 2 we looked at the vorticity equation for viscous fluids, which was + (u t ) = ( )u + 2 . (1) The second term on the left side of (1) is clearly recognizable as the advective term in the material time derivative, and describes the change in the vorticity of a fluid element as it's moved around by the velocity field. The first term on the right side of (1) describes the increase (or decrease) of vorticity due to stretching (or compression) of vortex lines. These terms are the mathematical representation of Helmholtz's first and second vortex theorems. The last term in (1) is the term describing the viscous diffusion of vorticity. For ideal fluids, this term isn't present. The Helmholtz theorems, and Kelvin's circulation theorem, were formulated for ideal fluids, so the viscous term somehow makes these theorems invalid. We discussed why this is in the context of Helmholtz's second theorem. That theorem led to the conclusion that stretching vortex tubes increases their vorticity proportionally, while compressing them decreases the vorticity proportionally. We pointed out that this is analogous to the conservation of angular momentum. Angular momentum is the moment of inertia times the angular velocity, where angular velocity is related to vorticity. By incompressibility, stretching a fluid element makes it thinner, which reduces its moment of inertia. To conserve the angular momentum, the vorticity increases proportionally to the decrease in moment of inertia. For a viscous fluid, however, we can't say that the angular momentum of a vortex tube is conserved, because vorticity can diffuse due to vorticity and thus be transferred out of or into a vortex tube. So Helmholtz's second theorem no longer holds. This viscous diffusion is like a frictional loss. Similar arguments explain the failure of the first Helmholtz theorem and Kelvin's theorem for viscous fluids. The book talks about the interesting example of the Burgers vortex. In Chapter 2 we worked the example of the viscous diffusion of a vortex line. The Burgers vortex is the superposition of a vortex line with a radially inward flow that is, by incompressibility, also extensional in the zdirection. The velocity field is given by: 1 ur =  r 2 uz = z 2 u = 1  er /4 . 2r The corresponding vorticity field is given by = r2 /4 ez . ^ e 4 (2) 1 The PrandtlBatchelor theorem Back in Chapter 1, we saw that for 2D, steady flow of an ideal fluid, the vorticity, , was constant along streamlines. We got this result by looking at the vorticity equation for ideal fluids. Another way is to remember that for 2D, incompressible flow, the velocity field can be written u= , y v= , x (3) 1 where , the stream function, is constant along streamlines. Since the vorticity is z = v/x  u/y, (3) means we can write = z (), that is, the vorticity is a function of the stream function value, which means the vorticity is constant on a given streamline. For the ubiquitous case of 2D flow around a wing, we used this result to conclude that the flow was irrotational everywhere because the flow was irrotational far upstream, and all streamlines started out far upstream. But what if there is a region of closed streamlines in the flow? If we consider ideal fluid, with = 0, then there's nothing we can say about the vorticity in such a region beyond that the vorticity is constant along streamlines. However, if we allow the viscosity to be nonzero but still small, i.e. 0, then we can get somewhere. In particular, we have the PrandtlBatchelor theorem, which states that in steady, 2D viscous flow the vorticity is constant throughout any region of closed streamlines in the limit 0. The book goes through the derivation of this. This is a good example of a phenomenon we've discussed occasionally, namely that considering a fluid with no viscosity is not the same as considering a fluid with vanishingly small viscosity. 2 ME 563  Intermediate Fluid Dynamics  Su Lecture 17  The NavierStokes equations: velocity gradient tensor
Reading: Acheson, 6.1, 6.2. Up to now we've discussed the Euler equations that describe the motion of ideal fluids, and the NavierStokes equations for viscous, incompressible fluids, without really explaining where they came from. In fact the NavierStokes equations in particular took a while to be `finalized,' and there are some important assumptions inherent in them. We're going to take some time to develop NavierStokes equations from the ground up. The discussion will not follow the book directly, but the notes and the text are intended to complement each other. Near the beginning of the course we introduced the velocity gradient tensor, u u v w x x x uj u= (1) = u v w . y y y xi u v w
z z z At the time we mentioned this because it represents a vector (the operator) operating on another vector (the velocity, u) and yielding a secondorder tensor. It turns out the velocity gradient tensor is central to the derivation of the NS equations. In this lecture we will discuss the physical implications of u. 1 Deformation, strain, and rotation As we've seen (and as you recall from your introductory fluids class), in a parallel flow of a Newtonian fluid, the velocity gradient normal to the flow direction is directly proportional to the imposed shear stress. That is, if the velocity is u = u(y)^x , the shear stress = (du/dy), where the e proportionality constant we've called the dynamic viscosity. Figure 1 illustrates this for the case of a flow between two parallel plates. Figure 1: Shear stress in a parallel flow. Physically (without being precise), we say that the viscosity of the fluid resists deformation because of internal friction. In this example, the `deformation' takes the form of the shear term du/dy. What about the more general case of nonparallel flow? For simplicity, we'll talk first about twodimensional flow, where the velocity field is u(x, y) = u(x, y)^x + v(x, y)^y , and the velocity gradient tensor is e e u= uj = xi
u x u y v x v y . (2) A priori we might expect that the offdiagonal terms of u, v/x and u/y, describe the shear stress, in analogy with the parallel flow case, but it's not quite so simple. Consider a twodimensional fluid element, with infinitesimal side lengths dx = dy (Fig. 2a), in 1 Figure 2: Shear deformation of a 2D fluid element in parallel flow. a twodimensional, steady (not timevarying) flow field. We are not interested now in expansion or compression terms, only shearing terms, so u(x, y) = u(y) and v(x, y) = v(x). (This is equivalent to looking only at the offdiagonal terms of the velocity gradient tensor.) First, consider the parallel velocity field u = u(y)^x . If the fluid element is square at time t = 0, as in Fig. 2a, then some e (infinitesimal) time later (t = dt) point 1 will have moved a distance u1 dt, and point 2 a distance (u1 + (du/dy) dy)dt, and the fluid element gets deformed as in Fig. 2b. (The side lengths dx = dy are infinitesimal, so we only have to worry about linear velocity variations.) The fluid element is deformed similarly if the velocity field is parallel in the ydirection, u = v(x)^y . In that case, the fluid element will look like Fig. 2c at time t = dt. Since the u and v e velocity components are orthogonal, for the 2D velocity field u = u(y)^x +v(x)^y we can superpose e e the two deformations to find the shape of the fluid element at time t = dt, as in Fig. 3. In the Figure 3: General shear deformation of a 2D fluid element. figure, the displacement vector for point 2 is labeled ds, and can be written ds = du dv dy ex + ^ dx ey dt = ^ dy dx du dv ex + ^ ey dt dx ^ dy dx making use of dx = dy. There are two special cases of this shear deformation. If du/dy = dv/dx, then the fluid element rotates, but remains square (Fig. 4). This corresponds to solid body rotation. In this case, there is actually no `deformation' of the fluid element and no viscous stress. The second special case is where du/dy = dv/dx (Fig. 5). In this case the diagonals of the fluid element maintain their angular orientation. This is called pure shear strain, in which there is no rotation. (Acheson has a useful definition of rotation, namely that there is local rotation in a flow if two initially perpendicular fluid line elements have a nonzero average angular velocity. For the cases in Figs. 4 and 5, this means that we just have to see if the diagonals in the fluid element are rotating.) 2 Figure 4: Pure rotation of a 2D fluid element. Figure 5: Pure strain of a 2D fluid element. Now let's go back and look at only the ondiagonal terms in the 2D velocity gradient tensor (Eq. 2). Let u = u(x)^x + v(y)^y . The effect of the velocity gradient terms du/dx and dv/dy on the e e square fluid element of Fig. 2 is obvious; positive du/dx stretches the element in the xdirection, and positive dv/dy stretches the element in the ydirection. Similarly, negative du/dx and dv/dy represent compressions in the x and ydirections. These are called normal strains, and obviously involve no rotation. Viscous fluids resist these normal strains also, but not in quite as intuitive a way as the shear strains (more about this later). We have now resolved the general `deformation' of a fluid element under the action of velocity gradients, and have shown that viscosity is relevant to those terms that are independent of rotation. For convenience, then, we want to separate the velocity gradient tensor into parts that describe the strain and rotation separately. The way to do this suggests itself if we appeal to the above discussion on shear deformations. We saw that there is no viscous stress for du/dy = dv/dx, and no rotation for du/dy = dv/dx. The logical thing to do, then, is to make the viscous shear stress proportional to (du/dy) + (dv/dx), and the rotation proportional to (du/dy)  (dv/dx). This is the same thing as separating the tensor u into its symmetric and antisymmetric parts. We'll call the symmetric part e and the antisymmetric part 1 u + ( u)T "rateofstrain tensor" 2 1 = (3) u  ( u)T "rotation rate tensor" 2 The factor of 1/2 is incorporated so e + = u. These tensors are called `rates' of strain and rotation because the terms have units (time)1 . Going back to three dimensions, e can be written (using symmetry) e11 e12 e13 1 ui uj e = e12 e22 e23 , where eij = (4) + 2 xj xi e13 e23 e33 e= 3 and can be written (using antisymmetry) 0 12 13 1 0 23 , where ij = = 12 2 13 23 0 ui uj  xi xj (5) The tensor e has six independent components, and three (not surprising, since u has nine). Since only has three independent components, we might expect that the rotation rate can be expressed as a vector. Recalling the definition of a curl from the first lecture, we write = u= = x1 e1 ^ u1 x2 e2 ^ u2 x3 e3 ^ u3 u3 u2 e1 + ^  x2 x3 u3  u2 x3 x2 = u1  u3 x3 x1 u2 u1 x1  x2 u1 u3  x3 x1 e2 + ^ u2 u1  x1 x2 e3 ^ (6) Comparing (6) with (5), we see that 2 ij = k ijk , or 0 3 2 1 = 3 0 1 . 2 2 1 0 (7) The quantity = 1 e1 + 2 e2 + 3 e3 is just the vorticity. We already discussed a little how the ^ ^ ^ vorticity is analogous physically to the angular velocity, or rotation rate. Equation (7) tells us the same thing mathematically. 2 Normal strains and viscosity
u1 x1 u2 = x2 u3 = . x3 Let's look back at the ondiagonal terms of (4): e11 = e22 e33 In Cartesian coordinates, these are just u/x, v/y, and w/z, describing expansion and compression. Notice that, for incompressible flows, eii = e11 + e22 + e33 = u = 0. The sum of the diagonal terms of a tensor is known as its trace. For incompressible flows, then, the trace of the rateofstrain tensor is zero. (This will become interesting later.) To summarize: the physical reason for separating u into the rateofstrain and rotation rate tensors in (3) is because of the effects of viscosity. The components of the rateofstrain tensor e describe motions that are resisted by viscosity; shearing, compression and expansion. The components of the rotation rate tensor describe motions that are not resisted by viscosity. We will incorporate this knowledge into the NS equations in the next lecture. 4 3 Further reading White (on reserve), Section 13.3. This treatment is not particularly logical mathematically. Batchelor (on reserve), Section 2.3. This treatment is very (maybe too) mathematical. Discussions similar to the above can be found in 5 ME 563  Intermediate Fluid Dynamics  Su Lecture 18  The NavierStokes equations: stress tensor
Reading: Acheson, 6.14. In the last class we learned that the deformation of a fluid element is described by the strain rate tensor, e, which is the symmetric part of the velocity gradient tensor, u 1 e= ( u+ 2 We can then write e as uT ). e11 e12 e13 1 e = e21 e22 e23 , where eij = eji = 2 e31 e32 e33 uj ui + xj xi . (1) The significance of this is that, in a viscous fluid, the straining or deforming motions described by e are precisely the motions that are resisted by viscosity. In order to derive the NavierStokes equations, we need to describe the forces that act on a fluid element, and to do that, we will need to incorporate this strain rate tensor. 1 The stress tensor The ideas we apply to derive the NavierStokes equations are pretty general and, in particular, similar ideas pop up in solid mechanics, especially the ideas behind the relationship between stresses, by which we mean forces acting over surfaces, and strains, by which we mean geometric deformations. What we want to do is describe the forces acting on the surfaces of a volume element of fluid. As usual, we will consider a cubic volume element in Cartesian coordinates, and also as usual, the side lengths will be infinitesimal, so any flow quantities (like velocity) will be essentially constant over each face. Obviously, there will be a pressure force on each face i equal to pi dSi , where pi Figure 1: Our usual cubic volume element. is the pressure on the face i and dSi is the area of the face. Since we are interested in the forces acting on the volume element by the fluid outside of the element, the pressure force points normal to each face and towards the inside of the element. There will also be viscous forces acting on each face i. Suppose we are talking about face 1 in Fig. 1. On face 1, there are three possible viscous stress terms, as depicted in Fig. 2. There is a normal viscous stress if du/dx = 0. If du/dx is positive, the force acting on the fluid element across face 1 is in the positive xdirection, as shown. There are also two shear viscous stresses, 1 Figure 2: Viscous forces acting on face 1 of the volume element in Fig. 1. shown in the figure as simple parallel flow shearing motions. If dv/dx is positive, there is a shear force acting on the fluid element across face 1 in the positive ydirection. If dw/dx is positive, the resulting shear force across face 1 is in the positive zdirection. (If these velocity gradient terms are negative, these forces are clearly in the appropriate negative directions.) Given this multitude of force terms on each face, the handy way to keep track of them is as follows. We define the stress tensor, T, such that T11 T12 T13 T = T21 T22 T23 , where T31 T32 T33 Tij is the icomponent of stress acting on a surface whose normal vector is in the jdirection. Thus, if the area of the surface is dS, the icomponent of force acting on the surface is Tij dS. We can make this definition a bit more general. For a planar surface with arbitrary orientation, ^ i.e. a planar surface whose normal vector n is not aligned with one of the coordinate axes, so that ^ n = nj ej , we can say: ^ Tij is the icomponent of stress acting on the surface that is associated with the jcomponent ^ of the normal vector n of that surface. Thus, the net icomponent of force acting on the surface is Tij nj dS (using the implied summation convention). As is evident from the above definition(s), T has two directions associated with it. This points up one of the really confusing aspects of (secondorder) tensors as opposed to vectors vectors have a single direction associated with them, but secondorder tensors have two. In (2), we are basically writing the ith row of the stress tensor T, which is the part of T that describes forces acting in the idirection, as a vector Ti , which is thus written Ti = Tij ej ^ so that (2) becomes ^ Tij nj dS = Ti n dS. Before we concern ourselves with the exact mathematical form of T, let's look at how T fits into the equation of motion for a fluid element. 2 (2) 2 Cauchy's equation of motion Consider an arbitrary fluid element, enclosed by a surface S. The net force in the idirection acting on the (entire) surface is ^ Ti n dS = Tij nj dS = Ti dV = Tij dV xj S S V V where we have applied the divergence theorem. As we've done before, we can let the arbitrary fluid element be of infinitesimal size, which lets us write Tij Tij dV = dV xj xj (3) V since the integrand is basically constant over the infinitesimal element. Meanwhile, the gravitational force on the fluid element is gi dV (4) where gi is the component of the gravity vector in the idirection. (We are not constraining the axes so that one of them is necessarily vertical.) The forces described by (3) and (4) are just equal to the mass of the fluid element times the icomponent of its acceleration, where the mass is dV and the acceleration is given by the material derivative of ui . Thus Dui Tij + gi , = Dt xj (5) where we've canceled the common dV terms. There is a little bit more precise derivation of (5) in the book (p.206). 3 Incompressible, Newtonian viscous fluids We said above that the surface of an element of viscous fluid experiences pressure forces, shear viscous forces, and normal viscous forces. Incorporating the pressure in the stress tensor T is easy. Obviously the pressure is a normal force, so we have Tij = p ij + (viscous terms). (6) The minus sign on the pressure ensures that the pressure force imposed on a fluid element always points into the fluid element. The Kronecker delta tells ensures that the pressure is a normal force only, i.e. it appears only when i = j. For the viscous forces, we already showed in the last class that the strain rate tensor (1) describes the motions resisted by viscosity. So we anticipate that the strain rate figures in T. On top of that, for Newtonian fluids, we said that the stresses were linearly proportional to strains, with a proportionality constant that we called the dynamic viscosity. This allows us to write (6) as Tij = p ij + 2eij = p ij + uj ui + xi xj . (7) There are some interesting properties of the form of this equation, and in fact you can work backwards and arrive at (7) by assuming these properties and then applying some tensor math (for the interested, the Batchelor book on reserve goes over this in detail). The properties, first enumerated by Stokes, are: 3 The stress tensor is a linear function of the velocity gradient components (this is essentially equivalent to the Newtonian assumption). The tensor is isotropic, meaning it doesn't care what coordinate system you choose to express it. If the fluid isn't deforming, that is if the strain rates are zero, T reduces to the pressure terms only. In deriving (7) we have also implicitly assumed that the flows described are incompressible. The isotropy property is interesting because, for starters, it tells us that the fluid resists stretching and compression in the same quantitative way that it resists shearing. In the equation this is expressed by there being only a single viscosity, . Thus the contribution of a normal strain, say u1 /x1 , to the stress is (u1 /x1 ), while the contribution of a shear strain, say u2 /x1 , to the stress is (u2 /x1 ). It seems kind of strange that the same coefficient should be valid for both normal and shear stresses. But then, it seems less strange if we look at a simple example. Consider a twodimensional velocity field defined as (it's as usual easier to conceive of these motions in two dimensions): u 1 = x2 , u2 = x1 which corresponds to a pure deformation like that on the left in Fig. 3. Notice that because of the Figure 3: Example of rotation of coordinate systems. symmetry of the strain rate tensor, the stress components T12 and T21 are equal (this is another key property of the stress tensor: think about what it implies for a simple parallel shear flow like u = (u(y), 0, 0)). Essentially, what we get is a stretching of the square fluid along the one diagonal, and a compression on the other. Now, what happens if our coordinate system happens to be aligned with those diagonals? Then, we get the situation depicted on the right in Fig. 3. The original shear terms T12 and T21 correspond to normal terms T1 1 and T2 2 in the new primed axes e1 and e2 . So ^ ^ it's clear why only one viscosity coefficient is used to describe shear and normal stresses; because by changes in coordinate systems, you can literally change shear terms into normal terms. 4 ME 563  Intermediate Fluid Dynamics  Su
Lecture 19  The NavierStokes equations: more on the stress tensor, and kinetic energy transport
Reading: Acheson, 6.4. At the end of the last lecture we saw saw why a single proportionality constant, the dynamic viscosity, , relates shear strain to shear stresses and also normal strain to normal stresses; the reason is that, just by rotating our coordinate systems, we can actually change shearing motions to compressional or extensional motions. We showed this schematically in terms of a stretching of a twodimensional, square fluid element along one of its diagonals, which turns into a combination of extension in one axis and compression in the other when we rotate our coordinate system by 45 . 1 Rotation of coordinate systems What exactly do we mean by rotating our coordinate systems? Consider Fig. 1, depicting a coordinate system with axes e1 and e2 and a second system, obtained by a counterclockwise ^ ^ rotation through , with axes e1 and e2 . As can be seen (with some thought) from the figure, for ^ ^ Figure 1: Relation between rotated coordinate systems. a given point in the plane, the coordinates (x1 , x2 ) relate to the coordinates in the rotated system, (x1 , x2 ), by x1 = x1 cos + x2 sin x2 = x1 sin + x2 cos which can be expressed in terms of vectors and tensors as x = Rx, where x = x1 x2 , R= cos sin ,  sin cos x= x1 x2 . We can also take the inverse of this rotation, defined as x = R1 x , 1 (1) where R1 R = RR1 = I, (2) with I being the identity matrix, with ones on the diagonal and zeroes everywhere else. How can we express tensors in these rotated coordinate systems? An easy way to see this is as follows. Often, a tensor will be multiplied by a vector to yield another vector, which can be written: Tx = y (3) where T is a tensor and x and y are vectors. The same relationship can be written in the rotated coordinate system x = Rx: Tx =y (4) where T is the original tensor T expressed in the rotated coordinate system. To relate T to T , we can apply (1) to (3), giving TR1 x = R1 y . Multiplying both sides by R, and using (2), we get RTR1 x = y . And, comparing (5) with (4), we see that RTR1 = T . (5) 2 Topology of the strain rate field There is a useful result from linear algebra and tensor math that says that, for any symmetric tensor T, there is a coordinate rotation represented by R such that the tensor expressed in the rotated coordinate system, T , has nonzero terms only on the diagonal. This has direct implications for the strain rate tensor, e. Since e is symmetric, it can be rotated into a coordinate system where it is diagonal, i.e. 0 0 e = 0 0 . (6) 0 0 Since the diagonal terms of the strain rate tensor describe normal strains (extensions and compressions) and the offdiagonal terms describe shears, (6) describes a strain field that has only normal strains and doesn't have shear strains. This formalizes what we discussed in the example at the end of the last lecture, that shear strains can become normal strains by rotating the coordinate system. The diagonal terms in the diagonalized strain rate tensor are called the principal strain rates. We can say even a little more about the strain field. We have also pointed out that the incompressibility condition, u = 0, means that the trace of the strain rate tensor (meaning the sum of the terms on the diagonal) has to be zero. This has to be true whatever Cartisian coordinate system we're looking at, independent of rotation. Let these diagonal terms be , and as in (6), and let us say that has the highest value and the lowest. Since + + = 0, this means that > 0, < 0, (if both and are zero, things are pretty uninteresting, since then the whole strain rate tensor is zero). The question then becomes, what is the sign of the intermediate principal strain rate, ? The three cases are: 2 Figure 2: Canonical strain rate topologies: left, positive intermediate principal strain rate, right, negative intermediate principal strain rate. 1. = 0: twodimensional flow, with one extensional principal strain rate and one compressional principal strain rate. 2. > 0: this is depicted on the left in Fig. 2. There are two extensional strains and one compressional. This can be viewed as pressing regions of fluid into `sheets.' 3. < 0: depicted on the right in Fig. 2. There is one extensional strain and two compressional. This can be viewed as stretching regions of fluid into `tubes.' The irrotational parts (that is, the parts described by the strain tensor) of any incompressible velocity field have to fit into one of these three categories. These socalled canonical topologies are widely studied in the context of turbulent flows. 3 Vector form of the NavierStokes equation
Tij Dui + gi , = Dt xj uj ui + xi xj We have gotten up to the following point with the NavierStokes equation: where Tij = p ij + . (8) (7) Let's complete the analysis and show how this gives us the form we've been using in class up to now. Combining (7) and (8), we get uj Dui ui p ij + + gi + = Dt xj xi xj uj p ui + gi = + + xi xj xi xj uj 2 ui p + = + + gi xi xj xi xj xj (9) Notice how we got rid of the Kronecker delta, ij , in the pressure term. Now, looking at the last line, notice that we can interchange the order of partial differentiation in the second term. That is, xj uj xi = 3 xi uj xj . But uj = xj u=0 since the flow is incompressible. Thus (9) becomes which, in vector notation, is just with which we are familiar. Du = P + Dt
2 Dui 2 ui p + + gi , = Dt xi xj xj u + g 4 Transport of kinetic energy Acheson talks about the transport of kinetic energy in an interesting fashion (one that I haven't seen in any other texts). The kinetic energy of the fluid in a volume element V enclosed by surface S is: 1 T = ui ui dV 2 V (Using T to denote the kinetic energy may seem confusing since we're already using T for the stress tensor, but the notation is standard in physics. Also, I'm not keen on how Acheson writes u2 instead of ui ui can you think why?) The time rate of change of this T is then i dT = dt
V ui Tij dV + xj V ui gi dV. In the course of the derivation Acheson has these intermediate steps: This step uses the product rule: Tij ui dV = xj V (ui Tij ) dV  xj ui dV xj V V Tij This step uses the divergence theorem: (ui Tij ) dV = V xj S ui Tij nj dS = S ui ti dS where the divergence theorem is applied using the j index as the direction index in the term ui Tij . That is, we treat ui Tij as the jcomponent of a vector u T. The last equality uses the definition we introduced earlier of the force acting on a surface. This step uses the fact of incompressibility: 1 Tij 2 ui uj + xj xi 1 ui ui uj uj =  p ij + + + 2 xj xi 2 xj xi ui = p + 2eij eij = 2eij eij xi
work done on fluid element 2 The result is:
gravitational potential energy viscous dissipation dT = dt V u g dV +
S t u dS  2 V eij eij dV The third term represents viscous losses, which are analogous to friction losses. 4 ME 563  Intermediate Fluid Dynamics  Su
Lecture 20  The NavierStokes equations: kinetic energy dissipation
Reading: Acheson, 6.5. At the end of the last lecture we looked at the kinetic energy of the fluid in a volume element, V , enclosed by a surface S, T = 1 2
V ui ui dV, and we saw that the time evolution of this quantity is given by
gravitational potential energy work done on fluid element viscous dissipation dT = dt V u g dV +
S t u dS  2 V eij eij dV . The first term on the righthand side represents changes in the gravitational potential energy; the second term represents work done on the surface of the fluid element, by pressure and viscous forces; and the third term represents viscous energy losses. This last term is, obviously, strictly negative, and expresses that whenever there is any deformation of fluid elements in a viscous fluid, energy will be dissipated due to viscosity. This viscous kinetic energy dissipation term is neither a conversion of kinetic energy to potential energy nor a work done by the fluid element on the surroundings, but is rather a conversion of kinetic energy to heat. Thus these viscous losses are analogous to friction losses. There is a very interesting result for the energy dissipation that was first shown by Helmholtz. Before we write that we will take a brief detour into physics. 1 Minimization principles One of the really attractive ways of thinking about physical problems is in terms of minimization principles. An example of these is the refraction of light rays at the interfaces of two media with different indices of refraction. Figure 1 shows a light ray traveling across a plane interface, from a n2 n1 1 2 Figure 1: Schematic for Snell's law of refraction. medium with index of refraction n1 into one with index n2 . Define 1 as the angle the ray makes in the first medium with the normal to the interface, and 2 as the same angle in the second medium. According to Snell's law of refraction, these angles will be related by: sin 1 n2 = . sin 2 n1 1 (1) This is useful enough, but a much more insightful description of refraction (which also works for reflection) is the principle of least time for light rays (due to Fermat, of Last Theorem fame), which states that A light ray always travels from one point to another in a medium by a path that requires the least time. We can get from this principle to Snell's law by remembering that the speed of light in a medium with refractive index n equals c/n, where c is the speed of light in a vacuum. Thus Snell's law and the principle of least time are equivalent, but clearly they have much different characters. In some sense, Snell's law is descriptive of the phenomenon of refraction, but the principle of least time is explanatory. Think of this another way: Snell's law is of no use to you if you don't remember the exact formula (1), but you can always derive the formula if you remember the principle of least time. One consequence of the distinction is that Snell's law is useless if the refractive index of the medium of interest varies continuously, while the principle of least time can still be used to predict light propagation paths. There are other wellknown minimization principles in physics. Physical bodies seek to minimize their gravitational potential energy, which tells you how to find the shape of a hanging chain. Then there is Hamilton's principle, which states that the path followed by a dynamical system minimizes the time integral of the difference between the kinetic and potential energies. You can actually use Hamilton's principle to calculate the equations of motion of a body without using any Newtonian theory, which is pretty remarkable. The question here is: is there any equivalent explanatory principle for fluids? 2 Minimum dissipation theorem It turns out that there is, but it's not as general as in the cases of light propagation or dynamical systems. As shown by Helmholtz, for incompressible viscous flow with negligible inertia forces, The actual flow has a smaller total rate of kinetic energy dissipation than any other incompressible flow in the same region having the same values of the velocity vector everywhere on the boundary of the region. That is, for given boundary conditions, viscous flows with negligible inertial terms arrange themselves so as to lose the least amount of energy to viscous dissipation. The constraint of negligible inertial terms here means that u and (u t )u the viscous and pressure terms which holds for steady parallel flows, for example (why?). Let's look at an example. Consider plane Poiseuille flow of a viscous fluid (Fig. 2). However, we're going to approach this y x y=1 ?
Figure 2: Plane Poiseuille flow. y = 1 as if we've never heard of the NavierStokes equations. Instead we will ask the question 2 Given a fixed mass flow rate per unit span, m, between the planes, what is the velocity profile, u(y)? We will make one assumption just for convenience: we will consider only velocity profiles of the form u(y) = A + By , (2) where A and B are constants, and can be any number. What we really want to know is the value of , which tells us the shape of the velocity profile. Because the fluid is viscous, the noslip condition pertains, which automatically tells us A = B, so we can write u(y) = A  Ay . Now let's apply the fixed mass flow assumption, which will give us our value A.
1 m = 2 0 (A  Ay ) dy, where is the fluid density. Integrating, we get m = 2 Ay  so that A= and u(y) = m( + 1) (1  y ). 2 (3) m( + 1) , 2 1 Ay +1 +1
1 0 = 2A 1  1 +1 Now we want to minimize the kinetic energy dissipation. The dissipation at any point is equal to 2eij eij , but here the only velocity gradient term is du/dy, so the dissipation is 2(du/dy)2 . We want to minimize this quantity integrated across the region between the planar boundaries. First write the total dissipation:
1 2
0 m( + 1) 1  y 2 2 m2 ( + 1)2 dy = 22 1 y 22 dy.
0 Integrating, we get m2 ( + 1)2 1 y 21 2 2 2  1
1 =
0 m2 ( + 1)2 . 22 (2  1) We want the value of that minimizes this. In particular we'll look for the value of that minimizes ( + 1)2 . 2  1 3 (4) So we want to set the derivative of (4) equal to zero: d ( + 1)2 22  2  4 2( + 1) 2( + 1)2 = . =  d 2  1 2  1 (2  1)2 (2  1)2 The numerator has roots = 1 and = 2. The = 1 solution is unphysical, because the velocity u would go to infinity at y = 0. So the proper value of is 2 (it is trivial to show that this represents a minimum of (4). This gives us a parabolic velocity profile, which is exactly what we found using the NavierStokes equations. What is the utility of this minimum dissipation formulation? Surprisingly (to me), this concept is never discussed in classes and appears in only one text I'm familiar with (the Batchelor text on reserve). I think it's very useful in that it gives you a clear sense of what the property of viscosity does in a flow. Minimizing dissipation is equivalent to minimizing deformation, which is similar to our prior description that viscosity resists deformations. From a practical standpoint, the minimum dissipation idea is useful in variety of problems including lubrication theory, flow in porous media, etc. (There is one major hurdle, which is that minimization ideas in general need lots of variational calculus....) 4 ME 563  Intermediate Fluid Dynamics  Su
Lecture 21  The energy equation
Reading: Acheson, 8.18.2. In this lecture we will finish up with the equations describing fluid flow. Next, we'll get into boundary layers (note the reading). 1 Energy equation Acheson doesn't cover it, but the energy equation (more specifically, the internal energy equation) is significant to a lot of engineering problems. To derive it, start by considering the fluid contained in a volume V , enclosed by S. This discussion will repeat some of the ideas we used in deriving the kinetic energy equation, but the distinctions will be clear. Work is done on the fluid element by body forces (gravity) and surface forces. The rate at which work is done on the fluid by gravity is ui gi dV. (1) V To interpret this, note that gi dV is the gravitational force acting on the fluid in an infinitesimal element dV . The work done by that force on a body (in this case, the fluid in dV ) is equal to the force times the distance traveled in the direction of the force. Thus the rate at which work is done on the fluid in dV , i.e. the work per unit time, is the force times the distance traveled per unit time = the force times the velocity in the direction of the force. The rate at which work is done on the fluid in V by surface forces is computed similarly. On ^ a given infinitesimal surface element dS, with normal vector n = nj ej , the surface forces in the ^ idirection are Tij nj dS (notice the implied summation over j). The rate at which work is done on the fluid in V by the surface forces on dS is thus ui Tij nj dS. (2) Again notice the implied summation over i, which means that (2) represents the sum of the work done by all of the components of the surface forces. Integrating over the entire surface, the rate at which work is done on the fluid in V by surface forces is ui Tij nj dS = ui Tij dV, xj (3) S V where we have applied the divergence theorem. To understand this, write the integral on the left side as
S ^ (u T) n dS and note that u T is a vector whose components are given by the j index of Tij . Combining (1) and (3), the rate at which work is done on the fluid in V by all forces is ui gi + ui Tij ui dV, + Tij xj xj (4) V where we have expanded (3) using the ubiquitous product rule. 1 Let's consider just the term in the square brackets in (4), which represents the rate at which work is done on the fluid in the infinitesimal element dV . Recall from a couple of lectures ago that the NavierStokes equation can be written Tij Dui + gi , = Dt xj so we can write the term in the square brackets in (4) as Rate of working on the fluid in dV = ui But the first term on the right can be written ui Dui D 1 = u i ui , Dt Dt 2 Dui ui . + Tij Dt xj (5) which is just the rate of change of kinetic energy (per unit volume). Thus the term Tij ui xj (6) in (5) is the rate of change of internal energy of the fluid in dV due to work. There will also be heat transferred into the fluid by molecular conduction, or diffusion. The rate of heat conduction per unit area is given by q = k , where k is the thermal conductivity (which we'll treat as constant) and is the temperature. Since q is a gradient, the rate of heat conduction per unit area in the direction of a vector s is q s. Using this and integrating, the rate of gain of heat by the fluid in V is k ni dS = xi k 2 dV. xi xi (7) S V Two observations. First, the integrand is positive although the q is defined with a minus sign, ^ because the normal vector n is defined positive out of V , and we want to know the heat flux into V . Second, we've applied the divergence theorem as usual. From (7), we can write that the rate of transfer of heat to the fluid in dV due to conduction is k 2 . xi xi (8) A brief appeal to thermodynamics now. We can write the change in internal energy (per unit mass), E, of a fluid as E = Q + W, where Q is heat transferred to the fluid and W is the work performed. Here, our Q is given by (8) and our W is given by (6). Thus, we can write the rate of change of internal energy for an infinitesimal fluid element as DE ui 2 +k . = Tij Dt xj xi xi (9) We can simplify this a bit. The first term on the right can be written 1 ui 1 uj + Tji Tij 2 xj 2 xi 2 since i and j are just dummy variables (we can treat them as such because we sum over them). But T is symmetric, so this becomes Tij We can write this out too: Tij eij = (pij + 2eij ) eij = p eii + 2eij eij = 2eij eij , since eii = u = 0 for incompressible flows. But this is just the kinetic energy dissipation term! Thus, the energy equation is DE 2 , = 2eij eij + k Dt xi xi (10) 1 ui 1 uj + 2 xj 2 xi = Tij eij . which says that the internal energy of the fluid increases because of the heat generated by viscous dissipation, and because of thermal conduction. Observe that we have not considered radiative heat transfer, or heat release by chemical reactions etc., which are among many effects that can modify this equation. 3 ME 563  Intermediate Fluid Dynamics  Su
Lecture 22  Boundary layers: physical insight
Reading: Acheson, 8.18.3. The study of boundary layers occupies a central place in fluid dynamics for a number of reasons. For one, boundary layer theory sits at the intersection of inviscid and viscous flow theory, which is one of the major conceptual divisions in fluids. For another, boundary layers are seen in some form in almost every flow of practical engineering interest. In fact, this is true even in the absence of solid surfaces; the relevant `boundary' can be between, say, high and lowspeed fluid streams. In general, we will define a boundary layer as follows: Boundary layer: a thin layer in which the effect of viscosity is important. Under what conditions will viscosity be important? We know from our discussion of the NavierStokes equations that viscosity comes into play when we deform fluid elements, so we can expect that viscosity will be important where velocity gradients are high. When will that happen? As an example, any time we introduce a stationary solid surface into a flow of viscous fluid, a boundary layer will form, because of the noslip condition. Near the solid surface, the flow needs to adjust itself from the main flow velocity to the velocity of the solid surface. This gives rise to a velocity gradient, which makes viscous effects important. Just because a boundary layer exists in a flow doesn't mean that the boundary layer is significant to the flow properties. Think back to the airfoil example we talked about near the beginning of class. We mentioned that at low angles of attack, the lift of a wing is well predicted by ideal flow theory, which ignores viscosity and thus ignores any possible boundary layers. However, at higher angles of attack, ideal flow theory is less able to predict lift. What causes the difference? The answer was boundary layer separation, also known as flow separation. When the boundary layer is attached to the surface of the wing, it remains thin, and the percentage of the flow volume in which viscous effects are significant is small. When the boundary layer separates, though, the percentage of the flow volume in which viscous effects are significant rises dramatically, and ideal flow theory is inadequate to describe the flow. It's worthwhile look at boundary layer separation in more detail. 1 Boundary layer separation A boundary layer can remain attached to a surface if the flow streamlines in the boundary layer have no difficulty following the surface contours. To understand what might cause the boundary layer to be unable to follow the surface contours, consider Fig. 1, depicting the boundary layer on a flat plate, where the mean flow is from left to right. In the leftmost cartoon, the flow is moving in the positive xdirection, pushed along by the negative pressure gradient, dp/dx < 0. The velocity gradient at the surface, du/dy (y = 0), is positive. Now suppose we bump up the pressure gradient at the surface. As the pressure gradient changes sign and becomes positive, it will want to push the flow in the negative xdirection. If the flow has some inertia, it will tend to push against this adverse pressure gradient. At some point, these effects will be in balance, and the situation in the middle cartoon in Fig. 1 will pertain, where du/dy at the surface is zero. As the pressure gradient is increased past this point, we get flow reversal near the surface (the cartoon on the right). When this happens, the lefttoright boundary layer flow can no longer follow the contour of the surface. Instead, the boundary layer is diverted upward to go around the region of reversed flow. This process is called boundary layer separation. The favorite example of fluids instructors and textbooks everywhere for illustrating boundary layer separation is flow around a sphere. Figure 2 shows schematically what happens if the fluid is 1 Figure 1: Representation of boundary layer separation. inviscid. As the flow moves around the front half of the sphere, it accelerates (basically, the fluid has to fit through a smaller space), and by the Bernoulli theorems, the pressure decreases. Then, as the flow moves past the back half of the sphere, the velocity decreases and the pressure increases until both return to their initial values. The velocity and pressure profiles are symmetric between Figure 2: Ideal fluid flow around a sphere. the front and back of the sphere, because (in the absence of viscosity) there is no frictiontype effect to sap energy from the flow. In particular, there is no pressure drag on the sphere (part of what D'Alembert noted in his paradox). Now let's consider what happens of the viscosity of the fluid is nonzero. This situation is depicted in Fig. 3. There is now a boundary layer at the sphere's surface. On the front half of the sphere, the pressure at the surface is decreasing as we move in the flow direction, so according to the examples of Fig. 1, the boundary layer flows happily along the surface of the sphere. However, as we move to the back half of the sphere, the pressure at the surface increases as we move in the flow direction. Eventually, the boundary layer will be unable to fight against this adverse pressure gradient, and the flow will separate. After that, the pressure and velocity profiles essentially freeze; that is, the profiles don't go on to recover to their initial values. The resulting pressure distribution looks something like that shown in Fig. 3. The net pressure on the front of the sphere is higher than that on the back of the sphere, meaning that there is a pressure drag on the sphere (as everyone is amply familiar with in the real world). A good illustration of this boundary layer separation phenomenon is given in Fig. 4. As is wellknown, the pressure drag on a sphere is lowered when the boundary layer flow is turbulent. This is because a turbulent boundary layer is better able to overcome the adverse pressure gradient. We will discuss the reasons for this when we discuss turbulent flows later in the semester. The boundary layer can be caused to become turbulent if the sphere's surface is 2 Figure 3: Viscous fluid flow around a sphere. Figure 4: Boundary layer separation on a sphere (from van Dyke, An Album of Fluid Motion.) somehow rough examples are the dimples on a golf ball or the raised seams on a baseball. In Fig. 5, the boundary layer is made turbulent by a wire placed around the sphere upstream of its equator. By comparing Fig. 5 with Fig. 4, it is clear that the turbulent boundary layer gets further around the sphere before separating. The result is that the pressure recovers better around the back of the sphere, so the pressure on the back of the sphere is higher for the turbulent boundary layer as compared with the laminar boundary layer. Thus, by `tripping' the boundary layer into turbulence, we can reduce the pressure difference between the front and back of the sphere, and we reduce the pressure drag on the sphere. We might reasonably expect that the same physical reasoning pertains when we talk about lift of bodies (wings, for example) in viscous flows instead of drag. After all, drag arises from the difference in pressure between the front and back of a body, and lift arises from the pressure difference between the top and bottom (Fig. 6). Let's see if our reasoning regarding drag on a sphere carries over to the lift generated by a wing. Consider the left cartoon in Fig. 7. In this case, the wing makes an angle of attack with the incoming flow, , of zero. The boundary layer stays attached over the whole wing, by design. The long, tapered rear to a wing lessens the adverse pressure gradient experienced by the boundary layer, by lengthening the distance over which the pressure recovers. The pressure profiles on the top surface and bottom surface of the wing are also shown in the figure. (Wings are designed so the pressure drop on the top of the wing is much higher than the pressure increase on the bottom surface.) 3 Figure 5: A turbulent boundary layer on a sphere (from van Dyke). Figure 6: Drag vs. lift. Now suppose we increase a little. As we learned in the first chapter of Acheson, as increases, the lift of the wing initially increases. How does this happen? The middle cartoon in Fig. 7 gives an idea. At increased , the flow on the wing's upper surface can't stay attached across the whole surface. (What might cause this?) A good photo of this is shown in Fig. 8. If we follow the reasoning we gave for sphere drag, what happens is that after the separation point on the upper surface, the pressure stops rising to reattain its original value, and so the pressure on top of the wing is lower than it would be if there were no separation. Lower pressure on top of the wing means higher lift, as we can observe experimentally. But...what happens if we keep bumping up the angle of attack, (the rightmost cartoon in Fig. 7)? This keeps moving the separation point forward, and what should happen if our reasoning is still valid is that the pressure on top should keep going down, and the lift should keep going up. But we know that it doesn't. What's going on? I have seen one text which explains that the separation region on the upper surface of the wing is a region of `dead air,' where the pressure is (allegedly) roughly equal to ther pressure in the in the undisturbed air stream infinitely far from the wing. According to this explanation, the drop in lift at high angles of attack is because this region of `dead air' dominates the upper surface of the wing, so the underpressure effect of the top of the wing is drastically rreduced. This sounds reasonable, but what about the case where the wing makes an angle of attack of = 90 ? This `dead air' explanation then predicts that there is no drag! The point is, neither intuitive explanation for the effect of the separation point above the wing is satisfactory. The idea that separation `freezes' the pressure profile doesn't work when the angle of attack gets high, and the idea that a region of `dead air' forms doesn't explain the behavior of a wing at low angles of attack and also contradicts our knowledge of drag behind bodies. Math is useful only to an extent here, too; what we really need are experiments. 4 Figure 7: Lift on an airfoil at different angles of attack Figure 8: Separation on the upper surface of a wing (from van Dyke). 5 ME 563  Intermediate Fluid Dynamics  Su
Lecture 23  Boundary layers: math basics
Reading: Acheson, 8.28.5. We now want to look at boundary layers in a more rigorous way. It will be helpful to keep in mind exactly what we mean when we speak of boundary layers. In the last lecture we gave the following definition: Boundary layer : a thin layer in a flow where viscous effects are important. Implicit in this definition is that viscous effects are not important outide of this boundary layer. This is, of course, exactly the situation that we considered with the regard to the problem of the flow of air over a wing. The viscosity of air is low enough that it can effectively be treated as inviscid, except near the surface of the wing, where the noslip condition ensures that viscous forces will be significant. To help us get a sense of the mathematical approach we'll use, let's look at a sample problem. 1 A sample differential equation
d2 u du + =1 dy 2 dy Consider the following problem (starting on p. 269 of the text) concerning a function u(y): with boundary conditions u(0) = 0, u(1) = 2. (2) (1) The variable is a small positive constant. The general solution to (1) with the given boundary conditions is u(y) = y + 1  ey/ . 1  e1/ (3) (We'll take this as given without bothering to derive it it's trivial to verify that (3) satisfies (1) and the boundary conditions.) Exponentials are very sensitive to the values of their arguments, and in particular the exponential term ey/ is very sensitive to the relative values of y and . The solution can be divided into two asymptotic cases. If y , then ey/ 1, and we get from (3) the main flow (Acheson: `mainstream') solution uM (y) = y + 1, while if y solution (4) , then ey/ is of order 1, and y 0, and we have from (3) the boundary layer uBL = 1  ey/ . (5) Graphically, the solution to (1), and the asymptotic limits (4) and (5), are as shown in Fig. 1 for = 0.05. The exact solution (3) matches (5) when y , which defines the `boundary layer' in the problem, and (3) matches (4) as y gets large. The boundary layer solution (5) depends on , while the main flow solution (4) does not. In a sense, then, the parameter is analogous to the viscosity, 1 1 0.8 BL solution main flow soln exact soln 0.6 0.4 0.2 0 0 0.5 1 1.5 2 Figure 1: Solution to (1), for = 0.05. because it is significant only in the boundary layer. In fact, the problem defined by (1) is quite a good representation for the NavierStokes equations, because this `viscosity' is the coefficient of the second derivative in (1), just as the kinematic viscosity is the coefficient of the Laplacian in the NS equations. In the above discussion, we proceeded by looking at the general solution (3) and identifying different solution regimes based on the relative values of y and . We can take an alternate approach, where we determine the asymptotic solutions by looking at the original equation (1). Let's look for an `outer' solution, uo (y) by ignoring the `viscous' term (d2 u/dy 2 ). This is the same thing as letting 0 for a given value of y. Our differential equation is then duo = 1, dy (6) which is a firstorder equation, so we get to use only one boundary condition from (2). The logical one is the boundary condition that u = 2 at y = 1, since it's the one in the main or outer part of the flow. The result is that uo = y + 1. (7) Now we want to find an `inner' solution, ui , where the viscosity is important. The way we proceed is by rescaling our y variable so that the viscous term (d2 u/dy 2 ) is of the same order as du/dy, the other term involving u in the equation (1), recognizing that is small. In other words, we want to rewrite the equation for a new variable Y = y/L, where we need to select the proper value of L. Plugging (8) into (1), we get 1 d2 ui 1 dui + = 1, 2 dY 2 L L dY (8) and obviously the two terms on the left side are of the same order if L = , and Y is of order one, which is the same thing as letting 0 while Y = y/ stays of order one. This is the point of 2 using the rescaled independent variable Y ; we want to be able to consider values of the independent variable that aren't at the limit of extremely small. One benefit to having the two terms on the left side be of the same order is computational if one term is much bigger than the other, but you don't want to neglect the small one, it's hard to do arithmetic on those terms because your computer only uses a finite number of significant digits. Anyway, with Y = y/, our equation becomes d2 ui dui + = 0, (9) dY 2 dY which is our boundary layer equation for this problem. The relevant boundary condition from (2) is u(0) = 0. By integrating (9), we first get dui = ui + C, dY and integrating again (separating variables) we get ui = C + DeY . Applying ui (0) = 0, we get ui = C(1  eY ) (10) as our `inner' solution for the problem. Notice that one free constant, C, still remains. As our last step, we have to ensure that the outer solution (7) matches with the inner solution (10) at the appropriate place. The appropriate place is for y 0 in the outer solution (i.e. as the outer solution approaches the boundary) and for Y in the inner solution (i.e. as the inner solution gets far from the boundary). These limits are uo 1, and ui C, so C = 1. So finally we have u(y) = y + 1 for fixed y as 0 1  ey/ for y/ of order one as 0, (11) as our solution, which agrees with the equations (4) and (5) we found earlier. This second method of arriving at the solutions for u(y) in the boundary layer, and in the main flow, is interesting because we didn't need to know the exact form of the solution to equation (1). Instead, we applied boundary layertype assumptions to simplify (1), and solved the simplified equations. Using this approach to look at fluid flow problems is very useful, because in general we can't find exact solutions to the NavierStokes equations. What we can do, instead, is assume the flow is inviscid away from boundaries, which gives us the Euler equations, which are easier to solve; near the boundary, we will apply boundary layer assumptions in the hope that we will get simplified equations that we can solve. 2 The 2D steady boundary layer equations
u u 1 p +v = + x y x v v 1 p u +v = + x y y u v + = 0. x y 3 2u 2u + 2 , x2 y 2v 2v + 2 , x2 y The exact NavierStokes equations for twodimensional flow are u (12) (13) (14) We can simplify these considerably for a boundary layer near a wall (say) at y = 0. (This is well covered in the text the discussion here will be sketchy.) Our basic simplifying assumption is that u and v vary much more rapidly in the ydirection than in the xdirection. One consequence of this can be seen by looking at (14), which is the 2D continuity equation u = 0. We need the terms u/x and v/y to be of the same order. Otherwise, u would not vary with x and v would not vary with y (meaning v = 0 everywhere), and the problem would be trivialized our purpose for deriving equations for the boundary layer is to look at their evolution in the xdirection. Since the xderivative is much smaller than the yderivative for both u and v, (14) tells us that v is much smaller than u. Now, look at (12) and (13) and compare them termbyterm. We get u since u v, and v again since u is that u y v v , y u x u v , x v. The same applies to the second derivative terms on the right sides. The result p x p , y (notice: not the same behavior as with the velocity components!) which means that p p(x), and we can concern ourselves with the uequation, (12), and basically ignore the vequation, (13). Finally, in the uequation, looking at the right side, we have 2u y 2 2u , x2 because the yderivative is much larger than the xderivative, and this carries over to the second derivatives. So finally we have u u u 1 dp 2u +v = + 2 x y dx y u v + =0 x y (15) as our 2D boundary layer equations. 4 ME 563  Intermediate Fluid Dynamics  Su
Lecture 24  Boundary layers: the boundary layer on a flat plate
Reading: Acheson, 8.28.5. In Lecture 23 we developed the boundary layer equations for 2D, steady flow, which were u u u 1 dp 2u +v = + 2 x y dx y u v + = 0. x y (1) (2) To derive these, our key assumption was that u and v vary much more rapidly in the ydirection than in the xdirection. In terms of derivatives, this meant u y v y u and x v , x (3) which allowed us to toss out various terms in the NavierStokes equations. We can get some interesting insight by looking further at (3). Let U0 be a characteristic value of u, and suppose that u varies by order U0 over a distance L in the xdirection, and suppose also that u varies by order U0 over a distance in the ydirection. Thus, L is a measure of the boundary layer length and is a measure of the boundary layer thickness. In terms of derivatives, we can write u =O y u =O x and by applying (3), we get L, (4) U0 U0 L and , which means that the boundary layer is much thinner than it is long, which we already anticipated earlier in our definition of a boundary layer as a thin layer in where viscous effects are important. Applying these estimates for the magnitude of u and its derivatives, we can also estimate the magnitude of v. Start with (2): v u = =O y x U0 L . (5) We can approximate v/y in the boundary layer as follows v v(y = )  v(y = 0) . y But v(y = 0) = 0 (the boundary is at y = 0), so combining (5) and (6), we get v=O U0 L , (6) which quantifies our qualitative conclusion from the last lecture that v was much smaller than u. 1 Finally, we can use a similar analysis to estimate the thickness of the boundary layer as a function of the Reynolds number. Looking at (1), we want to ensure that the viscous term on the right side of the equation is comparable in magnitude to the terms on the left. Remember the bigger picture; away from the boundary, the flow is inviscid, so the Euler equation pertains, which means we only have to worry about the left side of (1) and the pressure term on the right. In the boundary layer, viscosity is important, so we additionally have to consider the viscous term on the right side of (1). Applying the same estimates as above, the terms on the left side are each of order 2 U0 /L, while the viscous term is of order U0 /2 . Setting these magnitudes equal, we get 2 , L U0 L U0 L
1/2 so = Re1/2 , (7) where we have defined the Reynolds number as Re = U0 L/. Compare this with (4). Our initial assumption, that u and v vary more quickly in y than in x, led to the condition L. From (7), we can see that that assumption is equivalent to requiring that the Reynolds number be large. 1 The boundary layer on a flat plate Let's look at the example of a boundary layer on a flat plate (Fig. 1). In this situation, a uniform, parallel, incompressible flow with u = U encounters a flat plate, infinitely thin, aligned with the flow direction. The plate is located at y = 0 and begins at x = 0. The viscosity of the fluid is Figure 1: A uniform flow encountering a flat plate. very small, so the flow is essentially inviscid except near the flat plate, where the noslip condition pertains, so u(y = 0) = 0 for x > 0. We will consider the steadystate situation, so time derivatives are zero. We can make one simplification to the boundary layer equations for this problem. Consider a streamline that is located just `outside' of the boundary layer. By `outside,' we mean that the streamline is in that part of the u(y) distribution where u U . Because this streamline is outside of the boundary layer, we can treat the fluid as ideal on the streamline. Since the flow is steady, we p can apply the (first) Bernoulli theorem on the streamline, i.e. + 1 u2 is constant on the streamline. 2 Now, remember from Lecture 23 that we can split the flow into two parts, the main (inviscid) flow and the boundary layer, and get solutions for those parts separately. For this problem, the main flow is uniform, and since the main flow is inviscid, there is no noslip condition. Thus the main flow remains uniform even in the presence of the plate, and in particular u is constant on streamlines outside the boundary layer. Bernoulli's theorem then tells us that the pressure, p, is constant on those streamlines. Since the streamlines in the main flow are parallel, we have dp/dx = 0 in the main flow. But we saw in the last lecture that p is basically a function only of x in the boundary layer. The pressure distribution inside the boundary layer, then, is the same as the pressure distribution 2 in the main flow just outside the boundary layer. For this problem, this means dp =0 dx in the boundary layer. With the pressure term out, then, the boundary layer equations become u u u 2u +v = 2 x y y u v + = 0. x y (8) (9) As we've done in other problems, we're going to look for a similarity solution for u. Acheson justifies this approach somewhat in the text; remember (from the review notes for the first exam) that if a similarity solution doesn't exist, it will become obvious in the course of the calculation. The similarity solution we seek is of the form u = U h(), where = y , g(x) (10) and rather than assume a particular form for the normalizing function g, we'll let it come out of the calculation. The first thing we'll do is write u (and v) in terms of a stream function, , which will guarantee that our solutiion automatically satisfies (9). Thus u= , v= . y x (11) By integrating (10) with respect to y, we then get = U g(x)
0 h(s) ds + k(x). (12) Remember now that lines of constant are streamlines. We want the plate to be a streamline (it would be bad if velocity vectors pointed into or out of the plate). We can arbitrarily set = 0 on the plate. The surface of the plate is located at y = 0, which also means = 0. So (12) becomes = U g(x)
0 h(s) ds, which can also be written as = U g(x)f (), where f (0) = 0, so f () is now the similarity function we'll be looking for. From this form of , and using (11), we have for the velocity components u = U f () and v = U (f  f ) 3 dg dx (14) (13) where the primes mean d/d. From (14), the velocity derivative terms in (8) are u y dg = Uf = Uf  2 x x g dx u 1 = Uf = Uf y y g 2u 1 = Uf 2 , y 2 g and inserting into (8), we get U 2 f f which gives U 2 f f which finally simplifies to f + U dg g f f = 0. dx (16) 1 dg 1 dg 1 dg f + U 2f f  U 2f f = U 2 , g dx g dx g dx g y dg dg f f + U 2 (f  f ) = U 2 , 2 dx g dx g g , (15) Remember that in writing (10), we were making the assumption that the velocity, which was originally a function of the two (independent) variables x and y, could instead be expressed as the function of a single independent variable, the similarity variable, , that combines x and y. Equation (10) led to (13), where the similarity function was f (). In order for the assumption of similarity to be valid, the differential equation for f (), (16), has to have as the only independent variable. But g(x), which is a function of the independent variable x, appears in (16). So for the similarity assumption to hold, we have to eliminate x as an independent variable. The way to do this is to set g dg/dx constant. The choice of this constant is arbitrary. We will pick g dg = dx U (17) since this choice will also eliminate and U from (16). (Why can we seemingly arbitrarily pick g(x) to satisfy the similarity assumptions? The answer is that we are free to choose g anyway we like in this case, we're looking for a similarity solution with the similarity variable y/g(x) as long as the assumptions we make are consistent and our solution for u and v satisfies the boundary layer equations. Basically, the boundary layer equations have a solution. All we're doing is trying to find that solution. To make the job of finding the solution easier, we express the problem in a more convenient fashion, where in this case `convenient' means in terms of just one independent variable. By doing this we don't affect the solution is any way. If the solution can't actually be written in the form we want, we will know because the equations we get will be inconsistent with our assumptions.) To solve for g(x), first observe that g dg d = dx dx 1 2 g . 2 Substituting this in (17) and integrating, we get 1 2 x g = +d 2 U 4 Figure 2: Initial stage of the boundary layer. where d is a constant. To specify the value of d, look at what happens when g = 0. If g = 0, then by (15), u/y becomes undefined. But it is quite reasonable for u/y to be undefined at x = 0. Consider Fig. 2. Immediately upstream of the plate, at x 0, the velocity profile u(y) is uniform. Then, starting at x = 0, the noslip condition holds on the plate surface, y = 0. For x 0, however, the boundary layer is still vanishingly thin, so u has to vary from 0 to U over a distance 0. But the derivative u/y is of order U/, so at x = 0 we can say that u/y is undefined. We therefore set g = 0 at x = 0, which gives us 1 2 x g = 2 U or g(x) = Finally, we have = (2U x)1/2 f (), where = which satisfies the differential equation f with boundary conditions f (0) = f (0) f (0) = 0 f () = 1. The first boundary condition comes from setting v = 0 at y = 0 in (14), the second comes from setting u = 0 at y = 0, and the third is from setting u = U outside the boundary layer. There is no analytical solution for this problem. Instead, the problem can be solved numerically. The resulting velocity profile, f () = u/U , is shown in Fig. 3. + f f = 0, y , (2x/U )1/2 ) 2x U
1/2 . 5 Figure 3: The flat plate boundary layer profile, u/U = f () (Acheson's Figure 8.8). 6 ME 563  Intermediate Fluid Dynamics  Su
Reading: Acheson, 8.38.5. Lecture 25  Boundary layers: high Re flow in a converging channel Consider a converging channel, with plane walls located at y = 0 and y = x tan , between which a viscous fluid flows at high Reynolds number. The two walls intersect at the origin, (x, y) = (0, 0). We will assume that there is a slit at the origin through which fluid can either be extracted from the channel (we'll call this the inflow case) or injected into the channel (outflow). First, we want Figure 1: Flow in a converging channel at high Re (Acheson's Fig. 8.9). to find the inviscid, mainstream flow. Given the geometry of the system, we will look for a radial velocity field, where u = ur (r, )^r . The continuity condition u = 0 becomes: e u = which gives (rur ) = 0 r and integrating, we get rur = const. so (calling the integration constant Q as Acheson does) ur =  Q . r (1) 1 1 u uz 1 (rur ) + + = (rur ) = 0 r r r z r r Negative Q corresponds to inflow (velocity vectors in the negative r direction), and positive Q corresponds to outflow (velocity vectors in the positive r direction). With this result, our problem comes down to finding the velocity profile in the boundary layer, y 0. (The boundary layer on the other wall will, of course, be symmetric with this one.) With the coordinates defined as in Fig. 1, we can express the flow for y 0 in terms of Cartesian coordinates pretty easily. In particular, our velocity field, which we had expressed generally as u = ur (r, )^r , e can be written in the boundary layer as u = ux (x, y)^x e 1 We now want to establish the proper equations for the boundary layer. The equation for the ucomponent of velocity is u u u 1 dP 2u +v = + 2 x y dx y (2) We saw earlier that the pressure gradient in the boundary layer, dP/dx, can be written in terms of the velocity, U (x), in the inviscid, main flow just outside of the boundary layer. In Cartesian coordinates, this main flow is U (x) =  Q . x (3) (Think carefully about what it means to write (3) in place of (1). The coordinate x is exactly the same as the coordinate r for y = 0, and x r for y 0. To express the main flow in terms of x, as in (3), the main flow solution must be valid near y = 0, meaning the flow must become inviscid for y 0. In other words, the boundary layer in this problem is required to be especially thin. As we saw in the last lecture, the ratio /L that describes the thickness of the boundary layer is proportional to Re1/2 . For this converging channel problem we are therefore assuming that the Reynolds number is high.) To get the pressure gradient, we use the Bernoulli theorem for steady flow of ideal fluid, which says that p/ + 1 u2 is constant along streamlines. Thus, for a streamline just outside the boundary 2 layer, p 1 =  U (x)2 + const. 2 (4) which holds because the flow is parallel in the region of interest. Taking the xderivative of this equation, we get 1 dp dU = U dx dx which, substituting into the boundary layer equation (2), gives u u 2u u Q2 +v = 3 + 2 x y x y (6) (5) As with the flat plate problem of Lecture 24, we'll look for a similarity solution to this equation, which will be of the form u = U (x) h(), where = y/g(x). (By using the absolute value of the main flow velocity, the difference between the inflow and outflow cases is described by the sign of the similarity function.) Using a stream function formulation to satisfy the continuity condition, u = 0, we get for the stream function: = U (x)g(x)f () =  so the velocity components are u= and v= 2 . x (9) Q = U f () =  f () y x (8) Q g(x)f () x (7) (We will hold off on the particular form of v in terms of U and f until later, when it will be simplified considerably.) From (8), we get boundary conditions on f . Outside the boundary layer, i.e. as (which is not the same thing as y !), we need u U (x). For the inflow case U is negative, and for the outflow case U is positive, so f () = 1 in the inflow case 1 in the outflow case (10) Let's look now at (6). If we plug in the expressions for u and v ((8) and (9)) and the appropriate derivatives in terms of the similarity variable, , then a similarity solution will only be possible if we end up with an ordinary differential equation for f () (meaning an equation that involves only f , the derivatives of f , and ). But consider the first term on the right side of (6), what was originally the pressure term. It is unchanged by the substitution of the similarity forms. In order for our resulting differential equation to be independent of x, then, we have to multiply each term in (6) by x3 to eliminate the xdependence of the pressure term, then we need to require that each of the remaining three terms be independent of x. The simplest of these three terms to analyze is the second term on the right side of (6) (no messy product rule to deal with). Using (8), we have 2u Q 2 f () Q =  =  2 2 y x y x y df d y =  Q f x y g Q f =  . x g2 From the preceding paragraph, we need x3 u to be independent of x, which means g(x) has y 2 to be proportional to x, for there to be a possible similarity solution. We are free to choose the proportionality constant in the expression for g, and we choose g(x) = (/Q)x 2 because this gives us a `cleaner' differential equation (fewer constants floating around). In particular, with this choice of g the right side of (6) becomes:  Q2 2u Q2 Q2 + 2 = 3  3 f x3 y x x = Q2 (1 + f ) x3 (11) using Q2 = Q2 . This g also makes computing the terms on the left side of (6) relatively easy, because the stream function in (7) becomes = which gives v= = x Q f =  Qf x x =  Qf =  Qf = 3  y dg g2 dx Q/ x /Q Qf (),  Q f . x Meanwhile, u is given by (8), and the derivatives of u are u Q Q f Q3/2 = f = =  2f y x y x g x u Q Q Q Q = 2f  f = 2f  f x x x x x x Substituting all this into the left side of (6), we have u u u Q +v = f x y x = Q2 f x3 Q Q f + 2f  2 x x
2  y dg g2 dx = Q Q f + 2 f . 2 x x (12) Q f x Q2 2 f . x3 Q3/2  2f x (13) +f f =  Setting (13) and (11) equal, and canceling the common Q2 /x3 , we get f
2 =1+f (14) as our ordinary differential equation for f (). The boundary conditions are given by (10) for . For = 0 (which means y = 0), the boundary condition is the noslip condition, u = 0, so from (8), f (0) = 0. (15) We can make (14) a secondorder, rather than thirdorder, differential equation by making the substitution F () = f (). Then (14) becomes F  F2 + 1 = 0 and the boundary conditions given by (10) and (15) become F (0) = 0, F () = and (17) (16) 1 in the inflow case 1 in the outflow case. 1 The inflow case, F () = 1 The analytical solution for this case makes use of a variety of tricks that are nonintuitive. The first thing we do is multiply (16) by F': F F  F 2F + F = which, integrated, gives 1 F 2
2 d d 1 F 2 2  d d 1 3 F 3 + dF =0 d 1  F 3 + F = const. 3 (18) To determine the value of the constant, we just need to know F and F at one value of . We know that F () = 1 already. To find F () = f (), look back at u/y, given by (12) u Q3/2 Q3/2 =  2f =  2F . y x x 4 (19) For , u U (x), so u U (x) =0 y y and by (19), F () = 0. Substituting for F () and F () into (18), we get 1 F 2
2 1 2  F3 + F = . 3 3 (20) Now the sleight of hand starts. We can factor (20) as F Then the substitution G2 = 2 + F turns out to be convenient. First, taking the derivative of (22), we get 2GG = F and substituting into (21), we get 4G2 G 2 = which gives 1 3 G = (3  G2 ) = 6 6 Our next substitution is G = 3 tanh H. We have sinh H = 1 G2 3 . (23) 2 (3  G2 )2 G2 3 (22)
2 2 = (1  F )2 (2 + F ). 3 (21) eH  eH and 2 eH + eH cosh H = 2 (24) giving d (sinh H) = cosh H d d (cosh H) = sinh H d sinh2 H  cosh2 H = 1 dH d dH d and tanh2 H  1 = sech2 H so d d (tanh H) = d d sinh H cosh H = 1 5 sinh2 H cosh2 H dH = sech2 H H . d Thus, with this substitution and (23), we have 3 G2 G =  3 sech2 H H = 1 3 6 3 = (1  tanh2 H) 6 3 =  sech2 H 6 which gives 1 H = , 2 and, integrating with respect to , H = + C 2 where C is the constant of integration. Finally, putting all of these substitutions together: F = G2  2 = 3 tanh2 H  2 = 3 tanh2 +C 2  2, (27) (26) (25) where we've used only the positive solution for H from (26), becauase tanh2 (x) = tanh2 (x). From (24), tanh() = 1, so F from (27) satisfies the boundary condition F () = 1. To satisfy the other boundary condition, F (0) = 0, we need tanh C = and it turns out that the solution to this is C 1.14. Profiles of F vs. are plotted in Fig. 2, for C = 1.14 and C = 1.14. Recall from (8) that u = U f = U F . For positive C, we get the sort of boundary layer depicted in Fig. 1, while for negative C, the boundary layer has a region of reversed flow near the wall. Both of these solutions are obviously valid mathematically; the question of which one would be observed in an actual flow depends on the details of how the flow is set up (more on this in Chapter 9). Also of interest: remember that all the way at the beginning of this analysis we said that the Reynolds number had to be high so that the boundary layer would be thin. Acheson shows in the text how this condition comes down to Q/ 1 in this case. 2 , 3 2 The outflow case, F () = 1 As mentioned in the text (Exercise 8.6, which is explained in the hints), there is no solution to the boundary layer equations for this case. (Actually, all we can conclude from the book's analysis is that there is no similarity solution to the boundary layer equations for this case.) 6 8 7 6 5 C = 1.14 C = 1.14 4 3 2 1 0 1 0.5 0 0.5 1 1.5 2 F
Figure 2: Profiles of F = u/U  vs. for the inflow case. 7 ME 563  Intermediate Fluid Dynamics  Su
Lecture 26  Boundary layers: equations for rotating flow
Reading: Acheson, 8.58.6. We'll consider now boundary layers in rotating flows. Consider a system with viscous fluid between two rigid boundaries located at z = 0 and z = L, and suppose the lower boundary rotates around the zaxis with angular velocity , and the upper boundary rotates around the zaxis with angular velocity (1 + ), where is small. The characteristic length scale in this flow is clearly L, and we can define a velocity scale by L, so we can write a Reynolds number for this system as Re = L2 . Analogous to our prior examples, if this Reynolds number is large, we expect that viscous effects will be important only in thin layers near the boundaries, and that the flow will be essentially inviscid in the rest of the space. In our previous examples, though, we were able to proceed by determining the inviscid main flow separately, without worrying about the effect of the boundary. For each of those cases, the boundary layer just allowed the given main flow, which was presumed to be driven externally, to match the noslip boundary condition. In this rotating flow system, the situation is different, because the flow is driven by the moving boundaries themselves. 1 Basic equations As usual, we'll start by writing the equations for the system. If a fluid is rotating with nearly  uniform angular velocity , it is helpful to consider a coordinate system that also rotates with  angular velocity . The NavierStokes equations in such a system are u + (u t and u =0 where u, in this case, is the fluid velocity relative to the moving coordinate system. The third term on the left side of (1) is the Coriolis `force,' and the fourth term on the left side of (1) is the centrifugal `force.' These terms only arise because the coordinate system is rotating and thus aren't really true forces. (For example, to understand the centrifugal `force,' consider what happens to (1) when u is zero.) If we apply the vector identity (which we won't bother to prove)   ( x) =  1  x2 , 2 1    )u + 2 u + ( x) =  P + 2 u (1) (note: remember that when Acheson writes F2 for some vector F, he means F2 = F F = F2 ) then we can rewrite (1) in the somewhat simpler form u + (u t 1 1   )u + 2 u =  p+ ( x)2 + 2 1 = pR + 2 u 1
2 u (2) where we define the `reduced pressure' as 1  pR = p   x2 . 2 One of our key assumptions here is that u will be small relative to the overall rotation of the system. For the system we talked about at the start of the lecture, this will be guaranteed because the two boundaries have almost the same angular velocity. If U is a characteristic magnitude of u, and L is a characteristic length scale of the flow, then we can write (u and  2 u = O(U )  where =  . The condition that u is small relative to the overall system rotation means that U L. We can then write (u )u O(U 2 /L) =O =  O(U ) 2 u or (u )u  2 u, U L 1, )u = O U2 L that is, the inertial term in (2) is much smaller than the Coriolis term. Thus, the NavierStokes  equations for a small velocity u relative a uniform rotation with angular velocity are u 1  +2 u =  p+ t u = 0.
2 u (3) (4) 2 Steady, inviscid main flow Let's consider now Cartesian coordinates that are fixed in the rotating frame of reference, and let the zaxis be parallel to the rotation axis. Thus  = (0, 0, ). We are interested first in the steady, inviscid flow solution for this rotating system. Let this flow be denoted by uI = (uI , vI , wI ) so that the Coriolis term in (3) is ex ey ez ^ ^ ^  u = 0 0 = vI ex + uI ey + (0)^z . ^ ^ e uI vI wI 2 Because the flow is steady and inviscid, we can ignore the u/t and Considering the terms that remain, we get 2vI =  1 pI x 1 pI 2uI =  y 1 pI 0= z 2u terms in (3). (5) (6) (7) from the x, y and zcomponents of (3) respectively, and we get uI vI wI + + =0 x y z (8) from the continuity equation. Now here's some interesting reasoning. From (7), pI is independent of z. This means, then, that pI /x and pI /y are also independent of z, so from (5) and (6), uI and vI are independent of z. Also from (5) and (6), we find vI 1 pI 1 2 pI = = y y 2 x 2 yx uI 1 pI 1 2 pI = = x x 2 y 2 xy and since the order of partial differentiation can be exchanged, vI uI = . y x Plugging this into (8), we get wI = 0, z which means that wI is independent of z. Combining this with the result for uI and vI above, we have that uI is independent of z. This is known as the TaylorProudman theorem, and can be stated: For a system in uniform rotation with small velocities u relative to the rotating coordinate frame, the steady solution uI of the inviscid equations of motion will not vary along the axis of rotation. 3 ME 563  Intermediate Fluid Dynamics  Su
Lecture 27  Boundary layers: Ekman boundary layers
Reading: Acheson, 8.58.6. In the last lecture we began to consider the problem of the flow between parallel boundaries (at z = 0 and z = L) that are rotating steadily about the zaxis at slightly different angular velocities. We were able to determine that the velocity of the inviscid, main flow is independent of z (a result known as the TaylorProudman theorem), which holds generally for coordinate systems in uniform rotation where the fluid velocity relative to those coordinates is small. Now let's consider our particular geometry, taking into account the effect of the walls. Assuming that the Reynolds number defined by Re = L2 is large, we expect that the flow will be inviscid in the region away from the boundaries, and will have thin viscous layers on the walls. Let's focus on the z = 0 boundary. The equations for the  velocity field (valid for small departures from uniform rotation at angular velocity ) were u 1  +2 u =  p+ t u = 0.
2 u (1) (2) We will apply our usual boundary layer assumption, that the velocity components vary much more rapidly in the direction normal to the boundary than in the other directions. In this case that means that the zderivatives of velocity are much larger than the x or yderivatives, so (1) becomes (noting also that the flow is steady) 1 2u  2 u =  p+ 2, z but the incompressibility condition remains the 2D boundary layer equations before. The components of (3) work out as (3) u = 0, using the same reasoning we used in deriving 2v =  1 p 2u + 2 x z 1 p 2v 2u =  + 2 y z 1 p 2w 0= + 2. z z (4) and the continuity equation written out is u v w + + = 0. x y z From (5), we conclude that w u, v (5) 1 since all the terms in (5) need to be of the same order, and the x and yderivative components are much smaller than the zcomponent. Then, we can also conclude (from 4) that p z meaning p p(x, y). (6) p p , x y Again, analogously to the reasoning we used with the 2D boundary layers, (6) means that the pressure field in the boundary layers is basically the same as the pressure in the inviscid, main flow. From the last lecture, we have the result 2vI =  1 p x 1 p 2uI =  y where the subscript I denotes the velocity component values in the inviscid main flow, and we've removed the subscript from p because the pressure field in the inviscid main flow is also the pressure field in the boundary layers. Putting this result into (4), we get 2(v  vI ) = 2u z 2 2v 2(u  uI ) = 2 . z (7) The book works through the solution of these equations (p. 282), which gives u = uI  ez (uI cos z + vI sin z ) v = vI  ez (vI cos z  uI sin z ), (8) (z = / z) which is known as the Ekman boundary layer. It turns out that just at the outside of this boundary layer (z ), the zcomponent of velocity is important. First write using the definition of z , then note 1/2 1/2 w w = , z z w w = = z z = = u v + x y vI uI + ez sin z  x y vI uI + ez sin z , x y uI vI + x y (1  ez cos z ) (9) where the last step is true because uI /x + vI /y = 0 from the last lecture. We can integrate (9) from z = 0, where w = 0, to z = to get w(z = ), the zcomponent of the velocity at the edge of the Ekman boundary layer. The result is wE (x, y) = 1 2 1/2 vI uI  x y , 2 but the term in parentheses is just the zcomponent of flow. We can denote this vorticity component as I , so wE (x, y) = 1 2 1/2 uI , the vorticity in the inviscid main I , . (10) Finally, what happens if the boundaries are rotating at different angular velocities, as with our problem? If we let B be the angular velocity of the bottom boundary relative to the rotating coordinate system, and T be the angular velocity of the top boundary relative to the rotating coordinate system, we get wE (x, y) = just outside the bottom Ekman layer, and wE (x, y) = 1/2 1/2 1 I  B 2 (11) 1 T  I 2 (12) just outside the top Ekman layer. Now let's focus back on the inviscid main flow. The velocity components given by (11) and (12) have to match the inviscid main flow just outside the bottom and top boundary layers, respectively. But we saw in the last lecture that the inviscid main flow velocity is independent of z, which means that (11) and (12) have to match each other. This gives I = B + T . For our problem (Acheson's Fig. 8.10) the coordinate system and the bottom boundary have angular velocity , and the top boundary has angular velocity (1 + ), so B = 0 and T = . Thus I = . (13) But remember from our discussion of the NavierStokes equations a while back that the vorticity is equal to double the local angular velocity. Thus, (13) tells us that the inviscid main flow between our rotating boundaries rotates at an angular velocity 1 , 2 that is, the main flow has angular velocity that is halfway between the angular velocities of the two boundaries. 3 ME 563  Intermediate Fluid Dynamics  Su
Lecture 28  Waves: the basics
Reading: Acheson, 3.13.3. We're going to go back in the book briefly to touch on the subject of waves. The topic of waves is very extensive and covers surface waves, interfacial waves, sound or pressure waves, shocks, solitons and a variety of other phenomena that fill textbooks all by themselves. Our goal is just to cover some fundamental ideas about waves, both to get a feel for the subject, and to learn enough to look at the area of flow instability. We will consider a simple harmonic surface wave, in which the position of the surface is defined by y = (x, t) = A cos(kx  t) = (), (A is a constant), where we define (x, t) = kx  t. The driving force behind the wave motion will be gravity. Sinusoidal waves are interesting because Fourier theory tells us that we can decompose any function into sinusoidal components of different wavelengths. So individual waves defined by (1) are the building blocks of general waves and disturbances. To determine the speed, c, of the wave, we plant ourselves on one particular point on the wave, i.e. fixed . From time t to t + t, that point on the wave will have moved from spatial position x to x + x, where x, t, x and t are related by (x, t) = (x + x, t + t) = k(x + x)  (t + t). The speed of the wave is then given by c= As we will see later, is given by = where g is the gravitational acceleration, so c= g . k (3) gk, (2) x = . t k (1) This tells us that waves with different k travel at different speeds. To illustrate the meaning of k, we will determine the wavelength, , of the wave defined by (1), where is defined as the spatial distance between successive peaks at a given time. Recognizing that in terms of , successive peaks are separated by = 2, we can write (x, t) + 2 = (x + , t), which gives us kx  t + 2 = k(x + )  t. 1 Thus the wavelength of the wave can be written = 2 . k (We can follow a similar analysis to find the period, T , of the wave, defined as the time between successive peaks as seen by an observer at a fixed location.) Going back to (3), we see that waves of longer wavelength (smaller k) travel faster. The property whereby individual waves of different wavelengths travel at different speeds is known as dispersion. Dispersion explains an interesting property of groups of waves. By `group of waves' (some authors use the term `wave packet') we mean a pattern of waves that occupies a finite area in space, not bounded by walls. Acheson points out that while an observer viewing a group of waves, which consists of a multitude of different Fourier components (each defined by (1) with its own k and ), may count N peaks in that group at any particular time, that observer would count more than N peaks as the group passed a fixed location. This is because the velocity of the group of waves is different from the velocity of the individual peaks. The speed of the whole group of waves is given by cg = which, using (2), becomes cg = 1 2 g 1 = c, k 2 d , dk using (3). So, the group of waves travels half as fast as the individual wave peaks themselves. 1 Surface waves on deep water While there's a lot of interesting background material on waves in 3.1 in the text, we're going to skip the stuff on sound waves and solitary waves and go directly to the example of surface waves on deep water. We will treat the water as an ideal fluid (negligible viscous forces). We are interested in twodimensional surface waves, so the velocity field is u = [u(x, y, t), v(x, y, t), 0], and we will assume that the flow is irrotational, so z = v u  . x y (4) The irrotationality of the flow is ensured if the fluid is initially at rest, i.e. initially has no vorticity; we saw all the way back in the first chapter that in twodimensional flow of ideal fluids subject to gravity, the vorticity of each fluid element is conserved. Because the fluid is irrotational, there exists a velocity potential, (x, y, t), such that u= , and v = . x y Also, because the fluid is incompressible, we have u= u v 2 2 + 2 = + = x y x2 y 2
2 = 0. (5) The position of the free surface will be given by y = (x, t), and it is the deformation of this surface that will generate the fluid motion we're interested in. We will require that fluid particles on the surface remain there (this implicitly eliminates the possibility of large disturbances, for example). If we define F (x, y, t) = y  (x, t), then the free surface is defined as the surface where F = 0. Since fluid elements on the free surface have to stay there, that means F is constant (and = 0) for those fluid elements, so in particular we can write DF F = + (u )F = 0 on the surface. (6) Dt t Using the following relations, F = t t F u = u x x F v = v, y (6) becomes +u =v t x on the surface. (7) 1.1 Smallamplitude approximations Now, we'll impose the requirement that (x, t), the free surface displacement, and the velocity components u and v and their derivatives be `small'. (The `smallness' will amount to requiring that the surface displacement be , the wavelength of the disturbances.) Then, in (7), we will throw away terms that are the products of small quantities (this is called linearizing), leaving us with v(x, , t) = , (8) t where the arguments for v establish that we're talking about v at the surface. Since is small, we can expand v(x, , t) around = 0 in a Taylor series. In general, for a function f (x), the Taylor expansion around x = 0 is f (x) = f (0) + x For our problem, v(x, , t) expands as v(x, , t) = v(x, 0, t) + v (x, 0, t) + , y df x3 d3 f x2 d2 f (0) + (0) + . (0) + dx 2! dx2 3! dx3 and once again we toss out anything involving products of small quantities, so (8) becomes v(x, 0, t) = (x, 0, t) = y t using v = /y, where we can impose the condition y = 0 instead of y = because of the linearization. 3 ME 563  Intermediate Fluid Dynamics  Su
Lecture 29  Waves: more basics
Reading: Acheson, 3.13.3. In the last lecture, we started to consider the problem of twodimensional surface waves on deep water, for which the fluid motion is driven by gravity. We defined the freesurface position as y = (x, t), with the implicit assumption that y = 0 is the position of the undisturbed free surface. Thus (x, t) can be interpreted as the free surface displacement due to the wave motion; we will also refer to (x, t) as te `disturbances.' The following relation pertains between the components of the fluid velocity, u and v, on the surface, and the freesurface position y = (x, t): +u =v t x on the surface. (1) Imposing the condition that the disturbances, , the velocity components, u and v, and their derivatives be small, we linearized (1) to get v(x, 0, t) = (x, 0, t) = , y t (2) ( is the velocity potential) which says that the `speed' at which the free surface moves in the ydirection is equal to the vcomponent of the fluid velocity at y = 0. 1 Bernoulli's equation for unsteady irrotational flow We can get another freesurface condition relating the disturbances to the velocity field by looking at the following form of Euler's equation for the velocity of an ideal fluid u +( t u) u =  p 1 2 + u + , 2 (3) where is the velocity potential. (Recall that g =  , and since g = g^y for the present e problem, we have = gy.) since the flow in this problem is irrotational, we have also = u=0 u= which, inserted in (3), gives = t p 1 2 + u + , 2 which can be written, interchanging the spatial and time derivative on the left side, t We can get rid of the gradient operator be interpreted = p 1 2 + u + . 2 (4) in (4) by integrating. Remember that the gradient can d ds = ds ds 1 ds = where ds is an arbitrary vector with infinitesimal magnitude ds. To integrate (4), then, we first take the dot product of both sides with some vector ds. The integral of both sides can then be written d which gives = t p 1 2 + u + + G(t), 2 (6) t = d p 1 2 + u + , 2 (5) where the arbitrary `constant' of integration G is allowed to be a function of t because the integration is over spatial variables only. Rearranging this, we get p 1 2 + + u + = G(t), t 2 which is Bernoulli's equation for unsteady, irrotational flow. By (7), the quantity p 1 2 + + u + t 2 (8) (7) is constant in space. Now let's look specifically at our problem. At the surface, y = (x, t), the pressure is equal to the atmospheric pressure, p0 , which can be treated as constant. Then (7) gives us 1 2 + u + v 2 + g = G(t), t 2 (9) where we've incorporated the pressure term into G, and used = gy. Now, we also note that G is arbitrary, because we can write the velocity potential as = s (x, y, t) + t (t), (10) where now t (t) is arbitrary, because we are only interested in the velocity potential through the expression u= = s , (11) and t (t) thus has no effect on the resulting u. In (9), then, we are free to choose G(t) = 0, giving us 1 2 + u + v 2 + g = 0. t 2 (12) We can now apply our smallamplitude approximations to this. Specifically, we are treating the disturbances, , and the velocity components, u and v, as small. The condition that u and v are small means that the spatial derivatives of are small, and we assume here that we can extend this to say that the time derivative of is also small. Linearizing (12) gives us (x, 0, t) + g(x, t) = 0, t (13) where this condition is imposed at y = 0 rather than y = by the same arguments we used in the last lecture to derive (2). 2 2 Dispersion relation
(x, t) = A cos(kx  t). (14) We are going to look for a solution to (2) and (13) of the form Look at (13) first. Plugging in (14), we get (x, 0, t) = gA cos(kx  t), t and integrating over t, we get (x, 0, t) = gA sin(kx  t). We have to be a little careful in interpreting this because in this expression is explicitly evaluated at a fixed y. What we can hypothesize is that has the form (x, y, t) = f (y) sin(kx  t). We can make the same hypothesis by considering (2). The potential has to satisfy Laplace's equation, i.e.
2 (15) = 2 2 + 2, x2 y so f (y)k2 sin(kx  t) + f (y) sin(kx  t) = 0 which gives f  k2 f = 0. The general solution of this is (as seen in any textbook on differential equations) f = Ceky + Deky . (16) We can simplify this by considering what happens at y  (remember that we're considering the problem of surface waves on deep water). Look at the ucomponent of velocity, which is (using (15)) u= = f (y)k cos(kx  t) = Ceky + Deky k cos(kx  t). x f = Ceky , and we can write (15) as (x, y, t) = Ceky sin(kx  t). (17) In order for this to be finite as y , we must have D = 0. Thus (16) becomes Now we want to get some idea of how the constants in (17) relate to each other. Inserting (17) into (2), we get Ckek(0) sin(kx  t) = A sin(kx  t) 3 or Ck = A, and inserting (17) into (13) we get Cek(0) cos(kx  t) + gA cos(kx  t) = 0 or C + gA = 0. Equation (18) tells us C = A/k, so (x, y, t) = Also, (18) and (19) together give us 2 = gk (21) A ky e sin(kx  t). k (20) (19) (18) which we anticipated earlier. This is the dispersion relation that tells us how and k relate, i.e. (remembering that the wave speed is given by c = /k and the wavelength is 2/k) how the wave speed varies with the wavelength. 3 Particle paths
u = Aeky cos(kx  t) v = Aeky sin(kx  t). From our result (20) for the velocity potential, the velocity components are We will assume that the instantaneous position of any particle, which we'll denote by (x , y ), is very close to its average position, (x, y). We can then write u= x Aeky cos(kx  t) t y v= Aeky sin(kx  t) t which integrates to x = Aeky sin(kx  t) y = eky cos(kx  t) These trajectories describe circles (Fig. 1). Two observations. The radii of the trajectories is Aeky , which decreases exponentially as we go deeper (more negative y), just like the fluid velocities do. The other observation is that near the peak of a disturbance, there is forward motion of the water, and near the troughs, there is rearward motion; you can see this very readily by watching a floating object (say, yourself) as it encounters a series of wavefronts. 4 Figure 1: Particle paths (figure from Faber, Fluid dynamics for physicists). 5 ME 563  Intermediate Fluid Dynamics  Su
Lecture 30  Instability: introduction
Reading: Acheson, 9.19.2. The subject of flow instability (also referred to as hydrodynamic instability or laminar flow instability) is our bridge to one of the most interesting problems, if not the most interesting problem, in all of the physical sciences; that problem being turbulence. The classic experiment in turbulence is that of Osborne Reynolds that Acheson describes in 9.1. Reynolds's experiment concerned flow in a tube of circular cross section. To mark the flow pattern, Reynolds used thin streams of dye introduced into the flow in the tube. Reynolds's observations seem simple; for low velocities the stream of dye ran down the tube basically undisturbed (emblematic of laminar flow), while for higher velocities, there would appear a point in the tube where the dye pattern suddenly erupted into seeming randomness and disorder (emblematic of turbulence). To quantify the point of transition between laminar and turbulent flow, Reynolds used (logically enough) the Reynolds number: Re = UL where U is the characteristic velocity of the flow, L is a characteristic length scale, and is the kinematic viscosity of the fluid. There are numerous striking aspects about the transition from laminar to turbulent flow. A couple of them are: As Reynolds observed (and as anyone can see by looking at the plume of smoke rising from a cigarette), the transition from laminar flow to turbulence is very abrupt, occurring at a particular value of the Reynolds number (the critical Reynolds number, Rec ). Flows are either laminar or turbulent, they can't be slightly laminar or kind of turbulent. This can be thought of as a symptom of the nonlinearity of the NavierStokes equations and of turbulence. Think of the inputs to a physical system as lying on the xaxis of a plot and of the outputs as being on the yaxis, and for fluid flow systems consider the Reynolds number as quantifying the input and consider the amount of disorder as the output. If fluid flows were linear, then, say, doubling the input should double the output. However this is clearly not true when the increase in Reynolds number goes across the critical value Rec . Maybe more interesting, there is really no reason why we might predict the existence of turbulence. Think of our solutions for viscous flow in a pipe or between parallel plates. There was nothing in our solutions for those flows that indicated that they should only be valid for low Reynolds numbers only. So where does turbulence come from? Basically, turbulence can be thought to arise because the physical world isn't ideal. Walls are never completely still or smooth, flow profiles are never completely uniform, boundary conditions are never completely steady (in much the same way that friction is never zero, and vacuums are never completely free of particles, and so on). If walls were completely steady, etc., then we could just bump up the flow velocity in a tube without changing the flow qualitatively. But because there are always irregularities, however small, we have to worry about what happens to those irregularities as we change our flow parameters. In the case of Reynolds's experiment, when he increased the Reynolds number, the infinitesimal irregularities in his wall smoothness, etc. eventually grew until they made the drastic qualitative change in his flow that he observed as turbulence. To treat this phenomenon mathematically, we enter the area of stability (or instability) theory. 1 1 Linear stability theory and KelvinHelmholtz instability We will consider the twodimensional flow system in which one deep layer of inviscid fluid, with density 2 , flows at uniform speed U over another deep layer of inviscid fluid, with density 1 > 2 , which is at rest (Fig. 1). We permit there to be surface tension on the interface between the fluids. Now suppose that the interface (given by y = (x, t)) between the two fluids is subject to a 2 1 U
y=0 Figure 1: Flow system illustrating KelvinHelmholtz instability. sinusoidal, travellingwave disturbance, so that (x, t) = A Re ei(kxt) (1) where we recall that ei = cos + i sin , and the `Re' indicates that we are taking the real part of the complex exponential. We are interested in the temporal stability problem for this flow system (as opposed to the spatial stability problem, more on that later). That is, given an arbitrary initial disturbance in space, how does that disturbance evolve in time? Applying Fourier theory, we can express a general disturbance of the interface between our two fluids as (x, t) =
 a(k) Re ei(kxt) dk. (2) Written in this form, a general disturbance is made up of individual disturbances, each of the form (1), with different values of the wavenumber, k (recall 2/k = , where is the wavelength). We will make the following definition: A flow system is called stable when subject to a general disturbance of the form (2) if the magnitude of the disturbance remains finite for all time. This requires that the magnitudes of individual disturbances of the form (1) remain finite for all time, for all values of k. Our next question is, how do we determine the magnitude of a disturbance as a function of time? The key will be to find the dispersion relation, = (k). For the temporal instability problem, we will require that the wavenumber k be real, but the angular frequency, , will in general be a complex number, which we'll write as = R + iI . What is the significance of the complex angular frequency? If we plug (3) into (1), we get (x, t) = A Re ei(kxR tiI t) = A Re eI t ei(kxR t) . 2 (4) (3) The magnitude of the imaginary part of the exponential is always bounded by 1, because it's just sines and cosines. However, the term eI t is a different story. If I > 0, then this term clearly grows exponentially as time evolves. So our stability definition from above can be recast in the following form: A system is temporally stable if I 0 for all values of k. Also, and this is a significant point, we will consider infinitesimal disturbances, so that we can linearize the necessary equations, as we did with the wave problem of the past couple of lectures; with this condition the procedure we're describing is called linear stability theory. Now let's turn to the particular flow system in Fig. 1. In the homework you will show that, with the disturbances treated as small and linearizing the appropriate equations, the following relation holds:
a b c (1 + 2 ) 2 22 U k + 2 U 2 k2  k(k2 T + (1  2 )g) = 0, (5) which is a quadratic equation in the variable . The roots (solutions) to (5), of course, have the form b b2  4ac = = (k), (6) 2a where the form (k) expresses that the solutions to (5) represent the dispersion relation for the system. From (6), it is readily seen that (a) there will be an imaginary component to if b2  4ac < 0, and that (b) in that case, one of the roots of (5) will have imaginary part greater than zero, and the other will have imaginary part less than zero. So for the system to be temporally stable, the solutions to (5) cannot have an imaginary part, i.e. we must have b2  4ac 0. For this problem, this condition is written b2  4ac = 42 U 2 k2  4(1 + 2 )[2 U 2 k2  k(kT + (1  2 )g)] 0. 2 Playing with the algebra, we have 42 U 2 k2 [1 ] 4(1 + 2 )k(k2 T + (1  2 )g
"D" (7) U2 kT + (1  2 ) g k 1 1 + 1 2 . (8) So, we know that a disturbance , given by (1), with wavenumber k, is stable if (8) holds. But for the flow system to be stable when subject to a general disturbance (2), we need this condition (8) to hold for all values of k. To ensure this, we just need to make sure that (8) holds for the particular value of k that minimizes the quantity D. The find that k, first take the derivative of D with respect to k, then set that equal to zero d d D= dk dk which holds when k satisfies k = g (1  2 ) . T 3 kT + (1  2 ) g k = T  (1  2 ) g k2 = 0, Then, D is given by D = 2 (1  2 )gT , and, putting this back into (8), we find that the system in Fig. 1 is temporally stable if U 2 2 (1  2 )gT 1 1 + 1 2 . (10) (9) Some observations; from (10), the system becomes more stable (i.e. the critical value of velocity Uc below which the system is stable increases) if the density difference 1  2 between the fluids increases, or if the surface tension T on the interface increases. Also, gravity is a stabilizing influence. (This treatment dealt with temporal stability, which is sometimes also called absolute stability. The problem of spatial stability involves following a particular wave packet as it evolves in space. In that case, we would allow our disturbances (1) to have real values of the angular velocity, , but complex values of wavenumber, k. Spatial stability is also called convective stability.) Figure 2: An interesting example of KelvinHelmholtzlike instability (from van Dyke, An Album of Fluid Motion). A 3 inch square, horizontal plate starts at rest and is suddenly accelerated upward. The photo is taken 6.5 milliseconds after the start of the motion, when the speed has just reached the maximum 24 ft/s. The pins on the right are spaced 1/4 inch apart; the lowermost pin represents the starting position of the plate. 4 ME 563  Intermediate Fluid Dynamics  Su
Lecture 31  Instability: more on KelvinHelmholtz, and thermal convection
Reading: Acheson, 9.29.3. To recap what we did with KelvinHelmholtz instability in the last lecture We defined our undisturbed, twodimensional flow system: a deep layer of inviscid fluid, density 2 , for y > 0, and a deep layer of inviscid fluid of density 1 for y < 0, with 1 > 2 . The interface is permitted to have surface tension. Initially, the upper fluid moves at uniform speed U in the xdirection and the lower fluid is at rest (Fig. 1). 2 1 U
y=0 Figure 1: Undisturbed flow for the KelvinHelmholtz problem. We allowed the interface to be subject to sinusoidal, travelingwave disturbances with infinitesimal magnitude. We then wanted to see how the initially, spatiallydefined disturbances evolved with time; specifically, do they increase or decrease in magnitude with time? We call a flow system stable to a particular disturbance if that disturbance dies away in magnitude after it's introduced. For this temporal stability analysis, we considered disturbances with real wavenumbers k but complex angular frequencies . We showed that disturbances grow exponentially with time, and thus destroy the undisturbed flow pattern (Fig. 1), when the imaginary part of the complex angular frequency is greater then zero. What we wanted, then, was to find the dispersion relation for the system, that is, = (k). This would tell us whether a disturbance of a given wavenumber either grows or decays with time. Linearizing the appropriate equations, then solving for (k), we were able to show that a given disturbance decays with time (i.e. doesn't grow and eventually destroy the initial flow pattern) when the following relation holds: U2 kT + (1  2 ) g k 1 1 + 1 2 . (1) Recognizing that an arbitrary initial disturbance is made up of a combination of sinusoidal disturbances, each with different wavenumber, we said that the flow system of Fig. 1 could be called temporally stable if disturbances of all wavenumbers decay with time. Working from (1), we showed that this was true when U 2 2 (1  2 )gT 1 1 + 1 2 . (2) Acheson's Fig. 9.4(b) illustrates this mathematically argument. 1 1 Thermal convection Consider a system where viscous fluid sits at rest between two plane rigid boundaries, located at z = 0 and z = d. We hold the two boundaries at different temperatures Tl and Tu , respectively, where Tl  Tu = T > 0; that is, the lower boundary is hotter than the upper boundary. The fluid nearer the lower boundary will therefore have a lower density, and will, due to buoyancy, have a tendency to rise. Because the system is confined, though, it is not immediately obvious how it will be possible for the fluid near the lower boundary to exchange places with the fluid near the upper boundary. It turns out that the fluid will remain in the intiail state of no motion until T exceeds a critical value. We can demonstrate this using linear stability theory. In this problem the variation of density is significant, so we won't be making our usual incompressible, constant density assumptions. With the NavierStokes equations as we've been using them, we have an equation describing mass conservation and one describing the balance of momentum. Here, we will need two more equations: an equation for the (internal) energy, and also an equation of state relating the fluid density to temperature and pressure. To simplify things a little, we'll confine our interest to the case where our viscous fluid is a liquid. (This problem is not trivial to solve, as you'll see from reading pp. 309311. We'll concern ourselves mainly with setting up the problem, then discussing the properties of the solution. In setting up the problem, we will basically be linearizing like crazy.) For liquids, density varies with both temperature and pressure, but the temperature dependence is small and the pressure dependence much smaller still. So we can include in our state equation only a linear temperature dependence, or = [1  (T  T )], (3) where T is just an arbitrary reference temperature, and (T = T ) = . Note: while T is arbitrary, (3) only holds if the temperatures we're considering are close to T . For the mass conservation, the general expression we found back in the first homework assignment D + Dt u = 0. Because the variations in for liquids are very small, we can basically ignore the D/Dt term, leaving us with u = 0, with which we're quite familiar. The momentum equation is also familiar: Du =  p + Dt
2 (4) u + g, (5) except we need to recognize now that isn't a constant. We will, however, assume that is a constant, independent of temperature. The energy equation we will use is T +u t T =
2 T, (6) which, as you'll recognize from the last homework, means that temperature (which we think of here equivalently to internal energy) is a conserved scalar quantity. This means that we're neglecting any changes of the total internal energy in the system, e.g. work terms or viscous dissipation. 2 We want to know first what the temperature distribution is in the undisturbed system where the fluid is at rest. Call this temperature distribution T0 (z). In (6), the terms on the left side are zero, leaving us with 0= d2 T0 . dz 2 (7) Integrating twice and applying the boundary conditions, T (z = d) = Tu and T (z = 0) = Tl , we get z z T0 (z) = Tl  (Tl  Tu ) = Tl  T. d d We can plug (8) into (3) to get the density distribution in the undisturbed fluid 0 (z) = [1  (T0 (z)  T )], (9) (8) and using (5), we can find the equation for the pressure distribution in the undisturbed fluid by tossing all the terms with u in them, leaving us with 0= dp0  0 (z)g. dz (10) At this point, we want to disturb the system, so T = T0 (z) + T1 = 0 (z) + 1 p = p0 (z) + p1 u = u1 (11) where the subscript 1 denotes small quantities that are functions of x, y, z and t. Now we plug (11) into our four equations (3), (4), (5) and (6), and linearize. For the state equation, we get 0 + 1 = [1  (T0 + T1  T )], and using (9), we have 1 = T1 . The continuity equation is simple: u1 = 0. For the momentum equation, we have (0 + 1 ) u1 + (u1 t )u1 =  (p0 + p1 ) + (0 + 1 )
2 (12) (13) u1 + (0 + 1 )g and throwing out products of small quantities, we leave ourselves with 0 Applying (10), we have 0 u1 =  p 1 + 0 t 3
2 u1 =  (p0 + p1 ) + 0 t 2 u1 + (0 + 1 )g. u1 + 0 g. Finally, we can also linearize (9). Recognizing that T0  T is a small quantity, we can replace 0 in the above equation with , giving us u1 =  p1 + t
2 u1 + 1 g. (14) In the energy equation (6), we recognize that T0 is not a function of time and use (7). This leaves us with T1 + u1 t (T0 + T1 ) = 2 T1 . T0 = dT0 /dz, giving (15) Again, we ignore products of small quantities in the second term, and use T1 dT0 + w1 = t dz
2 T1 . We'll look at the solutions to the system (12), (13), (14) and (15) in the next lecture. 4 ME 563  Intermediate Fluid Dynamics  Su
Lecture 32  Instability: thermal convection, and parallel shear flow
Reading: Acheson, 9.3, 9.5, 9.8. In the last lecture we considered a viscous fluid at rest between parallel plates at z = 0 and z = d, where the lower plate is held at temperature Tl and the upper at Tu , and Tl > Tu . The equations describing the evolution of small disturbances we found as: 1 = T1 u1 = 0 u1 =  p1 + 2 u1 + 1 g t T1 dT0 + w1 = 2 T1 . t dz Courtesy of some manipulation covered in 9.3 of the text (which makes use of differential equation analysis that's a little beyond us), we can first reduce the problem to a problem in the zcomponent of the velocity disturbance, w1 . In particular, the form of w1 is w1 (x, y, t) = W (z)f (x, y)est where f (x, y) has to satisfy 2f 2f + 2 + a2 f = 0, x2 y (2) (1) with a being some constant (which can take on any value). It is perfectly intuitive that we'd be interested in the zcomponent of the velocity; when we talk about stability of this system, the issue is whether the buoyancy introduced by the temperature gradient between the boundaries is sufficient to induce large scale motion of the warm fluid upward and the cold fluid downward. From the form of (1), it is clear that the magnitude of the disturbance w1 increases with time if s > 0. To simplify the problem somewhat Acheson gives the solution for the case where the boundaries are assumed to be stressfree (essentially ignoring viscosity at the boundaries). It turns out that s > 0 if the following holds: gT 1 > 2 d a N 2 2 a + d2
2 3 , (3) where N is any positive integer. In order for us to declare the system to be stable for a particular value of the temperature difference T , (3) needs to hold for all infinitesimal disturbances, i.e. all possible values of N and a. Basically, we want to find the minimum value of the right side of (3). This defines Tc , the critical value of the temperature difference. If T < Tc , then all infinitesimal disturbances decay with time. The value of N that minimizes the right side of (3) is obvious. Some simple calculus suffices to find the minimizing value of a. We get N = 1, and ac = . 2d Defining the Rayleigh number as Ra = gT d3 , 1 (5) (4) and using (4) in (3), we find that the critical Rayleigh number (above which the system is unstable and below which it's unstable)is Rac = 27 4 658. 4 If we don't impose the stressfree condition on the boundary, Acheson gives the critical Rayleigh number as Rac = 1708. Note what it means that the Rayleigh number governs the stability or instability. Increasing the viscosity or thermal conductivity stabilizes the system (raises the value of T before instability sets in), while increasing (which measures the response of the fluid density to temperature changes) and the spacing of the boundaries destabilizes the system. Acheson discusses the relation of these analytical results to experimental observations in the text. 1 Parallel shear flow: inviscid theory The last problem in instability we'll consider is that of the stability of parallel shear flow. Reynolds looked at this problem in some detail in the early 1880's, as part of the experimental effort that also involved the tube flow problem we discussed earlier; in the decade or so prior to Reynolds's efforts, Helmholtz, Kelvin and Rayleigh had looked at the problem from a theoretical standpoint. The theoretical work concentrated on the inviscid problem, which is much more tractable analytically than the viscous problem. The viscous analysis really became workable thanks to Orr and Sommerfeld in the early 1900's. We will get to the very influential OrrSommerfeld equation for viscous stability of parallel flows in the next lecture. First, the inviscid case. The problem setup is simple. Suppose we have two plane boundaries at y = L and y = L. We are only interested in twodimensional flow, so the velocity field has only uand vcomponents, and these components vary spatially only in x and y. The instability problem is inherently unsteady, of course, so we have to allow time variation of u and v. The governing equations are the Euler equations: u u u 1 p +u +v = t x y x v v v 1 p +u +v = t x y y u v + = 0. x y For our undisturbed flow conditions, we consider the parallel shear flow given by u0 = (U (y), 0, 0). (7) (6) Substituting (7) into the Euler equations, it is clear that (7) gives zero left hand side for each of the equations (6. Thus, the parallel shear flow (7) is a solution to (6) for any choice of the velocity profile U (y), provided that the pressure is given by p0 = const. (8) There are no boundary conditions on U (y), because the noslip condition doesn't pertain for inviscid flows. 2 Now we'll impose twodimensional disturbances to the undisturbed system. (small) disturbances by the subscript 1, we write u = (U (y) + u1 (x, y, t), v1 (x, y, t), 0) p = p0 + p1 (x, y, t). Plugging (9) into (6) gives us u1 d(U + u1 ) u1 1 p1 + (U + u1 ) + v1 = t x dy x v1 v1 1 p1 + (U + u1 ) = t x y u1 v1 + =0 x y and linearizing, we get u1 dU u1 1 p1 +U + v1 = t x dy x v1 v1 1 p1 +U = t x y u1 v1 + = 0. x y Denoting the (9) (10) Our experience with surface waves might suggest that we look for solutions for u1 , v1 and p1 that are sinusoidal in form. In fact, we do, except that because the equations (10) have coefficients of the derivatives that are functions of y (in particular, U = U (y)), we will look for solutions that have the following form (never mind how we know to do this; it's in the theory of differential equations) u1 = Re u(y)ei(kxt) ^ v1 = Re v (y)ei(kxt) ^ p1 = Re p(y)ei(kxt) ^ (11) where k is taken to be real, but the and the (^) quantities are complex, similar to what we did in considering the temporal stability of a given form of spatial disturbance: u, v , p C, ^ ^ ^ k Putting (11) into (10), we get 1 i(  U k)^ + U v =  ikp u ^ ^ 1 i(  U k)^ =  p v ^ ik^ + v = 0 u ^ where we've taken out the common exponential terms and where ( ) means d( )/dy. (12) (13) (14) , and C. 3 We can manipulate this to get an equation for v only. First, rewrite (14) as ^ u= ^ and plug this into (12), giving 1 1 ^ (  U k)^ + U v =  ikp. v ^ k Now differentiate this in y, and compare with (13):
cancels cancels 1 v ^ ik 1 1 U v + (  U k)^ + U v + U v =  ikp = ik[i(  U k)] = k(  U k), ^ ^ ^ v ^ k giving us 1 ^ (  U k)^ + U v = k(  U k). v k Multiply through by (  U k)/k, and we get v + ^ kU ^  k2 v = 0,  Uk (15) where we recognize that the boundary condition on v1 requires v = 0 at y = L, ^ (16) which expresses that even in inviscid flows, fluid can't go into or come out of boundaries. Rayleigh was able to take the problem defined by (15) and (16) to get a condition for the stability of this system without actually solving the equation (15). First, note that since v is ^ complex, we can write v = vR + i^I , ^ ^ v and we write the complex conjugate of v as ^ v = vR  i^I . ~ ^ v Then, v and its complex conjugate v satisfy ^ ~ v v = ^2 . ~^ v What Rayleigh did is to multiply (15) by v , then integrate across the domain in the ydirection: ~
L L v v dy + ~^
L L kU v  k2 ^2 dy = 0.  Uk (17) We attack the first integral by integrating by parts:
L use (16 L L v ~ L v v dy = [^ v ]L  ~^
L L ^ 2 dy =  v
L ^ 2 dy v 4 so (17) becomes
L L 
L ^  dy + v
L 2 kU v  k2 ^2 dy = 0.  Uk (18) Now, remember that since is complex, we can write = R + iI and using the complex conjugate and noting that U and k are real, we also have (  U k)(~  U k) =   U k2 , which allows us to write (  U k) = Then, plugging (19) into (18), we have
L L   U k2 . R  iI  U k (19) 
L ^  dy + v
L 2 (R  iI  U k)kU v  k2 ^2 dy = 0.   U k2 (20) We can see why we wrote (19) by noting that in (20), the real and imaginary parts can easily be separated, which wouldn't be the case if there were complex terms in the denominator of the fractional term. In (20), the real and imaginary parts of the left side have to be zero individually. Look at the imaginary part:
L I k
L U ^2 v dy = 0.   U k2 (21) Recall from our discussion of surface waves that disturbances like (11), where k is real and complex, grow exponentially in time if I > 0. If this is true, then the integral in (21) has to be zero. But everything in the integrand is necessarily positive everywhere in the domain of integration except for U . So for the integral to be zero, U either has to be zero everywhere, or it has to be positive in part of the domain and negative in another part. That is, if we consider profiles where U is continuous, then for I > 0 to hold, the velocity has to have an inflection point (where U = 0). This is Rayleigh's theorem: Rayleigh's Inflection Point Theorem. A necessary condition for the linear instability of an inviscid shear flow U (y) is that U (y) must change sign somewhere in the flow, i.e. the velocity profile U (y) must have an inflection point. This is only a necessary condition for linear instability, so just because a profile has an inflection point doesn't mean that particular parallel inviscid flow is unstable. But we can use Rayleigh's theorem to get a necessary condition for stability (meaning that disturbances don't grow in time); in particular, if a flow doesn't have an inflection point, then it is stable. 5 ME 563  Intermediate Fluid Dynamics  Su
Lecture 33  Instability: parallel shear flow, transition to turbulence
Reading: Acheson, 9.5, 9.8. In the last lecture we investigated the stability of parallel shear flows of ideal fluids, and proved Rayleigh's theorem, which says that a necessary condition for the linear instability of an inviscid shear flow U (y) is that the U (y) profile must have an inflection point. What happens when we consider viscous fluids? As before, suppose we have two plane boundaries at y = L and y = L, and we consider twodimensional flows, u = (u, v, 0), for which the NavierStokes equations take the following form: u u u 1 p +u +v = + t x y x v v v 1 p +u +v = + t x y y u v + = 0, x y u = 0 at y = L. As in the inviscid case, we will consider an undisturbed flow given by u0 = (U (y), 0, 0), and we will look at small perturbations to this flow given by u1 = (u1 (x, y, t), v1 (x, y, t), 0) . We'll do one thing differently from the viscous case; that is, we'll express the disturbances u1 in terms of a stream function , specifically u1 = and we look at ^ = Re (y)ei(kxt) ^ where k is real and and are complex. As before when we consider this type of periodic disturbance, the imaginary part of = R + iI determines the (temporal) stability of the problem. If I > 0 then the disturbance grows exponentially in time (i.e. the flow is unstable to the disturbance), if I < 0 the disturbance decays exponentially in time (the flow is stable to the disturbance), and if I = 0, the disturbance remains constant in magnitude (this condition is known as marginal stability. The boundary conditions on the disturbances are that the vcomponent at the boundaries needs to be zero, as before, and also that the noslip condition needs to be observed at the boundaries, since the fluid is viscous. Thus ^ v1 = 0 = 0 at y = L, and ^ u1 = 0 = 0 at y = L. 1 and v1 =  y x 2u 2u + 2 x2 y 2v 2v + 2 x2 y (1) and because the fluid is viscous, we also have the noslip condition (2) Figure 1: Stability curves for plane Poiseuille flow (a), from Acheson (Fig. 9.11). See the text for an explanation of curve (b). Inserting these forms for the velocity and pressure disturbances in the NavierStokes equations (1), we eventually get the OrrSommerfeld equation ^ i( ^ ^ ^ ^ ^  2k2 + k4 ) + (U k  )(  k2 )  U k = 0. (3) Acheson discusses the particular case of plane Poiseuille flow, in which the undisturbed flow between the boundaries is driven by a constant pressure gradient. In that case, the undisturbed velocity profile, U (y) is: U (y) = Umax 1  y2 L2 . (4) Acheson then presents the stability results for this case (the math itself is, again, beyond us) as a plot (given as Fig. 1) of wavenumber, k, vs. the Reynolds number, which is defined for this problem as Re = Umax L/. (5) The key observation from the plot is that is that there is a region of wavenumbers and Reynolds numbers for which the flow is unstable. In the inviscid case, the profile (4) is stable, because it has no inflection point (the second derivative is never zero). So in the viscous case, the existence of viscosity can actually destabilize the flow! This sort of runs counter to our intuition, which is that the the ability of a viscous flow to dissipate kinetic energy through friction ought to help damp out disturbances, and thus stabilize the flow. There are other observations we can make from Fig. 1. For one, we can see one sense in which viscosity does have a stabilizing effect. There is a critical value of the Reynolds number, which Acheson gives as Rec = Uc L/nu = 5772 for plane Poiseuille flow, below which the flow is stable to disturbances of all wavenumbers. Since the Reynolds number is defined by (5), increasing the viscosity, , stabilizes the flow in the sense that the critical velocity, Uc , is increased. Also, we can see from the figure that as we increase the Reynolds number, i.e. as viscosity becomes less important relative to inertial forces, the range of wavenumbers k for which the flow is unstable becomes smaller, eventually disappearing as Re . Thus, the inviscid stability solution from our previous lecture can be treated as the limiting case of the viscous stability solution for 0. 1 Turbulence The ideas of instability are understood to explain the onset of turbulence, at least in a certain sense. In the context of Reynolds's experiment of flow in a tube, we know that the flow is subject 2 to background disturbances. For low enough Reynolds numbers, these disturbances are unable to affect the flow structure, but as the Reynolds number increases, the disturbances grow and cause the flow structure to make the transition to turbulence. Of course, there are complicating factors. For starters, we have dealt only with infinitesimal disturbances, for which we can apply linear approximations; specifically, we have been making the assumption that products of perturbation terms are negligible relative to terms that are linear in perturbations. We then tracked the time evolution of these disturbances when applied to a particular background flow (the same principles could have been applied to tracking the spatial evolution of the disturbances). In the event that the amplitude of the disturbance grew with time, we said that the background flow was unstable relative to that disturbance, while in the event that the disturbance either decayed or remained of constant amplitude with time, we said that the flow was stable relative to the disturbance. We can see immediately that there are limitations to this view of things. For example, it is invalid to apply linear theory once the disturbance magnitudes are not infinitesimal, so technically, we can't say with certainty that a disturbance that initially grows in time will continue to do so for all time. As a disturbance grows, we start to have to consider nonlinear effects; and a priori, we have no way to predict if the disturbance will continue to grow with time when nonlinearity is considered. There has been a great deal of work in nonlinear instability theory (see, for example, the book by Drazin and Reid that's on reserve). However, the particular link between flow instability and turbulence is still not well understood. Actually there have been a multitude of attempts to `explain' turbulence theoretically, and still we search for a satisfactory `solution.' `Solution' appears in quotation marks because it is not, in fact, easy to say exactly what it would mean to `solve' turbulence. Consider the following analysis. Recognizing that turbulence is in some sense a random process superimposed on some orderly flow arrangement, it is customary in discussing turbulence to separate flow quantities in turbulence into their mean and fluctuating parts. For the velocity field, this means u(x, t) = u(x, t) + u (x, t), (6) where denotes some averaging process (which, depending on the flow of interest, could be in time or space, for now it doesn't matter which) and the primed term is the fluctuating term. Clearly, the fluctuating term describes the deviation of the actual value of the velocity vector from the mean value, at a particular point in space and time (x, t). The form (6) is known as the Reynolds decomposition. We point out for now some properties of the mean. The mean of a sum is equal to the sum of the means, or a+b = a + b , (where a and b are functions that can vary in space and time); taking the mean of a mean leaves it unchanged, or a = a, and the mean commutes with both temporal and spatial differentiation, so a a = and t t a a = . xj xj (7) What happens to the NavierStokes equations when we write the velocity in the decomposed form (6)? For the continuity condition, u = 0, we get, taking the average u = 3 u =0 (8) since averaging commutes with differentiation. We can also show easily from this that u = 0. (9) So (8) and (9) tell us that the mean and fluctuating parts of the velocity field each satisfy the continuity condition. More interesting is what happens to the equation for the momentum, which is in its original form Dui 1 p + = Dt xi Now write the material derivative as Dui ui ui (ui uj ) ui = = + uj + Dt t xj t xj where the second step uses the product rule and the continuity condition. Now take the mean of this: Dui Dt and substitute (6), which leads us to Dui Dt = u i uj ui ui + . + uj t xj xj (11) = ui u i uj , + t xj
2 ui . (10) Finally, taking the mean of (10) and substituting (11), we get ui ui 1 p = + + uj t xj xi
2 ui  u i uj xj . (12) Now suppose you're someone whose job is to compute turbulence numerically. As we will see, even when the mean terms vary slowly in space (and time), the fluctuating terms can vary over very small length and time scales. In fact, in general, it is unrealistic (in terms of computational resources) to compute all of the length scales in turbulence directly, because the required computational grid would have an unmanageable number of points. The problem, relative to the decomposed NavierStokes equation (12), becomes, how can we relate the term ui uj , xj which is known as the Reynolds stress, and which varies on the smallest length and time scales of the flow, to quantities like ui that we can compute on reasonably large scales? This problem (or some closely related form of it) has occupied a lot of people for a long time now. Ultimately, this problem (and, in fact, most of the difficulty with studying turbulent flows) arises from the wide disparity in the length scales of interest in turbulence. To get some feel for this aspect of turbulent flows, consider Figs. 2 and 3. 4 Figure 2: A mixing layer at high Reynolds number. The upper stream is moving at 100 m/s and the lower at 38 m/s, both from left to right in the image. Two things: (1) the transition to turbulence is evident on the left side of the image, where the initially smooth roller structures suddenly develop smallscale detail, and (2) the largescale organization of the flow is evident as the flow moves downstream, even though the flow has a lot of smallscale activity. Figure 3: The same flow arrangement as above, but at double the Reynolds number. There is more smallscale activity, but the largescale organization is not greatly affected. 5 ME 563  Intermediate Fluid Dynamics  Su
Lecture 34  Turbulence: basic ideas
We're not going to have time to discuss turbulence in any level of detail, so we'll just touch on a few significant ideas. When we speak of turbulence, certain things should come to mind: `Randomness' and disorder at small flow scales Coherence and order at larger flow scales Energy cascade: transfer of energy between length scales. Similarity and scaling of flow properties. We'll talk about thise things briefly in the time remaining. 1 Randomness and disorder  statistical tools When we speak of `randomness,' we have a particular definition in mind. Saying for example that the velocity, u, is random doesn't mean that it can take on any value; instead, it means that any time we run an experiment under the same conditions, u will take on a different value. It may seem strange that u can be a random variable when the time evolution of u is described by a deterministic system of equations, the NavierStokes equations. The NavierStokes equations have the form u = (something), t and, in principle, we should be able to integrate this equation in time, subject to initial and boundary conditions, to get the solution, u(x, t), for the velocity at all later times. So where does the randomness come from? Well, as we discussed in the context of stability theory, there will always be perturbations in the initial and boundary conditions of the flow walls are not perfectly still or flat, pressure is never exactly uniform, etc. Then on top of that, turbulent flows are particularly sensitive to perturbations (think of the exponential growth of disturbances in unstable flows). To describe the randomness of variables like the velocity component u, we use something called the probability density function, or PDF, which is itself definable in terms of the cumulative distribution function, or CDF. We start this discussion (which follows the notation of the book Turbulent Flows by S. B. Pope) with the following definition. If we run an experiment and monitor the velocity at a particular point, the probability, P , that u will be less than the value V is F (V ) = P (u < V ) (1) (V is a dummy variable called the samplespace variable, not the ycomponent of velocity). The function F is the CDF. There are three basic properties to the CDF: F () = 0, F () = 1, Vb > Va F (Vb ) F (Va ). (2) The first two of these hese tell us that u can't be less than , and that u has to be less than , respectively. To understand the last property, observe that, based on the definition (1), the 1 probability that u is between Va and Vb > Va is given by the probability that u < Vb minus the probability that u < Va , or P (Va < u < Vb ) = F (Vb )  F (Va ). Since the probability of an event can't be less than zero, this last expression means that F (Vb ) F (Va ). Now we define the probability density function (PDF), as the derivative of the CDF f (V ) = dF (V ) . dV (3) From the properties (2) of the CDF, we can derive the following properties for the PDF. First, f (V ) using the third property in (2). Second, F (V + dV )  F (v) 0 dV f (V )dV = f (V = )  f (V = ) = 1.
 Finally,
Vb P (Va < u < Vb ) = F (Vb )  F (Va ) =
Va f (V )dV. (4) If we let Vb and Va be separated by the infinitesimal interval dV in (4, we get P (V < u < V + dV ) = F (V + dV )  F (V ) = dF dV = f (V )dV, dV (5) which expresses that f (V ) is the probability per unit (whatever the random variable is, in this case velocity), which is why we call it the probability density function. The form of the PDF, f , completely characterizes the random variable u. In practice we are often interested in quantities like the mean. The mean of the random variable u is written u =
 V f (V )dV (6) i.e. it is the probabilityweighted average of all possible u values. In general, if G(u) is some function of u, we can write the mean of G as G(u) =
 G(V )f (V )dV. (7) It's worth thinking a little bit about how to interpret the mean defined by (6). u can be thought of as the expected value of u. For example, if we run our experiment N times with the same conditions, and measure u at the same spatial (and temporal) location, then as N gets very large, the average of our observations of u will approach u . (So it's not really a spatial average or a time average, although for certain flow conditions, it could be either). 2 The form of (7) makes it easy to see why the mean of a sum of random variables is equal to the sum of the means. Suppose G and H are both functions of the random variable u (which makes G and H also random variables). Then the mean of G + H is G(u) + H(u) =
 (G(V ) + H(V ))f (V )dV =
 G(V )f (V )dV +
 H(V )f (V )dV = G + H . Also, it is easy to see that the mean of a random variable isn't a random variable. Again, look at (7). Clearly, G(u) is not a function of u (and thus not a random variable), because the udependence has been integrated out. So we can write G(u) =
 G f (V )dV = G
 f (V )dV = G . In the last lecture we touched on the idea of the decomposing u using the mean: u= u +u. (8) One of the interesting things to look at in turbulence is the magnitude of u , which tells us in some sense just how random turbulence is. To quantify the magnitude of u , we might try just to look at its mean value. But if we take the mean of (8), we get u = u + u (using u = u ), which gives us u =0 i.e. the mean value of the fluctuating term in (8) is zero, which should be no surprise. To measure the magnitude of u , we instead use the variance, which is u2 =
 (V  u )2 f (V )dV. The square root of this is called the standard deviation, sometimes seen as u , and often called the `r.m.s.' (rootmeansquare, referring to how it's computed) of u. We can quantify the intensity of the turbulence by looking at the value of u . In a pipe flow, for example, we might compare u to the mean velocity in the pipe, u . If u is on the order of 10% of the mean flow velocity, that's pretty turbulent. Finally, using the definition of the mean, we can show easily that the mean commutes with differentiation. For time differentiation, first write du(t) u(t + dt)  u(t) = lim dt0 dt dt then write du(t) dt = u(t + dt)  u(t) dt0 dt 1 = lim ( u(t + dt)  u(t) ) dt0 dt u(t + dt)  u(t) = lim dt0 dt d u(t) = . dt lim 3 (9) The same argument pertains to spatial differentiation. It is important to keep in mind that in turbulence, as we saw in the figures in Lecture 33, the randomness is most evident at small scales. Therefore the fluctuation terms like u tend to vary on small flow scales, while mean terms like u tend to vary over much larger scales. This insight motivates virtually all theoretical approaches to understanding turbulence. 4 ME 563  Intermediate Fluid Dynamics  Su
Lecture 35  Turbulence: the energy cascade, scaling concepts
Now that we've seen that turbulent flows tend to be characterized simultaneously by largescale organization and a high degree of smallscale disorder, we might reasonably ask what sort of connection there might be between these scales. This question was first addressed by Richardson in the 1930s, and was subsequently taken up also by the mathematician Kolmogorov in the following decade. In this discussion we will be interested in very high Reynolds number turbulence, also called fullydeveloped turbulence, for which the notion of a wide disparity between the large and small scales is well justified. 1 The energy cascade Richardson conceived of turbulence as being composed of `eddies' of different sizes. What is an eddy? Though it isn't strictly correct, we can imagine eddies as being individual vortices of some measurable diameter, l. (As an alternative, Pope offers the following in his book Turbulent Flows; "An `eddy' eludes precise definition, but it is conceived to be a turbulent motion, localized within a region of size l, that is at least moderately coherent over this region.") Richardson then had the idea that large eddies are prone to break into smaller eddies, which break up into smaller eddies, etc. etc. In each breakup step, the larger eddy transfers its energy to the smaller ones, without dissipation, meaning that the energy transfer process in turbulence is inviscid. This makes perfect sense, since turbulence requires high Reynolds numbers, which in turn mean that viscous forces are negligible (with respect to inertial forces). To complete this picture, we need to say what happens at the ends of the process. Energy can't just be transferred from large scales to smaller scales forever. Somehow energy has to be put into the system, and somehow, energy has to be taken out. We know that energy dissipation is a viscous process, and we know that the energy transfer process is inviscid, so this tells us that energy is only taken out of the system at the smallest length scales, and that viscosity is important only at those scales. Meanwhile, the process is kicked off by the breakup of the largest eddies. These large eddies have the same characteristic velocity scale, U , and length scale, which we'll call , as the overall flow does. Using only U and , we can form a rate of energy transfer as U 3 / (to see why, verify the units). This represents the rate at which energy is put into the transfer process by the breakup of the large eddies. Figure 1 outlines this energy cascade process. We denote the rate at which energy is put into Energy dissipation, Energy transfer Energy production, P (length scale) Figure 1: The energy cascade process in fullydeveloped turbulence. The length scale represents a characteristic scale of energydissipating eddies we'll discuss this later. the system by P, and the rate at which energy is dissipated as . The energy input P determines . They are not, however, generally equal, though you might think they need to be by an energy balance. Keep in mind that the energy we're talking about is the energy contained in the turbulent motions only. There is, in addition, energy in the flow that is contained in the mean part of the 1 flow. The production term P really represents energy that is taken in some sense from the mean flow. This last point gives us insight into exactly how energy enters the system. The scales U and represent a level of mean shear, if we view as a distance over which the mean velocity can vary by U . Since U and are sufficient to define the rate of energy input to the system, we can say that energy is put into the system by the mean shear of the flow. There is an interesting analogy to the turbulent transfer process that exploits this observation. Consider a system of two parallel plates, between which sit a bunch of rocks with characteristic size roughly equal to the distance between the plates. Now suppose we allow the plates to move tangentially and in opposite directions to each other, generating a shear on the rocks. Eventually the big rocks break to form smaller ones, until finally you're left with dust that can't be made any smaller (Fig. 2). (Descriptions of this process are called breakage models.) Considering that energy Figure 2: The rock breakup process. is put into the system through the shearing of the rocks, one can ask questions like, what is the distribution of rock sizes at any point in time? 2 Kolmogorov's hypotheses Kolmogorov came to the subject of turbulence partly because of an interest in associated problems like the rock breakage problems. In 1941, he wrote a paper, barely more than two pages long, in a Russian mathematics journal that has gone on to spawn an amazing amount of work. That original paper and its followups addressed questions like: What is the size of the smallest length scales in turbulence, i.e. the ones that are responsible for dissipating energy? What's the relation between the size, l, of an eddy and the characteristic velocity, u(l), associated with it? In attacking these problems, Kolmogorov came up with three hypotheses. Kolmogorov's hypothesis of local isotropy. For high enough Reynolds number Re = U /, the smallscale turbulent eddies (with l ) are statistically isotropic. What this means is, while the large scales of the flow are permitted to look different in different directions (i.e., to be anisotropic), with properties defined by the mean flow field and the boundary conditions, the small scales will look the same in all directions. ('Local isotropy' means that the isotropy is confined to small scales.) Kolmogorov argued that, in proceeding down the energy cascade, all of the directional information in the largescale motions is lost. 2 Kolmogorov's first similarity hypothesis. For high enough Re, the statistics of the smallscale eddies (l ) take on a universal form determined by and , the energy dissipation rate. We can label these length scales as l (where we will define later). This is basically a stronger version of the isotropy hypothesis. Just as the directional information in the large scales is lost in the energy cascade process, so will all information about the geometry of the large scales, which are set by the mean flow and the boundary conditions. We are left with small scales that don't know anything about the large scales. This means, in other words, that the small scales are somehow universal and should be similar in every turbulent flow with high enough Reynolds number. Since energy is dissipated at these small scales, the important parameters are (1) the rate at which energy comes to the small scales, which (it turns out) is given by , and (2) the viscosity , since energy dissipation is a viscous process. Kolmogorov's second similarity hypothesis. For high enough Re, the statistics of motion in the range between and the dissipation scale (specifically, for length scales l such that l ) have a universal form that is determined by and is independent of . This hypothesis comes from the idea that if the Reynolds number is high enough, so the length scales and are really widely separated, there will be eddies of a size sufficiently smaller than that they will retain no memory of the large scale properties of the flow, yet they will be sufficiently larger than that the viscosity will still be unimportant. 3 Aside: shear flow scaling We'll take a moment to discuss an interesting way that we can use some of these ideas, along with some ideas we've seen earlier in the class. Let's look at an axisymmetric turbulent jet (Fig. 3). In this flow, there is a point source of momentum with a momentum flux J0 , which has units of uc(x) (x)
x r J0
Figure 3: An axisymmetric turbulent jet. momentum per unit time, in the xdirection. The jet fluid and ambient fluid both have density, 0 , and viscosity, 0 . If the momentum flux is high enough, the jet will be turbulent. Suppose now that we're only interested in the largescale, mean flow properties. In an earlier homework assignment we worked the problem of a laminar jet using the boundary layer equations. A key aspect of that solution was that the velocity profile was selfsimilar. We can apply the same principle here, and by the geometry of the problem we will also assume the the flow is axisymmetric. Specifically, we have u(x, r) = uc (x)f (), where = r , (x) 3 where u(x, r) is the mean axial velocity. (We're using cylindrical coordinates, but with x instead of z as the axial coordinate.) What we want to know is, how do uc and depend on x? Finding the flow width, , is easy. The major assumption is taken from the discussion of the previous section. In turbulent flows, the importance of viscosity is basically confined to the smallest scales, and since we're interested in the largescale properties, this means we can ignore viscosity. So we can say that depends on x, J0 and 0 . Considering the units of these: length, x length, J0 mass length mass , 0 . 2 time length3 The units of don't involve time, so can't depend on J0 ; but then also doesn't involve mass, so we can't use 0 either. Thus we just end up with = Cx. To find the xdependence of uc , we make use of the conservation of momentum, since there aren't any forces in our problem (why can't we use the conservation of mass?). The momentum flux through any arbitrary x = const. surface has to equal J0 , or
2 J0 = J(x)0
0 d
0 dr u2 r dr = 2u2 c
0 d f 2 = 2 u2 2 c
0 f 2 d. Looking at the last equality, the integral is constant, since f is a function only of . So u2 2 has c to be a constant, which means that uc = D . x 4 ME 563  Intermediate Fluid Dynamics  Su
Lecture 36  Turbulence: more on scaling, and the Reynolds stress
In the last class we looked at the Kolmogorov similarity hypotheses, which expressed the idea that the smaller scales of turbulence approach a universal state for high enough Reynolds number. (The Reynolds number we're interested in is Re = u0 l0 /, where u0 is a characteristic velocity of the flow large scale, and l0 is a characteristic largescale flow dimension.) That is, while the large scales will obviously be different, being set by the different flow boundary conditions, etc., the statistics of the small scales of the flow will be the same between flows if the Reynolds numbers are high enough. Also, if you look at the smallest turbulent scales, the turbulence at those scales will be isotropic, meaning that all of the spatial directions will look the same. Based on Kolmogorov's hypotheses, we can also extimate the size of the smallest turbulence length, time and velocity scales. From Kolmogorov's first similarity hypothesis, the universal form of the smallest turbulence scales is dependent only on the viscosity, and the rate of energy dissipation, . The units of these terms are = length2 velocity2 length2 , = = time time time3 (1) where we'll note that is really the specific rate of energy dissipation, i.e. the rate per unit mass. From these, we can form a unique length scale, , time scale, , and velocity scale, u , as: = 3 1/2 = u = ()1/4 .
1/4 (2) These are collective called the Kolmogorov scales, though if you see reference to `the Kolmogorov scale,' usually this means just the length scale, . From the scales (2), we can form the Reynolds number Re = u = 1, (3) which tells us that at the Kolmogorov scales, inertial effects and viscous effects are of roughly equal importance. This is consistent with the idea that the the scales (2) represent the smallest scales of turbulence, and with the notion of the energy cascade we also discussed earlier. Remember that we said that in a turbulent flow, energy is transferred without dissipation to successively smaller scales until, at the smallest scales, viscosity becomes important and energy is dissipated. The scales defined by (2) give us useful insight, but in reality, the dissipation is very difficult to measure experimentally. It would be more convenient if we could express the Kolmogorov scales in terms of largescale quantities. In particular, we'd like to express , and u in terms of the largescale length and velocity scales l0 and u0 . To do this, we remember from our prior discussion (Lecture 34) that the dissipation rate scales with l0 and u0 as u3 0 , l0 (4) which told us that the energy dissipated by a turbulent flow is proportional to the energy that's 1 put into the cascade process at the large scales. Using (4) in (2), we get the relations Re3/4 l0 u Re1/4 u0 Re1/2 (l0 /u0 ) (5) where Re = u0 l0 /. From (5), it is clear that as the Reynolds number increases, the disparity between the large and small scales increases. For example, if we take a given flow system and bump up the Reynolds number by increasing the velocity, the small length scale will get smaller even as the large scale l0 stays the same. 1 The Reynolds stress Recall from Lecture 34 that if we decompose the velocity field in turbulent flows as u i = ui + ui , and also write p = p +p, then the evolution equation for ui becomes ui ui 1 p = + + uj t xj xi
2 (6) (7) ui  u i uj xj , (8) which we compare to the NavierStokes equation for the raw ui : ui ui 1 p = + + uj t xj xi
2 ui . (9) The difference is obviously the last term in (8), which we'll write as Rij , where Rij = ui uj . xj (10) (Note: calling this term Rij is not standard notation.) It goes without saying that this term Rij /xj is significant. From (8) and (9), the raw field ui and the mean field ui would be identical except for Rij /xj . This means that all of the disorder that a turbulent flow manifests can be traced to Rij /xj . This is a pretty strong statement. (Consider Fig. 1, which shows the scalar concentration of a turbulent water jet. (The scalar in this case is a laser fluorescent dye in the jet fluid.) The mean scalar field is nice and smooth, basically Gaussian in profile. All of the disorder that you see is due to the scalar equivalent of the term Rij /xj .) What exactly does this Rij /xj represent? From our discussion of the NavierStokes equations, remember that (ignoring gravity) Dui 1 Tij , = Dt xj (11) 2 Figure 1: The scalar concentration of a turbulent water jet. where Tij was the stress tensor, given by Tij = pij + uj ui + xi xj . (12) Supposing that we're thinking about a small fluid element, what (11) and (12) tell us is that the acceleration of the fluid element is caused by pressure forces acting on the element, and by the transfer of momentum between the element and the surrounding fluid. Now let's look at (8). With (11) and (12) in mind, rewrite (8) as ~ D ui 1 Tij , = Dt xj where ~ Tij =  p ij + uj ui + xi xj  u i uj , (14) (13) ~ and Tij is, effectively, the stress tensor that's relevant for the mean velocity field, ui (why can't we call this Tij ?). If we're again considering some fluid element, (13) and (14) describe the mean acceleration of the fluid element; and, looking at the terms in (14), the first term tells us that there's acceleration due to the mean pressure, the second term tells us that there's acceleration due to transfer of mean momentum, and the third term tells us that there's acceleration due to momentum transfer by the fluctuating velocity field. This last point is very significant, because one of the major features of turbulent flows is that they show increased transport of quantities like momentum, scalar concentration etc. In other words, the disorder and smallscale activity in turbulent flows makes it easier for quantities like momentum and scalar concentration to mix. 3 ME 563  Intermediate Fluid Dynamics  Su
Lecture 37  Turbulence: some thoughts
The following is from A First Course in Turbulence, by Tennekes and Lumley: "Most flows occurring in nature and in engineering applications are turbulent. The boundary layer in the earth's atmosphere is turbulent (except possibly in very stable conditions); jet streams in the upper troposphere are turbulent; cumulus clouds are in turbulent motion. The water currents below the surface of the oceans are turbulent; the Gulf Stream is a turbulent walljet kind of flow. The photosphere of the sun and the photospheres of similar stars are in turbulent motion; interstellar gas clouds (gaseous nebulae) are turbulent; the wake of the earth in the solar wind is presumably a turbulent wake. Boundary layers growing on aircraft wings are turbulent. Most combustion processes involve turbulence and even depend on it; the flow of natural gas and oil in pipelines is turbulent. Chemical engineers use turbulence to mix and homogenize fluid mixtures and to accelerate chemical reaction rates in liquids or gases. The flow of water in rivers and canals is turbulent; the wakes of ships, cars, submarines and aircraft are in turbulent motion. The study of turbulence is clearly an interdisciplinary activity, which has a very wide range of applications. In fluid dynamics laminar flow is the exception, not the rule: one must have small dimensions and high viscosities to encounter laminar flow." We have discussed turbulence without offering any actual definition for it, becuase the turbulence is resistant to precise definition. Instead, we have mentioned the following properties. High Reynolds number. Turbulence occurs at high Reynolds number. Remember that we said that instability of laminar flow was one way that flows become turbulent; instability of a particular flow requires that the Reynolds number exceed come minimum value. The high Reynolds number also means that at the large scales of the flow, inertial forces dominate over viscous forces. Randomness and disorder. Turbulent flows exhibit a high degree of randomness and disorder, particularly at small scales. We describe these random fluctuations by decomposing the flow variables into mean and fluctuating parts, e.g. ui = ui + ui . Disparity in length scales. A wide range of length scales are relevant in turbulent flows. The spread between the largest and smallest length scales in a flow increases with increasing Reynolds number. This makes computation extremely difficult. Energy cascade. In a turbulent flow, energy generally gets transferred from the large scales to the small scales in an inviscid fashion, then gets dissipated at the small scales by the action of viscosity. The presence of dissipation in turbulent flows is significant; even though the large flow scales may be essentially inviscid, there is viscous dissipation of kinetic energy going on in the flow, at the small scales. Increased transport of momentum, scalars etc. The random fluctuations of turbulent flows provide another mechanism by which quantities can get transferred from one portion of a flow to another. In a laminar flow, transport of momentum mostly occurs through viscosity (diffusion), while in a turbulent flow, transport can occur through the random motion embodied by the fluctuating terms. Think of a parallel flow, with all the (mean) velocity vectors pointing in the same direction (say, the xdirection). In the laminar case, momentum only gets transferred between layers by viscous drag. In the turbulent case, however, there can be velocity fluctuations in the ydirection to carry momentum across the mean streamlines. 1 ME 563  Intermediate Fluid Dynamics  Su HW 1  due 4 February Problem 1  Cylindrical coordinates. Up to now we've focused on the Cartesian coordiate system. We will also make extensive use of cylindrical coordinates. In cylindrical coordinates, the xy plane is represented by the radial coordinate r and the angular coordinate , and the zaxis is unchanged, as in Fig. 1. The appendix in Acheson's book goes over the vector operations in cylindrical coordinates (Section A.6.) Here we will familiarize ourselves a little with the cylindrical coordinate system. (a) Compute the volume of a right circular cylinder of height h and radius R using cylindrical coordinates. Hint: the coordinate does not have dimensions of length. (b) The divergence theorem introduced in Lecture 1 is valid for any coordinate system. Using an analysis similar to that in Lecture 1, show that the divergence in cylindrical coordinates is as given in Section A.6 in the text. (c) Unlike Cartesian coordinates, where the unit vectors ex , ey ^ ^ and ez have the same orientation at all points in space, the ^ vectors er and e in cylindrical coordinates vary in space. ^ ^ Let s be the displacement of a point particle relative to the origin. Express s in cylindrical coordinates. Then, Figure 1: The cylindrical coordiexpress the velocity of the particle, ds/dt, and its accelernate system. ation, d2 s/dt2 , in cylindrical coordinates. Hint: you may want to express er and e in Cartesian coordinates as a ^ ^ convenience. Problem 2  Vector identities. (a) Using the list of vector identities in Appendix A in the text, show that when a flow is irrotational, the velocity can be expressed in Cartesian coordinates as u= , , x y z , where is a scalar function known as the velocity potential. (b) Define a vortex tube as a thin tube whose surface is parallel to the vorticity vector at every point. Is it possible for a vortex tube to come to an end in a flow, i.e. for the vorticity to go to zero along a surface that cuts through the vortex tube? Why or why not? Is the result limited to ideal fluids? Problem 3  Acheson, Exercise 1.1. Problem 4  Acheson, Exercise 1.3. Hints: refer to the discussion of Exercise 1.2, and when analyzing the term (u )u, be careful how you handle the unit vectors. Problem 5  Acheson, Exercise 1.6. Problem 6  Acheson, Exercise 1.7, but also show how the `Lagrangian' description given in the problem is derived. ME 563  Intermediate Fluid Dynamics  Su HW 2  due 13 February Problem 1  Vorticity in viscous fluids. (a) Derive the vorticity equation, Eq. (2.39) in the text, starting with Eq. (2.3). Recall that gravity is a conservative force. (b) Is the Reynolds number, Re = U L/, the proper parameter to describe the relative importance of the inertial and viscous terms in the vorticity equation? Show why or why not. Problem 2  Solutions of the NavierStokes equations. Consider the steady flow of an incompressible, viscous fluid between two parallel plates, where the the upper plate is moving at speed U in the xdirection, and the lower plate is fixed. The pressure is uniform everywhere in the fluid. (a) Solve for the velocity field u. Clearly state all assumptions. (b) Consider the point x = 0, y = h/2. What is the shear stress, y = u , and what is the magnitude of the viscous term y in the NavierStokes equation, 2 u, there? Are the results consistent with each other? Explain why or why not. (c) Show that the vorticity is nonzero in the fluid. Explain why this is when the flow system is clearly not rotating. Problem 3  Acheson, Exercise 2.1. Explain how you estimated the various parameters. Problem 4  Acheson, Exercise 2.4. Make sure you understand the relevant example in 2.3 in the text. Problem 5  Acheson, Exercise 2.6. Assume the pressure is constant between the plates. You will need to use separation of variables to solve the differential equation  and be careful with your arbitrary constants. Problem 6  Acheson, Exercise 2.14. Follow the reasoning in the book's hints carefully. Also, there are four boundary conditions, covering u and v at y = 0 and y . The condition on v at y is expressed as the similarity function f in terms of . You will need this to get the expression for the pressure in the hints, which you will need to get the differential equation for f given in the problem. ME 563  Intermediate Fluid Dynamics  Su HW 3  due 25 February Problem 1  Induced velocity fields. Suppose that we have two line vortices, each with circulation and each oriented in the positive zdirection, separated by a distance 2a in the x  y plane (see the figure). Assume that the two line vortices are fixed in space. Compute the velocity vectors at the points A through I indicated. Don't work any integrals unless you feel you need to. Figure 1: Geometry for Problem 1. Problem 2. The stream function. (a) Write the stream function for the inviscid stagnation point flow u = x, v = y. Compute the gradient of the stream function, . Show graphically how the direction of compares to the streamline directions. (b) Write the stream function for the viscous stagnation point flow given in Exercise 2.14. Problem 3. Acheson, Exercise 4.1. Note: a `simply connected' region is one whose border is connected, i.e. a simply connected region has no holes. Problem 4. Acheson, Exercise 5.1. For the circulation computation, you may want to use these trig identities: sin2 + cos2 = 1 sin 2 = 2 sin cos Problem 5. Acheson, Exercise 5.7. ME 563  Intermediate Fluid Dynamics  Su HW 4  due 15 March Problem 1  Isotropy of the stress tensor. Consider the example given at the end of Lecture 18. Ignore pressure. Compute the forces acting on a square fluid element in both coordinate systems, the original one (^1 , e2 ) and the one rotated by 45 . How do the forces expressed in the e ^ different coordinate systems relate? (This is similar to the book example starting on p. 210, but our velocity field is different, and I want you to work out all forces.) Problem 2. The Burgers vortex. (a) Acheson, Exercise 5.19. (b) Compute the terms in Eq. 5.38 in the book for the Burgers vortex of part (a) and the decaying line vortex given by Eq. 2.37 in the book. Explain what the terms in 5.38 are describing for both cases. Problem 3. Acheson, Exercise 5.18. Problem 4. Acheson, Exercise 6.3. Problem 5. Acheson, Exercise 6.4. Problem 6. Acheson, Exercise 6.6. ME 563  Intermediate Fluid Dynamics  Su HW 5  due 5 April Problem 1 Kinetic energy dissipation. Consider an infinitely long cylindrical tube with radius R, inside of which there is a flow of viscous Newtonian fluid. The flow is parallel, so u = u(r)^z . e For a fixed mass flow rate in the tube, m, find u(r). Assume that u(r) is of the form u(r) = A+Br . Do not use the NavierStokes equations. Problem 2 Scalar transport. Let (x, t) (`zeta') represent the concentration of a diffusive, conserved scalar. The units of are (1/volume). At any point, the flux of the scalar due to diffusion is q = D , where D is the diffusivity of the scalar. D has units (length2 /time) and is a constant. Derive the differential equation describing the time evolution of . As a starting point, you may want to apply control volume principles. Problem 3  Acheson, Exercise 8.2. Problem 4 (20 pts)  Acheson, Exercise 8.4. This is a long problem. Note the following. 1. When Acheson says one of the boundary conditions is "u 0 as we leave the jet," he means u 0 as y . 2. In the derivation of Eq. (8.65), keep in mind what symmetry the jet has. 3. In the hints you are told that g1/2 g has to be constant. Why? 4. In solving for f , these will be useful: sin2 + cos2 = 1. tanh x = d(sinh x) dx ex ex ex +ex d(cosh x) dx = cosh x, = sinh x. 5. Don't use Acheson's reasoning to estimate the width, , of the jet. Think instead about what the similarity solution means, and show graphically how the curves of u vs. y differ as x increases, and how they are similar. 6. You don't need to do anything with the material inside the square brackets on p. 295. 7. If you can't solve any part of the problem, just assume the result and pick up the next part (notice this problem is worth 20 points). ME 563  Intermediate Fluid Dynamics  Su HW 6  due 26 April Problem 1 Harmonic waves. (5 pts) Consider a sinusoidal disturbance defined by (x, t) = A sin(kx  t). (1) Define the period of this disturbance, T , as the time between the passage of adjacent peaks of the disturbance as seen by an observer at fixed x. Write (1) in terms of T and the wavelength, , instead of k and . Problem 2 Surface tension. (5 pts) Before doing this problem, read 3.4 up to Eq. 3.44, then look at the figure below y x T pu pl x T Consider the element of the free surface between the dots, with length x on the xaxis. The surface tension vectors at the endpoints have magnitude T (which is treated as a constant) and directions as indicated. Show that the vertical component of the surface tension at each endpoint is T /x, and then verify the first equation on p. 75 in the text by using Taylor expansions of the terms on the left side. Finally, give the expression relating the pressure above and below the surface element, pu and pl . On which side of the surface is the pressure higher? Problem 3 (10 pts)  Acheson, Exercise 3.2. In this problem you don't have to consider the surface tension between the interfaces. Problem 4 (10 pts)  Acheson, Exercise 3.6. Here you do need the surface tension, as covered in 3.4. You don't have to do the last part of the problem, deducing the condition for instability. 530.621: Graduate Fluid Mechanics I
Homework Set 4, due Wednesday Oct 15 with exam Reading: Panton: Chapter 6, 7.2. Panton does not discuss Bernoulli equation at length, so you should complement this material by reviewing your undergraduate course notes, read other books of the recommended list, etc... Problems: And*: 2. A nozzle is connected to a venturi and to a tank as indicated in the figure. Water flows in a steady fashion from the tank through the nozzle. The crosssectional areas at indicated sections are A1 =100 cm2 , A2 =10 cm2 and A3 =8 cm2 . The legs of the mercury manometer are filled with water and the reading is ha=300mm. (a) Calculate the volume flow rate Q, (b) the heigh hb .
3 From Panton: 6.4 Water hb venturi
1 2 nozzle r Free surface ha h 3. Not too close to the center of a whirlpool, the free surface takes on a shape qualitatively shown in the figure. For this region, it is known that particle paths are concentric cirles and that the velocity varies with radius r according to V = G/(2pr). Neglecting friction, determine the shape of the free surface in terms of r, G and fluid properties. What happens near r=0? 4. A scrubber is a device which removes pollutants from air by allowing the air to pass through a spray of water droplets. The pollutants are absorbed by the droplets, and the air that passes through is clean. An inventor designs the scubber shown in the figure. A highspeed water jet with velocity Vj and mass flow rate mw enters at section 1 and quickly breaks up into droplets which become dispersed in the tube. At station 2, the water droplets and the air are moving at essentially the same uniform speed Ve. The device automatically draws air from a source at atmospheric pressure. Calculate the mass flow rate of air, mair. Polluted, pa Polluted, pa Air Air Vj Water R Vj Air and, pa water spray Water
1 R
2 Ve h Polluted, pa Air 5. A hovercraft is supported by a highpressure curtain jet. Air is compressed from atmospheric pressure pa to stagnation pressure po by means of a compressor, and then issues from a narrow slot of width t that extens around the periphery of the base. The whole apparatus has a weight W and floats at a distance h above the ground. Assuming incompressible, inviscid, and steady flow, and that the mean pressure acting on the bottom surface is much smaller than po (pi <<po), and that t<<D and t<<h, compute h/D in terms of driving pressure popa, W, t and D.
pa Airflow pa 1 p0 Slot detail pi pi Ground D a h h pi Ground t 6. Consider a lawn sprinkler with radius R2 . The jet exits at an angle f to the plane of revolution. At station 1, the pressure is p1 and the velocity is negligible. In the jet discharge, the pressure p2 is atmospheric. Neglecting gravity and friction, and assuming that the sprinkler rotates at constant speed, (a) derive an expression for angular velocity in terms of p1 ,p2 ,r,f, and R2 . (b) Find upwards force Ry exerted by the lower flange on the upper flange near station 1, in terms of p1 ,p2 ,A1 ,A2 and f.
R2 w
2 2 f Water
1 1 7. Water flows in a 2D channel whose bottom surface is inclined locally at an angle a. The flow is steady and there is a free surface at the top. Local water velocity is V and local depth of channel is h. The angle a and its rate of change are small. Show that the slope of the free surface is given by [F2 /(F2 1)] tana, where F = V/(gh)1/2 is the socalled Froude number.
y V h a x Free surface 8. An incompressible inviscid fluid escapes through a plane 2D nozzle from a large reservoir where the pressure is po, to a region where the pressure is maintained at the pressure pa above the jet and at pb below the jet. The flow is steady and gravity effects are negligible. At some distance from the nozzle, the jet assumes an equilibrium configuration with a thickness t, a radius of curvature Ro at the center of the jet sheet. Show that the radius Ro is related to the volume flow rate Q and the pressure drops by Ro = [(1+N1/2)/ln(N1)] Q / [2 (popb) /r]1/2, where N = (popb)/(popa). p0 pa > pb y Q R0 pb t 9. Consider a bubble of highpressure gas expanding in an incompressible liquid in a spherically symmetric fashion. The gas is not soluble in the liquid, and the liquid does not evaporate into the gas. At any instant, R(t) is the radius of the bubble, pg is the gas pressure, r is the liquid density and p is the liquid pressure at a great distance from the bubble and s is the surface tension. Neglecting gravity, show that the radius of the bubble is described by the (RayleighPlesset) equation: R d2 R/dt2 + (3/2) [dR/dt]2 + 2s/(rR) = (pg p )/r. (Hint: First find the water velocity at any radial distance in terms of R(t)).
* Problems selected and adapted from `Advanced Fluid Mech.' course notes (MIT, 1984), courtesy of Prof. O. Knio. 530.621: Graduate Fluid Mechanics I
Homework Set 5, due Monday Nov. 3, 2003 Reading: Problems: Panton: Chapter 8 From Panton: 8.6, 8.15 And*: 3. The figure shows a "water cannon" to produce slugs of water at very high speed. Initially, with the nozzle exit temporarily closed by a light "cork", the cylinder and nozzle are loaded with water. Air at high pressure is the supplied behind the piston; the cork pops out and a jet of water is accelerated out of the nozzle at a volume flow rate Qe(t). The crosssectional area of this particular nozzle follows: A = Ao  (x/L) (Ao Ae). It is agreed to assume that (1) the flow is 1D, (2) water is incompressible, (3) viscous effects are negligible, (4) the crosssectional area of the cylinder is extremely large compared to that of the nozzle, (5) the exit pressure pe=patm, (6) h<<L and b<<h, and (7) the driving pressure difference pd(t) is known, increasing function of time. (a) Calculate (pd  po ) by making appropriate approximations. (b) Derive the governing equation for Qe, in terms of (pdpo ), L, Ao , Ae, and relevant fluid properties. (c) Assuming Qe(t) has been determined, find the total tensile force F(t) holding the nozzle to the cylinder. (d) Express F(t) and its dependence on relevant parameters in dimensionless form. Piston Cylinder, Ac Ac Water nozzle Qe(t) Ae, x pe Air pd (t) A0 , p0 b h L 4. The velocity of propagation c of a surface wave in deep water is assumed to depend on the wavelength l, the fluid density r, acceleration of gravity g and the surface tension s. Viscous effects can be ignored. (a) Using dimensional analysis, obtain a functional relationship for the dimensionless propagation velocity. (b) Use dimensional analysis to obtain an expression for c assuming that surface tension forces are negligible compared to gravitational forces. (c) Same, but assuming gravitational forces are negligible compared to surface tension forces. (d) Formulate a dimensionless group which indicates relative importance of gravity and surface tension effects. (e) On the basis of (bd), construct a schematic chart showing the variation of c with l, all other parameters remaining fixed. 5. It is proposed to evaluate characteristics of a journal bearing by testing a similar model 1/6th the size. (a) Assuming inertial forces in the oil are negligible (density is irrelevant), find the dimensionless form of an equation connecting the friction torque T, the load supported W, the diameter of the bearing D, the angular speed w, the oil viscosity m, and the clearance h (difference in diameter of bearing and shaft). (b) If the prototype is designed to carry a load of 15 tons at 120 rpm, and the model is to be tested at 180 rpm using the same oil as the propotype, what should be the model load? (c) With the load of part (b), what would be the relationship between the frictional torques in the model and prototype? 6. The vibrations of a long thin cantilevered rod of diameter d and length L are damped because of the drag exerted by the surrounding fluid which has density r, viscosity m, and bulk modulus b=r (dp/dr). To characterize the decay of periodic motion, we are interested in how the dissipation to the fluid, PD, depends on the oscillation frequency f, the amplitude of the motion at a particular instant, z, and any pertinent fluid properties. (a) Define a dimensionless power dissipation and find the dimensionless groups it depends on. (b) The independent groups you obtained in (a) may or may not have a clear physical meaning. Combine them as necessary to form a set of dimensionless groups which most clearly portray the important physics of the problem. Give a physical interpretation of each. L d 7. A pair of plates of length L is hinged at O such that the angle q(t) contained between the pates may vary in time. An incompressible fluid of density r and viscosity m fills the triangular gap between the plates and the space surrounding the plates. The angle q is always small, so that sinq ~q and cosq~1. When the angular velocity is positive, fluid is pulled into the gap from the surroundings, and viceversa. The plates are very wide in the direction normal to the paper, and the flow may be considered as twodimensional. It is agreed that p/y is negligible. (a) Express the pressure difference p(x)p(L) in terms of dimensionless parameters. It is assumed that q(t) is a known function of time. Explain the reason for including each independent variable you choose. In particular, explain whether or not the successive time derivatives dq/dt, d2 q/dt2 etc.. should be included. (b) Let U(x,t) be the average velocity in the xdirection, at the location x and time t. Solve (dimensionally) for U as function of q, dq/dt and x. (c) Show that the orderofmagnitude estimate v/u ~ O(q), where v is the local vertical velocity and u the local horizontal velocity in the xdirection, holds. (d) From an orderofmagnitude analysis of the NavierStokes equation in the xdirection, determine the criteria which would allow one to consider the local flow as inertiafree. (e) What are the criteria for the entire flow to be inertiafree? Show that even if the entire flow cannot be considered inertiafree, at least some part of it may be so considered. y q (t) x x L 8. An engineer tests an aircraft of length 20m and crosssectional area 6m2. He tests a scale model of length 1.2m in a water tunnel. It is found that the drag force of the model at cruising conditions depends on the flow speed according to D = 20V + 0.4 V2 (N), where V is in m/s. It the test, the model stalls when V decreases to 10 m/s. (a) Predict the stalling speed of the real aircraft. (b) Obtain an expression for the drag force on the real aircraft in terms of its flying speed. (c) Calculate the power output required by the real aircraft to maintain a steady cruising speed of 100 m/s.
* Problems selected and adapted from `Advanced Fluid Mech.' course notes (MIT, 1984), courtesy of Prof. O. Knio. 530.621: Graduate Fluid Mechanics I
Homework Set 6, due Monday Nov 10. Reading: Problems: 7.4). And: 4. A piston moves in a cylinder whose dimensions are shown below. The piston moves in at a speed U. Assume that the gap is much smaller than the piston length. Solve for the velocity profile in the gap. Find the force that must be applied to the piston. Panton: Chapter 7 and 11 (in depth reading) From Panton: This is like a single problem: 7.5,7.6,7.7,7.8 (they all refer to 7.5, not h (<<D) m D U L 5.* A horizontal water stream covered with a very thin layer of oil impinges on a flat plate which is inclined at an angle q with the horizontal direction. Upsteam, the water depth is hw1 and the oil film has depth ho 1 (ho 1 << hw1). The two fluids move at the same speed V1 . The water density is larger than that of oil (rw > ro ). The plate elevation is given by h(x) = h1  qx. The oil thickness is ho (x) << h(x). Viscous forces dominate in the oil layer, but the water stream may be considered inviscid throughout. We may assume 2D, steady flow, neglect gravity and endeffects near edges of plate. (a) If the water flow is inviscid, what force drives the thin film forward? (b) Find the pressure distribution p(x) in the water (c) Find the oil layer thickness ho (x) on the plate. oil water V 1 x h(x)
ho(x) q * Problem adapted from `Advanced Fluid Mech.' course notes (MIT, 1984), courtesy of Prof. O. Knio. 530.621: Graduate Fluid Mechanics I
Homework Set 7, due Wednesday Nov. 19 Reading: Problems: 1. Find the pressure difference p1 p2 required to drive a volume flowrate Q across the axisymmetric channel with parabolic shape as shown. Use the lubrication approximation. Express limits of validity of the approach in terms of given variables. Panton: Chapter 21 p1 m z=0 z r(z) = Ro+ a z 2 p2 2. (previous midterm #2 problem): Consider twodimensional laminar flow of Newtonian gas (density r1 , viscosity m1 ) in a channel, that has a groove with Newtonian liquid. Gas and liquid are separated by a horizontal interface. The upper wall of the channel moves from left to right at speed U. The lower wall contains a very long (in x and z directions) and narrow (in ydirection) groove filled with a liquid of density r2 , viscosity m2 . H and h are of the same order of magnitude. In the central region A away from the ends of the groove, the streamlines are essentially horizontal (v=0). y U gas r1 , m1 liquid x A>> h L >>> h r 2, m2 h H (a) (5 pt) Assuming that the upper wall has been in motion for a very long time, sketch the streamlines in the gas and in the liquid region. (b) (10 pt) Continue assuming that in region A, the flow can be considered as horizontal (v=0) and steady, write down the equations governing the motion in the gas (velocity ug (y)) and in the liquid (velocity uL(y)). In this problem, you may assume that, due to very high surface tension, the interface remains horizontal as shown above, even when a pressure difference develops between liquid and gas, and between left and right edges of the groove. (c) (5pt) Write down all necessary boundary conditions. (d) (5pt) Write down the appropriate condition from which you can solve for the horizontal pressure gradient p/x that exists in the liquid groove. (e) (5pt) Solve for ug (y), uL(y), and p/x. (f) (5pt) If at t<0 the upper wall and fluids were stationary, and the upper wall impulsively starts moving at speed U=0.1 m/s at t=0, estimate (no exact solution required) the order of magnitude of the time required before the gas near the liquid surface (y=h) achieves a significant fraction of its final velocity. Other numerical parameters are as follows: h=1cm, H=2cm, n1 =1 104 m2 /s. (g) (5pt) Still assuming fullydeveloped conditions in the xdirection and a flat interface at y=h, clearly write out (do not solve) the complete mathematical formulation which would allow you to solve for unsteady velocity profiles in both fluids. (general notation, do not use numerical parameters of part (f).) 3. (previous exam problem) A piston with a conical shape (segment of a cone) is pushed into a liquidfilled cylindrical container as shown in the figure. The incompressible fluid inside region C is squeezed through a narrow passage of linearly decreasing height. We are interested in calculating the force required to push the piston at a constant speed U. The analysis will be based on the following assumptions: (i) The angle a (approximately equal to (h1 h2 )/L) is very small (but not negligible). (ii) The gap height h1 is much smaller than D1 . Also, h1 << Lc. Lc L Detail of gap: h1 h(x) h2 a
D2 region C pc fluid m, r D1 D D(x) patm = 0 (say) U, F (?) (a) (4p) Without solving partial differential equations, calculate the average horizontal velocity (in the gap) of the fluid which is being squeezed out, in terms of U, D1 , h1 . (b) (4p) Give an estimate of the ratio of inertia and viscous terms in the momentum (NS) equation in the horizontal direction for the flow in the gap. In terms of a, U, D1 , h1 and fluid properties, under what condition can we neglect the inertia terms in the equations? Now, assume that this condition is met. Therefore, we can assume mixed Poiseulle+Couette flow at any location x in the gap, but with a varying (unknown) pressure gradient (lubrication approximation). (c) (4p) Is it important to keep the moving piston wall (Couette part) when solving for the velocity profile inside the gap? To answer this question, construct a dimensionless parameter b (velocity U divided by mean velocity in gap). Consider the limit when h/D1 > 0. (d) (12p) Neglecting the Couette part in the flow inside the gap, find the pressure in region C. (I am assuming you are able to integrate (a+bx)2). (e) (12p) Find the force F. (I am assuming you are able to integrate (a+bx)1). There are three basic contributions to the force. Compute each, and determine which is dominant (assume L~O(D)). When integrating, you may also assume that D(x) ~ D1 . 4. Solve axisymmetric Stokes flow of a liquid spherical droplet of viscosity mi and radius a, in another external fluid with viscosity me. Some results are exhibited in Panton, but you are asked to derive and show all details, following approach used in class. (a) Using the general expression derived in class (which is valid for both external flow and inetrnal circulation) you must now solve for 8 coefficients Ai, Ae, Bi, Be, ...etc.. using appropriate boundary conditions. Obtain the velocity field ur and uq both inside and outside the sphere in terms of the parameter g=me/mi . (b) Find the pressure distribution at r=a. (c) Find the distribution of tangential and normal viscous stresses (note: now the normal viscous stress need not be zero) (d) Integrate to find the drag force in terms of U, me, a and g. (e) Confirm that you recover the case of a solid sphere if g > 0. Obtain an expression for drag force for the case of a bubble in a liquid of viscosity me, by letting me/mi > infinity. 5. Using any computer plotting package, plot the steamlines for Oseen's solution of low Renumber flow over a sphere in a frame of reference moving with the sphere. Consider Re=0.001, Re=0.1 and Re=1. Comment on qualitative differences (e.g. frontaft symmetry, etc..). Homework Set 8 Due by 5pm Monday Dec 8 (in mailbox box Latrobe 127) Reading: Problems: Panton: Chapter 16.4, all of chapter 20. From Panton: 20.3, 20.6. 530.621: Graduate Fluid Mechanics I 3. Consider the steady laminar incompressible flow in entry of a 2D channel of constant crosssection. Boundary layers on the walls begin at x=0 and it is agreed to assume that the velocity profile in the boundary layer is linear with distance from the wall. Derive a formula relating the pressure drop coefficient [p(x)  p1 ] /(rV2) (h/x) to the length and breadth Reynolds numbers, defined respectively as rVx/m and rVh/m. p1 V p(x) h d x 4. A 2D channel has a convergingdiverging crosssection, whose height is given by h(x) = h0 for x<0 and x>L, and h(x)=h0 [1+ a (1cos(2px/L) )]1 for 0<x<L. The bottom wall is straight and only the top wall has the cosshape. Inlet velocity is U0. To a first approximation, neglect the effect of displacement thickness on the outer flow. A laminar boundary layer has been developing on the bottom surface. Assume that at x=0, the displacement boundarylayer thickness is d*(0). (a) Using the momentum integral formulation with an appropriate parabolic testprofile between y=0 and y=d(x), write down an equation for d(x). (b) Using the values U0=1 m/s, n=105 m2/s, h0=0.03m, a=0.2, L=0.3m and d*(0)=1mm, solve this ODE numerically (by writing a short program). Integrate between x=0 and x=L using a simple firstorder Euler integration: dd/dx = g(d,x) leading to di+1 = di + (xi+1  xi) g(di,xi), starting from the initial (boundary) condition at i=0, up to i=100. xi=(i1) L /100. Plot the resulting d*(x) in dimensional units. ME 563  Intermediate Fluid Dynamics  Su Exam 1  1 March 2002 50 points Problem 1  Short answers. (Total 20 pts) Name (a) (5 pts) Consider the velocity field defined in Cartesian coordinates as: u(x, y, z) = 1 1 1 x2 y  yz 2 ,  xy 2  xz 2 , xyz . 2 2 2 Is it possible to write a velocity potential, , for this flow? Why or why not? (Even if it is possible, you don't need to find .) (b) (5 pts) Draw the streamline patterns for the following stream functions: (1) = 2r 3 , (2) = 4/ r. Think before you compute any derivatives. 1 (c) (10 pts) Consider the discussion of D'Alembert's paradox in 4.13 in the book, which makes use of the figure below. A body sits in a flow of ideal fluid in a long, straight channel with uniform crosssection the paradox is that it feels no drag force, D. The discussion in the text hinges on the assumptions that the flow is uniform and has speed U0 both far upstream and far downstream of the body. The assumption that the flow is uniform far upstream is easily justified as part of the definition of the problem. It is also easy to justify that the flow is parallel (i.e. u = (u, 0, 0)) far upstream and downstream, since we've assumed that the channel has a uniform crosssection. However, can you think of a way to justify the assumption that the flow is uniform far downstream of the body? Figure 1: Sketch for illustrating D'Alembert's paradox. 2 Problem 2. (Total 15 pts) Two viscous, incompressible, unmixable fluids of the same density, , but different viscosities 1 and 2 , flow in separate layers between parallel boundaries located dp at y = h, driven by a constant pressure gradient P = dx < 0. (The lower fluid has viscosity 1 , and the upper has viscosity 2 .) l y
2, y=h y =h x dp/dx =P <0 1, a) (10 pts) Find the velocity field between the boundaries. b) (5 pts) Consider a section of the flow of length l in the xdirection. What is the total force exerted on the fluid by the plates (per unit span in the zdirection) over that length? Compute this in the simplest way possible. 3 4 Problem 3. (Total 15 pts) Consider two infinitely long, concentric circular cylinders, of radii a (the inner cylinder) and b (the outer cylinder). The space between the cylinders is filled with a viscous, incompressible fluid with density and viscosity . Suppose both cylinders rotate with angular velocity (i.e. the linear velocity at the surface of the inner cylinder is a, and at the surface of the outer cylinder is b.) Suppose also that the inner cylinder is moving in the outofpage, zdirection at a constant velocity U . a b The relevant NavierStokes equations are (2.22 in the book): u2 ur ur 2 u 1 p 2 ur  2  2 + (u )ur  =  + t r r r r u u r u 2 ur 1 p u 2 u + 2 + (u )u + = +  2 t r r r r uz 1 p + (u )uz =  + 2 uz t z where we also have (u 2 u + + uz , and r r z 1 2 1 2 = r + 2 2 + 2. r r r r z ) = ur a) (7 pts) Assuming steady flow, find the velocity field between the cylinders. There is no pressure gradient p/z and you can ignore gravity. Assume that the velocity is of the form u = u (r)^ + uz (r)^z . The book example on page 44 is useful. e e b) (8 pts) Suppose the pressure at the surface of the inner cylinder is pa . Find the pressure field between the cylinders. 5 (More space for Problem 3) 6 Bonus questions from the world of the arts.
1. Which Mahler symphony uses sleigh bells? a) 2nd b) 4th c) 9th d) 256th 2. Who of the following is not an Impressionist (could be more than one)? a) Vincent van Gogh b) Jackson Pollock c) Claude Monet d) Anselm Kiefer 3. The Beach Boys' hit records from the 1960s were originally in mono, not stereo, because: a) Brian Wilson was deaf in one ear b) Stereo recording equipment was expensive in those days c) Capitol Records thought stereo was a fad d) They were originally in stereo, what are you talking about? 4. The first line of MobyDick is: a) "It was a dark and stormy night." b) "It was the best of times, it was the worst of times." c) " `Thar she blows!' " d) "Call me Ishmael." 5. My favorite color is: a) Black b) Gray c) Blue d) None, I'm colorblind 7 ME 563  Intermediate Fluid Dynamics  Su Exam 2  15 April 2002 50 points Problem 1  Short answers. (Total 20 pts)
(a) (10 pts) For a Newtonian fluid, the stress tensor is given by Tij = p ij + uj ui + xj xi , Name Solutions and the momentum equation (in index notation, and ignoring gravity) is ui ui 1 p 2 ui = + , + uj t xj xi xj xj and the incompressibility condition is uj = 0. xj Now suppose we have a fluid (nonNewtonian) with stress tensor given by Tij = p ij + ui uj + xj xi + ui uj . ( is a constant.) Write the momentum equation (ignoring gravity) and incompressibility condition in index notation for this fluid. 1 (b) (10 pts) In working the problem of the steady, twodimensional flat plate boundary layer, we said that the (inviscid) main flow velocity, U (x), was constant, so that the pressure gradient was zero, i.e. dp/dx = 0. Then, we found the velocity in the boundary layer subject to the noslip condition. But consider the following points: For the main flow velocity to be constant despite no driving pressure gradient, it seems that there must be no friction losses in the system. But there is friction, because the shear stress at the surface of the flat plate is nonzero. It appears that these two points are inconsistent with each other. Are they? Why or why not? 2 Problem 2  Boundary layers. (30 pts)
a) (5 pts) We've seen how to find velocity fields by minimizing the kinetic energy dissipation. Why can't we use this approach for boundary layers? b) (3 pts) (The rest of this problem is based on Acheson's Exercise 8.5). Consider the flow downstream of a twodimensional streamlined body at high Reynolds number. The flow upstream of the body is constant and in the xdirection, i.e. u = (U, 0); assume that the body is centered at y = 0 and is symmetric about y = 0. There is a thin wake downstream of the body in which the velocity varies much more rapidly with y than with x. Assume that we are far enough downstream of the body that we can write the velocity field as u = (U + u1 (x, y), v1 (x, y)), where u1 < 0, u1 /U  = O( ), and 1. First, draw what the u vs. y profiles might look like for a couple of different positions downstream of the body. 3 c) (7 pts) Let be the characteristic width of the wake (in the ydirection) and L be a characteristic length (in the xdirection). Show by explicitly estimating the magnitudes of the terms in the boundary layer equations that the momentum equation reduces to U Can you relate and L? u1 2 u1 . = x y 2 d) (5 pts) Show that u1 dy = const.
 (1) 4 e) (5 pts) Using control volume arguments, show how (1) relates to the drag on the body. f) (5 pts) Assume a similarity solution for the velocity of the form u1 = F (x)f (), where = For this assumption to be valid, how must g vary with x? y . g(x) 5 Bonus questions from the world of sports.
1. Which three drivers were the first to qualify at over 200 mph for the Indy 500? a) Johnny Rutherford, Al Unser Sr., Rick Mears b) Tom Sneva, Al Unser Sr., Johnny Rutherford c) Bobby Unser, A.J. Foyt, Gordon Johncock d) Rick Mears, Tom Sneva, Danny Ongais 2. Which of these alpine skiing stars represented Luxembourg? a) Pirmin Zurbriggen b) Hanni Wenzel c) Franck Piccard d) Marc Girardelli 3. The infield fly rule pertains in which of these situations (could be more than one)? a) Bases loaded, two out b) Runners on first and third, no outs c) Runners on first and second, one out d) Runners on first and second, two out 4. Wayne Gretzky's singleseason NHL record for goals is: a) 83 b) 86 c) 88 d) 92 5. Which style of play was Martina Navratilova famous for? a) Baseline b) Serve and volley c) Selective attacking d) Mash the ball senseless, then hit the weights 6 ME 563  Intermediate Fluid Dynamics  Su Final Exam  13 May 2002 120 points Problem 1  Short answers. (Total 50 pts) Name (a) (25 pts) In class we saw a demonstration of vortex rings, generated when fluid was forced out of a box through a circular orifice. Vortex rings are able to propagate a long distance through the mechanism illustrated in the figure. We can approximate the vortex ring as a circular vortex tube, and we know that in the inviscid case, vortex lines (and thus tubes) move at the velocity of the fluid. At any part of the ring, the rotation induces the portion of the ring on the opposite side of the circle to move in a direction perpendicular to the plane of the ring. In this way, the ring can propagate indefinitely without change in size or shape, in the absence of viscosity. When two identical vortex rings are produced successively, the second vortex ring catches up with the first. Ignoring viscous effects, explain how this happens. You can do this without equations, but you will need to use Helmholtz's vortex theorems. 1 (Space for part (a).) 2 (b) (10 pts) Consider the transport equation for a conserved scalar quantity, : +u t
inertial diffusive 2 D = 0. Let U be a characteristic velocity, L be a characteristic length scale, and Z be a characteristic scalar concentration. Derive the nondimensional expression that describes the relative importance of the inertial and diffusive terms for the scalar transport equation (analogous to the Reynolds number for the NavierStokes equations). 3 (c) (15 pts) Consider an axisymmetric, turbulent jet flowing in the xdirection. The mean velocity in the xdirection is given by u(x, r) = uc (x)f (), where = r , (x) and, as we saw in class, = Cx and uc = D/x, where C and D are constants. Now define m(x) as the mean mass flux across a surface of fixed x. How does m depend on x? Is the mass flux a constant? Whether it is or isn't, explain what's happening. 4 Problem 2. (Total 35 pts) Consider a viscous, incompressible fluid flowing with uniform velocity U in the zdirection. The fluid encounters a thinwalled, circular tube of radius R = 1 whose axis is aligned with the zaxis. We are interested in the steady flow that develops inside this tube. Working in cylindrical coordinates, the appropriate NavierStokes equations are (u (u u2 ur 2 u 1 p 2 ur  2  2 = + r r r r u r u 2 ur 1 p u 2 )u + u + 2 = +  2 r r r r 1 p (u )uz =  + 2 uz z 1 1 u uz u= (rur ) + + =0 r r r z )ur  (1) (2) (3) (4) making use of the following: (u 2 u + + uz , and r r z 1 2 1 2 = r + 2 2 + 2. r r r r z ) = ur a) (8 pts) Since the boundary conditions are axisymmetric, it is reasonable to assume that the velocity field inside the tube will also be axisymmetric. Rewrite Eqs. (1), (2), (3) and (4) for an axisymmetric velocity field. 5 b) (15 pts) Immediately after the fluid enters the tube, a boundary layer will begin to develop on the tube's interior walls. In the boundary layer, let U be a characteristic value of uz , and let uz change by an amount of order U over a distance L in the zdirection. Also, let be a characteristic thickness of the boundary layer, where L. Working from the equations you developed in part a), derive the boundary layer equations for this flow. 6 c) (12 pts) If we go far enough downstream in the tube, we expect that the flow inside will become parallel in the zdirection, and that uz will have a parabolic profile. Estimate just how far we have to go in the tube for this socalled `fullydeveloped flow' to pertain. Explain your reasoning. 7 Problem 3. (Total 35 pts) This problem concerns the steady, parallel flow in the xdirection of a viscous, incompressible fluid between plane boundaries at y = 0 and y = h. Assume that the flow is fully developed, so none of the flow properties depends on x.
a) (5 pts) Consider the internal energy of the fluid. Is there heat generated by the flow? If so, what happens to this heat, and if not, why not? b) (5 pts) Suppose that the thermal conductivity of the fluid, k, is constant; also assume that the internal energy can be written E = cT , where c is the (constant) specific heat of the fluid and T is the temperature. Write the internal energy equation for the particular flow in this problem. 8 c) (15 pts) Show that there is no solution to the internal energy equation (from part b) for this problem in the case where the boundaries are thermally insulated, i.e. when T /y = 0 at y = 0, h. 9 d) (10 pts) Suppose now that the boundary at y = 0 is stationary, the boundary at y = h moves at a speed U in the xdirection, and instead of being thermally insulated, both boundaries are held at a temperature T0 . What's the temperature field, T (y), between the boundaries? 10 Bonus questions from film and television.
1. Before becoming a private detective, Pierce Brosnan's title character on "Remington Steele" was a) A secret agent b) A police officer c) A classical musician d) An art thief 2. Which school is Ashley Judd's alma mater? a) University of Tennessee b) Vanderbilt University c) Caltech d) University of Kentucky She's a wellknown UK basketball fan 3. Who first played Jack Ryan, the recurring hero in Tom Clancy's novels, in the movies? a) Harrison Ford b) Wesley Snipes c) Alec Baldwin In `The Hunt for Red October' d) Kurt Russell 4. In the movie "Random Hearts", Kristin ScottThomas plays a) A Senator b) A member of the House of Representatives c) A Cabinet secretary d) A casino blackjack dealer 5. Which one of these actresses did not star on "The Facts of Life"? a) Lisa Whelchel b) Molly Ringwald appeared in the first season only c) Dana Plato d) Nancy McKeon 11 530.621: Graduate Fluid Mechanics I
Midterm 1, Fall 1997 Time: 2 hrs. Open books & notes Please write and sketch in legible fashion. 1a. (5 points) The enstrophy per unit volume of a flow is defined as w2 /2 = wiwi /2, where wi is the . vorticity vector w =  x u. Consider a divergencefree velocity field,  u = 0. Using index notation, prove the following identity: w2 = 2 (u2 /2)  u 2 u  () : (u u), where () : (u u) = i j ( ui uj ), and i = /xi. 1b. (5 points) Consider an incompressible, Newtonian, viscous fluid but with spatially varying viscosity m(x). Prove, using index notation, the following statement: "The total viscous force fv obtained by integrating the differential force (generated by the viscous stress) over a closed surface S bounding a control volume V, is given by: . fv = (m 2 u ) dV +  m ( u +  uT) dV. " 2. (20 points) A sprinkler system with two vertical exit pipes is externally driven. The pressure at the inlet is p1 (known), and the angular velocity of the system is constant, equal to W (also known). All pipe diameters are constant, equal to D. Lengths of the two segments are L and 2 L, as shown in figure below. Incompressible fluid, with density r. W A patm 2 Dh patm 3 g
. p1 L 2L Distances between tank entrance, axis of rotation and gate, and the vertical heights at 2 and 3, are all very small compared to L. (cont'd on next page). (a) Initially, at t<0, the gate valve at A is closed and there is no flow. Determine the height difference between fluid levels at exits 2 and 3. (b) At t=0, the gate is suddenly opened completely. For any t>0, formulate appropriate equations that will allow you to obtain the velocity u2 and u3 at exits 2 and 3. Assume uniform, inviscid flow inside the pipes, and neglect the length of the vertical portions as compared to the length of the horizontal portions. Also, for this and next parts (b & c) neglect gravity. (c) Find expressions from which to determine the horizontal (parallel to the pipe) and vertical reaction forces at any time t>0, acting at the weld near A. Neglect moments of inertia of the pipe material, but not of the fluid. 3. (10 points) Consider a porous pipe of constant crosssectional area A and perimeter P. The flow rate into the pipe across the pipe wall is a constant b (per unit area of pipe surface, units of (m3 /s)/m2 ). It does not depend on x. The flow enters the pipe normal to the streamwise direction. Assume the wall shear stress is a known function of x, namely t(x). We wish to find the streamwise velocity distribution u(x) and pressure gradient p/x(x). The flow is not inviscid, is steady, has constant density r, and inside the pipe the flow can be assumed to be unidirectional in the streamwise xdirection. Use a thin control volume as shown in the figure (entire crosssection by dx) to: (a) (5 points) Derive a differential equation relating u(x) to x, A,P,b. Solve, using u(x=0)=U. (b) (5 points) Derive another differential equation relating the pressure gradient to u(x), the known t(x), and the other parameters of the problem. (You need not solve this differential equation). b etc... P A u(x) p(x) t(x) A dx x Computer Assignment 2 StepbyStep Instructions
1 Introduction
As in the first project, GAMBIT will be used to generate the mesh and FLUENT will be used to solve the corner flow. You may use GAMBIT & FLUENT on either the Windows PCs in the Kreiger computer lab or on apserv1.hcf.jhu.edu (Unix, for remote use). Instructions for both systems are given here, although unless you are comfortable with Unix, it is recommended that you follow the Windows instructions. 1.1(a) Starting GAMBIT in Windows To start GAMBIT:
Start > Programs Engineering Applications > Fluent & Gambit > Gambit 2.0 On the main menu bar of GAMBIT:
File > New... Enter the id of your project (these instructions assume you enter corner) and title. Click Accept then Yes . All files will be created and saved in the directory C:/SaveFilesHere. The files can be moved using Windows Explorer later. To open Windows Explorer in Kreiger lab:
Start > Accessories & Media > Command Prompt Type explorer at the prompt. 1.1(b) Starting GAMBIT in Unix From your local system, open (secure) a connection to apserv1 (you must already have an account). Make sure that your system will display graphics windows remotely requested by apserv1. If you use ssh to login to apserv1, the graphics will probably be automatically forwarded to your local machine. At the shell prompt, type (nozzle will be the id of the project):
gambit nozzle & The & symbol runs GAMBIT in the background and gives you back control of your interactive shell. 1.2 GUI Layout The Graphical User Interface layout for GAMBIT is shown below. Sometimes the windows in this field are too long and the bottom part might be hidden under the taskbar at the bottom of the screen (or in Unix, it may be drawn off the bottom edge of the screen)this will hide important buttons such as Apply and Create . You can move the window up by clicking on the window name bar (in dark blue on Windows PCs in Kreiger) and dragging upwards.
Main Menu Bar Operation Toolpad Graphics Window Global Control Toolpad Transcript Window Description Window 1.3 Help/Undo While using GAMBIT, any steps can be undone with the undo button found in the Global control toolpad. The delete button is available most of the time under most options in the operation "Mesh command". Use this as an alternative to the UNDO button to remove unwanted objects. 2 Grid Generation Using GAMBIT
The first step in solving this fluid flow problem is to generate a grid that approximates the true shape of the corner. We begin by setting up a coordinate system, and then defining important vertices for the corner shape. These vertices are then used to define edges, which in turn are used to define surfaces. 2.1 Domain Borders Begin by setting up the coordinate system and defining the borders of the flow domain. Select the TOOLS COMMAND BUTTON Select the COORDINATE SYSTEM COMMAND BUTTON. Select the DISPLAY GRID BUTTON. Under Axis, make sure the diamond next to X is filled in. Set the number in the Minimum box to 0. Set the number in the Maximum box to 10. Update the list. (The name of this button might not show completely, it is directly to the right of the Maximum field) Under Axis, select the diamond next to Y. Set the number in the Minimum box to 0. Set the number in the Maximum box to 10. Update the list. Apply the changes (may need to move Display Grid window up to see the button). Close the Display Grid window Near the bottom right corner of the screen, in the Global Control Toolpad, click FIT TO WINDOW. The coordinate system is complete. Now your graphics window should look something like this: 2.2 Flow Geometry Now we specify the corner geometry. Select the GEOMETRY COMMAND BUTTON. Select the VERTEX COMMAND BUTTON. The Create Real Vertex window should appear. Under Global, make sure the x,y, and z entries are all set to 0. Then press This defines the point (0,0,0) as a vertex of the nozzle and a + should appear at the point showing that it is a vertex. Now we define other vertices: Under Global, set (x,y,z)=(10, 0, 0) and press Repeat with (x,y,z)=(10, 10, 0) and press Repeat with (x,y,z)=(6.67, 10, 0) and press Repeat with (x,y,z)=(6.67, 5 ,0) and press Repeat with (x,y,z)=(0, 5, 0) and press 2.3 Edge Generation Select the EDGE COMMAND BUTTON. The Create Straight Edge window should appear. This is for defining edges using the vertices that were already created. Holding the shift key down, click point (0, 0) in the graphics window and, while keeping the shift key down, click the point (10 , 0) and then press The line between the two points you just clicked should become highlighted. This process of holding shift and clicking on a point will be referred to as shiftclick from now on. Shiftclick points (10, 0) and (10, 10) and press Shiftclick points (10, 10) and (6.67, 10) and press Shiftclick points (6.67, 10) and (6.67, 5) and press Shiftclick points (6.67, 5) and (0, 5) and press Shiftclick points (0, 5) and (0, 0) and press The Edge generation is complete. 2.4 Surface Definition Now, the edges must be combined to define a surface (that will become the corner interior): In the upper right, click the FACE COMMAND BUTTON. Click the arrow next to the Edges field, and move all to the picked column (use All>) and close the edge selection window Press The edge outlines should become highlighted and you should be able to see something like this in the graphics window: 2.5 Mesh Creation Now that we have the corner interior as a surface, we must "mesh" it. Meshing defines the discrete points at which our flow variables will be defined. We start by meshing the edges of the surface: Near the top right, click the MESH COMMAND BUTTON (yellow color). Click the EDGE COMMAND BUTTON. Click the arrow to the right of the Edges field, move edge.1 to the picked column and close the Edge List window Under Spacing, change Interval size to Interval count, enter Nx given to you in the problem statement sheet, and press Repeat, but with edge.2 and an interval count of Ny Repeat, but with edge.3 and an interval count of Nx/3 Repeat, but with edge.4 and an interval count of Ny/2 Repeat, but with edge.5 and an interval count of 2Nx/3 Repeat, but with edge.6 and an interval count of Ny/2 Now you should see something like this: Don't forget for each edge to press Now click the FACE COMMAND BUTTON. The Mesh Faces window will appear. Under Faces, pick face.1 then close the face selection window Apply the mesh face We have now constructed a computational mesh for our domain. You should see something like this: 2.6 Edge Specification We must now specify the type of each edge (inlet, wall, outlet, etc.) From the main menu bar click
Solver > FLUENT 5/6 Now click the ZONES COMMAND BUTTON (blue color) Click the SPECIFY BOUNDARY TYPES COMMAND BUTTON (white color). Next to Entity:, change Faces to Edges. Click on the arrow next to the Edges button and pick edge.6. Change Type: to VELOCITY_INLET and press Apply. If the VELOCITY_INLET option is hidden under the taskbar, select SYMMETRY instead but do not Apply, and then click Type again, the menu will appear higher than before and you can select VELOCITY_INLET. Repeat, but with edges 1 and 2 (move them both to the picked column) and with type WALL and press Repeat, but with edge.3 and type OUTFLOW and press Repeat, but with edges 4 and 5 and type WALL and press The mesh is now ready. 2.7 Mesh Export You now have to save the mesh and export it: From the main menu bar click
File > Save From the main menu bar click
File > Export > Mesh... IMPORTANT: make sure to press the button Export 2d Mesh. The default name of the mesh file will be corner.msh, you can keep it and then press Accept . From the main menu bar click
File > Exit If you do not exit GAMBIT properly you might have problems restarting the program for a short time (due to licensing issues). 3 Corner Flow Computation with FLUENT Now that the grid has been generated, we use FLUENT to specify the velocity at the inlet and then calculate flow on this grid. The flow is assumed to be fully developed channel flow at the inlet, with a specified flow rate Q. To define the inlet profile, a short program written in C++, must be used. 3.1 Inlet Profile Program Listing The program listing is given below, and the first step is to copy this into a file (call it something like "profileu.cpp"). It is NOT necessary to compile the program, since FLUENT can interpret this file. NOTE: YOU WILL HAVE TO CHANGE THE FLOW RATE Q IN THE FILE TO THE ONE GIVEN TO YOU IN THE ASSIGNMENT SHEET. If you have problems email me at yili@pegasus.me.jhu.edu and I can email you the file. Make sure the file is in your working directory.
#include "udf.h" DEFINE_PROFILE(inlet_x_velocity, /* function name */ thread, /* thread */ nv) /* variable number */ { real Q; /* flow rate (m^3/s) */ real height; /* height of channel inlet (m) */ real depth; /* depth of channel (m) */ face_t f; real x[ND_ND]; Q=1.e5; /* enter your flow rate here */ height=5.e5; depth=0.01; /* loop over each of the faces of this zone */ begin_f_loop (f,thread) { F_CENTROID(x,f,thread); /* make sure you get the ()'s right here!! */ F_PROFILE(f,thread,nv) = 6.*Q/(pow(height,3)*depth)*x[1]*(heightx[1]); } end_f_loop (f,thread) } 3.2(a) Start FLUENT in Windows To start FLUENT:
Start > Programs > Engineering Applications > Fluent & Gambit > FLUENT6.0 The FLUENT Version should appear. Select 2d and press Run . 3.2(b) Start FLUENT in Unix At the shell prompt, type:
fluent & A new window should appear. Enter 2d next to version> 3.3 Flow Computation On the menu bar click:
File > Read > Case... This will open the Case submenu, select corner.msh (or whatever you called your mesh) under Files and press OK . On the menu bar click:
Grid > Scale... This opens the Scale submenu. Change the scale factors to
x: y: 1e5 1e5 then press Scale (make sure you only press this button once). This corresponds to a scale of 100 m, which is the dimension of our micromachine. Now Close . We must now define the material; we are using air at 20C which is the default material. On the menu bar click:
Define > Materials... Set the properties to 20C, i.e. Density = 1.205 kg/m3 and Viscosity = 1.81e5 kg/ms Close the Materials window Close . On the menu bar click:
Define > UserDefined > Functions > Interpreted... The Interpreted UDFs window should appear. Type in "profileu.cpp" (or whatever you called the .cpp file) under Source File Name. Press Compile In the main FLUENT window you should see a long line starting with cpp and ending profileu.cpp, make sure the line does not have any error messages. An error message indicates that there is probably a typo in the .cpp file. Close the Interpreted UDFs window Close On the menu bar click: Define > Boundary Conditions Select the velocity_inlet entry then press the Set... Make sure the Velocity Specification Method is Magnitude, Normal to Boundary. Leave the value in the Velocity Magnitude window at 0. In the field to the right of Velocity Magnitude change it from Constant to udf inlet_x_velocity. Press OK Close the Boundary Conditions window Close . On the menu bar click: Solve > Initialize > Initialize... Press Init
Close On the menu bar click: Solve > Iterate... Set the number of iterations to 20 and press Iterate This will start the solver and execute 20 iterations. Some results will be displayed in the main FLUENT window. Keep iterating in steps of 20 (just repress the Iterate button) until the main FLUENT window displays the message solution is converged. Close the iterate window. Now the computational part is done. 3.3 PostProcessing Postprocessing must be performed to extract the desired data. First we have to change the function of the mouse buttons. On the menu bar click: Display > Mouse buttons... Change the function of the right mouse button to mousezoom. Click OK 3.3.1 Grid Display Now we plot the grid On the menu bar click: Display > Grid... Then press Display . A window containing the grid should appear. Zoom in by windowing the grid while holding down the right mouse button. Start the window at the upper left and end at lower right. If you need to zoom out, hold the right mouse button and draw a window going from lower left to upper right. 3.3.1.1 Save the Plot On the main FLUENT window menu bar click: File > Hardcopy ... Select the JPEG format (or PostScript, etc.) for the plot. Select Gray scale (if you want to print on a grayscale printer) or color for the coloring option. Leave all other options as is and click Save... . Enter the name you want to call your plot and press OK . Make sure you ad the extension jpg to your filename. If you don't you can edit the filename later in Windows Explorer or just treat the file as a jpg. The figure will have the same size as the plot on the screen, you can modify the size of your plot window before saving the hardcopy. 3.3.2 Velocity Vectors On the menu bar click: Display > Vectors... The Velocity... and Velocity Magnitude options will be selected by default; keep them. You might need to change the scale to get a nice looking plot, start by testing a Scale of 1. Press Display Once you get a nice plot, save the plot as you did in section 3.3.1.1. The velocity vector plot should look like this: 3.3.3 Streamlines While keeping the velocity vector plot open: On the menu bar click: Display > Contours... Change the Contours of to the Velocity... option. Change the second box to Stream Function option. Set Levels to 30. Click Display . Save the plot as you did in section 3.3.1.1. The plot should look like: 3.3.4 XVelocity Profiles To output xvelocity profiles, we first define a few stations (xlocation) for the profiles: On the menu bar click: Surface > Isosurface... In the Surface of Constant box, select Grid... in the upper box. In the second box select XCoordinate. In the IsoValues box, enter the x coordinate where you want a velocity profile and then press Create . For this problem it would be good to create stations at: X=2e5; X=4e5; X=6e5; X=8e5. And with YCoordinate selected in the second box: Y=2e5; Y=4e5; Y=6e5; Y=8e5. Note the names of the lines that are being created in the New Surface Name box. Press Create each time you change the value in the IsoValues box.
Close the Isosurface window. We now want to display the velocity profiles along the lines we just defined. On the menu bar click: Display > Contours... In Contours of select Velocity... In the second box select X Velocity. Check the box next to Draw Profiles . The Profile Options window will pop up. Select a scale factor of 1e6 (You may need to play with the scale factor to get a nice plot, depending on Q and the number of vertical lines you created). After closing the draw profiles window, it can be opened again by deselecting and reselecting the Draw Profiles option. Press Apply. And then Close the profile option window. Under surfaces, select velocity inlet and ALL of your vertical lines (with names like xcoordinate5, etc.). Press Display . Save the plot as you did in section 3.3.1.1. To generate the profiles of the yvelocity, repeat, but: (i) Change XVelocity to YVelocity. (ii) Change selected surfaces to: OUTFLOW, and ALL of your horizontal lines (with names like ycoordinate5, etc). (iii) In the Profile Options window where you select the Scale, change the Scale to 5e7 and change projection direction to (0,1,0) (previously it was (1,0,0)). Press Display . The plots should look like this: xvelocity profiles yvelocity profiles 3.3.5 Forces on Wall One way to find the force on the wall is to begin by extracting the velocity near the wall, and use this to estimate the derivative of the velocity profile at the wall. This information can then be used to calculate the wall shear stress and hence the total force on the wall. FLUENT does have the capability to calculate derivatives and forces automatically, and you may want to use these as a rough check of your results. The first steps are to generate isosurfaces of constant xcoordinate and constant ycoordinate near the wall. We will extract the velocity along these lines. Follow the steps for creating these isosurfaces given above, and create the constant ycoordinate line at y = 50 (Dy/2) mm, where Dy is the vertical spacing of your mesh. Likewise, create the constant xcoordinate line at x = 66.67+(Dx/2) mm, where Dx is the horizontal spacing of your mesh. Note that Dx = Lx/(Nx1), and similarly for Dy. 3.3.5.1 Horiztonal Line We'll start with the line along the horizontal part of the wall: On the menu bar click: Plot > XY Plot... The Solution XY Plot window will popup. In the upperright corner under Y Axis Function select Velocity... in the first box, and Xvelocity in the second box. Under Surfaces, select only the line you created at y = 50 (Dy/2) mm and press Plot . Make sure the values look reasonable and make a note of the units (since these won't be in your text file). Note that in calculating the derivatives and forces, you won't be using the points past where the wall turns, i.e. in this case we won't be using x > 66.7 m. Keep the plot open. We will now save the data in a text file. 3.3.5.2 Write to File In the Solution XY Plot window check the box next to the option write to file . Click the Write... button that will appear at the bottom of the window. Call the file something like horizontalline.xy and press OK . You can view this file with any text editor or import it into a spreadsheet program (e.g. Excel). 3.3.5.3 Vertical Line Now for the line near the vertical wall. On the menu bar click: Plot > XY Plot... The Solution XY Plot window will popup. Deselect the writing option: p write to file Under Plot Direction change (X,Y)=(1,0) to (0,1). In the upperright corner under Y Axis Function select Velocity... then Yvelocity. Under Surfaces, select only the line you created at x = 66.67+(Dx/2) mm and press Plot . Make sure the values look reasonable and make a note of the units (since these won't be in your text file). Keep in mind you'll only use the values that are actually near a wall (as above). 3.4 Finalizing Now we are done, we just have to save our work: On the menu bar click: File > Write > Case & Data... The default name corner.cas should be fine. Press OK On the menu bar click: File > Exit 3.5 Important Notes It would be best if you take all of the files you created in C:/SaveFilesHere (on Windows only), and either save them on a floppy diskette or ftp them somewhere else and then erase them from the PC. At least move them into a new folder that you create under C:/SaveFilesHere. If you have money on your JHUID card, you can print in BW or color in the Krieger lab. BW printouts are 5c/page, color printouts are 25c/page. Good luck! ME 563  Intermediate Fluid Dynamics  Su
Some review materials for the first exam 3/1/02
Note: Review session Wednesday 2/27, 5:307 pm, 318 ME (normal class room). 1 Framework of ideas With all of the equations and theorems we've covered, it is important that we be able to organize the various ideas into a clear framework, so that we know what assumptions we need to make to apply which results, and so we know under what circumstances we can apply which assumptions. The property that is common to all flows we have discussed in this class is that they are incompressible. The major distinction in the class so far is between inviscid and viscous fluids. To recap: Ideal Fluids Euler's equations The only force acting across a surface element in the fluid is the pressure Flows are driven by pressure gradients (p 0), or gravity Incompressible: u = 0 Kelvin circulation theorem, Helmholtz vortex theorems, CauchyLagrange theorem Steady flow: Bernoulli streamline theorem  If also irrotational: Bernoulli theorem for entire flow 2D flow: vorticity conserved for fluid elements  If also steady: vorticity constant on streamlines Besides pressure, there are also viscous forces acting across surface elements in the fluid Flows observe the noslip condition at solid boundaries Viscous, incompressible fluids NavierStokes equations Viscous forces: shear forces (resist `sliding'), normal forces (resist compression/expansion) Viscous forces = additional flow drivers Because of noslip condition, viscous effects most strongly felt near walls viscousaffected region = boundary layer Viscosity diffuses momentum and vorticity  analogous to diffusion of heat 1 As we have discussed, we are generally free to use ideal fluid concepts until boundary layer effects become significant. Obviously, this can happen if there is a solid boundary in the problem. In the wing example, Acheson points out that ideal fluid theory is pretty good at predicting the lift for small angles of attack, , but that ideal fluid theory fails when the angle of attack gets large. Even though real fluids (e.g. air) are viscous, at small ideal fluid theory works OK because the boundary layer, where viscous effects are important, is thin and occupies a very small portion of the flow. At higher , the wing experiences boundary layer separation, so the boundary layer occupies a much larger portion of the flow and viscous effects can no longer be neglected. Another instructive example is the one in the book about the impulsively moved plane boundary (2.3). In that case, the boundary and fluid above it are initially at rest. Suddenly, the boundary starts moving at a fixed speed, which it then maintains for all time. As time passes, fluid that lies further and further away from the boundary gets accelerated, as viscosity allows the momentum (and vorticity) which is introduced to the fluid at the boundary (by the noslip condition) to get diffused upward into the regions of the fluid that are essentially still at rest. 2 Similarity solutions The impulsively moved plane boundary problem also introduced the idea of similarity solutions. In their most basic form, similarity solutions allow us to reduce the number of independent variables in a problem from two to one, by introducing a similarity variable that relates the original two independent variables. Then, the original partial differential equation, which involves two variables, can now be written as an ordinary differential equation, which involves just one variable. The value of the similarity solution approach arises because ordinary differential equations are much easier to solve, analytically or numerically. Let's review what we did for the impulsively moved plane boundary, emphasizing the similarity aspects. Our original equation was u 2u = 2, t y (1) which is a partial differential equation with independent variables t and y. In the notes, we argued for a similarity solution on physical grounds. The mathematical argument goes as follows. The variable y enters (1) as a secondorder term, so if we define a new variable y = y, where is some constant, then (1) becomes u 2u . = 2 t u 2 (2) The variable t enters as a firstorder term, so if we define a new variable t = t, where is some constant, (2) becomes u 2u = 2 . t u 2 (3) By inspection, if = 2 , then (3) is the same equation as (1), except that the independent variables are t and y . This means that if we do the transformation y = y t = 2 t (4) the mathematics of the problem defined by 1 is the same. We can also interpret the transformation (4) as meaning that the problem is unchanged if y/ t is unchanged, i.e. if y / t = y/ t. This then 2 suggests that the solution of the equation (1) can be expressed in terms of the variable = y/ t. We are free to multiply this resulting variable by any constant we like (except zero), since that won't affect the dependence of the problem on y and t. We can choose, as the book does, to let = y/ t, because this removes the constant from (1). This is also convenient because defined this way is dimensionless. With this new independent variable, we then look for a solution of the form u = U f () or, as the book does, u = f (). There is no difference in these two approaches, since U is a constant in the problem. With either form for u, (1) becomes f + f = 0, 2 (5) which is an ordinary differential equation involving the single independent variable, . As we saw in class, there's no problem solving (5), subject to the proper boundary conditions. This process may seem kind of arbitrary, since it seems like we're guessing the form of the solution. What's to stop us from guessing some other solution? To which the response is twofold. First, if you guess a solution and it's right, then you're done, because you can't have different solutions to the same problem if the problem is fully specified1 . Secondly, let's look what happens if we guess a different similarity solution for (1). Say, let's look for a solution of the form u = f () Then we have u y df = = f () 2 t d t t u 1 df = = f () y d y t 2u 1 = f () 2 . y 2 t Inserting these into (1), we get f y 1 = f 2 , which gives 2 t t y f + f = 0, where = y . t (6) (7) which isn't an ordinary differential equation in a single independent variable; in fact it's a lot more complicated, because we have a free y floating around. This means that the solution f to (7) will depend on y, which is inconsistent with (6), where we assumed that f was only a function of . Thus (6) is not a valid solution of the differential equation (1). So you see that what we did in class to solve the problem wasn't quite so arbitrary after all, because if you guess a solution and it isn't valid, then the equations will tell you. Things aren't quite this simple really there are some complications with higherorder differential equations, etc. but for us this statement is valid. 1 3 Turbines  water  power generation
Impulse Turbines  Pelton wheel (kinetic energy only  nozzle  at patm) Reaction Turbines (both kinetic energy and delta P  energy conversion takes place in enclosed space above patm) Francis Kaplan Pelton wheel Francistype reaction turbine Kaplantype reaction turbine Jet engine:  gas turbine General Electric CF650 Gas turbine
ANALYSIS: Velocity Triangles Thermal aspects: blade cooling Pumps: Centifugal Axial Pump Compressors AXIAL General Electric CF650 Complex Flow in Compressors Fans
This 4stage ATEGG Compressor Rig, the highest loaded compressor ever built with acceptable stall margin and efficiency characteristics, meets the Phase III stage loading goals. JTAGG III Advanced Concept Centrifugal Impeller with independent inducer and exducer will provide higher pressure ratio and efficiency. JTDE Forward Swept Fan Blisk will be the first large forward swept rotor to be tested in a demonstrator engine. Integrating Forward Sweep and Splitter Technology are key features to achieving high efficiency and high pressure ratio. This stage will be utilized as the JTAGG III low pressure compressor. A Forged Orthorhombic Transformed Super Alpha2 Billet will be bonded to Gamma TiAl to form a novel, dual alloy impeller. Centrifugal compressors Specific speed: N ...
View
Full
Document
This note was uploaded on 04/18/2008 for the course CE 320 taught by Professor Kf during the Spring '08 term at NJIT.
 Spring '08
 kf

Click to edit the document details