14

Applications of the Theory of Partial Differential Equations to Problems of Fluid Mechanics

PAUL R. GARABEDIAN

PROFESSOR OF MATHEMATICS
NEW YORK UNIVERSITY

14.1    Introduction

In this chapter we shall be concerned with Cauchy’s problem for a linear or quasi-linear partial differential equation in two independent variables. We shall be interested in applications of this subject to fundamental, and also to less familiar, problems in mathematical physics.

The usual physical context in which one encounters Cauchy’s problem is that of a mechanical configuration, with known initial position and velocity, for which the subsequent motion is governed by a partial differential equation of the second order. The classical example is the case of the vibrating string, with behavior described by the wave equation^2,4

u_xx − u_tt = 0

Another example is the motion of a perfect inviscid gas in a pipe, which leads to a nonlinear partial differential equation. In these familiar cases, one of the two independent variables is the time t, and the Cauchy data, or initial data, are simply the values of the displacements and velocities of a mechanical system at the time t = 0. Thus it is truly an initial-value problem to which such physical situations lead, and, furthermore, the partial differential equations involved are of the hyperbolic type.² This is as it should be, since we can then deduce that the Cauchy problem, or initial-value problem, is properly set in the sense that the solution depends in a continuous manner on the given physical data, so that a small error in measurement of the latter would not lead to significant errors in the results.

Although the great bulk of the theory of Cauchy’s problem has been developed for hyperbolic partial differential equations, and although the direct physical applications of the theory involve hyperbolic equations, it is nevertheless true that Cauchy’s problem is of significance in the elliptic case as well. Hadamard has shown, of course, as stated in Chap. 4, that the problem is not properly set for an elliptic equation, in the sense that the solution does not in general depend continuously on the data. To illustrate this, one usually cites (cf. pages 111 to 112 of Ref. 10) the case of the Laplace equation

u_xx + u_yy = 0

with Cauchy data

u(x,0) = 0u_y(x, 0) = n⁻¹ sin nx

The expression

u(x, y) = n⁻² sin nx sh ny

is easily verified to be the solution of this particular Cauchy problem. The solution, however, does not approach zero as n → ∞, despite the fact that the Cauchy data do tend toward zero in this limiting process.

We have thus established without difficulty the incorrectness, in Hadamard’s sense, of the Cauchy problem for an elliptic partial differential equation, but we can also present physical applications in which it would be quite useful to find the solution even in this case. Such situations do not usually arise from the direct physical formulation of a problem, but they can occur in the construction, by inverse methods, of solutions possessing prescribed properties that would not ordinarily be part of a given set of physical data. Thus, with suitable idealizations, if we set out to describe an electromagnetic field guiding an electron beam in a prescribed channel, then we must discuss a Cauchy problem for Laplace’s equation. The direct, as opposed to the inverse, physical problem here would have been to locate the paths of the electrons when the field is given. Another question that results in Cauchy’s problem for an elliptic equation concerns the determination of flows of an incompressible, inviscid fluid possessing prescribed free streamlines. Finally, the inverse construction of steady gas flows through a given detached shock wave in front of a blunt body can be effected by solving Cauchy’s problem for a nonlinear partial differential equation of mixed elliptic-hyperbolic type. A common feature of all these examples is the appearance of an inverse solution to what is essentially a free-boundary problem.

From the list of applications we have mentioned, it becomes evident that a successful analytical treatment of the Cauchy problem is called for even in the case of an elliptic equation. The Cauchy-Kowalewski solution by means of power-series expansions clearly is not satisfactory here any more than it is for hyperbolic problems. What is actually required is a method that can be adapted for numerical integration based on a finite-difference scheme. We shall present such an approach in this chapter. We emphasize some of the numerical aspects of this topic because of the frequent occurrence in practice of situations that are so complicated that only such an analysis proves to be adequate. A treatment of Cauchy’s problem that lends itself to numerical computation is, of course, as important for the hyperbolic as for the elliptic case, and we shall bear this in mind in presenting our material.

We begin with a discussion of the method of characteristics for a hyperbolic partial differential equation in two independent variables. This is followed by the description of a finite-difference scheme for the numerical solution of Cauchy’s problem in a characteristic coordinate system. Next we take up elliptic equations and study them by means of analytic extension into the complex domain. This material is followed by an elementary application to a free-streamline flow problem concerned with a bubble that rises under the influence of gravity. We close the chapter with a brief treatment of an inverse solution of the detached-shock problem for supersonic flow past a blunt body.

14.2    Cauchy’s Problem for a Hyperbolic Partial Differential Equation in Two Independent Variables

We consider the quasi-linear partial differential equation

where the coefficients a, b, c, and d are given functions of the two independent variables x and y and of the three quantities u, u_x, and u_y. For the moment, we are interested primarily in the hyperbolic case of a solution u in a region where

Cauchy’s problem for Eq. (14.1) consists in finding a solution u that assumes on a given curve Γ prescribed values of u, u_x, and u_y. Of course, in order that u_x and u_y should be the first partial derivatives of the function u, the Cauchy data u, u_x, and u_y must fulfill along Γ the familiar condition

so that we actually prescribe there only two independent functions of a single variable.

A natural first step in any attempt to solve the Cauchy problem just described would be to try to calculate the second derivatives u_xx, u_xy, and u_yy of u along Γ. The partial differential equation (14.1) itself yields one linear relationship among these three unknowns, and two further linear equations of the form

can be found for them along Γ by differentiation there of the Cauchy data. We are assured of the existence of a unique solution of this system of three simultaneous linear equations in three unknowns if and only if the determinant

