Fluid Mechanics

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 Classiﬁcation (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 diﬀraction by a pair of cylinders, by John Bell, Phillip Colella, William Crutchﬁeld, 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 ﬂuid 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 ﬂuid mechanics, nor to assess the engineering value of various approximation procedures. The goals were: • to present some of the basic ideas of ﬂuid 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 diﬃcult 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 ﬁrst 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 ﬂow, 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 ﬂow 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 speciﬁc 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. Birkhoﬀ [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 Diﬀerential 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 ﬁrst edition. We also thank O. Hald and P. Arminjon for a careful proofreading of the ﬁrst 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 ﬂuid 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 ﬂuid. These assumptions are relaxed in the third section to allow for viscous eﬀects 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 ﬁlled with a ﬂuid. Our object is to describe the motion of such a ﬂuid. Let x ∈ D be a point in D and consider the particle of ﬂuid 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 ﬂuid; this particle traverses a well-deﬁned trajectory. Let u(x, t) denote the velocity of the particle of ﬂuid that is moving through x at time t. Thus, for each ﬁxed time, u is a vector ﬁeld on D, as in Figure 1.1.1. We call u the (spatial ) velocity ﬁeld of the ﬂuid. For each time t, assume that the ﬂuid has a well-deﬁned mass density ρ(x, t). Thus, if W is any subregion of D, the mass of ﬂuid 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 ﬂowing 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 ﬂuid 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 ﬁxed 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 deﬁned at points of ∂W ; and let dA denote the area element on ∂W . The volume ﬂow rate across ∂W per unit area is u · n and the mass ﬂow 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 diﬀerential 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 diﬀerential 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 ﬂuid particle, so that the velocity ﬁeld 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 ﬂow).1 u(x(t), t) = The acceleration of a ﬂuid 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: ﬁrst, 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 ﬁnal 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 ﬂuid is moving and that the positions of ﬂuid 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 ﬁeld, 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 deﬁne an ideal ﬂuid as one with the following property: For any motion of the ﬂuid there is a function p(x, t) called the pressure such that if S is a surface in the ﬂuid 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 ﬂuid as a mathematical deﬁnition 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 ﬂuids 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 ﬂuid, nor, if it is there at the beginning, to stop. This idea will be ampliﬁed in the next section. However, even here we can detect physical trouble for ideal ﬂuids because of the abundance of rotation in real ﬂuids (near the oars of a rowboat, in tornadoes, etc.). If W is a region in the ﬂuid at a particular instant of time t, the total force exerted on the ﬂuid 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 ﬁxed 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 ﬂuid material, force per unit volume = −grad p + ρb. By Newton’s second law (force = mass × acceleration) we are led to the diﬀerential 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 ﬁrst as a deduction from the diﬀerential form and second from basic principles. From balance of momentum in diﬀerential 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 ﬁxed 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 ﬁxed 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 ﬂux 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 diﬀerential law. With an eye to assuming as little diﬀerentiability as possible, it is useful to proceed to the integral law directly and, as with conservation of mass, derive the diﬀerential form from it. To do this carefully requires us to introduce some useful notions. As earlier, let D denote the region in which the ﬂuid 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 ﬁxed t, ϕ is an invertible mapping. Let ϕt denote the map x → ϕ(x, t); that is, with ﬁxed t, this map advances each ﬂuid particle from its position at time t = 0 to its position at time t. Here, of course, the subscript does not denote diﬀerentiation. We call ϕ the ﬂuid ﬂow map. If W is a region in D, then ϕt (W ) = Wt is the volume W moving with the ﬂuid . 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 ﬂuid 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 ﬂuid in W ﬂow for time t.

where J(x, t) is the Jacobian determinant of the map ϕt . Because the volume is ﬁxed at its initial position, we may diﬀerentiate 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 diﬀerentiation following the ﬂuid.) Next, we learn how to diﬀerentiate 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 deﬁnition of the velocity ﬁeld of the ﬂuid. The determinant J can be diﬀerentiated by recalling that the determinant of a matrix is multilinear in the columns (or rows). Thus, holding x

1.1 Euler’s Equations

9

ﬁxed 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 diﬀerential 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 ﬂow incompressible if for any ﬂuid 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 ﬂuid 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 ﬂuid is incompressible if and only if Dρ/Dt = 0, that is, the mass density is constant following the ﬂuid . If the ﬂuid is homogeneous, that is, ρ = constant in space, it also follows that the ﬂow is incompressible if and only if ρ is constant in time. Problems involving inhomogeneous incompressible ﬂow occur, for example, in oceanography. We shall now “solve” the equation of continuity by expressing ρ in terms of its value at t = 0, the ﬂow 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 ﬂuid that is homogeneous at t = 0 but is compressible will generally not remain homogeneous. However, the ﬂuid 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 ﬂow 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 ﬁve functions: u, ρ, and p. Thus, one might suspect that to specify the ﬂuid motion completely, one more equation is needed. This is in fact true, and conservation of energy will supply the necessary equation in ﬂuid 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 conﬁne ourselves to two special cases here, and later we shall treat another case for an ideal gas. For ﬂuid moving in a domain D, with velocity ﬁeld 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 ﬂuid 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 ﬂuid or if we allow the ﬂuid to do work, Etotal will change. The rate of change of kinetic energy of a moving portion Wt of ﬂuid 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 ﬂuid 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 ﬂuid 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 ﬂow 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 ﬂow: 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 ﬁrst law is a statement of conservation of energy; a statement equivalent to (TD1) is, as is readily veriﬁed, d = T ds + p dρ. ρ2 (TD2)

