Ockendon J., Howison S., Lacey A., Movchan A. Applied Partial Differential Equations

Подождите немного. Документ загружается.

398

MISCELLANEOUS TOPICS

If, instead, we specify the initial data u = uo(x') for xl = 0, where x' _

(x2i ... , xm)', and f = 0, our integration provides the standard formula

uW =

J _

GA1uodx'. (9.15)

What is not standard is the fact that, for m-dimensional problems, even for

x1 < t1, G is a distribution with `mass' concentrated on the (m - 1)-dimensional

characteristic surfaces through x = t. Hence the integrals in (9.14) and (9.15) in-

clude contributions from these surfaces and their intersection with x1 = 0, as well

as from the region inside. This is a crude generalisation of the discussion about

propagation of discontinuities for the Riemann functions that we gave at the end

of §4.2.2.

An exactly similar calculation can be carried out for the Dirichlet problem for

an elliptic system of the form (9.11) with, say,

u=(tf.l,...ru2n)T,

and

(ul,...,un)T =0

on the boundary 3D of a closed region D. We again define G to satisfy (9.13), but

now its last n columns vanish on 3D. The asymptotic behaviour of G near x = f

is now less easy to discern than it is for the ode case, as can be illustrated with

the following model.

9.2.2 Linear elasticity

In §4.7.1, we have already encountered the famous Navier-Lame equations of linear

elasticity. In the case of static equilibrium in three dimensions, they are

Cu = pV2u + (A + µ)VV u = -f(x),

(9.16)

where u is the displacement, f is the body force per unit volume and A and is are

the so-called Lame constants which characterise the material.

In view of the identities

and

it is easy to see that, when suitable boundary conditions are prescribed for u, the

operator G is self-adjoint.

There are two remarks to be made before we start. First, the system is second

order rather than first; although it could, with labour, be recast in the form (9.11),

there is no reason why any of the discussion above should be confined only to

linear systems of first-order pdes, and we will thus proceed directly with (9.16).

Secondly, (9.16) has constant coefficients and hence could, by cross-differentiation,

be reduced to a scalar equation for, say, any of the components of u. This is also

an unnecessarily tiresome procedure because, although the outcome is relatively

LINEAR SYSTEMS REVISITED 399

simple (see Exercise 9.4), the physical significance of the Green's matrix is soon

obscured.

For simplicity, let us just consider the solution of (9.16) in the case when

the elastic continuum extends to infinity in all directions, and let us assume

that Jul grows no faster than lxi at infinity and that the strain components

8u;/8x; + au;/ax; decay at least as fast as

O(1/IxI3).175 We

note that the rigid

body displacement u = c. where c is constant, satisfies Cu = 0 and the decay

conditions above, and hence the Alternative implies that176

for all c, so that

11(x) dx = 0.

Also Cu = 0, with the decay conditions, is satisfied by u = w A x for constant w,

which corresponds to a rigid body rotation (since we are only considering small

displacements u, a rotation is

rby a vector product). Hence

that

f exerts zero net force and moment on

the elastic continuum.

Now, again using the fact that C is self-adjoint, we are guided by (9.13) to

simply define the Green's matrix G to be the solution of

CG(x - 4) = -6(x - 4)I. (9.17)

with

G(x - 4) -> 0

as lxI -> oo. (9.18)

each side being interpreted entry by entry. It is easy to see that G is symmetric

and depends only on x - 4; our usual procedure gives at once that

f ((CG) u - GCu) dx = -

J(o(x

- 4)u(x) - G(x - 4)f(x)) dx, (9.19)

and hence, since the left-hand side vanishes, that

u(4) =

G(x - 4)f(x) dx.

It is interesting to note that, were we to carry out this calculation ab initio without

relying on the validity of (9.19) for distributions (cf. the comments after (4.13)),

175These conditions ensure that the elastic strain energy is finite.

176For the rest of this section we denote integrals over the whole of three-dimensional space by

f . dx.

400

MISCELLANEOUS TOPICS

then the application of Green's theorem yields a contribution from a small sphere

centred at t that can be interpreted in terms of the forces applied to that sphere.

Thus G can be interpreted physically in terms of the response to arbitrarily di-

rected point forces applied to the medium at x = 4 (and hence we do not expect

its entries to decay as rapidly at infinity as we insisted for u above).

In order to find G explicitly, we simply look at the Fourier transform of (9.17).

Defining

G(k, t) = r

G(x - t)e'k'x dx,

where the transform variable k is now a vector, it is easy to see that177

(9.20)

ik'f

kkT I

I -

G(k

) ( .

(9 21).

j 2

p(A + 2p) IkI4

Some care is needed with the Fourier inversion of this singular function (see Ex-

ercise 9.5), but it transpires that the entries in G are

G+,(x -

4) =

A + p

A +

3pai} +

(xi - G)(x,Z fit)

(9.22)

87rp(A+2p)

Ix - fl

The presence of off-diagonal elements would have been difficult to spot directly

from (9.17), and they illustrate the statement made above, that each component

of f influences all the components of u. The columns of G are precisely the above-

mentioned displacement vectors associated with the `point forces' along the coor-

dinate axes.

When boundaries are present, (9.22) still describes the local behaviour of G

near x = , but, not surprisingly, the application of the concept of images is now

much more complicated than in §5.6.

9.2.3 Linear inviscid hydrodynamics

We now consider what can happen to singular solutions of systems that have a

real characteristic but are not hyperbolic, thereby again illustrating the surprising

way in which the components of f in (9.11) can influence the components of u.

We will study a very simple model of inviscid hydrodynamics which shows how

the fundamental concepts of lift and drag on a body moving through a stationary

liquid can be understood in terms of distributional solutions of a system of linear

pdes. The model is a truncated version of the incompressible Euler equations, to

which (2.3)-(2.7) are related, and to which we will refer again in §9.4.4. When the

fluid is being driven by a localised force that is so weak that only linear terms need

be retained, these equations take the dimensionless form17e

S +Vp=f,

at t = 0, and we now study the response to various forms of f which

might describe different kinds of small bodies moving through the fluid.

177 By kkT we mean the matrix (kikt).

1781n this chapter only, we revert to the conventional use of u as the fluid velocity in hydrody-

namics.

LINEAR SYSTEMS REVISITED 401

In all cases we assume the fluid extends to infinity in all directions and is at

rest there. Also, we expect that much of the flow is irrotational because, in regimes

where f = 0, the vorticity V A u is independent of time. Hence, if the vorticity is

zero initially, then we expect it to stay zero, but we will encounter some surprises

in this respect.

A lift force in two dimensions

Consider a small lifting body, such as an aerofoil moving along the x axis. We

assume that the lift force Lj, in the y direction, is prescribed, along with the

aerofoil position X(t)i. When we take axes moving with the aerofoil, we expect u

and p to satisfy

8 - X(t) 8 +Vp =

V u = 0,

(9.23)

where l; = x - X (t). When L and X are constant., so that we can set Ou/Ot = 0,

we find that the Fourier transform

u(k, y) =

o0 oc

ueiktf+ik2Y d{ dy

Z. f

is given by

with

u =

(k2i - k1j),

(9.24)

ik2L

kl +

It is a simple matter to invert these transforms using the fact that the transform

of logr is -27r/(kl + ka) because V2logr = 21r&(x)6(y); the answer is

L 1 Ly

u =

2rX

+ y 2 "

P = 2x02 + y2) ,

(9.25)

which is the velocity field of a vortex centred at x = X (t), y = 0. Hence a small

moving aerofoil can be identified with such a vortex.

A quite different result sometimes emerges when we allow the velocity and the

lift on the aerofoil to vary. As long as L = I'X, where I' is constant, we find

d ik, .r

I'Xeik, X

(ue

) = k4(-k2i+k1j),

(9.26)

and hence we retrieve (9.25). However, as soon as L s PX, (9.26) cannot be

integrated directly and u acquires a `history'. For example, if L = rxH(t), where

H is the Heaviside function, (9.25) is replaced by

I yi-xj\\yi-&j_

27r

f 2

+ y2 x2 + y2 1 ,

(9.27)

402 MISCELLANEOUS TOPICS

which shows that a `starting vortex' is shed at the origin x = y = 0.179 Our

expectation that the flow would be irrotational away from the aerofoil is unjustified

in this case; the single real characteristic of (9.23) has propagated information

about the motion of the aerofoil.

Drag in two dimensions

This is an interesting problem to pose because it is well known that d'Alembert's

paradox [27J precludes the existence of a drag force on a body moving with constant

velocity in an irrotational flow. However, if we attempt to model a drag force D

by naively replacing (9.23) by

X a + Vp

and X are constant and 8u/8t = 0,

A(k2 k2)

i-ks.l),

kDk ,

(9.29)

and inversion leads to the appearance of a delta function in the components of u

as we approach the drag-producing body at f = y = 0. In fact, it is well known

that drag on an accelerating body, usually referred to as added mass, must be

modelled by

o - X8u + Vp = 0,

V- u =M61 wgy).