of the coefficients does not vanish along Γ. Curves Γ along which this determinant vanishes identically must be treated as inappropriate for the solution of the Cauchy problem, and they are known as the characteristics of the partial differential equation (14.1). Along a characteristic curve Γ, we can form a linear combination of Eqs. (14.4) and (14.5) that must reduce to Eq. (14.1), so that an ordinary differential equation is obtained there for the Cauchy data themselves. Even when the Cauchy data fulfill this additional condition along a characteristic, they do not determine uniquely the second or higher derivatives of u, so that in this case more than one solution of Cauchy’s problem can exist. A consequence of this remark is that a curve across which the higher derivatives of a solution u of Eq. (14.1) exhibit discontinuities must be a characteristic. Thus, for example, the characteristics of the equation of one-dimensional gas dynamics can represent sound waves.

The relationship

which defines the characteristics of Eq. (14.1), can be interpreted as an ordinary differential equation for a family of curves in the xy plane. Of course, when a, b, and c depend on u, u_x, or u_y, the same can be true of the characteristics. Under our assumption (14.2) of hyperbolicity, the quadratic equation (14.7) for the slope dy/dx has two distinct roots σ and τ, so that actually two separate one-parameter families of characteristics are defined. Incidentally, there is no real loss of generality here if we assume for the sake of simplicity that a ≠ 0, since when this is not true, it can easily be achieved by rotating the xy plane. Now it would seem natural to try to introduce our two families of characteristic curves as a new coordinate network. This can be done by defining two functions α(x,y) and β(x,y) such that one family of characteristics is represented in implicit form by the equation α(x,y) = const, while the other family is represented by the equation β(x,y) = const. Of course, α becomes a natural parameter to use along each curve of the latter family, whereas β is the parameter to be used along the characteristics of the first family.

Let us rewrite the quadratic ordinary differential equation (14.7) for the characteristics in the form

where σ and τ depend in an obvious fashion on a, b, and c, and let us suppose that it is the first factor in Eq. (14.8) that vanishes on the characteristics along which α is the parameter. If we think of α and β as new independent variables, it follows that

whereas by symmetry

If our partial differential equation (14.1) is linear, so that a, b, and c, and hence σ and τ, do not depend on u, u_x, or u_y, then the pair of partial differential equations (14.9) and (14.10) can be solved for x and y independently of Eq. (14.1). We are interested, however, in the quasi-linear case when σ and τ do involve u, u_x, and u_y, and therefore we shall reformulate the whole problem in terms of the characteristic coordinates α and β, and we shall study Eq. (14.1) simultaneously with Eqs. (14.9) and (14.10).

Our next step is to express the original partial differential equation (14.1) in terms of the new independent variables α and β. This is achieved by following through on our earlier remark that along a characteristic we can combine Eqs. (14.4) and (14.5) with Eq. (14.1) to derive an ordinary differential equation there for the Cauchy data. Let the characteristic in question be one along which α is the parameter, and introduce the familiar notation p = u_x, q = u_y. We notice that along our characteristic the determinant (14.6) vanishes, while Eqs. (14.1), (14.4), and (14.5) express the fact that a linear combination of the columns of this determinant with coefficients chosen to be u_xx, u_xy, and u_yy is equal to the triple of quantities −d, p_α, and q_α. Hence any three-by-three determinant from the matrix

must vanish. We consider, in particular, the determinant composed of the first, third, and fourth columns, and we divide it by x_α and use Eq. (14.9) to obtain

This is the ordinary differential equation to which we have been referring for the Cauchy data p and q along a characteristic with the parameter α variable. In a symmetric fashion we can derive along a characteristic of the opposite family, with β variable, the analogous equation

If u itself does not appear as an argument of the coefficients a, b, c, and d, then Eqs. (14.9) to (14.12) comprise an elegant canonical system of four first-order quasi-linear partial differential equations for the new unknowns x, y, p, and q as functions of the new independent variables α and β. The system is equivalent to Eq. (14.1), but it has the advantage that in each equation partial derivatives with respect to only one of the two independent variables appear. In the new formulation, the characteristics are simply the parallels to the coordinate axes. Thus, the characteristics and the solution of Eq. (14.1) are found simultaneously in the xy plane when we solve the canonical system (14.9) to (14.12) and then invert to introduce x and y once more as the independent variables.

If u occurs as an argument of the given functions a, b, c, and d, then an additional artifice is required in our derivation of the canonical system. This consists in choosing either of the two equations

or

as a fifth equation of the system. These latter equations are readily obtained from the requirement that p and q are the first derivatives of u with respect to x and y, and they are, indeed, equivalent to this requirement. It is an instructive exercise for the reader to verify that either of the two relationships (14.13) and (14.14) implies, in conjunction with the remainder of the system (14.9) to (14.12), that the other must hold for a solution of a Cauchy problem with data satisfying Eq. (14.3) along a non-characteristic initial curve Γ. Thus it is not difficult to obtain a canonical system of five first-order partial differential equations in characteristic coordinates for the five unknowns x, y, p, q, and u that can be used in place of the original Eq. (14.1) in the general quasi-linear case, but the final results are somewhat unattractive because of the lack of symmetry made necessary by our choice of either Eq. (14.13) or Eq. (14.14) as the final equation.

A considerable amount of freedom remains in our choice of a characteristic coordinate system, for we can always replace α by an arbitrary function of α, and β by an arbitrary function of β, in any of our canonical equations (14.9) to (14.14) without changing their form. For a given initial curve that is not characteristic, we can always pick the characteristic coordinates so that the initial curve corresponds to the line α + β = 0 in the αβ plane and so that on this line x = α − β, whereupon no further transformations or normalizations are possible. It is also instructive to use s = α − β and t = α + β as the independent variables and to reformulate our canonical system in terms of matrix notation. To do this, we set

and

where

It then turns out that the system (14.9) to (14.13) assumes the simple matrix form

The reader can easily verify as an exercise that the square B² of the matrix B is the identity matrix I. Thus also

If we differentiate Eq. (14.16) with respect to t and Eq. (14.17) with respect to s and subtract the latter from the former, we find

This second-order system of partial differential equations helps to motivate our reduction to a canonical form, since the wave-equation operator now appears on the left. One application that this formulation suggests, for example, is a solution of the Cauchy problem based on a sequence of successive approximations found by inserting an established approximation into the terms on the right in Eq. (14.18) and calculating the next approximation by substituting it into the operator on the left. The Cauchy data here would be given at t = 0.

