Generalized Eigenfunction Expansion for Water Waves for a Fixed Body

This is a complete page

Introduction

The generalized eigenfunction method goes back to the work of Povzner 1953 and Ikebe 1960. The generalized eigenfunction method has been applied to water-wave problems by Friedman and Shinbrot 1967, Hazard and Lenoir 2002 (for the case of a rigid body in water of infinite depth) and Meylan 2002 (for a thin plate on water of shallow draft). We will present here the theory for a rigid body in water of finite depth. rms of this solutions (which we call the generalised eigenfunctions) because they solve for

The generalized eigenfunction method is based on an inner product in which the evolution operator is self-adjoint. It follows from the self-adjointness that we can expand the solution in the eigenfunctions of the operator. These eigenfunctions are nothing more than the single-frequency solutions. The main difficulty is that the eigenfunctions are associated with a continuous spectrum, and this requires that they be carefully normalized. Once this is done, we can derive simple expressions which allow the time-domain problem to be solved in terms of the single-frequency solutions. The mathematical ideas are discussed in detail in Hazard and Loret 2007.

We consider a two-dimensional fluid domain of constant depth, which contains a finite number of fixed bodies of arbitrary geometry. We denote the fluid domain by [math]\displaystyle{ \Omega }[/math], the boundary of the fluid domain which touches the fixed bodies by [math]\displaystyle{ \partial\Omega }[/math], and the free surface by [math]\displaystyle{ F. }[/math] The [math]\displaystyle{ x }[/math] and [math]\displaystyle{ z }[/math] coordinates are such that [math]\displaystyle{ x }[/math] is pointing in the horizontal direction and [math]\displaystyle{ z }[/math] is pointing in the vertical upwards direction (we denote [math]\displaystyle{ \mathbf{x}=\left( x,z\right) ). }[/math] The free surface is at [math]\displaystyle{ z=0 }[/math] and the sea floor is at [math]\displaystyle{ z=-h }[/math]. The fluid motion is described by a velocity potential [math]\displaystyle{ \Phi }[/math] and free surface by [math]\displaystyle{ \zeta }[/math].

The equations of motion in the time domain are Laplace's equation through out the fluid

[math]\displaystyle{ \Delta\Phi\left( \mathbf{x,}t\right) =0,\ \ \mathbf{x}\in\Omega. }[/math]

At the bottom surface we have no flow

[math]\displaystyle{ \partial_{n}\Phi=0,\ \ z=-h. }[/math]

At the free surface we have the kinematic condition

[math]\displaystyle{ \partial_{t}\zeta=\partial_{n}\Phi,\ \ z=0,\ x\in \partial\Omega_{F}, }[/math]

and the dynamic condition (the linearized Bernoulli equation)

[math]\displaystyle{ \partial_{t}\Phi = -g\zeta ,\ \ z=0,\ x\in \partial\Omega_{F}. }[/math]

The body boundary condition for a fixed body is

[math]\displaystyle{ \partial_{n}\Phi=0,\ \ \mathbf{x}\in\partial\Omega, }[/math]

The initial conditions are

[math]\displaystyle{ \left.\zeta\right|_{t=0} = \zeta_0(x)\,\,\,\mathrm{and}\,\,\, \left.\partial_t\zeta\right|_{t=0} = \partial_t\zeta_0(x). }[/math]

Equations in the Frequency Domain

We also assume that Frequency Domain Problem with frequency [math]\displaystyle{ \omega }[/math] and we assume that all variables are proportional to [math]\displaystyle{ \exp(-\mathrm{i}\omega t)\, }[/math]

The water motion is represented by a velocity potential which is denoted by [math]\displaystyle{ \phi\, }[/math] so that

[math]\displaystyle{ \Phi(\mathbf{x},t) = \mathrm{Re} \left\{\phi(\mathbf{x},\omega)e^{-\mathrm{i} \omega t}\right\}. }[/math]

The equations are the following

[math]\displaystyle{ \begin{align} \Delta\phi &=0, &-h\lt z\lt 0,\,\,\mathbf{x} \in \Omega \\ \partial_z\phi &= 0, &z=-h, \\ \partial_z \phi &= \alpha \phi, &z=0,\,\,\mathbf{x} \in \partial \Omega_{\mathrm{F}}, \end{align} }[/math]

