## Abstract

An analysis of wave motion under stochastic excitation is presented. The starting point is a Hamiltonian model of surface wave motion. This model is augmented with linear damping and stochastic forcing terms. An asymptotic scaling parameter is then introduced to show that the problem has three time scales. Integrability of the system is established for the case of two wave modes near 1:1 resonance. Stochastic averaging is used to reduce the dimensions of the system from four to two and the evolution of surface waves is now described, over long time scales and small amplitude stochastic forcing, by the evolution of the integrals of motion, as a random process. The domain of the generator of this random process is characterized and it includes a ‘gluing’ boundary condition. The adjoint of the generator yields the Fokker–Planck partial-differential equation (FPE). This equation governs the time evolution of the joint probability density of the two integrals of motion. The coefficients of the FPE are calculated numerically and the steady-state solution is found using the finite-element method. Results of the analysis show a distinct peak in the probability density along one edge of the reduced domain.

## 1. Introduction

The history of the problem we consider dates back to wave motion experiments performed in the nineteenth century by Faraday (1831). We focus on the motion of waves that form at the free surface of a liquid where gravity is the dominant restoring force, thus *surface-gravity* waves. The work we present is fundamental in nature, but certain applications seem pertinent, for example, in the context of ‘g-jitter’—the small fluctuations in gravitational force introduced by machinery on spacecraft (Walter 1987). Another possible area of application is in the field of liquid sloshing (Ibrahim 2005), where issues of stability can arise. An online report, Schlee (2007), gives an example of the relevance of sloshing in fuel tanks on the stability of spacecraft during their launch.

Our analysis starts with the formulation by Miles (1976), where a Hamiltonian form for the equations governing the motion of surface waves is presented (in a closed basin). We focus on the motion of the two dominant wave modes in 1:1 resonance and excited by horizontal motion, in a sense our analysis is the stochastic counterpart to the work by Miles (1984*b*). We hope that stochastic forcing will lead to more realistic modelling for the application areas mentioned above.

The mathematical analysis we perform falls under the rubric of stochastic averaging methods. Seminal work in stochastic averaging was performed in Khas’minskii (1968). This classical approach has limitations when the system being averaged has more than just one centre fixed point, as is the case in our problem. To perform stochastic averaging in the presence of multiple fixed points, we follow the method of Freidlin & Wentzell (2004), where the underlying approach is based on the martingale problem. The problem we study has three time scales rather than two, thus we extend stochastic averaging by introducing a second averaging operator.

This article is organized as follows. In §2, we present the Hamiltonian model of surface wave motion. In §3, the formalism and application of the stochastic analysis method is presented. In §4, we calculate steady-state probability densities of the Fokker–Planck partial-differential equation (FPE) using the finite-element method (FEM). We conclude with a discussion of our results, in §5.

## 2. Surface-gravity wave model

### (a) Governing equations

Although our ultimate goal is to study surface waves, we begin by considering the ‘internal’ motion of the fluid at whose surface the waves form. The first portion of our presentation reviews the development presented by Miles (1976). For an incompressible fluid whose velocity is described by the vector ** u**, the Boussinesq approximation (Kundu 1990) holds and conservation of mass is expressed by the continuity equation
2.1
For a fluid particle with density

*ρ*, pressure

*p*and kinematic viscosity

*ν*(and with gravity represented by the vector

**), conservation of momentum is expressed by the Navier–Stokes equation. We are working towards a Hamiltonian formulation, however, and therefore discard viscous terms. This means the Navier–Stokes equation is replaced by the Euler equation**

*g*Note, however, that we will reintroduce viscous effects (§2*c*) in the form of linear damping. Next, the fluid motion is assumed irrotational: ∇×** u**=0. This simplification can be inferred from Kelvin’s circulation theorem; it states that for inviscid fluids, rotational motion is conserved. Therefore, if we take the initial state of the fluid to be irrotational, it will remain irrotational. Irrotational fluid motions can be described by a velocity potential

**=∇**

*u**ϕ*. The Euler equation becomes In the above, it has been assumed that gravity acts along the vertical direction, therefore |

**|=**