We are now in possession of the more formal aspects of the theory of characteristics for the hyperbolic equation (14.1). We turn our attention to the matter of examples and applications. Perhaps the most instructive case that we could treat here is that of the isentropic one-dimensional flow of a perfect inviscid gas. Such a flow can be visualized as the motion of a gas in a very narrow rectilinear pipe. Let x denote a space coordinate measured along the pipe, let t denote the time, let ρ denote the density, and let v denote the velocity of the gas; v is actually a scalar, since the flow is only one-dimensional. The pressure P is given explicitly as a function of the density ρ alone by an equation of state

Conservation of mass provides us with the first-order partial differential equation

connecting the quantities ρ and v, while conservation of momentum yields the further equation

When the initial distributions of density ρ and velocity v are given in the pipe, the subsequent motion of the gas can be determined by solving the system of equations (14.19) to (14.21).

We can easily reformulate the initial-value problem just described as a Cauchy problem for a second-order partial differential equation of the type (14.1). Indeed, the first-order equation (14.20) asserts the existence of a function u of x and t such that

and according to Eqs. (14.19) and (14.21) this function must satisfy the second-order hyperbolic partial differential equation

u_x²u_tt − 2u_tu_xu_tx + [u_t² − u_x²P′ (u_x)]u_xx = 0

The coefficients in this equation do not depend on u itself, whence our equivalent canonical system in characteristic coordinates α and β reduces to the four first-order equations

where p = u_t and q = u_x and where for clarity we have replaced u_x and u_t in the coefficients by their values from Eqs. (14.22). The last two Eqs. (14.25) and (14.26) can also, of course, be written in the form

From Eqs. (14.23) and (14.24) we see that at each point in the tx plane the two slopes of the characteristics through that point are given by

We recall that the characteristics can be interpreted in our time-space diagram as the possible trajectories of infinitesimal disturbances or sound waves. Equation (14.29) shows that, relative to the motion of the gas with speed v, these trajectories advance along the pipe with the velocity ± P′(ρ)^½. Thus we recognize that, at each point, the quantity P′ (ρ)^½ represents the local speed of sound.

If we introduce the expression

we can put Eqs. (14.27) and (14.28) in the form

Consider a situation in which the initial distributions of velocity v and density ρ in the pipe are so arranged that v − R has a constant value there at the time t = 0. By integrating Eq. (14.30), we conclude that v − R must then be constant throughout the subsequent flow. But if we now advance along any specific characteristic on which β varies, Eq. (14.31) tells us that v + R also remains constant, whence both v and ρ must have constant values there. Hence from Eq. (14.24) the slope of this characteristic does not vary, and it must actually be a straight line. A flow of this special type, for which one family of the characteristics in the tx plane consists exclusively of straight lines, is called a simple wave. Consider now the case in which the slopes dt/dx of these lines increase as we move along our pipe in the direction of increasing x. As time goes on, the lines must reach an envelope, and this envelope must form a singularity in the flow. Thus by a consideration of the geometry of the characteristics alone we have been able to establish that discontinuities must eventually appear in the flow; see Ref. 2, pages 101 to 105. These discontinuities are shock waves; as is well known, their theory³ is extremely important in fluid mechanics.

A useful exercise for the reader is to show that, if ρ and v are introduced as the independent variables, the resulting partial differential equations obtained for t and x are linear. Riemann⁹ has used this remark to derive an elegant explicit solution of Cauchy’s problem for the equation of one-dimensional gas dynamics in the case of a polytropic gas with the equation of state

P = Aρ^γ

where A is a positive constant and γ, the ratio of specific heats, is a constant exceeding 1. His results provide one of the few examples where an exact solution of a nonlinear problem in the mechanics of continua can be achieved. The details of his method, however, would carry us too far from our present discussion.

14.3    The Method of Finite Differences

Cauchy’s problem for the hyperbolic partial differential equation (14.1) reduces, through our transformation to the characteristic coordinates α and β, to an initial-value problem for the canonical system (14.9) to (14.13) in which the values of x, y, p, q, and u are prescribed, for example, along the line α + β = 0. We are interested in setting up an approximate solution of this problem by means of the method of finite differences. This method consists in laying down a grid of mesh points over the αβ plane, in replacing the derivatives in the canonical system by suitable differences of function values at the mesh points in order to obtain a discrete set of equations, and in solving these equations step by step for the function values involved. There are many possible formulations of such a technique that differ from one another in the type of grid that is used and in the type of differences that are selected to approximate the derivatives. We shall describe here some of the simpler difference schemes that are in common use for the solution of initial-value problems.

First, let us consider a square grid over the αβ plane, having mesh-point coordinates of the form

α = mhβ = nh

where m and n are arbitrary integers and where h is a fixed positive constant called the mesh size of the grid. We have been given the values of the unknowns x, y, p, q, and u at the grid points such that m + n = 0, and our objective is to calculate step by step, from the values of the unknowns at grid points with a fixed choice of m + n, the corresponding values at grid points with m + n one unit greater. We use the subscripts m, n to indicate quantities evaluated at the mesh point (mh, nh); thus, for example, we write

x_m,n = x(mh,nh)

With this notation, a useful difference approximation to the canonical system (14.9) to (14.13) is given by the equations

This set of equations is so arranged that it is possible to solve for the unknowns x, y, p, q, u at the grid point (mh, nh) directly in terms of their values at the previous pair of grid points ((m − 1)h, nh) and (mh, (n − 1)h). Such a formulation enables us to march forward, calculating the values of the unknowns at successive levels of mesh points in the desired fashion. The rules of the procedure are altogether elementary and straightforward, but in principle one should also justify the convergence of the answer to the solution of the initial-value problem for the actual differential equations (14.9) to (14.13) as the mesh size h tends to the limit zero. We shall not go into this more difficult question in the present short chapter, where we must content ourselves with a heuristic discussion of the order of magnitude of the error in the difference scheme when the convergence is granted.