(note that the last expression can be obtained from combining the expressions:

[math]\displaystyle{ \begin{align} \partial_z \phi &= -\mathrm{i} \omega \zeta, &z=0,\,\,\mathbf{x} \in \partial \Omega_{\mathrm{F}}, \\ \mathrm{i} \omega \phi &= g\zeta, &z=0,\,\,\mathbf{x} \in \partial \Omega_{\mathrm{F}}, \end{align} }[/math]

where [math]\displaystyle{ \alpha = \omega^2/g \, }[/math]) The body boundary condition for a rigid body is just

[math]\displaystyle{ \partial_{n}\phi=0,\ \ \mathbf{x}\in\partial\Omega_{\mathrm{B}}, }[/math]

The equation is subject to some radiation conditions at infinity. We assume the following. [math]\displaystyle{ \phi^{\mathrm{I}}\, }[/math] is a plane wave travelling in the [math]\displaystyle{ x }[/math] direction,

[math]\displaystyle{ \phi^{\mathrm{I}}(x,z)=A \phi_0(z) e^{\mathrm{i} k x} \, }[/math]

where [math]\displaystyle{ A }[/math] is the wave amplitude (in potential) [math]\displaystyle{ \mathrm{i} k }[/math] is the positive imaginary solution of the Dispersion Relation for a Free Surface (note we are assuming that the time dependence is of the form [math]\displaystyle{ \exp(-\mathrm{i}\omega t) }[/math]) and

[math]\displaystyle{ \phi_0(z) =\frac{\cosh k(z+h)}{\cosh k h} }[/math]

In two-dimensions the Sommerfeld Radiation Condition is

[math]\displaystyle{ \left( \frac{\partial}{\partial|x|} - \mathrm{i} k \right) (\phi-\phi^{\mathrm{{I}}})=0,\;\mathrm{{as\;}}|x|\rightarrow\infty\mathrm{.} }[/math]

where [math]\displaystyle{ \phi^{\mathrm{{I}}} }[/math] is the incident potential.

Time domain calculations

The solution in the frequency domain can be used to construct the solution in the time domain. This is well-known for the case of a plane incident wave which is initially far from the body, and in this case the solution can be calculated straightforwardly using the standard Fourier transform. However, when we consider an initial displacement which is non-zero around the bodies, we cannot express the time-dependent solution in terms of the single frequency solutions by a standard Fourier transform.

We begin with the equations in the time domain subject to the initial conditions given by Denoting the potential at the surface by

[math]\displaystyle{ \psi(x,t)=\left.\Phi\left(\mathbf{x},t\right)\right|_{z=0}, }[/math]

we introduce the operator [math]\displaystyle{ \mathbf{G} }[/math] which maps the surface potential to the potential throughout the fluid domain. The operator [math]\displaystyle{ \mathbf{G}\psi }[/math] is found by solving

[math]\displaystyle{ \begin{matrix} \Delta\Psi\left( \mathbf{x}\right) & = 0,\ \ \mathbf{x}\in\Omega,\\ \partial_{n}\Psi & = 0,\ \ \mathbf{x}\in\partial\Omega,\\ \partial_{n}\Psi & = 0,\ \ z=-h,\\ \Psi & =\psi,\ \ z = 0,\ x\in F, \end{matrix} }[/math]

and is defined by [math]\displaystyle{ \mathbf{G}\psi=\Psi. }[/math] The operator [math]\displaystyle{ \partial_{n}\mathbf{G} }[/math], which maps the surface potential to the normal derivative of potential at the surface (called the Dirichlet-to-Neumann map) is given by

[math]\displaystyle{ \partial_{n}\mathbf{G}\psi=\left. \partial_{n}\Psi\right\vert _{z=0}. }[/math]

Therefore the equations in the time domain can be written as

[math]\displaystyle{ \partial_{t}^2 \zeta + \partial_{n}\mathbf{G} \zeta = 0. }[/math]

where we can recover the potential using the operator [math]\displaystyle{ \mathbf{G} }[/math]. The evolution operator [math]\displaystyle{ \partial_{n}\mathbf{G} }[/math] is symmetric in the Hilbert space given by the following inner product

[math]\displaystyle{ \left\langle \zeta,\eta \right\rangle _{\mathcal{H}}= \int_{F}\zeta \eta ^{*}\,\mathrm{d} x, }[/math]

where [math]\displaystyle{ * }[/math] denotes complex conjugate, and we assume that this symmetry implies that the operator is self-adjoint. We can prove the symmetry by using Green's second identity