*g**g*.

In order to relate the motion of the surface waves to the internal motion of the fluid, boundary conditions are needed. First, there is the kinematic boundary condition. As explained in Whitham (1974), derivation of this boundary condition, the kinematic boundary condition, is based on the surface of the fluid being defined by the property that no fluid crosses it. We denote the surface height by *η*
(*x*_{1} and *x*_{2} refer to the coordinates of the horizontal plane, in actual calculations we will use cylindrical coordinates, thus (*x*_{1},*x*_{2})=(*r*,*θ*)), and the kinematic boundary condition is given by
In the above, we have made use of the notation ** u**=(

*u*

_{x1},

*u*

_{x2},

*u*

_{z}). In our treatment, we take the fluid to be confined to a cylindrical tank of radius

*a*and depth

*d*(at rest, the fluid is contained between

*z*=−

*d*and

*z*=0). This leads to boundary conditions at the lateral and bottom walls of the tank that impose zero flow across these walls

Miles showed that the conservation of mass equation, the kinematic boundary condition and the boundary conditions at the lateral and bottom walls can be obtained from a variational formulation. Specifically, the first variation of the integral
must vanish (*S*=*πa*^{2}). Note that this type of variational formulation was first presented by Zakharov (1968).

Miles’s analysis proceeds with separation of variables
and
The relationship between *ψ*_{n} and *χ*_{n} is found from the conservation of mass equation, the boundary condition at the container’s walls and a linearization of the kinematic boundary condition
The relation between *ϕ*_{n}(*t*) and *q*_{n}(*t*) is found using the variational formulation. Details are given by Miles (1976, §2).

Specifically, for a cylindrical container, as given by Miles (1984*a*), the eigenfunctions *ψ*_{n} contain Bessel functions
where *N*_{ij} is a normalization factor and *J*_{i} is a Bessel function and the eigenvalues satisfy

### (b) The Hamiltonian of surface-gravity waves

Averaging techniques reduce the dimensions of a system. The lower dimensional system’s long-term evolution is given in terms of quantities that are constants of motion of the original system, over short time spans. This motivates our interest in determining the Hamiltonian for the surface-gravity wave system. First, we write expressions for the kinetic and potential energies. These give the Lagrangian from which the Hamiltonian follows. The kinetic energy is
The second line makes use of the results by Miles (1976), where it is shown that
Denoting the natural frequency of the *n*th mode by *ω*_{n}, . The calculations of the correlation integrals defining *a*_{lmn} and *a*_{jlmn}, specific to cylindrical basins, are provided in appendix A of Miles (1984*a*). These quantities depend on the depth and radius of the cylindrical basin. The potential energy of the free surface, displaced from its equilibrium position, is
2.2
The Lagrangian, normalized by *ρS*, is
and the Hamiltonian is
where
Note that this Hamiltonian contains quadratic terms, but also higher order terms as the coefficient *h*_{mn} depends on *q*_{n}
where and equations for *h*_{lmn} and *h*_{jlmn} are available in §4 of Miles (1976).

### (c) Damping and forcing effects

We are interested in the effects of stochastic forcing. To balance forcing effects, damping is necessary. Physically, damping stems from the viscosity of the fluid, but for simplicity we introduce damping ad hoc. We use linear damping coefficients, *α*_{k}, to supplement the momentum equation.

To account for forcing effects, terms are added to the potential energy given in equation (2.2), thus
In the above, *ξ*_{x1} and *ξ*_{x2} specify the imposed horizontal velocities of the tank while *ξ*_{z} is the vertical velocity and
Note that the potential energy could be modified to include surface tension. Details are provided in Miles & Henderson (1990).

With forcing and damping, the equations governing the motion of the surface-gravity waves are, for , and

## 3. Stochastic analysis

In this section, we analyse the equations of motion presented in the previous section under the influence of stochastic forcing. We concentrate on the motion of two wave modes near 1:1 resonance and with stochastic forcing imposed in the horizontal direction. Mathematically, we will define the domain of the generator of the stochastic process that is central to our analysis.

### (a) Integrals of motion