There are two kinds of errors that should be distinguished; see page 486 of Ref. 7. The first of these is the truncation error caused by approximating derivatives by finite differences, and the second is the rounding error that occurs at each step in the computation when we retain only a fixed number of decimal places in our arithmetic. As to the truncation error, we can estimate it crudely by noting that, as h → 0, a forward difference such as (x_{m, n}, − x_m−1,n)/h approaches the corresponding derivative x_α[(m − 1)h, nh] with a remainder of the order of magnitude h, provided we assume that the second derivative x_αα is uniformly bounded. Thus, substitution of such finite differences into our system of partial differential equations in place of the actual derivatives can be viewed as equivalent to the introduction of additional terms of the order of magnitude h. Such new terms might be expected to alter the solution by an equivalent amount, and therefore we can assume that the truncation error itself is of the order of magnitude h. It is significant in this connection that, if we replace the coefficients in the difference equations (14.32) to (14.36), such as σ_{m−1, n}, by the appropriate averages, such as (σ_{m, n} + σ_{m−1, n})/2, we achieve a symmetric arrangement more or less equivalent to the use of a central-difference approximation to the derivatives. An expansion in Taylor series, which the reader can easily carry out as an exercise, shows that this symmetric approximation is equivalent to introducing into the differential equations extra terms only of the order of magnitude h², whence the truncation error of the modified scheme should be proportional only to h² instead of h; since, however, the coefficients now have the unknowns at the advanced grid point

α = mhβ = nh

among their arguments, the new system of difference equations is in implicit form and can be solved only by a sequence of iterations; this fact, of course, is a disadvantage of the refinement.

The magnitude of the rounding error can be estimated by a similar technique. Let denote the maximum contribution we can omit in any single arithmetic operation by retaining only a certain fixed number of significant figures or decimal places in the computation. It is then clear that an error of the order of magnitude could appear in each of the difference equations (14.32) to (14.36) because of rounding. Notice, however, that these equations must be divided by h before the differences appearing in them occur in a form that approximates the derivatives in our system of differential equations. The rounding error should therefore be reviewed as arising from insertion in the differential equations of additional terms proportional to /h. We conclude that the rounding error itself is of the order of magnitude /h. It follows that, to obtain a significant numerical solution of the Cauchy problem, we must always take much smaller than h. Therefore, as the mesh size h is taken smaller and smaller, it is necessary to retain more and more decimal places in the calculations.

In our numerical treatment of the Cauchy problem by finite differences, we notice that a certain segment of Cauchy data along the initial line α + β = 0 determines the solution in a triangular region bounded by this segment and by two intersecting characteristics α = const and β = const through its end points; see Fig. 14.1. By letting the mesh size h → 0, we see that the same result is valid for our canonical system (14.9) to (14.13) or for the original hyperbolic partial differential equation (14.1) itself. This analysis of the precise amount of data upon which the solution of the Cauchy problem at any specific point must depend is one of the most important features of our theory. It is interesting and instructive now to manipulate with difference schemes that do not take this phenomenon into account but are designed instead to yield the solution in too large a triangle. We shall indicate how such a procedure leads to unstable results as the mesh size h diminishes.

Fig. 14.1 Finite-difference grid for Cauchy’s problem.

For this purpose we refer to the canonical system of partial differential equations (14.18). In the ts plane, we introduce the grid of mesh points

t = mhs = nk

with h and k not necessarily equal, and the corresponding difference equation is seen to be

in an obvious notation. For our analysis of the stability of this difference scheme, we shall confine our attention to the simplest case B ≡ 0 of the wave equation. Here we can find explicit solutions of Eq. (14.37) of the form

where T and S are appropriately selected constants. Now if |T| ≤ 1 whenever |S| = 1, we shall say that the difference scheme (14.37) is stable, since then a bounded initial choice of U_m,n at m = 0 corresponds to a bounded solution as m increases. On the other hand, if there exists a T with |T| > 1 for |S| = 1, we shall say that the difference equation (14.37) is unstable, since we can find a solution with bounded initial values that grows exponentially as m increases.

If we set , we see that Eq. (14.38) is a solution of Eq. (14.37) if and only if

a result that the reader can check for himself as an exercise. When h ≤ k, this shows that |T| = 1 for all real choices of θ, and Eq. (14.37) is therefore stable. The restriction h ≤ k on our grid in the ts plane means that a triangle of mesh points used in a step-by-step solution of Cauchy’s problem for the difference equation (14.37) lies within the corresponding triangle bounded by a pair of characteristics emanating from the ends of the segment of Cauchy data involved. Thus all is as it should be and we have not attempted to compute the solution of Eq. (14.18) from data insufficient to determine it. On the other hand, when h > k, there is a real choice of θ such that one of the two values of T given by Eq. (14.39) has a modulus larger than unity, |T| > 1. Hence the difference scheme (14.37) is unstable; indeed, in this case the triangle of mesh points where the difference equation has a solution depending on given Cauchy data reaches beyond the characteristic triangle in which the solution of the original partial differential equation (14.18) is determined by these data. Thus we must never use an h > k in numerical applications; actually, the best choice would be h = k, which corresponds more or less to the formulation of the finite-difference method that we presented relative to the αβ plane.

Of course, the difficulty with the difference scheme (14.37) in the unstable case h/k > 1 is that, when we pass to the limit h → 0, k → 0 with the mesh sizes in such a ratio, our solution can exhibit large oscillations in spite of any restrictions on its initial data. It is interesting at this point to remark that this phenomenon occurs for any ratio h/k of mesh sizes, no matter how small, if we attempt to treat by these methods Laplace’s equation

U_tt + U_ss = 0

instead of the hyperbolic equation (14.18). For this elliptic case, indeed, our whole discussion could proceed along the same lines except for a difference in sign in one place, but this change in sign would bring Eq. (14.39) into the new form