[math]\displaystyle{ \begin{matrix} \int_{F}\left( \partial_{n}\mathbf{G}\zeta \right) \left( \eta\right) ^{*}\,\mathrm{d} x &=\int_{F}\left( \partial_{n}\mathbf{G}\partial_t\phi \right) \left( \partial_t\psi\right) ^{*}\,\mathrm{d} x \\ &=\int_{F}\left( \partial_t\phi \right) \left( \partial_{n}\mathbf{G} \partial_t\psi\right) ^{*}\,\mathrm{d} x \\ &=\int_{F}\left( \zeta \right) \left( \partial_{n}\mathbf{G}\eta\right) ^{*}\,\mathrm{d} x, \end{matrix} }[/math]

where [math]\displaystyle{ \phi }[/math] and [math]\displaystyle{ \psi }[/math] are the surface potentials associated with [math]\displaystyle{ \zeta }[/math] and [math]\displaystyle{ \eta }[/math] respectively.

Eigenfunctions of [math]\displaystyle{ \partial_{n}\mathbf{G} }[/math]

The eigenfunctions of [math]\displaystyle{ \partial_{n}\mathbf{G} }[/math] satisfy

[math]\displaystyle{ \partial_{n}\mathbf{G} \zeta = \omega^2 \zeta. }[/math]

To solve for the eigenfunctions of [math]\displaystyle{ \partial_{n}\mathbf{G} }[/math] we need to solve the frequency-domain equations, and the radian frequency [math]\displaystyle{ \omega }[/math] is exactly the eigenvalue. To actually calculate the eigenfunctions of [math]\displaystyle{ \partial_{n}\mathbf{G} }[/math] we need to specify the incident wave potential, and for each frequency we have two eigenfunctions (waves incident from the left and from the right). It is possible for there to exist point spectra for this operator which correspond to the existence of a trapped mode McIver 1996 and the presence of a trapped mode requires that the generalized eigenfunction expansion we derive must be modified.

Normalization of the Eigenfunctions

The eigenfunctions of [math]\displaystyle{ \partial_{n}\mathbf{G} }[/math] (with eigenvalue [math]\displaystyle{ \omega }[/math]) are denoted by [math]\displaystyle{ \zeta_{\kappa}(x,k\left( \omega\right) ) }[/math], where [math]\displaystyle{ \kappa = \pm 1 }[/math] for waves incident from the left or right respectively. We must also specify radiations conditions to normalize the eigenfunctions and they are given by

[math]\displaystyle{ {\zeta}_1 = \left( \mathrm{e}^{\mathrm{i} kx} + R_1 e^{- \mathrm{i} k x} \right)\frac{\cos k_0\left( z+h\right) }{\cos k_0h}, \,\,\,\mathrm{as}\,\,x\to-\infty, }[/math] [math]\displaystyle{ {\zeta}_1= T_1\mathrm{e}^{\mathrm{i} k x}\frac{\cos k_0\left( z+h\right) }{\cos k_0h}, \,\,\,\mathrm{as}\,\,x\to\infty, }[/math] [math]\displaystyle{ {\zeta}_{-1} = T_{-1}\mathrm{e}^{-\mathrm{i} k x} \frac{\cos k_0\left( z+h\right) }{\cos k_0h}, \,\,\,\mathrm{as}\,\,x\to\infty, }[/math] [math]\displaystyle{ {\zeta}_{-1}= \left( \mathrm{e}^{-\mathrm{i} k x} + R_{-1} e^{\mathrm{i} k x} \right)\frac{\cos k_0\left( z+h\right) }{\cos k_0 h}, \,\,\,\mathrm{as}\,\,x\to-\infty. }[/math]

Note that [math]\displaystyle{ R_{\kappa} }[/math] and [math]\displaystyle{ T_{\kappa} }[/math] are the reflection and transmission coefficients respectively and that we have normalized so that the amplitude (in displacement) is unity. Also [math]\displaystyle{ k = - \mathrm{i} k_0 }[/math].

and we will consider both [math]\displaystyle{ k(\omega) }[/math] and [math]\displaystyle{ \omega(k) }[/math] in what follows as required. The solution of the single-frequency equation may be computationally challenging and for the generalized eigenfunction expansion the major numerical work is to determine the single-frequency solutions. As mentioned previously, determining [math]\displaystyle{ \zeta_\kappa }[/math] is the major computation of the generalized eigenfunction method, but we simply assume that they are known. We know that the eigenfunctions are orthogonal for different [math]\displaystyle{ \omega }[/math] (from the self-adjointness of [math]\displaystyle{ \partial_n\mathbf{G} }[/math]), and that the waves incident from the left and right with the same [math]\displaystyle{ \omega }[/math] are orthogonal from the identity