In order to study the effects of small amplitude noise over long time scales, we introduce a scaling parameter *ϵ* (0<*ϵ*≪1) and a canonical change of variables
with *ω*_{1,2}=*ω*+*ϵσ*_{1,2}. This transformation eliminates terms of order *ϵ* in the Hamiltonian. The terms of order *ϵ*^{3/2} vanish under the averaging operator
3.1
Terms of order *ϵ*^{3/2} can contribute at order *ϵ*^{2}; however, these are second-order effects that do not qualitatively change the underlying 1:1 resonant normal form. In our analysis, these contributions are not considered.

The dynamics of the system acquire the form
3.2
where ** x**=(

*x*

_{1},

*y*

_{1},

*x*

_{2},

*y*

_{2}). Equation (3.2) is helpful in showing the three time scales present in this problem. Owing to their small amplitude (i.e.

*ϵ*

^{2}), the effects of noise and damping have an influence only over long times, of order 1/

*ϵ*

^{2}. The Hamiltonian dynamics, associated with the components of

*b*^{1}, are of order

*ϵ*, and thus they fluctuate over a faster time scale, of order 1/

*ϵ*. Finally, the fastest fluctuations occur owing to periodicity of some of the coefficients that constitute

*b*^{1}. The coefficients fluctuate over a period

*ω*/2

*π*, which is of order 1.

We rescale time such that leading order dynamics become of order 1. With , equation (3.2) becomes
With the equation in this form, we are able to apply a result from Khas’minskii (1966*b*). This result tells us that as *ϵ*→0, converges in probability to
In the above, denotes the symplectic gradient
The averaged Hamiltonian, , is our first integral of motion. Written out explicitly,
The second integral of motion will be denoted by
Note that this integral is found by using a canonical transformation that uses the sum and difference of two angles. In this way, *I* is interpreted as an angular momentum. The Poisson bracket of and *I* vanishes, making these two quantities independent; we have therefore obtained the two integrals of motion of our system.

### (b) Structure of the unperturbed system

As stated in the previous section, the surface wave system in 1:1 resonance has two constants of motion, and *I*. We now introduce a canonical transformation that takes advantage of this fact by making *I* one of the canonical variables. The canonical transformation from (*x*_{1},*y*_{1},*x*_{2},*y*_{2}) to (*X*,*Y*,*θ*,*I*) is
Introducing the variables
the averaged Hamiltonian becomes
3.3
and equations of motion are
This canonical transformation decouples the equations for *X*_{t} and *Y* _{t} from those for *θ*_{t} and, consistent with results of §3*a*, *I*_{t} is a constant. Our calculations will therefore be performed in the *X*–*Y* plane and *I* will act as a bifurcation parameter. Defining the critical value
when *I*<*I*_{c} phase portraits in the *X*,*Y* plane have a single elliptic fixed point at the origin. For *I*>*I*_{c}, the fixed point at the origin becomes a saddle point and two fixed points on the *X* axis appear. The coordinates of these fixed points are

The stochastic diffusion process that we will study is described in terms of the variables and *I*. The –*I* domain is shown in figure 1. While the figure shows a maximum value of *I*=0.5, this value could be made arbitrarily large. The range of is not arbitrary. The maximum value, , occurs at the elliptic fixed points. There is also a minimum value, , as the canonical transformation introduces the restriction . This leads to a Hamiltonian surface that has two hill-tops, at each of the elliptic fixed points. Our equations yield hill-tops rather than valleys, but also negative damping, making the system consistent. Note that the domain between and exits for both elliptic fixed points; thus, the entire –*I* domain may be described as a three-leaved ‘open-book’, a nomenclature used in Freidlin & Wentzell (2004).

### (c) Evolution of the integrals of motion