Here for any positive pair of mesh sizes h and k there is a real choice of θ such that |T| > 1, so that the difference scheme for the solution of Cauchy’s problem is always unstable. This result is consistent with the comments in Sec. 14.1 to the effect that Cauchy’s problem is incorrectly set for Laplace’s equation. Despite all this, however, we shall have occasion later on to develop a procedure that yields a valid numerical solution of Cauchy’s problem for elliptic equations in those cases that are of practical significance for physical applications.

14.4    Cauchy’s Problem in the Elliptic Case

We shall present here a theory of Cauchy’s problem for an elliptic partial differential equation in two independent variables with analytic coefficients and analytic data. Our treatment will be suitable for numerical analysis by the method of finite differences when the Cauchy data have an analytic continuation into the complex domain in the large. In Sec. 14.1, we have described applications of this technique to problems of plasma physics, stationary shock waves, and flow with free streamlines, but we shall begin by discussing a much simpler special example.

Consider the linear elliptic equation

with Cauchy data

where f and g are assumed to be real analytic functions of the real argument y. This problem has a solution u in the form of a convergent power series, and thus u can be extended analytically to a complex-valued function of the complex arguments

x = x₁ + jx₂y = y₁ + jy₂

by direct substitution of these variables into its series representation. By holding x₂ and y₁ fixed, we find from Eq. (14.40) and the Cauchy-Riemann equations that u = u₁ + ju₂ satisfies the telegraph equation

which is of hyperbolic type. In particular, we set x₂ = 0 and find that u = u(x₁, y₁ + jy₂) has as a function of x₁ and y₂ the Cauchy data

along the initial line x₁ = 0, where f(y₁ + jy₂) and g(y₁ + jy₂) are the analytic continuations into the complex y plane of the two functions f and g appearing in Eq. (14.41). Thus u satisfies the properly set Cauchy problem defined by Eqs. (14.42) and (14.43) when it is considered as a complex-valued function of the pair of real variables x₁ and y₂. This Cauchy problem can be treated either numerically by the method of finite differences or explicitly by means of Bessel functions. The original elliptic Cauchy problem (14.40) and (14.41) can therefore be solved in two stages consisting of the analytic continuation of the data into the complex y plane and of the solution by any standard procedure of the hyperbolic problem (14.42) and (14.43). To obtain the final real answer, it suffices to set y₂ = 0 once more and to take u = u₁; but, of course, to calculate this result in a full region of the real xy plane, we must repeat the solution of the hyperbolic problem for a complete interval of values of the parameter y₁.

In the foregoing analysis, the unstable dependence of the solution of the elliptic equation (14.40) on the real Cauchy data (14.41) appears exclusively in the step where we extend the data to complex values of the independent variable y. Therefore, in cases when this analytic continuation of f and g can be performed in an elementary way—e.g., by direct substitution of complex values of y into explicit expressions for f and g—we encounter no difficulty in solving the elliptic Cauchy problem in the manner described. Furthermore, it is evident from inspection of our construction (see Fig. 14.2) that the region of existence of the solution in the real xy plane is identical with the region of regularity of the analytic functions f(y) and g(y) in the y₂y₁ plane.

Fig. 14.2 Geometry of Cauchy’s problem in the complex domain.

Our next concern will be the more complicated situation of Cauchy’s problem for a quasi-linear elliptic equation (14.1) with analytic coefficients a, b, c, and d satisfying

For analytic data, a convergent power-series solution still exists, and we can again extend it into the complex domain by direct substitution of complex arguments into the series expansion. In this fashion we have at hand an analytic function u of two complex variables x and y, and there is no reason why we cannot perform with it all the calculations that led us earlier from the partial differential equation (14.1) to the canonical systems (14.9) to (14.13) or (14.16). All the manipulations required here are just as valid for analytic functions of several complex variables as they are for differentiable real functions of real arguments. This reduction to canonical form will prove to be an essential feature of our treatment of Cauchy’s problem in the complex domain, where it no longer makes sense to speak of the sign of the discriminant b² − ac and where, in general, the roots σ and τ of the quadratic equation (14.7) are complex numbers.

If we think of u as known, we can determine a general solution of the ordinary differential equation (14.7) for the characteristics in the implicit form

solved for the arbitrary constant of integration. This can be done, for example, by working out the power-series solution of Eq. (14.7) that assumes the value y = y₀ at x = x₀ and then expressing y₀ as a function of the remaining variables in the result. Since

according to Eq. (14.45), the function α = α₁ + jα₂ satisfies the first-order partial differential equation

and this is true in particular for strictly real values of x and y. In this real domain, a, b, and c are real, whence also satisfies Eq. (14.46) there. Thus if we define

β(x,y) = α₁(x,y) − jα₂(x,y)

for real values of x and y and extend it afterward to complex values by analytic continuation, we obtain a second solution of Eq. (14.46) valid in the entire complex domain. Evidently, the implicit equation

β(x,y) = const

defines another family of solutions of the quadratic ordinary differential equation (14.7). The pair of functions α and β that we have introduced in this manner form a possible choice of the characteristic coordinates in the canonical system (14.9) to (14.13), and they have the special property that for real x and y the expression t = α + β is real, while s = α − β is pure imaginary.

We set s = ξ + jη and introduce t and η as new real independent variables, whereupon the canonical system (14.16), which is equivalent to Eq. (14.1), takes the form

The point here is that for ξ = 0, and for real choices of t and η, the components x, y, p, q, and u of the solution U of Eq. (14.47) are all real. Since the radical (b² − ac)^½ appearing in our definition (14.15) of B is pure imaginary in the real domain, it even turns out that all the elements of the matrix jB except those in the last row are real there. If we had elected to take the average of Eqs. (14.13) and (14.14) instead of Eq. (14.13) alone as our fifth equation in the canonical system (14.16), we would actually have been led to a matrix equation of the type (14.47) that would have been entirely real in the real domain. The choice we have made, however, has the advantage that by cross differentiation we obtain the result