If the pressure is a function of ρ only, then the ﬂow 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 satisﬁes 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 ﬂows with p a function of ρ, the integral form of energy balance reads as follows: The rate of change of energy in a portion of ﬂuid 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 ﬂow 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 ﬂuid 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 ﬂuid, 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 ﬂuid ﬂow with velocity ﬁeld u(x, t), a streamline at a ﬁxed 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 satisﬁes dx = u(x(s), t), t ﬁxed. ds We deﬁne a ﬁxed 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 diﬀerential 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 ﬂow is called stationary. Bernoulli’s Theorem In stationary isentropic ﬂows 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 ﬂow 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 ﬂow 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 ﬂuid-ﬁlled 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 ﬂow in a channel.

Suppose that the pressure p1 at x = 0 is larger than that at x = L so the ﬂuid 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 ﬂow with a constant pressure gradient increases indeﬁnitely. Of course, this cannot be the case in a real ﬂow; 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 ﬂow without any body force. Show that for a ﬁxed volume W in space (not moving with the ﬂow). d dt
1 2ρ

u

2

+ρ

dV = −
∂W

ρ

W

1 2

u

2

+ w u · n dA.

Use this to justify the term energy ﬂux 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 ﬁeld of a ﬂuid is u = (u, v, w), then its curl,

is called the vorticity ﬁeld of the ﬂow. We shall now demonstrate that in a small neighborhood of each point of the ﬂuid, u is the sum of a (rigid ) translation, a deformation (deﬁned later ), and a (rigid ) rotation with rotation vector ξ/2. This is in fact a general statement about vector ﬁelds u on R3 ; the speciﬁc features of ﬂuid 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 ﬁxed, an orthonormal basis e1 , e2 , e3 in which D is diagonal:   d1 0 0 D =  0 d2 0  . 0 0 d3 Keep x ﬁxed and consider the original vector ﬁeld as a function of y. The motion of the ﬂuid is described by the equations dy = u(y). dt If we ignore all terms in (1.2.1) except D · h, we ﬁnd dy = D · h, dt i.e., dh = D · h. dt

20

1 The Equations of Motion

This vector equation is equivalent to three linear diﬀerential 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 ﬁeld 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 conﬁrms the fact proved in §1.1 that volume elements change at a rate proportional to div u. Of course, the constant vector ﬁeld u(x) in formula (1.2.1) induces a ﬂow that is merely a translation by u(x). The other term, 1 2 ξ(x) × h, induces a ﬂow dh = 1 ξ(x) × h, 2 dt (x ﬁxed).

The solution of this linear diﬀerential 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 ﬂuid at t = 0. Let Ct be the contour carried along by the ﬂow. 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 ﬂuid ﬂow map discussed in §1.1 (see Figure 1.2.1). The circulation around Ct is deﬁned to be the line integral

ΓC t =
Ct

u · ds.

Kelvin’s Circulation Theorem For isentropic ﬂow without external forces, the circulation, ΓCt is constant in time. For example, we note that if the ﬂuid 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 ﬁeld of a ﬂow and C a closed loop, with Ct = ϕt (C) the loop transported by the ﬂow (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 deﬁnition 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 ﬁrst 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 ﬂow is isentropic and without external forces), we ﬁnd 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 ﬂux of vorticity across a surface moving with the ﬂuid 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 ﬂow.