The ultimate goal of our analysis is to understand the effects of stochastic forcing. The amplitude of the forcing is small and its effects are seen by examining the evolution of the two conserved quantities, and *I*, over long times. Thus, we define the map by
As ,
The evolution of is given by
where
As and *I* are integrals of motion, . Thus, to see the fluctuations of and *I*, we need to look on a time scale of order 1/*ϵ*^{2}. Thus, we rescale time again, setting and . This yields
3.4

### (d) Martingale problem

In order to carry out our stochastic analysis, we must rigorously characterize the dynamics of equation (3.4). The classical analysis of such equations dates back to Khas’minksii. As shown in §3*b*, the reduced domain (i.e. ) of the surface wave equations contains a saddle fixed point. For such a case, as opposed to a domain free of fixed points in its interior, showing that is a Markov process (where it must be stressed that by this we mean alone, without being coupled to the equations) is done with a martingale problem approach.

The martingale approach begins by considering rather than . The martingale problem is stated in terms of a canonical process *X*_{t}. This canonical process is made of deterministic parts (the Hamiltonian dynamics and the damping in our case) and stochastic parts (forcing.) The stochastic part has a probability measure and this is also the ‘induced’ probability for . *X*_{t} on the other hand has a probability measure and the correspondence between the two probability measures is given, in terms of a set *A*, as
This transfers the dependence on *ϵ* from the process to the measure.

The martingale problem states that the measure is unique and that it has a martingale associated with it. The process ‘lives’ in a reduced two-dimensional space. Denoting the reduced-space canonical process corresponding to by and the probability measure of by the novel requirements for the class of problems within which the stochastic analysis of the surface waves falls are to show the existence of the limit

The majority of the results needed to prove this are given by Freidlin & Wentzell (2004); however, their results do not account for averaging. Here, we simply point out that in order to treat a domain with fixed points, the martingale problem approach is required and using this approach, classical results are augmented with gluing conditions.

### (e) Generator of the reduced Markov process

In the previous section, we stated the existence of a reduced process with a probability measure . In this section, we characterize the Markov process in detail.

Let us denote by the three-leaf space depicted in figure 1. To begin, we decompose into three disjoint regions: two inner leaves for the regions enclosed by the homoclinic orbit and the third outer leaf for the region outside the homoclinic orbit

Within each leaf, , classical stochastic averaging results (e.g. Khas’minskii 1968) hold, and thus the generator of the reduced Markov process is, for a function ,
where is the drift coefficient and the diffusion coefficient. In the present problem, because the equations of motion involve time derivatives of noise, we apply a ‘real’ noise with a specified spectral density rather than white noise (as the latter cannot be differentiated.) The results of Khas’minskii (1966*a*) yield the following:
3.5and
3.6
where
3.7
and
3.8

Note that Khas’minskii’s definition of includes an additional term known as the ‘second-averaging’ term. We do not include this term here because it may change our results quantitatively but not qualitatively.

is a spatial averaging operator. It is applied along orbits *S* on which and *I* are constant
3.9
is a measure of speed, thus the denominator corresponds to the period, , associated with a Hamiltonian orbit. At fixed points and on the homoclinic orbit, velocity is zero, thus the definition above does not hold. On the homoclinic orbit, while at the elliptic fixed point the period is finite and equal to the period of the system linearized about the fixed point. At the saddle point, we will need to evaluate drift and diffusion coefficients without dividing by the (infinite) period. For notational convenience, we introduce

Thus, each has a distinct generator . The generators are linked along the edge by a ‘gluing’ boundary condition. For a process that does not ‘stick’ to any of the orbits, which is true for the surface waves, the gluing condition is
The + sign is taken if on the leg and the − sign is taken otherwise. *ν*_{1} and *ν*_{2} are unit vectors in the directions of and *I*, respectively. corresponds to the gluing vertex where the three leaves meet (and ).

With the gluing condition specified, the domain of the generator of the process evolving on can be defined
3.10
The first condition specifies behaviour at the boundary where *X*^{2}+*Y* ^{2}=2*I*, the outer leaf boundary. The limit as specifies behaviour at the two centre fixed points. The third condition holds in the limit .

### (f) Calculation of the drift and diffusion coefficients

The drift vector, , and diffusion matrix, , are defined by equations (3.5) and (3.6). These equations involve two averaging operators, averaging, which is defined by equation (3.1) and averaging, defined by equation (3.9). Up to and including the averaging operator, we are able to carry out our calculations analytically, whereas the final operation, averaging, is performed numerically. In this section, we explain how and are calculated.

The expectation operator used in equations (3.7) and (3.8) leads to the introduction of correlation functions. In our case, with horizontal forcing in a single direction (i.e. *ξ*_{z}=*ξ*_{x2}=0),
Furthermore, using a Fourier transform
Applying these two simplifications gives formulas for and in terms of *X*,*Y* and *I*. These formulas are given in appendix A.

We turn to numerical integration. The integrands of the drift and diffusion coefficients have been calculated in Cartesian coordinates (*X*,*Y*). Standard numerical integrators are programmed for integrals of the form
To cast integration over Hamiltonian orbits, as expressed in equation (3.9) into this form, we parametrize the Hamiltonian orbits by the angle variable of polar coordinates. This parametrization by *θ* yields two cases:

— Orbits for which

*I*<*I*_{c}or . These orbits are centred at the origin, as illustrated in figure 2.— Orbits for which . These orbits are centred at either the positive or negative elliptic fixed point, illustrated in figure 3.

Simplifying the integrands appearing in equation (3.9) to *g*(*X*,*Y*) and denoting the coordinate of the centre fixed point by (*x*_{c},0), transformation to polar coordinates and parametrization by *θ* gives
3.11
(For case 1, *x*_{c}=0.) To apply equation (3.11), *r*(*θ*) is calculated with a numerical root solver by transforming equation (3.3) to polar coordinates.

The procedure described above holds in the –*I* domain wherever the period is not infinite. The period is infinite along the edges and . Along the edge , the diffusion coefficients have a finite value whereas the drift coefficients do not, but only the former need to be evaluated to impose the conservation of probability flux condition. Along the edge , the drift and diffusion coefficients are found by linearization about the coordinates of the elliptic fixed points. This allows us to calculate and everywhere needed. Sample representations of the drift and diffusion coefficients are shown in figure 4. The figures show that the drift vector can be both negative and positive whereas the elements of the diffusion matrix on the diagonal are non-negative. Note also that the diffusion matrix is positive-definite.

## 4. Steady-state probability densities

### (a) The Fokker–Planck equation

Having defined the reduced Markov process and its generator for our problem, we apply our results to calculate steady-state probability densities.

The time evolution of probability density is governed by the FPE; this can be demonstrated by requiring that the generator be self-adjoint, see Sri Namachchivaya & Sowers (2001). Furthermore, self-adjointness dictates what boundary conditions must be imposed on the FPE equation.

We restrict our analysis to the steady-state FPE. This is justified because our formulation is for the dynamics of two wave modes in resonance, a regime observed after transients have died out. The steady-state FPE is
We introduce the *probability flux*, as it allows us to write some of the equations that follow more compactly. The probability flux is
With *J*, the FPE equation can be written
showing its similarity to the conservation of mass equation in continuum mechanics.

Boundary conditions consistent with the domain of the generator given in equation (3.10) follow.

— At the vertex , conservation of probability flux 4.1

— Along the edges and , zero probability flux is imposed to create a reflection boundary

— At the upper edge

*I*_{max}, the boundary condition requires special care. The next section discusses this.

### (b) The boundary condition at *I*_{max}

Establishing the correct boundary condition to impose at *I*_{max} must be done carefully. When characterizing the generator of the Markov process in the context of a martingale problem (§3*e*), proofs are carried out by first selecting a finite value for *I*_{max} and killing the process at that value of *I*. After completing the proofs, *I*_{max} is allowed to increase to infinity.

In the case of a steady-state numerical solution, imposing a killing condition for a steady-state solution at finite value of *I* would simply lead to the trivial solution (i.e. *p*(*h*,*i*)=0) which is not satisfactory. Given this, we choose to impose a reflecting boundary condition at *I*_{max}. In doing so, we must keep in mind that the value of *I*_{max} may have an effect on the solution of the FPE.

Note that were we to solve the time-dependent FPE, applying the killing condition for a finite *I*_{max} (and obtaining a non-trivial solution) may be possible.

### (c) Finite-element solution of the FPE

In certain cases, it is possible to solve the steady-state FPE analytically (Soize 1994), but this is not expected here; the geometry within which our problem is stated is not particularly simple and we only know how to evaluate the drift and diffusion coefficients for our problem numerically. Therefore, we look for a solution using the FEM. The two texts consulted to implement the finite-element solution were Becker *et al.* (1981) and Pozrikidis (2005). Our implementation is programmed in Octave (Eaton 2002).

If *v*(** y**) represents an arbitrary test function (this is part of the FE formulation), the weak-form solution of the FPE is
We simplify this equation by discarding the time derivative as we seek only steady-state solutions. In the FEM, the second-order derivatives are changed to first-order derivatives using integration by parts. Thus, we obtain
It is worth noting that only Neumann boundary conditions are imposed in our problem; a unique solution is obtained as a consequence of the normalization condition.

Following the FEM, *v* and *p* are approximated by a finite number of basis functions *ϕ*_{1},…,*ϕ*_{N} so that
The finite element form of the surface integral is, denoting, ** y**=(

*x*,

*y*) The

*β*

_{i}, being coefficients of arbitrary test functions, the FEM proceeds by choosing the set of

*N*test functions such that the first test function has

*β*

_{1}=1 and all the other

*β*

_{i}=0, the second test function has only

*β*

_{2}=1 as the non-zero

*β*

_{i}and so forth, thus the equation we obtain is To find steady-state solutions to the FPE, the null space of the matrix

*K*needs to be determined. Numerically, the null space is found by calculating the singular value decomposition of

*K*The column of

*V*associated with the smallest singular value is normalized and this is the steady-state solution to the FPE.

Note that the conservation of probability flux condition, equation (4.1), implies *C*^{1} continuity across the gluing vertex. This degree of continuity can be satisfied automatically by treating all the leaves on which the FPE is posed as a single domain.

## 5. Discussion and conclusions

Figure 5 shows a solution to the FPE produced with the FEM. The solution shows that the probability density function (PDF) is highest along the edge . This feature of the solution remains present when the aspect ratio of the cylindrical basin is varied from *d*/*a*=1, as in figure 5 down to 0.3 (solution not illustrated). In the X–Y phase plane, the edge is situated outside the homoclinic orbit and our results suggest that the wave motions most likely to be observed would be those corresponding to circular motion where for a given value of *I*, is near its maximum.