analogous to Eq. (14.18), in which the highest derivatives appear only as a Laplacian. The reader should verify as an exercise that in this form the canonical system is invariant under conformal mappings¹ of the tη plane.

Cauchy’s problem for the elliptic system (14.48) can be handled by the same technique that we used for studying the special case (14.40), and thus the general quasi-linear elliptic equation (14.1) from which Eq. (14.48) arose can also be treated. For this purpose, we have only to exploit the analytic continuation that has led us to Eq. (14.48). First, however, we perform a preliminary conformal transformation of the tη plane, mapping the initial curve onto the η axis. In this connection, it should be emphasized again that in our change of coordinates the real tη plane corresponds to the real xy plane. Now our method is to allow s = ξ + jη to have arbitrary complex values once more and to pick t and ξ as the independent variables, while η is held fixed as a parameter. The system (14.47) then reverts to its hyperbolic form

for which the Cauchy problem can be solved by the classical procedures that we described earlier. The Cauchy data for Eq. (14.49) at t = 0 now appear, to be sure, as the analytic extensions into the full complex s plane of the Cauchy data for the elliptic system (14.47), which were originally defined only as real analytic functions of the real argument η. Thus Cauchy’s problem for the general quasi-linear elliptic partial differential equation (14.1) is solved in three steps consisting of a reduction to the canonical form (14.47), an analytic continuation of the data, and, finally, a solution of the hyperbolic problem (14.49) by, for example, the method of finite differences.

We should point out at this stage that further conformal transformations of the tη plane are still at our disposal, provided that they preserve the relevant segment of the η axis where the Cauchy data are prescribed. Locally, such a conformal transformation can be obtained by analytic continuation to complex arguments t + jη of an arbitrary purely imaginary analytic function of the purely imaginary variable jη. Thus, at the beginning, any analytic change of scale along the η axis can be introduced. This can be arranged, in fact, so that the initial values of any one of the unknowns x, y, p, q, and u, say of x, are given by a quite arbitrary analytic expression in η. Consider what becomes of this freedom in our choice of x as a function of η after we have continued the Cauchy data analytically into the complex s plane. Our remark then implies that for a fixed value of η we can pick the initial values of x to be any complex-valued analytic expression in the real variable ξ that assumes at −ξ the conjugate of its value at ξ. But these data are to be substituted into a stable Cauchy problem for the hyperbolic system (14.49) in the tξ plane, and therefore any convergent limiting process performed on them still leaves us with valid results. In particular, we can approximate any continuously differentiable function x = x(ξ) arbitrarily well by analytic expressions in ξ, and hence such a function is admissible as our initial choice of x in terms of ξ, even when it is not analytic. The remaining unknowns y, p, q, and u, however, must still be related initially to x by the analytic conditions of our original Cauchy problem in the xy plane.

That our initial choice of x = x(ξ) need not be analytic is a somewhat surprising feature of our procedure for solving Cauchy’s problem for an elliptic equation by analytic extension into the complex domain; the underlying explanation of this phenomenon, however, is that the function x = x(ξ) essentially determines what might be thought of as a path of integration in the real xy plane. Indeed, if we solve the hyperbolic Cauchy problem (14.49) in the tξ plane for a fixed value of η, the real portion of the answer is found along the line ξ = 0, which corresponds to a curve in the real xy plane along which the values of p, q, and u are determined. This curve would feature as a path of integration if we were to solve the hyperbolic problem (14.49) by the classical method of successive approximations. Furthermore, with each initial value of x is associated a specific initial value of y, and through the corresponding point (x,y) in the complex domain we can pass a complex characteristic of the partial differential equation (14.1) that intersects the real xy plane in a unique point lying on this same path. Therefore, for each fixed η our choice of the initial function x = x(ξ) serves to determine in a specific manner the path in the real domain along which the results of the solution of the Cauchy problem will appear. The path need not be analytic because x(ξ) need not be analytic; it could, for example, have several corners. It is of interest that, in the actual execution of a finite-difference scheme, the characteristic coordinates do not occur explicitly, so that the freedom in the selection of x(ξ) allows in practice for a situation in which we can assign an arbitrary initial column of complex values of x and subsequently can march forward with a valid solution by fixed elementary rules.

We have used the method of characteristics and of analytic continuation into the complex domain to define a stable finite-difference scheme for the solution of Cauchy’s problem in the elliptic case. There are many other significant applications of this theory. For example, Lewy⁸ has used it to establish that a sufficiently differentiable solution of an analytic elliptic partial differential equation in two independent variables is necessarily itself analytic. It is also possible to show by this procedure that a free streamline in the steady flow of an incompressible, inviscid fluid must always be an analytic curve.⁶ Such results concerning the analyticity of solutions of elliptic boundary-value problems are of interest for our discussion of inverse methods because they show that we start in the right class of functions when we assume analyticity of our data. We shall not devote space to a more detailed presentation of such theoretical matters here, however, since our principal concern is with the applications.

14.5    Flow around a Bubble Rising under the Influence of Gravity

We shall consider the plane, steady, irrotational motion of an incompressible, inviscid liquid that falls over a free streamline bounding a bubble of gas held at constant pressure. The flow is described (see pages 375, 384 of Ref. 1) by a velocity potential ϕ and a stream function ψ that can be combined to form a complex potential ζ = ϕ + jψ that is analytic as a function of the complex variable z = x + jy in the physical plane. The velocity of the motion is the gradient of ϕ, and the pressure P is given by Bernoulli’s equation

in which we have taken the density ρ to be identically 1. We assume that the fluid has significant weight, and g denotes the acceleration of gravity.

Along a free streamline bounding a bubble of the type described above, two boundary conditions are required, stating that ψ and P are constant. There is no loss of generality if we assume that the constant value of ψ along the free streamline is

and by translating the coordinate system up or down appropriately we can write the constant-pressure boundary condition in the form