[math]\displaystyle{ R_1 T_{-1}^{*} + R_{-1}^{*} T_{1}=0, }[/math]

Mei 1983. It therefore follows that

[math]\displaystyle{ \left\langle \left( {\zeta}_{\kappa}(x,k\left( \omega_{1}\right) )\right) ,{\zeta}_{\kappa^{\prime}}(x,k\left( \omega_{2}\right) )\right\rangle _{\mathcal{H} }=\Lambda_{n}\left( \omega_{1}\right) \delta\left( \omega_{1}-\omega _{2}\right) \delta_{\kappa\kappa^{\prime}}, }[/math]

but we need to determine the normalizing function [math]\displaystyle{ \Lambda_{n}\left( \omega_{n}\right) }[/math]. This is achieved by using the result that the eigenfunctions satisfy the same normalizing condition with and without the scatterers present. This result, the proof of which is quite technical, is well-known and has been shown for many different situations. The original proof was for Schr\"odinger's equation and was due to Povzner 1953 and Ikebe 1960. A proof for the case of Helmholtz equation was given by Wilcox 1975. Recently the proof was given for water waves by Hazard and Lenoir 2002 and Hazard and Loret 2007.

Since the eigenfunctions satisfy the same normalizing condition with and without the scatterers, we normalize with the scatterers absent. This means that the eigenfunctions are simply the incident waves, and the free surface [math]\displaystyle{ F }[/math] is the entire axis. This allows us to derive

[math]\displaystyle{ \begin{matrix} \left\langle \left( {\zeta}_{\kappa}(x,k\left( \omega_{1}\right) )\right) ,{\zeta}_{\kappa^{\prime}}(x,k\left( \omega_{2}\right) )\right\rangle _{\mathcal{H}} &= \int_{\mathbb{R}}\left( e^{\kappa\mathrm{i} k_{1} x}\right) \left( e^{\kappa^{\prime}\mathrm{i} k_{2}x}\right) ^{*}\,\mathrm{d} x \\ & =2\pi \delta_{\kappa\kappa^{\prime}}\delta\left( k_{1}-k_{2}\right) \\ & =2\pi\delta_{\kappa\kappa^{\prime}}\delta\left( \omega_{1}-\omega_{2}\right) \left. \frac{\mathrm{d}\omega}{\mathrm{d}k}\right\vert _{\omega=\omega_{1}}. \end{matrix} }[/math]

This result allows us to calculate the time-dependent solution in the eigenfunctions (or single-frequency solutions).

Expansion in Eigenfunctions

We expand the solution for the displacement in the time domain as

[math]\displaystyle{ \zeta\left( x,t\right) =\int_{\mathbb{R}^{+}} \sum_{\kappa\in\left\{ -1,1\right\}} \left\{ f_{\kappa}\left( \omega\right) \cos(\omega t)+ g_{\kappa}\left( \omega\right) \frac{\sin(\omega t)}{\omega}\right\} \zeta_{\kappa}(x,k) \mathrm{d}\omega, }[/math]

where [math]\displaystyle{ f_\kappa }[/math] and [math]\displaystyle{ g_\kappa }[/math] will be determined from the initial conditions. Note that here, and in subsequent equations, we are assuming that [math]\displaystyle{ k }[/math] is a function of [math]\displaystyle{ \omega }[/math] or that [math]\displaystyle{ \omega }[/math] is a function of [math]\displaystyle{ k }[/math] as required. If we take the inner product with respect to the eigenfunctions [math]\displaystyle{ \zeta_\kappa }[/math] we obtain

[math]\displaystyle{ \left\langle \zeta_0\left( x\right) ,\zeta_{\kappa}(x,k)\right\rangle _{\mathcal{H}}=2\pi f_{\kappa}\left( \omega\right) \frac{\mathrm{d}\omega}{\mathrm{d}k}, }[/math]

and

[math]\displaystyle{ \left\langle \partial_t\zeta_0\left( x\right) ,\zeta_{\kappa}(x,k)\right\rangle _{\mathcal{H}}=2\pi g_{\kappa}\left( \omega\right) \frac{\mathrm{d}\omega}{\mathrm{d}k}. }[/math]

We can therefore write, changing the variable of integration to [math]\displaystyle{ k }[/math], as