(9.30)

The last term is a `mass dipole' and it leads to the velocity field

u = 2

(8F82

a(logr)i +

80Y

(logr)j)

r2 = f2 +y2,

(9.31)

which differs from the inverse transform of (9.29) by the aforementioned delta

function in u. Hence (9.28) is not a good representation of a drag force. The

correct representation (9.30) models the flow generated by a small non-lifting body

for arbitrary k and the exact solution of the Euler equations reveals that M is

proportional to X.

Three-dimensional flow

The situation becomes even more interesting in three dimensions. When a small

aerofoil at (Vt, 0, 0) is modelled by

-V

+ Vp = L6(06(y)b(x)j,

V u = 0,

we find the three-dimensional Fourier transform

179When viscosity is taken into account, this vortex diffuses through the whole flow field, which

is eventually steady and such that f u ds around the aerofoil is equal to r, which is called the

circulation.

LINEAR SYSTEMS REVISITED 403

Fig. 9.1 Horseshoe vortex.

iLk2

Vkl(k2 + k22 +

k2)k.

Hence u = V O, where

L y

x - Vt

-47rVy2+z2

((x-Vt)2+y2+z2)1/2

so now, even in steady flow, a wake is shed along the degenerate characteristic

y = z = 0, x < Vt. This is the famous `horseshoe vortex' [26] (see Fig. 9.1).