By deﬁnition, 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 ﬂow of an isentropic ﬂuid 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 ﬂux 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 ﬂow (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 ﬂow (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 ﬂow map.

and ω(ϕ(x, t), t) = ∇ϕt (x) · ω(x, 0), (1.2.7) where ϕt is the ﬂow 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 ﬁrst-order diﬀerential 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 ﬂow using (1.2.7) and (1.2.10). For two-dimensional ﬂow, 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 ﬂuid, 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 ﬂow. Employing (1.2.10) and the change of variables theorem gives (1.2.11) as a special case. In three-dimensional ﬂows, 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 ﬂuid mechanics it is customary to be sloppy about this deﬁnition and make tacit assumptions to the eﬀect that the tube really “looks like” a tube. More precisely, we assume S is diﬀeomorphic to a disc (i.e., related to a disc by a one-to-one invertible diﬀerentiable transformation) and that the resulting tube is diﬀeomorphic 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 ﬂuid 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 ﬂuid. 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 ﬂuid. It either forms a ring (such as the smoke from a cigarette), extends to inﬁnity, 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 ﬂuid. 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 ﬂuid is clearly false if we allow ξ to have zeros and probably is false even if ξ has no zeros (an orbit of a vector ﬁeld 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 diﬀerence 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 ﬂuid 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 ﬂow. For two-dimensional homogeneous incompressible ﬂow, the vorticity equation is Dξ (1.2.12) = ∂t ξ + (u · ∇)ξ = 0, Dt where ξ = ξ(x, y, t) = ∂x v −∂y u is the (scalar) vorticity ﬁeld of the ﬂow and u, v are the components of u. Assume that the ﬂow is contained in some plane domain D with a ﬁxed 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 ﬁxed 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 ﬂow 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 ﬂow. 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 ﬂow 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 deﬁned as above is a solution of the two-dimensional stationary incompressible equations of ideal ﬂow.

30

1 The Equations of Motion

is

For three-dimensional incompressible ideal ﬂow, 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 ﬁeld 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 ﬁeld. Exercise 1.2-2 Couette ﬂow. Let Ω be the region between two concentric cylinders of radii R1 and R2 , where R1 < R2 . Deﬁne 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 ﬂow on the two cylinders is ω1 and ω2 .

1.3

The Navier–Stokes Equations

In §1.1 we deﬁned an ideal ﬂuid as one in which forces across a surface were normal to that surface. We now consider more general ﬂuids. To understand the need for the generalization beyond the examples already given, consider the situation shown in Figure 1.3.1. Here the velocity ﬁeld 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 ﬂuid 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 diﬀuse across S and impart momentum to the ﬂuid below, and, likewise, slower molecules from below S will diﬀuse across S to slow down the ﬂuid above S. For reasonably fast changes in velocity over short distance, this eﬀect is important.5
B' u u B
Figure 1.3.1. Faster molecules in B can diﬀuse across S and impart momentum to B.

S

We thus change our previous deﬁnition. 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 deﬁnite 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 ﬂuid 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 σ modiﬁes 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 ﬂuid undergoes a rigid body rotation, there should be no diﬀusion 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 deﬁnes 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 ﬁrst coeﬃcient of viscosity , and ζ = λ + 2 µ is the 3 second coeﬃcient 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 ﬂow of a compressible viscous ﬂuid. In the case of incompressible homogeneous ﬂow ρ = ρ0 = constant, the complete set of equations becomes the Navier–Stokes equations for incompressible ﬂow, Du = − grad p + ν∆u Dt div u = 0

(1.3.4)

where ν = µ/ρ0 is the coeﬃcient of kinematic viscosity, and p = p/ρ0 . These equations are supplemented by boundary conditions. For Euler’s equations for ideal ﬂow we use u · n = 0, that is, ﬂuid 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 ﬂuid 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 ﬂow9 . In any case, adding the tangential boundary condition is crucial for viscous ﬂow. The physical need for the extra boundary conditions comes from simple experiments involving ﬂow past a solid wall. For example, if dye is injected into ﬂow 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 diﬀusion. 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 ﬂuid. Another crucial feature of the boundary condition u = 0 is that it provides a mechanism by which a boundary can produce vorticity in the ﬂuid. 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 eﬀect of viscosity on the ﬂow. 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 ﬂow past a sphere, L could be either the radius or the diameter of the sphere, and U could be the magnitude of the ﬂuid velocity at inﬁnity. L and U are merely reasonable length and velocity scales typical of the ﬂow 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 deﬁne the Reynolds number R to be the dimensionless number LU R= . ν For example, consider two ﬂows past two spheres centered at the origin but with diﬀering radii, one with a ﬂuid where U∞ = 10 km/hr past a sphere of radius 10 m and the other with the same ﬂuid 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 inﬁnity, then the Reynolds number is the same for each ﬂow. The equations satisﬁed by the dimensionless variables are thus identical for the two ﬂows. Two ﬂows with the same geometry and the same Reynolds number are called similar . More precisely, let u1 and u2 be two ﬂows 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 ﬂow, and let the viscosities be ν1 and ν2 respectively. If L1 U1 L2 U2 = , R1 = R2 , i.e., ν1 ν2 then the dimensionless velocity ﬁelds 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 ﬂows 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 ﬂuid ﬂow 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, coeﬃcient of viscosity, and so on, such that the Reynolds number for the ﬂow in our experiment matches that of the actual ﬂow. We can then expect the results of our experiment to be relevant to the actual ﬂow 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 eﬀects 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 ﬂow, the pressure p in incompressible viscous ﬂow is determined through the equation div u = 0. We now shall explore the role of the pressure in incompressible ﬂow 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 ﬁeld 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 deﬁned 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, deﬁne 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 ﬁeld 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 deﬁned. 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 satisﬁed 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 ﬂows is conceptually diﬀerent than in incompressible ﬂows just as it was in ideal ﬂow. If we think of viscous ﬂow as ideal ﬂow with viscous eﬀects 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 deﬁned here (through equation (1.3.1)) is identical to p as deﬁned in that other science. Not all quantities called p are equal. The use of expressions from equilibrium thermodynamics requires an additional physical justiﬁcation, which is indeed often available, but which should not be forgotten. According to the analysis given earlier, the pressure p in incompressible ﬂow 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 ﬂow with p = p(ρ), where p (ρ) > 0. If ﬂuid ﬂows into a given ﬁxed 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 ﬂuid to ﬂow 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 ﬂuid, 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 diﬀusion or dissipation term,

which are the Stokes’ equations for incompressible ﬂow. 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 ﬂows 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 eﬀect, 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 diﬀerence between the ideal and viscous ﬂow with regard to the energy of the ﬂuid. 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 ﬂow, we should have d Ekinetic ≤ 0. dt (1.3.9)

40

1 The Equations of Motion

We calculate (d/dt)Ekinetic for incompressible viscous ﬂow 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 ﬂow 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 ﬂow in a channel leads to unreasonable results. We now reconsider this example with viscous eﬀects. Example Consider stationary viscous incompressible ﬂow 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 ﬂuid 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 ﬂuid is pushed from left to right and correspondingly, p1 > p2 .

Because each side depends on diﬀerent 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 proﬁle 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 ﬂuid to achieve a stationary state. We saw that this was not possible for ideal ﬂow. Next we consider the vorticity equation for (homogeneous) viscous incompressible ﬂow. In the two-dimensional case we proved in §1.2 (see equation (1.2.12)) that for isentropic ideal plane ﬂow, Dξ/Dt = 0. The derivation is readily modiﬁed to cover viscous incompressible ﬂow; the result is Dξ 1 = ∆ξ. (1.3.11) Dt R This shows that the vorticity is diﬀused by viscosity as well as being tranported by the ﬂow. 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 ﬂow 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 ﬂow, the vorticity equation is Dξ 1 − (ξ · ∇)u = ∆ξ. (1.3.13) Dt R Thus, vorticity is convected, stretched, and diﬀused. (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 ﬂow, 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 ﬂow allows for the generation of vorticity.

1.3 The Navier–Stokes Equations

43

This is possible because of the diﬀerence in boundary conditions between ideal and viscous ﬂows. The mechanism of vorticity generation is related to the diﬃculties 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 ﬂows for simplicity. Assume that we have an equation of state p = p(ρ), Deﬁne c= p (ρ). For reasons that will become clear later, c is called the sound speed of the ﬂuid. Thus, we have c2 dρ = dp. (1.3.14) Let u = u be the ﬂow speed. One calls M = u/c the (local) Mach number of the ﬂow; it is a function of position in the ﬂow. From Bernoulli’s theorem proved in §1.1, u2 + 2 dp = constant on streamlines. ρ(p) (1.3.15) p (ρ) > 0.

Also, diﬀerentiating 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 ﬂow map. Combining (1.3.14), (1.3.15), and (1.3.16) we get dJ du = −M . J c The ﬂow will be approximately incompressible if J changes only by a small amount along streamlines. Thus, a steady ﬂow can be viewed as incompressible if the ﬂow 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 ﬂow will be approximately incompressible if γ is very large. For nonsteady ﬂow one also needs to know that l τ c,

where l is a characteristic length and τ is a characteristic time over which the ﬂow picture changes appreciably.13 The presence of viscosity does not alter these conclusions signiﬁcantly.

Exercises
Exercise 1.3-1 Find a stationary viscous incompressible ﬂow 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 inﬁnite pipe.

(i) Poiseuille ﬂow . 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 ﬂow 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 ﬂow 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 ﬂow 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 ﬂows. We begin with a more detailed study of inviscid irrotational ﬂows, that is, potential ﬂows. Then we go on to study boundary layers, where the main diﬀerence between slightly viscous and inviscid ﬂows originates. This is further developed in the third section using probabilistic methods. For most of this chapter we will study incompressible ﬂows. A detailed study of some special compressible ﬂows is the subject of Chapter 3.

2.1

Potential Flow

Throughout this section, all ﬂows 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 ﬂow with zero vorticity is called irrotational . For ideal ﬂow, this holds for all time if it holds at one time by the results of §1.2. An inviscid, irrotational ﬂow is called a potential ﬂow . 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 ﬂow 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 ﬂow 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 ﬂow 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 ﬂow, note that w = p/ρ0 from the deﬁnition of w. For potential ﬂow 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-deﬁned and is constant in time. Next, consider incompressible potential ﬂow 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 ﬂow is potential in Σ.

Let the velocity of ∂D be speciﬁed 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 ﬂow (satisfying (2.1.4)) in D if and only if ∂D V · n dA = 0; this ﬂow is the minimizer of the kinetic energy function Ekinetic = 1 2 ρ u
D 2

dV,

among all divergence-free vector ﬁelds 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 ﬂow in a bounded region with ﬁxed boundary is the trivial ﬂow u = 0. For unbounded regions this is not true without a careful speciﬁcation of what can happen at inﬁnity; the above uniqueness proof is valid only if the integration by parts (i.e., use of the divergence theorem) can be justiﬁed. 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 ﬂow on the simply connected

2.1 Potential Flow

51

y D streamlines

x

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

region D shown in Figure 2.1.2. This ﬂow may be arrived at by the methods of complex variables to which we will now turn. Incompressible potential ﬂow is very special, but is a key building block for understanding complicated ﬂows. For plane ﬂows 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 deﬁne an incompressible (stationary) potential ﬂow. 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 ﬂow 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 ﬁxed 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 ﬂow 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 identiﬁed with a complex number in the standard way; i.e., (x, y) is identiﬁed with z = x + iy. Proof If dz = dx + i dy represents an inﬁnitesimal 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 ﬂow exterior to a region B. Let the velocity ﬁeld approach the constant value (U, V ) = U at inﬁnity. 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 inﬁnity, 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 inﬁnity, 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 ﬂow is normal to the direction of ﬂow and is proportional to the circulation around the body. In any case, the body experiences no drag (i.e., no force opposing the ﬂow) in contradiction with our intuition and with experiment. The diﬃculty, 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 ﬂow. 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 ﬁeld is u = (U, V ). This is two-dimensional ﬂow 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 ﬂow 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 ﬂow 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 ﬂow around a disc.

From (2.1.12) with z = aeiθ , that is, z on ∂B, we ﬁnd 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 deﬁned in the whole plane, then ˜ W (z) = W (z) + W a2 , z |z| ≥ a

is a potential describing a ﬂow exterior to the disc of radius a > 0, but possibly with a more complicated behavior at inﬁnity. 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 ﬂow that is incompressible and has vorticity ξ = −∆ψ. If we can arrange for ψ to be the imaginary part of an analytic function, then the ﬂow 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 ﬁeld is zero at inﬁnity. For incompressible potential ﬂow about a disc of radius a centered at z0 , we need only choose W (z) = Γ log(z − z0 ). 2πi

The boundary conditions are satisﬁed because ψ is constant on any circle centered at z0 (see Figure 2.1.6). The incompressible potential ﬂow 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 ﬂow 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 ﬂow separately, it is also true for W given by (2.1.15). Thus, we get an incompressible potential ﬂow on the exterior of the disc |z| ≤ a with circulation Γ around the disc. The velocity ﬁeld is (U, 0) at inﬁnity (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 deﬁned 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 ﬂow around an obstacle in three dimensions with constant velocity U at inﬁnity, not only can there be no drag, there can be no lift either.

58

2 Potential Flow and Slightly Viscous Flow

The diﬀerence 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 ﬁnite 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 ﬁrst term in the expansion in powers of 1/r is now missing. For an incompressible potential ﬂow there will be a potential ϕ, that is, u = grad ϕ (because the exterior of the body is simply connected). The potential satisﬁes ∆ϕ = 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 outﬂow 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 ﬂow 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 ﬂow 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 inﬁnity. We leave the detailed veriﬁcation to the reader.1 Next we will discuss a possible mechanism, ultimately to be justiﬁed by the presence of viscosity, by which one can avoid d’Alembert’s paradox. An eﬀort to resolve the paradox is of course prompted by the fact that real bodies in ﬂuids do experience drag. By an almost potential ﬂow , we mean a ﬂow in which vorticity is concentrated in some thin layers of ﬂuid; the ﬂow is potential outside these thin layers, but there is a mechanism for producing vorticity near boundaries. For example, one can postulate that the ﬂow past the obstacle shown in Figure 2.1.8 produces an almost potential ﬂow with vorticity produced at the boundary and concentrated on two streamlines emanating from the body. We image diﬀerent potential ﬂows in the two regions separated by these streamlines with the velocity ﬁeld discontinuous across them. For such a model, the Kutta–Joukowski theorem does not apply and the drag may be diﬀerent from zero. There are a number of situations in engineering where real ﬂows 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 ﬂow 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 ﬂow around them are outside the scope of this book. Next we shall examine a model of incompressible inviscid ﬂow inspired by the idea of an almost potential ﬂow and Example 3. We imagine the vorticity in a ﬂuid is concentrated in N vortices (i.e., points at which the vorticity ﬁeld 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 ﬂow generated by point vortices in the plane.

As the ﬂuid moves according to Euler’s equations, the circulations Γj associated with each vortex will remain constant. The vorticity ﬁeld 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 satisﬁes

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 ﬁeld 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 ﬁeld
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 ﬁeld 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 ﬂows 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). Deﬁne 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 satisﬁed. If C is a contour containing l vortices at x1 , x2 , . . . , xl , then l ΓC = − i=1 Γi and ΓC is invariant under the ﬂow. However, the relationship between these solutions and bona ﬁde solutions of Euler’s equations is not readily apparent. Such a relationship can, however, be established rigorously and such vortex systems do contain signiﬁcant 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. Deﬁne 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 inﬁnite. 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 ﬂuids, 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 ﬂow 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 ﬂow 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 reﬂection principle: G(x, xj ) = 1 ˆ (log x − xj + log x − xj ) , 2π

ˆ where xj = (xj , −yj ) is the reﬂection 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 diﬀerential 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 ﬂow 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 speciﬁc 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 ﬁeld 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 deﬁned by the above integral satisﬁes 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 ﬁeld ξ is concentrated on C only, that is, the ﬂow is potential outside the ﬁlament 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 ﬂow induced by a vortex ﬁlament.

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 ﬁrst 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 ﬁnd a formula for potential ﬂow 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 ﬂuid in R3 . For slow velocities, assume Stokes’ equations apply. Take the point of view that the object is stationary and the ﬂuid streams by. The setup for the boundary value problem is as follows: given U = (U, 0, 0), U constant, ﬁnd 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 inﬁnity. 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 outﬂow at inﬁnity (an inﬁnite wake). Exercise 2.1-5 Because of the diﬃculties 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 ﬂow 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 diﬀerence 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 diﬀerent; 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 inﬁnite! 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 diﬀerent a ﬂow governed by (2.2.1) is from one governed by the Euler equations for incompressible ideal ﬂow:  ∂t u + (u · )u = − grad p,   div u = 0, (2.2.2)   u · n = 0 on ∂D. Imagine that both ﬂows coincide at t = 0 and, say, are irrotational, that is, ξ = 0. Thus, under (2.2.2) the ﬂow stays irrotational, and thus is a potential ﬂow. However, we claim that the presence of the (small) viscosity term (1/R) u and the diﬀerence in the boundary conditions have the following eﬀects: 1. The ﬂow governed by (2.2.2) is drastically√ modiﬁed near the wall in a region with thickness proportional to 1/ R. 2. The region in which the ﬂow is modiﬁed may separate from the boundary. 3. This separation will act as a source of vorticity.
5 See Birkhoﬀ’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 eﬀects via 2 and 3 may persist and not diminish in the limit R → ∞. A model for the ﬁrst eﬀect 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) diﬀers 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 veriﬁed 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 diﬀerent is conﬁned 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 diﬀerence 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 ﬂow 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 ﬂat plate.

We seek a ﬂow 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 deﬁned 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 inﬁnity, 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 ﬁxed t > 0 is obtained by rescaling the above graph. The region in which u departs signiﬁcantly from the constant Euler ﬂow is again called the boundary layer and is illustrated in the ﬁgure. By √ √ the preceding scaling argument, this thickness is proportional √ to √ νt = t/ R . Thus, for ﬁxed time, the boundary layer decreases as 1/ R , which was our ﬁrst contention. Our second contention was that the modiﬁed ﬂow in the boundary layer may separate from the body. We can give a heuristic argument justifying this contention by considering potential ﬂow 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 eﬀect, points of high potential energy. As a ﬂuid 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 ﬂow. Within this layer the ﬂuid is slowed down by friction. Therefore, ﬂuid 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 ﬂuid 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 ﬂows really can undergo boundary layer separation. Our third contention was that the boundary layer produces vorticity. To see this, imagine ﬂow 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 speciﬁcally 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 ﬂow 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 ﬂow 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 eﬀects 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 ﬂow should be eﬀectively inviscid by deﬁnition (we are focusing on the region where the viscous eﬀects are important). We claim that v is small at a distance δ from the wall. Indeed, at the wall, v = 0 and so inside the ﬂuid, 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 eﬀects 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 eﬀects 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 eﬀects 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 ﬁelds. 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 diﬀused vertically into the ﬂuid. 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 ﬂow, 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 ﬂow past the ﬂat 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 ﬂow ξ = −∂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 ﬂat. 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 modiﬁed 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 ﬁxed curvature of the wall, the pressure change across the boundary layer is negligible. However, (2.2.10)2 may still be signiﬁcant in aﬀecting 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 ﬁeld 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.

Diﬃculties On the line y = δ the ﬂows may not be parallel, so matching u or some other ﬂow 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. Diﬃculties The calculations are diﬃcult, and the results are insuﬃciently 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.

Diﬃculties 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. Diﬃculties Same matching problems as before, but also modiﬁed 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 diﬃculties 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 proﬁle near the point of separation. Contemplation of this velocity ﬁeld 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 ﬂow 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 ﬂow on a ﬂat plate of inﬁnite width. Consider a ﬂat 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 ﬂuid lies in the domain given by y > 0 and the ﬂow is two dimensional. We then let the plate be washed by a steady uniform ﬂow of velocity U in the positive x-direction. We assume that a steady boundary layer is developed after a short distance oﬀ 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 ﬂow 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 ﬁnd 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 ﬂuid 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 ﬂow over a ﬂat 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 ﬂow described by the Euler equations near the boundary to simulate the eﬀects 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 eﬀect 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 deﬁned 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 deﬁned 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 deﬁnition, 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 Tchebysheﬀ ’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 inﬁnite 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 diﬀerential 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 ﬁnd 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 ﬁnal 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 ﬁxed. 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 reﬂected 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 halﬂine 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∗ satisﬁes 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 ﬂow past an inﬁnite ﬂat plate at rest. In §2.2 (see Figure 2.2.3), we saw that the ﬂow 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 deﬁne a vortex sheet of intensity ξ as follows. As in Figure 2.3.2, it consists of a ﬂow 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 ﬂow 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 ﬂow (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 “reﬂect” 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 eﬀort.) 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 ﬂow will be approximated at t = 0 by a collection of N vortex sheets of ﬁnite 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 ﬂow ∂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 ﬂow 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 ﬂat plate example (including reﬂections) 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 ﬂow. The velocity ﬁeld satisﬁes (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 ﬁeld 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 ﬁnd 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 ﬁrst 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 ﬁeld is now determined by the same vortex sheets, but at their new positions. This new velocity ﬁeld satisﬁes 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 diﬀerent 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 ﬂow. 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 ﬂow together with a random y component simulating viscous diﬀusion. Vortex sheets newly created to accommodate the boundary conditions diﬀuse out from the boundary by means of the random component and then get swept downstream by the main ﬂow. 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 ﬂow past a semi-inﬁnite ﬂat 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 ﬂow 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 ﬂow 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 ﬂow is stationary; we may regard t as a parameter and we may eliminate it to ﬁnd vertical displacement ∼ horizontal displacement.

This shows that the structure of the boundary on a semi-inﬁnite ﬂat 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 diﬀerential equations; we can add a random term to take care of the diﬀusion 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 ﬁlaments, 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 ﬁrst a diﬀerential 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 eﬃciency 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 deﬁne a map Ht : Rn → Rn , called the ﬂow map. Let Kt and Lt

2.4 Remarks on Stability and Bifurcation

95

be the corresponding ﬂow 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 ﬂow of the Stokes equation and Lt the ﬂow 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 ﬂows 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 modiﬁed 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 ﬁeld, u = 0, on the boundary ∂Ω of the region Ω in question; • Lt is the ﬂow of the Euler equation (boundary conditions u parallel to ∂Ω), which may be solved approximately by a vortex method; • Kt is the ﬂow 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 ﬁeld to u whose backﬂow cancels u on ∂Ω; and • Ht is the ﬂow 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 diﬀerence 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 ﬂow 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 ﬂow. The other extreme, the limit R → ∞ concerns very fast or only

96

2 Potential Flow and Slightly Viscous Flow

slightly viscous ﬂow. The situation for intermediate values of R involves many interesting transitions, or changes of ﬂow 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 ﬂow of a diﬀerential 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 ﬁxed point of X. It follows from the initial condition F0 (x0 ) = x0 and uniqueness of solutions that x0 is a ﬁxed point of X, i.e., Ft (x0 ) = x0 . Deﬁnition A point x0 is an asymptotically stable ﬁxed point of the vector ﬁeld 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] Diﬀerential 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 ﬁxed point of a smooth autonomous vector ﬁeld 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 suﬃciently close to x0 , the term DX(x0 ) dominates the behavior of the ﬂow 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 ﬁeld Xµ (x, y) = (y, µ(1 − x2 )y − x) on R2 . Determine if the ﬁxed 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 ﬂow 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 ﬁxed point of a vector ﬁeld 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 ﬁelds satisfying the appropriate boundary conditions. The Navier–Stokes equations determine the dynamics on this space P of velocity ﬁelds. 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 diﬀerential 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 diﬀerentiation and DX(v0 ) is a linear diﬀerential operator. Its spectrum will consist of inﬁnitely 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, ﬂow in a pipe and Couette ﬂow are stable (in fact globally stable) if the Reynolds number is not too big. Sometimes one is interested in the loss of stability of ﬂows as the Reynolds number R is increased. In general, any such qualitative change in the nature of a ﬂow 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 ﬁeld (depending on a parameter µ) possessing a ﬁxed point x0 for all µ. Assume the eigenvalues of DXµ (x0 ) are in the left half-plane of C for µ < µ0 , where µ0 is ﬁxed. 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 ﬂuid mechanics this means that a stationary solution (ﬁxed point) of the Navier– Stokes equations is replaced by a periodic solution (stable closed orbit). For example, consider the case of ﬂow 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 ﬂuid 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 ﬁrst-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 ﬂow, 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 ﬂows from this approach, but for some complex ﬂows like TaylorCouette ﬂow the approach has been very successful.

14 See Marsden and Hughes [1994]; M. Golubitsky, I. Stewart and D. Schaeﬀer 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 ﬂow in one dimension. In the ﬁrst 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 ﬂow 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 ﬁnal section we generalize the discussion to the ﬂow of a gas that allows chemical energy release, such as occurs in combustion.

3.1

Characteristics

One-dimensional isentropic ﬂow 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. Deﬁne c = p (ρ) (assuming p (ρ) > 0) and note that px = p (ρ)ρx by the chain rule. Thus, the ﬁrst 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 ﬁrst 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 diﬀerentiating the ﬁrst 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 deﬁne 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 ﬁrst-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 reﬂected 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 diﬀerential equations. Let us verify that, as in the ﬁrst example, solutions are constant along characteristics. Indeed, if x = x(s), t = t(s) is a characteristic and v satisﬁes (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 coeﬃcient 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 deﬁne 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 oﬀ 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 deﬁned 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 diﬀerent points can indeed intersect (see Figure 3.1.3). It is also easy to see in principle why the nonlinear case diﬀers from the linear case. In the linear case, the characteristics are determined by two ordinary diﬀerential 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 deﬁne 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 deﬁne a characteristic to be a curve with the following property: if data are given on that curve, the diﬀerential 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 deﬁnition.)

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 deﬁnition. 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 oﬀ 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 ﬁnd 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 diﬀerentiate 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 ﬁnd 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 ﬁnds 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 diﬃculties with two C+ characteristics crossing. There may also be boundary conditions that have to be taken into account. We now discuss a diﬀerent 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 coeﬃcient 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 deﬁnes the characteristics; that is, we recover the Riemann invariants obtained earlier. As we saw earlier, it may be diﬃcult to determine the characteristics in a general ﬂow. 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 deﬁned analogously. In a Γ+ simple wave, both Γ− and Γ+ are constant along C− characteristics by deﬁnition 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 inﬁnitely long tube extending along the x-axis. We assume that gas ﬁlls 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 rareﬁed 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 ﬂow at every point (x, t). Two constant states, S1 and S2 , are said to be connected by a Γ− rarefaction wave if the conﬁguration shown in Figure 3.1.8 is a solution of the diﬀerential equations for gas ﬂow. 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 conﬁguration 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 ﬁnd ρ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 ﬂow.1 We assume there is a speciﬁc internal energy function (“speciﬁc” 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 ﬂuid domain from its boundaries, the total energy can only be reduced when the ﬂuid does work. The work done by a ﬂuid 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 diﬀerential 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 satisﬁed, as long as u and ρ are smooth. The energy equation (3.2.1) is called the ﬁrst 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 speciﬁc 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 diﬀerentiability to justify the diﬀerential 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 justiﬁed 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 deﬁne a weak solution of (3.2.3) to be a function u that satisﬁes (3.2.4) for all smooth ϕ with compact support. A slightly diﬀerent 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 diﬀerential 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 diﬀerentiable. Note, however, that a function that satisﬁes the integral form of the equation of motion does not have to be diﬀerentiable either, and the integral forms in ﬂuid 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, aﬃrmative 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 ﬁnd 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 ﬁgure. 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 diﬀerential form of the equation is satisﬁed whenever u is smooth, and (3.2.8) is satisﬁed across any discontinuity. It is, of course, assumed that all discontinuities are jump discontinuities. A function u that satisﬁes the diﬀerential equations where possible and (3.2.8) across a jump discontinuity satisﬁes 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 ﬂow: 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 ﬁnd 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 ﬁrst two equations in (3.2.11) are called the mechanical jump relations. The equations of ﬂuid 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 ﬁxed 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 diﬀerent 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 ﬂuid. 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 speciﬁc 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. Deﬁne 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 deﬁned 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 ﬁll 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 ﬁgure shows two ways to ﬁnd a globally deﬁned 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 satisﬁes 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. Speciﬁcally, 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 eﬀect 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 deﬁne 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 ﬁnite 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 ﬂow, 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 satisﬁes 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 satisﬁes 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 satisﬁed.) 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 ﬂuid particle crosses the shock. We next show that for a γ-law gas a shock is compressive if and only if it satisﬁes 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 deﬁnition, ∗ ρ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 ﬂuid 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 conﬁguration 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 deﬁnition. 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 diﬀerentiation, 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 conﬁguration 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 satisﬁed. 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 ﬁrst 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 conﬁrms 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. Deﬁne 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 ﬁnd pl , ul such that Sl is connected to Sr by a centered wave. If pl > pr , we can ﬁnd a straight line shock connecting the two states, and if pl ≤ pr , we can ﬁnd 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 ﬁnite 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 deﬁnitions 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 ﬁnd 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 deﬁnition 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. Deﬁne 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 ﬁrst 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 speciﬁc 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 deﬁnition 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 deﬁne 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, deﬁne Λ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 deﬁne 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 deﬁned in §3.2) of the initial data f (x) is suﬃciently 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 ﬁxed 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 ﬂow 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 ﬂow equations, other methods of proof are available.9 A major component in the Glimm construction is the solution of the Riemann problem. We have oﬀered 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 modiﬁed 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 speciﬁc system. We examine in four steps the cases where the graph of f is straight, concave up, concave down, and, ﬁnally, 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. Diﬀerential 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 conﬁguration shown in Figure 3.4.3. We get a shock conﬁguration 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 , deﬁne 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 coeﬃcient matrix A(u) when our equations are written in hyperbolic form (3.1.3). For gas dynamics we have deﬁned the families C+ , C− , C0 of characteristics. We say that a family of characteristics crosses a discontinuity that satisﬁes 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 satisﬁes 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 modiﬁcation 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 eﬀects, heat conduction, and radiative heat transfer. Of these the most serious is the exclusion of the eﬀects 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 simpliﬁed situation. We shall furthermore assume that the region in which combustion occurs (and in which q changes from q0 to q1 ) is inﬁnitely 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 ﬂow with ﬁnite 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 ﬂow. 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 deﬁne 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 diﬀers 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 deﬂagration 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 deﬂagration branch, the part IV below CJ2 the strong deﬂagration 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 ﬂow 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 ﬁgure 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 suﬃcient 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 ﬂow 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 ) ﬁrst 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 ﬂow 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 deﬂagration 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 deﬂagration 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 ﬁnite 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 deﬂagration that can occur is one across which the pressure is constant and there is no mass ﬂow, that is, this deﬂagration 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 deﬂagration. Case 5 S lies in the strong deﬂagration 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 deﬂagration branch.

Then the gas ﬂow relative to the reaction front is subsonic in the front side and supersonic in the back side. If this strong deﬂagration moves in the positive x direction, it separates the C+ characteristics. However, this separation does not satisfy the geometric entropy condition, although the strong deﬂagration is consistent with the conservation laws. Thus, we must exclude strong deﬂagrations. 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 ﬂow,” 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 deﬁne 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 ﬁrst 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 ﬂow, 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 ﬂow 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 ﬁelds 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 ﬁeld ∇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 deﬁne (F × ∇) × G = U × G

Vector Identities

159

where we deﬁne 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 ﬂuid particle, 4 algorithm, 94 almost potential ﬂow, 59 asymptotically stable, 97 autonomous, 97 back, 124 balance of momentum, 4, 6, 12 diﬀerential 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 ﬂow, 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 coeﬃcients of viscosity, 33 combustion, 103, 152 front, 152 wave, 145 complex potential, 51 variables methods, 51 velocity, 51, 53 compressible

162

Index

ﬂow, 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 ﬂow, 31 Courant–Friedrichs–Lewy condition, 141 cylindrical coordinates, 45 d’Alembert’s paradox, 54 in 3d, 58 decompostion theorem, 37 deﬂagration branch, 154 deformation, 18 tensor, 19, 31 degenerate linearly, 151 density, 1, 14 derivative material, 5 detonation branch, 154 diﬀerential form, 120 form of mass conservation, 3 diﬀusion, 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 ﬂux, 18 internal, 12, 117 kinetic, 12 per unit volume, 117 enthalpy, 14 entropy, 14, 118, 158 condition, 127, 130, 133, 157 equation continuity, 3 diﬀerential 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 ﬁeld velocity, 1 vorticity, 18 ﬁlament, 65 Filon’s paradox, 67

Index

163

ﬁrst coeﬃcient of viscosity, 33 law of thermodynamics, 14, 118 ﬁxed point, 97 ﬂat plate, 69 ﬂow, 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 ﬁlament, 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 ﬂuid ﬂow map, 7 ideal, 5 particle, 4 velocity, 1 viscous, 33 ﬂux, 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 ﬂow, 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 ﬂow, 69 Hamiltonian system, 62 heat equation, 84, 86 Helmholtz decomposition theorem, 37 theorem, 26, 37 Hodge theorem, 37 hodograph transformation, 110 homogeneous, 11 ﬂow, 48 Hopf bifurcation theorem, 99 Hugoniot curve, 126, 153 equation, 125, 139 function, 125, 153 hyperbolic, 104 ideal ﬂow, 48 ﬂuid, 5, 31

164

Index

gas, 44, 118, 125, 131, 139 ignition temperature, 152 incompressible, 11 approximately, 44 ﬂow, 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 ﬂow, 14, 15 ﬂuids, 14 gas ﬂow, 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 ﬂow rate, 46 matching solutions, 78 material derivative, 5, 18 mean, 82 mechanical jump relations, 124 momentum balance of, 6 ﬂux, 7 moving with the ﬂuid, 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 ﬂow 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 ﬂow, 45 piston, 111 plate, 80 ﬂow between, 41 ﬂow over, 66, 69

Index

165

ﬂow past, 87 point vortices, 60 Poiseuille ﬂow, 45 potential complex, 51 ﬂow, 47, 48, 51, 56, 58 almost, 59 ﬂow around a disk, 55 velocity, 48 vortex, 56 vortex ﬂow, 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 coeﬃcient 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 ﬂows, 36 simple wave, 111 simply connected, 47 skin drag, 80 slightly viscous ﬂow, 47 slip line, 124, 138, 140 sound speed, 44, 104, 131 space-time divergence, 119 spatial velocity ﬁeld, 1 speed discontinuity, 147 sphere ﬂow 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 ﬂow, 16, 58 ﬂow criterion, 29 steady ﬂow, 58 Stokes equation, 40, 67, 96 ﬂow, 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 deﬂagration, 154 detonation, 154 subsonic, 131, 157 supersonic, 131 symmetry, 101 tangential boundary condition, 34 forces, 5 Tchebysheﬀ’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 ﬁrst 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 ﬂow in, 51 variance, 83, 143 variation, 129 velocity characteristic, 35 complex, 51 ﬁeld, 1 potential, 48, 51 proﬁle, 42 viscosity, 129 coeﬃcients, 33 kinematic, 34 viscous ﬂuid, 33 vortex, 61 ﬁlament, 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 deﬂagration, 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 modiﬁed from the third corrected printing version typeset on 10-28-99.

Modiﬁcations:
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 ﬁgures 2.1.x.eps, x=3,4,5,6, and save as Illustrator 8 eps ﬁles. Apr-03-2000: Redo the frontmatter to include a Series Preface (Springer) page as page (v) and start the preface on (vii)

Similar Documents

Fluid Mechanics

Fluid Mechanics

Fluid Mechanics

Cfd in Fluid Mechanics

Written Report

3m Hvlp Spray Gun

Nt1310 Unit 4 Lab

Bible Among the Myths

Smooth Bore vs. Fog Nozzles

Fluid Dynamics

Civil

Fluids

Fluids

Decompression Speed Lab Report

System Design for Transport of Limestone

Popular Essays