## Abstract

Internal gravity waves are generated in a stratified fluid by arbitrary forced oscillations of a horizontal disc. The wave fields are calculated in both the time domain and the frequency domain. In the time domain, an initial-value problem is solved using Laplace transforms; causality is imposed. In the frequency domain (time-harmonic oscillations), a radiation condition is imposed: a plane-wave (Fourier) decomposition is used in which waves with outgoing group velocity are selected. It is shown that both approaches lead to the same solution, once transient effects are ignored. Then, a method is given for calculating the far-field, using asymptotic approximations of double integrals. It is shown that the total energy flux is outwards, for arbitrary forcings of the disc. Further investigations of energy transport are made with a view to clarifying the nature of radiation conditions in the frequency domain.

## 1. Introduction

It is well known that internal gravity waves exist in a density stratified fluid, and that they can be generated by oscillating immersed bodies or by scattering from topography. The motivating problem here is to calculate these wave fields, in a three-dimensional unbounded fluid domain. It is also known that the energy produced by the oscillations is confined to characteristic conical wave beams. See, for example Lighthill (1978, chapter 4) or LeBlond & Mysak (1978, §8).

Under the Boussinesq approximation, with an inviscid fluid, a uniform density stratification and small-amplitude motions, we can write down the governing equations (see §2 for details). In particular, for time-harmonic motions, the pressure is Re{*p*e^{−iωt}}, where *p*(*x*,*y*,*z*) solves
1.1*z* is the vertical coordinate and *N* is the constant Brunt–Väisälä frequency. We are interested in the case where 0<*ω*<*N*, so that equation (1.1) is hyperbolic. This equation has to be solved subject to a boundary condition (prescribed normal velocity) and a far-field (radiation) condition.

There are several ways to specify the radiation condition in the frequency domain; for a review, see Voisin (1991, §3.2). One way is to replace *ω* by *ω*+*iε*, letting *ε*→0 at the end of the calculation (Lighthill 1978, §4.9). Another way uses analytical continuation in the complex *ω*-plane; see below for more details. When plane-wave decompositions are available, we can examine the corresponding group velocity; note that, for plane internal waves, group and phase velocities are orthogonal (Lighthill 1978, §4.4). Barcilon & Bleistein (1969) used a local plane-wave approximation at each point in the fluid but their approach was criticized by Baines (1971).

Underlying radiation conditions in the frequency domain is the requirement of causality in the time domain. Therefore, because of the variety of frequency-domain radiation conditions, we start by treating an initial-value problem. The problem we have chosen concerns an oscillating horizontal disc: there is no motion for *t*<0, and then time-harmonic forcing is switched on at *t*=0. We solve the resulting problem using Laplace transforms in time and Fourier transforms in the horizontal coordinates *x* and *y*. This leads to solvable dual-integral equations. Inverting the Laplace transform, imposing causality, gives a transient contribution plus a time-harmonic forced response; we are interested in the latter.

Next, we formulate the analogous time-harmonic problem. A radiation condition is imposed in the Fourier domain: we ensure that only those plane-wave contributions with group velocity directed away from the disc are retained. This approach has been used by others: it is limited to those problems where plane-wave decompositions are available throughout the fluid domain. We find that the pressure field obtained is exactly the same as the forced response obtained in the initial-value problem. (This is an example of the so-called Limiting Amplitude Principle.)

Having confirmed that our time-harmonic solution does respect causality, we go on to estimate far-field quantities (§7). This requires the asymptotic estimation of certain double integrals; a form of the method of stationary phase is used, in which the stationary phase points occur along a line in the domain of integration. Estimates for both the pressure, *p*, and the velocity, ** v**, within the wave beams are obtained. Then, the time-averaged energy transport vector,

**I**, defined by 1.2can be calculated, where the overbar denotes complex conjugation and

*ρ*

_{0}is the constant background density. The equations of motion imply that div

**I**=0. Consequently, if

*S*is a closed surface enclosing the oscillating plate, then we should have , where

**is the unit normal on**

*n**S*pointing away from the plate. We verify this (global) inequality in the far field, for arbitrary forcings of the plate (§7

*c*).

It is clear that **I** is arbitrary in the sense that it could be replaced by **I**+**I**_{1} provided that div **I**_{1}=0 (Longuet-Higgins 1964). This means that there is no reliable physical interpretation of **I****⋅**** n**, although Lighthill (1978, p. 14) states that it gives the energy ‘being transported in the direction of

**across a small plane element at right angles to**

*n***’. Nevertheless, we might expect that**

*n***I**

**⋅**

**>0, pointwise. However, we show, by example, that this is not the case. This confirms that examining the sign of**