Our approach can even shed light on the well-known phenomenon of vortex

ring propagation, as produced, say, by a suitably exhaled puff of cigarette smoke.

Using the technique described above, we can show that the solution of a steady

axisymmetric flow modelled by

-V 8 + Vp =

(z)i, V V. u = 0

D(g2 _

- zz, 3t b, 2z)

for (t, y, z) 0 0,

41r

y2 + z2)

which is the far-field of a vortex ring (see [39] and (7.114)). This velocity field

differs by a delta function at the origin from that produced by a mass dipole as in

(9.30), but, more importantly, it shows that we can regard the motion produced

by a smoke ring as equivalent to that generated by a suitable point force in the

equations of motion.

9.2.4

Wave propagation and radiation conditions

For linear systems of pdes which model wave propagation, it is possible to proceed

directly in the time domain by analysing the equation with three independent

space variables and time t, say by using Riemann matrix ideas. However, we have

often remarked that, if the coefficients in the system are independent of time,

404

MISCELLANEOUS TOPICS

then it is easier and frequently of more practical interest to restrict attention to

monochromatic waves in the frequency domain by writing the dependent variables

as functions proportional to a-i' t, say 180 We did that for the scalar wave equa-

tion in §§5.1.5 and 8.1, and some of the ideas in those sections carry over to the

vector case. However, many fascinating problems arise, such as the question of

propagation through periodic media. As mentioned in §4.5.4, in one dimension

this leads to the notions of pass and stop bands bounded by eigenvalues of a pe-

riodic Sturm-Liouville problem, but this framework relies on the well-developed

Floquet theory for odes with periodic coefficients. It is much harder to understand

multidimensional configurations of this kind because of the geometric complexity

of the waves as they reflect from each periodic cell boundary. It is a great pity

that we cannot discuss this further here because of the fundamental importance

of wave propagation in crystal lattices.

Another issue that arises when we consider problems in the frequency domain

for unbounded wave-bearing media that are uniform at large distances is the ques-

tion of the radiation conditions that generalise the Sommerfeld condition (5.75)

which was so vital in ensuring uniqueness for scalar problems. For example, for

Maxwell's equations (see §4.7.2) we may suppose that the leading-order far-field

radiation due to any bounded source of non-zero intensity decays with distance

in proportion to e'kr/r (or eikr/f in two dimensions), where k = w/c is the

wavenumber. Now all electromagnetic waves are transverse waves in the sense that

the directions of the fields E and H are perpendicular to the direction of propaga-

tion of any wave; this is a trivial consequence of (4.82) in the frequency domain,

because, when we seek such solutions in which E = E(k r) and H = H(k r),

both E and H are clearly perpendicular to r. Hence we may write that, at large

distances, where k and r are nearly parallel,

eikr eikr

EH--h,

r r

where e and h are both azimuthal and are independent of r to lowest order.

Then (4.82), suitably scaled, gives at once the radiation conditions

O= VAE-iwH.i(kAe-wh) =ol 1 1

/\l

O= V AH+iwE.i(kAh+we) = o 11)

asr-aoo.

Similarly, in elasticity, the wave equations (4.73) imply that, far from any

bounded region of sources, u decays like uoeikr/r, where again k is radial. As

expected, we find that either

k = k

where

uo k, = 0

and w2 =

k°µ

'80Such waves can always be superimposed to generate solutions in the time domain, but this

may be easier said than done.

COMPLEX CHARACTERISTICS AND CLASSIFICATION BY TYPE 405

corresponding to transverse shear waves (S-waves), or

k = kp,

where

uo A ky = 0 and

k»(,\+ 2µ)

corresponding to longitudinal compressive waves (P-waves). Hence the radiation

condition is that, at large distances,

elks eik,r

u-u8-+up-,

r r

where u8 k8 =0 and uy A kp =0.

9.3

Complex characteristics and classification by type

We have highlighted in many places the dramatic differences between ellipticity,