[math]\displaystyle{ \zeta\left( x,t\right) =\frac{1}{2\pi}\int_{\mathbb{R}^{+}} \sum_{\kappa\in\left\{ -1,1\right\}} \Big\{ \left\langle \zeta_0\left( x\right) ,\zeta_{\kappa}(x,k)\right\rangle _{\mathcal{H}} \cos(\omega t) + \left\langle \partial_t \zeta_0\left( x\right) ,\zeta_{\kappa}(x,k)\right\rangle _{\mathcal{H}} \frac{\sin(\omega t)}{\omega}\Big\} \zeta_{\kappa}(x,k) \mathrm{d}k, }[/math]

If we take the case when [math]\displaystyle{ \partial_t\zeta_0( x) =0 }[/math] and write the integral given by the inner product explicitly, we obtain

[math]\displaystyle{ \zeta\left( x,t\right) =\int_{\mathbb{R}^{+}}\Big\{ \sum_{\kappa\in\left\{ -1,1\right\} }\left( \frac{1}{2\pi}\int_{F}\zeta_0\left( x^{\prime}\right) \zeta_{\kappa}(x^{\prime},k) ^{*}\,\mathrm{d} x^{\prime}\right) \zeta_{\kappa}(x,k)\Big\}\cos(\omega t)\mathrm{d}k. }[/math]

Case of a trapped mode

If the scattering structure supports a trapped mode at a particular frequency [math]\displaystyle{ \omega_{0} }[/math] then the expression for the free-surface elevation becomes

[math]\displaystyle{ \zeta\left( x,t\right) =\int_{\mathbb{R}^{+}}\Big\{ \sum_{\kappa\in\left\{ -1,1\right\} }\left( \frac{1}{2\pi}\int_{F}\zeta_0\left( x^{\prime}\right) \zeta_{\kappa}(x^{\prime},k) ^{*}\,\mathrm{d} x^{\prime}\right) \zeta_{\kappa}(x,k)\Big\}\cos(\omega t)\mathrm{d}k + \left(\frac{\int_{F}\zeta_{0}(x^{\prime})\tilde\zeta(x^{\prime})^{*}\mathrm{d}x^{\prime} } {\int_{F}\tilde\zeta(x^{\prime})\tilde\zeta(x^{\prime})^{*}\mathrm{d}x^{\prime} }\right)\tilde\zeta(x)\cos(\omega_{0} t) }[/math] [math]\displaystyle{ + \int_{\mathbb{R}^{+}}\Big\{ \sum_{\kappa\in\left\{ -1,1\right\} }\left( \frac{1}{2\pi}\int_{F}\partial_t\zeta_0\left( x^{\prime}\right) \zeta_{\kappa}(x^{\prime},k) ^{*}\,\mathrm{d} x^{\prime}\right) \zeta_{\kappa}(x,k)\Big\}\frac{\sin(\omega t)}{\omega} \mathrm{d}k + \left(\frac{\int_{F} \partial_t\zeta_{0}(x^{\prime})\tilde\zeta(x^{\prime})^{*}\mathrm{d}x^{\prime} }{\int_{F} \tilde\zeta(x^{\prime})\tilde\zeta(x^{\prime})^{*}\mathrm{d}x^{\prime} }\right)\tilde\zeta(x)\frac{\sin(\omega_{0} t)}{\omega_0} }[/math]

where [math]\displaystyle{ \tilde\zeta(x) }[/math] is the trapped mode free-surface elevation. This formula can be extend to the case when there is more than one trapped mode.

An identity linking waves from the left and right

A consequence of the requirement that the displacement be real, if the initial displacement and initial derivative of displacement is real, is that

[math]\displaystyle{ \sum_{\kappa\in\left\{ -1,1\right\}} \left\langle \zeta_0\left( x\right) ,\zeta_{\kappa}(x,k)\right\rangle _{\mathcal{H}} \zeta_{\kappa}(x,k), }[/math]

must be purely real. This can only be true if

[math]\displaystyle{ \Im \left\{ \zeta_{1}(x^\prime,k)^{*} \zeta_{1}(x,k) \right\} = - \Im \left\{ \zeta_{-1}(x^\prime,k)^{*} \zeta_{-1}(x,k) \right\}, \,\,\,x,x^{\prime} \in F. }[/math]

Matlab Code

• A program to solve two small docks gem_two_small_dock.m

with output shown on the right

The evolution of an initial surface displacement, with two docks as shown

Additional code

This program requires

• Calculate Fourier integral k_integrate.m
• A program to calculate the coefficients for the finite dock problems can be found here two_finite_docks_symmetry.m
• dispersion_free_surface.m
• gif_add_frame.m