*n***I**

**⋅**

**at each point in the far-field cannot be used as a radiation condition: doing so could exclude physically meaningful solutions.**

*n*Let us conclude this introduction by mentioning previous work on the generation of three-dimensional internal gravity waves by oscillating bodies. Most of it is concerned with axisymmetric motions of spheres and spheroids (see, for example, Hendershott 1969; Sarma & Krishna 1972; Lai & Lee 1981; Appleby & Crighton 1987; Voisin 1991, 2003; Flynn *et al.* 2003; Voisin *et al.* 2010). The approach of these authors starts by solving equation (1.1) when *ω*>*N* using a scaling in *z*. This leads to a boundary value problem for Laplace's equation exterior to a spheroid, a problem that can be solved by separation of variables. The solution for *ω*<*N* is then obtained by analytical continuation. In particular, Sarma & Krishna (1972) and Lai & Lee (1981) obtained the solution for vertical oscillations of a rigid disc; we find agreement with their solutions (§6*a*).

For disc problems, there is a paper by Gabov & Pletner (1988), in which the pressure (not the normal velocity) is prescribed on the disc. There is a recent study by Davis & Llewellyn Smith (2010), where the plate oscillates in its own plane but the fluid is viscous. There are also several papers by Chashechkin and his collaborators, in which internal waves are generated in a viscous fluid by a circular piston in a horizontal rigid plane, *z*=0. Mathematically, this problem is simpler than a disc problem because the velocity is prescribed over the whole plane, *z*=0. See, for example, Chashechkin *et al.* (2004) and Bardakov *et al.* (2007); these papers include both theoretical and experimental results.

## 2. Mathematical formulation

We consider a variable density fluid without rotation. We use the linearized form of the Boussinesq equations (Vallis 2006, §2.4.2),
Here, we have Cartesian coordinates *Oxyz*, with *z* pointing upwards; is a unit vector in the *z*-direction. The velocity is and is a rescaled pressure: the actual excess pressure is , where *ρ*_{0} is the constant background density. The buoyancy frequency *N* is a positive constant and is the buoyancy. Henceforth, tildes indicate (real) quantities in the physical domain.

We remove the time-dependence by Laplace transformation. Thus, we define, for example,
Assuming zero initial data, and eliminating *b*, we obtain
The first of these gives
Then, div ** v**=0 gives
2.1This is a partial differential equation for

*p*. For time-harmonic problems (with a time-dependence of e

^{−iωt}), the same equations (for

*p*,

*u*,

*v*and

*w*) are obtained but with

*s*=−i

*ω*(see §5).

## 3. A horizontal plate in a stratified fluid

We have a thin flat plate, *Ω*, in the *xy*-plane. The rest of the *xy*-plane is denoted by *Ω*′. The motion is forced by oscillating the plate, with prescribed normal velocity on *Ω*, . Thus,
which becomes, after Laplace transformation,
In fact, this holds on both sides of *Ω*. It follows that *p* must be an odd function of *z*, so the problem can be reduced to one in the half-space *z*>0. As we have split the problem into two half-space problems, we must also impose continuity of *p* across *Ω*′.