Similarly to our analysis, in Miles (1984*c*), a system of two surface-gravity wave modes in a cylindrical basin subjected to horizontal excitation is studied. Unlike our problem, the two wave modes are set to be in exact resonance with one another, and a detuning parameter is introduced to represent the frequency difference between the natural wave frequency and the forcing frequency. A significant conclusion reached in that study is that for cylindrical basins with an aspect ratio between 0.3 and 0.5, limit cycles and chaotic motion are not possible, making the wave dynamics observed in that interval of aspect ratios qualitatively different from the dynamics seen outside that interval. In our case, with stochastic excitation, the probability densities obtained inside and outside the interval do not exhibit any qualitative differences; in both cases, the probability density is highest near the edge .

In conclusion, we started with a Hamiltonian model of surface wave motion. We then analysed the stochastic dynamics of this system by applying the method of stochastic averaging, paying special attention to the three time scales present in the original system as well as the bifurcations in the phase space of the system. Finally, we demonstrated that the system dynamics can be described with the FPE and we solved that equation numerically, with the FEM.

We believe the work we have presented lays the foundation for physically motivated investigations of wave motion. For example, our results may be extended so as to relate them to research on wave evolution in systems with discrete spectra (Lvov *et al.* 2006; Kartashova *et al.* 2008).

## Acknowledgements

The authors would like to acknowledge the support of the AFOSR under grant no. FA9550-08-1-0206 and the National Science Foundation under grant no. CMMI 07-58569. Any opinions, findings and conclusions or recommendations expressed in this paper are those of the authors and do not necessarily reflect the views of the National Science Foundation.

## Appendix A. Drift and diffusion coefficient integrands

A1 A2 A3 A4 A5 A6and A7

## Footnotes

- Received December 22, 2009.
- Accepted February 8, 2010.

- © 2010 The Royal Society