Free Essay

Fluid Mechanics

In: Computers and Technology

Submitted By shegzydee
Words 50231
Pages 201
This is page i Printer: Opaque this

A Mathematical Introduction to Fluid Mechanics
Alexandre Chorin
Department of Mathematics University of California, Berkeley Berkeley, California 94720-3840, USA

Jerrold E. Marsden
Control and Dynamical Systems, 107-81 California Institute of Technology Pasadena, California 91125, USA

ii

iii

A Mathematical Introduction to Fluid Mechanics

iv

Library of Congress Cataloging in Publication Data Chorin, Alexandre A Mathematical Introduction to Fluid Mechanics, Third Edition

(Texts in Applied Mathematics) Bibliography: in frontmatter Includes. 1. Fluid dynamics (Mathematics) 2. Dynamics (Mathematics) I. Marsden, Jerrold E. II. Title. III. Series. ISBN 0-387 97300-1 American Mathematics Society (MOS) Subject Classification (1980): 76-01, 76C05, 76D05, 76N05, 76N15 Copyright 1992 by Springer-Verlag Publishing Company, Inc. All rights reserved. No part of this publication may be reproduced, stored in a retrieval system, or transmitted, in any or by any means, electronic, mechanical, photocopying, recording, or otherwise, without the prior written permission of the publisher, Springer-Verlag Publishing Company, Inc., 175 Fifth Avenue, New York, N.Y. 10010. Typesetting and illustrations prepared by June Meyermann, Gregory Kubota, and Wendy McKay The cover illustration shows a computer simulation of a shock diffraction by a pair of cylinders, by John Bell, Phillip Colella, William Crutchfield, Richard Pember, and Michael Welcome.

The corrected fourth printing, April 2000.

v

Series Preface Page (to be inserted)

vi

blank page

This is page vii Printer: Opaque this

Preface

This book is based on a one-term course in fluid mechanics originally taught in the Department of Mathematics of the University of California, Berkeley, during the spring of 1978. The goal of the course was not to provide an exhaustive account of fluid mechanics, nor to assess the engineering value of various approximation procedures. The goals were: • to present some of the basic ideas of fluid mechanics in a mathematically attractive manner (which does not mean “fully rigorous”); • to present the physical background and motivation for some constructions that have been used in recent mathematical and numerical work on the Navier–Stokes equations and on hyperbolic systems; and • to interest some of the students in this beautiful and difficult subject. This third edition has incorporated a number of updates and revisions, but the spirit and scope of the original book are unaltered. The book is divided into three chapters. The first chapter contains an elementary derivation of the equations; the concept of vorticity is introduced at an early stage. The second chapter contains a discussion of potential flow, vortex motion, and boundary layers. A construction of boundary layers using vortex sheets and random walks is presented. The third chapter contains an analysis of one-dimensional gas flow from a mildly modern point of view. Weak solutions, Riemann problems, Glimm’s scheme, and combustion waves are discussed. The style is informal and no attempt is made to hide the authors’ biases and personal interests. Moreover, references are limited and are by no

viii

Preface

means exhaustive. We list below some general references that have been useful for us and some that contain fairly extensive bibliographies. References relevant to specific points are made directly in the text.
R. Abraham, J. E. Marsden, and T. S. Ratiu [1988] Manifolds, Tensor Analysis and Applications, Springer-Verlag: Applied Mathematical Sciences Series, Volume 75. G. K. Batchelor [1967] An Introduction to Fluid Dynamics, Cambridge Univ. Press. G. Birkhoff [1960] Hydrodynamics, a Study in Logic, Fact and Similitude, Princeton Univ. Press. A. J. Chorin [1976] Lectures on Turbulence Theory, Publish or Perish. A. J. Chorin [1989] Computational Fluid Mechanics, Academic Press, New York. A. J. Chorin [1994] Vorticity and Turbulence, Applied Mathematical Sciences, 103, Springer-Verlag. R. Courant and K. O. Friedrichs [1948] Supersonic Flow and Shock Waves, WileyInterscience. P. Garabedian [1960] Partial Differential Equations, McGraw-Hill, reprinted by Dover. S. Goldstein [1965] Modern Developments in Fluid Mechanics, Dover. K. Gustafson and J. Sethian [1991] Vortex Flows, SIAM. O. A. Ladyzhenskaya [1969] The Mathematical Theory of Viscous Incompressible Flow , Gordon and Breach. L. D. Landau and E. M. Lifshitz [1968] Fluid Mechanics, Pergamon. P. D. Lax [1972] Hyperbolic Systems of Conservation Laws and the Mathematical Theory of Shock Waves, SIAM. A. J. Majda [1986] Compressible Fluid Flow and Systems of Conservation Laws in Several Space Variables, Springer-Verlag: Applied Mathematical Sciences Series 53. J. E. Marsden and T. J. R. Hughes [1994] The Mathematical Foundations of Elasticity, Prentice-Hall, 1983. Reprinted with corrections, Dover, 1994. J. E. Marsden and T. S. Ratiu [1994] Mechanics and Symmetry, Texts in Applied Mathematics, 17, Springer-Verlag. R. E. Meyer [1971] Introduction to Mathematical Fluid Dynamics, Wiley, reprinted by Dover. K. Milne–Thomson [1968] Theoretical Hydrodynamics, Macmillan. C. S. Peskin [1976] Mathematical Aspects of Heart Physiology, New York Univ. Lecture Notes. S. Schlichting [1960] Boundary Layer Theory, McGraw-Hill. L. A. Segel [1977] Mathematics Applied to Continuum Mechanics, Macmillian. J. Serrin [1959] Mathematical Principles of Classical Fluid Mechanics, Handbuch der Physik, VIII/1, Springer-Verlag. R. Temam [1977] Navier–Stokes Equations, North-Holland.

Preface

ix

We thank S. S. Lin and J. Sethian for preparing a preliminary draft of the course notes—a great help in preparing the first edition. We also thank O. Hald and P. Arminjon for a careful proofreading of the first edition and to many other readers for supplying both corrections and support, in particular V. Dannon, H. Johnston, J. Larsen, M. Olufsen, and T. Ratiu and G. Rublein. These corrections, as well as many other additions, some exercises, updates, and revisions of our own have been incorporated into the second and third editions. Special thanks to Marnie McElhiney for typesetting the second edition, to June Meyermann for typesetting the third edition, and to Greg Kubota and Wendy McKay for updating the third edition with corrections. Alexandre J. Chorin
Berkeley, California

Jerrold E. Marsden
Pasadena, California

Summer, 1997

x

Preface

This is page xi Printer: Opaque this

Contents

Preface 1 The 1.1 1.2 1.3 Equations of Motion Euler’s Equations . . . . . . . . . . . . . . . . . . . . . . . . Rotation and Vorticity . . . . . . . . . . . . . . . . . . . . . The Navier–Stokes Equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

vii 1 1 18 31 47 47 67 82 95 101 101 115 135 143 161

2 Potential Flow and Slightly Viscous Flow 2.1 Potential Flow . . . . . . . . . . . . . . . 2.2 Boundary Layers . . . . . . . . . . . . . . 2.3 Vortex Sheets . . . . . . . . . . . . . . . . 2.4 Remarks on Stability and Bifurcation . . 3 Gas 3.1 3.2 3.3 3.4 Index Flow in One Dimension Characteristics . . . . . . . Shocks . . . . . . . . . . . . The Riemann Problem . . . Combustion Waves . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

xii

Contents

This is page 1 Printer: Opaque this

1
The Equations of Motion

In this chapter we develop the basic equations of fluid mechanics. These equations are derived from the conservation laws of mass, momentum, and energy. We begin with the simplest assumptions, leading to Euler’s equations for a perfect fluid. These assumptions are relaxed in the third section to allow for viscous effects that arise from the molecular transport of momentum. Throughout the book we emphasize the intuitive and mathematical aspects of vorticity; this job is begun in the second section of this chapter.

1.1

Euler’s Equations

Let D be a region in two- or three-dimensional space filled with a fluid. Our object is to describe the motion of such a fluid. Let x ∈ D be a point in D and consider the particle of fluid moving through x at time t. Relative to standard Euclidean coordinates in space, we write x = (x, y, z). Imagine a particle (think of a particle of dust suspended) in the fluid; this particle traverses a well-defined trajectory. Let u(x, t) denote the velocity of the particle of fluid that is moving through x at time t. Thus, for each fixed time, u is a vector field on D, as in Figure 1.1.1. We call u the (spatial ) velocity field of the fluid. For each time t, assume that the fluid has a well-defined mass density ρ(x, t). Thus, if W is any subregion of D, the mass of fluid in W at time t

2

1 The Equations of Motion

trajectory of fluid particle

D

x

u(x,t)

Figure 1.1.1. Fluid particles flowing in a region D.

is given by m(W, t) =
W

ρ(x, t) dV,

where dV is the volume element in the plane or in space. In what follows we shall assume that the functions u and ρ (and others to be introduced later) are smooth enough so that the standard operations of calculus may be performed on them. This assumption is open to criticism and indeed we shall come back and analyze it in detail later. The assumption that ρ exists is a continuum assumption. Clearly, it does not hold if the molecular structure of matter is taken into account. For most macroscopic phenomena occurring in nature, it is believed that this assumption is extremely accurate. Our derivation of the equations is based on three basic principles: i ii iii mass is neither created nor destroyed ; the rate of change of momentum of a portion of the fluid equals the force applied to it (Newton’s second law ); energy is neither created nor destroyed.

Let us treat these three principles in turn.

i Conservation of Mass
Let W be a fixed subregion of D (W does not change with time). The rate of change of mass in W is d d m(W, t) = dt dt ρ(x, t) dV =
W W

∂ρ (x, t) dV. ∂t

1.1 Euler’s Equations

3

Let ∂W denote the boundary of W , assumed to be smooth; let n denote the unit outward normal defined at points of ∂W ; and let dA denote the area element on ∂W . The volume flow rate across ∂W per unit area is u · n and the mass flow rate per unit area is ρu · n (see Figure 1.1.2).

u n portion of the boundary of W

Figure 1.1.2. The mass crossing the boundary ∂W per unit time equals the surface integral of ρu · n over ∂W.

The principle of conservation of mass can be more precisely stated as follows: The rate of increase of mass in W equals the rate at which mass is crossing ∂W in the inward direction; i.e., d dt

ρ dV = −
W ∂W

ρu · n dA.

This is the integral form of the law of conservation of mass. By the divergence theorem, this statement is equivalent to ∂ρ + div(ρu) dV = 0. ∂t

W

Because this is to hold for all W , it is equivalent to ∂ρ + div(ρu) = 0. ∂t The last equation is the differential form of the law of conservation of mass, also known as the continuity equation. If ρ and u are not smooth enough to justify the steps that lead to the differential form of the law of conservation of mass, then the integral form is the one to use.

4

1 The Equations of Motion

ii Balance of Momentum
Let x(t) = (x(t), y(t), z(t)) be the path followed by a fluid particle, so that the velocity field is given by u(x(t), y(t), z(t), t) = (x(t), y(t), z(t)), ˙ ˙ ˙ that is, dx (t). dt This and the calculation following explicitly use standard Euclidean coordinates in space (delete z for plane flow).1 u(x(t), t) = The acceleration of a fluid particle is given by a(t) = d2 d x(t) = u(x(t), y(t), z(t), t). dt2 dt

By the chain rule, this becomes a(t) = Using the notation ux = and u(x, y, z, t) = (u(x, y, z, t), v(x, y, z, t), w(x, y, z, t)), we obtain a(t) = uux + vuy + wuz + ut , which we also write as a(t) = ∂t u + u · ∇u, where ∂t u = ∂u ∂t and u · ∇ = u ∂ ∂ ∂ +v +w . ∂x ∂y ∂z ∂u , ∂x ut = ∂u , ∂t etc., ∂u ∂u ∂u ∂u x+ ˙ y+ ˙ z+ ˙ . ∂x ∂y ∂z ∂t

1 Care must be used if other coordinate systems (such as spherical or cylindrical) are employed. Other coordinate systems can be handled in two ways: first, one can proceed more intrinsically by developing intrinsic (i.e., coordinate free) formulas that are valid in any coordinate system, or, second, one can do all the derivations in Euclidean coordinates and transform final results to other coordinate systems at the end by using the chain rule. The second approach is clearly faster, although intellectually less satisfying. See Abraham, Marsden and Ratiu [1988] (listed in the front matter) for information on the former approach. For reasons of economy we shall do most of our calculations in standard Euclidean coordinates.

1.1 Euler’s Equations

5

We call

D = ∂t + u · ∇ Dt the material derivative; it takes into account the fact that the fluid is moving and that the positions of fluid particles change with time. Indeed, if f (x, y, z, t) is any function of position and time (scalar or vector), then by the chain rule, d Df f (x(t), y(t), z(t), t) = ∂t f + u · ∇f = (x(t), y(t), z(t), t). dt Dt

For any continuum, forces acting on a piece of material are of two types. First, there are forces of stress, whereby the piece of material is acted on by forces across its surface by the rest of the continuum. Second, there are external, or body, forces such as gravity or a magnetic field, which exert a force per unit volume on the continuum. The clear isolation of surface forces of stress in a continuum is usually attributed to Cauchy. Later, we shall examine stresses more generally, but for now let us define an ideal fluid as one with the following property: For any motion of the fluid there is a function p(x, t) called the pressure such that if S is a surface in the fluid with a chosen unit normal n, the force of stress exerted across the surface S per unit area at x ∈ S at time t is p(x, t)n; i.e., force across S per unit area = p(x, t)n. Note that the force is in the direction n and that the force acts orthogonally to the surface S; that is, there are no tangential forces (see Figure 1.1.3).

S n force across S = pn

Figure 1.1.3. Pressure forces across a surface S.

Of course, the concept of an ideal fluid as a mathematical definition is not subject to dispute. However, the physical relevance of the notion (or mathematical theorems we deduce from it) must be checked by experiment. As we shall see later, ideal fluids exclude many interesting real physical

6

1 The Equations of Motion

phenomena, but nevertheless form a crucial component of a more complete theory. Intuitively, the absence of tangential forces implies that there is no way for rotation to start in a fluid, nor, if it is there at the beginning, to stop. This idea will be amplified in the next section. However, even here we can detect physical trouble for ideal fluids because of the abundance of rotation in real fluids (near the oars of a rowboat, in tornadoes, etc.). If W is a region in the fluid at a particular instant of time t, the total force exerted on the fluid inside W by means of stress on its boundary is S∂W = {force on W } = −
∂W

pn dA

(negative because n points outward). If e is any fixed vector in space, the divergence theorem gives e · S∂W = −
∂W

pe · n dA = −
W

div(pe) dV = −
W

(grad p) · e dV.

Thus, S∂W = −
W

grad p dV.

If b(x, t) denotes the given body force per unit mass, the total body force is B= ρb dV.
W

Thus, on any piece of fluid material, force per unit volume = −grad p + ρb. By Newton’s second law (force = mass × acceleration) we are led to the differential form of the law of balance of momentum: ρ Du = − grad p + ρb. Dt (BM1)

Next we shall derive an integral form of balance of momentum in two ways. We derive it first as a deduction from the differential form and second from basic principles. From balance of momentum in differential form, we have ρ ∂u = −ρ(u · ∇)u − ∇p + ρb ∂t

and so, using the equation of continuity, ∂ (ρu) = − div(ρu)u − ρ(u · ∇)u − ∇p + ρb. ∂t

1.1 Euler’s Equations

7

If e is any fixed vector in space, one checks that e· ∂ (ρu) = − div(ρu)u · e − ρ(u · ∇)u · e − (∇p) · e + ρb · e ∂t = − div(pe + ρu(u · e)) + ρb · e.

Therefore, if W is a fixed volume in space, the rate of change of momentum in direction e in W is e· d dt ρu dV = −
W ∂W

(pe + ρu(e · u)) · n dA +
W

ρb · e dV

by the divergence theorem. Thus, the integral form of balance of momentum becomes: d dt ρu dV = −
W ∂W

(pn + ρu(u · n)) dA +
W

ρb dV.

(BM2)

The quantity pn+ρu(u·n) is the momentum flux per unit area crossing ∂W , where n is the unit outward normal to ∂W . This derivation of the integral balance law for momentum proceeded via the differential law. With an eye to assuming as little differentiability as possible, it is useful to proceed to the integral law directly and, as with conservation of mass, derive the differential form from it. To do this carefully requires us to introduce some useful notions. As earlier, let D denote the region in which the fluid is moving. Let x ∈ D and let us write ϕ(x, t) for the trajectory followed by the particle that is at point x at time t = 0. We will assume ϕ is smooth enough so the following manipulations are legitimate and for fixed t, ϕ is an invertible mapping. Let ϕt denote the map x → ϕ(x, t); that is, with fixed t, this map advances each fluid particle from its position at time t = 0 to its position at time t. Here, of course, the subscript does not denote differentiation. We call ϕ the fluid flow map. If W is a region in D, then ϕt (W ) = Wt is the volume W moving with the fluid . See Figure 1.1.4. The “primitive” integral form of balance of momentum states that d dt ρu dV = S∂Wt +
Wt Wt

ρb dV,

(BM3)

that is, the rate of change of momentum of a moving piece of fluid equals the total force (surface stresses plus body forces) acting on it. These two forms of balance of momentum (BM1) and (BM3) are equivalent. To prove this, we use the change of variables theorem to write d dt ρu dV =
Wt

d dt

(ρu)(ϕ(x, t), t)J(x, t) dV,
W

8

1 The Equations of Motion

moving fluid

D W t =0 t Wt

Figure 1.1.4. Wt is the image of W as particles of fluid in W flow for time t.

where J(x, t) is the Jacobian determinant of the map ϕt . Because the volume is fixed at its initial position, we may differentiate under the integral sign. Note that ∂ (ρu)(ϕ(x, t), t) = ∂t D ρu (ϕ(x, t), t) Dt

is the material derivative, as was shown earlier. (If you prefer, this equality says that D/Dt is differentiation following the fluid.) Next, we learn how to differentiate J(x, t). Lemma ∂ J(x, t) = J(x, t)[div u(ϕ(x, t), t)]. ∂t

Proof Write the components of ϕ as ξ(x, t), η(x, t), and ζ(x, t). First, observe that ∂ ϕ(x, t) = u(ϕ(x, t), t), ∂t by definition of the velocity field of the fluid. The determinant J can be differentiated by recalling that the determinant of a matrix is multilinear in the columns (or rows). Thus, holding x

1.1 Euler’s Equations

9

fixed throughout, we have  ∂ ∂ξ  ∂t ∂x   ∂  ∂ ∂ξ J =  ∂t ∂y ∂t   ∂ ∂ξ ∂t ∂z ∂ξ  ∂x    ∂ξ +  ∂y   ∂ξ ∂z Now write 

∂η ∂x ∂η ∂y ∂η ∂z ∂η ∂x ∂η ∂y ∂η ∂z

  ∂ζ ∂ξ ∂x   ∂x     ∂ζ   ∂ξ + ∂y   ∂y   ∂ζ   ∂ξ ∂z ∂z  ∂ ∂ζ ∂t ∂x    ∂ ∂ζ  . ∂t ∂y   ∂ ∂ζ  ∂t ∂z

∂ ∂η ∂t ∂x ∂ ∂η ∂t ∂y ∂ ∂η ∂t ∂z

 ∂ζ ∂x    ∂ζ   ∂y   ∂ζ  ∂z

∂ ∂ξ ∂ ∂ξ ∂ = = u(ϕ(x, t), t), ∂t ∂x ∂x ∂t ∂x ∂ ∂ξ ∂ ∂ξ ∂ = = u(ϕ(x, t), t), ∂t ∂y ∂y ∂t ∂y . . . ∂ ∂ζ ∂ ∂ζ ∂ = = w(ϕ(x, t), t). ∂t ∂z ∂z ∂t ∂z The components u, v, and w of u in this expression are functions of x, y, and z through ϕ(x, t); therefore, ∂ ∂u ∂ξ ∂u ∂η ∂u ∂ζ u(ϕ(x, t), t) = + + , ∂x ∂ξ ∂x ∂η ∂x ∂ζ ∂x . . . ∂ ∂w ∂ξ ∂w ∂η ∂w ∂ζ w(ϕ(x, t), t) = + + . ∂z ∂ξ ∂z ∂η ∂z ∂ζ ∂z When these are substituted into the above expression for ∂J/∂t, one gets for the respective terms ∂u ∂v ∂w J+ J+ J = (div u)J. ∂x ∂y ∂z

10

1 The Equations of Motion

From this lemma, we get d dt ρu dV =
Wt W

D ρu (ϕ(x, t), t) + (ρu)(div u)(ϕ(x, t), t) Dt × J(x, t) dV D (ρu) + (ρ div u)u dV, Dt

=
Wt

where the change of variables theorem was again used. By conservation of mass, D ∂ρ ρ + ρ div u = + div(ρu) = 0, Dt ∂t and thus d dt ρu dV =
Wt Wt

ρ

Du dV. Dt

In fact, this argument proves the following theorem. Transport Theorem d dt For any function f of x and t, we have Df dV. Dt

ρf dV =
Wt Wt

ρ

In a similar way, one can derive a form of the transport theorem without a mass density factor included, namely, d dt ∂f + div(f u) dV. ∂t

f dV =
Wt Wt

If W , and hence, Wt , is arbitrary and the integrands are continuous, we have proved that the “primitive” integral form of balance of momentum is equivalent to the differential form (BM1). Hence, all three forms of balance of momentum—(BM1), (BM2), and (BM3)—are mutually equivalent. As an exercise, the reader should derive the two integral forms of balance of momentum directly from each other. The lemma ∂J/∂t = (div u) J is also useful in understanding incompressibility. In terms of the notation introduced earlier, we call a flow incompressible if for any fluid subregion W , volume(Wt ) =
Wt

dV = constant in t.

1.1 Euler’s Equations

11

Thus, incompressibility is equivalent to 0= d dt dV =
Wt

d dt

J dV =
W W

(div u)J dV =
Wt

(div u) dV

for all moving regions Wt . Thus, the following are equivalent: (i) the fluid is incompressible; (ii) div u = 0;

(iii) J ≡ 1. From the equation of continuity ∂ρ + div(ρu) = 0, ∂t i.e., Dρ + ρ div u = 0, Dt

and the fact that ρ > 0, we see that a fluid is incompressible if and only if Dρ/Dt = 0, that is, the mass density is constant following the fluid . If the fluid is homogeneous, that is, ρ = constant in space, it also follows that the flow is incompressible if and only if ρ is constant in time. Problems involving inhomogeneous incompressible flow occur, for example, in oceanography. We shall now “solve” the equation of continuity by expressing ρ in terms of its value at t = 0, the flow map ϕ(x, t), and its Jacobian J(x, t). Indeed, set f = 1 in the transport theorem and conclude the equivalent condition for mass conservation, d ρ dV = 0 dt Wt and thus, ρ(x, t)dV =
Wt W0

ρ(x, 0) dV.

Changing variables, we obtain ρ(ϕ(x, t), t)J(x, t) dV =
W0 W0

ρ(x, 0) dV.

Because W0 is arbitrary, we get ρ(ϕ(x, t), t)J(x, t) = ρ(x, 0) as another form of mass conservation. As a corollary, a fluid that is homogeneous at t = 0 but is compressible will generally not remain homogeneous. However, the fluid will remain homogeneous if it is incompressible. The example ϕ((x, y, z), t) = ((1 + t)x, y, z) has J((x, y, z), t) = 1 + t so the flow is not incompressible, yet for ρ((x, y, z), t) = 1/(1 + t), one has mass conservation and homogeneity for all time.

12

1 The Equations of Motion

iii Conservation of Energy
So far we have developed the equations ρ and Dρ + ρ div u = 0 Dt (conservation of mass). Du = − grad p + ρb (balance of momentum) Dt

These are four equations if we work in 3-dimensional space (or n + 1 equations if we work in n-dimensional space), because the equation for Du/Dt is a vector equation composed of three scalar equations. However, we have five functions: u, ρ, and p. Thus, one might suspect that to specify the fluid motion completely, one more equation is needed. This is in fact true, and conservation of energy will supply the necessary equation in fluid mechanics. This situation is more complicated for general continua, and issues of general thermodynamics would need to be discussed for a complete treatment. We shall confine ourselves to two special cases here, and later we shall treat another case for an ideal gas. For fluid moving in a domain D, with velocity field u, the kinetic energy contained in a region W ⊂ D is Ekinetic = 1 2 ρ u
W 2

dV

where u 2 = (u2 + v 2 + w2 ) is the square length of the vector function u. We assume that the total energy of the fluid can be written as Etotal = Ekinetic + Einternal where Einternal is the internal energy , which is energy we cannot “see” on a macroscopic scale, and derives from sources such as intermolecular potentials and internal molecular vibrations. If energy is pumped into the fluid or if we allow the fluid to do work, Etotal will change. The rate of change of kinetic energy of a moving portion Wt of fluid is calculated using the transport theorem as follows: d d 1 Ekinetic = dt dt 2 = =
Wt

ρ u
Wt

2

dV

1 2

ρ
Wt

D u Dt

2

dV ∂u + (u · ∇)u ∂t dV.

ρ u·

1.1 Euler’s Equations

13

Here we have used the following Euclidean coordinate calculation 1 D u 2 Dt
2

=

∂ 1 ∂ 2 1 (u + v 2 + w2 ) + u (u2 + v 2 + w2 ) 2 ∂t 2 ∂x ∂ ∂ + v (u2 + v 2 + w2 ) + w (u2 + v 2 + w2 ) ∂y ∂z ∂u ∂v ∂w ∂u ∂v ∂w =u +v +w +u u +v +w ∂t ∂t ∂t ∂x ∂x ∂x ∂u ∂v ∂w ∂u ∂v ∂w +v u +v +w +w u +v +w ∂y ∂y ∂y ∂z ∂z ∂z ∂u =u· + u · (u · ∇)u). ∂t

A general discussion of energy conservation requires more thermodynamics than we shall need. We limit ourselves here to two examples of energy conservation; a third will be given in Chapter 3. 1 Incompressible Flows Here we assume all the energy is kinetic and that the rate of change of kinetic energy in a portion of fluid equals the rate at which the pressure and body forces do work: d Ekinetic = − dt pu · n dA +
∂Wt Wt

ρu · b dV.

By the divergence theorem and our previous formulas, this becomes ρ u·
Wt

∂u + u · ∇u ∂t

dV = −
Wt

(div(pu) − ρu · b) dV (u · ∇p − ρu · b) dV
Wt

=−

because div u = 0. The preceding equation is also a consequence of balance of momentum. This argument, in addition, shows that if we assume E = Ekinetic , then the fluid must be incompressible (unless p = 0). In summary, in this incompressible case, the Euler equations are: ρ Du = − grad p + ρb Dt Dρ =0 Dt div u = 0

with the boundary conditions u·n=0 on ∂D.

14

1 The Equations of Motion

2 Isentropic Fluids A compressible flow will be called isentropic if there is a function w, called the enthalpy , such that grad w = 1 grad p. ρ

This terminology comes from thermodynamics. We shall not need a detailed discussion of thermodynamics concepts in this book, and so it is omitted, except for a brief discussion of entropy in Chapter 3 in the context of ideal gases. For the readers’ convenience, we just make a few general comments. In thermodynamics one has the following basic quantities, each of which is a function of x, t depending on a given flow: p = pressure ρ = density T = temperature s = entropy w = enthalpy (per unit mass) = w − (p/ρ) = internal energy

(per unit mass).

These quantities are related by the First Law of Thermodynamics, which we accept as a basic principle:2 dw = T ds + 1 dp ρ (TD1)

The first law is a statement of conservation of energy; a statement equivalent to (TD1) is, as is readily verified, d = T ds + p dρ. ρ2 (TD2)

If the pressure is a function of ρ only, then the flow is clearly isentropic with s as a constant (hence the name isentropic) and ρ w=

p (λ) dλ, λ

which is the integrated version of dw = dp/ρ (see TD1). As above, the internal energy = w − (p/ρ) then satisfies d = (pdρ)/ρ2 (see TD2) or, as a function of ρ, ρ ∂ε p(λ) p = ρ2 , or dλ. = ∂p λ2
2 A. Sommerfeld [1964] Thermodynamics and Statistical Mechanics, reprinted by Academic Press, Chapters 1 and 4.

1.1 Euler’s Equations

15

For isentropic flows with p a function of ρ, the integral form of energy balance reads as follows: The rate of change of energy in a portion of fluid equals the rate at which work is done on it: d d Etotal = dt dt =
Wt 1 2ρ

u

2



dV (BE) pu · n dA.

Wt

ρu · b dV −
∂Wt

This follows from balance of momentum using our earlier expression for (d/dt)Ekinetic , the transport theorem, and p = ρ2 ∂ /∂ρ. Alternatively, one can start with the assumption that p is a function of ρ and then (BE) and balance of mass and momentum implies that p = ρ2 ∂ /∂ρ , which is equivalent to dw = dp/ρ, as we have seen.3 Euler’s equations for isentropic flow are thus ∂u + (u · ∇)u = −∇w + b, ∂t ∂ρ + div(ρu) = 0 ∂t in D, and u·n=0 on ∂D (or u · n = V · n if ∂D is moving with velocity V). Later, we will see that in general these equations lead to a well-posed initial value problem only if p (ρ) > 0. This agrees with the common experience that increasing the surrounding pressure on a volume of fluid causes a decrease in occupied volume and hence an increase in density. Gases can often be viewed as isentropic, with p = Aργ , where A and γ are constants and γ ≥ 1. Here, ρ w=

γAsγ−1 γAργ−1 ds = s γ−1

and

=

Aργ−1 . γ−1

Cases 1 and 2 above are rather opposite. For instance, if ρ = ρ0 is a constant for an incompressible fluid, then clearly p cannot be an invertible function of ρ. However, the case ρ = constant may be regarded as a limiting case p (ρ) → ∞. In case 2, p is an explicit function of ρ (and therefore
3 One can carry this even further and use balance of energy and its invariance under Euclidean motions to derive balance of momentum and mass, a result of Green and Naghdi. See Marsden and Hughes [1994] for a proof and extensions of the result that include formulas such as p = p2 ∂ε/∂p amongst the consequences as well.

16

1 The Equations of Motion

depends on u through the coupling of ρ and u in the equation of continuity); in case 1, p is implicitly determined by the condition div u = 0. We shall discuss these points again later. Finally, notice that in neither case 1 or 2 is the possibility of a loss of kinetic energy due to friction taken into account. This will be discussed at length in §1.3. Given a fluid flow with velocity field u(x, t), a streamline at a fixed time is an integral curve of u; that is, if x(s) is a streamline at the instant t, it is a curve parametrized by a variable, say s, that satisfies dx = u(x(s), t), t fixed. ds We define a fixed trajectory to be the curve traced out by a particle as time progresses, as explained at the beginning of this section. Thus, a trajectory is a solution of the differential equation dx = u(x(t), t) dt with suitable initial conditions. If u is independent of t (i.e., ∂t u = 0), streamlines and trajectories coincide. In this case, the flow is called stationary. Bernoulli’s Theorem In stationary isentropic flows and in the absence of external forces, the quantity
1 2

u

2

+w

is constant along streamlines. The same holds for homogeneous (ρ = constant in space = ρ0 ) incompressible flow with w replaced by p/ρ0 . The conclusions remain true if a force b is present and is conservative; i.e., b = −∇ϕ for some function ϕ, with w replaced by w + ϕ. Proof has From the table of vector identities at the back of the book, one
1 2 ∇(

u 2 ) = (u · ∇)u + u × (∇ × u). (u · ∇)u = −∇w

Because the flow is steady, the equations of motion give

and so ∇
1 2

u

2

+ w = u × (∇ × u). x(s2 )

Let x(s) be a streamline. Then
1 2

u

2

+w

x(s2 ) x(s1 )

= x(s1 ) x(s2 )



1 2

u

2

+ w · x (s) ds

= x(s1 )

(u × (∇ × u)) · x (s) ds = 0

1.1 Euler’s Equations

17

because x (s) = u(x(s)) is orthogonal to u × (∇ × u). See Exercise 1.1-3 at the end of this section for another view of why the combination 1 u 2 + w is the correct quantity in Bernoulli’s theorem. 2 We conclude this section with an example that shows the limitations of the assumptions we have made so far. Example Consider a fluid-filled channel, as in Figure 1.1.5. y flow direction

channel with

p1 > p2 x

pressure =

p1 L

pressure =

p2

0

Figure 1.1.5. Fluid flow in a channel.

Suppose that the pressure p1 at x = 0 is larger than that at x = L so the fluid is pushed from left to right. We seek a solution of Euler’s incompressible homogeneous equations in the form u(x, y, t) = (u(x, t), 0) and p(x, y, t) = p(x).

Incompressibility implies ∂x u = 0. Thus, Euler’s equations become ρ0 ∂t u = 2 −∂x p. This implies that ∂x p = 0, and so p(x) = p1 − p1 − p2 L x.

Substitution into ρ0 ∂t u = −∂x p and integration yields u= p1 − p 2 t + constant. ρ0 L

This solution suggests that the velocity in channel flow with a constant pressure gradient increases indefinitely. Of course, this cannot be the case in a real flow; however, in our modeling, we have not yet taken friction into account. The situation will be remedied in §1.3.

Exercises
Exercise 1.1-1 Prove the following properties of the material derivative

18

1 The Equations of Motion

(i) (ii) (iii)

D Df Dg (f + g) = + , Dt Dt Dt D Dg Df (f · g) = f +g Dt Dt Dt (Leibniz or product rule),

D Dg (h ◦ g) = (h ◦ g) (chain rule). Dt Dt Exercise 1.1-2 Use the transport theorem to establish the following formula of Reynolds: d dt f (x, t) dV =
Wt Wt

∂f (x, t) dV + ∂t

f u · n dA.
∂Wt

Interpret the result physically. Exercise 1.1-3 Consider isentropic flow without any body force. Show that for a fixed volume W in space (not moving with the flow). d dt
1 2ρ

u

2



dV = −
∂W

ρ

W

1 2

u

2

+ w u · n dA.

Use this to justify the term energy flux vector for the vector function ρu 1 u 2 + w and compare with Bernoulli’s theorem. 2

1.2

Rotation and Vorticity ξ = ∇ × u = (∂y w − ∂z v, ∂z u − ∂x w, ∂x v − ∂y u)

If the velocity field of a fluid is u = (u, v, w), then its curl,

is called the vorticity field of the flow. We shall now demonstrate that in a small neighborhood of each point of the fluid, u is the sum of a (rigid ) translation, a deformation (defined later ), and a (rigid ) rotation with rotation vector ξ/2. This is in fact a general statement about vector fields u on R3 ; the specific features of fluid mechanics are irrelevant for this discussion. Let x be a point in R3 , and let y = x + h be a nearby point. What we shall prove is that u(y) = u(x) + D(x) · h + 1 ξ(x) × h + O(h2 ), 2 (1.2.1)

where D(x) is a symmetric 3 × 3 matrix and h2 = h 2 is the squared length of h. We shall discuss the meaning of the several terms later. Proof of Formula (1.2.1) Let  ∂x u ∇u =  ∂x v ∂x w  ∂z u ∂z v  ∂z w

∂y u ∂y v ∂y w

1.2 Rotation and Vorticity

19

denote the Jacobian matrix of u. By Taylor’s theorem, u(y) = u(x) + ∇u(x) · h + O(h2 ), (1.2.2)

where ∇u(x) · h is a matrix multiplication, with h regarded as a column vector. Let D = 1 ∇u + (∇u)T , 2 where
T

denotes the transpose, and S=
1 2

∇u − (∇u)T . (1.2.3)