parabolicity and hyperbolicity. One attempt that has been made to reconcile these

two concepts is Garabedian's theory of complex characteristics [21). The basic idea

relies on the fact that the solution of an elliptic partial differential equation is

usually an analytic function of the independent variables. Hence, for the elliptic

equation

020 820

0{2 + 9y2

°,

we could seek solutions analytic in

_ + irl and note that this implies that

020

82¢

3 2

X12

= 0.

Subtracting the two equations, we see that

820

020

aye

2 =

which is a hyperbolic equation in (y, 7)). However, this procedure inevitably involves

analytic continuation of the boundary data which, as we have already remarked,

is a dangerous procedure, and it is only justified as long as no singularities are

encountered in the continuation into rl $ 0 (see Exercise 9.7).

Armed with our knowledge of ray theory from Chapter 8, we can gain further

insight into the pitfalls encountered in trying to unify ellipticity, parabolicity and

hyperbolicity. Consider, for example, the following problems in y > 0:

ate- +8-=0,

(9.32)

with

0(x, 0) =

e-: /2`,

e > 0,

(9.33)

and either

(i)

a>0

and asy-oo

406

MISCELLANEOUS TOPICS

(ii) a<0 and

ony=0.

It is easy to see that the decay condition is sufficient to ensure uniqueness in the

elliptic case (i); the Cauchy data, which is equivalent to saying that 0 is a function

of x - Vr--ay, ensures uniqueness in case (ii). (When a = 0, the `parabolic' case,

the only solution that is bounded at infinity is O(x, y) = e z2/2i.)

Although the parameter a can be scaled out of either problem, we might expect

that sensible limits could be retrieved as a -> 0. Indeed, it is possible that in this

limit 0 would be closely related to the Green's function or the Riemann function,

respectively, and we will see the precise relation later. Now a ray theory ansatz in

which Ae°/[ yields the eikonal equation

a MY

MU), =0,

and Charpit's method then gives u = -s2/2, where, in case (i),

x=s(1-t), y=fib, u

-_(xf

and so there is no possibility of exponential decay as y -+ oo. In case (ii), however,

we have

x=s(1-t), y= -

Thus in case (ii) the characteristics are x = s + y//&, and we have retrieved

the hoped-for propagation of the data along the characteristic x = ys, but

case (i) has woefully failed to reveal anything like the singularity associated with

the Green's function for (9.32). The ray ansatz has worked well when the equation

is hyperbolic but leads to solutions with unbounded oscillations and growth in the

elliptic case.

Now we can solve either case by Fourier transforms. The answer is

dk, (9.34)

(i) O(x, Y) _

F00

dk; (9.35)

- 2

irv

in each case the coefficient of y in the exponent is determined uniquely by the

data.

We now have to rely on some results from the theory of asymptotic expansions

of integrals [22]. As a - 0, (9.35) is always approximated by the saddle point

contribution that emerges when we write k = rc/e to obtain

1 e-(K2/2+1 Udsc

2e _.

QUASILINEAR SYSTEMS WITH ONE REAL CHARACTERISTIC 407

Fig. 9.2 Inversion contours for case (i) and case (ii); we have taken x = 0 in both cases.

The steepest descents direction at the saddle is shown dashed.

The saddle point is at r. = i(x - yvr--a), leading to a term e-2in 4.

However, the integral for case (i) must be written as

dK,

a 2e

which is always approximated by its end-point contribution from near r. = 0,

namely

eyV/Q-/(ay2 +x2), which is proportional to the y derivative of the Green's

function for (9.32) with Dirichlet data on y = 0.

Thus we see that the switch from elliptic to hyperbolic can be identified with the

presence or absence of a saddle point on the inversion contour for a Fourier integral.

In case (i), the inversion path cannot be deformed into a 'steepest descents' contour

(see Fig. 9.2(i)), but in case (ii) it can (Fig. 9.2(ii)). Note that the function that

emerges from the ray ansatz has no Fourier transform in the usual sense and is

hence legislated against when we write down (9.34).

9.4

Quasilinear systems with one real characteristic

Annoyingly, practical problems often give rise to systems of pdes which, like that

in §9.2.3, are neither elliptic nor hyperbolic in the sense of Chapter 3 and are

often nonlinear as well. Hence they usually need to be treated individually on

their merits and here we will simply list three examples.

9.4.1 Heat conduction with ohmic heating

In models of many electric devices, such as thermistors, the equations of conductive

heat transfer must be coupled with those of electromagnetism via the resistance

or 'ohmic' heating of the device. In one of the simplest configurations, the current

density j is related to the electric potential 0 by Ohm's law

j = -a(T)VO, (9.36)

where the electric conductivity a is a positive function of the temperature T. In

quasi-steady conditions, conservation of charge gives