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.
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 (1984b). 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 g), 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
Note, however, that we will reintroduce viscous effects (§2c) 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 η (x1 and x2 refer to the coordinates of the horizontal plane, in actual calculations we will use cylindrical coordinates, thus (x1,x2)=(r,θ)), and the kinematic boundary condition is given by In the above, we have made use of the notation u=(ux1,ux2,uz). 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=πa2). 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 qn(t) is found using the variational formulation. Details are given by Miles (1976, §2).
Specifically, for a cylindrical container, as given by Miles (1984a), the eigenfunctions ψn contain Bessel functions where Nij is a normalization factor and Ji 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 nth mode by ωn, . The calculations of the correlation integrals defining almn and ajlmn, specific to cylindrical basins, are provided in appendix A of Miles (1984a). 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 hmn depends on qn where and equations for hlmn and hjlmn 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=(x1,y1,x2,y2). 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 b1, 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 b1. 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 (1966b). 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 (x1,y1,x2,y2) 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 Xt and Yt from those for θt and, consistent with results of §3a, It 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<Ic phase portraits in the X,Y plane have a single elliptic fixed point at the origin. For I>Ic, 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 §3b, 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 Xt. 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 . Xt 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 (1966a) 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 X2+Y2=2I, 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<Ic 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 (xc,0), transformation to polar coordinates and parametrization by θ gives 3.11 (For case 1, xc=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 Imax, the boundary condition requires special care. The next section discusses this.
(b) The boundary condition at Imax
Establishing the correct boundary condition to impose at Imax must be done carefully. When characterizing the generator of the Markov process in the context of a martingale problem (§3e), proofs are carried out by first selecting a finite value for Imax and killing the process at that value of I. After completing the proofs, Imax 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 Imax. In doing so, we must keep in mind that the value of Imax 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 Imax (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 C1 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 (1984c), 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).
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
- Received December 22, 2009.
- Accepted February 8, 2010.
- © 2010 The Royal Society