Thus, ∇u = D + S. It is easy to check that the coordinate expression for S is   ξ2 0 −ξ3 1 ξ3 0 −ξ1  S= 2 −ξ2 ξ1 0 and that S · h = 1 ξ × h, 2 (1.2.4) where ξ = (ξ1 , ξ2 , ξ3 ). Substitution of (1.2.3) and (1.2.4) into (1.2.2) yields (1.2.1). Because D is a symmetric matrix, D(x) · h = gradh ψ(x, h), where ψ is the quadratic form associated with D; i.e., ψ(x, h) =
1 2

D(x) · h, h ,

where , is the inner product of R3 . We call D the deformation tensor. We now discuss its physical interpretation. Because D is symmetric, there ˜ ˜ ˜ is, for x fixed, an orthonormal basis e1 , e2 , e3 in which D is diagonal:   d1 0 0 D =  0 d2 0  . 0 0 d3 Keep x fixed and consider the original vector field as a function of y. The motion of the fluid is described by the equations dy = u(y). dt If we ignore all terms in (1.2.1) except D · h, we find dy = D · h, dt i.e., dh = D · h. dt

20

1 The Equations of Motion

This vector equation is equivalent to three linear differential equations that ˜ ˜ ˜ separate in the basis e1 , e2 , e3 : ˜ d hi ˜ = d i hi , dt i = 1, 2, 3.

˜ The rate of change of a unit length along the ei axis at t = 0 is thus di . The vector field D · h is thus merely expanding or contracting along each ˜ of the axes ei —hence, the name “deformation.” The rate of change of the ˜ ˜ ˜ ˜ ˜ ˜ volume of a box with sides of length h1 , h2 , h3 parallel to the e1 , e2 , e3 axes is ˜ ˜ ˜ d ˜ ˜ ˜ d h1 ˜ ˜ ˜ ˜ d h2 h 3 + h 1 h2 d h 3 ˜ ˜ (h1 h2 h3 ) = h2 h 3 + h 1 dt dt dt dt ˜ ˜ ˜ = (d1 + d2 + d3 )(h1 h2 h3 ). However, the trace of a matrix is invariant under orthogonal transformations. Hence, d1 + d2 + d3 = trace of D = trace of
1 2

(∇u) + (∇u)T = div u.

This confirms the fact proved in §1.1 that volume elements change at a rate proportional to div u. Of course, the constant vector field u(x) in formula (1.2.1) induces a flow that is merely a translation by u(x). The other term, 1 2 ξ(x) × h, induces a flow dh = 1 ξ(x) × h, 2 dt (x fixed).

The solution of this linear differential equation is, by elementary vector calculus, h(t) = R(t, ξ(x))h(0), where R(t, ξ(x)) is the matrix that represents a rotation through an angle t about the axis ξ(x) (in the oriented sense). Because rigid motion leaves volumes invariant, the divergence of 1 ξ(x) × h is zero, as may also be 2 checked by noting that S has zero trace. This completes our derivation and discussion of the decomposition (1.2.1). We remarked in §1.1 that our assumptions so far have precluded any tangential forces, and thus any mechanism for starting or stopping rotation. Thus, intuitively, we might expect rotation to be conserved. Because rotation is intimately related to the vorticity as we have just shown, we can expect the vorticity to be involved. We shall now prove that this is so. Let C be a simple closed contour in the fluid at t = 0. Let Ct be the contour carried along by the flow. In other words, Ct = ϕt (C),

1.2 Rotation and Vorticity

21

D C Ct

Figure 1.2.1. Kelvin’s circulation theorem.

where ϕt is the fluid flow map discussed in §1.1 (see Figure 1.2.1). The circulation around Ct is defined to be the line integral

ΓC t =
Ct

u · ds.

Kelvin’s Circulation Theorem For isentropic flow without external forces, the circulation, ΓCt is constant in time. For example, we note that if the fluid moves in such a way that Ct shrinks in size, then the “angular” velocity around Ct increases. The proof of Kelvin’s circulation theorem is based on a version of the transport theorem for curves. Lemma Let u be the velocity field of a flow and C a closed loop, with Ct = ϕt (C) the loop transported by the flow (Figure 1.2.1). Then d dt u · ds =
Ct Ct

Du ds. Dt

(1.2.5)

Proof Let x(s) be a parametrization of the loop C, 0 ≤ s ≤ 1. Then a parameterization of Ct is ϕ(x(s), t), 0 ≤ s ≤ 1. Thus, by definition of the line integral and the material derivative, d dt u · ds =
Ct

d dt
1 0

1

u(ϕ(x(s), t), t) ·
0

∂ ϕ(x(s), t) ds ∂s

= +

Du ∂ (ϕ(x(s), t), t) · ϕ(x(s), t) ds Dt ∂s
1

u(ϕ(x(s), t), t) ·
0

∂ ∂ ϕ(x(s), t) ds. ∂t ∂s

22

1 The Equations of Motion

Because ∂ϕ/∂t = u, the second term equals
1

u(ϕ(x(s), t), t) ·
0

∂ u(ϕ(x(s), t), t) ds ∂s = 1 2
1 0

∂ (u · u)(ϕ(x(s), t), t) ds = 0 ∂s

(since Ct is closed). The first term equals Du ds, Dt

Ct

so the lemma is proved. Proof of the Circulation Theorem Using the lemma and the fact that Du/Dt = −∇w (the flow is isentropic and without external forces), we find d d ΓC = dt t dt u ds =
Ct Ct

Du ds Dt ∇w · ds = 0 (since Ct is closed).
Ct

=−

We now use Stokes’ theorem, which will bring in the vorticity. If Σ is a surface whose boundary is an oriented closed oriented contour C, then Stokes’ theorem yields (see Figure 1.2.2) ΓC =
C

u · ds =
Σ

(∇ × u) · n dA =
Σ

ξ · dA.

dA = n dA

Σ

C

Figure 1.2.2. The circulation around C is the integral of the vorticity over Σ.

Thus, as a corollary of the circulation theorem, we can conclude that the flux of vorticity across a surface moving with the fluid is constant in time.

1.2 Rotation and Vorticity

23

ξ x flow x

ξ

u

S vortex sheet

L vortex line

Figure 1.2.3. Vortex sheets and lines remain so under the flow.

By definition, a vortex sheet (or vortex line) is a surface S (or a curve L) that is tangent to the vorticity vector ξ at each of its points (Figure 1.2.3). Proposition If a surface (or curve) moves with the flow of an isentropic fluid and is a vortex sheet (or line) at t = 0, then it remains so for all time. Proof Let n be the unit normal to S, so that at t = 0, ξ · n = 0 by hypothesis. By the circulation theorem, the flux of ξ across any portion ˜ S ⊂ S at a later time is also zero, i.e., ξ · n dA = 0.

˜ St

It follows that ξ · n = 0 identically on St , so S remains a vortex sheet. One can show (using the implicit function theorem) that if ξ(x) = 0, then, locally, a vortex line is the intersection of two vortex sheets. Next, we show that the vorticity (per unit mass), that is, ω = ξ/ρ, is propagated by the flow (see Figure 1.2.4). This fact can also be used to give another proof of the preceding theorem. We assume we are in three dimensions; the two-dimensional case will be discussed later. Proposition For isentropic flow (in the absence of external forces) with ξ = ∇ × u and ω = ξ/ρ, we have Dω − (ω · ∇)u = 0 Dt (1.2.6)

24

1 The Equations of Motion

ω at time 0

ω is dragged by the flow

x ϕ(x,t )

ω at time t

Figure 1.2.4. The vorticity is transported by the Jacobian matrix of the flow map.

and ω(ϕ(x, t), t) = ∇ϕt (x) · ω(x, 0), (1.2.7) where ϕt is the flow map (see §1.1) and ∇ϕt is its Jacobian matrix. Proof Start with the following vector identity (see the table of vector identities at the back of the book)
1 2 ∇(u

· u) = u × curl u + (u · ∇)u.

Substituting this into the equations of motion yields ∂u 1 + 2 ∇(u · u) − u × curl u = −∇w. ∂t Taking the curl and using the identity ∇ × ∇f = 0 gives ∂ξ − curl(u × ξ) = 0. ∂t Using the identity (also from the back of the book) curl(F × G) = F div G − G div F + (G · ∇)F − (F · ∇)G for the curl of a vector product, gives ∂ξ − [ (u(∇ · ξ) − ξ(∇ · u) + ξ · ∇)u − (u · ∇)ξ ] = 0, ∂t that is, Dξ − (ξ · ∇)u + ξ(∇ · u) = 0, Dt since ξ is divergence free. Also, Dω D = Dt Dt ξ ρ = 1 Dξ ξ + (∇ · u) ρ Dt ρ (1.2.8)

(1.2.9)

by the continuity equation. Substitution of (1.2.8) into (1.2.9) yields (1.2.6).

1.2 Rotation and Vorticity

25

To prove (1.2.7), let F(x, t) = ω(ϕ(x, t), t) and G(x, t) = ∇ϕt (x) · ω(x, 0). By (1.2.6), ∂F/∂t = (F · ∇)u. On the other hand, by the chain rule: ∂G ∂ϕ =∇ (x, t) · ω(x, 0) = ∇(u(ϕ(x, t), t)) · ω(x, 0) ∂t ∂t = (∇u) · ∇ϕt (x) · ω(x, 0) = (G · ∇)u Thus, F and G satisfy the same linear first-order differential equation. Because they coincide at t = 0 and solutions are unique, they are equal. The reader may wish to compare (1.2.7) with the formula ρ(x, 0) = ρ(ϕ(x, t), t)J(x, t) (1.2.10)

proved in §1.1. As an exercise, we invite the reader to prove the preservation of vortex sheets and lines by the flow using (1.2.7) and (1.2.10). For two-dimensional flow, where u = (u, v, 0), ξ has only one component; ξ = (0, 0, ξ). The circulation theorem now states that if Σt is any region in the plane that is moving with the fluid, then ξ dA = constant in time.
Σt

(1.2.11)

In fact, one can say more using (1.2.7). In two dimensions, (1.2.7) specializes to ξ ξ (ϕ(x, t), t) = (x, 0), (1.2.7) ρ ρ that is, ξ/ρ is propagated as a scalar by the flow. Employing (1.2.10) and the change of variables theorem gives (1.2.11) as a special case. In three-dimensional flows, the relation (1.2.7) allows rather complicated behavior. We shall now discuss the three-dimensional geometry a bit further. A vortex tube consists of a two-dimensional surface S that is nowhere tangent to ξ, with vortex lines drawn through each point of the bounding curve C of S. These vortex lines are integral curves of ξ and are extended as far as possible in each direction. See Figure 1.2.5. In fluid mechanics it is customary to be sloppy about this definition and make tacit assumptions to the effect that the tube really “looks like” a tube. More precisely, we assume S is diffeomorphic to a disc (i.e., related to a disc by a one-to-one invertible differentiable transformation) and that the resulting tube is diffeomorphic to the product of the disc and the real line. This tacitly assumes that ξ has no zeros (of course, ξ could have zeros!).

26

1 The Equations of Motion

S

C

vortex line

Figure 1.2.5. A vortex tube consists of vortex lines drawn through points of C.

Helmholtz’s Theorem

Assume the fluid is isentropic. Then

(a) If C1 and C2 are any two curves encircling the vortex tube, then u · ds =
C1 C2

u · ds.

This common value is called the strength of the vortex tube. (b) The strength of the vortex tube is constant in time, as the tube moves with the fluid. Proof (a) Let C1 and C2 be oriented as in Figure 1.2.6.
V = region enclosed

S2 C2

S

S1
C1

Figure 1.2.6. A vortex tube enclosed between two curves, C1 and C2 .

The lateral surface of the vortex tube enclosed between C1 and C2 is denoted by S, and the end faces with boundaries C1 and C2 are denoted by S1 and S2 , respectively. Since ξ is tangent to the lateral surface, S is a

1.2 Rotation and Vorticity

27

vortex sheet. Let V denote the region of the vortex tube between C1 and C2 and Σ = S ∪ S1 ∪ S2 denote the boundary of V . By Gauss’ theorem, 0=
V

∇ · ξ dx =
Σ

ξ · dA =

S1 ∪S2

ξ · dA +
S

ξ · dA.

By Stokes’ theorem u · ds =
C1 S1

ξ · dA

and
C2

u · ds = −
S2

ξ · dA,

so (a) holds. Part (b) now follows from Kelvin’s circulation theorem. Observe that if a vortex tube gets stretched and its cross-sectional area decreases, then the magnitude of ξ must increase. Thus, the stretching of vortex tubes can increase vorticity, but it cannot create it. A vortex tube with nonzero strength cannot “end” in the interior of the fluid. It either forms a ring (such as the smoke from a cigarette), extends to infinity, or is attached to a solid boundary. The usual argument supporting this statement goes like this: suppose the tube ended at a certain cross section S, inside the fluid. Because the tube cannot be extended, we must have ξ = 0 on C1 . Thus, the strength is zero—a contradiction. This “proof” is hopelessly incomplete. First of all, why should a vortex tube end in a nice regular way on a surface? Why can’t it split in two, as in Figure 1.2.7? There is no a priori reason why this sort of thing cannot happen unless we merely exclude it by tacit assumption.4 . In particular, note that the assertion often made that a vortex line cannot end in the fluid is clearly false if we allow ξ to have zeros and probably is false even if ξ has no zeros (an orbit of a vector field can wander around forever without accumulating at an endpoint—as with a line with irrational slope on a torus) Thus, our assertion about vortex tubes “ending” is correct if we interpret “ending” properly. But the reader is cautioned that this may not be all that can happen, and that this time-honored statement is not at all a proved theorem. The difference between the two-dimensional and three-dimensional conservation laws for vorticity is very important. The conservation of vorticity (1.2.7) in two dimensions is a helpful tool in establishing a rigorous theory of existence and uniqueness of the Euler (and later Navier–Stokes) equations. The lack of the same kind of conservation in three dimensions is a major obstacle to the rigorous understanding of crucial properties of the solutions of the equations of fluid dynamics. The main point here is to get existence theorems for all time. At the moment, it is known only in two dimension that all time smooth solutions exist.
4 H. Lamb [1895] Mathematical Theory of the Motion of Fluids, Cambridge Univ. Press, p. 149.

28

1 The Equations of Motion

a zero of ξ

C S

C1 P

C2

this vortex line ends at P
Figure 1.2.7. Can this be a vortex tube generated by S? Is the circulation around C1 equal to that around C2 ?

Our last main goal in this section is to develop the vorticity equation somewhat further for the important special case of incompressible flow. For two-dimensional homogeneous incompressible flow, the vorticity equation is Dξ (1.2.12) = ∂t ξ + (u · ∇)ξ = 0, Dt where ξ = ξ(x, y, t) = ∂x v −∂y u is the (scalar) vorticity field of the flow and u, v are the components of u. Assume that the flow is contained in some plane domain D with a fixed boundary ∂D, with the boundary condition u·n=0 on ∂D, (1.2.13)

where n is the unit outward normal to ∂D. Let us assume D is simply connected (i.e., has no “holes”). Then, by incompressibility, ∂x u = −∂y v, and so from vector calculus there is a scalar function ψ(x, y, t) on D unique up to an additive constant such that u = ∂y ψ and v = −∂x ψ. (1.2.14)

The function ψ is the stream function for fixed t; streamlines lie on level curves of ψ. Indeed, let (x(s), y(s)) be a streamline, so x = u(x, y) and y = v(x, y). Then d ψ(x(s), y(s), t) = ∂x ψ · x + ∂y ψ · y = −vu + uv = 0. ds In particular, by (1.2.13), ∂D lies on a level curve of ψ, and we can adjust the constant so that ψ(x, y, t) = 0 for (x, y) ∈ ∂D.

This convention and (1.2.14) determine ψ uniquely. (∂D need not be a whole streamline, but can be composed of streamlines separated by zeros of u, that is, by stagnation points.) The scalar vorticity is now given by
2 2 ξ = ∂x v − ∂y u = −∂x ψ − ∂y ψ = −∆ψ,

1.2 Rotation and Vorticity

29

2 2 where ∆ = ∂x + ∂y is the Laplace operator in the plane. We can summarize the equations for ξ for two-dimensional incompressible flow as follows:  Dξ  Dt ≡ ∂t ξ + (u · ∇)ξ = 0,     ∆ψ = −ξ,     with (1.2.15)  ψ = 0 on ∂D,      and with    u = ∂y ψ and v = −∂x ψ.

These equations completely determine the flow. Note that given ξ, the function ψ is determined by ∆ψ = −ξ and the boundary conditions, and hence u by the last equations in (1.2.15). Thus, ξ completely determines ∂t ξ and hence the evolution of ξ and, through it, ψ and u. Another remark is useful: (u · ∇)ξ = u∂x ξ + v∂y ξ = (∂y ψ)(∂x ξ) − (∂x ψ)(∂y ξ) = det ∂x ξ ∂x ψ ∂y ξ ∂y ψ = J(ξ, ψ),

the Jacobian of ξ and ψ. Thus, the flow is stationary (time independent) if and only if ξ and ψ are functionally dependent. (If functional dependence holds at one instant it will hold for all time as a consequence.) Example Suppose at t = 0 the stream function ψ(x, y) is a function only of the radial distance r = (x2 + y 2 )1/2 . Thus, the streamlines are concentric circles. Write ψ(x, y) = ψ(r) and assume ψr > 0. The velocity vector is given by u = ∂y ψ = ∂r ψ∂y r = y ∂r ψ, r x v = −∂x ψ = −∂r ψ∂x r = − ∂r ψ, r (1.2.16) (1.2.17)

that is, u is tangent to the circle of radius r with magnitude |∂r ψ| and oriented clockwise if ψr > 0 and counterclockwise if ψr < 0. Next, observe that 1 ∂ ∂ψ ξ = −∆ψ = − r , r ∂r ∂r a function of r alone. Because ψr = 0, r is a function of ψ so ξ is also a function of ψ. Thus, J(ξ, ψ) = 0. Hence, motion in concentric circles with u defined as above is a solution of the two-dimensional stationary incompressible equations of ideal flow.

30

1 The Equations of Motion

is

For three-dimensional incompressible ideal flow, the analogue of (1.2.15)  Dξ   − (ξ · ∇)u = 0,   Dt (1.2.18) ∆A = −ξ, div A = 0,     u = ∇ × A.

Here we used ∇·u = 0 to write u = ∇×A, where div A = 0. (This requires not that D be simply connected, but that it not have any “solid holes” in it; for instance, if D is convex, this will hold.) Then ξ = curl u = curl(curl A) = −∆A + ∇(div A) = −∆A. One of the troubles with (1.2.18) is that given ξ, the vector field A is not uniquely determined (we cannot impose boundary condition such as A = 0 on ∂D because A need not be constant on ∂D as was the case with ψ).

Exercises
Exercise 1.2-1 Derive a formula akin to the transport theorem and Kelvin’s circulation theorem for d v · n dA, dt St where St is a moving surface and v is a vector field. Exercise 1.2-2 Couette flow. Let Ω be the region between two concentric cylinders of radii R1 and R2 , where R1 < R2 . Define v in cylindrical coordinates by vr = 0, vz = 0, and vθ = where A=− Show that (i) v is a stationary solution of Euler’s equations with ρ = 1; (ii) ω = ∇ × v = (0, 0, 2B); (iii) the deformation tensor is D=− A r2 0 1 1 0
2 2 R1 R2 (ω2 − ω1 ) 2 1 R2 − R2

A + Br, r and B = −
2 2 R1 ω 1 − R2 ω 2 . 2 − R2 R2 1

and discuss its physical meaning;

1.3 The Navier–Stokes Equations

31

(iv)

the angular velocity of the flow on the two cylinders is ω1 and ω2 .

1.3

The Navier–Stokes Equations

In §1.1 we defined an ideal fluid as one in which forces across a surface were normal to that surface. We now consider more general fluids. To understand the need for the generalization beyond the examples already given, consider the situation shown in Figure 1.3.1. Here the velocity field u is parallel to a surface S but jumps in magnitude either suddenly or rapidly as we cross S. If the forces are all normal to S, there will be no transfer of momentum between the fluid volumes denoted by B and B in Figure 1.3.1. However, if we remember the kinetic theory of matter, we see that this is actually unreasonable. Faster molecules from above S will diffuse across S and impart momentum to the fluid below, and, likewise, slower molecules from below S will diffuse across S to slow down the fluid above S. For reasonably fast changes in velocity over short distance, this effect is important.5
B' u u B
Figure 1.3.1. Faster molecules in B can diffuse across S and impart momentum to B.

S

We thus change our previous definition. Instead of assuming that force on S per unit area = −p(x, t)n, where n is the normal to S, we now assume that force on S per unit area = −p(x, t)n + σ(x, t) · n, (1.3.1)

where σ is a matrix called the stress tensor, about which some assumptions will have to be made. The new feature is that σ·n need not be parallel to n. The separation of the forces into pressure and other forces in (1.3.1) is somewhat ambiguous because σ · n may contain a component parallel to n. This issue will be resolved later when we give a more definite functional form to σ.
5 For more information, see J. Jeans [1867] An Introduction to the Kinetic Theory of Gases, Cambridge Univ. Press.

32

1 The Equations of Motion

As before, Newton’s second law states that the rate of change of any moving portion of fluid Wt equals the force acting on it (balance of momentum): d ρu dV = − (p · n − σ · n) dA dt Wt ∂Wt (compare (BM3) in §1.1). Thus, we see that σ modifies the transport of momentum across the boundary of Wt . We will choose σ so that it approximates in a reasonable way the transport of momentum by molecular motion. One can legitimately ask why the force (1.3.1) acting on S should be a linear function of n. In fact, if one just assumes the force is a continuous function of n, then, using balance of momentum, one can prove it is linear in n. This result is called Cauchy’s Theorem.6 Our assumptions on σ are as follows: 1. σ depends linearly on the velocity gradients ∇u that is, σ is related to ∇u by some linear transformation at each point. 2. σ is invariant under rigid body rotations, that is, if U is an orthogonal matrix, σ(U · ∇u · U−1 ) = U · σ(∇u) · U−1 . This is reasonable, because when a fluid undergoes a rigid body rotation, there should be no diffusion of momentum. 3. σ is symmetric. This property can be deduced as a consequence of balance of angular momentum.7 Since σ is symmetric, if follows from properties 1 and 2 that σ can depend only on the symmetric part of ∇u; that is, on the deformation D. Because σ is a linear function of D, σ and D commute and so can be simultaneously diagonalized. Thus, the eigenvalues of σ are linear functions of those of D. By property 2, they must also be symmetric because we can choose U to permute two eigenvalues of D (by rotating through an angle π/2 about an eigenvector), and this must permute the corresponding eigenvalues of σ. The only linear functions that are symmetric in this sense are of the form σi = λ(d1 + d2 + d3 ) + 2µdi , i = 1, 2, 3,

where σi are the eigenvalues of σ, and di are those of D. This defines the constants λ and µ. Recalling that d1 + d2 + d3 = div u, we can use property 2 to transform σi back to the usual basis and deduce that σ = λ(div u)I + 2µ D,
6 For 7 Op.

(1.3.2)

a proof and further references, see, for example, Marsden and Hughes [1994]. cit.

1.3 The Navier–Stokes Equations

33

where I is the identity. We can rewrite this by putting all the trace in one term: σ = 2µ[D − 1 (div u)I] + ζ(div u)I (1.3.2) 3 where µ is the first coefficient of viscosity , and ζ = λ + 2 µ is the 3 second coefficient of viscosity. If we employ the transport theorem and the divergence theorem, as we did in connection with (BM3), balance of momentum yields the Navier– Stokes equations, Du = −∇p + (λ + µ)∇(div u) + µ∆u Dt ∂2 ∂2 ∂2 + 2+ 2 ∂x2 ∂y ∂z

ρ where

(1.3.3)

∆u =

u

is the Laplacian of u. Together with the equation of continuity and an energy equation, (1.3.3) completely describes the flow of a compressible viscous fluid. In the case of incompressible homogeneous flow ρ = ρ0 = constant, the complete set of equations becomes the Navier–Stokes equations for incompressible flow, Du = − grad p + ν∆u Dt div u = 0

(1.3.4)

where ν = µ/ρ0 is the coefficient of kinematic viscosity, and p = p/ρ0 . These equations are supplemented by boundary conditions. For Euler’s equations for ideal flow we use u · n = 0, that is, fluid does not cross the boundary but may move tangentially to the boundary. For the Navier– Stokes equations, the extra term ν∆u raises the number of derivatives of u involved from one to two. For both experimental and mathematical reasons, this is accompanied by an increase in the number of boundary conditions. For instance, on a solid wall at rest we add the condition that the tangential velocity also be zero (the “no-slip condition”), so the full boundary conditions are simply u = 0 on solid walls at rest. The mathematical need for extra boundary conditions hinges on their role in proving that the equations are well posed; that is, that a unique solution exists and depends continuously on the initial data. In three dimensions, it is known that smooth solutions to the incompressible equations

34

1 The Equations of Motion

exist for a short time and depend continuously on the initial data.8 It is a major open problem in fluid mechanics to prove or disprove that solutions of the incompressible equations exist for all time. In two dimensions, solutions are known to exist for all time, for both viscous and inviscid flow9 . In any case, adding the tangential boundary condition is crucial for viscous flow. The physical need for the extra boundary conditions comes from simple experiments involving flow past a solid wall. For example, if dye is injected into flow down a pipe and is carefully watched near the boundary, one sees that the velocity approaches zero at the boundary to a high degree of precision. The no-slip condition is also reasonable if one contemplates the physical mechanism responsible for the viscous terms, namely, molecular diffusion. Our opening example indicates that molecular interaction between the solid wall with zero tangential velocity (or zero average velocity on the molecular level) should impart the same condition to the immediately adjacent fluid. Another crucial feature of the boundary condition u = 0 is that it provides a mechanism by which a boundary can produce vorticity in the fluid. We shall describe this in some detail in Chapter 2. Next, we shall discuss some scaling properties of the Navier–Stokes equations with the aim of introducing a parameter (the Reynolds number) that measures the effect of viscosity on the flow. For a given problem, let L be a characteristic length and U a characteristic velocity . These numbers are chosen in a somewhat arbitrary way. For example, if we consider flow past a sphere, L could be either the radius or the diameter of the sphere, and U could be the magnitude of the fluid velocity at infinity. L and U are merely reasonable length and velocity scales typical of the flow at hand. Their choice then determines a time scale by T = L/U . We can measure x, u, and t as fractions of these scales, by changing variables and introducing the following dimensionless quantities u x t u = , x = , and t = . U L T The x component of the (homogeneous) incompressible Navier–Stokes equation is ∂u ∂u ∂u ∂u 1 ∂p ∂2u ∂2u ∂2u + 2 + 2 . +u +v +w =− +ν ∂t ∂x ∂y ∂z ρ0 ∂x ∂x2 ∂y ∂z
8 For a review of much of what is known, see O. A. Ladyzhenskaya [1969] The Mathematical Theory of Viscous Incompressible Flow , Gordon and Breach. See also R. Temam [1977] Navier–Stokes Equations, North Holland. 9 Op. cit. and W. Wolibner, Math. Zeit. 37 [1933], 698–726; V. Judovich, Mat. Sb. N.S. 64 [1964], 562–588; and T. Kato, Arch. Rational Mech. Anal. 25 [1967], 188–200.

1.3 The Navier–Stokes Equations

35

The change of variables produces ∂(u U ) ∂t ∂(u U ) ∂x ∂(u U ) ∂y ∂(u U ) ∂w + Uu + Uv + Uw ∂t ∂t ∂x ∂x ∂y ∂y ∂z ∂z 1 ∂p ∂x ∂ 2 (u U ) ∂ 2 (u U ) ∂ 2 (u U ) , =− +ν + + ρ0 ∂x ∂x ∂(Lx )2 ∂(Ly )2 ∂(Lz )2 U2 L ∂u ∂u ∂u ∂u +u +v +w ∂t ∂x ∂y ∂z 2 2 U ∂(p/(ρ0 U )) ∂2u U ∂2u ∂2u =− ν . + 2 + 2 + L ∂x L2 ∂x ∂y ∂z 2

Similar equations hold for the y and z components. If we combine all three components and divide out by U 2 /L, we obtain ∂u ν + (u · ∇ )u = − grad p + ∆u, ∂t LU where p = p/(ρ0 U 2 ). Incompressibility still reads div u = 0. The equations (1.3.5) are the Navier–Stokes equations in dimensionless variables. We define the Reynolds number R to be the dimensionless number LU R= . ν For example, consider two flows past two spheres centered at the origin but with differing radii, one with a fluid where U∞ = 10 km/hr past a sphere of radius 10 m and the other with the same fluid but with U∞ = 100 km/hr and radius = 1 m. If we choose L to be the radius and U to be the velocity U∞ at infinity, then the Reynolds number is the same for each flow. The equations satisfied by the dimensionless variables are thus identical for the two flows. Two flows with the same geometry and the same Reynolds number are called similar . More precisely, let u1 and u2 be two flows on regions D1 and D2 that are related by a scale factor λ so that L1 = λL2 . Let choices of U1 and U2 be made for each flow, and let the viscosities be ν1 and ν2 respectively. If L1 U1 L2 U2 = , R1 = R2 , i.e., ν1 ν2 then the dimensionless velocity fields u1 and u2 satisfy exactly the same equation on the same region. Thus, we can conclude that u1 may be obtained from a suitably rescaled solution u2 ; in other words, u1 and u2 are similar. (1.3.5)

36

1 The Equations of Motion

This idea of the similarity of flows is used in the design of experimental models. For example, suppose we are contemplating a new design for an aircraft wing and we wish to know the behavior of a fluid flow around it. Rather than build the wing itself, it may be faster and more economical to perform the initial tests on a scaled-down version. We design our model so that it has the same geometry as the full-scale wing and we choose values for the undisturbed velocity, coefficient of viscosity, and so on, such that the Reynolds number for the flow in our experiment matches that of the actual flow. We can then expect the results of our experiment to be relevant to the actual flow over the full-scale wing. We shall be especially interested in cases where R is large. We stress that one cannot say that if ν is small, then viscous effects are unimportant, because such a comment fails to consider the other dimensions of the problem, that is, “ν is small” is not a physically meaningful statement unless some scaling is chosen, but “1/R is small” is a meaningful statement. As with incompressible ideal flow, the pressure p in incompressible viscous flow is determined through the equation div u = 0. We now shall explore the role of the pressure in incompressible flow in more depth. Let D be a region in space (or in the plane) with smooth boundary ∂D. We shall use the following decomposition theorem. Helmholtz–Hodge Decomposition Theorem can be uniquely decomposed in the form w = u + grad p, A vector field w on D

(1.3.6)

where u has zero divergence and is parallel to ∂D; that is, u · n = 0 on ∂D. Proof First of all, we establish the orthogonality relation u · grad p dV = 0.
D

Indeed, by the identity div(pu) = (div u)p + u · grad p, the divergence theorem, and div u = 0, we get u · grad p dV =
D D

div(pu) dV =
∂D

pu · n dA = 0,

because u · n = 0 on ∂D. We use this orthogonality to prove uniqueness. Suppose the w = u1 + grad p1 = u2 + grad p2 . Then 0 = u1 − u2 + grad(p1 − p2 ).

1.3 The Navier–Stokes Equations

37

Taking the inner product with u1 − u2 and integrating, we get 0=
D

u1 − u2

2

+ (u1 − u2 ) · grad(p1 − p2 ) dV =
D

u1 − u2

2

dV

by the orthogonality relation. It follows that u1 = u2 , and so, grad p1 = grad p2 (which is the same thing as p1 = p2 + constant). If w = u + grad p, notice that div w = div grad p = ∆p and that w · n = n · grad p. We use this remark to prove existence. Indeed, given w, let p be defined by the solution to the Neumann problem ∆p = div w in D, with ∂p = w · n on ∂D. ∂n

It is known10 that the solution to this problem exists and is unique up to the addition of a constant to p. With this choice of p, define u = w−grad p. Then, clearly u has the desired properties div u = 0, and also u · n = 0 by construction of p. The situation is shown schematically in Figure 1.3.2. gradient fields

vector fields that are divergence free and parallel to the boundary

Figure 1.3.2. Decomposing a vector field into a divergence-free and gradient part.

It is natural to introduce the operator P, an orthogonal projection operator, which maps w onto its divergence-free part u. By the preceding theorem, P is well defined. Notice that by construction P is a linear operator and that w = Pw + grad p. (1.3.7) Also notice that Pu = u
10 See R. Courant and D. Hilbert [1953], Methods of Mathematical Physics, Wiley. The equation ∆p = f, ∂p/∂n = g has a solution unique up to a constant if and only if D f dV = ∂D g dA. The divergence theorem ensures that this condition is satisfied in our case.

38

1 The Equations of Motion

provided div u = 0 and u · n = 0, and that P(grad p) = 0. Now we apply these ideas to the incompressible Navier–Stokes equations (1.3.5). If we apply the operator P to both sides, we obtain P(∂t u + grad p) = P −(u · ∇)u + 1 ∆u . R

Because u is divergence-free and vanishes on the boundary, the same is true of ∂t u (if u is smooth enough). Thus, by (1.3.7), P∂t u = ∂t u. Because P(grad p) = 0, we get ∂t u = P −u · ∇u + 1 ∆u . R (1.3.8)

Although ∆u is divergence free, it need not be parallel to the boundary and so we cannot simply write P∆u = 0. This form (1.3.8) of the Navier– Stokes equations eliminates the pressure and expresses ∂t u in terms of u alone. The pressure can then be recovered as the gradient part of −u · ∇u + 1 ∆u. R

This form (1.3.8) of the equations is not only of theoretical interest, shedding light on the role of the pressure, but is of practical interest for numerical algorithms.11 The pressure in compressible flows is conceptually different than in incompressible flows just as it was in ideal flow. If we think of viscous flow as ideal flow with viscous effects added on, it is not unreasonable to assume that p is still a function of ρ. A note of caution is appropriate here. The expressions for p(ρ) used in practical situations are often borrowed from the science of equilibrium thermodynamics. It is not obvious that p as defined here (through equation (1.3.1)) is identical to p as defined in that other science. Not all quantities called p are equal. The use of expressions from equilibrium thermodynamics requires an additional physical justification, which is indeed often available, but which should not be forgotten. According to the analysis given earlier, the pressure p in incompressible flow is determined by the equation of continuity div u = 0. To see why this
11 See, for instance, A. J. Chorin, Math. Comp. 23 [1969], 341-353 for algorithms, and D. Ebin and J. E. Marsden, Ann. of Math. 92 [1970], 102–163 for a theoretical investigation of the projection operator and the use of material coordinates.

1.3 The Navier–Stokes Equations

39

is physically reasonable, consider a compressible flow with p = p(ρ), where p (ρ) > 0. If fluid flows into a given fixed volume V , the density in V will increase, and if p (ρ) > 0, then p in V will also increase. If either the change in ρ is large enough or p (ρ) is large enough, −grad p at the boundary of V will begin to point away from V , and through the term −grad p in the equation for ∂t u, this will cause the fluid to flow away from V . Thus, the pressure controls and moderates the variations in density. If the density is to remain invariant, this must be accomplished by an appropriate p, that is, div u = 0 determines p. In the Navier–Stokes equations for a viscous incompressible fluid, namely, ∂t u + (u · ∇)u = −∇p + we call and (u · ∇)u, the inertia or convective term. The equations say that u is convected subject to pressure forces and, at the same time, is dissipated. Suppose R is very small. If we write the equations 1 in the form ∂t u = P(−u · ∇u + R ∆u), we see that they are approximated by ∂t u = P that is, ∂t u = − grad p + 1 ∆u R 1 ∆u , R and div u = 0, 1 ∆u, R 1 ∆u, R

the diffusion or dissipation term,