Take the Fourier transform of equation (2.1) with respect to *x* and *y*, with
3.1the inverse Fourier transform is
3.2(Henceforth, we write when the integration limits are as in equation (3.1) or (3.2).) Then, with *κ*^{2}=*ξ*^{2}+*η*^{2}, and writing *P*′′=∂^{2}*P*/∂*z*^{2}, we obtain *P*′′=(*κ*/*λ*)^{2}*P*. Hence,
3.3We require *p*→0 as , so we take *C*=0, assuming that Re *λ*≥0.

Let us define the discontinuity in *p* across *z*=0 by
3.4Clearly, we must have *δ*=0 on *Ω*′, so the Fourier transform of *δ* is
3.5using equations (3.3) and (3.4). Also, the Fourier transform of *w* is given by

Applying the boundary conditions gives
3.6and
3.7Substituting for *W*(*ξ*,*η*,0) and *P*(*ξ*,*η*,0) gives
3.8and
3.9where
3.10Equations (3.8) and (3.9) are a pair of dual-integral equations for *D*. They are of a standard kind, so that various methods are available.

Another option is to substitute for *D* in equation (3.8) using equation (3.5). This would lead to a hypersingular boundary integral equation over *Ω* for the discontinuity function, *δ*.

To make analytical progress, we now suppose that *Ω* is circular.

## 4. Circular disc

### (a) Dual-integral equations

For a disc of radius *a*, use polar coordinates,
4.1and suppose for simplicity that *δ* has a Fourier expansion in the form
4.2where *ϵ*_{0}=1 and *ϵ*_{n}=2 for *n*≥1. Using equation (3.5) and
we find that
4.3Then,
4.4with a similar expansion for the integral in equation (3.8). Therefore, if *q* has an expansion
4.5the dual-integral equations (3.8) and (3.9) become, for *n*=0,1,2,…,
4.6and
4.7

### (b) Expansion methods

A standard method for solving integral equations (such as equations (4.6) and (4.7)) begins by expanding the unknown function in a series of basis functions. We should choose the expansion functions with two criteria in mind. First, *d*_{n}(*r*/*a*) gives a radial component of *δ*, and this quantity is known to have a square-root zero at the edge of the plate, *r*=*a*; this behaviour (which can be inferred from a local expansion near any point on the edge) should be incorporated in the expansion functions. Second, it would be useful if the integrals relating *d*_{n} to *D*_{n}, equation (4.3), could be evaluated analytically for the chosen expansion functions. Of course, we may choose to expand *D*_{n} directly (this is essentially Tranter's method (Tranter 1954, 1966)) but we prefer to begin with a quantity that has known physical properties.

Define a function by
4.8where is a Gegenbauer polynomial. As is of the form multiplied by a polynomial in *ρ*, the square-root zero at the plate edge can be incorporated automatically. It can also be shown that is proportional to both
where is an associated Legendre function and is a Jacobi polynomial.

As Gegenbauer polynomials are orthogonal, so are . Perhaps more importantly, we have 4.9where is a spherical Bessel function. This formula is due to Tranter (1954, eqn (5), 1966, eqn (8.6); see also Krenk 1979, eqn (72)).

We note in passing that the polynomials have been advocated by Boyd (2001) (§18.5.1) for use as radial basis functions in spectral methods for problems posed on a disc, 0≤*r*<1. For us, the functions , defined by equation (4.8), are preferable because of the availability of Tranter's integral, equation (4.9), and because of the known behaviour at *r*=1.

### (c) Solution of the dual-integral equations

Expand *δ*(*r*) as equation (4.2) with
4.10Then, from equation (4.3),
4.11using equation (4.9). Substitution in equation (4.7) leads to the integral
this Weber–Schafheitlin integral (Watson 1944, §13.4) vanishes, so that equation (4.7) is satisfied identically.

If we denote the left-hand side of equation (4.6) by *E*_{n}(*r*/*a*), we have
where (Martin 1986, p. 278)
which is a polynomial in *ρ*. Hence, it follows that the right-hand side of equation (4.6) should have a similar polynomial expansion,
4.12and then
4.13This determines from ; in general, these coefficients depend on *s*.

### (d) The pressure field

The Laplace-transformed pressure is given by equations (3.2) and (3.3) (with *C*=0). If we expand *D* as equation (4.3), we obtain
Substituting equations (4.11) and (4.13) gives
4.14where

To conclude, we invert, using
4.15As usual, *p* is supposed to be analytic to the right of Re *s*=*c*>0, ensuring that *p*≡0 for *t*<0. From equation (2.1), we have
4.16with branch points at *s*=±i*N*. We require that Re *λ*≥0 for Re *s*≥0 and we take the cuts to the left, parallel to the real *s*-axis. In addition, for time-harmonic forcing at a frequency *ω*<*N*, the quantities will have poles at *s*=±i*ω* (see equation (4.22) below). Then, moving the inversion contour to the left gives residue contributions from the poles and branch line integrals; the latter decay algebraically, they give the transient behaviour and they are henceforth ignored as we are interested in the ultimate time-harmonic behaviour.

For the residue contributions, we note that
4.17and
4.18where defines *θ*_{c}, 0<*θ*_{c}<*π*/2 and
4.19The substitutions *ξ*→−*ξ* and *η*→−*η* imply that *β*→*β*+*π* hence
4.20where the overbar denotes complex conjugation. Further discussion of the integrals and their far-field behaviour can be found in §7.

To be more specific, suppose that
where *w*_{p} does not depend on the transform variable, *s*. Then, from equation (3.10),
So, if we expand the given function *w*_{p} as (cf. equations (4.5) and (4.12))
4.21with (real) constants , then
4.22Hence, equations (4.14) and (4.16) give
where . Using equations (4.17), (4.18) and (4.20), the sum of the residues of at *s*=±i*ω* is found to be
Hence, inverting gives
4.23where
4.24

Before considering specific forcings and their consequences, we examine the associated time-harmonic problem instead of solving an initial-value problem.

## 5. The time-harmonic disc problem

The pressure is supposed to have the form Re{*p*(*x*,*y*,*z*)e^{−iωt}}, where *p* satisfies
5.1

As in §3, we consider *z*>0, with the plane *z*=0 partitioned into *Ω* (the plate) and *Ω*′, with *p*=0 on *Ω*′ and prescribed normal velocity on *Ω*,

Introduce Fourier transforms,
5.2and
5.3Then, with *κ*^{2}=*ξ*^{2}+*η*^{2}, the Fourier transform of equation (5.1) is *P*′′=−(*κ*/*γ*)^{2}*P*. Hence,
5.4

Equations (5.3) and (5.4) can be interpreted as giving a plane-wave decomposition of *p*. Thus, equation (5.1) has solutions in the form
where . The corresponding group velocity, defined by *c*_{g}=(∂*ω*/∂*k*_{1},∂*ω*/∂*k*_{2},∂*ω*/∂*k*_{3}), is given by
So, in the half-space *z*>0, we may argue that the group velocity should be upwards, away from the disc, implying that we should take *k*_{3}<0, that is, we should take *C*=0 in equation (5.4). Thus, in effect, we have just applied a radiation condition. This idea has been used by Bell (1975), Nycander (2006) and others.

We now proceed as in §3. We define the discontinuity in *p* across *z*=0, *δ*, by equation (3.4); its Fourier transform, *D*, is given by equation (3.5). The Fourier transform of *w* is given by
Applying the boundary conditions gives equations (3.6) and (3.7), involving *W*(*ξ*,*η*,0)=(*γ*/*ω*)*κD* and *P*(*ξ*,*η*,0)=*D*. Hence, we obtain the dual-integral equations, (3.8) and (3.9), for *D*, in which .

For a circular disc, we expand *δ* and *w*_{p} using equations (4.2), (4.10) and (4.21), with known coefficients and unknown coefficients . We have , with related to by equation (4.13):
5.5Then, we compute the pressure, using equations (5.3) and (5.4) with *C*=0, equations (4.3), (4.11) and (5.5): we obtain precisely the formula (4.24).

Summarizing, we solved the time-harmonic problem with the radiation condition of upward-directed group velocity for each constituent plane wave. This led to exactly the same pressure field that we obtained by solving an initial-value problem, enforcing causality and ignoring transient effects.

Henceforth, we focus on the time-harmonic problem.

## 6. Some examples

We consider three simple examples. Two are axisymmetric, including a heaving rigid disc. The third is a rigid disc oscillating about a diameter. For each example, we solve the problem. Later, we shall calculate the far-field radiated field within the conical wave beams.

### (a) The heaving disc

This is the simplest problem. It has been solved by Sarma & Krishna (1972) and by Lai & Lee (1981). We suppose that , a given constant, so that we are only concerned with *n*=0. As , and so equation (4.21) gives , all other being zero.

The pressure is given by equation (4.24) as
with defined by equation (4.19). The pressure discontinuity across the disc is 2*ρ*_{0}*δ*, where
and equation (5.5) gives . The force on the disc is in the vertical direction with magnitude
which agrees with Lai & Lee (1981, eqn (47)).

### (b) Another axisymmetric problem

As , we have . Therefore, from equation (4.21), all other being zero; here, is a constant. The pressure is The force on the disc is readily calculated.

### (c) The rolling disc

We consider a rigid disc oscillating about a diameter taken as the *y*-axis. For this ‘rolling’ motion, we have , where is a given constant. Thus, *n*=1, and as , , all other being zero. The pressure is
and the pressure discontinuity is
where . Thus, there is no net force on the disc but there is a moment about the *y*-axis with magnitude
We are unaware of the previous results on rolling discs.

## 7. The far field

In the far field, is large. However, we are expecting characteristic conical wave beams, so it is useful to introduce conical polar coordinates, *σ* and *ζ*, defined by
7.1equivalently, and . There are two characteristic (double) cones that intersect the edge of the disc at *r*=*a*. Restricting to *z*≥0, *r*≥0, the lower cone is with *r*≥*a*, and the upper cone is . We expect to find energy propagation between these two cones, away from the disc. The lower cone corresponds to and the upper cone corresponds to ; the width of the beam is 2*b* with .

### (a) The far-field pressure

The pressure is given by equation (4.24) in terms of the integral , defined by equation (4.19). As is an even function of *y*, we can assume that *y*≥0.

To start, we use equation (4.1) in equation (4.19), with 0≤*ϕ*≤*π*, giving
7.2The integration with respect to *β* can be evaluated but it turns out to be more convenient to retain the double integral.

In the far field, but *σ* is finite; note that *ζ*∼*R* in this limit. Then, equation (7.2) becomes
7.3where
and
The form of equation (7.3) suggests using the two-dimensional method of stationary phase. We have
These both vanish when *β*=*ϕ*, giving a line of stationary points. On this line, *f*=0, *f*_{κκ}=0, and . Then, from Wong (2001, p. 454, theorem 1), we obtain the estimate
Hence,
7.4where
This integral can be evaluated. From Watson (1944, p. 405, eqns (2) and (3))
7.5Let and . For 0≤*σ*≤*b*, we take the plus in the first of equation (7.5) with . For −*b*≤*σ*≤0, we take the minus in the first of equation (7.5) with . Hence,
7.6with
Formula (7.6) shows the phase variation of *h*_{n}(*σ*) across the wave beam. Outside the beam, we use the second of equation (7.5); for example,

Combining equations (7.4) and (7.6) gives, within the upper beam,
Hence,
7.7where
7.8Notice that *p* decays as *ζ*^{−1/2} for all forcings of the disc, whereas the behaviour across the beam (via *ϑ*) depends on those values of *n* and *m* present in the forcing (see equation (4.21)).

We remark that integrals of the form (7.2) have been considered by Voisin (2003, equation (4.21)). He obtained far-field estimates (see his equation (4.16)) but our method is much more direct.

### (b) The far-field velocity

In Cartesian coordinates, the (time-harmonic) velocity field is ** v**=(

*u*

_{x},

*u*

_{y},

*w*), where

*γ*=

*C*/

*S*, , and

*ω*=

*NC*. We use cylindrical polar coordinates (

*r*,

*ϕ*and

*z*) and conical polar coordinates (

*ζ*,

*σ*and

*ϕ*), defined as follows (see equation (7.1)): Corresponding orthogonal unit vectors are Applications of the chain rule give Within the conical wave beams,

*ζ*is large and

*p*decays as

*ζ*

^{−1/2}(see equation (7.7)), so that Thus, using equation (7.7), 7.9where

*G*

_{n}(

*σ*)=(i

*NS*)

^{−1}

*F*

_{n}′(

*σ*). As , and 7.10

### (c) The far-field energy transport

The time-averaged energy transport vector is (see, for example, Voisin 2003, eqn (4.39))
Let us calculate the total energy flux along the upper conical wave beam,
using (Voisin 2003, eqn (4.40)). Substituting equations (7.7) and (7.9), and integrating over *ϕ* gives
Use of equations (7.8) and (7.10) gives
The remaining integral has the form
where *δ*_{ij} is the Kronecker delta. Hence,
Evidently, this is positive: the total energy flux is always outwards. However, the energy transport vector, **I**, need not point outwards for all *σ* and for all *ϕ*, as we show next using examples.

### (d) Some examples

Define by We calculate for several examples.

#### (i) Axisymmetric forcings

Start with axisymmetric forcing (*n*=0). Then, we have
where

For a heaving rigid disc (§6*a*), and
Similarly, for the forcing described in §6*b*, and
For these two examples, for |*σ*|<*b*. However, suppose we combine the two forcings, retaining both and . Then,
Suppose, further, that , where *α* is a real parameter. Then
where varies from −*π* to *π* as *σ* varies from −*b* to *b*, across the beam. We can easily choose *α* so that is negative in part of the beam; a simple choice is *α*=1.

#### (ii) Non-axisymmetric forcings

For a plate rolling about a diameter (§6*c*), we have *n*=1, and
which is positive in all directions (*ϕ*) and across the beam (|*σ*|<*b*).

Suppose, now, that we combine heaving and rolling; the non-zero coefficients are and . We have
If we suppose, further, that , where *β* is a real constant, then
The sign of this quantity varies with *β* and with position within the beam (via *ϕ* and *ϑ*).

## 8. Discussion

We have solved a class of problems for oscillating horizontal discs. As there is uncertainty about the radiation condition in the frequency domain, we began by solving an initial-value problem in the time domain, where causality can be imposed. We also solved a boundary-value problem in the frequency domain, using plane-wave (Fourier) decompositions with waves selected so as to have outgoing group velocity. Both methods led to the same solution (once transient effects are discarded in the time domain).

We used our time-harmonic solution to solve specific problems, such as those for heaving and rolling rigid discs. It is expected that similar problems for elliptical discs may also be solved (Martin 1986).

We showed how to calculate the far field, using asymptotic methods. In particular, we estimated the pressure and velocity within the conical wave beams, and we showed that the total energy flux is outwards.

- Received March 25, 2011.
- Accepted June 17, 2011.

- This journal is © 2011 The Royal Society