according to Bernoulli’s equation (14.50). We shall study the inverse problem in which the free streamline is prescribed to be the principal arc of the cycloid

corresponding to the interval −π ≤ θ ≤ π. This is a Cauchy problem for the Laplace equation

with Cauchy data defined by Eqs. (14.51) and (14.52) along the initial curve (14.53). We shall solve it in closed form by using the fundamental procedures we have already outlined.

A first remark is that Eq. (14.54) is already in canonical form. We can make a conformal transformation that brings the initial curve into a segment of the η axis in the tη plane by interpreting Eq. (14.53) as a conformal mapping of the z plane onto the complex θ plane and then setting

Of course, ψ is still harmonic as a function of t and η,

From Eqs. (14.51) and (14.52), we easily find for ψ the initial conditions

at t = 0. Now we can extend ψ analytically to complex values η − jξ of its second argument and consider it as a function of t and ξ with η held fixed as a parameter. Under this complex substitution of variables, the Cauchy problem (14.56) and (14.57) assumes the hyperbolic form

with

for t = 0.

We can use the general solution

ψ = α(ξ + t) + β(ξ − t)

of the wave equation (14.58) in terms of two arbitrary functions α and β in order to solve this problem explicitly. We find from Eqs. (14.59) that

α(ξ) + β(ξ) = 0

and

α′(ξ) − β′(ξ) = −2g^½ sin (η − jξ)

whence

ψ = jg^½ cos (η − jξ − jt) − jg^½ cos (η − jξ + jt)

Specializing to the case ξ = 0, we obtain the real solution

of the original elliptic Cauchy problem (14.56) and (14.57).

In view of Eq. (14.55), it is an easy matter to derive, from Eq. (14.60), the relationship

ζ = −2g^½ cos θ

for the complex potential ζ = ϕ + jψ; therefore, according to our interpretation of Eq. (14.53) as a conformal transformation of the complex θ plane onto the complex z plane, we have

Although it is not possible to invert this formula in order to express ζ explicitly in terms of z, it is easy to read off the essential properties of the resulting flow directly. The liquid moves downward with symmetry about the y axis in the upper half plane. The central streamline divides at the origin into two forks along the principal arc of the cycloid (14.53), which forms the free boundary of a gas bubble underneath. At the ends of this arc of the cycloid, the two forks of the central streamline change over into vertical segments that descend indefinitely. Far to the right and to the left the velocity approaches a constant downward value. Thus the flow represents the motion of a liquid falling around a semi-infinite channel of gas bounded at the top by a smoothly forked free streamline in the shape of a cycloid. See Fig. 14.3. If we think of our coordinate system in the xy plane as moving upward uniformly at a speed just sufficient to cancel out the velocity of the liquid at infinity, the plane flow we have described can be interpreted as the motion of a gas bubble rising in a liquid of infinite extent. The part of the gas under the bubble, which is bounded by vertical walls, can be understood in this connection to represent a wake.

Fig. 14.3 Flow over a gas bubble.

We have gone out of our way to solve the elementary Cauchy problem defined by Eqs. (14.51), (14.52), and (14.54) by using general procedures applicable to a wider variety of cases. This particular example, however, could be treated more directly by means of analytic continuation alone. It suffices for this purpose to derive a relationship between ζ and z along the cycloid (14.53) and then continue it analytically over the whole z plane. As an exercise, the reader should carry out the details of such a solution based exclusively on the theory of analytic functions of a complex variable.

14.6    The Detached-shock Problem

We are finally in a position to develop with reasonable facility an inverse method for calculating steady plane or axially symmetric flows of a compressible inviscid fluid behind a prescribed shock wave.⁵ In our presentation, we shall confine our attention to the axially symmetric case, which we can study in a meridian xy plane. Here y plays the role of the radius in a cylindrical coordinate system. We denote by u and v the components of velocity parallel, respectively, to the x axis and the y axis. Although the flow we shall consider is rotational, the law of conservation of mass permits us to describe it by means of a stream function ψ = ψ(x,y) such that

ψ_x = −yρvψ_y = yρu

where ρ is the density. The fundamental laws of conservation of momentum and energy in fluid mechanics can be used to show that ψ satisfies the second-order partial differential equation

of the quasi-linear type (14.1). Here A is an arbitrary function of ψ, C is the local speed of sound defined by

C² = γAρ^γ−1

where γ is the ratio of specific heats, and

is a constant called the stagnation enthalpy.

We are interested in the flow behind a shock wave that occurs along a prescribed curve

We assume the flow in front of the shock wave to be uniform and in the direction of the x axis, with a Mach number M > 1. Across the shock wave, of course, the various flow quantities have jumps that are determined by the laws of conservation of mass, momentum, and energy, and that can be expressed at each point of the wave in terms of the slope

there alone. In particular, conservation of energy implies that the stagnation enthalpy H is continuous across the shock, and this fact explains why it is to be constant throughout the flow. Similarly, the law of conservation of mass shows that the stream function ψ has no jump. Conservation of momentum in the directions of the x axis and the y axis can be used to calculate the arbitrary function A explicitly in terms of the given function F and to determine the normal derivative of ψ at each point of the shock wave. If we normalize so that the pressure P and density ρ in front of the shock wave are identically 1, and if we call the speed there Q, then the formulas derived from the above shock conditions for flow quantities evaluated just behind the shock are

and

where

Somewhat more complicated expressions would feature in these relationships if we were to consider a more involved, but occasionally useful, formulation in which the constant γ merely plays the role of an effective ratio of specific heats behind the shock and differs from the actual ratio of specific heats γ₁ in front.

The conditions (14.63) and (14.64) along the prescribed curve (14.62) present us with a Cauchy problem for the partial differential equation (14.61). The coefficients of this equation depend in an essential way on the quantities y, ψ_x, ψ_y, and ψ, and therefore we must use our full theory to obtain a solution. Furthermore, the discriminant of Eq. (14.61) has the form

b² − ac = C²(u² + v² − C²)