which are the Stokes’ equations for incompressible flow. These are linear equations of “parabolic” type. For small R (i.e., slow velocity, large viscosity, or small bodies), the solution of the Stokes equation provides a good approximation to the solution of the Navier–Stokes equations. Later, we shall mostly be interested in flows with the large R; for these the inertial term is important and in some sense is dominant. We hesitate and say “in some sense” because no matter how small (1/R)∆u may be, it still produces a large effect, namely, the change in boundary conditions from u · n = 0 when (1/R)∆u is absent to u = 0 when it is present. There is a major difference between the ideal and viscous flow with regard to the energy of the fluid. The viscous terms provide a mechanism by which macroscopic energy can be converted into internal energy. General principles of thermodynamics state that this energy transfer is one-way. In particular, for incompressible flow, we should have d Ekinetic ≤ 0. dt (1.3.9)

40

1 The Equations of Motion

We calculate (d/dt)Ekinetic for incompressible viscous flow using the transport theorem, as we did in §1.1. We get d d Ekinetic = dt dt =
D 1 2

ρ u 2 dV =
D D

ρu ·

Du dV Dt dV,

1 −u · ∇p + u · ∆u R

by (1.3.3) and div u = 0. Because u is orthogonal to grad p, we get d 1 Ekinetic = dt R u · ∆u dV.
D

The vector identity div(f V) = f div V + V · ∇f gives ∇ · (u∇u + v∇v + w∇w) = ∇u · ∇u + ∇v · ∇v + ∇w · ∇w + u∆u + v∆v + w∆w. This equation, the divergence theorem, and the boundary condition u = 0 on ∂D enable us to simplify the expression for (d/dt)Ekinetic to d Ekinetic = −µ dt ∇u
D 2

dV,

(1.3.10)

where ∇u 2 = ∇u · ∇u = ∇u 2 + ∇v 2 + ∇w 2 . Notice that (1.3.9) and (1.3.10) are compatible exactly when µ ≥ 0 (or, equivalently, ν ≥ 0 or 0 < R ≤ ∞). In other words, there is no such thing as “negative viscosity.” A similar analysis for compressible flow and making use of (1.3.2)’ leads to the inequalities µ ≥ 0 and λ + 2 µ ≥ 0 3 and with σ given by (1.3.2).12 At the end of §1.1 we noted that ideal flow in a channel leads to unreasonable results. We now reconsider this example with viscous effects. Example Consider stationary viscous incompressible flow between two stationary plates located at y = 0 and y = 1, as shown in Figure 1.3.3. We seek a solution for which u(x, y) = (u(x, y), 0) and p is only a function of x, with p1 = p(0), p2 = p(L), and p1 > p2 , so the fluid is “pushed” in the positive x direction. The incompressible Navier–Stokes equations are ∂x u = 0 and 0 = −u ∂x u − ∂x p + 1 2 2 ∂ u + ∂y u R x (incompressibility)

12 See, for example, S. Chapman and T. G. Cowling, The Mathematical Theory of Non-uniform Gases, Cambridge University Press, 1958.

1.3 The Navier–Stokes Equations

41

with boundary conditions u(x, 0) = u(x, 1) = 0. Because ∂x u = 0, u is only a function of y and thus, writing u(x, y) = u(y), we obtain p = 1 u . R

y y =1 flow direction y =0 pressure = p 1 x =0 pressure = p 2 x =L fixed wall

x fixed wall

Figure 1.3.3. Flow between two parallel plates; the fluid is pushed from left to right and correspondingly, p1 > p2 .

Because each side depends on different variables, p = constant, Integration gives p(x) = p1 − and u(y) = y(1 − y)R ∆p x, L ∆p = p1 − p2 , ∆p . 2L 1 u = constant. R

Notice that the velocity profile is a parabola (Figure 1.3.4). The presence of viscosity allows the pressure forces to be balanced by 1 the term R u (y) and allows the fluid to achieve a stationary state. We saw that this was not possible for ideal flow. Next we consider the vorticity equation for (homogeneous) viscous incompressible flow. In the two-dimensional case we proved in §1.2 (see equation (1.2.12)) that for isentropic ideal plane flow, Dξ/Dt = 0. The derivation is readily modified to cover viscous incompressible flow; the result is Dξ 1 = ∆ξ. (1.3.11) Dt R This shows that the vorticity is diffused by viscosity as well as being tranported by the flow. Introduce the stream function ψ(x, y, t) by means of

42

1 The Equations of Motion

y y=1

u(y) x

Figure 1.3.4. Viscous flow between two plates.

(1.2.15)2 and (1.2.15)3 as before. We saw that we could impose the boundary condition ψ = 0 on ∂D. Now, however, the no-slip condition u = 0 on ∂D implies that ∂x ψ = 0 = ∂y ψ on ∂D by (1.2.15)3 . Because ψ = 0 on ∂D implies that the tangential derivative of ψ on ∂D vanishes, we get the extra boundary condition ∂ψ =0 ∂n on ∂D

This extra condition should be somewhat mystifying; certainly we cannot impose it when we solve ∆ψ = −ξ, ψ = 0 on ∂D, because this problem already has a solution. Thus, it is not clear how to get the system  Dξ 1   = ∆ξ,   Dt R  (1.3.12) ∆ψ = −ξ, ψ = 0 on ∂D,      u = ∂y ψ, v = −∂x ψ to work. We shall study this problem in §2.2. For three-dimensional viscous incompressible flow, the vorticity equation is Dξ 1 − (ξ · ∇)u = ∆ξ. (1.3.13) Dt R Thus, vorticity is convected, stretched, and diffused. (The left-hand side of (1.3.13) is called the Lie derivative. It is this combination, rather than each term separately, that makes coordinate independent sense.) Here the problems with getting a system like (1.3.12) are even worse; even in the isentropic case we had trouble with (1.2.16) because of boundary conditions. For viscous flow, circulation is no longer a constant of the motion. One might suspect from (1.3.13) that if ξ = 0 at t = 0, then ξ = 0 for all time. However, this is not true: viscous flow allows for the generation of vorticity.

1.3 The Navier–Stokes Equations

43

This is possible because of the difference in boundary conditions between ideal and viscous flows. The mechanism of vorticity generation is related to the difficulties with the boundary conditions in equations (1.3.12) and will be discussed in §2.2. For many of our discussions we have made the assumption of incompressibility. We now give a heuristic analysis of when such an assumption will be reasonable and when, instead, the compressible equations should be used. We shall do this in the context of isentropic stationary flows for simplicity. Assume that we have an equation of state p = p(ρ), Define c= p (ρ). For reasons that will become clear later, c is called the sound speed of the fluid. Thus, we have c2 dρ = dp. (1.3.14) Let u = u be the flow speed. One calls M = u/c the (local) Mach number of the flow; it is a function of position in the flow. From Bernoulli’s theorem proved in §1.1, u2 + 2 dp = constant on streamlines. ρ(p) (1.3.15) p (ρ) > 0.

Also, differentiating the continuity equation in the form (1.2.10) along streamlines gives 0 = Jdρ + ρ dJ, (1.3.16) where J is the Jacobian of the flow map. Combining (1.3.14), (1.3.15), and (1.3.16) we get dJ du = −M . J c The flow will be approximately incompressible if J changes only by a small amount along streamlines. Thus, a steady flow can be viewed as incompressible if the flow speed is much less than the sound speed, u c, i.e., M 1,

or if changes in the speed along streamlines are very small compared to the sound speed. For example, for equations of state of the kind associated with ideal gases, p = Aργ , γ > 1,

44

1 The Equations of Motion

we have c= γp ρ

so the flow will be approximately incompressible if γ is very large. For nonsteady flow one also needs to know that l τ c,

where l is a characteristic length and τ is a characteristic time over which the flow picture changes appreciably.13 The presence of viscosity does not alter these conclusions significantly.

Exercises
Exercise 1.3-1 Find a stationary viscous incompressible flow in a circular pipe with radius a > 0 and with pressure gradient ∇p. Exercise 1.3-2 Show that the incompressible Navier–Stokes equations in cylindrical coordinates are (i) ρ Dvr v2 − θ Dt r Dvθ vr vθ + Dt r = ρfr − ∂p vr 2 ∂vθ + µ ∆vr − 2 − 2 ∂r r r ∂θ vθ 1∂p 2∂v2 + µ ∆vθ + 2 − 2 r∂θ r ∂θ r .

(ii) ρ

= ρfθ −

.

(iii) ρ

Dvz ∂p = ρfz − + µ∆vz , Dt ∂z ∆= 1 ∂ r ∂r r ∂ ∂r + ∂2 1 ∂2 + 2 r2 ∂θ2 ∂z

where and Exercise 1.3-3

∂ ∂ D ∂ vθ ∂ = + vr + + vz . Dt ∂t ∂r r ∂θ ∂z Flow in an infinite pipe.

(i) Poiseuille flow . Work in cylindrical coordinates with a pipe of radius a aligned along the z-axis. The no-slip boundary condition is
13 Theoretical work on the limit c → ∞ is given by D. Ebin, Ann. Math. 141 [1977], 105, and S. Klainerman and A. Majda, Comm. Pure Appl. Math., 35 [1982], 629. Algorithms for solving the equations for incompressible flow by exploiting the regularity of the limit c → ∞ can be found in A. J. Chorin, J. Comp. Phys.12 [1967], 1.

1.3 The Navier–Stokes Equations

45

v = 0 when r = a. Assume the solution takes the form p = Cz, C constant, vz = vz (r), and vr = vθ = 0. Using Exercise 1.3-2, obtain C = µ∆vz = µ Integration yields vz = − C 2 r + A log r + B, 4µ 1 ∂ r ∂r r ∂vz ∂r .

where A, B are constants. Because we require that the solution be bounded, A must be 0, because log r → −∞ as r → 0. Use the no-slip condition to determine B and obtain vz = C 2 (a − r2 ). 4µ

(ii) Show that the mass flow rate Q = s ρvz dA through the pipe is Q = ρπCa4 /8µ. This is the so-called fourth-power law. (iii) Determine the pressure on the walls. Exercise 1.3-4 Compute the solution to the problem of stationary viscous flow between two concentric cylinders and determine the pressure on the walls. (Hint: Proceed as above, but retain the log term.)

46

1 The Equations of Motion

This is page 47 Printer: Opaque this

2
Potential Flow and Slightly Viscous Flow

The goal of this chapter is to present a deeper study of the relationship between viscous and nonviscous flows. We begin with a more detailed study of inviscid irrotational flows, that is, potential flows. Then we go on to study boundary layers, where the main difference between slightly viscous and inviscid flows originates. This is further developed in the third section using probabilistic methods. For most of this chapter we will study incompressible flows. A detailed study of some special compressible flows is the subject of Chapter 3.

2.1

Potential Flow

Throughout this section, all flows are understood to be ideal (i.e., inviscid); in other words, either incompressible and nonviscous or isentropic and nonviscous. Although we allow both, our main emphasis will be on the incompressible case. A flow with zero vorticity is called irrotational . For ideal flow, this holds for all time if it holds at one time by the results of §1.2. An inviscid, irrotational flow is called a potential flow . A domain D is called simply connected if any continuous closed curve in D can be continuously shrunk to a point without leaving D. For example, in space, the exterior of a solid sphere is simply connected, whereas in the plane the exterior of a solid disc is not simply connected.

48

2 Potential Flow and Slightly Viscous Flow

For irrotational flow in a simply connected region D, there is a scalar function ϕ(x, t) on D called the velocity potential such that for each t, u = grad ϕ. In particular, it follows that the circulation around any closed curve C in D is zero. Using the identity (u · ∇)u = 1 ∇ 2 u
2

− u × (∇ × u),

(2.1.1)

we can write the equations for isentropic potential flow in the form ∂t u + 1 ∇( u 2 ) = − grad w, 2 where w is the enthalpy, as in §1.1. Substituting u = grad ϕ, we obtain grad ∂t ϕ + and thus ∂t ϕ +
1 2 1 2 1 2

u

2

+ w = 0, (2.1.2)

u

2

+ w = constant in space.

In particular, if the flow is stationary, u
2

+ w = constant in space.

For the last equation to hold, simple connectivity of D is unnecessary. The version of Bernoulli’s theorem given in §1.1 concluded that 1 u 2 + w 2 was constant on streamlines. The stronger conclusion here results from the additional irrotational hypothesis, ξ = 0. For homogeneous incompressible ideal flow, note that w = p/ρ0 from the definition of w. For potential flow in nonsimply connected domains, it can occur that the circulation ΓC around a closed curve C is nonzero. For instance, consider u= −y x , x2 + y 2 x2 + y 2

outside the origin. If the contour C can be deformed within D to another contour C , then ΓC = ΓC ; see Figure 2.1.1. The reason is that basically C ∪ C forms the boundary of a surface Σ in D. Stokes’ theorem then gives ξ · dA =
Σ C

u · ds −
C

u · ds = ΓC − ΓC

and because ξ = 0 in D, we get ΓC = ΓC . (A more careful argument proving the invariance of ΓC under deformation is given in books on complex variables.) Notice that from §1.2, the circulation around a curve is constant in time. Thus, the circulation around an obstacle in the plane is well-defined and is constant in time. Next, consider incompressible potential flow in a simply connected domain D. From u = grad ϕ and div u = 0, we have ∆ϕ = 0.

2.1 Potential Flow

49

C C' Σ

Figure 2.1.1. The circulations about C and C are equal if the flow is potential in Σ.

Let the velocity of ∂D be specified as V, so u · n = V · n. Thus, ϕ solves the Neumann problem: ∆ϕ = 0, ∂ϕ = V · n. ∂n (2.1.3)

If ϕ is a solution, then u = grad ϕ is a solution of the stationary homogeneous Euler equations, i.e., ρ(u · ∇)u = − grad p, div u = 0, u · n = V · n on ∂D,

(2.1.4)

where p = −ρ u 2 /2. This follows from the identity (2.1.1). Therefore, solutions of (2.1.3) are in one-to-one correspondence with irrotational solutions of (2.1.4) (with ϕ determined only up to an additive constant) on simply connected regions. This observation leads to the following. Theorem Let D be a simply connected, bounded region with prescribed velocity V on ∂D. Then i ii there is exactly one potential homogeneous incompressible flow (satisfying (2.1.4)) in D if and only if ∂D V · n dA = 0; this flow is the minimizer of the kinetic energy function Ekinetic = 1 2 ρ u
D 2

dV,

among all divergence-free vector fields u on D satisfying u ·n = V·n.

50

2 Potential Flow and Slightly Viscous Flow

Proof i The Neumann problem (2.1.3) has a solution if and only if the obvious necessary condition ∂D V · n dA = 0 holds, as was mentioned earlier. We can demonstrate the uniqueness of u directly as follows: Let u and u be two solutions, and let v = u − u , ψ = ϕ − ϕ . Then ∆ψ = 0, ∂ψ/∂n = 0, and v = grad ψ. Hence, div(ψv) dV =
D D

v · grad ψ dV +
D

ψ div v dV =
D

v · v dV.

On the other hand, div(ψv) dV =
D ∂D

ψv · n dA = 0.

Thus, ii

D

v 2 dV = 0 and v = 0, that is, u = u .

Let u solve (2.1.4) and let u be divergence free and u = n = V · n. Let v = u − u ; then div v = 0 and v · n = 0 on ∂D. Therefore, Ekinetic − Ekinetic =
1 2

ρ( u
D

2

− u
2

2

) dV ρ(u − u ) · u dV
D

= −1 2 ≤
D

ρ u−u
D

dV +

ρv · grad ϕ dV = 0.

The last equality follows by the orthogonality relation proved in §1.3. Thus, Ekinetic ≤ Ekinetic as claimed.

Notice, in particular, that the only incompressible potential flow in a bounded region with fixed boundary is the trivial flow u = 0. For unbounded regions this is not true without a careful specification of what can happen at infinity; the above uniqueness proof is valid only if the integration by parts (i.e., use of the divergence theorem) can be justified. For example, in polar coordinates in the plane, ϕ(r, θ) = r+ 1 r cos θ

solves (2.1.3) with ∂ϕ/∂n = 0 on the unit circle and on the x-axis. It represents a nontrivial irrotational potential flow on the simply connected

2.1 Potential Flow

51

y D streamlines

x

Figure 2.1.2. Potential flow in the upper half-plane outside the unit circle.

region D shown in Figure 2.1.2. This flow may be arrived at by the methods of complex variables to which we will now turn. Incompressible potential flow is very special, but is a key building block for understanding complicated flows. For plane flows the methods of complex variables are useful tools. Let D be a region in the plane and suppose u = (u, v) is incompressible, that is, ∂u ∂v + =0 (2.1.5) ∂x ∂y and is irrotational, that is, ∂u ∂v − = 0. ∂y ∂x Let F = u − iv, (2.1.7) which is called the complex velocity. Equations (2.1.5) and (2.1.6) are exactly the Cauchy-Riemann equations for F , and so F is an analytic function on D. Conversely, given any analytic function F , u = Re F and v = −Im F define an incompressible (stationary) potential flow. If F has a primitive, F = dW/dz, then we call W the complex potential . (If one allows multivalued functions, W will always exist, but such a convention could cause confusion.) Write W = ϕ + iψ. Then (2.1.7) is equivalent to u = ∂x ϕ = ∂y ψ and v = ∂y ϕ = −∂x ψ, (2.1.6)

that is, u = grad ϕ and ψ is the stream function. In what follows, however, we do not and must not assume a (single-valued) W exists. Consider a flow in the exterior of an obstacle B (Figure 2.1.3).

52

2 Potential Flow and Slightly Viscous Flow

B D n

Figure 2.1.3. Flow around an obstacle.

The force on the body equals the force exerted on ∂B by the pressure, that is, F =−
∂B

pn ds,

(2.1.8)

which means that for any fixed vector a, F ·a=−
∂B

pn · a ds.

Formula (2.1.8) was already discussed at length in §1.1. We next prove a theorem that gives a convenient expression for F. Blasius’ Theorem For incompressible potential flow exterior to a body B (with rigid boundary) and complex velocity F , the force F on the body is given by −iρ F= F 2 dz (2.1.9) 2 ∂B where the overbar denotes complex conjugation and where the vector F is identified with a complex number in the standard way; i.e., (x, y) is identified with z = x + iy. Proof If dz = dx + i dy represents an infinitesimal displacement along the boundary curve C = ∂B, then (1/i)dz = dy − i dx represents a normal displacement. Thus, by (2.1.8) F =−
C

p dy + i
C

p dx = i
C

p(dx + i dy).

As we observed in (2.1.4), p= −ρ(u2 + v 2 ) , 2 and therefore F = −iρ 2 (u2 + v 2 ) dz.
C

2.1 Potential Flow

53

On the other hand, F 2 = (u − iv)2 = u2 − v 2 − 2iuv, and because u is parallel to the boundary, we get u dy = v dx. Thus, F 2 dz = (u2 − v 2 − 2iuv)(dx + i dy) = (u2 + v 2 )(dx − i dy), and because u2 + v 2 is real, F 2 dz = (u2 + v 2 ) dz. This formula will be used to prove the following (Figure 2.1.4):
F U U B

U

Figure 2.1.4. The Kutta–Joukowski theorem gives the force exerted on B.

Kutta–Joukowski Theorem Consider incompressible potential flow exterior to a region B. Let the velocity field approach the constant value (U, V ) = U at infinity. Then the force exerted on B is given by F = −ρΓC U n, (2.1.10)

where ΓC is the circulation around B and n is a unit vector orthogonal to U. Proof By assumption, the complex velocity F is an analytic function outside the body B. It may therefore, be expanded in a Laurent series. Because F tends to a constant U at infinity, no positive powers of z can occur in the expansion. In other words, F has the form a1 a3 a2 F = a0 + + 2 + 3 + ··· z z z valid outside any disc centered at the origin and containing B. Because U is the velocity at infinity, a0 = U − iV . By Cauchy’s theorem, F dz = 2πa1 i,
C

where C = ∂B. (The integral is unchanged if we change C to a circle of large radius.) However, F dz =
C C

(u − iv)(dx + i dy) =
C

u dx + v dy =
C

u · ds = ΓC

54

2 Potential Flow and Slightly Viscous Flow

because u dy = v dx, i.e., u is parallel to ∂B. Therefore, a1 = ΓC . 2πi

Squaring F gives the Laurent expansion F 2 = a2 + 0 2a0 a2 + a2 2a0 a1 1 + + ··· . z z2

By Blasius’ theorem and Cauchy’s theorem, F =− iρ 2 F 2 dz = −
C

iρ · (2πi · 2a0 a1 ) = ρΓC (V − iU ) 2

which proves the theorem. Notice that the force exerted on the body B by the flow is normal to the direction of flow and is proportional to the circulation around the body. In any case, the body experiences no drag (i.e., no force opposing the flow) in contradiction with our intuition and with experiment. The difficulty, of course, stems from the fact that we have neglected viscosity. (We shall remedy this in the succeeding two sections.) Even worse, if ΓC = 0, there is no net force on the body at all, a fact hard to reconcile with our intuition even for ideal flow. This result is called d’Alembert’s paradox. Example 1 For a complex number α = U − iV , let W (z) = αz. Thus, F (x) = α, so the velocity field is u = (U, V ). This is two-dimensional flow moving with constant velocity in the direction (U, V ). Example 2 Let B be the disc of radius a > 0 centered at the origin in the complex plane, and let W (z) = U z+ a2 z (2.1.11)

for a positive constant U . The complex velocity is F (z) = W (z) = U 1− a2 z2 , (2.1.12)

which approaches U at ∞. The velocity potential ϕ and the stream function ψ are determined by W = ϕ + iψ . To verify that the flow is tangent to the circle |z| = a, we need only to show that ψ = constant when |z| = a. In fact, for |z|2 = z z = a2 , we have from (2.1.1), ¯ W (z) = U (z + z ), ¯ so W is real on |z| = a, that is, ψ = 0 on |z| = a. The flow is shown in Figure 2.1.5.

2.1 Potential Flow

55

y ψ = constant ϕ = constant velocity = u B A B D C x

Figure 2.1.5. Potential flow around a disc.

From (2.1.12) with z = aeiθ , that is, z on ∂B, we find F (z) = U 1− a2 a2 e2iθ = U (1 − cos 2θ + i sin 2θ).

Thus, the velocity is zero at A and C; that is, A and C are stagnation points. The velocity reaches a maximum at B and D. By Bernoulli’s theorem, we can write ρ p = − u 2 + constant; 2 thus, the pressure at A and C is maximum and is a minimum at B and D. The disc has zero circulation because F = W and W is single-valued. If W is any analytic function defined in the whole plane, then ˜ W (z) = W (z) + W a2 , z |z| ≥ a

is a potential describing a flow exterior to the disc of radius a > 0, but possibly with a more complicated behavior at infinity. This is proved along the same lines as in the argument just presented. Example 3 In §1.2 we proved that choosing ψ to be an arbitrary increasing function of r alone yields a flow that is incompressible and has vorticity ξ = −∆ψ. If we can arrange for ψ to be the imaginary part of an analytic function, then the flow will be irrotational as well. The function W (z) = Γ log z 2πi (2.1.13)

56

2 Potential Flow and Slightly Viscous Flow

has this property, because log z = log |z| + i arg z. Of course, W (z) is not single-valued, but the complex velocity F (z) = Γ 2πiz (2.1.14)

is analytic and single-valued outside z = 0. The circulation is indeed Γ. Note that the velocity field is zero at infinity. For incompressible potential flow about a disc of radius a centered at z0 , we need only choose W (z) = Γ log(z − z0 ). 2πi

The boundary conditions are satisfied because ψ is constant on any circle centered at z0 (see Figure 2.1.6). The incompressible potential flow with W (z) = Γ log(z − z0 ) 2πi

will be called a potential vortex at z0 .

ϕ = constant

a

B z0

ψ = constant

Figure 2.1.6. Potential vortex flow centered at z0 .

Example 4

We combine Examples 2 and 3 by forming W (z) = U z+ a2 z + Γ log z, 2πi (2.1.15)

where |z| ≥ a. Because ψ is constant on the boundary for each flow separately, it is also true for W given by (2.1.15). Thus, we get an incompressible potential flow on the exterior of the disc |z| ≤ a with circulation Γ around the disc. The velocity field is (U, 0) at infinity (therefore, the

2.1 Potential Flow

57

Kutta–Joukowski theorem applies). On the surface of the disc the velocity u = grad ϕ is tangent to the disc and is given in magnitude by velocity = Here ϕ = Re W , so that ϕ(r, θ) = U cos θ r + and thus velocity = a2 r + Γθ , 2π Γ . 2πa 1 ∂ϕ r ∂θ . r=a 1 ∂ϕ a ∂θ

= −2U sin θ + r=a If |Γ| < 4πa U , there are two stagnation points A and C defined by sin θ = Γ 4πaU

on the boundary, where the pressure is highest. See Figure 2.1.7.

Jerry’s Books

B Γ U

A

D

C

Figure 2.1.7. Flow around a disc with circulation.

This example helps to explain the Kutta–Joukowski theorem; note that the vertical lift may be attributed to the higher pressure at A and C. The symmetry in the y-axis means that there is no drag. D’Alembert’s Paradox in Three Dimensions In the case of steady incompressible potential flow around an obstacle in three dimensions with constant velocity U at infinity, not only can there be no drag, there can be no lift either.

58

2 Potential Flow and Slightly Viscous Flow

The difference with the two-dimensional case is explained by the fact that the exterior of a body in three-space is simply connected, whereas this is not true in two dimensions. We will not present the detailed proof of d’Alembert’s paradox here, but we can give the idea. Recall that the solution of ∆ϕ = −ρ in space is ϕ(x) = 1 4π ρ(y) dV (y), x−y

that is, ϕ is the potential due to a charge distribution ρ. Notice that if ρ is concentrated in a finite region, then ϕ(x) = O where r = x , that is, constant r for large r. In fact, as we know physically, ϕ(x) ≈ Q/4πr for r large, where Q = ρ(y) dV (y) is the total charge. If Q = 0, then ϕ(x) = O(1/r2 ) because the first term in the expansion in powers of 1/r is now missing. For an incompressible potential flow there will be a potential ϕ, that is, u = grad ϕ (because the exterior of the body is simply connected). The potential satisfies ∆ϕ = 0, ∇ϕ = U at ∞, |ϕ(x)| ≤ ∂ϕ = 0 on the boundary of the obstacle. ∂n The solution here can then be shown to satisfy ϕ(x) = U · x + O 1 r and 1 r ,

as in the potential case above. However, there is an integral condition analogous to Q = 0, namely, the net outflow at ∞ should be zero. This means ϕ(x) = U · x + O Hence, 1 r2 .

u = U + O(r−3 ).

(2.1.16)

Because p = −ρv 2 /2, we also have p = p0 + O(r−3 ). (To see that this is true, write u 2 = U 2 + (u − U) · (u + U)). The force on the body B is F =−
∂B

pn dA.

2.1 Potential Flow

59

Let Σ be a surface containing B. Because u · n = 0 on ∂B and the flow is steady, equation (BM3) from §1.1 applied to the region between B and Σ shows that F = − (ρ(u · n)u + pn) dA.
Σ

We are free to choose Σ to be a sphere of large radius R enclosing the obstacle. Then F =−
Σ

(p0 n + ρ(U · n)U) dA + (area Σ) · O(R−3 )

= 0 + O(R−1 ) → 0 as R → ∞. Hence, F = 0. One may verify d’Alembert’s paradox directly for flow past a sphere of radius a > 0. In this case ϕ=− where n = x x

a3 U · n + x · U, 2r2

, and u=− a3 [3n(U · n) − U] + U, 2r2

where U is the velocity at infinity. We leave the detailed verification to the reader.1 Next we will discuss a possible mechanism, ultimately to be justified by the presence of viscosity, by which one can avoid d’Alembert’s paradox. An effort to resolve the paradox is of course prompted by the fact that real bodies in fluids do experience drag. By an almost potential flow , we mean a flow in which vorticity is concentrated in some thin layers of fluid; the flow is potential outside these thin layers, but there is a mechanism for producing vorticity near boundaries. For example, one can postulate that the flow past the obstacle shown in Figure 2.1.8 produces an almost potential flow with vorticity produced at the boundary and concentrated on two streamlines emanating from the body. We image different potential flows in the two regions separated by these streamlines with the velocity field discontinuous across them. For such a model, the Kutta–Joukowski theorem does not apply and the drag may be different from zero. There are a number of situations in engineering where real flows can be usefully idealized as “nearly potential.” These situations arise in particular when one considers “streamlined” bodies, that is, bodies
1 See L. Landau and E. Lifschitz [1959] Fluid Mechanics, Pergamon, p. 34 for more information.

60

2 Potential Flow and Slightly Viscous Flow

potential flow

potential flow

ξ≠0 ξ≠0

Figure 2.1.8. Almost potential flow has vorticity concentrated on two curves.

so shaped as to reduce their drag. The discussion of such bodies, their design, and the validity of potential approximation to the flow around them are outside the scope of this book. Next we shall examine a model of incompressible inviscid flow inspired by the idea of an almost potential flow and Example 3. We imagine the vorticity in a fluid is concentrated in N vortices (i.e., points at which the vorticity field is singular), located at x1 , x2 , . . . , xN in the plane (Figure 2.1.9). The stream function of the jth vortex, ignoring the other vortices for a moment, is by Example 3, ψj (x) = − Γj log x − xj . 2π (2.1.17)

x3

x4

x2

x1

Figure 2.1.9. The flow generated by point vortices in the plane.

As the fluid moves according to Euler’s equations, the circulations Γj associated with each vortex will remain constant. The vorticity field produced by the jth vortex can be written as ξj = −∆ψj = Γj δ(x − xj ), where δ is the Dirac δ function. This equation arises from the fact, which we just accept, that the Green’s function for the Laplacian in the plane

2.1 Potential Flow

61

is G(x, x ) = that is, G satisfies

1 log x − x , 2π

∆x G(x, x ) = δ(x − x ). The solution of ∆ψ = −ξ is then given by the superposition ψ(x) = − ξ(x )G(x − x ) dx
N j=1

which in our case reduces to ψ(x) = ψj (x) = −

ψj (x), where

1 Γj log x − xj . 2π

The velocity field induced by the jth vortex (again ignoring the other vortices) is given by uj = (∂y ψj , −∂x ψj ) = − Γj 2π y − yj r2 , Γj 2π x − xj r2 , (2.1.18)

where r = x − xj . Let the vortices all move according to the velocity field
N

u(x, t) = j=1 uj (x, t),

where uj is given by (2.1.18) except we now allow, as we must, the centers of the vortices xj , j = 1, . . . , N to move. Each one ought to move as if convected by the net velocity field of the other vortices. Therefore, by (2.1.18), xj moves according to the equations dxj 1 =− dt 2π Γi (yj − yi ) 2 rij and 1 dyj = dt 2π Γi (xj − xi ) , 2 rij (2.1.19)

i=j

i=j

where rij = xi − xj . Let us summarize the construction of the flows we are considering: choose constants Γ1 , . . . , ΓN and initial points x1 = (x1 , y1 ), . . . , xN = (xN , yN ) in the plane. Let these points move according to the equations (2.1.19). Define uj by (2.1.18) and let
N

u(x, t) = − j=1 uj (x, t).

This construction produces formal solutions of Euler’s equation (“formal” because the meaning of δ-function solutions of nonlinear equations is not obvious). These solutions have the property that the circulation theorem

62

2 Potential Flow and Slightly Viscous Flow

is satisfied. If C is a contour containing l vortices at x1 , x2 , . . . , xl , then l ΓC = − i=1 Γi and ΓC is invariant under the flow. However, the relationship between these solutions and bona fide solutions of Euler’s equations is not readily apparent. Such a relationship can, however, be established rigorously and such vortex systems do contain significant information about the solutions of Euler’s equations under a wide variety of conditions.2 An important property of the equations is that they form a Hamiltonian system. Define H=− 1 4π Γi Γj log xi − xj . i=j (2.1.20)

First of all, it is easy to check that (2.1.19) is equivalent to Γj dxj ∂H , = dt ∂yj Γj dyj ∂H , =− dt ∂xj (2.1.21)

where j = 1, . . . , N (there is no sum on j). Introduce the new variables xi = |Γi |xi , yi = |Γi | sgn(Γi )yi , i = 1, . . . , N,

where sgn(Γi ) is 1 if Γi > 0, and is −1 otherwise. Then (2.1.19) is equivalent to the following system of Hamiltonian equations dxi ∂H , = dt ∂yi dyi ∂H , =− dt ∂xi i = 1, . . . , N, (2.1.22)

with Hamiltonian H and generalized coordinates (xi , yi ). As in elementary mechanics, dH = dt
N

i=1

∂H dxi ∂H dyi + = 0, ∂xi dt ∂yi dt i=1

N

that is, H is a constant of the motion. A consequence of this fact is that if the vortices all have the same sign they cannot collide. If xi −xj = 0, i = j at t = 0, then this remains so for all time because if xi − xj → 0, H becomes infinite. This Hamiltonian system is of importance in understanding how vorticity evolves and organizes itself.3
2 See C. Anderson and C. Greengard, On Vortex Methods, SIAM J. Sci. Statist. Comput. 22 [1985], 413. 3 The Euler equations themselves form a Hamiltonian system (this is explained, along with references, in R. Abraham and J. E. Marsden, Foundations of Mechanics, 2nd Edition [1978]), and the Hamiltonian nature of the vortex approximation is consistent with this. See also J. E. Marsden and A. Weinstein, Coadjoint orbits, vortices and Clebsch variables for incompressible fluids, Physica 7D [1983], 305–323.

2.1 Potential Flow

63

Let us now generalize the situation and imagine our N vortices moving in a domain D with boundary ∂D. We can go through the same construction as before, but we have to modify the flow uj of the jth vortex in such a way that u · n = 0, that is, the boundary conditions appropriate to the Euler equations hold. We can arrange this by adding a potential flow vj to uj such that uj · n = −vj · n on ∂D. In other words, we choose the stream function ψj for the jth vortex to solve ∆ψj = −ξj = −Γj δ(x − xj ) with ∂ψj = 0 on ∂D. ∂n

This is equivalent to requiring ψj (x) = −Γj G(x, xj ) where G is the Green’s function for the Neumann problem for the region D. This procedure will appropriately modify the function (1/2π) log x−xj and allow the analysis to go through as before. For example, suppose D is the upper half-plane y ≥ 0. Then we get G by the reflection principle: G(x, xj ) = 1 ˆ (log x − xj + log x − xj ) , 2π

ˆ where xj = (xj , −yj ) is the reflection of xj across the x-axis (see Figure 2.1.10). For the Neumann-Green’s functions for other regions the reader may consult textbooks on partial differential equations. y (x,y )

(xi ,yi ) Γ

motion of the vortex

D :y >0

x –Γ

(xi ,–yi )

Figure 2.1.10. The stream function at (x, y) is the superposition of those due to vortices with opposite circulations located at (xi , yi ) and (xi , −yi ).

64

2 Potential Flow and Slightly Viscous Flow

Consider again Euler’s equations in the form ∆ψ = −ξ, One can write ψ=− ξ(x )G(x, x ) dx , u = ∂y ψ, v = −∂x ψ, Dξ = 0. Dt

where G(x, x ) = 1 π log x − x , and set u = ∂y ψ, v = −∂x ψ. The re2 sulting equations resemble the equations just derived for a system of point vortices. The integral for ψ here resembles the formula for ψ in the point vortex system somewhat as an integral resembles one of its Riemann sum approximations. This suggests that an ideal incompressible flow can be approximated by the motion of a system of point vortices. There are in fact theorems along these lines.4 Vortex systems provide both a useful heuristic tool in the analysis of the general properties of the solutions of Euler’s equations, and a useful starting point for the construction of practical algorithms for solving these equations in specific situations. One can ask if there is a similar construction in three dimensions. First of all, one can seek an analogue of the superposition of stream functions from point potential vortices. Given u satisfying div u = 0, there is a vector field A such that div A = 0 and such that u = curl A, and therefore ∆A = −ξ. In three dimensions, Green’s function for the Laplacian is given by G(x, x ) = − 1 1 , 4π x − x x=x.

Then we can represent A in terms of ξ by A=− 1 4π ξ(x ) dV (x ), s

where s = x − x , and where dV (x ) is the usual volume element in space. It is easy to check that A defined by the above integral satisfies the normalization condition div A = 0. Thus, because u = curl A, we obtain u(x) = 1 4π s×ξ dV (x ), s3

where s = x − x and ξ = ξ(x ). Suppose that we have a vortex line C in space with circulation Γ (see Figure 2.1.11) and we assume that the
4 The discrete vortex method is discussed in L. Onsager, Nuovo Cimento 6 (Suppl.) [1949], 229; A. J. Chorin, J. Fluid Mech. 57 [1973], 781; and A. J. Chorin, SIAM J. Sci. Statist. Comput. 1 [1980], 1. Convergence of solutions of the discrete vortex equations to solutions of Euler’s equations as N → ∞ is discussed in O. H. Hald, SIAM J. Numer. Anal. 16 [1979], 726; T. Beale and A. Majda, Math. Comp. 39 [1982], 1–28, 29–52; and K. Gustafson and J. Sethian, Vortex Flows, SIAM Publications, 1991.

2.1 Potential Flow

65

vorticity field ξ is concentrated on C only, that is, the flow is potential outside the filament C. Then u(x) can be written as u(x) = 1 4π s × Γ ds s3

C

where ds is the line element on C.

Figure 2.1.11. The flow induced by a vortex filament.

Exercises
Exercise 2.1-1 If f : D → D is a conformal transformation (an analytic function that is one to one and onto), it can be used to transform one complex potential to another. Use f (z) = z +a2 /z (which takes the exterior of the disc of radius a in the upper half-plane to the upper half-plane) and the complex potential in the upper half-plane to generate formula (2.1.11). Exercise 2.1-2 Let F (z) = z 2 be a complex potential in the first quadrant. Sketch some streamlines and the curves φ = constant, ψ = constant, where F = φ + iψ . What is the force exerted on the walls? Exercise 2.1-3 Use conformal maps to find a formula for potential flow over the plate in Figure 2.1.12. What is the force exerted on this plate? Exercise 2.1-4 Let a spherical object move through a fluid in R3 . For slow velocities, assume Stokes’ equations apply. Take the point of view that the object is stationary and the fluid streams by. The setup for the boundary value problem is as follows: given U = (U, 0, 0), U constant, find u and p such that Stokes’ equation holds in the region exterior to a sphere of radius R, u = 0 on the boundary of the sphere and u = U at infinity. The solution to this problem (in spherical coordinates centered in the object)

66

2 Potential Flow and Slightly Viscous Flow

y

(0,L)

x

Figure 2.1.12. Flow over a vertical plate.

is called Stokes’ Flow: 3 U + n(U · n) 1 3 U − 3n(U · n) u=− R + U, − R 4 r 4 r3 3 U·n p = p0 − ν 2 R, 2 r where p0 is constant and n = r/r . (a) Verify this solution. (b) Show that the drag is 6πRνU and there is no lift. (c) Show there is net outflow at infinity (an infinite wake). Exercise 2.1-5 Because of the difficulties in Exercise 2.1-4, Oseen in 1910 suggested that Stokes’ equations be replaced by −ν∆u + (U · ∇)u = − 1 grad p, ρ

(2.1.23)

with div u = 0, where u represents the true velocity minus U. This amounts to linearizing the Navier–Stokes equations about U, whereas Stokes’ equations may be viewed as a linearization about 0. One would thus conjecture that Oseen’s equations are good where the flow is close to the free stream velocity U (away from the sphere) and that Stokes’ equations are good where the velocity is 0 (near the sphere). The solution of Oseen’s equations in the region exterior to a sphere in R3 can be found in Lamb’s book. Show that drag on the sphere for the Oseen solution is F = 6πRU ν(1+ 3 R), 8 where R = U R/ν is the Reynolds number. Thus, there is a difference of the order R in the Stokes and Oseen drag forces. Notes on Exercise 2.1-4 and Exercise 2.1-5: If D is bounded with smooth boundary, then there exists at most one solution to Stokes’ equations. See Ladyzhenskaya’s book listed in the Preface. In the exterior of a bounded region in R3 there exists a unique solution to Stokes’ equations. The situation in R2 is different; in fact, we have the following strange situation:

2.2 Boundary Layers

67

Stokes’ Paradox There is no solution to Stokes’ equations in R2 in the region exterior to a disc (with reasonable boundary conditions).5 Stokes’ paradox does not apply to the Oseen or Navier–Stokes equations in R2 or R3 . However, Filon in 1927 pointed out that for other reasons Oseen’s equations also lead to unacceptable results. The example he gives is a skewed ellipse in a free stream. Computation of the moment exerted on the ellipse reveals that it is infinite! This is not so surprising in view of the fact that Oseen’s equations represent linearization about the free stream. One cannot expect them to give good results around the obstacle because the equations contain errors there of order U 2 .

2.2

Boundary Layers
1 ∂t u + (u · )u = − grad p + R div u = 0, u = 0 on ∂D,  u,     (2.2.1)

Consider the Navier–Stokes equations

and assume the Reynolds number R is large. We ask how different a flow governed by (2.2.1) is from one governed by the Euler equations for incompressible ideal flow:  ∂t u + (u · )u = − grad p,   div u = 0, (2.2.2)   u · n = 0 on ∂D. Imagine that both flows coincide at t = 0 and, say, are irrotational, that is, ξ = 0. Thus, under (2.2.2) the flow stays irrotational, and thus is a potential flow. However, we claim that the presence of the (small) viscosity term (1/R) u and the difference in the boundary conditions have the following effects: 1. The flow governed by (2.2.2) is drastically√ modified near the wall in a region with thickness proportional to 1/ R. 2. The region in which the flow is modified may separate from the boundary. 3. This separation will act as a source of vorticity.
5 See Birkhoff’s book, and J. Heywood, Arch. Rational Mech. Anal. 37 [1970], 48–60, and Acta Math. 129 [1972], 11–34.

68

2 Potential Flow and Slightly Viscous Flow

Whereas the region referred to in 1 gets arbitrarily small as R increases, its effects via 2 and 3 may persist and not diminish in the limit R → ∞. A model for the first effect is given in the following: Example 1 Consider the equation dy = a, dx y(1) = 1, (2.2.3)

where a is a constant and x ranges between x = 0 and x = 1. The solution of (2.2.3) is clearly y = a(x − 1) + 1. Now add a “viscosity” term and a boundary condition to (2.2.3) as follows: d2 y dy + = a, with y(0) = 0, y(1) = 1. (2.2.4) dx2 dx As with the Navier–Stokes and Euler equations, (2.2.4) differs from (2.2.3) by the addition of a small constant times a higher-order term, as well as by the necessary addition of an extra boundary condition. The solution of (2.2.4) is verified to be y= 1−a 1 − e−1/ 1 − e−x/ + ax.

For 0 < a < 1, y is graphed in Figure 2.2.1. y y =1 solution with ε = 0

y =1–a solution with ε > 0

ε

x =1

x

Figure 2.2.1. Comparing solutions of (2.2.3) and (2.2.4).

For small and x > , the two solutions are close together. The region where they are drastically different is confined to the interval [0, ], which

2.2 Boundary Layers

69

we can refer to as the boundary layer . Note that as → 0, the boundary layer shrinks to zero, but the maximum difference between the solutions remains constant. Next we consider a special case in which (2.2.1) and (2.2.2) can be solved and the boundary layer can be seen explicitly. Example 2 Consider two-dimensional flow in the upper half-plane y ≥ 0, and suppose the boundary y = 0 (“the plate”) is rigid and the velocity at y = ∞ is parallel to the x-axis and has magnitude U . See Figure 2.2.2. y U rigid plate x

Figure 2.2.2. Flow over a flat plate.

We seek a flow satisfying (2.2.1) that is parallel to the plate and independent of x, that is, u has the form u = (u(y, t), 0), and, with constant pressure, so grad p = 0. The appropriate solution of Euler’s equations (2.2.2) is the solution u(x, y, t) = (U, 0). The boundary conditions for (2.2.1) are u(0, t) = 0 and u(∞, t) = U. Under the preceding conditions, u·∇u = u∂x u + v∂y u = 0, and so equation (2.2.1) reduces to ∂u ∂2u (2.2.5) = ν 2, ∂t ∂y
1 where ν = R . We can eliminate one variable by the following scaling argument: If L and T are length and time scales, respectively, the transformation y = y/L, t = t/T brings (2.2.5) to the form

∂u νT ∂ 2 u = 2 . ∂t L ∂y 2

(2.2.6)

If L2 /T = 1, then (2.2.5) and (2.2.6) are the same equation and the boundary conditions are unaltered as well. Thus, if equation (2.2.6) with the

70

2 Potential Flow and Slightly Viscous Flow

preceding boundary conditions has a unique solution, we must have y t , = u(y, t) L T √ Picking T = t and L = t , we obtain u u y √ ,1 t if L2 = t.

= u(y, t).

√ Thus, u depends only on y and t through the combination y/ t. Set η = √ y/(2 νt) and let u(y, t) = U f (η). Then, substitution in (2.2.5) yields the equation f (η) + 2ηf (η) = 0 with f (∞) = 1, f (0) = 0. (2.2.7)

Integration of f + 2ηf = 0 gives f (η) = ce−η ,
2

c a constant.

The solution of (2.2.7), shown in Figure 2.2.3, is u = U erf(η), where the error function is defined by 2 erf(η) = √ π u U 0.8U 0.6U 0.4U 0.2U 0 boundary layer 1 2 3 4 η η 0

e−s ds.
2

Figure 2.2.3. The boundary layer has thickness proportional to (t/R)1/2 .

In matching the boundary condition at infinity, we have used the following basic identity for the error function: √ ∞ 2 π e−s ds = . 2 0

2.2 Boundary Layers

71

The graph of u as a function of y for fixed t > 0 is obtained by rescaling the above graph. The region in which u departs significantly from the constant Euler flow is again called the boundary layer and is illustrated in the figure. By √ √ the preceding scaling argument, this thickness is proportional √ to √ νt = t/ R . Thus, for fixed time, the boundary layer decreases as 1/ R , which was our first contention. Our second contention was that the modified flow in the boundary layer may separate from the body. We can give a heuristic argument justifying this contention by considering potential flow around a cylinder (Figure 2.2.4). By Bernoulli’s theorem, y B A D C x

Figure 2.2.4. Illustration for the discussion of separation.

the highest pressure on the boundary occurs at the stagnation points A and C. These points are, in effect, points of high potential energy. As a fluid particle moves from A to B, its pressure drops and its velocity increases. As it moves from B to C it recoups its pressure and loses velocity. This is similar to the following situation: As a perfect frictionless bicycle goes down a hill into a valley, its momentum will enable it to climb another hill of equal height. Imagine a boundary layer is created in the flow. Within this layer the fluid is slowed down by friction. Therefore, fluid traversing the arc AC near the boundary may not have enough kinetic energy left to reach point C. It will leave the boundary somewhere between A and C. This argument only hints at what is, in reality, very complicated. For instance, the presence of viscosity near the boundary not only slows the fluid down by frictional forces, but it destroys the validity of Bernoulli’s theorem! Nevertheless, the preceding argument probably does have some

72

2 Potential Flow and Slightly Viscous Flow

merit; in fact, experiments show that flows really can undergo boundary layer separation. Our third contention was that the boundary layer produces vorticity. To see this, imagine flow near the boundary, as shown in Figure 2.2.5. y A circulation D

B C

x

Figure 2.2.5. Vorticity is produced near the boundary.

The circulation around the contour ABCD is u · ds =
ABCD DA

v dy +
AB

u dx −
CB

v dy −
DC

u dx.

Because u = 0 on the boundary, ∂x u = 0; so, by div u = 0, ∂y v = 0. Thus, because v = 0 on the boundary, v is small near the boundary. Let us, in fact, assume more specifically that v is small compared to the value of u along AB and that u is nearly zero along DC. (Near points of separation these assumptions might not be valid.) Thus, u · ds ∼ =
ABCD AB

u dx > 0

and so the flow has circulation and therefore vorticity. Next we shall make some approximations in the Navier–Stokes equations that seem to be intuitively reasonable near the boundary, away from points of separation. This will yield the Prandtl boundary layer equations. We now consider two-dimensional incompressible homogeneous flow in the upper half-plane y > 0. We write the Navier–Stokes equations as  1 ∂t u + u ∂x u + v ∂y u = −∂x p + ∆u,     R    1  ∂t v + u ∂x v + v ∂y v = −∂y p + ∆v,  R (2.2.8)     ∂x u + ∂y v = 0,      u = v = 0 on ∂D.

2.2 Boundary Layers

73

Let us assume that the Reynolds number R is large and that a boundary layer of thickness δ develops. From our previous discussions, we expect that √ δ ∼ 1/ R , where ∼ means “of the same order as.” See Figure 2.2.6. y main flow

boundary layer

δ x

Figure 2.2.6. The boundary layer thickness δ is of order 1/R1/2 .

We also expect, as discussed earlier, v to be small in the boundary layer. We want to perform a change of variables that will focus attention on the thin layer near the wall in which viscous effects are important. The layer is thin, of order δ. Write y = y/δ, so that when y varies between 0 and δ, y varies between 0 and 1. We expect that u is of order 1 at a distance of order δ from the wall, because at the edge of the layer, the flow should be effectively inviscid by definition (we are focusing on the region where the viscous effects are important). We claim that v is small at a distance δ from the wall. Indeed, at the wall, v = 0 and so inside the fluid, v(x, y) ≈ v(x, 0) + y ∂v ∂v =y ; ∂y ∂y

y is of order δ at most, and ∂v/∂y is of order 1 because ∂v ∂u + =0 ∂y ∂x and ∂u is of order 1. ∂x Thus, we perform the change of variables x = x, y = y , δ t = t, u = u, v = v , δ p = p.

74

2 Potential Flow and Slightly Viscous Flow

Equations (2.2.8) become ∂t u + u ∂x u + v ∂y u = −∂x p + δ ∂t v + δu ∂x v + δv ∂y v 1 1 = − ∂y p + δ R
2 δ ∂x v +

1 R

2 ∂x u +

1 2 ∂ u δ2 y δ 2 ∂ v δ2 y

,

                 ,                 (2.2.9)

∂x u + ∂y v = 0, u =v =0 on the x-axis.

We now collect the largest terms in each equation and eliminate the √ lower-order ones. If we had not known that δ ∼ 1/ R, we would have been able to deduce it from the scaled equation. Indeed, the purpose of the change of variables y = y/δ, etc., is to focus attention on the region in which the viscous wall effects are important and in which the transition from √ wall region to the outer region takes place. If δ is smaller than the O(1/ R), then only the viscous effects survive as the dominant terms in our equation (2.2.9)1 , and so we are looking too close to the wall. If δ is √ altogether and we are larger than O(1/ R) the viscous effects disappear √ looking at a region that is too large. Thus, δ = O(1/ R) is the only choice. From these arguments, it is plausible to believe that in the boundary layer, a good approximation to the Navier–Stokes equations is achieved by the following equations:  1 2 ∂t u + u ∂x u + v ∂y u = −∂x p + ∂y u,     R    ∂y p = 0, (2.2.10)   ∂x u + ∂y v = 0,      u = 0 = v on the x-axis. These are called the Prandtl boundary layer equations. One usually hopes, or assumes, that the following is true: If u(x, y, t) is a solution of the Navier–Stokes equations and up (x, y, t) is the corresponding solution of the boundary layer equations with the same initial conditions, then for some constant α > 0 and some constant C, u(x, y, t) − up (x, y, t) ≤ C Rα

for 0 ≤ y ≤ δ as R → ∞. Here · is a norm (to be found) on the space of velocity fields. The boundary layer equations have been quite successful and so one expects some estimate like this to be valid, at least under

2.2 Boundary Layers

75

some reasonable conditions. While some partial results verify this,6 a full mathematical proof under satisfactory hypotheses is not available. Let us derive a few consequences of the boundary layer equations (2.2.10). First of all, (2.2.10)2 implies that p is a function of x alone. This implies that if the pressure has been determined outside the boundary layer, (2.2.10)2 determines it inside the boundary layer. Secondly, let us derive an equation for the propagation of vorticity in the boundary layer. From the derivation of the equations, we see that ∂x v is to be neglected in comparison to ∂y u, and thus we have ξ = ∂x v − ∂y u = −∂y u. Therefore, from (2.2.10)1 , ∂t ξ = −∂y ∂t u = −∂y −∂x p + = 1 2 ∂ u − u ∂x u − v ∂y u R y

1 2 ∂ ξ − v ∂y ξ − u ∂x ξ + (∂x u + ∂y v)∂y u R y 1 2 = ∂y ξ − v ∂y ξ − u ∂x ξ, R because ∂x u + ∂y v = 0. Thus, Dξ 1 2 = ∂y ξ Dt R (2.2.11)

describes the propagation of vorticity in the boundary layer. We can interpret (2.2.11) by saying that the vorticity is convected downstream and, at the same time, is diffused vertically into the fluid. For the full Navier–Stokes equations in two dimensions, recall that Dξ 1 = ∆ξ Dt R (see equation (1.3.11) of Chapter 1). In the boundary layer approximation (2.2.11) only the y part of the Laplacian has survived. Notice also that in the full equations we recover u from ξ by the equations ∆ψ = −ξ and u = ∂y ψ, v = −∂x ψ with ψ = 0 on ∂D,

6 See O. A. Oleinik, Soviet Math. Dokl. [1968], for existence and uniqueness theorems. P. C. Fife (Arch. Rational Mech. Anal. 28 [1968], 184) has proved the above estimate for α = 1/2 assuming a forward pressure gradient, steady flow, and some other technical conditions. The boundary layer equations for curved boundaries are more complicated (see below and S. L. Goldstein, Modern Developments in Fluid Mechanics (two volumes), Dover [1965]). To the authors’ knowledge, no theorems of existence or approximation have been proved in this case.

76

2 Potential Flow and Slightly Viscous Flow

(see equation (1.3.12) of Chapter 1). In the boundary layer approximation the situation is much simpler: y ξ = −∂y u,

so

u=−
0

ξ dy.

Notice that the flow past the flat plate described at the beginning of this section is not only a solution of the Navier–Stokes equations but is a solution of the Prandtl equations as well. That solution was u(y, t) = U erf y √ 2 νt .

The vorticity for this flow ξ = −∂y u is shown in Figure 2.2.7. The vorticity is everywhere negative corresponding to clockwise circulation, which agrees with our intuition. y δ

ξ

Figure 2.2.7. Vorticity in the boundary layer.

Our derivation of the boundary layer equations assumed the wall was flat. For curved walls one can imagine using coordinate systems with y = 0 on the wall and x = constant giving the normal to the wall. One can show (see the Goldstein reference in the preceding footnote) that (2.2.10)1 , (2.2.10)3 and (2.2.10)4 are unaltered, but (2.2.10)2 should be modified to ku2 = ∂y p, (2.2.10)2

where k is the curvature of the wall. See Figure 2.2.8. This means that p is no longer a function of x alone. Of course, for δ small (R large) and fixed curvature of the wall, the pressure change across the boundary layer is negligible. However, (2.2.10)2 may still be significant in affecting the dependence of p on x. We have not yet discussed what boundary conditions should (or could) be imposed on the solution of the Prandtl equations beyond those at the

2.2 Boundary Layers

77

y

boundary layer wall

δ

x

Figure 2.2.8. Boundary layer for a curved wall.

wall, u = v = 0; these boundary conditions are essential for the derivation of the equations to make sense. We shall see later that far from the body, we may impose u but not v. Boundary conditions may also be required on the left and on the right. For a discussion, see the next section. Now we want to discuss how a knowledge of boundary layer theory might be useful in constructing solutions of the Navier–Stokes equations. The overall method is to use the Prandtl equations in the boundary layer, the inviscid equations outside the boundary layer, and to try to match these two solutions to produce an approximate solution to the Navier–Stokes equations. To carry out such a program, there are several alternative procedures, four of which are listed here. Strategy 1 Match the actual velocity field u of the two solutions at a √ distance δ ∼ 1/ R from the wall (Figure 2.2.9). y u Euler solution match u here

u

Prandtl solution δ x

Figure 2.2.9. Matching solutions of the Euler and Prandtl equations.

Difficulties On the line y = δ the flows may not be parallel, so matching u or some other flow quantity may not be easy.

78

2 Potential Flow and Slightly Viscous Flow

Strategy 2 Require that the Euler solution outside the layer and the Prandtl solution match as R → ∞. As R → ∞, the layer becomes thinner; we can impose v at the boundary for the Euler equations, but we cannot impose u. We can impose u at the edge of the layer and obtain a solution of the Prandtl equations, but we cannot impose v. We may wish to identify the vertical velocity v at the edge of the layer obtained from the Prandtl equations with the boundary condition required for the Euler equations, and conversely, use the tangential velocity v obtained from the Euler equations as the velocity at the edge of the layer imposed on the Prandtl equations. Difficulties The calculations are difficult, and the results are insufficiently informative when separation occurs (see below). Strategy 3 We can match the inner (boundary layer) and outer (Euler) solutions over a region whose height is of order 1/Ra , 0 < a < 1/2 (Figure 2.2.10). The matching can be done by some smooth transition between the Euler and Prandtl solutions.

inviscid region 1/Rα solutions agree 1/R1/2 boundary layer x

Figure 2.2.10. Matching method in strategy 3.

Difficulties The calculations are enormously complicated and require choices (of a for example) for which a rationale is not always readily available. Strategy 4 As both the Euler and Prandtl equations are solved and matched by some method, check (by means of a computer) whether or not the assumptions underlying the derivation of the Prandtl equations are valid. If not, modify the equations. Difficulties Same matching problems as before, but also modified boundary layer equations might be needed.

2.2 Boundary Layers

79

Techniques such as these are in fact fairly extensively developed in the literature.7 The scope of the difficulties involved, both theoretical and practical, should be evident. It is observed experimentally that when a boundary layer develops near the surface of an obstacle, the streamlines may break away from the boundary. This phenomenon is called boundary layer separation. Figure 2.2.11 shows a typical velocity profile near the point of separation. Contemplation of this velocity field shows that it is plausible to identify the point of separation C with points where the circulation vanishes: ∂u ∂n = ξ = 0.
C

There are, however, no known theorems that apply to this situation that guarantee such points must correspond to separation points. n u u u u

A

B

C separation point

back flow D separation line

Figure 2.2.11. Separation of the boundary layer.

The phenomenon of separation can be used to justify our model of almost potential flow introduced in the previous section (see Figure 2.1.8). The two curves of vorticity springing from the boundary are imagined to arise from boundary layer separation. One would hope that boundary layer theory would predict the strength of these vortex lines and their points of separation. Near the boundary, one should also be able to compute the pressure and hence the form drag − boundary pn ds.

This is the drag due to actual pressure forces on the obstacle. There will also be a skin drag due to frictional forces in the boundary layer. These
7 See, for instance, M. Van Dyke, Perturbation Methods in Fluid Mechanics, Parabolic Press.

80

2 Potential Flow and Slightly Viscous Flow

two give the total drag. If one computes the drag from a nearly potential approximation, the skin drag turns out to be zero, because Euler’s equations do not take into account the viscous forces. For large R, the skin drag is very often negligible compared to the form drag. Example This example concerns steady boundary layer flow on a flat plate of infinite width. Consider a flat plate on the xz plane with its leading edge coinciding with x = 0, as shown in Figure 2.2.12. Let the y-axis be normal to the plate so that the fluid lies in the domain given by y > 0 and the flow is two dimensional. We then let the plate be washed by a steady uniform flow of velocity U in the positive x-direction. We assume that a steady boundary layer is developed after a short distance off the leading edge. Because the pressure is uniform outside the boundary layer, it is reasonable to seek a solution satisfying ∂x p = 0 in the boundary layer. Because ∂x p = 0, p is a constant everywhere. Hence, the steady boundary layer equations are 2 u ∂x u + v ∂y u = ν ∂y u, (2.2.12) ∂x u + ∂y v = 0. (Here we have written ν instead of 1/R.) The boundary conditions are given by u(x, 0) = 0 v(x, 0) = 0 u(x, y) → U for x > 0, for x > 0, as y → ∞.

Because the flow is incompressible, there is a stream function related to u by the usual equations u = ∂y ψ, We assume ψ to be of the form ψ(x, y) = where U . νx This form may be derived through the use of a scaling argument similar to the one used earlier in §2.1. Thus, η=y u = U f (η), v = 1 2 Uν (ηf (η) − f (η)). x (2.2.13) √ v = −∂x ψ.

νU x f (η),

Substituting u, v from (2.2.13) into (2.2.12) and using the boundary conditions, we find the equations f f + 2f =0 (2.2.14)

2.2 Boundary Layers

81

with boundary conditions f (0) = f (0) = 0, f (∞) = 1.

Equation (2.2.14) has a numerically computable solution. In the fluid solution given by (2.2.13), notice that the velocity u depends only on η. The level curves of u are exactly the lines where η = constant. But level curves of η are of the form y 2 = (constant)x, that is, parabolas, as shown in Figure 2.2.12. y velocity profile for u

U main flow U 0 leading edge

level curve for u

plate

x

Figure 2.2.12. Boundary layer flow over a flat plate.

Exercises
Exercise 2.2-1 Show that the solution of y −y =0 in the interval (−1, 1) with the boundary conditions y(−1) = 0, is


y(1) = 1


e(1+x)/ − e−(1+x)/ √ √ y(x) = e2/ − e−2/ Plot the graph of y(x) for = 0.1, 0.01 and 0.001.

.

Exercise 2.2-2 Derive the solution in Exercise 2.2-1 from that for (2.2.4) in the text by a scaling argument. Exercise 2.2-3 Plot the shape of the graph of the function f (η) that solves (2.2.14) with f (0) = f (0) = 0, f (∞) = 1.

82

2 Potential Flow and Slightly Viscous Flow

2.3

Vortex Sheets