so that this partial differential equation is of the hyperbolic type in a region of supersonic flow with u² + v² > C² and is of the elliptic type in a region of subsonic flow with u² + v² < C². As we have seen, Cauchy’s problem is straightforward for a hyperbolic equation, but can be handled in the elliptic case only after an extension into the complex domain. We can treat both situations by means of our reduction to canonical form, but the details of the analysis are really quite different for the two types of flow. In particular, we must assume that any arc of the shock wave (14.62) bounding a subsonic region is analytic.

We have given a complete description of the method of characteristics and of finite differences for the solution of a Cauchy problem such as the one with which we are confronted in the present situation. In this special application of our procedure, all the steps to be taken are quite routine except for the analytic continuation, into the complex domain, of the Cauchy data along a subsonic arc of the shock wave (14.62). Even here, the only difficulty is with the extension of the function F(y²) to complex values of the variable y, since the remaining expressions that appear in the data (14.63) and (14.64) or in the canonical system (14.49) are elementary in the sense that they involve only additions, multiplications, divisions, or possibly square roots, and hence their extension is merely a matter of direct substitution of complex arguments. Thus, in case F(y²) is also defined by an explicit formula that can be continued analytically over the y plane by straight-forward substitution of complex values of y, the Cauchy problem (14.61), (14.63), and (14.64) is completely solved in a satisfactory fashion. Since a wide variety of shock waves (14.62) of an elementary form like this can be found, we can obtain a large class of flows with detached shock waves by applying our inverse method based on prescribing F. The body that generates the shock wave (14.62) appears, of course, as the level curve ψ = 0 of the stream function ψ. It is not too hard in practice to determine a choice of the shock wave that corresponds to a blunt body within reasonable approximation to a desired shape. This can be achieved with moderate success by studying the behavior of the analytic function F(y²) in the complex y plane.

Fig. 14.4 Detached shock in front of a blunt body in hypersonic flow.

In Fig. 14.4 we have sketched an axially symmetric flow of the type under consideration here, for which the shock wave was selected to be a hyperbola. In this calculation, the Mach number in front of the shock wave was taken to be M = 20, and the ratio of specific heats there was γ₁ = 1.4. Behind the shock, however, an effective value γ = 1.17 of the ratio of specific heats was assumed in order to approximate in some degree the effects of dissociation. It is this small choice of γ that brings the shock wave so close to the body, which turned out to be very nearly spherical.

Certain qualitative properties of flows with a detached shock wave can be deduced from our analysis based on analytic continuation. We observe, for example, that, in order for the flow constructed by our inverse method to exist in the large and to reach as far as the body without exhibiting a singularity, it is necessary that the portion of the Cauchy data prescribed along the subsonic section of the shock wave should have analytic continuations that are regular in the large in the complex domain. It follows that the subsonic arc of the shock wave lies in the interior of the region of regularity of the stream function ψ in the four-dimensional domain of the two complex variables x and y. For this reason, the values of ψ or of its derivatives on this arc should not change very much when the boundary values of ψ are perturbed along the body, since in the closed interior of its region of regularity an analytic function is not greatly affected by alterations of its boundary values. In other words, the shape of the shock wave in front of a subsonic portion of flow should not be very sensitive to changes in the shape of the body that generates it. Those readers familiar with experimental work on detached shock waves will recognize that the phenomenon that we can predict in this manner does actually occur in nature.

A final comment about our construction of flows through prescribed shock waves is that it applies to the case of attached shocks, too. The difficulty here is to formulate the correct asymptotic form for the equation of the shock wave near its point of attachment. This is easy if the flow is entirely supersonic, but if it is subsonic in the neighborhood of the point of attachment, we are led to a problem concerning the behavior of the solution of an elliptic partial differential equation near the confluence of two analytic boundary conditions. The method of analytic extension into the complex domain indicates that the appropriate asymptotic expansion (see Chap. 5) here should proceed in fractional powers of z and of z log z when we normalize so that the point of attachment lies at the origin. However, we must close with the admission that this interesting question is far too involved for further discussion in the present chapter.

REFERENCES

1.   Beckenbach, Edwin F., Conformal Mapping Methods, chap. 15 in “Modern Mathematics for the Engineer,” First Series, edited by E. F. Beckenbach, McGraw-Hill Book Company, Inc., New York, 1956.

2.   Courant, Richard, Hyperbolic Partial Differential Equations and Applications, chap. 5 in “Modern Mathematics for the Engineer,” First Series, edited by E. F. Beckenbach, McGraw-Hill Book Company, Inc., New York, 1956.

3.   —— and K. O. Friedrichs, “Supersonic Flow and Shock Waves,” Interscience Publishers, Inc., New York, 1948.

4.   —— and D. Hilbert, “Methoden der mathematischen Physik,” vol. 2, Springer-Verlag, Berlin, Vienna, 1937.

5.   Garabedian, P. R., and H. M. Lieberstein, On the Numerical Calculation of Detached Bow Shock Waves in Hypersonic Flow, J. Aero. Sci., vol. 25, pp. 109–118, 1958.

6.   ——, H. Lewy, and M. Schiffer, Axially Symmetric Cavitational Flow, Ann. Math., vol. 56, pp. 560–602, 1952.

7.   Lehmer, Derrick H., High-speed Computing Devices and Their Applications, chap. 19 in “Modern Mathematics for the Engineer,” First Series, edited by E. F. Beckenbach, McGraw-Hill Book Company, Inc., New York, 1956.

8.   Lewy, H., Neuer Beweis des analytischen Charakters der Lösungen elliptischer Differentialgleichungen, Math. Ann., vol. 101, pp. 609–619, 1929.

9.   Riemann, B., “Gesammelte mathematische Werke,” Leipzig, 1876.

10.   Schiffer, Menahem M., Boundary-value Problems in Elliptic Partial Differential Equations, chap. 6 in “Modern Mathematics for the Engineer,” First Series, edited by E. F. Beckenbach, McGraw-Hill Book Company, Inc., New York, 1956.