This section describes a simple version of an algorithm for modifying a flow described by the Euler equations near the boundary to simulate the effects of boundary layers and thereby obtain an approximation to the Navier–Stokes equations. The method utilizes some concepts from probability theory, so we begin by summarizing some relevant facts, mostly without proofs.8 The possible outcomes of an experiment (such as throwing a die) form the points in a sample space S. A subset E of S is called an event. We assume that to each event is assigned a probability P (E), a number between 0 and 1 which intuitively represents the fraction of times an outcome in E will occur if the experiment is repeated many times. We assume, therefore, that P (S) = 1. Moreover, if two events, E1 and E2 , are disjoint, that is, E1 ∩ E2 = ∅, then P (E1 ∪ E2 ) = P (E1 ) + P (E2 ). (Technically, we assume P is a measure on S with total mass 1.) Two events, E1 and E2 , are called independent if P (E1 ∩ E2 ) = P (E1 ) · P (E2 ). Intuitively, two events are independent if the occurrence of one of them has no effect on the probability of the occurrence of the other one. (For instance, in the toss of two dice marked #1 and #2, the events “a two on #1” and “a three or a four on #2” are independent.) A random variable is a mapping η : S → R, interpreted as a number attached to the outcome of an experiment. The expectation or mean of η is defined by E(η) =
S

η dP.

For instance, if S = {s1 , . . . , sN } and the probability of si occurring is pi , then n E(η) = i=1 η(si )pi

Suppose there is a function f on the real line such that the probability b of η lying between a and b is a f (x) dx, that is, b P ({s ∈ S | η(s) ∈ [a, b]}) = a f (x) dx.

8 For details see, for instance, J. Lamperti [1966] Probability: A Survey of the Mathematical Theory, Benjamin.

2.3 Vortex Sheets

83

Then we say that η has the probability density function f . Clearly, ∞ f (x) dx = 1. Also, one can show that −∞


E(η) =
−∞

xf (x) dx.

The variance of η is defined by Var(η) = E [η − E(η)] and the standard deviation by σ(η) = Var(η).
2

= E(η 2 ) − [E(η)]2

Two random variables, η1 and η2 , are called independent if for any two sets A1 , A2 in the real line R, the events {s ∈ S | η1 (s) ∈ A1 } and {s ∈ S | η2 (s) ∈ A2 }

are independent. For independent random variables, one has E(η1 η2 ) = E(η1 )E(η2 ) and Var(η1 + η2 ) = Var(η1 ) + Var(η2 ). (From the definition, E(η1 + η2 ) = E(η1 ) + E(η2 ) is always true.) The law of large numbers states that if η1 , η2 , . . . , ηn are random variables that are independent and have the same mean and variance as η, then n 1 E(η) = lim ηi . n→∞ n i=1 Part of the theorem is that the right-hand side is a constant. This result sheds light on our intuition that E(η) is the average value of η when the experiment is repeated many times. The meaning of the standard deviation is illuminated by Tchebysheff ’s inequality : If σ is the standard deviation of η, 1 P ({s ∈ S | |η(s) − E(η)| ≥ kσ}) ≤ 2 k for any number k > 0. For example, the probability that η will deviate from its mean by more than two standard deviations is at most 1/4. If a random variable η has the probability density function f (x) = √ 1 2πσ 2 e−(x−a)
2

/2σ 2

,

84

2 Potential Flow and Slightly Viscous Flow

we say that η is Gaussian. One can check that E(η) = a and Var(η) = σ 2 . If η1 and η2 are independent Gaussian random variables, then η1 + η2 is Gaussian as well. The central limit theorem states that if η1 , η2 , . . . are independent random variables with mean 0 and variance σ 2 , then the √ probability density function of (1/ n)(η1 + · · · + ηn ) converges to that of a Gaussian with mean 0 and variance σ 2 as n → ∞. For instance, if measurement errors arise from a large number of independent random variables, the errors can be expected to have a Gaussian distribution. In many problems, however, one does not get a Gaussian distribution for the errors because they arise from nonindependent sources. Next we show how Gaussian random variables can be used in the study of the heat equation for an infinite rod: vt = νvxx , −∞ < x < ∞, t ≥ 0. (2.3.1)

Here v represents the temperature in the rod as a function of x and t, and ν represents the conductivity of the rod. If v is given at t = 0, then (2.3.1) determines it for t ≥ 0. If initially v(x, 0) = δ(x), a delta function at the origin, then the solution of (2.2.1) is given by G(x, t) = √ 1 exp 4πνt −x2 4νt . (2.3.2)

This is the Green’s function for the heat equation (see any textbook on partial differential equations). We can interpret the function (2.3.2) from a probabilistic point of view in two ways as follows. Method 1 Fix time at t, and place N particles at the origin. Let each of the particles “jump” by sampling the Gaussian distribution with mean zero and variance 2νt. Thus, the probability that a particle will land between x and x + dx is 1 −x2 √ dx. exp 4νt 4πνt If we repeat this with a large number of particles, we find number of particles between 1 x and x + dx at time t exp =√ N →∞ N dx 4πνt lim −x2 4νt . (2.3.3)

It is convenient to think of each particle as having mass 1/N , so the total mass is unity. Method 2 We split up the time interval [0, t] into l pieces, each with length t = t/l, and carry out the procedure in a step-by-step manner. Again, place N particles at the origin. Let them undergo a random walk,

2.3 Vortex Sheets

85

that is, the position of the ith particle at time m∆t (i = 1, . . . , N ; m = 1, . . . , l) is m xm+1 = xm + ηi , i i 0 xi = 0,

(2.3.4)

m where ηi are independent Gaussian random variables, each with mean 0 and variance 2ν∆t. The final displacement of the ith particle is the sum of its displacements, and has a Gaussian distribution with mean 0 and variance l × 2ν∆t = 2νt. Thus, the distribution of particles at time t is given again by formula (2.3.3), and methods 1 and 2 are equivalent. Next consider the solution v(x, t) of the heat equation with given initial data v(x, 0) = f (x). The solution is ∞

v(x, t) =
−∞

G(x, x , t)f (x ) dx ,

(2.3.5)

−(x − x )2 1 exp . 4νt 4πνt This general solution has a probabilistic interpretation as well. Instead of starting N particles at the origin, start N randomly spaced particles on the line, at positions, say, x0 , i = 1, . . . , N , and assign to the ith particle i the mass f (x0 ) i . N Let these particles perform a random walk as in (2.3.4)1 , keeping their mass fixed. Then after l steps as in method 2, the expected distribution of mass on the real line approximates (2.3.5). In this process the total mass of the particles remains constant. This corresponds to the fact that G(x, x , t) = √
∞ ∞

where

∂t

v(x, t) dx = ν
−∞ −∞

vxx (x, t) dx = 0

(assuming vx → 0 as x → ±∞). Of course, one’s intuitive feeling that the solutions of the heat equation decay is also correct. Indeed,


∂t

v 2 (x, t) dx =

∞ −∞

−∞

2νvvxx dx = −2ν

∞ −∞

(vx )2 dx < 0.

The decay of v 2 dx (which occurs while v dx remains constant) is accomplished by spreading. As time advances, the maxima of the solution decay and the variation of the solution decreases. To see intuitively why the integral of v 2 decreases, consider the two functions v1 = 2, − 1 ≤ x ≤ 2 0, elsewhere
1 2

and v2 =

1, −1 ≤ x ≤ 1 . 0, elsewhere

86

2 Potential Flow and Slightly Viscous Flow

The function v2 is more “spread out” than v1 . This is reflected by the calculations v2 dx = but
2 v2 dx = 2 and 2 v1 dx = 4.

v1 dx = 2,

Note that as time unfolds, the variance of the random walk that is used to construct the solution increases, whereas the integral of v 2 , which is related to the variance of v, decreases. The variance of the random walk increases as the solution spreads out, whereas the integral of v 2 decreases because the range of values assumed by v decreases. Next consider the heat equation on the halfline x ≤ 0 with boundary condition v(0, t) = 0. The Green’s function for this problem is G∗ (x, x , t) = G(x, x , t) − G(x, −x , t), where G(x, x , t) is given by (2.3.5)2 . That is, G∗ satisfies G∗ (0, x , t) = 0, and The solution of vt = νvxx , with v(x, 0) = f (x) and v(0, t) = 0 is given in terms of the initial data and the Green’s function by


(2.3.6)

G∗ (x, x , 0) = δ(x − x ),

2 ∂t G∗ (x, x , t) = ν∂x G∗ (x, x , t).

x ≤ 0, t > 0,

v(x, t) =
−∞

G∗ (x, x , t)f (x ) dx .

(2.3.7)

The random walk interpretation for (2.3.6) is obtained by superposing the random walks that generate Green’s function for the whole line; namely, the distribution G∗ (x, x , t) is obtained by starting N particles with weight 1/N at x and N with weight −1/N at −x , and letting them all walk by, say, method 2. There is an analogous interpretation for (2.3.7), as shown in Figure 2.3.1. Random walk methods will now be applied to vortex sheets. Consider flow past an infinite flat plate at rest. In §2.2 (see Figure 2.2.3), we saw that the flow is u = u(y, t), v = 0,

2.3 Vortex Sheets

87

G(x, x', t) x x –G(x, –x', t)

x'

–x'

x'

–x'

x G*(x, –x', t)

Figure 2.3.1. Green’s function G∗ is obtained by superimposing two random walks.

where u = U erf(η), and the vorticity is

with η = √ −y 2 4νt

y , 4νt

(2.3.8)

U ξ = −√ exp 4νt The velocity and vorticity satisfy ut = νuyy , u(0, t) = 0,

.

(2.3.9)

u(∞, t) = U,

ξt = νξyy .

(2.3.10)

The boundary conditions for the vorticity are not explicit; they must be determined from the fact that the velocity vanishes on the boundary. To reconstruct this solution using random walks, we define a vortex sheet of intensity ξ as follows. As in Figure 2.3.2, it consists of a flow parallel to the x-axis such that u jumps by the amount ξ as y crosses a line, say y = y0 ; that is, u(y0 +) − u(y0 −) = −ξ. The flow described by equation (2.3.8) near the boundary looks as shown in Figure 2.3.3. As t → 0+, the solution approaches the constant value U for y > 0 and yet u = 0 on the x-axis; that is, as t → 0+, the solution approaches a vortex sheet on the x-axis with intensity −U . We replace this vortex sheet by N “smaller” vortex sheets, each of intensity −2U/N . Allow each of these smaller vortex sheets to undergo a random walk in the y-direction of the form m+1 m m 0 yi = yi + ηi , yi = 0, m where ηi is chosen from a Gaussian distribution with mean zero and variance 2ν∆t, where, as in method 2, ∆t = t/l, m = 1, 2, . . . , l. We claim

88

2 Potential Flow and Slightly Viscous Flow y + u(y0 )

y0
– u(y0 )

x

Figure 2.3.2. A vortex sheet in the plane: u jumps across the line y = y0 .

y u=U

boundary layer x

Figure 2.3.3. The flow (2.3.8).

that for N large the distribution of vorticity constructed this way and the resulting velocity


u(y, t) = U + y ξ(y, t) dy

(2.3.11)

satisfy the correct equations (2.3.10). That ξ and hence u satisfy the heat equation follows from our random walk solution of the heat equation. What requires explanation is why u should vanish on the boundary. To see this, note that on the average, half the vortex sheets are above the x-axis and half below, because they started on the x-axis and the Gaussian distribution has mean zero. Thus,


u(0, t) = U +
0

ξ(y, t) dy,

or, in our discrete version, we have on the average
N/2

u(0, t) = U + i=1 ξi .

2.3 Vortex Sheets

89

But the intensity of the ith vortex sheet is ξi = −2U/N , and therefore u(0, t) = 0. (By symmetry in the x-axis the random walks below the x-axis can be ignored; whenever a vortex sheet crosses the x-axis from above, we can “reflect” it back. On the average, this will balance those vortex sheets that would have crossed the x-axis from below during the course of their random walk. This device can save half the computational effort.) The random walk method based on vortex sheets will now be generalized to the solution of the Prandtl boundary layer equations: ξt + uξx + vξy = νξyy , ∂x u + ∂y v = 0, with u = 0 = v on the x-axis (the boundary of the region), y ≥ 0, and u(∞) = U. (We shall see shortly why v cannot be prescribed at ∞.) The overall flow will be approximated at t = 0 by a collection of N vortex sheets of finite width h, extending from xi − h/2 to xi + h/2 with ycoordinate yi and with strength ξi . This approximation is in the same spirit as the point vortex approximation discussed in §2.1 (see Figure 2.3.4). y h yi xi ξi x

Figure 2.3.4. A collection of vortex sheets.

To move these sheets, we divide the time [0, t] into l parts of duration ∆t = t/l and proceed in a step-by-step manner. The advancement in time from time t to t + ∆t takes place by means of the following algorithm: i The vortices move according to a discrete approximation to the Euler flow ∂t ξ + uξx + vξy = 0, ∂x u + ∂y v = 0. Vorticity is added by placing new vortex sheets on the boundary so that the resulting flow has u = 0 = v on the boundary (creation of vorticity).

ii

90

2 Potential Flow and Slightly Viscous Flow

iii

The vortex sheets undergo a random walk as described in the flat plate example (including reflections) to approximate the solution of the equation ξt = νξyy and preserve the boundary conditions u = 0 = v.

iv

Time is advanced by ∆t and one goes back to step i of the procedure, etc., until time t is reached.

Notice that the number of vortex sheets will increase in time; this corresponds to the fact that vorticity is created in boundary layers. We now give some more details on this procedure. First, we discuss how the vortex sheets move according to the Euler flow. The velocity field satisfies (see §2.2): u(x, y) = u(∞) − y ∞

∂u dy = u(∞) + ∂y



ξ dy, y or in a discrete version, the velocity of the ith vortex is u(xi , yi ) = u(∞) + j ξj ,

(2.3.12)

where the sum is over all vortex sheets such that yj > yi and |xi −xj | < h/2; that is, all vortex sheets whose “shadow” (on the x-axis) engulfs (xi , yi ); see Figure 2.3.5. Equation (2.3.12) determines the u-component of the velocity field produced by the vortex sheets. y (xj,yj) (xk,yk) shadow

(xi,yi)

x

Figure 2.3.5. Vortex sheets and their shadows.

From incompressibility and the boundary condition v(x, 0, t) = 0, we find y y

v(x, y, t) = v(x, 0, t) +
0

ux (x, s, t) ds = ∂x
0

u(x, s, t) ds.

2.3 Vortex Sheets

91

This equation determines v in terms of u (and shows why we could not prescribe v at ∞). A discrete evaluation of y may be given by v(x, y, t) = − 1 h y u x+
0

h , s, t ds 2 y −
0

u x−

h , s, t ds . (2.3.13) 2

A more useful approximation is given by rewriting (2.3.13) directly in terms of the vortex strengths ξj , namely, vi (xi , yi , t) = where i I+ (xi , yi ) = +

1 i i (I − I− ), h +
∗ yj ξj , ∗ yj ξj .

(2.3.14)

and i I− (xi , yi ) = −

with Here

∗ yj = min(yj , yi ). +

means that one is to sum over all j = i for which xj − xi + h 2 < h , 2

and



means sum over all j for which xj − xi − h 2 < h . 2

The expression (2.3.14) comes about by remembering that as we move vertically from the x-axis, u remains constant and then jumps by an amount −ξj when a vortex sheet of strength ξj is crossed. This leads us directly from (2.3.13) to (2.3.14). We can summarize all this by saying that in our first step of the algorithm the ith vortex sheet is moved by xm+1 = xm + ui ∆t i i m+1 m yi = yi + vi ∆t,

(2.3.15)

where ui is given by (2.3.12) and vi by (2.3.14). The new velocity field is now determined by the same vortex sheets, but at their new positions. This new velocity field satisfies v = 0 on the x-axis by construction, and u(∞) = U . However, because the vortex sheets move, even if u = 0 on the x-axis at the beginning of the procedure, it need not remain so. This is another aspect of the main fact we are dealing with: The boundary

92

2 Potential Flow and Slightly Viscous Flow

conditions for the Euler equations, u · n = 0, are different from those for the Navier–Stokes or Prandtl equations, u = 0. The purpose of the second step in our procedure is to correct the boundary condition. This may be done as follows: Divide the x-axis into segments of length h, and suppose that at the center Pl of one of these segments u = ul = 0. Then at Pl place one (or more) vortex sheets, the sum of whose intensities is 2ul . This will guarantee that (on the average) u = 0 on the x-axis in the new flow. In the third step in our procedure, we add a random y component to positions (xi , yi ) of the existing (as well as of the newly created) vortex sheets: xm+1 = xm + ui ∆t i i m+1 m m yi = yi + vi ∆t + ηi ,

(2.3.16)

m where ηi are Gaussian random variables with mean 0 and variance 2ν∆t. Intuitively, the vortex sheets move about in ideal flow together with a random y component simulating viscous diffusion. Vortex sheets newly created to accommodate the boundary conditions diffuse out from the boundary by means of the random component and then get swept downstream by the main flow. This mechanism then has all the features a boundary layer should have as we discussed at the beginning of §2.2.9 Consider the problem of flow past a semi-infinite flat plate, situated on the positive x-axis (Figure 2.3.6). Far enough along the positive x-axis, say at x = x2 , most of the vortex sheets that move downstream are replaced by vortex sheets coming from the left. Thus, there is little need to create much vorticity at x2 . However, at x1 (called the leading edge) any vortex sheet is immediately swept downstream by the flow with no replacement. Thus, we are constantly forced to create more vortex sheets at x = x1 to satisfy the noslip condition. One can see that most of the vorticity in the flow is created at the leading edge. In a time ∆t, how far does a vortex sheet move? In the x direction the displacement is proportional to (∆t); in the y direction the displacement m is (∆t)vi + ηi . The standard deviation tells us the length of the average “jump.” Thus, √ √ average jump in y direction ∼ var ∼ ∆t. =
9 The generation of vorticity is discussed is G. K. Batchelor [1967] An Introduction to Fluid Mechanics, Cambridge Univ. Press; J. Lighthill [1963] “Introduction to Boundary Layer Theory”, in Laminar Boundary Layers, edited by L. Rosenhead, Oxford Univ. Press. The vortex sheet model is due to A. J. Chorin, J. Comp. Phys. 27 [1978], 428. Some theoretical aspects are found in A. J. Chorin, T. J. R. Hughes, M. J. McCracken, and J. E. Marsden, Comm. Pure. Appl. Math. 31 [1978], 205; C. Anderson, J. Comp. Phys., 80 [1989], 72; C. Marchioro and M. Pulvirenti, Vortex Methods in 2-Dimensional Fluid Mechanics, Springer-Verlag, [1984]; K. Gustafson and J. Sethian, Vortex Flows, SIAM Publications [1991].

2.3 Vortex Sheets

93

y U

x x1 x2

Figure 2.3.6. Vortex sheets get created and swept downstream.

√ The term (∆t)vi is also proportional to ∆t, as can be seen from equation (2.3.14) and the fact that I+ and I− contain ys as factors. Suppose the flow is stationary; we may regard t as a parameter and we may eliminate it to find vertical displacement ∼ horizontal displacement.

This shows that the structure of the boundary on a semi-infinite flat plate is parabolic, as was derived earlier. Finally, we ask, are there more general constructions available, similar to this vortex sheet construction? For Euler’s equations, we have already described the model of point vortices satisfying dxi ∂H , = dt ∂yi dyi ∂H , =− dt ∂xi

where 4πH = ΣΓi Γj log rij . This is a system of ordinary differential equations; we can add a random term to take care of the diffusion as with the boundary layer equations. Indeed, this construction has been carried out and one can conceivably use it to study the convergence of the Navier– Stokes equations to Prandtl’s equations, boundary layer separation, and other questions of interest. In three dimensions a possible construction might use vortex filaments, but a discussion of these constructions is outside the scope of this book. These ideas are the starting points for the construction of practical numerical algorithms. There is a way of writing the above general scheme for solving the Navier– Stokes equations that sheds light on the mathematical structure of this and related methods.10 The methods are related to some basic facts about
10 Some additional references, in addition to those in the preceding footnote, are A. J. Chorin, J. Fluid Mech. 57 [1973], 785–796, J. E. Marsden, Bull. Amer. Math. Soc. 80 [1974], 154–158, G. Benfatto and M. Pulvirenti, Convergence of the Chorin-Marsden product formula in the half-plane, Comm. Math. Phys. 106 [1986], 427–458, and J. E. Marsden, Lectures on Mechanics, Cambridge University Press, [1992].

94

2 Potential Flow and Slightly Viscous Flow

algorithms and our purpose is just to provide a hint about them. Consider first a differential equation x = X(x) ˙ (2.3.17)

on Rn (eventually replaced by a more general space). To solve this equation, a computer would move a point forward during a time step ∆t = τ according to some rule. More precisely, an algorithm is a collection of maps Fτ : Rn → Rn . The associated iterative process is denoted xk+1 = Fτ (xk ). To be consistent with (2.3.17) one requires d Fτ (x) dτ = X(x). τ =0

(2.3.18)

Under some additional hypotheses one can establish convergence of the algorithm, namely, lim (Ft/n )n x0 = x(t), (2.3.19) n→∞ where x(t) is the solution of (2.3.17) with initial condition x0 . Of course, questions of rate of convergence and efficiency of the algorithm are important in practice, but are not discussed here. An important simple case is the following. Consider solving the linear system x = Ax + Bx, ˙ (2.3.20) where A and B are n × n matrices. The solution is x(t) = et(A+B) x0 , where eC = I + C + 1 C 2 + 2
1 3 3! C

(2.3.21) + ···

However, it may be easier to solve the equation by breaking it into two pieces: x = Ax and x = Bx and solving these successively. Doing so leads ˙ ˙ to the algorithm Fτ x = eτ A eτ B x, and yields the basic formula et(A+B) = lim (etA/n etB/n )n n→∞ (2.3.22)

(sometimes called the Lie–Trotter product formula). More generally, consider solving a nonlinear equation x = X(x) + Y (x). ˙ (2.3.23)

Let Ht (x0 ) denote the solution of (2.3.23) with initial condition x0 , so we thereby define a map Ht : Rn → Rn , called the flow map. Let Kt and Lt

2.4 Remarks on Stability and Bifurcation

95

be the corresponding flow maps for x = X(x) and x = Y (x). Then (2.3.22) ˙ ˙ generalizes to Ht (x0 ) = lim (Kt/n ◦ Lt/n )n (x0 ), (2.3.24) n→∞ where the power n indicates repeated composition. It is tempting to apply (2.3.24) to the Navier–Stokes equations, where Kt is the flow of the Stokes equation and Lt the flow of the Euler equations. First of all, consider the case of regions with no boundary (e.g., this occurs for the case of spatially periodic flows in the plane or space). In this case, the method works quite well and was used by Ebin and Marsden in 1970 to show the convergence in a suitable sense of the solutions of the Navier– Stokes equations to the Euler equations as the viscosity ν → 0. If a boundary is present, we need to modify the scheme because of the boundary conditions, as has been explained. Translated into a product formula, the modified scheme reads as follows. Ht (u) = lim (Kt/n ◦ Φt/n ◦ Lt/n )n u. n→∞ (2.3.25)

Here, • u is a divergence-free vector field, u = 0, on the boundary ∂Ω of the region Ω in question; • Lt is the flow of the Euler equation (boundary conditions u parallel to ∂Ω), which may be solved approximately by a vortex method; • Kt is the flow of the Stokes equation (boundary conditions u = 0 on ∂Ω), solved for example by a random walk procedure on vorticity; • Φt is a “vorticity creation operator” that maps a u ∂Ω to a u = 0 on ∂Ω by adding on a suitably constructed vorticity field to u whose backflow cancels u on ∂Ω; and • Ht is the flow of the Navier–Stokes equations. Convergence of (2.3.25) in reasonable generality is not yet proven, although special cases are given in the cited references. What formula (2.3.25) does is to make explicit the intuitive idea that vorticity is created on the boundary because of the difference in the boundary conditions between the Euler and Navier–Stokes equations, and that if the Reynolds number is high this vorticity is swept downstream by the Euler flow to form a structured or turbulent wake. The references should be consulted for details.

2.4

Remarks on Stability and Bifurcation

If R is the Reynold’s number, then limit R → 0 corresponds to slow or viscous flow. The other extreme, the limit R → ∞ concerns very fast or only

96

2 Potential Flow and Slightly Viscous Flow

slightly viscous flow. The situation for intermediate values of R involves many interesting transitions, or changes of flow pattern, and is the subject of much research by mathematicians and engineers alike. Here we shall give a few informal comments using the subject of nonlinear dynamics, complementing topics already discussed in this chapter. The concept of stability plays a central role in this discussion, so we begin with this concept. Let P be Rn (more generally a Banach or Hilbert space of functions), and let Ft denote a flow of a differential equation x = X(x) on P ; that is, ˙ d Ft (x) = X(Ft (x), t), dt where x ∈ P and t ≥ 0. Assume that we have a good existence and uniqueness theorem. If X is independent of t, we say it is autonomous, and if x0 is such that X(x0 ) = 0, we call x0 a fixed point of X. It follows from the initial condition F0 (x0 ) = x0 and uniqueness of solutions that x0 is a fixed point of X, i.e., Ft (x0 ) = x0 . Definition A point x0 is an asymptotically stable fixed point of the vector field X if there is a neighborhood U ⊂ P of x0 such that if x ∈ U , then Ft (x) → x0 as t → ∞. In the sequel we shall refer to such a point simply as a stable point. If one can choose U = P , we say that x0 is globally stable. See Figure 2.4.1. y p1

p2

p3

x

Figure 2.4.1. A system with two asymptotically stable points.

Next we give a basic stability result of Liapunov11 that is based on computing the spectrum of the linearization.
11 For the proof, see for example, M. Hirsch and S. Smale [1974] Differential Equations, Dynamical Systems, and Linear Algebra, Academic Press, or the book by Abraham, Marsden, and Ratiu listed in the Preface.

2.4 Some Remarks on Stability and Bifurcation

97

Liapunov Stability Theorem Let x0 be a fixed point of a smooth autonomous vector field X on P = Rn . Then x0 is a stable point if the eigenvalues of the n × n matrix DX(x0 ) of partial derivatives of X have real parts < 0. Example Suppose we have a linear autonomous system, that is, X = A, an n×n real matrix, which we assume for simplicity is diagonalizable. Then the solution of d Ft (x) = A · Ft (x), F0 (x) = x dt is Ft (x) = etA x, where etA = i=o ti Ai /i!. This series is absolutely convergent. Transforming to coordinates that diagonalize A, we get A = diag(λ1 , . . . , λn ) and so etA = diag(etλ1 , . . . , etλn ) where λ1 , . . . , λn are the eigenvalues of A (they occur in conjugate pairs because A is real). This indicates why the stability theorem is true, namely, if Re(λi ) < 0, i = l, . . . , n, then the matrix etA converges to the zero matrix as t → ∞. Observe that for the linear case, consideration of neighborhoods is unnecessary. The preceding example is helpful in understanding the nonlinear case as well. To see this, expand X by Taylor’s formula about x0 : X(x) = X(x0 ) + DX(x0 ) · (x − x0 ) + θ(x − x0 ), where X(x0 ) = 0, θ(x − x0 ) is o( x − x0 ), and · is the norm on P . Thus, for x sufficiently close to x0 , the term DX(x0 ) dominates the behavior of the flow of X. The stability theorem requires that the spectrum of the matrix DX(x0 ) lies entirely in the strict left half-plane of the complex plane C for x0 to be stable. Exercise Consider the vector field Xµ (x, y) = (y, µ(1 − x2 )y − x) on R2 . Determine if the fixed point (0, 0) is stable or unstable for various values of µ. We shall now examine how these concepts enable us to study the stability of solutions of the Navier–Stokes equations. Consider the example of flow in a pipe and assume v0 is a stationary solution of the Navier–Stokes equations, for example, the Poiseuille solution (Exercise 1.3-3). We are interested in what happens when v0 is perturbed, that is, when v0 → v0 + δv.


98

2 Potential Flow and Slightly Viscous Flow

Note that the disturbed solution of the Navier–Stokes equations, v0 + δv, will not be a stationary solution. This is part of our notion of stability, that is, v0 is stationary corresponding to v0 being a fixed point of a vector field on the space of all possible solutions of the Navier–Stokes equations, which we shall denote by P . For the present example, P would consist of the set of all divergence-free velocity fields satisfying the appropriate boundary conditions. The Navier–Stokes equations determine the dynamics on this space P of velocity fields. We write dv = X(v), dt v a curve in P,

where X does not depend on time unless the boundary conditions or body force are time dependent. In this abstract notation, the Navier–Stokes equations, boundary conditions, divergence-free condition, and so on, look like an ordinary differential equation. In the present circumstances X is an unbounded operator. However, we can still apply the stability concepts for dynamical systems on Rn , due to the following conditions: (1) The preceding equation really does determine the dynamics on P . This follows from the short t-interval existence and the uniqueness theorem for the Navier–Stokes equations. (2) The stability theorem really works for the Navier–Stokes equations.12 For condition (2), DX is obtained via differentiation and DX(v0 ) is a linear differential operator. Its spectrum will consist of infinitely many eigenvalues or it may be continuous. For v0 to be stable the entire spectrum must lie in the left half-plane of C. For example, flow in a pipe and Couette flow are stable (in fact globally stable) if the Reynolds number is not too big. Sometimes one is interested in the loss of stability of flows as the Reynolds number R is increased. In general, any such qualitative change in the nature of a flow is called a bifurcation. In this regard we consider X to be parametrized by R and study the behavior of the spectrum of DX(v0 ) as a function of this parameter. As R is increased, we anticipate that conjugate pairs of eigenvalues of DX(v0 ) may drift across the imaginary axis. In this case stability is lost and an oscillation develops. A major result dealing with this situation is given by the following theorem of Poincar´, Andronov, and Hopf.13 e
12 See for example, D. Joseph, Stability of Fluid Motion, I, II, Springer-Verlag and for the Euler equations, V. Arnold and B. Khersin, Topological Methods in Hydrodynamics, Springer-Verlag, 1997, and references therein. 13 References are J. E. Marsden and M. McCracken The Hopf Bifurcation, Springer Applied Mathematics Series, Vol. 19 [1976]; D. Sattinger Lectures on Stability and Bifurcation Theory, Springer Lecture Notes 309 [1973]. Books that approach this subject

2.4 Some Remarks on Stability and Bifurcation

99

Hopf Bifurcation Theorem (Nontechnical version) Let Xµ be a vector field (depending on a parameter µ) possessing a fixed point x0 for all µ. Assume the eigenvalues of DXµ (x0 ) are in the left half-plane of C for µ < µ0 , where µ0 is fixed. Assume that as µ is increased, a single conjugate pair λ(µ), λ(µ) of eigenvalues crosses the imaginary axis with nonzero speed at µ = µ0 (Figure 2.4.2 ). Then there is a family of closed orbits of X near µ0 . If the point x0 is stable for Xµ0 , then the closed orbits appear for µ > µ0 and are stable. For each µ > µ0 , near µ0 , there is one corresponding stable closed orbit and its period is approximately equal to Im(λ(µ0 )/2π). spectrum of the linearization at equilibrium

Im λ C

Re λ

Figure 2.4.2. A pair of eigenvalues crosses the imaginary axis in the Hopf bifurcation.

Roughly speaking, what this means is that when stability is lost, a stable point is replaced by a stable closed orbit. Translated to the case of fluid mechanics this means that a stationary solution (fixed point) of the Navier– Stokes equations is replaced by a periodic solution (stable closed orbit). For example, consider the case of flow around a cylinder. As R is increased, the stationary solution becomes unstable and goes to a stable periodic solution—the wiggly wake depicted in Figure 2.4.3. Such a development of periodic oscillations is the main content of the Hopf bifurcation theorem. It applies to a wide variety of physical and biological phenomena, and this ubiquity leads one to strongly suspect it is the mechanism underlying even complex fluid dynamical phenomena like the singing of transmission lines in a strengthening wind. We say that a bifurcation has occurred when R reaches a critical value where stability is lost and is replaced by oscillations. The term is used generally when sudden qualitative changes occur. In Hopf’s theorem the from a more physical point of view are C. C. Lin, The Theory of Hydrodynamic Stability, Cambridge [1955], and S. Chandrasekhar Hydrodynamic and Hydromagnetic Stability, Oxford [1961].

100

2 Potential Flow and Slightly Viscous Flow

Figure 2.4.3. The development of patterns like this is related to bifurcation, stability loss and symmetry breaking.

stability analysis fails at µ0 and stability must be determined by “higherorder” analysis; that is, the stability theorem, based on first-order analysis just fails at µ = µ0 . How to do this is discussed in the cited references. One should also note that in highly symmetrical situations such a Couette flow, the hypotheses of Hopf’s theorem (simplicity of the eigenvalues, or the nonzero speed assumption) can fail. In this context the theory of bifurcation with symmetry is needed.14 As R is increased, additional bifurcations and more complex motions can occur. It is still not clear how one might set useful information about turbulent flows from this approach, but for some complex flows like TaylorCouette flow the approach has been very successful.

14 See Marsden and Hughes [1994]; M. Golubitsky, I. Stewart and D. Schaeffer Symmetry and Groups in Bifurcation Theory, Vol. II, Springer-Verlag [1988]; and P. Chossat and G. Iooss, The Taylor-Couette Problem, Springer Applied Math. Sciences Series [1991] for further information and references.

This is page 101 Printer: Opaque this

3
Gas Flow in One Dimension

In this chapter we discuss compressible flow in one dimension. In the first section we develop the geometry of characteristics and in the second we introduce the notion of a weak solution and the entropy condition for shocks. In the third section we discuss the Riemann problem, that is, a flow problem with particular discontinuous initial data. A general construction, due to Glimm, which uses the solution of Riemann problems to produce solutions of arbitrary problems, is then presented. This construction is the basis of both some existence proofs and some methods of numerical computation in gas dynamics. In the final section we generalize the discussion to the flow of a gas that allows chemical energy release, such as occurs in combustion.

3.1

Characteristics

One-dimensional isentropic flow with an equation of state p = p(ρ) is described by the following equations derived in §1.1:  px  ut + uux = − , ρ (3.1.1)  ρt + uρx + ρux = 0. Define c = p (ρ) (assuming p (ρ) > 0) and note that px = p (ρ)ρx by the chain rule. Thus, the first equation in (3.1.1) becomes ut + uux = −c2 ρx . ρ

102

3 Gas Flow in One Dimension

Before proceeding, let us explain why c is called the sound speed . Sound can be viewed as a small disturbance in an otherwise motionless gas, due to density changes. Let u=u and ρ = ρ0 + ρ ,

where u and ρ are small and ρ0 is constant. Correspondingly, neglecting products of small quantities, we have p = p0 + ∂p ρ = p0 + c2 ρ ∂ρ and c2 = ∂p (ρ0 ) = constant. ∂ρ

Substitution of these expressions in (3.1.1) gives ut + u ux = −c2 ρx ρ and ρt + uρx + ρux = 0.

To first order (i.e., neglecting again squares of small quantities) these become the linearized equations ρ0 ut = −c2 ρx and ρt + ρ0 ux = 0.

Eliminating u by differentiating the first equation in x and the second in t, we get ρtt = c2 ρxx , which is the wave equation. A basic fact about this equation is that the general solution is ϕ(x + ct) + ψ(x − ct), that is, small disturbances (“sound waves”) propagate with speed c. If we define u ρ A= , c2 /ρ u then (3.1.1) becomes ρ u +A t ρ u

= 0. x (3.1.2)

Notice that the eigenvalues of A are u + c and u − c. Equation (3.1.2) is a special case of a first-order quasilinear hyperbolic system in one dimension, that is, a special case of the system ut + A(x, t, u)ux = B(x, t, u), (3.1.3)

where u = u(x, t) is an n-component vector function of x and t, and A is an n × n matrix function of x, t, and u that has n distinct real eigenvalues and hence n linearly independent eigenvectors (i.e., as a real matrix, A is diagonalizable). The nonlinearity of (3.1.3) is reflected in the dependence of A or B on the unknown vector u. We shall now examine systems of the form (3.1.3).

3.1 Characteristics

103

Example 1

Consider the single linear equation vt − vx = 0. (3.1.4)

Given initial data v(x, 0) = f (x), the solution is v(x, t) = f (x + t). Thus, v is a wave, with shape f , traveling with unit speed to the left on the x-axis. On the lines x + t = constant, v is constant. Thus, we can think of information from the initial data as being propagated along these lines. These lines are called the characteristics of equation (3.1.4). If initial data are given on a curve C that is transverse to the characteristics (i.e., nowhere tangent to them), then (3.1.4) is solved by setting v(x0 , t0 ) = the value of the initial data on the curve C at the point where C intersects the characteristic through (x0 , t0 ). See Figure 3.1.1. If initial data are given on t (x0 , t0) C

x

Figure 3.1.1. Characteristics intersecting the curve C.

a curve not transverse to the characteristics, such as on a characteristic itself, then the equation does not necessarily have a solution. Example 2 Consider the general single linear homogeneous equation vt + a(x, t)vx = 0. (3.1.5)

This time the characteristics are the curves t = t(s), x = x(s), satisfying dt = 1, ds dx = a(x, t). ds (3.1.6)

If a(x, t) is smooth, these curves exist (at least locally) and never intersect without being coincident, by the existence and uniqueness theorem for

104

3 Gas Flow in One Dimension

ordinary differential equations. Let us verify that, as in the first example, solutions are constant along characteristics. Indeed, if x = x(s), t = t(s) is a characteristic and v satisfies (3.1.5), then by the chain rule, d dx dt v(x(s), t(s)) = vx + vt = vt + a(x, t)vx = 0. ds ds ds Thus, if initial data is given on a curve C that is transverse to the characteristics, we can assert as before that v(x0 , t0 ) = the initial data at point P in Figure 3.1.2. If the data are given on a curve that is not transverse to the characteristics, a solution does not necessarily exist. t characteristics (x0,t0)

P x

Figure 3.1.2. Characteristics of a variable coefficient linear equation are nonintersecting curves.

Example 3 equation

Consider the generalization of (3.1.5) to the inhomogeneous vt + a(x, t)vx = b(x, t), (3.1.7)

and define the characteristics by the same equations (3.1.6). This time d v(x(s), t(s)) = b(x, t), ds and so, as the right-hand side is known, s v(x(s), t(s)) = v(x(s1 ), t(s1 )) + s1 b(x(α), t(α)) dα,

which again determines v off any curve C if initial data is given on C and C is transverse to the characteristics.

3.1 Characteristics

105

Example 4

As a nonlinear example, consider the equation ut + uux = 0. (3.1.8)

Let us search for curves x(s), t(s) along which u is constant. By the chain rule, d dx dt u(x(s), t(s)) = ux + ut . ds ds ds Thus, we should choose dt = 1, ds dx = u. ds

Thus, the characteristics depend on u. Along the curve defined by dt/ds = 1 and dx/ds = u, we have du/ds = 0; that is, u is constant. The situation is similar to the one in the linear case, except for the following crucial fact: characteristics can now intersect. This is easy to see directly; if u = constant on a characteristic, then dt/ds = 1, dx/ds = u is a straight line, and such lines issuing from different points can indeed intersect (see Figure 3.1.3). It is also easy to see in principle why the nonlinear case differs from the linear case. In the linear case, the characteristics are determined by two ordinary differential equations dt = 1, ds dx = a(x, t). ds

If a(x, t) is reasonably smooth, the solution of these equations through a given point (x0 , t0 ) is unique, and thus the characteristics cannot intersect. In the nonlinear case, the characteristics are determined by three equations (in the special case above, dt/ds = 1, dx/ds = u, du/ds = 0). The solution of these equations through (x0 , t0 , u0 ) is still unique, but the characteristics in the (x, t) plane are the projections of the three-dimensional solutions of the system on that plane, and thus they may indeed intersect. When such intersection occurs, our method of solution by following characteristics breaks down and the solution is no longer uniquely determined. We shall resolve this problem in the next section. Now we return to the system (3.1.3) and define its characteristics. Setting B = 0 for simplicity, the rate of change of u along a curve (x(s), t(s)) is du dt dx dt dx = ut + ux = −A(x, t, u) + ux . ds ds ds ds ds (3.1.9)

We define a characteristic to be a curve with the following property: if data are given on that curve, the differential equation does not determine the solution at any point not on the characteristic. (The data may also fail to be consistent, but we do not make this fact a part of the definition.)

106

3 Gas Flow in One Dimension

t

(x0,t0)

characteristics

(x1,0)

(x2,0)

x

Figure 3.1.3. The characteristics of a nonlinear equation may intersect.

One can readily check that all the examples of characteristics we have given meet this definition. If the characteristic is not parallel to the x-axis, this means we cannot determine ux from knowledge of u on the curve. From (3.1.9) this will happen if −A(x, t, u) is a singular matrix, that is, if dt = 1 and ds dx = λ(x, t, u), ds (3.1.10) dt dx + I ds ds

where λ(x, t, u) is an eigenvalue of A(x, t, u). Notice that there are n characteristics through each point and that the characteristics depend on u. If a change of variables is made on the dependent variable u that has a nonzero Jacobian, then A is changed by a similarity transformation, and thus the eigenvalues and hence the characteristics are unchanged. If a change of variables is made on the independent variables (x, t) (with nonzero Jacobian), then this change maps characteristics onto characteristics of the transformed problem. If C is a curve transverse to all the characteristics, then −A dt dx + ds ds

will be invertible, and from (3.1.9), ux = −A dx dt + ds ds
−1

us ,

3.1 Characteristics

107

so we can propagate u off C. Thus, it is reasonable to expect (and indeed it is true) that initial data on C uniquely determine a solution near C. (We have to say “near” because, as in Example 4, characteristics for the same eigenvalue can intersect.) In this general case we cannot expect u to be constant along characteristics. However, we might seek functions f1 , . . . , fn of u that are constant along the characteristics associated with the eigenvalues λ1 , . . . , λn . Indeed, suppose we can find a function f associated with an eigenvalue λ with the property that AT ∂f ∂f =λ , ∂u ∂u (3.1.11)

that is, ∂f /∂u is an eigenvector of the transpose of A. We claim that f is constant along the corresponding characteristic. Such functions f are called the Riemann invariants of the equation. To prove our contention, write out (3.1.11) in components: n Aki k=1 ∂f ∂f =λ , ∂uk ∂ui

(3.1.12)

where u = (u1 , . . . , un ). Now differentiate f along a curve satisfying (3.1.10): ∂f = ∂s = k=1 n n

k=1 n

∂f ∂uk = ∂uk ∂s

n

∂f ∂uk ∂ui − +λ Aki ∂uk ∂x ∂x i=1 n k=1 n

∂f ∂uk

∂uk ∂uk dx + ∂t ∂x ds

= k=1 − i=1 Aki

∂f ∂ui + ∂uk ∂x

n

λ k=1 ∂f ∂uk . ∂uk ∂x

This is zero by virtue of (3.1.12). Thus, if n functions f = (f1 , . . . , fn ) are found that are constant along characteristics, we may invert them to express u in terms of the f s and hope in this way to be able to determine explicitly the characteristics in terms of the initial data alone. An example will be given later to show that this can sometimes be done. Let us now see these ideas work for the equations of gas dynamics (3.1.1), written in the form (3.1.2) with u = (ρ, u). Because the eigenvalues of A are u ± c, the characteristics are the curves C+ : dx =u+c dt and C− : dx = u − c. dt (3.1.13)

To find the Riemann invariants we use (3.1.11). We seek eigenvectors of AT with eigenvalues u ± c; that is, vectors (h± , k± ) such that u ρ c2 /ρ u h± k± = (u ± c) h± k± .

108

3 Gas Flow in One Dimension

This is easy to do: one finds h± k± = ±c/ρ 1 .

Thus, the Riemann invariants are found by inspection to be Γ± = u ± that is, ∂Γ± = h± , ∂ρ ∂Γ± = k± , ∂u c(ρ) dρ, ρ (3.1.14)

as required by (3.1.11). Thus, Γ± are constant along the ± characteristics, respectively. The fact that Γ+ (resp. Γ− ) is constant on C+ (resp. C− ) does not in itself enable us to integrate the characteristic equations. The trouble is that Γ+ need not be constant on C− , and Γ− need not be constant on C+ ; if they were, then by inverting (3.1.14), u and ρ would be constant on characteristics and the characteristics would be straight lines. If, by some device, the characteristics can be found, then the equations can be solved as follows: From (x, t) follow the ± characteristics back to the curve on which initial data are prescribed to determine Γ± via (3.1.14). Using these values, solve (3.1.14) for (u, ρ) and the result will be the value of (u, ρ) at (x, t). Of course, as in Example 4 there may be difficulties with two C+ characteristics crossing. There may also be boundary conditions that have to be taken into account. We now discuss a different method for obtaining the Riemann invariants. Consider the change of variables u = u(x, t), ρ = ρ(x, t) from the (x, t) plane to the (ρ, u) plane and assume it is invertible (with nonzero Jacobian) so x and t may be regarded as functions of ρ and u. This is called the hodograph transformation. Because the transformations (x, t) → (ρ, u) and (ρ, u) → (x, t) are inverses, their Jacobian matrices are inverses: ρx ux ρt ut = xρ tρ xu tu
−1

=

1 J

tu −tρ

−xu xρ

,

where J = xρ tu − xu tρ = 0. Substitution into (3.1.1) yields a linear hyperbolic system with independent variables (ρ, u) and dependent variables (x, t) as follows: xρ − utρ + c2 tu = 0, ρ ρtρ + xu − utu = 0,

3.1 Characteristics

109

that is, x t + ρ 1 0

u/ρ ρ−1

0 c2 /ρ 1 −u

x t

= 0. u The eigenvalues of the coefficient matrix 1 0 u/ρ ρ−1 0 c2 /ρ 1 −u = 1 ρ u c2 − u 2 1 −u

are easily computed to be λ = ±c/ρ . Thus, the characteristics in the (ρ, u) plane are du c =± , dρ ρ that is, u± c(ρ) dρ = constant. ρ

Thus, we recover the fact that Γ± = constant defines the characteristics; that is, we recover the Riemann invariants obtained earlier. As we saw earlier, it may be difficult to determine the characteristics in a general flow. However, in special circumstances, the characteristics are straight lines. We now introduce some terminology relevant to this situation. A state of a gas is a pair of values (ρ, u). A constant state is a region in the (x, t) plane in which ρ and u are constant. A simple wave is a region in the (x, t) plane in which one of Γ+ , Γ− is constant. If it is Γ+ that is constant, we refer to the region as a Γ+ simple wave. A Γ− simple wave is defined analogously. In a Γ+ simple wave, both Γ− and Γ+ are constant along C− characteristics by definition of Γ− . From (3.1.14), u and c are constant too, and so C− is a straight line. Similarly, in a Γ− simple wave, the C+ characteristics are straight lines. Consider a constant state S that is bounded by a smooth curve σ in the (x, t) plane, as shown in Figure 3.1.4. (We implicitly assume that S is the largest region on which (ρ, u) are constant, that is, the region across σ exterior to S is not a constant state.) Because we are assuming c = 0, then the C+ and C− characteristics are not parallel at any point. Therefore, at any point P of σ, one of C+ or C− must cross σ. Suppose that C+ crosses σ. Now Γ+ is constant along each C+ characteristic and is constant throughout S. Therefore, assuming u and ρ are smooth enough, Γ+ will be constant on a region exterior to S. Thus there is a portion of the (x, t) plane adjoining S that is a Γ+ (or Γ− ) simple wave. If at some point along σ, the other characteristic C− was not equal to σ, it would cross it (again assuming smoothness) and so Γ− would be constant on some region exterior to S. That region would thus be a constant state. Therefore, the boundary of a constant state is composed of straight lines.

110

3 Gas Flow in One Dimension

σ C+ P S

constant state

Figure 3.1.4. Characteristics crossing the boundary of a constant state at P .

We will now give an example that shows the usefulness of the preceding observations. Consider an infinitely long tube extending along the x-axis. We assume that gas fills the half-tube x > 0 with a piston situated at x = 0, and that the gas is in a motionless state with a constant density ρ0 (Figure 3.1.5). At t = 0, we start pulling the piston to the left so that the tube (cut-out)

piston

gas in constant state u = 0, ρ = ρ0 x=0 x>0

Figure 3.1.5. Gas in a tube with a piston at t = 0.

piston follows a path x = x(t) for t > 0 in the (x, t) plane (Figure 3.1.6). As a result the gas is set into motion. Because the gas is originally motionless with uniform density ρ0 , there is a constant state I in the (x, t) plane where the C± characteristics have constant slope ±c(ρ0 ). The Riemann invariants in I are Γ0 = ± ± ρ0 a

c(s) ds = constant. s

The boundary of I is a straight C+ characteristic that must emanate from the origin. This is readily seen from the fact that two distinct C+ characteristics in a Γ− wave cannot intersect at t > 0. If they did, they would coincide. The C− characteristics of I penetrate into the adjoining region II, and they can be traced at most as far as the piston path. Hence,

3.1 Characteristics

111

particle path piston path x = x(t)

t C+

II I

slope = c(ρ0)

0 C–

constant state u = 0, ρ = ρ0

x

Figure 3.1.6. Dynamics of the gas in the piston.

the region must be a Γ− simple wave (part of it may also be a constant state). The C+ characteristics in the region II are all straight lines. In particular, if the piston is withdrawn with a constant velocity U (which is negative), that is, x(t) = U t, then we claim that the density of the gas is constant along the piston path (Figure 3.1.7). To prove this fact, pick any left boundary C+

constant state St = III piston path x = Ut C–

II

right boundary

I = Sr constant state u = 0, ρ = ρ0 particle path x

Figure 3.1.7. Piston pulled with constant (negative velocity).

point B on the piston path. The gas particle at B moves with velocity U .

112

3 Gas Flow in One Dimension

Because Γ− (B) = Γ0 , we have − ρ(B) a

c(ρ) dρ = U − Γ0 = constant. − ρ

(3.1.15)

This is possible only if ρ(B) is constant along the piston path, because c(ρ) and ρ are positive. Hence, the Riemann invariant Γ+ , which is equal to U + (U − Γ0 ) on the piston path, is a constant independent of any − C+ characteristics issuing from the piston path. Thus, there is a constant state III adjacent to the piston path. The other boundary of III must be a straight C+ characteristic, which issues from the origin. The constant states I and III enclose a region II which must be a Γ− simple wave. Because the boundary of II consists of straight C+ characteristics issuing from the origin, all C+ characteristics of II intersect at the origin. A simple wave whose straight characteristics intersect at a point on the x-axis is called a centered rarefaction wave. We justify this terminology by noting the following fact: From (3.1.15), ρ(B) ρ0

c(ρ) dρ = ρ

ρ(B) a

c(ρ) dρ + Γ0 = U < 0. − ρ

Because c(ρ) and ρ are positive, we must have ρ(B) < ρ0 (i.e., in Figure 3.1.7 the gas is rarefied when the gas state passes through II from I to III ). We can carry out the explicit construction of the solution to the piston problem at any point (x, t) in the special case of “γ law gas” whose equation of state is given by p = Aργ , A, γ = constants, γ > 1. The sound speed in such a gas is c= dp = dρ Aγργ−1 = γp , ρ (3.1.16)

and the Riemann invariants are Γ± = u ± c(ρ) 2 dρ = u ± c. ρ γ−1

We call the constant state I the right state Sr . In this state ur = 0, ρr = ρ0 , pr = Aργ , 0 and cr = γpr . ρ0

We call state III the left state Sl and in it (u, ρ, p, c) have constant values, say (ul , ρl , pl , cl ). Obviously, ul = U . From Γ− = U − 2 2 2 cl = ur − cr = − cr , γ−1 γ−1 γ−1

3.1 Characteristics

113

we can determine cl and hence the entire left state from (3.1.16). Because we know cl , the left boundary of the rarefaction wave is given by x = (U + cl )t. Finally, if (x, t) is a point inside the rarefaction, then u+c= x , t (3.1.17)

since the C+ characteristic through (x, t) is a straight line issuing from the origin. On the other hand, the rarefaction is a Γ− simple wave, and Γ− = u − 2 2 c = ur − cr . γ−1 γ+1 (3.1.18)

Relations (3.1.17) and (3.1.18) yield u(x, t) = c(x, t) = 2 x − cr , γ+1 t γ−1 x 2 + cr . γ+1 t γ+1

We now use (3.1.16) and p = Aργ to obtain ρ(x, t) and p(x, t). Thus we have constructed the gas flow at every point (x, t). Two constant states, S1 and S2 , are said to be connected by a Γ− rarefaction wave if the configuration shown in Figure 3.1.8 is a solution of the differential equations for gas flow. t S2

S1

x
Figure 3.1.8. Two constant states connected by a rarefaction wave.

Remember that in a Γ− simple wave, the C+ characteristics are straight lines. The configuration in Figure 3.1.8 depicts two constant states, S1 and S2 , connected by a Γ− rarefaction wave. We ask the question: Given a state S1 , that is, given values ρ1 and u1 , what are the possible states S2 that can be connected to S1 by means of a centered rarefaction wave? We are only interested in states S2 whose density ρ2 is less than ρ1 ; this corresponds to the situation created when a piston moves to the left. For this discussion, we assume the relation p = Aργ . Choose p2 such that 0 < p2 < p1 . We claim that p2 will uniquely determine the state S2 that can be connected to S1 by a centered rarefaction

114

3 Gas Flow in One Dimension

wave. Given p2 , we find ρ2 from the relation p2 = Aργ . Given S1 , we know 2 ρ1 and u1 and thus can compute Γ− = u 1 − 2 c(ρ1 ), γ−1

which is constant in S1 . Because Γ− is constant throughout the centered rarefaction, we have u2 − 2 2 c(ρ2 ) = u1 − c(ρ1 ). γ−1 γ−1

This determines u2 and hence the constant state S2 . How far can these C+ characteristics in Figure 3.1.8 fan out? The answer is seen from the relation Γ− = u − 2 c(ρ); γ−1

that is, when ρ is zero we must stop. The methods developed so far cannot deal with the case of a piston being pushed in. In such a case, our C+ characteristics would soon collide, as shown in Figure 3.1.9. t piston path

collision S1 x

Figure 3.1.9. The characteristics can collide if the piston is pushed in.

When the characteristics collide, the solutions become multiple valued and thus, as presented so far, meaningless. To deal with this situation we will modify our approach and introduce the concept of a “weak solution” in the next section.

3.2 Shocks

115

3.2

Shocks

To deal with the crossing of characteristics and the formation of shocks, we introduce some concepts from thermodynamics and the notion of a weak solution. We begin with a brief review of relevant parts of thermodynamics for ideal flow.1 We assume there is a specific internal energy function (“specific” means “per unit mass”). For example, as we saw in §1.1, for ideal gases, p 1 = . ρ γ−1 The total energy per unit volume is e = 1 ρu2 + ρ . 2 Assuming that heat does not enter the fluid domain from its boundaries, the total energy can only be reduced when the fluid does work. The work done by a fluid volume W per unit time is −
∂w

pu · n dA

where n is the outward unit normal to ∂W and dA is the area element. This must equal the rate of change of energy in W : ∂t
W

e dV = −
∂W

pu · n dA = −
W

div(pu) dV.

By the transport theorem in §1.1, this integral form is equivalent to the differential equation ∂t e + div(eu) + div(pu) = 0. that is, ∂t e + div((e + p)u) = 0, In one dimension this becomes ∂t e + ((e + p)u)x = 0. (3.2.1)

Let us check that (3.2.1) indeed holds for an ideal gas. Here we choose p = Aργ and = p 1 ; ργ−1

1 For an “axiomatic” treatment, see Marsden and Hughes [1994]. See also L. Malvern [1969] Introduction to the Mechanics of a Continuous Medium, Prentice-Hall.

116

3 Gas Flow in One Dimension

then (3.2.1) reads ∂t
2 1 2ρ u

+

p γ−1

+

2 1 2ρ u

+

p +p u γ−1

= 0. x Substituting p = Aργ and the equation of motion ut + uux = −px /ρ, we see that this equation is indeed satisfied, as long as u and ρ are smooth. The energy equation (3.2.1) is called the first law of thermodynamics. We regard it as one of the basic equations. Our system, thus enlarged, becomes  ρt + uρx + ρux = 0,    px ut + uux + = 0, (3.2.2) ρ    et + ((e + p)u)x = 0. When shocks form, discontinuities develop and so the meaning of equations (3.2.2) has to be carefully considered. Later we shall examine in detail why it is desirable to work with (3.2.2). For instance, for a gas (for which e = 1 ρ u2 + p/(γ − 1)) we shall see that (3.2.2) can be meaningful even 2 when p = Aργ breaks down. Physically, it is reasonable to concentrate on the system (3.2.2) because it expresses conservation of mass, momentum, and energy. We remark that the system (3.2.2) has no Riemann invariants in general. The second law of thermodynamics in the form of a mathematically reasonable postulate about shock waves will be introduced later. In thermodynamics, the second law asserts the existence of an “entropy” η that has the property Dη/Dt ≥ 0; η increases when energy is converted from kinetic to internal energy. In works on thermodynamics, η is related to the other variables, for example, to specific heats. For an ideal gas, one can show that η = [ constant ] log(pρ−γ ). If η is increasing, we cannot have p = Aργ , where A is a constant. When shocks form, η increases and p = Aργ must be abandoned. However, a suitable formulation of equations (3.2.2) will continue to make sense. We turn to this formulation next. We now introduce the concept of a “weak solution” that will allow discontinuous solutions. The basic idea behind weak solutions is to go back to the integral formulation of the equations; this may remain valid when there is not enough differentiability to justify the differential form.2
2 For a general mathematical reference, see P. D. Lax [1973] Hyperbolic Systems of Conservation Laws and the Mathematical Theory of Shock Waves, Conference Board of the Mathematical Sciences Regional Conference Series in Applied Mathematics, No. 11. Society for Industrial and Applied Mathematics, Philadelphia, Pa., v+48 pp.

3.2 Shocks

117

Let us begin by considering general nonlinear general equations in conservation form, that is, of the form ut + (f (u))x = 0. (3.2.3)

Example The equation ut + uux = 0 was discussed in the previous section. It may be written in the conservation form ut + 1 u2 x = 0. However, 2 it can also be written as
1 2 2u t

+

1 3 3u x

= 0.

The results one gets may depend on which conservation form is chosen. Usually, physical considerations will dictate which one is appropriate. Let F = (f (u), u) and note that (3.2.3) is equivalent to Div F = 0, where Div(f1 , f2 ) = (f1 )x + (f2 )t is the space-time divergence. Let ϕ be a smooth function with compact support in the (x, t) plane (i.e., ϕ vanishes outside a compact set). Then (3.2.3) is equivalent to ϕ · Div F dx dt = 0, for all such test functions ϕ. Integrating by parts (which is justified because ϕ has compact support), Grad ϕ · F dx dt = 0. (3.2.4)

where Grad ϕ = (ϕx , ϕt ) is the space-time gradient of ϕ. If u is smooth (3.2.3) and (3.2.4) are equivalent, but if u is not smooth, (3.2.4) can make sense even when (3.2.3) does not. We define a weak solution of (3.2.3) to be a function u that satisfies (3.2.4) for all smooth ϕ with compact support. A slightly different formulation of (3.2.4) may be used to take into account the initial conditions, say u(x, 0) = q(x). Namely, we can replace the integration over all of space-time by integration over the half-space t ≥ 0. Then during integration by parts we pick up a boundary term: Grad ϕ · F dx dt + t≥0 x ∞

ϕ(x, 0) q (x) dx = 0.
−∞

(3.2.5)

If ϕ has support away from the x-axis (i.e., ϕ = 0 on the x-axis) then (3.2.5) reduces to (3.2.4). So far we have the differential form of equation (3.2.3) and the weak form (3.2.4). We also have an integral form corresponding to the integral form of

118

3 Gas Flow in One Dimension

the equations of motion discussed in §1.1. Here it expresses “conservation of u.” Indeed, consider an interval W = [a, b] on the x-axis; from (3.2.3), d dt Thus, d dt b a b

u dx =
W W

ut dx = −
W

(f (u))x dx.

u dx = −f (u) . a ∞ −∞

In particular, note that if u → 0 at ±∞, then literal sense: d ∞ u dx = 0. dt −∞

u dx is conserved in the

A weak solution does not have to be differentiable. Note, however, that a function that satisfies the integral form of the equation of motion does not have to be differentiable either, and the integral forms in fluid mechanics are in fact the basic laws we are working with, while weak solutions are merely convenient mathematical tools. A natural question to ask is: Does a weak solution necessarily satisfy the integral form of the equation? The answer is yes, and, therefore, weak solutions are in fact the objects we are looking for. (Note that the answer is, in general, affirmative only if the conserved quantities are the same in the integral forms and in the conservation forms of the equations.) The proof of the fact that weak solutions satisfy the integral form of the equation is reasonably straightforward. Let ψD be the characteristic function of the domain D over which the integral form of the equation is used; i.e, ψD (x, t) = 1 if (x, t) ∈ D, while ψD (x, t) = 0 if (x, t) is not in D. Substitute ψD instead of ϕ into equation (3.2.4). Some straightforward manipulation of improper functions will yield the desired integral form. However, ψD is not an acceptable test function ϕ because it is not smooth. Despite this, one can find sequences of acceptable test functions ϕn (x, t) such that ϕn (x, t) div F dx dt → ψD (x, t) div F dx dt,

as n → ∞ for all piecewise smooth F. During this process, one must account carefully for lines on which F is not smooth. We omit the details. We now investigate properties of weak solutions of the conservation law ut + f (u)x = 0 near a jump discontinuity. Let u be a weak solution with a jump discontinuity across a smooth curve Σ in the (x, t) plane. Let ϕ be a smooth function vanishing outside the region S shown in Figure 3.2.1.

3.2 Shocks

119

t S Σ S1 S2 n

x

Figure 3.2.1. A weak solution may have a jump discontinuity across the curve Σ.

Write S = S1 ∪ S2 as shown in the figure. Then 0= =
S1

Grad ϕ · F dx dt Grad ϕ · F dx dt +
S2

Grad ϕ · F dx dt.

(3.2.6)

Let us assume that u is smooth in each of regions S1 and S2 . Thus grad ϕ · F dx dt =
S1 S1

div(ϕF) dx dt −
S1

ϕ Div F dx dt ϕ Div F dx dt.

=
Σ

ϕF · n ds dt −
S1

In region S1 where u is smooth, Div F = 0, and therefore Grad ϕ · F dx dt =
S1 Σ

ϕF1 · n ds,

where in the integral over Σ, F1 means u is evaluated by taking the limit from region S1 . Similarly, one has Grad ϕ · F dx dt = −
S2 Σ

ϕF2 · n ds.

The minus sign occurs because the outward normal n for S1 is the inward normal for S2 . Substitution into (3.2.6) yields ϕ(F1 − F2 ) · n ds = 0.
Σ

120

3 Gas Flow in One Dimension

Because this is valid for all ϕ, [F · n] = 0 on Σ, (3.2.7) where [F · n] = F1 · n − F2 · n denotes the jump in F · n across Σ. Let Σ be parametrized by x = x(t) so s = dx/dt is the speed of the discontinuity. Then (1, −s) n= √ s2 + 1 thus, (3.2.7) becomes −s[u] + [f (u)] = 0 on Σ (3.2.8) and F = (f (u), u);

where again [ ] denotes the jump in a quantity produced when (x, t) crosses from S1 to S2 . Equation (3.2.8) is the constraint that the assumed weak form of the equations (3.2.4) imposes on the values of u on both sides of a discontinuity. Also, u dx is conserved if the differential form of the equation is satisfied whenever u is smooth, and (3.2.8) is satisfied across any discontinuity. It is, of course, assumed that all discontinuities are jump discontinuities. A function u that satisfies the differential equations where possible and (3.2.8) across a jump discontinuity satisfies the integral form and the weak form of the equations. What we have just done works just as well for systems of conservation laws, that is, equations of the form (ui )t + (fi (u1 , . . . , un ))x = 0, i = 1, . . . , n

for a vector unknown u = (u1 , . . . , un ). Setting Fi = (fi (u), ui ), these may be rewritten as Div Fi = 0, i = 1, . . . , n. The comments on the single equation case then carry over verbatim, equation by equation, to the case of systems. In particular, each component ui represents a conserved quantity of the system in the sense spelled out earlier. As an example of a system of conservation laws, consider the equations of isentropic gas flow: px ρt + (ρu)x = 0, ut + uux + = 0. ρ They can be written in conservation form if one writes the second equation in terms of the momentum m = ρu. One gets  ρt + mx = 0,   (3.2.9) 2 m  mt + = 0.  +p ρ x

3.2 Shocks

121

In this system we assume p = p(ρ). We can drop this assumption p = p(ρ) if we add the energy equation et + ((e + p)u)x = 0, where e = ρ + 1 ρu2 and 2 is given, for example, by = p 1 . ργ−1

In summary, using the energy equation, we have the system of conservation laws  ρt + mx = 0,       2  m  mt + +p = 0, (3.2.10) ρ x      m  et + (e + p) = 0.   ρ x The system (3.2.9) together with p = Aργ will not have the same weak solutions as (3.2.10) together with e = 1 ρu2 + ρ , 2 = p 1 . ργ−1

It is believed on physical grounds that the condition of conservation of energy is more fundamental than p = Aργ and that, indeed, A may depend on the entropy and thus may not be constant. Therefore, we adopt the system (3.2.10) as our basic system of conservation laws. There are problems, however (e.g., in the theory of water waves), where systems such as the 2 × 2 system with p = Aργ are appropriate. Physical considerations dictate which quantities must be conserved. From the jump conditions (3.2.8) applied to the system of conservation laws (3.2.10) we find that the jump relations across a discontinuity Σ in the (x, t) plane with velocity dx/dt = s are given as follows:  s[ρ] = [m],      2 m (3.2.11) s[m] = +p ,  ρ     s[e] = [(e + p)u]. These are called the Rankine–Hugoniot relations. The first two equations in (3.2.11) are called the mechanical jump relations. The equations of fluid mechanics, like those of classical particle mechanics, are Galilean invariant. Therefore, it is legitimate to transform variables to a coordinate system moving with uniform velocity. It is convenient to

122

3 Gas Flow in One Dimension

t

shock

back(1)

front(0) x

Figure 3.2.2. The shock in moving coordinates.

pick a coordinate system whose velocity at some fixed time t0 , say t0 = 0, equals that of the discontinuity and such that the discontinuity is at x = 0 at t = 0. See Figure 3.2.2. In this coordinate system, the Rankine–Hugoniot relations at the origin become  ρ 0 u 0 = ρ1 u 1 ,   (3.2.12) ρ0 u2 + p0 = ρ1 u2 + p1 , 0 1   (e0 + p0 )u0 = (e1 + p1 )u1 , where the subscripts indicate the different sides (0) or (1) of the discontinuity. Let M = ρ0 u0 = ρ1 u1 . If M = 0 we call the discontinuity a contact discontinuity or a slip line. Because u0 = u1 = 0, these discontinuities move with the fluid. From (3.2.12)2 , p0 = p1 , but in general ρ0 = ρ1 . If M = 0, we call the discontinuity a shock . Because u0 = 0, u1 = 0, gas is crossing the shock, or, equivalently, the shock, is moving through the gas. The side that consists of gas that has not crossed the shock is called the front of the shock, and we identify it by a subscript 0. The other side, denoted by a subscript 1, is called the back . Some simple algebraic identities for a shock may be derived from the Rankine–Hugoniot relations. From (3.2.12)1 and (3.2.12)2 , M u0 + p0 = M u1 + p1 , i .e., M = − p0 − p1 . u0 − u1 (3.2.13)

Substituting u0 = M/ρ0 = M τ0 and u1 = M/ρ1 = M τ1 , where τ = 1/ρ is the specific volume, into (3.2.13) gives M2 = − p0 − p1 . τ0 − τ1 (3.2.14)

3.2 Shocks

123

If we write M 2 = M · M = (ρ0 u0 ) (ρ1 u1 ) and τ = 1/ρ , (3.2.14) becomes u0 u1 = p 0 − p1 . ρ 0 − ρ1 (3.2.15)

The identities (3.2.13), (3.2.14), and (3.2.15) are consequences of the mechanical jump relations only. To bring in the energy, combine (3.2.12)1 with (3.2.12)3 to give e0 τ0 − e1 τ1 = p1 τ1 − p0 τ0 . However, e0 τ0 − e1 τ1 =
2 1 2 ρ0 u 0

+ ρ0

0

τ0 −

2 1 2 ρ1 u 1

+ ρ1

1

τ1

= 1 (u0 − u1 )(u0 + u1 ) + 0 − 1 2 p0 − p 1 (M τ0 + M τ1 ) + ( 0 − 1 ) =− 2M (from (3.2.12)1 and (3.2.13)) p 0 − p1 = ( 0 − 1) − (τ0 + τ1 ). 2 p0 + p1 (3.2.16) (τ1 − τ0 ) = 0. 2 This relation is called the Hugoniot equation for the shock. Notice that it depends only on p and τ and not on u. Define the Hugoniot function with center (τ0 , p0 ) to be
1

Thus,



0

+

H(τ, p) = (τ, p) − (τ0 , p0 ) +

p0 + p (τ − τ0 ) 2

so that (3.2.16) may be written as H(τ1 , p1 ) = 0. For a “γ-law gas” for which = 1 pτ, γ−1

it is easily checked that H(τ, p) is given by 2µ2 H(τ, p) = (τ − µ2 τ0 ) p − (τ0 − µ2 τ ) p0 , where µ2 = γ−1 . γ+1

In this case, the curve H(τ, p) = 0 is the hyperbola shown in Figure 3.2.3. The Hugoniot equation states that (τ1 , p1 ) lies on this hyperbola. The hyperbola represents all possible states that can be connected to the state (τ0 , p0 ) through a shock. Notice from (3.2.14) that −M 2 is the slope of the

124

3 Gas Flow in One Dimension

p H(τ,p) = 0

(τ1, p1)

(τ0, p0) µ2τ0 –µ2p
0

τ

Figure 3.2.3. The Hugoniot curve H(τ, p) = 0.

line through (τ1 , p1 ) and (τ0 , p0 ) and that this then determines u0 and u1 γ via (3.2.12)1 . Because the curve pρ−γ = p0 ρ−γ , that is, pτ γ = p0 τ0 does 0 −γ not coincide with the curve H(τ, p) = 0, we will not have pρ = constant across a shock, in general. Similarly, if we had insisted on (3.2.9) and pρ−γ = constant, we would have obtained shocks for which energy conservation would have been violated. Further conditions must be imposed on weak solutions of conservation laws to select a unique, physically correct solution. The next example demonstrates nonuniqueness with ut + uux = 0. Example Consider the conservation law ut + u2 2 = 0. x (3.2.17)

Recall from §3.1 that its characteristics are straight lines and that u is constant along them. Consider the initial data u(x, 0) = 0, if x ≥ 0; 1, if x < 0.

The corresponding characteristics are shown in Figure 3.2.4(a). To keep the characteristics from crossing, we introduce a shock with propagation speed 1 2 u 1 s= 2 = [u] 2 (see equation (3.2.8)). We thus get a globally defined weak solution by letting u = 1 to the left (behind) of the shock and u = 0 to the right (in front) of the shock. See Figure 3.2.4(b).

3.2 Shocks

125

t

t

shock; s = 1/2

u=1 x (a) characteristics (b)

u=0 x

Figure 3.2.4. Characteristics for the initial data u = 0 if x ≥ 0 and u = 1 if x < 0.

Next consider the initial data u(x, 0) = 1, if x ≥ 0; 0, if x < 0.

The characteristics, shown in Figure 3.2.5(a), do not fill out the (x, t) plane. t t s = 1/2 u=1

u=0

u=0

x (a) u=0 t rarefaction fan u=1 (b)

x

x (c)
Figure 3.2.5. Characteristics for the initial data u = 1 if x ≥ 0 and u = 0 if x < 0.

126

3 Gas Flow in One Dimension

The figure shows two ways to find a globally defined weak solution. We shall introduce a condition that excludes the solution (b) of Figure 3.2.5. Consider the characteristics of our problem and consider a shock (i.e., a discontinuity that satisfies the jump condition −s[u] + 1 u2 = 0) . It 2 is easy to see that the following is true: Through every point on the path of the shock in the (x, t) plane one can draw two characteristics, one on each side of the shock; and either both of them can be traced back to the initial line or both of them can be traced upwards to the “future,” that is, toward larger t. In either case, the shock is needed to avoid having characteristics intersect and create a multivalued solution. We say that this shock separates the characteristics. (We shall work out an example more general than 3.2.17 in greater detail later.) A shock is said to obey the entropy condition if the two characteristics that intersect on each point can be traced backward to the initial line, as in Figure 3.2.4(b). A shock that does not obey the entropy condition is called a rarefaction shock . We shall allow only shocks that do obey the entropy condition and exclude solutions such as those in Figure 3.2.5(b). This restriction will make the weak solution of the problem unique. The reason for the name “entropy condition” will appear later. The entropy condition can be viewed as a causality condition: The shock is determined by the given data, and not by future events.3 We shall formulate shortly the general entropy condition for shocks in systems of conservation laws. Before doing so, we list some of the reasons rarefaction shocks are excluded in gas dynamics. 1. If a rarefaction shock is allowed, then the problem will not have a unique solution. 2. A solution that includes rarefaction shocks need not depend continuously on the initial data. Specifically, in the neighborhood of a rarefaction shock, we can specify u on the left, u on the right, and s; the only required relation is the jump condition; u is not constrained by the initial data. Thus, we can alter the solution without altering the initial data. 3. When we write the equations of gas dynamics in hyperbolic form, we neglect the viscosity. The hidden assumption is that the effect of viscosity should be small. As an example, we consider the viscous equation ut + u2 2 = ν uxx , x ν>0

(3.2.18)

3 See P. D. Lax [1973] Hyperbolic Systems of Conservation Laws and the Mathematical Theory of Shock Waves, Conference Board of the Mathematical Sciences Regional Conference Series in Applied Mathematics, No. 11. Society for Industrial and Applied Mathematics, Philadelphia, Pa., v+48 pp.

3.2 Shocks

127

corresponding to the conservation law (3.2.17). Let the solution of (3.2.18) be uν and let u0 be the solution of (3.2.17) with the same initial data. Then the hidden assumption is lim uν = u0 . ν→0+ In fact, one can show that limν→0+ uν is a weak solution of the inviscid equation (3.2.17), and one can also show that only solutions of (3.2.17) satisfying the entropy condition are the limits of the corresponding viscous solutions. We will omit the proof.4 4. We have seen that the quantity u(x, t) dx is conserved for a solution of the conservation law (3.2.17). One can show that the “energy,” that is, u2 (x, t) dx, cannot increase for a weak solution whose shocks satisfy the entropy condition. To make this plausible, consider a solution u(x, t); we define its variation at t to be var(u(x, t)) = sup n |u(xn , t) − u(xn−1 , t)|

where the sup is taken over all the possible partitions {x1 , . . . , xN } of the x-axis. (A partition is a finite division of the x-axis, −∞ < x1 < · · · < xN < ∞.) We assume that var(u(x, 0)) < ∞ initially. Then, because a smooth solution is propagated along characteristics from the initial data, var(u) is invariant for smooth solutions (see Figure 3.2.6). t shock satisfying the entropy condition characteristics t0 t

t0

x y2 0 y1 var (u(x,t0)) < var (u(x,0)) (a)

x 0 var(u(x,t0)) = var(u(x,0)) (b)

Figure 3.2.6. The variation is decreasing for a shock satisfying the entropy condition.

However, if u(x, t) has a shock satisfying the entropy condition, u(x, t) loses a part of its initial variation at any time when this shock is present. Therefore, the variation var(u) cannot increase. This fact implies that
4 See

E. Hopf, Comm. Pure Appl. Math. 4 [1950], 201.

128

3 Gas Flow in One Dimension

u2 (x, t) dx cannot increase. This is not the case with a rarefaction shock. A similar situation holds for the equations of gas dynamics. 5. In the case of the equations of gas dynamics, rarefaction shocks violate basic thermodynamic principles. We now formulate the entropy condition for systems of conservation laws. Consider an n-component system ut + (f (u))x = 0. Using the chain rule, we can write ut + A(u)ux = 0, where A(u) is the Jacobian matrix of f . We assume this is hyperbolic in the sense of the previous section so A(u) has n real eigenvalues λ1 , λ2 , . . . , λn . Correspondingly, we have n families of characteristics given by the curves dxi = λi , dt i = 1, . . . , n.

For example, in the one-dimensional case we have considered, n = 1 and λ1 = λ = u. In the two-component case of isentropic flow, we have two families of characteristics, C+ and C− , associated with λ1,2 = u ± c. In the 3 × 3 system of gas dynamics, λ1 = u + c, λ2 = u − c, λ3 = u, and we have three corresponding families, C+ , C− , and C0 , respectively. A shock is said to separate characteristics of a given family if i ii it satisfies the jump conditions; through every point of the trajectory of the shock in (x, t) plane one can draw two characteristics of the family, one on each side of the shock; and either both characteristics can be traced back to the initial line or they can both be traced upwards toward increasing t. (We shall see later that it is a property of the gas dynamic equations that each shock separates either the C+ or the C− family.)

iii

The entropy condition for systems of conservation laws is the following: A shock satisfies the entropy condition if, when it separates characteristics of one family, the characteristics on each side can be traced back to the initial data. We shall allow only shocks that satisfy the entropy condition. (Some authors call a discontinuity a shock only if the entropy condition is satisfied.) A shock for gas dynamics is called compressive if the pressure behind the shock is larger than the pressure in front of the shock. Thus, in a compressive shock the pressure is raised as a fluid particle crosses the shock. We next show that for a γ-law gas a shock is compressive if and only if it satisfies the entropy condition. To demonstrate this, we proceed in a number of steps.

3.2 Shocks

129

Step 1 The velocity of a shock is larger than the sound speed (i.e., is supersonic) on one side and is smaller than the sound speed (i.e., is subsonic) on the other. To see this, consider the Hugoniot relation e0 τ0 + p0 τ0 = e1 τ1 + p1 τ1 derived in the course of establishing equation (3.2.16). We can rewrite it as
1 2 2 u0

+

0

+ p 0 τ0 = 1 u 2 + 2 1 pτ , γ−1

1

+ p1 τ 1 .

Substituting µ2 = we get µ2 u2 + (1 − µ2 )c2 = µ2 u1 + (1 − µ2 )c2 . 0 0 1 Let c2 ∗ denote the common value of both sides so that (1 − µ2 )(u2 − c2 ) = u2 − c2 0 0 0 ∗ and (1 − µ2 )(u2 − c2 ) = u2 − c2 . 1 1 1 ∗ However, µ2 < 1 because γ > 1, and so |u0 | > c0 if and only if |u0 | > c∗ and |u1 | > c1 if and only if |u1 | > c∗ Rewrite the expression for c2 as follows: by definition, ∗ ρ1 c2 = ρ1 [µ2 u2 + (1 − µ2 )γp1 τ1 ] = ρ1 µ2 u2 + (1 − µ2 )γp1 . ∗ 1 1 Because (1 − µ2 )γ = 1 + µ2 , ρ1 c2 = ρ1 µ2 u2 + (1 + µ2 )p1 . ∗ 1 From (3.2.12)2 we get ρ1 c2 = µ2 P + p1 ∗ (3.2.20) (3.2.19) γ−1 , γ+1 = and c2 = γ p τ,

where P = ρ0 u2 + p0 = ρ1 u2 + p1 . Similarly, ρ0 c2 = µ2 P + p0 . Eliminating ∗ 0 1 P yields p1 − p0 c2 = . ∗ ρ1 − ρ0 From (3.2.15) we then get c2 = u0 u1 . (3.2.21) ∗ This is called the Prandtl relation. It follows that |u1 | > c∗ implies |u0 | < c∗ , and vice versa. Together with (3.2.20), this proves our contention in Step 1.

130

3 Gas Flow in One Dimension

Step 2 shock.

Determination of the family of characteristics separated by a

Consider a shock facing to the right, that is, its front is on its right, its back is on its left, and fluid crosses it from right to left. The characteristics for our system (3.2.10) are given by dx = u + c (the C+ characteristic) dt dx = u − c (the C− characteristic) dt and dx =u dt (the C0 characteristic)

The shock cannot separate C0 characteristics. Indeed, by our conventions, the velocity u0 in front of the shock is negative; by the jump relation ρ0 u0 = ρ1 u1 , we see that u1 is negative as well. Therefore, the configuration of the C0 characteristics is such that on the left (labeled with a (1)) the characteristics go to the future and on the right (labeled with a (0)) they go to the past (see Figure 3.2.7), and so the shock does not separate them. The right-facing shock cannot separate C− characteristics either. Indeed, because u0 < 0, we have u1 < 0 as well and u0 − c < 0, u1 − c < 0; thus, one has a picture qualitatively similar to that shown in Figure 3.2.7 for C0 characteristics. t shock

C0 characteristics

x

Figure 3.2.7. The appearance of C0 characteristics separated by a shock.

Thus, a right-facing shock can only separate C+ characteristics. (Similarly, a left-facing shock separates C− characteristics.) Step 3 A shock is compressive if and only if ρ0 < ρ1 . Indeed, a shock is compressive when p1 > p0 by definition. From the jump condition we get ρ0 u0 = ρ1 u1 , and thus u1 and u0 are both negative

3.2 Shocks

131

for a right-facing shock. From (3.2.19), c2 = ∗ p1 − p0 > 0; ρ1 − ρ0

therefore, p1 − p0 and ρ1 − ρ0 have the same sign, as required. A similar argument yields the same conclusion for a left-facing shock. Step 4 A noncompressive shock violates the entropy condition.

A shock is noncompressive when ρ1 < ρ0 , from Step 3. From the jump condition ρ0 u0 = ρ1 u1 we have u1 < u0 < 0. Subject to the Hugoniot constraint H(τ, p) = 0, and using implicit differentiation, one sees that c= γp ρ

is an increasing function of ρ. From ρ1 < ρ0 and u1 < u0 < 0, we get the inequality u1 + c(ρ1 ) < u0 + c(ρ0 ). From Step 1, it follows that |u1 | > c(ρ1 ), |u0 | < c(ρ0 ), u1 + c(ρ1 ) < 0, and u0 + c(ρ0 ) > 0.

Thus, the configuration of the slopes of the C+ characteristics have the form shown in Figure 3.2.8. This clearly violates the entropy condition. t shock

dx = u1 + c(ρ1) < 0 dt

dx = u0 + c(ρ0) > 0 dt x

Figure 3.2.8. The C+ characteristics for a noncompressive shock.

Step 5

A compressive shock obeys the entropy condition.

132

3 Gas Flow in One Dimension

Here ρ1 > ρ0 and u0 < u1 < 0. Thus, c(ρ1 ) > c(ρ0 ) and u1 + c(ρ1 ) > u0 + c(ρ0 ). By Step 1 and Prandtl’s relation we must have |u0 | > c(ρ0 ) and |u1 | < c(ρ1 ) so that u1 + c(ρ1 ) > 0 and u0 + c(ρ0 ) < 0. Thus, the characteristics have the appearance shown in Figure 3.2.9, and therefore the entropy condition is satisfied. t shock dx = u0 + c0 < 0 dt

dx = u1 + c1 > 0 dt

x

Figure 3.2.9. The appearance of the characteristics for a compressive shock.

We note again that if the front side were the left side, then the shock would separate the C− characteristics, but the same conclusion would be reached. This completes the proof of our original contention that compressive shocks for a γ–law gas obey the entropy condition. We remarked earlier that the thermodynamic entropy can be shown to be an increasing function of pρ−γ and that the entropy should increase across a shock. This is easy to see by consulting Figure 3.2.10. Indeed, the curve H(τ, p) = 0 through (τ0 , p0 ) is drawn as well as the curves pρ−γ = constant; one sees that pρ−γ increases as we move along H(τ, p) = 0 in the direction of increasing p. Thus, for a γ-law gas, we see that the geometric entropy condition we have seen so far and the condition on entropy that one can obtain from thermodynamics, happen to coincide. We shall see in §3.4 that this is not always the case. The geometric condition is a stronger condition and may be needed for the construction of a unique solution even under circumstances where thermodynamics has little to say. We already saw such an example in the equation ut + uux = 0. We close this section with two remarks.

3.2 Shocks

133

p

increasing pρ–γ

pρ–γ = constant

compressive shocks

(τ0,p0) τ

Figure 3.2.10. Entropy is an increasing function of pρ−γ and increases across a shock.

Remark 1 As long as the solution of our system of equations is smooth, the initial value problem is reversible. For example consider the conservation law ut + (u2 /2)x = 0. If the solution is smooth for all t satisfying 0 ≤ t ≤ T , and if we know the solution u(x, t), then by the change of variables u = −u, t = −t we may solve backward to recover the initial data u(x, 0). As soon as our solution becomes discontinuous, then the reversibility of the solution is lost. Consider a solution with prescribed initial data u(x, 0), where the characteristics first collide at time T1 ; see Figure 3.2.11. At a later t T1 x

A

B

Figure 3.2.11. If characteristics collide at time T1 , the system loses its reversibility.

134

3 Gas Flow in One Dimension

time T > T1 , it will be impossible to uniquely reconstruct the initial data prescribed between A and B, because the characteristics containing these initial data have been “swallowed up” by the shock. Thus, the information between A and B is “lost,” which confirms the notion that an increase in entropy means that information has been lost.

Remark 2 In the case of nonconstant pρ−γ , we say that a state of a gas is a set of three values S = (ρ, u, p). We ask the question: What constant states can be connected to a given state by a shock? (see Figure 3.2.12.)

p

S = (ρr,ur,pr)

S = (ρl,ul,pl) τ
Figure 3.2.12. Can the left and right states be connected by a shock?

The Hugoniot relation yields a curve on which lie states that can be connected to Sr = (ρr , ur , pr ) by a shock; the curve represents the set of possible transitions. Furthermore, a shock that separates two constant states must be a straight line in the (x, t) plane, because both [f (u)] and [u] are constant in time and s[u] + [f (u)] = 0, thus s = constant. Define a centered wave to be either a straight line shock or an isentropic centered rarefaction. Then, given a right state Sr , and a pressure pl ≥ 0 in a left state Sl , we can find pl , ul such that Sl is connected to Sr by a centered wave. If pl > pr , we can find a straight line shock connecting the two states, and if pl ≤ pr , we can find a centered rarefaction wave that connects them.

3.3 The Riemann Problem

135

3.3

The Riemann Problem

The conservation laws for a γ-law gas were shown in the previous section to be ρt + (ρ u)x = 0, (ρ u)t + (ρ u2 + p)x = 0, et + ((e + p)u)x = 0, where e = 1 ρ u2 + p/(γ − 1). The state of the gas is a vector 2   ρ u =  u . e The Riemann problem is the initial value problem for (3.3.1) with special initial data of the form u(x, 0) = where  ρr ur =  ur  er  ur , x ≥ 0, ul , x ≤ 0,  ρl and ul =  ul  el  (3.3.1)

are two constant states of the gas. The special case in which the velocity components vanish, that is, ur = 0 and ul = 0, is called a shock tube problem. The main objective of this section is to show how to solve the Riemann problem and how to use the solution to construct the solution of (3.3.1) with general initial data. Consider a change of variables x = Lx, t = Lt,

where L > 0. Clearly this leaves the form of equations (3.3.1) unchanged and the initial data are unchanged as well. Thus, if we assume that the solution is unique, then u(x, t) = u(x , t ) = u (xL, tL) = u x , t t > 0.

Thus, the solution of the Riemann problem is constant along straight lines issuing from the origin in the (x, t) plane. Because hyperbolic equations have a finite speed of propagation of data, we conclude that at any instant,

136

3 Gas Flow in One Dimension

the state of the gas is ur far enough to the right and is ul far enough to the left. Let the corresponding regions in the (x, t) plane be denoted Sr and Sl , respectively. From the results of §3.1, the boundary of Sr is either a C+ characteristic of (3.3.1) or a jump discontinuity emanating from the origin. A similar statement holds for Sl . We will now give a plausible reason why Sr can be connected to Sl through a right-centered wave, a constant state I, a slip line, a constant state II , and a left-centered wave (see Figure 3.3.1). Recall from the end of the previous section that a centered wave is either a centered isentropic rarefaction wave or a shock. We now sharpen our definitions slightly. A right-centered wave is a wave that is either a shock wave facing to the right or an isentropic centered rarefaction wave in which Γ− is constant. left-centered wave constant state II t slip line dx = u* dt right-centered wave Sr ur x

Figure 3.3.1. The constant states Sr and Sl connected through constant states, a slip line, and a left-centered wave.

Similarly, a left-centered wave is a wave that is either a shock wave facing to the left or a Γ+ isentropic rarefaction. From the discussion at the end of the previous section, given a right state Sr , we can find a oneparameter family of constant states I (parametrized by the pressure, for example), which can be connected to Sr by a right-centered wave. Note that because Sr is a constant state, pρ−γ is a constant in Sr , and if the right wave is a rarefaction, pρ−γ is also constant in that rarefaction; thus, it is consistent to allow only isentropic rarefaction waves in the definition of a right-centered wave. Because the density is the only quantity that can be discontinuous across a slip line, we have a two-parameter family of constant states II that can be connected through a slip line, a constant state I, and a right-centered wave to Sr . Continuing once more, we get a three-parameter family of states that can be connected to Sr by a left-centered wave, a constant state II , a slip line, a constant state I, and a right-centered wave. The question is whether or not we can choose the parameters so that we end up with the desired constant state ul . If we can, we have a solution of the Riemann problem. If ul and ur are close enough, one can demonstrate the result by means of the implicit function theorem. To demonstrate the

{
Sl ul

constant state I 0

3.3 The Riemann Problem

137

result for general ul and ur (for a γ-law gas), let p∗ and u∗ denote the pressure and velocity of the constant states I and II , respectively. Define the quantities pr − p∗ pl − p∗ Mr = − , Ml = − ur − u∗ ul − u∗ (see formula (3.2.13)). We claim that these quantities are functions of p∗ only, that is, Mr = Mr (p∗ ) and Ml = Ml (p∗ ). To prove this claim, we assume first that the right wave is a shock. Then p∗ > pr . From (3.2.14),
2 Mr = −

p r − p∗ τr − τl

(3.3.2)

where τl is the specific volume of the constant state I. The Hugoniot equation (3.2.16) gives (τl − µ2 τr )p∗ = (τr − µ2 τl )pr where µ2 = (γ − 1)/(γ + 1). Therefore, τl = τ r Thus, using (3.3.2),
2 M r = ρr

pr + µ2 p∗ . µ2 pr + p∗

µ2 1 p + p∗ . 2 r 1−µ 1 − µ2

This shows that Mr is a function of p∗ only, because the state Sr is known from the initial data. If the right wave is not a shock, it is a centered rarefaction wave. Of course, the same argument applies to Ml . On the other hand, from the definition of Mr and Ml , u∗ = u∗ (Mr ) = ur + pr − p∗ , Mr pl − p∗ u∗ = u∗ (Ml ) = ul + , Ml

Because u∗ is continuous across the slip line, we must have u∗ (Mr ) = u∗ (Ml ). This is an equation for p∗ . Some elementary algebra shows that this equation has a unique real solution p∗ (for a γ-law gas). Hence, the Riemann problem can be solved. In particular, in the case of the shock tube problem, one of the waves must be a shock and the other a rarefaction wave. The position of the slip line determines which one is the shock. If the slip line moves with positive velocity, that is, in the positive x-direction, then it acts

138

3 Gas Flow in One Dimension

as if it were a piston pushing into the gas in the positive x-axis. Therefore, the right wave must be a shock, and the left wave is a rarefaction.5 Suppose now that we have arbitrary initial data u(x, 0) = f (x) for the equation (3.3.1) where we assume f to be of compact support for convenience. We will use the preceding solution of the Riemann problem to construct a family of approximate solutions that converge to a solution of our initial value problem for (3.3.1). The x-axis and the t-axis are divided into intervals of length h and k, respectively. The approximate solution is to be computed at the mesh points (ih, jk) and i + 1 h, j + 1 k . We 2 2 let j+1/2 uj and ui+1/2 i approximate u(ih, jk) and u respectively. Thus, initially, u0 = f (ih), i (See Figure 3.3.2.) t 1 u–1

i+

1 2

h, j +

1 2

k ,

|i| = 0, 1, 2, . . . .

u1 0

u1 1 t=k

u1/2 –1/2
– 3h 2 –h 2

Pi u1/2 1/2 h 2 3h 2

t = k/2

0 u–2 –2h

0 u–1 –h

u0 0 0

u0 1 h

u0 2 2h x

Figure 3.3.2. Construction of the solution of (3.3.1) for general initial data.

The method described next permits us to advance by a time step k. Suppose we have already constructed the values of uj at time t = jk for i j+1/2 each i. To compute ui+1/2 , we consider the Riemann problem for (3.3.1) with the initial data u(x, 0) = x< i+ uj , i uj , x ≥ i + i+1
1 2 1 2

h; h.

5 This solution is given in S. K. Godunov, Mat. Sbornik 47 [1957], 537. See also P. D. Lax, Comm. Pure Appl. Math. 10 [1957], 537.

3.3 The Riemann Problem

139

Let v(x, t) be the solution to this Riemann problem. We will define ui+1/2 j+1/2 to be a randomly chosen value of v(x, k/2) on the interval ih ≤ x ≤ (i+1)h. For this purpose, define Λj = max i uj + cj i i

where uj is the velocity component of uj and cj is the corresponding sound i i i speed obtained from the formula c = γp/ρ . We require that the time step k satisfy the Courant–Friedrichs–Lewy condition Λj k ≤ h for all j. (3.3.3) We pick a random variable θ equidistributed in [−1/2, 1/2], that is, θ has the probability density function that takes the value 1 in [−1/2, 1/2], and is zero otherwise. Then we define ui+1/2 = v(Pi ), where Pi = ((i + 1/2) h + θh, j + k/2). Therefore, we have obtained the approximate values ui+1/2 j+1/2 j+1/2

at time t = j +

1 2

k.

With exactly the same construction (which involves another independent random variable θ equidistributed in [−1/2, 1/2]), we can advance a further half-time step to get the approximate values uj+1 i at time t = (j + 1)k.

Thus, we have constructed an approximate solution. We note that condition (3.3.3) is necessary to ensure that, at each half-time step construction, the centered waves generated by each Riemann problem will not interact. When the total variation (as defined in §3.2) of the initial data f (x) is sufficiently small, it can be shown that the time step k can always be chosen to satisfy (3.3.3), and that the approximate solution tends to a weak solution of our initial value problem as h tends to zero.6 Computational experience seems to support the belief that the approximate solutions converge to the solution for all data.7
6 See J. Glimm, “Solution in the large of hyperbolic conservation laws,” Comm. Pure Appl. Math. 18 [1965], 69. 7 A. J. Chorin, “Random choice solution of hyperbolic systems,” J. Comp. Phys. 22 [1976], 519.

140

3 Gas Flow in One Dimension

We will now give some simple examples to illustrate this construction. Consider the single linear equation vt = vx , with v(0) = f (x).

We already know that the solution is v(x, t) = f (x + t). The initial data propagate along the characteristics with speed 1. Thus, Λj = 1 in each step of our construction, and the Courant–Friedrichs–Lewy condition (3.3.3) becomes k ≤ h. We will carry out the construction under this condition. To compute vi+1/2 j+1/2 j j from vi and vi+1 ,

we consider the Riemann problem and j v(x, 0) = vi j+1 v(x, 0) = vi

for x < i + for x ≥ i +

1 2 1 2

h, h,

The initial discontinuity propagates along the characteristic passing through the point (i + 1 ) h, 0 , as in Figure 3.3.3. 2 characteristic with slope –1 Pi t

vji + 1/2 + 1/2

t = k/2

vji ih x = ih + (1/2)h

vji +1 (i+1)h

x

Figure 3.3.3. The Riemann problem for vt = vx .

Hence, if the randomly chosen point Pi lies to the right of the characteristic, then j+1/2 j vi+1/2 = vi+1 , that is, the solution moves h/2 to the left. Because this characteristic intersects the line t = k/2 at the point i+ 1 2 h− k k , 2 2

3.3 The Riemann Problem

141

and θ is uniformly distributed, Pi lies to the right of the characteristic with probability (h + k)/(2h) and lies to the left with probability (h − k)/(2h). n Therefore, after 2n steps, the approximate value vi at x = ih and t = nk comes from the initial disturbance f at the position x − Xn , where
2n

Xn = j=1 ηj ,

and ηj are independent, identically distributed random variables with the probability distribution Prob ηj = − and Prob ηj = h 1 = (h + k), 2 2h

h 1 = (h − k). 2 2h

The expectation and the variance of η = ηj are E[η] = and Var[η] = E[η 2 ] − (E[η])2 = Hence,
2n 1 4

1 h 1 h k (h + k) − + (h − k) = − , 2h 2 2h 2 2 h2 − k 2 .

E[Xn ] = j=1 E[ηj ] = −nk = −t,

and Var[Xn ] =

2n

Var[ηj ] = j=1 k h2 n 2 − 1 t. h − k2 = 2 2 k2

If we keep the ratio h/k constant, Var[Xn ] tends to zero for fixed t as h tends to zero. Thus, we have shown that the solution constructed by this random choice method tends to the correct solution as h tends to zero. As another example, we apply this construction directly to a Riemann problem. Consider a Riemann problem for (3.3.1) with initial data u(x, 0) = ur ul for x ≥ 0, for x < 0.

We assume that corresponding regions Sr and Sl can be connected by a shock propagating with speed s. When we divide the x-axis into pieces of

142

3 Gas Flow in One Dimension

t shock State Sl State Sr t = k/2

x –2h –h 0 h 2h

Figure 3.3.4. Breaking a Riemann problem into smaller ones.

length h, we have only one nontrivial Riemann problem at each half-time step construction (see Figure 3.3.4). This nontrivial Riemann problem is also solved by connecting Sr to Sl with the appropriate connecting shock. Hence, the argument given for a single linear equation shows that the discontinuity of the approximate solution is at a position whose expected value is correct, and that the variance of the position also tends to zero as h tends to zero, that is, the approximate solution tends to the correct solution. A similar conclusion applies to the case where Sr can be connected to Sl either by a centered rarefaction wave or by a slip line. Thus, we obtain the conclusion that, if we apply this random choice method directly to the simplest Riemann problems, the approximate solutions tend to the correct solution as h tends to zero. Note that even if the flow is not isentropic, that is, pρ−γ = constant, the construction uses Riemann problems in which pρ−γ = constant on either side of the slip line in the absence of a shock. The construction just presented is the basis of Glimm’s existence proof for hyperbolic systems.8 For the 2×2 isentropic flow equations, other methods of proof are available.9 A major component in the Glimm construction is the solution of the Riemann problem. We have offered a construction of a solution of the Riemann problem but have said nothing about its uniqueness. The existence and uniqueness of the solution of the Riemann problem for gas dynamics subject to an appropriate formulation of the entropy condition have been
Glimm, Loc. cit. for example, R. DiPerna, “Convergence of the Viscosity Method for Isentropic Gas Dynamics,” Comm. Math. Phys., 91 [1983], 1.
9 See, 8 J.

3.4 Combustion Waves

143

established.10 It is worth noting that not very much is known rigorously about the solutions of the equation of gas dynamics in more than one space variable. (See Majda’s book listed in the Preface for some information.) Glimm’s construction is also the basis for numerical algorithms.11 The strategy for picking the numbers θ is the determining factor in the accuracy of these algorithms.

3.4

Combustion Waves

This section examines some additional features of conservation laws and applies them to a modified system of gas dynamic equations that allows for combustion. We shall begin by studying a single conservation law ut + (f (u))x = 0 (3.4.1)

and concentrate on the information the shape of the graph of f gives about discontinuities. Later, we generalize the results to systems and apply them to our specific system. We examine in four steps the cases where the graph of f is straight, concave up, concave down, and, finally, neither concave up nor down. Case 1 f is linear ; that is, f (u) = au, a = constant. Here (3.4.1) becomes ut + aux = 0

(3.4.2)

whose characteristics are the straight lines dx/dt = a. Because the characteristics do not intersect, smooth initial data propagate to a smooth solution; in fact, the solution with u(x, 0) = ϕ(x) is u(x, t) = ϕ(x − at). However, discontinuities in the initial data are propagated along characteristics. This result is also a consequence of our formula for the speed of the discontinuity (see equation (3.2.8)) s= f (ul ) − f (ur ) = a, ul − ur

where ul is the value of u to the left of the discontinuity and ur that to the right. In particular, a discontinuity must be a characteristic. See Figure 3.4.1.
10 See T. P. Liu, “The Riemann Problem for General Systems of Conservation Laws,” J. Differential Equations, 18 [1975], 218. 11 See A. J. Chorin, “Random Choice Solution of Hyperbolic Systems,” J. Comp. Phys., 22 [1976], 517.

144

3 Gas Flow in One Dimension

t

f

f (u) = au (ur, f (ur)) (ul, f (ul))

ul

Σ: speed = a ur x x

Figure 3.4.1. If f is linear, discontinuities move along characteristics.

Case 2 f is concave up; that is, fuu > 0. (An example is ut + uux = 0, whose characteristics are again straight lines, as we showed in §3.2, and where f (u) = u2 /2.) The characteristics of (3.4.2) are given by dx = f (u) dt because (3.4.1) is equivalent to ut + f (u)ux = 0 in the region where u is smooth. We still have s= f (ul ) − f (ur ) ul − ur

which, by the mean value theorem, gives s = f (ξ) for some ξ between ul and ur . See Figure 3.4.2. Let us examine the two possibilities ur < ul and ul < ur separately. Because fuu > 0, f is increasing. Thus, if ul < ur , f (ul ) < s < f (ur ) and 1 1 1 > > . f (ul ) s f (ur ) Thus, the slopes of the characteristics have the configuration shown in Figure 3.4.3. We get a shock configuration that violates the entropy condition. Thus, if ul < ur these states should not be connected by a shock; a rarefaction fan must be used. On the other hand, if ur < ul , then f (ur ) < s < f (ul ) and 1 1 >s> . f (ur ) f (ul ) This time the characteristics enter the discontinuity as in Figure 3.4.4, which is consistent with the entropy condition.

3.4 Combustion Waves

145

f

slope of chord = speed of discontinuity

ur

ul

x

Figure 3.4.2. The discontinuity speed in the case f is concave up.

t Σ

slope = 1 f '(ul) slope = 1 f '(ur) x

Figure 3.4.3. The slopes of the characteristics for f concave up and ul < ur .

Case 3 f is concave down, that is, fuu < 0. An argument similar to that of case 2 shows that for a shock we must have ul < ur and for a rarefaction fan, ur < ul . In either case 2 or case 3, note that the discontinuity separates the characteristics. If either fuu < 0 or fuu > 0, f is called convex. Case 4 This is the general case in which f is neither linear nor concave up nor concave down. Given ul and ur , define the linear function l(u) by l(u) = f (ul ) − f (ur ) (u − ul ) + f (ul ), ul − ur

that is, the graph of l is the straight line through the pair of points (ur , f (ur )) and (ul , f (ul )). If the graph of f is concave up on an inter-

146

3 Gas Flow in One Dimension

f (u) slope = 1 f '(ul)

slope =

1 f '(ur) x

Figure 3.4.4. The slopes of the characteristics for f concave up and ur < ul .

val [a, b] containing both ur and ul , then by elementary calculus, l(u) > f (u) whereas l(u) < f (u) if f is concave down. For general f , one has the following theorem:12 Solutions exist, are unique, and depend continuously on the initial data (in a certain function space) if one allows only discontinuities that satisfy: (a) if ur > ul , then l(u) ≤ f (u) for all u ∈ [ul , ur ]; and (b) if ur < ul , then l(u) ≥ f (u) for all u ∈ [ur , ul ]. This is known as Oleinik’s condition (E). In summary, if f is convex (concave up or concave down), elementary calculus shows that a discontinuity allowed by the condition separates the characteristics; that is, two characteristics cross the graph of the discontinuity at each point of the (x, t) plane, and either they both point forward in time or they can both be traced backward in time. The entropy condition rules out the former possibility. One can readily verify that Oleinik’s condition (E) rules out the same shocks as the entropy condition; for example, if fuu > 0, condition (E) requires that ur < ul across a shock. If f is not convex, condition (E) is a generalization of the entropy condition; indeed, if f is not convex, it is not obvious exactly what our earlier entropy condition does or does not allow. Condition (E) forbids in particular shocks that move faster or slower than the characteristics on both of their sides.
12 O. A. Oleinik, “Existence of Solutions for a Simple Hyperbolic Equation,” Amer. Math. Soc. Transl. Ser. 2, 285 [1965], 33.

3.4 Combustion Waves

147

Notice that a discontinuity can satisfy part (a) of Oleinik’s theorem without f being concave down on the whole interval [ul , ur ], as in Figure 3.4.5. f (u)

(ur, f (ur))

(ul, f (ul)) u

Figure 3.4.5. A situation obeying (a) of Oleinik’s condition.

Consider next a situation in which the line joining the point (ul , f (ul )) to (ur , f (ur )) is neither wholly above nor wholly below the graph of f , so Oleinik’s condition (E) is violated; see Figure 3.4.6. f (u)

(ur, f (ur)) (ul, f (ul))

u

Figure 3.4.6. A situation violating Oleinik’s condition.

Thus, we cannot connect ur and ul by means of a single shock. One can also show13 that they cannot be connected by a single rarefaction fan. One can prove, however, that the gap between ur and ul can be bridged by a
13 See

Oleinik, Loc. cit.

148

3 Gas Flow in One Dimension

compound wave consisting of shocks, rarefaction fans, and slip lines, all consistent with the entropy condition. We now consider systems of conservation laws. Recall from §3.1 that we have a family of characteristics associated with each eigenvalue of the coefficient matrix A(u) when our equations are written in hyperbolic form (3.1.3). For gas dynamics we have defined the families C+ , C− , C0 of characteristics. We say that a family of characteristics crosses a discontinuity that satisfies the jump conditions if through every point of the graph of the discontinuity in the (x, t) plane, one can draw only one characteristic of that family, with that characteristic being traceable backward in time on one side and forward in time on the other. Intuitively, a family of characteristics crosses a discontinuity that is not needed to prevent characteristics of that family from intersecting. A family of characteristics is called linearly degenerate if a discontinuity that satisfies the jump conditions can coincide with a member of that family. For example, the C0 family is linearly degenerate, because a slip line is also a C0 characteristic. We said earlier that a discontinuity separated a family of characteristics if through each point of the graph of the discontinuity in the (x, t) plane there exist two characteristics, with the two of them either traceable back in time or forward in time. (Slip lines do not separate the C0 characteristics, nor are they crossed by them.) A family of characteristics is called convex if, whenever a discontinuity is allowed by the algebraic jump conditions, it is either crossed by the family or it separates the family. If a family of characteristics is convex, the entropy condition rules out discontinuities crossed by characteristics that can be traced forward in time. The Prandtl relation u0 u1 = c2 ∗ derived in §3.2 shows that the C+ and C− characteristics of gas dynamics for a γ-law gas are convex families. Thus, these equations admit only two kinds of discontinuities: discontinuities that separate either the C+ or the C− family of characteristics, and linearly degenerate slip lines. A system of conservation laws is said to be convex if all the discontinuities allowed by the jump conditions separate one (and only one) of the families of characteristics and if all the families are convex. A single equation can be readily seen to be convex in this sense if and only if f (u) is convex. The system of isentropic gas dynamic equations (i.e., pρ−γ = constant) is thus convex, because it has only the C+ and C− families of characteristics, both of which are convex. A system of conservation laws is said to be linearly degenerate if at least one of its families of characteristics is linearly degenerate and the others are convex. Thus, the 3×3 system of gas dynamics is linearly degenerate because the C0 family is linearly degenerate. In other cases, the system is called nonconvex , and compound waves are needed to connect states and construct a solution of the Riemann problem

3.4 Combustion Waves

149

(and thus of the initial value problem). Examples of nonconvex systems are provided by the equations of gas dynamics with combustion, to which we turn next.14 The equations of gas dynamics with combustion are ρt + (ρ u)x = 0, (ρ u)t + (ρ u2 + p)x = 0, et + ((e + p)u)x = 0, where e = 1 ρ u2 + ρ 2 and = pτ + q. γ−1 (3.4.3)

The only modification from the gas dynamic equations for a γ-law gas that we have already studied is the presence of the term q in the energy density; we have yet to specify q. The term q represents chemical energy that can change due to burning. We assume here for simplicity that q takes only two values: q0 , if the gas is unburned, q= (3.4.4) q1 , if the gas is burned, where q0 and q1 are constants. To complete our model, we assume that the burning reaction occurs when some threshold is reached as follows: Let T = p/ρ measure the temperature. We assume that the gas burns when T ≥ Tc for a constant Tc called the ignition temperature. Let ∆ = q1 − q0 ; the reaction is called exothermic if ∆ < 0, that is, q1 < q0 , and is endothermic if ∆ > 0, that is, q1 > q0 . The third equation in (3.4.3) asserts that the “total” energy is conserved, where the “total” energy is the sum of the kinetic energy, the internal energy, and a chemical energy, part of which can be released through a process we shall call “burning” or “combustion.” In this model we have neglected viscous effects, heat conduction, and radiative heat transfer. Of these the most serious is the exclusion of the effects of heat conduction.15 We also assume that the “combustion” only occurs once; after the gas is burned, it stays burned. We have already assumed in writing (3.4.3) that we are dealing with a γ-law gas where γ is the same for both burned and unburned gas. This is not normally realistic, but the mathematics in the more general situation is essentially the same as in our simplified situation. We shall furthermore assume that the region in which combustion occurs (and in which q changes from q0 to q1 ) is infinitely thin, and that the transition occurs instantaneously.
14 Another interesting nonconvex physical system is analyzed in A. J. Chorin, Comm. Math. Phys. 91 [1983], 103. 15 For the equations for flow with finite conduction and reaction rates, see F. A. Williams [1965] Combustion Theory, Addison-Wesley.

150

3 Gas Flow in One Dimension

We can picture the gas as spread along the x-axis; at any instant some sections of the gas are burned (q = q1 ) and others are not (q = q0 ). The lines between them will be associated with discontinuities Σ in the flow. A discontinuity across which ∆ = 0 is called a reaction front or a combustion front. The jump relations across a discontinuity Σ are given by the relations described previously. In a frame moving with the discontinuity, they read: [ρu] = 0, ρu + p = 0, [(e + p)u] = 0. Letting M = ρ u, we derive M2 = p1 − p1 τ0 − τ1 (3.4.5)
2

as before (see equation (3.2.14)). We can eliminate u1 and u2 from the equation to obtain a Hugoniot relation by the methods of §3.2. We still define the Hugoniot function with center (τ0 , p0 ) to be H(τ, p) = −
0

+

τ0 − τ (p0 + p). 2

(3.4.6)

so the energy jump condition becomes H(τ1 , pl ) = 0. We can rewrite (3.4.6) as 2µ2 H(τ, p) = p(τ − µ2 τ0 ) − p0 (τ0 − µ2 τ ) + 2µ2 ∆, where µ2 = γ−1 . γ+1

This differs from our expression in §3.2 by the term 2µ2 ∆. In particular, note that a state cannot be connected to itself by a reaction front because H(τ0 , p0 ) = ∆. There can be discontinuities that are not reaction fronts (∆ = 0), in which case our analysis in §3.2 applies. For a given ∆ = 0, we still call the curve H(τ, p) = 0 the Hugoniot curve. It is the locus of all possible states that can be connected to the given state (τ0 , p0 ). We shall now discuss reaction fronts for exothermic reactions; that is, ∆ < 0. In this case the state (τ0 , p0 ) lies below the hyperbola representing the Hugoniot curve, as Figure 3.4.7 shows. One portion of the Hugoniot curve (from A to B) is omitted because the states in it correspond to negative M 2 (by (3.4.5)), which is an impossibility. The remaining part of the curve is divided into two branches. The

3.4 Combustion Waves

151

p

the Hugoniot curve I = strong detonation CJ1 = Chapman-Jouguet detonation II = weak detonation A CJ2 = Chapman-Jouguet deflagration B III = weak deflagration

p0

Rayleigh line (τ0,p0)

Rayleigh line τ0

IV = strong deflagration

τ

Figure 3.4.7. For exothermic reactions, (τ0 , ρ0 ) lies below the Hugoniot curve.

upper branch, which corresponds to a compressive reaction front (p > p0 ), is called the detonation branch. The lower branch, which corresponds to an expansive reaction front (p < p0 ), is called the deflagration branch. The lines through (τ0 , p0 ) tangent to the Hugoniot curve are called the Rayleigh lines, and the points of tangency CJ1 and CJ2 are called the Chapman–Jouguet points. These points divide the Hugoniot curve further into four subbranches. The part I above CJ1 is called the strong detonation branch, the part II between CJ1 and A the weak detonation branch. Similarly, the part III between B and CJ2 is called the weak deflagration branch, the part IV below CJ2 the strong deflagration branch. We will discuss each branch separately, and thus exclude some branches on the basis of an analogue of Oleinik’s condition (E), that is, by a geometrical entropy condition. Consider a constant state S corresponding to a point (τ1 , p1 ) that lies on the Hugoniot curve, that is, it is possible to connect S and (τ0 , p0 ) by a chemical reaction front. Depending on the position of S, we consider the following cases: Case 1 S lies in the strong detonation branch, that is, p1 > pCJ1 (See Figure 3.4.8). We can connect S to (τ0 , p0 ) by a strong detonation front. One can show by the argument in §3.2 that the gas flow relative to the reaction front is supersonic in the front and subsonic in the back, that is, |u0 | > |u1 | < c1 .

152

3 Gas Flow in One Dimension

t state S

C0 C–

strong detonation C+ state (τ0,p0) C+ 0 C– x

C0

Figure 3.4.8. A state S on the strong detonation branch.

Therefore, if the strong detonation front faces in the positive x direction, it separates the C+ characteristics, as the figure shows. The other two families of characteristics cross the reaction front. Thus, the geometry of the characteristics in the strong detonation case is similar to what it is in the case of a shock. Hence, the pressure in the back state is sufficient to determine the back state from the front state. Case 2 S = CJ1 , that is, p1 = pCJ1 (See Figure 3.4.9). In this case the velocity relative to the reaction front is supersonic in the front and sonic in the back, that is, |u0 | > c0 and |u1 | = c1 .

Thus, the Chapman–Jouguet reaction front, when observed from the burned gas state, is one of the characteristics (of C+ or C− family). The condition |u1 | = c1 enables us to determine the back state from the front state without further assumption. Case 3 S lies in the weak detonation branch, that is, pCJ1 > p1 ≥ pA (See Figure 3.4.10). In this case one can show that the gas flow relative to the reaction front is faster than the sound speed on either side.

3.4 Combustion Waves

153

t C0 C+ C+ C– state S Chapman–Jouguet detonation state (τ0,p0) C+ 0 C0 C– x C+

Figure 3.4.9. A state S corresponding to the Chapman Jouguet point CJ1 .

t C–

C0

C+

state S weak detonation

state (τ0,p0)

C+ 0

C0

C–

x

Figure 3.4.10. A state S on the weak detonation branch.

Each family of characteristics crosses the reaction front, as Figure 3.4.10 shows. This violates the analogue of Oleinik’s condition (E). Instead of connecting S and (τ0 , p0 ) by a single wave front, we use a compound wave to connect them. We connect (τ0 , p0 ) first to the Chapman–Jouguet detonation and followed by an isentropic-centered rarefaction wave to reach the

154

3 Gas Flow in One Dimension

state S. This compound wave is possible because the Chapman–Jouguet detonation moves with sound speed with respect to the gas flow in its back side. If we adopt this compound wave to connect them, the pressure is S and the state (τ0 , p0 , u0 ) will determine S completely. If weak detonations are not excluded, the solution of the initial value problem is not unique; indeed, the state S can then be connected to a state with a pressure p0 by either a weak detonation or a CJ detonation followed by a rarefaction. The exclusion of the weak detonation is analogous to the exclusion of discontinuities that move faster than the characteristics on both their sides by means of condition (E). Case 4 S lies in the weak deflagration branch, that is, p0 ≥ p1 > pCJ2 (See Figure 3.4.11). t C–

C0 weak deflagration

state S C+ state (τ0,p0) C+

C0 0

C–

x

Figure 3.4.11. A state S on the weak deflagration branch.

One can show that the reaction front moves with respect to the gas slower than the sound speed on both sides. This is an indeterminate case. One can determine the solution uniquely only by taking into account heat conduction and a finite reaction rate. In the limiting case of no heat conduction

3.4 Combustion Waves

155

(which we have assumed here), it can be shown16 that the only weak deflagration that can occur is one across which the pressure is constant and there is no mass flow, that is, this deflagration is indistinguishable from a slip line. Thus, only in the case where p1 = p0 can the state S be connected to (τ0 , p0 ) by a weak deflagration. Case 5 S lies in the strong deflagration branch, that is, p1 < pCJ2 (See Figure 3.4.12). t C–

C0

C+ state S strong deflagration C+ state (τ0,p0)

C0 0

C– x

Figure 3.4.12. A state S on the strong deflagration branch.

Then the gas flow relative to the reaction front is subsonic in the front side and supersonic in the back side. If this strong deflagration moves in the positive x direction, it separates the C+ characteristics. However, this separation does not satisfy the geometric entropy condition, although the strong deflagration is consistent with the conservation laws. Thus, we must exclude strong deflagrations. Note that we have excluded weak detonations and strong detonations by means of plausible geometrical entropy conditions; neither one could be
16 See A. J. Chorin, “Random choice methods with applications to reacting gas flow,” J. Comp. Phys., 25 [1977], 253.

156

3 Gas Flow in One Dimension

excluded by the requirement that the physical entropy should increase as the gas crosses the discontinuity. They could, of course, be excluded by other considerations of a plausible physical nature. (The theory of such waves is still in a state of development.)17 These considerations permit us to propose a solution for the Riemann problem for the system (3.4.3). We define a centered wave to be either a shock, a centered rarefaction, a strong detonation, or the compound wave consisting of the Chapman–Jouguet detonation followed by a centered rarefaction. A Riemann problem for the system (3.4.3) is an initial value problem for (3.4.3) with initial data of the form   ρ  u     e  = Sr , if x ≥ 0; q and   ρ  u     e  = Sl , q

if x < 0,

where Sr and Sl are two constant states of the gas. We claim that the Riemann problem is solvable by connecting Sr to Sl through a right-centered wave, a constant state I , a slip line, a constant state II , and a left-centered wave (Figure 3.4.13) The argument is similar to the case of an inert gas (§3.3). If we can determine the constant state I from Sr and a given pressure p∗ in I, we can proceed as in §3.3 to get an algebraic equation for p∗ through the continuity of the velocity across the slip line, and hence obtain a solution to the Riemann problem. For this purpose, we assume first that Sr is in an unburned state. If constant state I is in an unburned state, we are back to the case of inert gas. If constant state I is in a burned state, we can connect I to Sr by a centered wave depending on the position of p∗ in the Hugoniot curve with center Sr . The critical criterion is whether the temperature T computed in the constant state I exceeds the ignition temperature Tc or not. One can actually show that there is a consistent way of solving the Riemann problem.18 If Sr contains burned gas, the constant state I also contains burned gas, because the burning can occur only once.
17 A. Bourlioux, A. Majda, and V. Roytburd, “Theoretical and numerical structures for unstable one dimensional detonations,” SIAM J. Appl. Math., 51,[1991], 303–343. P. Colella, A. Majda, and V. Roytburd, “Theoretical and numerical structure of reactive shock waves,” SIAM J. Sci. Comput., 1, [1980], 1059–1080. 18 See A. J. Chorin, Random choice methods with applications to reacting gas flow, J. Comp. Phys., 25 [1977], 253 for the details.

Vector Identities

157

slip line left-centered wave constant state II Sl

t

right-centered wave constant state I Sr

0

x

Figure 3.4.13. The Riemann problem for the gas flow with combustion.

Having solved the Riemann problem, we can use the random choice method described in §3.2 and §3.3 to obtain the approximate solutions for the general initial value problem for the system (3.4.3). Results similar to those in §3.3 can be expected. We need not repeat the construction.

Vector Identities
The following gives some general formulas that are useful for calculations with vector fields in R3 .

158

Vector Identities

1. ∇(f + g) = ∇f + ∇g 2. ∇(cf ) = c∇f , for a constant c 3. ∇(f g) = f ∇g + g∇f 4. ∇ f g = g∇f − f ∇ g2

5. div(F + G) = div F + div G 6. curl(F + G) = curl F + curl G 7. ∇(F · G) = (F · ∇)G + (G · ∇)F + F × curl G + G × curl F 8. div(f F) = f div F + F · ∇f 9. div(F × G) = G · curl F − F · curl G 10. div curl F = 0 11. curl(f F) = f curl F + ∇f × F 12. curl(F × G) = F div G − G div F + (G · ∇)F − (F · ∇)G 13. curl curl F = grad div F − ∇2 F 14. curl ∇f = 0 15. ∇(F · F) = 2(F · ∇)F + 2F × (curl F) 16. ∇2 (f g) = f ∇2 g + g∇2 f + 2(∇f · ∇g) 17. div(∇f × ∇g) = 0 18. ∇ · (f ∇g − g∇f ) = f ∇2 g − g∇2 f 19. H · (F × G) = G · (H × F) = F · (G × H) 20. H · ((F × ∇) × G) = ((H · ∇)G) · F − (H · F)(∇ · G) 21. F × (G × H) = (F · H)G − H(F · G) Notes. In identity 7, V = (F·∇)G has components Vi = F·(∇Gi ), for i = 1, 2, 3, where G = (G1 , G2 , G3 ). In identity 13, the vector field ∇2 F has components ∇2 Fi , where F = (F1 , F2 , F3 ). In identity 20, (F×∇)×G means ∇ is to operate only on G in the following way: To calculate (F × ∇) × G, we define (F × ∇) × G = U × G

Vector Identities

159

where we define U = F × ∇ by: i U=F×∇= F1 ∂ ∂x j F2 ∂ ∂y k F3 ∂ ∂z .

160

Vector Identities

This is page 161 Printer: Opaque this

Index

acceleration of a fluid particle, 4 algorithm, 94 almost potential flow, 59 asymptotically stable, 97 autonomous, 97 back, 124 balance of momentum, 4, 6, 12 differential form, 6 integral form, 7 Bernoulli’s theorem, 16, 48, 55, 71 bifurcation, 99, 100 Blasius’ theorem, 52 body force, 6 force on, 52 boundary condition, 34, 41 layer, 67, 68 approximation, 76 equation, 75 separation, 79 thickness, 71 vorticity in, 76 layer separation, 71

burning, 152 Cauchy’s theorem, 32 centered rarefaction wave, 114 wave, 137, 138, 158 central limit theorem, 84 channel flow, 17 Chapman–Jouguet points, 154 characteristic, 103, 105, 108 intersecting, 106, 108 length, 35 linearly degenerate, 150 velocity, 35 chemical energy, 151 circulation, 21, 48, 57 coefficients of viscosity, 33 combustion, 103, 152 front, 152 wave, 145 complex potential, 51 variables methods, 51 velocity, 51, 53 compressible

162

Index

flow, 14, 33, 103 compressive shock, 131, 134 concave up, 146 conformal transformation, 65 connected state, 115 wave, 115 conservation law, 122, 129, 145 of energy, 12 of mass, 2, 12 of vorticity, 28 consistency of an algorithm, 94 constant state, 111 contact discontinuity, 124 continuity equation, 3 continuum assumption, 2 convective term, 40 convergence of an algorithm, 94 convex characterisitics, 151 Couette flow, 31 Courant–Friedrichs–Lewy condition, 141 cylindrical coordinates, 45 d’Alembert’s paradox, 54 in 3d, 58 decompostion theorem, 37 deflagration branch, 154 deformation, 18 tensor, 19, 31 degenerate linearly, 151 density, 1, 14 derivative material, 5 detonation branch, 154 differential form, 120 form of mass conservation, 3 diffusion, 39 discontinuity, 121, 146, 148 separating, 148 disjoint, 82 event, 82

dissipation term, 39 divergence free part, 38 space-time, 119 downstream, 93 drag, 54, 57, 60, 66 form, 80 skin, 80 dynamics, 96 endothermic, 152 energy, 12 equation, 118 flux, 18 internal, 12, 117 kinetic, 12 per unit volume, 117 enthalpy, 14 entropy, 14, 118, 158 condition, 127, 130, 133, 157 equation continuity, 3 differential form, 120 Euler, 13, 15 heat, 84 Hugoniot, 125 Navier–Stokes, 34 of state, 44 Prandtl, 75 Stokes’, 40 vorticity, 24 weak form, 120 error function, 70 Euler equation, 13, 15, 49, 78, 94, 96 event, 82 disjoint, 82 exothermic, 152 expectation, 82, 143 field velocity, 1 vorticity, 18 filament, 65 Filon’s paradox, 67

Index

163

first coefficient of viscosity, 33 law of thermodynamics, 14, 118 fixed point, 97 flat plate, 69 flow, 14 almost potential, 59 around a disk, 55 around a half circle, 51 around an obstacle, 52 between plates, 41 between two plates, 42 compressible, 14, 33 Couette, 31 gas, 103 homogeneous, 11, 48 ideal, 48 in a channel, 17 in a pipe, 45 in the half-plane, 69 in the upper half-plane, 51 incompressible, 10 induced by a vortex filament, 65 irrotational, 47 isentropic, 14, 15 map, 7, 95 over a plate, 66 past a sphere, 59 Poiseuille, 45 potential, 47, 48 potential around a disk, 55 potential vortex, 56 stationary, 29, 49 supplementary region half-plane, 51 fluid flow map, 7 ideal, 5 particle, 4 velocity, 1 viscous, 33 flux, 7 of vorticity, 22

force, 5, 53 on a body, 52 form drag, 80 fourth power law, 46 front, 124 function error, 70 Green’s, 61, 86 gamma law gas, 114, 118, 125, 131, 139 simple wave, 111 gas dynamics, 103, 111 flow, 103 ideal, 118, 125, 131, 139 Gaussian, 84 generation of vorticity, 43 geometric entropy condition, 135 Glimm’s existence proof, 145 global stability, 97 gradient part, 38 Green’s function, 61, 63, 64, 86 half-plane flow, 69 Hamiltonian system, 62 heat equation, 84, 86 Helmholtz decomposition theorem, 37 theorem, 26, 37 Hodge theorem, 37 hodograph transformation, 110 homogeneous, 11 flow, 48 Hopf bifurcation theorem, 99 Hugoniot curve, 126, 153 equation, 125, 139 function, 125, 153 hyperbolic, 104 ideal flow, 48 fluid, 5, 31

164

Index

gas, 44, 118, 125, 131, 139 ignition temperature, 152 incompressible, 11 approximately, 44 flow, 10, 13, 34 independent, 82 random variables, 83 integral form, 120 form of balance of momentum, 7 form of mass conservation, 3 intensity of a vortex sheet, 88 internal energy, 12, 40, 117 invariant Riemann, 109 irrotational, 47 isentropic flow, 14, 15 fluids, 14 gas flow, 122 Jacobian, 8 matrix, 24 jump, 122 discontinuity, 121, 122 relations, 124 Kelvin’s circulation theorem, 21 kinematic viscosity, 34 kinetic energy, 12, 40, 49 Kutta–Joukowski theorem, 53 law conservation, 129, 145 of large numbers, 83 layer boundary, 67, 71 leading edge, 93 left -centered wave, 138 state, 115 length characteristic, 35 Liapunov stability theorem, 97

Lie derivative, 43 Lie–Trotter product formula, 95 line vortex, 22 linearly degenerate, 150, 151 Mach number, 44 mass conservation, 11 density, 1 flow rate, 46 matching solutions, 78 material derivative, 5, 18 mean, 82 mechanical jump relations, 124 momentum balance of, 6 flux, 7 moving with the fluid, 7 Navier–Stokes equation, 31, 33, 38, 67, 77, 94, 95 Neumann problem, 37, 49, 63 Newton’s second law, 2, 6 no-slip condition, 34, 43 noncompressive shock, 133 nonconvex, 151 nonlinear dynamics, 96 obstacle, 58 flow around, 52 Oleinik’s condition, 149, 154 orthogonal projection, 38 oscillations, 100 Oseen’s equation, 66 paradox d’Alembert, 54, 58 Filon, 67 Stokes, 67 pipe flow, 45 piston, 111 plate, 80 flow between, 41 flow over, 66, 69

Index

165

flow past, 87 point vortices, 60 Poiseuille flow, 45 potential complex, 51 flow, 47, 48, 51, 56, 58 almost, 59 flow around a disk, 55 velocity, 48 vortex, 56 vortex flow, 56 Prandtl boundary layer equation, 73, 75 equation, 78, 94 relation, 132, 151 pressure, 5, 14 probability, 82 density function, 83 theory, 82 product formula, 95 projection operator, 38 quasilinear, 104 random choice method, 144, 159 variable, 82, 141 Gaussian, 84 walk, 85, 88, 96 Rankine–Hugoniot relations, 124 rarefaction fan, 148 shock, 127, 129 wave, 114 Rayleigh lines, 154 reaction endothermic, 152 exothermic, 152 front, 152 reversibility, 136 Reynolds number, 36, 96 Riemann invariant, 109, 110, 113

problem, 103, 137, 139, 144, 158, 159 right -centered wave, 138 state, 115 rigid rotation, 18 rotation, 18 sample space, 82 scaling argument, 81 second coefficient of viscosity, 33 law, 118 separate characteristics, 130 separated, 150 separation boundary, 68, 71, 79 shadow of a vortex sheet, 90 sheet vortex, 22 shock, 117, 124 back, 124 compressive, 131, 133 front, 124 noncompressive, 133 rarefaction, 127, 129 separating, 130 tube problem, 137 similar flows, 36 simple wave, 111 simply connected, 47 skin drag, 80 slightly viscous flow, 47 slip line, 124, 138, 140 sound speed, 44, 104, 131 space-time divergence, 119 spatial velocity field, 1 speed discontinuity, 147 sphere flow past, 59, 66 stability, 96 stable point, 97 stagnation point, 29, 55 standard deviation, 83

166

Index

state, 111, 136, 137 connected, 115 equation of, 44 stationary, 49 flow, 16, 58 flow criterion, 29 steady flow, 58 Stokes equation, 40, 67, 96 flow, 66 paradox, 67 stream function, 29, 43, 54, 80 streamline, 16, 29, 44 strength of a vortex tube, 26 stress tensor, 32 stretched, 27 strong deflagration, 154 detonation, 154 subsonic, 131, 157 supersonic, 131 symmetry, 101 tangential boundary condition, 34 forces, 5 Tchebysheff’s inequality, 84 temperature, 14, 152 test functions, 119 theorem Bernoulli’s, 16 Blasius’, 52 Cauchy’s, 32 central limit, 84 Helmholtz, 26 Helmholtz–Hodge, 37 Kelvin circulation, 21 Kutta–Joukowski, 53 transport, 10 thermodynamics, 14, 118 first law, 14, 118 thickness boundary layer, 71, 73 total force, 8 trajectory, 16

transfer of momentum, 31 transformation hodograph, 110 translation, 18 transport theorem, 10, 18, 117 tube vortex, 26 upper half-plane flow in, 51 variance, 83, 143 variation, 129 velocity characteristic, 35 complex, 51 field, 1 potential, 48, 51 profile, 42 viscosity, 129 coefficients, 33 kinematic, 34 viscous fluid, 33 vortex, 61 filament, 65 line, 22 potential, 56 sheet, 22, 82, 87 sheet intensity, 88 tube, 26 tube, strength of, 26 vortices point, 60 vorticity, 18, 63 conservation, 28 creation operator, 96 equation, 28 generation of, 43 in boundary layer, 76 transport, 23 wave centered, 137 connected, 115 left-centered, 138

Index

167

rarefaction, 114 right-centered, 138 simple, 111 weak deflagration, 154 detonation, 154 solution, 118, 119

168

Index

This is page 169 Printer: Opaque this

For authors use only DO NOT PUBLISH THIS PAGE
A Mathematical Introduction to Fluid Mechanics
(The corrected fourth printing, April 2000) Prepared on April 3, 2000 & retypeset August 14, 2003.

This document was modified from the third corrected printing version typeset on 10-28-99.

Modifications:
Mar-24-2000: Move the Preface to page (v) before the Table of contents. Mar-24-2000: Replace the ZapfChancery Fonts with the mathcal fonts in figures 2.1.x.eps, x=3,4,5,6, and save as Illustrator 8 eps files. Apr-03-2000: Redo the frontmatter to include a Series Preface (Springer) page as page (v) and start the preface on (vii)…...

Similar Documents

Free Essay

Fluids

...   It   is generally  known  that  water  boils  at   100°C,  but  really  that  is  only  true  if  the  water  is  at  a pressure  of If  the  water  pressure  increases  so  too  does  its  boiling  point, and  vice  versa.  The  relationship  between  pressure  and  temperature  of boiling  water  can  be plotted  on  a  so-called  transition  or phase  diagram  (Figure  1).  In  this  lab,  the  boiling  of water will be studied  and  plotted on a transition  diagram. Liquid 9 s \ 1- Solid A J    Vapour / i! ! 0.01 Temperature Figure  1:   Phase  Diagram  of Water TM PDF Editor Lab Worksheet Apparatus The  apparatus comprises  a  boiler and  a  heating element  (Figure  2).  The  boiler (and  the fluid it  contains)  is  separated  from  the  surroundings  by  isolating valves. Pressure Saturated vapour \ valve T1 Filling point Isolating valve Viewing  port Boiler Heating  elements Figure  2 TH3  Saturation  Pressure  Apparatus In   setting  up   the   laboratory,   the   water   is   boiled   with  the   isolating  valves  open.   Vapour leaves  the  system,  pushing  out  any  atmospheric  gases  (oxygen  and  nitrogen).  During  your experimental   work,   you   will   heat   the   vessel   and   then   allow   it   to   cool,   monitoring   the temperature and  pressure as you do so. Measurement  Devices Temperature  measurement Platinum   resistance  thermometers  are......

Words: 806 - Pages: 4

Free Essay

Fluid Mechanics

...Fluid Mechanics Learning Objectives Outcomes • Explain the pressure-depth relationship. Pressure increases with depth. • Define Pascal’s Principle. Pascal's Principle states that the pressure is transmitted evenly through a liquid. • Describe how to use Pascal’s Principle in practical application. When you inflate a balloon with air, it expands evenly in all directions, this is an example. • Describe Archimedes Principle. States that the mass of a liquid displaced by a floating body is equal to the mass of that body. • Determine if an object will float in a fluid based on its relative densities. So if you fill a tumbler up with water to the brim, put an object into it, weigh the water that has been pushed out of the tumbler, and compare that with the weight of the object, you'll know whether it floats or not. • Use the continuity equation and Bernoulli’s equation to explain common effects of ideal fluid flow. The pressure in a fluid moving steadily without friction or outside energy input decreases when the fluid velocity increases Assignment Requirements 3. Mass is the same, so if the whale is taking up less volume, the density must have increased. The whale has displaced a greater mass of water at the depth, so the buoyant force is greater. 20. Ice cubes float in water, and sink in alcohol. Anything with less density than the liquid that it's in will float. 22. It will increase 35. It would be harder on the top of a mountain because the pressure of the......

Words: 342 - Pages: 2

Premium Essay

Engg Mechanics

...SRI SHANMUGHA COLLEGE OF ENGINEERING AND TECHNOLOGY Pullipalayam,Morur(P.O),Sankari(T.k),Salem(D.T). Two Mark Questions Unit I – Basics 1. What is meant by mechanics? Mechanics is a branch of physical science which deals with the study of a body or bodies such as machines and structures at rest or in motion subjected to external mechanical disturbances such as forces, moments etc. What is meant by Engineering mechanics? Application of the principles of science of mechanics to the practical engineering problems is known as Engineering Mechanics. State the different types of mechanics? Depending upon the nature of the body involved, Mechanics can be classified into two types * Mechanics of Solids * Mechanics of Fluids Define Statics The study of a body which is in motion is known as statics Define Dynamics. The study of a body which is in motion is known as dynamics. Define Kinematics. It is the branch of dynamics which deals with the relationship between displacement, velocity, acceleration and time of a given motion, without considering the forces that cause the motion. Define Kinetics It is the branch of dynamics which deals with the relationship between the forces acting on a body, the mass of the body and the motion of the body. What do you understand from the concept of “Law of dimensional homogeneity”? Law of dimensional homogeneity states that “basic equation representing physical phenomenon must be valid for all systems of units”. State Parallelogram law. It states......

Words: 5082 - Pages: 21

Free Essay

Technology in Mechanics

...Technology in Mechanics Nowadays, there has been a significant advancement of technology in mechanics. These advancement made it possible for human being to walk on Moon (kurianjoseph1, 2013). This feature, in turn had created a fact as the powerful of technology.Various technologies surrounding us are helping people to improve their life with more luxury. So the using of technology in mechanics impacts to us positively. There are one cause and two effects why the technology is using. Mechanics is the basis of early and modern technology. According to a study by (B.Cotterell & J.Kamminga, 1992), the advancement of technology in mechanics plays a fundamental role in the study of material culture. In simple term, mechanics can explain how artifacts operate technologically and how efficient or the effective they are. For example, a cooling system is a technology in a mechanics. In relation to this, (V.A.W.Hillier, 2004) stated that the function of cooling system is to remove heat from the mechanics’ engine for keeping their temperature within safe limits and avoid any problems. With the engine combustion reaching extremely high temperature, the heat needs to be dissipated. A cooling system is essential to prevent the engine from burning. Another function of the cooling system is to regulate the temperature inside the passenger compartment and ensure the passengers’ comfort. So, it can be seem that the use of technology in mechanics has become an...

Words: 613 - Pages: 3

Premium Essay

Mechanics

...UNIVERSITI TUNKU ABDUL RAHMAN |Centre |: Centre for Foundation Studies (CFS) |Unit Code |: FHSC1014 | |Course |: Foundation in Science |Unit Title |: Mechanics | |Year/ Trimester |: Year 1 / Trimester 1 | | | |Session |: 201401 | | | Tutorial 4: Application of Newton’s Laws. 1. The distance between two telephone poles is 50.0 m. When a 1.00 kg bird lands on the telephone wire midway between the poles, the wire sags 0.200 m. Draw a free-body diagram of the bird. How much tension does the bird produce in the wire? Ignore the weight of the wire. [614 N] 2. A 40 kg crate rests on a horizontal floor, and a 75 kg person is standing on the crate. Determine the magnitude of the normal force that (a) the floor exerts on the crate and (b) the crate exerts on the person. [(a) 1.13 x 103 N, (b) 735 N] 3. A worker stands still on a roof sloped at an angle of 45° above the horizontal. He is prevented from slipping by a static frictional force of 450 N. Find the mass of the worker. [85 kg] 4. A 4.0-kg bucket of water is raised from a well by a rope. If the upward acceleration of the......

Words: 937 - Pages: 4

Premium Essay

Fluid Mechanics

...Variation of Pressure vertically in a fluid under gravity P_2 A P_1 Force due to P_1 on area A acting up = P_1 A Force due P_2 on area A acting down = P_2 A Force due to the weight of the element =mg= ρA(z_2- z_1 )g Since the fluid is at rest, there can be no shear forces and hence no vertical forces on the side of element due to surrounding fluid. Considering upward as positive, sum of all forces = 0, We have P_1 A-P_2 A - ρA(z_2- z_1 )g=0 Or P_2- P_1= - ρ(z_2- z_1 )g So, under the influence of gravity , pressure decreases with the increase height. Equality of Pressure at the same level in a static fluid Since the fluid is at rest, there are no horizontal shear stresses on the sides of the element. For static equilibrium the sum of horizontal forces must be zero, P_1 A=P_2 A P_1=P_2 Therefore the pressure at any two points at the same level in a body of fluid at rest will be the same If the (x,y) is the horizontal plane then ∂p/∂x=0 and ∂p/∂y=0 Variation of pressure due to gravity from point to point in a static fluid Since the element is in equilibrium, the resultant of the forces in any direction should be zero Resolving forces in the axial direction, pA-(p+δp)A-ρgAδs cos⁡θ=0 δp=-ρgδs cos⁡θ Putting above equation in differential form, dp/ds=-ρgcos......

Words: 609 - Pages: 3

Premium Essay

Fluid Management

...Fluid Management Leadership Learning Experience A1. Problem or Issue Research shows that dialysis patients who have problems with fluid management have an increase in hospitalizations, disease processes, and poor clinical outcomes.  Research has also proven that fluid is a strong predictor of mortality and morbidity. A1a. Explanation of Problems or Issues We have noticed in our clinic a trend of increasing fluid overloaded patients over the last few months.  This trend has also resulted in multiple hospitalizations for congestive heart failure, pulmonary edema, and respiratory distress resulting from fluid overload.  This problem was selected after a patient expired in our clinic with cardiac arrest after a long code attempting to revive.  This topic of discussion made a big impact in our dialysis clinic because it is an open area that can be seen by all the patients.  Seeing a desperate attempt to save a life, it is not soon forgotten in the minds of a room full of people who all share the same disease.  The fear of what they have just seen could almost be read on their faces, most thinking that could have been me.  After every code, there is a review of what could have been done differently, as far as the code was concerned the process was correct.  However, looking back you wonder if there were something more you could have done to prevent the code altogether.  Not only was the patient someone whom I was given the responsibility of caring for, he was also......

Words: 2110 - Pages: 9

Premium Essay

Fluid

...Nurs 2820 Fluid, Electrolyte, and Acid-Base Balance Name: _______________________________________ Case Study Jimmy Lewis is brought to the hospital emergency room by some friends. He had been vomiting for several days and was complaining of heart palpitations. Mr. Lewis is a 58-year-old white male who is homeless. He has not had any health care for at least 10 years. He is an alcoholic and drinks a quart of gin or vodka every day. He does not have a job, and his family is all out of state. The emergency physician does an initial assessment and transfers him to a hospitalist, who admits him to a medical-surgical unit for further evaluation and treatment. Mr. Lewis has lab work drawn. His electrolytes are as follows: sodium 138 mEq/L, potassium 3.1 mEq/L (low), chloride 104 mEq/L, and magnesium 1.5 mEq/L (low). His arterial blood gas measurements are as follows: pH 7.48 (high), PaCO2 40 mm Hg, HCO3 29 (high). Jamie Taylor, a 22-year-old nursing student, is assigned to Mr. Lewis. She reviews Mr. Lewis’ medical record before going in to assess him. 1. After reviewing his chart and lab work, what fluid and electrolyte imbalances would Jamie determine? (Select all that Apply) A. Fluid volume deficit B. Hypokalemia C. Hypermagnesemia D. Hyperkalemia E. Hypomagnesemia 2. What acid-base imbalance is Mr. Lewis experiencing? A. Metabolic acidosis B. Respiratory acidosis C. Metabolic alkalosis D. Respiratory alkalosis 3. The hospitalist orders an IV of D5NS......

Words: 338 - Pages: 2

Free Essay

Mechanic

...UNIVERSITI TUNKU ABDUL RAHMAN CentreCourseYear/ Trimester | : Centre for Foundation Studies (CFS): Foundation in Science: Year 1 / Trimester 1 | Session Course CodeCourse Title | : 201605: FHSC1014: Mechanics | Tutorial 1: Introduction 1. What are the dimensions for the physical quantity expressed as , where m, g, d, T and l are the mass, acceleration due to gravity, distance, period and length respectively. A MT2 B ML2 C MLT2 D MLT–2 [Answer: B] 2. Suppose A = BC, where A has the dimension LM–1 and C has the dimension LT–1. What is the dimension of B? A TM–1 B L2T–1M–1 C TML–2 D L2TM–1 [Answer: A] 3. (a) State the dimension for (i) velocity (ii) force (b) Suppose an equation related with force F, viscosity ƞ, speed v, radius r is given by this formula , where k is a dimensionless constant; x, y, and z are the power to which ƞ, v and r. Given that the dimensions of viscosity ƞ is ML–1T–1. Use dimensional analysis to determine the value of x, y and z. [Answer: (a) (i) LT–1, (ii) MLT–2; (b) 1, 1, 1] 4. (a) Discuss briefly whether in general the method of dimensional checking can positively confirm that an equation is correct. (b) Newton’s law of universal gravitation is represented by where F is the gravitational force, M and m are masses, and r is a length. Given that the dimensions of force is MLT–2. What are the dimensions and SI unit of the proportionality...

Words: 770 - Pages: 4

Premium Essay

Mechanics

...UNIVERSITI TUNKU ABDUL RAHMAN |Centre |: Centre for Foundation Studies (CFS) |Unit Code |: FHSC1014 | |Course |: Foundation in Science |Unit Title |: Mechanics | |Year/ Trimester |: Year 1 / Trimester 1 |Lecturer | | |Session |: 201605 | | | Additional Tutorial 2: Vector and translational kinematics. 1. Find the x and y-components of: (a) a displacement of 200 km, at 30.0o. (b) a velocity of 40.0 km/h, at 120o; and (c) a force of 50.0 N at 330o. [(a) 173 km, 100 km, (b) -20.0 km/h, +34.6 km/h, (c) 43.3 N, -25.0 N] 2. Three forces are applied to an object, as indicated in the drawing. Force [pic] has a magnitude of 21.0 Newton (21.0 N) and is directed 30.0° to the left of the + y axis. Force [pic] has a magnitude of 15.0 N and points along the + x axis. What must be the magnitude and direction (specified by the angle ( in the drawing) of the third force [pic] such that the vector sum of the three forces is 0 N? [18.7 N, 76o] 3. A 200 N block rests on a 30o inclined plane as shown in figure below. If the weight of the block acts vertically downward, what are the components of......

Words: 989 - Pages: 4

Free Essay

Fluid Dynamics

...FLUID DYNAMICS In physics, fluid dynamics is a subdiscipline of fluid mechanics that deals with fluid flow—the natural science of fluids (liquids and gases) in motion. Fluid dynamics is "the branch of applied science that is concerned with the movement of liquids and gases," according to the American Heritage Dictionary. Fluid dynamics is one of two branches of fluid mechanics, which is the study of fluids and how forces affect them. (The other branch is fluid statics, which deals with fluids at rest.)  Scientists across several fields study fluid dynamics. Fluid dynamics provides methods for studying the evolution of stars, ocean currents, weather patterns, plate tectonics and even blood circulation. Some important technological applications of fluid dynamics include rocket engines, wind turbines, oil pipelines and air conditioning systems. FLOW The movement of liquids and gases is generally referred to as "flow," a concept that describes how fluids behave and how they interact with their surrounding environment — for example, water moving through a channel or pipe, or over a surface. Flow can be either steady or unsteady. In his lecture notes, "Lectures in Elementary Fluid Dynamics" (University of Kentucky, 2009) J. M. McDonough, a professor of engineering at the University of Kentucky, writes, "If all properties of a flow are independent of time, then the flow is steady; otherwise, it is unsteady." That is, steady flows do not change over time. An example......

Words: 1253 - Pages: 6

Free Essay

Fluids

...shall begin by studying incompressible flow problems. Of course all fluids are, to some extent, compressible but under steady flow conditions we may assume that the effects of changes in fluid density are small. In fact, it is the velocity of the fluid that dictates whether changes in density are significant and must be accounted for. In Chapter 6 we shall quantify the velocity limit, below which may assume that the fluid is incompressible; however, the majority of fluid flow problems that you are likely to encounter may be assumed to be incompressible. We shall focus in this chapter on incompressible flow, and on problems in which the fluid is bounded by a surface (we shall call this internal flow); the next chapter will focus on unbounded (or external) fluid flow problems. Both chapters will study real fluid flows and do this by taking into account the effects of viscosity. To do this we must examine how fluids interact with boundaries and here the concept of zero fluid velocity on a surface (boundary) is important. Once we have an understanding of how real fluid flows behave – and see how difficult it is to analyses turbulent flows –then in Chapters 4 and 5 we shall turn our attention to modelling techniques useful for examining simple fluid flow problems commonly found in engineering. Accordingly, this chapter will look at laminar and turbulent bounded fluid flows. We shall focus on pipe flow as this represents a classic......

Words: 270 - Pages: 2

Free Essay

Fluid Mechanics

...Fundamentals of Fluid Mechanics Fluid mechanics is the study of fluids and the forces on them. (Fluids include liquids, gases, and plasmas.) Fluid mechanics can be divided into fluid kinematics, the study of fluid motion, and fluid dynamics, the study of the effect of forces on fluid motion, which can further be divided into fluid statics, the study of fluids at rest, and fluid kinetics, the study of fluids in motion. Fluid mechanics is very important to engineers when observing flow in pipes, viscous effects of fluids, and the forces that act on a fluid. As a student, I am suppose to demonstrate an adequate understanding of many properties involved fluid mechanics. Some learning outcomes that must be accomplished by taking this class are: * Demonstrate understanding of fluid mechanics fundamentals, fluid and flow properties such as compressibility, viscosity, buoyancy, hydrostatic pressure and forces on surfaces * Apply Bernoulli equation to solve problems in fluid mechanics * Solve fluid mechanics problem using control volume analysis using conservation of mass, energy equation and irreversible flow * Use differential analysis of fluid flow, potential flow theory, viscous flow, Navier Stokes equations to solve problems * Perform modeling and similitude using Buckingham Pi theorem, correlation of experimental data. * Analyze flow in pipes to determine laminar and turbulent flow behaviors. * Apply energy and momentum equations to......

Words: 1681 - Pages: 7

Free Essay

Fluids

...CVEN 3100: Fluid Mechanics Fluid Properties: Review Questions 1. What is the definition of a fluid? A substance that deforms continuously when acted on by a shearing stress of any magnitude. 2. Normal force per unit area in a fluid is called what? Pressure 3. True or False: - Static fluids are not subjected to shear force at any time. T 4. True or False: - Normal forces can occur in a fluid whether it is static or Moving T 5. What is the relation between absolute pressure and gage pressure? Absolute pressure can be found from the gage pressure by adding the value of the atmospheric pressure. 6. What formula is used to calculate density of gases? Identify the parameters in the formula p/RT 7. Define specific weight. What is relation to density? Weight per unit volume. Multiply by gravity 8. Because of viscosity, what happens when a fluid tries to flow? It resists and does not flow quick 9. What is kinematic viscosity? The ratio of absolute viscosity to density 10. State the Newton’s law of viscosity and express it mathematically. Change in velocity over distance which velocity changes. Du/dy 11. What is the purpose of lubricating metal hinges? 12. Why does viscosity of liquids decrease with temperature? Molecules are spread further apart 13. Why does viscosity of gases increase with temperature? Molecular activity increases 14. What is an ideal fluid? A fluid that......

Words: 543 - Pages: 3

Premium Essay

Cfd in Fluid Mechanics

...the earth is roughly 7000˚C. The two reasons why the earth is hot is when the earth was formed, the interior was heated rapidly due to the gravitational forces being converted into heat and radioactive isotopes within the earth liberate heat as they continue to decay. Our ability to drill into the earth is restricted to the upper few kilometers of the earth’s crust. We must look for a location where the earth’s interior heat is brought within our reach. This is most common at plate boundaries. Direct heating systems are designed to supply hot water only with no electricity generation. Borehole drilled to depth of 1800 meters beneath the city of Southampton, UK. Near the bottom of the hole is water at 70˚C. The fluid contains dissolved salts. The fluid is more accurately described as Brine. The Brine is pressurized and so rises unaided to within 100 m of the surface. A turbine pumps the Brine up to the surface through a heat exchanger that transfers heat to clean water. There are four main types of installations, Dry Steam Power Plant, Single Flash Steam Power Plant, Binary Cycle Power Plant, and Double Flash Power Plant....

Words: 343 - Pages: 2