Floating Elastic Plate on Variable Bottom Topography

From WikiWaves
Jump to navigationJump to search


Introduction

We present here the solution of Wang and Meylan 2004. The solution method for the linear wave forcing of a two dimensional floating plate on water of variable depth is derived from Liu and Liggett 1982, Hazard and Lenoir 1993 and Meylan and Squire 1994. The method is based on dividing the water domain into two semi-infinite domains and a finite domain. The finite domain contains both the plate and the region of variable water depth. Laplace's equation in the finite domain is solved by the boundary element method. Laplace's equation in the semi-infinite domains is solved by separation of variables. The solution in the semi-infinite domains gives an integral equation relating the normal derivative of the potential and the potential on the boundary of the finite and semi-infinite domains. The thin plate equation is expressed as an integral equation relating the normal derivative of the potential and the potential under the surface of the plate. The boundary element equations and the integral equations are solved simultaneously using the same discretisation.

Problem Formulation

We consider a thin plate, floating on the water surface above a sea bed of variable depth, which is subject to an incoming wave. The plate is approximated as infinitely long in the [math]y[/math] directions which reduces the problem to two dimensions, [math]x[/math] and [math]z[/math]. The [math]x[/math]-axis is horizontal and the [math] z [/math]-axis points vertically up with the free water surface at [math]z=0.[/math] The incoming wave is assumed to be travelling in the positive [math]x[/math]-direction with a single radian frequency [math]\omega [/math]. The theory which will be developed could be extended to oblique incident waves using the standard method Namba and Ohkusu 1999. However, to keep the treatment straightforward, this will not be done. We assume that the wave amplitude is sufficiently small that the problem can be approximated as linear. From the linearity and the single frequency wave assumption it follows that all quantities can be written as the real part of a complex quantity whose time dependence is [math]e^{-i\omega t}[/math].

The linear boundary value problem for the water is the following,

[math] \left. \begin{matrix} \nabla ^{2}\phi =0, \\ -\rho g\,w+i\omega \rho \phi =p,\qquad z=0, \\ \phi _{n}=0,\qquad z=d\left( x\right) , \end{matrix} \right\} (1) [/math]

where [math]\rho [/math] is the density of the water, [math]p[/math] is the pressure on water surface, [math]w[/math] is the displacement of the water surface, [math]g[/math] is the gravitational acceleration, [math]d(x)[/math] is the water depth and [math]\phi _{n}[/math] is the outward normal derivative of the potential. We assume that the water depth is constant outside a finite region, [math]-l\lt x\lt l[/math], but allow the depth to be different at either end. Therefore

[math] d\left( x\right) =\left\{ \begin{matrix} -H_{1},\;\;x\lt -l, \\ d\left( x\right) ,\;\;-l\lt x\lt l, \\ -H_{2},\;\;x\gt l, \end{matrix} \right. [/math]

where [math]H_{1}[/math] and [math]H_{2}[/math] are the water depths in the left and right hand domains of constant depth.

A thin plate, of negligible draft, floats on the surface of the water and occupies the region [math]-L\leq x\leq L\ [/math]as is shown in Figure \ref{fig_region}. Without loss of generality, we assume that [math]l[/math] is sufficiently large that [math] L\lt l.[/math] For any point on the water surface not under the plate the pressure is the constant atmospheric pressure whose time-dependent part is zero. Under the plate, the pressure and the displacement are related by the Bernoulli-Euler equation

[math] D\frac{\partial ^{4}w}{\partial x^{4}}-\rho ^{\prime }a\omega ^{2}\,w=p,\qquad -L\leq x\leq L\;\;\; \mathrm{and}\;\;\; z=0, (2) [/math]

where [math]\rho ^{\prime }[/math] is the density of the plate, [math]a[/math] is the plate thickness, and [math]D[/math] is the bending rigidity of the plate. We assume that the plate edges are free so the bending moment and shear must vanish at both ends of the plate, i.e.

[math] \frac{\partial ^{2}w}{\partial x^{2}}=\frac{\partial ^{3}w}{\partial x^{3}} =0,\;\;\; \mathrm{at}\;\;\; x=-L \;\;\; \mathrm{and}\;\;\;x=L. [/math]

The kinematic boundary condition at the surface allows us to express the displacement as a function of the outward normal derivative of the potential,

[math] w=\frac{i\phi _{n}}{\omega },\qquad z=0. (3) [/math]

Substituting equations (2) and (3)\ into equation (1), we obtain the following boundary value problem for the potential only,

[math] \left. \begin{matrix} \nabla ^{2}\phi =0, \\ -\rho \left( g\phi _{n}-\omega ^{2}\,\phi \right) =\left\{ \begin{matrix} 0,\;x\notin \left[ -L,L\right] ,\;z=0, \\ D\,\frac{\partial ^{4}\phi _{n}}{\partial x^{4}}-\rho ^{\prime }a\,\omega ^{2}\phi _{n},\;x\in \left[ -L,L\right] ,\;z=0, \end{matrix} \right. \\ \phi _{n}=0,\qquad z=d\left( x\right) , \end{matrix} \right\} (4) [/math]

together with the free plate edge conditions

[math] \frac{\partial ^{2}\phi _{n}}{\partial x^{2}}=\frac{\partial ^{3}\phi _{n}}{ \partial x^{3}}=0,\;\;\; \mathrm{at}\;\;\; x=\pm L, \ \ = z=0. (5) [/math]

Radiation Boundary Conditions

Equation (4) is subject to radiation conditions as [math] x\rightarrow \pm \infty .[/math] We assume that there is a wave incident from the left which gives rise to a reflected and transmitted wave. Therefore the following boundary conditions for the potential apply as [math]x\rightarrow \pm \infty [/math]

[math] \lim_{x\rightarrow -\infty }\phi \left( x,z\right) =\cosh \left( k_{t}^{\left( 1\right) }\left( z+H_{1}\right) \right) e^{ik_{t}^{\left( 1\right) }x}+R\,\cosh \left( k_{t}^{\left( 1\right) }\left( z+H_{1}\right) \right) e^{-k_{t}^{\left( 1\right) }x}, (6) [/math]

and

[math] \lim_{x\rightarrow \infty }\phi \left( x,z\right) =T\cosh \left( k_{t}^{\left( 2\right) }\left( z+H_{2}\right) \right) e^{k_{t}^{\left( 2\right) }x}, (7) [/math]

where [math]R[/math] and [math]T[/math] are the reflection and transmission coefficients respectively and [math]k_{t}^{\left( j\right) }\,\left( j=1,2\right) [/math] are the positive real solutions of the Dispersion Relation for a Free Surface

[math] gk_{t}^{\left( j\right) }\tanh \left( k_{t}^{\left( j\right) }\,H_{j}\right) =\omega ^{2}. [/math]

Non-dimensionalisation

Non-dimensional variables are now introduced. We non-dimensionalise the space variables with respect to the water depth on the left hand side, [math] H_{1},[/math] and the time variables with respect to [math]\;\sqrt{g/H_{1}}[/math]. The non-dimensional variables, denoted by an overbar, are

[math] \bar{x}=\frac{x}{H_{1}},\;\bar{z}=\frac{z}{H_{1}},\;\bar{t}=t\sqrt{\frac{g }{H_{1}}},\;\mathrm{and} \;\bar{\phi}=\frac{1}{H_{1}\sqrt{H_{1}g}}\,\phi . [/math]

Applying this non-dimensionalisation to equation (4) we obtain

[math] \left. \begin{matrix}{c} \nabla ^{2}\bar{\phi}=0, \\ \left( \bar{\phi}_{n}-\nu \bar{\phi}\right) =\left\{ \begin{matrix}{c} 0,\qquad \bar{x}\notin \left[ -L,L\right] ,\;\bar{z}=0, \\ -\beta \,\frac{\partial ^{4}\bar{\phi}_{n}}{\partial \bar{x}^{4}}+\gamma \nu \,\bar{\phi}_{n},\qquad \bar{x}\in \left[ -L,L\right] ,\;\bar{z}=0, \end{matrix} \right. \\ \bar{\phi}_{n}=0,\qquad \bar{z}=\bar{d}\left( \bar{x}\right) , \end{matrix} \right\} (9) [/math]

where

[math] \beta =\frac{D}{\rho gH_{1}^{4}},\;\gamma =\frac{\rho ^{\prime }a}{\rho H_{1}},\;\mathrm{and} \;\nu =\frac{\omega ^{2}H_{1}}{g}. (10) [/math]

We will refer to [math]\beta [/math] as the stiffness, [math]\gamma [/math] as the mass and [math]\nu [/math] as the wavenumber. The non-dimensional water depth is

[math] \bar{d}\left( \bar{x}\right) =\left\{ \begin{matrix}{c} -1,\;\;\bar{x}\lt -\bar{l}, \\ \bar{d}\left( \bar{x}\right) ,\;\;-\bar{l}\lt \bar{x}\lt \bar{l}, \\ -H,\;\;\bar{x}\gt \bar{l}, \end{matrix} \right. [/math]

where [math]H=H_{2}/H_{1}.[/math] Equation (9) is also subject to the non-dimensional free edge conditions of the plate (5) and the radiation conditions (6) and (7). With the understanding that all variables have been non-dimensionalised, from now onwards we omit the overbar.

Reduction to a Finite Domain

We solve equation (9) by reducing the problem to a finite domain which contains both the region of variable depth and the floating thin plate. This finite domain [math]\Omega =\left\{ -l\leq x\leq l,\;d\left( x\right) \leq z\leq 0\right\} [/math] is shown in Figure \ref{fig_region}. We will solve Laplace's equation in [math]\Omega [/math] using the boundary element method. To accomplish this we need to express the normal derivative of the potential on the boundary of [math]\Omega \;(\partial \Omega )[/math] as a function of the potential on the boundary.

Green's Function Solution for the Floating Thin Plate

We begin with the boundary condition under the plate which is the following

[math] \beta \frac{\partial ^{4}\phi _{n}}{\partial x^{4}}-\left( \gamma \nu -1\right) \phi _{n}=\nu \phi ,\;-L\leq x\leq L\;\;\; \mathrm{and}\;\;\; z=0, (11) [/math]

together with the free edge boundary conditions (5). Following Meylan and Squire 1994 we can transform equations (5) and (11) to an integral equation using the Green function, [math]g,[/math] which satisfies

[math] \left. \begin{matrix}{c} \beta \frac{\partial ^{4}}{\partial x^{4}}g\left( x,\xi \right) -\left( \gamma \nu -1\right) g\left( x,\xi \right) =\nu \delta \left( x-\xi \right) , \\ \frac{\partial ^{2}}{\partial x^{2}}g\left( x,\xi \right) =\frac{\partial ^{3}}{\partial x^{3}}g\left( x,\xi \right) =0,\;\;\; \mathrm{at}\;\;\; x=-L,\;x=L. \end{matrix} \right\} [/math]

This gives us the following expression for [math]\phi _{n}[/math] as a function of the potential under the plate [math]\phi ,[/math]

[math] \phi _{n}\left( x\right) =\int_{-L}^{L}g\left( x,\xi \right) \;\phi \left( \xi \right) \;d\xi . (12) [/math]

We will write this in operator notation as [math]\phi _{n}=\mathbf{g}\phi [/math] where [math]\mathbf{g}[/math] denotes the integral operator with kernel [math]g\left( x,\xi \right) .[/math]

Solution in the Semi-infinite Domains

We now solve Laplace's equation in the semi-infinite domains [math]\Omega ^{-}=\left\{ x\lt -l,\;-1\leq z\leq 0\right\} [/math] and [math]\Omega ^{+}=\left\{ x\gt l,\;-H\leq z\leq 0\right\} [/math] which are shown in Figure (\ref{fig_region}). Since the water depth is constant in these regions we can solve Laplace's equation by separation of variables. The potential in the region [math]\Omega ^{-} [/math] satisfies the following equation

[math] \left. \begin{matrix} \nabla ^{2}\phi =0,\;\;\mathbf{x}\in \Omega ^{-}, \\ \phi _{n}-\nu \phi =0,\;\;z=0, \\ \phi _{n}=0,\;\;z=-1, \\ \phi =\tilde{\phi}\left( z\right) ,\;\;x=-\,l, \\ \lim\limits_{x\rightarrow -\infty }\phi \left( x,z\right) =\cosh \left( k_{t}^{\left( 1\right) }\left( z+1\right) \right) e^{ik_{t}^{\left( 1\right) }x} \\ +R\,\cosh \left( k_{t}^{\left( 1\right) }\left( z+1\right) \right) e^{-ik_{t}^{\left( 1\right) }x}, \end{matrix} \right\} (13) [/math]

where [math]\mathbf{x}=\left( x,z\right) \ [/math]and [math]\tilde{\phi}\left( z\right) [/math] is an arbitrary continuous function[math].[/math]\ Our aim is to find the outward normal derivative of the potential on [math]x=-l[/math] as a function of [math]\tilde{\phi}\left( z\right) [/math].

We solve equation (13) by separation of variables \cite{Liu82, Hazard} and obtain the following expression for the potential in the region [math] \Omega ^{-},[/math]

[math]\begin{matrix} \phi \left( x,z\right) &=&\cosh \left( k_{t}^{\left( 1\right) }\left( z+1\right) \right) \,e^{ik_{t}^{\left( 1\right) }x}+R\cosh \left( k_{t}^{\left( 1\right) }\left( z+1\right) \right) \,e^{-ik_{t}^{\left( 1\right) }x} \\ &&+\sum_{m=1}^{\infty }\left\langle \tilde{\phi}\left( z\right) ,\tau _{m}^{\left( 1\right) }\left( z\right) \right\rangle \tau _{m}^{\left( 1\right) }\left( z\right) e^{k_{m}^{\left( 1\right) }\left( x+l\right) }. (14) \end{matrix}[/math]

The functions [math]\tau _{m}^{\left( 1\right) }\left( z\right) [/math] ([math]m\geq 1)[/math] are the orthonormal modes given by

[math] \tau _{m}^{\left( 1\right) }\left( z\right) =\left( \frac{1}{2}+\frac{\sin \left( 2k_{m}^{\left( 1\right) }\,\right) }{4k_{m}^{\left( 1\right) }} \right) ^{-\frac{1}{2}}\cos \left( k_{m}^{\left( 1\right) }\left( z+1\right) \right) ,\;\;m\geq 1. [/math]

The evanescent eigenvalues [math]k_{m}^{\left( 1\right) }[/math] are the positive real solutions of the Dispersion Relation for a Free Surface

[math] -k_{m}^{\left( 1\right) }\,\tan \left( k_{m}^{\left( 1\right) }H_{j}\right) =\nu ,\;\;m\geq 1, (15) [/math]

ordered by increasing size. The inner product in equation (13) is the natural inner product for the region [math]-1\leq z\leq 0[/math] given by

[math] \left\langle \tilde{\phi}\left( z\right) ,\tau _{m}^{\left( j\right) }\left( z\right) \right\rangle =\int_{-1}^{0}\tilde{\phi}\left( z\right) ,\tau _{m}^{\left( j\right) }\left( z\right) dx. (16) [/math]

The reflection coefficient is determined by taking an inner product of equation (14) with respect to [math]\cosh \left( k_{t}^{\left( 1\right) }\left( z+1\right) \right) .[/math] This gives us the following expression for [math]R[/math],

[math] R=\frac{\left\langle \tilde{\phi}\left( z\right) ,\cosh \left( k_{t}^{\left( 1\right) }\left( z+1\right) \right) \right\rangle }{\frac{1}{2 }+\sinh \left( 2k_{t}^{\left( 1\right) }\,\right) /4k_{0}^{\left( 1\right) }} e^{-ik_{t}^{\left( 1\right) }l}-e^{-2ik_{t}^{\left( 1\right) }l}. (17) [/math]

The normal derivative of the potential on the boundary of [math]\Omega ^{-}[/math] and [math] \Omega [/math] [math]\left( x=-l\right) [/math] is calculated using equation (14) and we obtain,

[math] \left. \phi _{n}\right| _{x=-l}=\mathbf{Q}_{1}\tilde{\phi}\left( z\right) -2ik_{t}^{\left( 1\right) }\cosh \left( k_{t}^{\left( 1\right) }\left( z+1\right) \right) \,e^{-ik_{t}^{\left( 1\right) }l}, [/math]

where the outward normal is with respect to the [math]\Omega [/math] domain. The integral operator [math]\mathbf{Q}_{1}[/math] is given by

[math]\begin{matrix} \mathbf{Q}_{1}\tilde{\phi}\left( z\right) &=&\sum_{m=1}^{\infty }k_{m}^{\left( 1\right) }\left\langle \tilde{\phi}\left( z\right) ,\tau _{m}^{\left( 1\right) }\left( z\right) \right\rangle \,\tau _{m}^{\left( 1\right) }\left( z\right) (18) \\ &&+ik_{t}\frac{\left\langle \tilde{\phi}\left( z\right) ,\cosh \left( k_{t}^{\left( 1\right) }\left( z+1\right) \right) \right\rangle \cosh \left( k_{t}^{\left( 1\right) }\left( z+1\right) \right) }{\frac{1}{2}+\sinh \left( 2k_{t}^{\left( 1\right) }\,\right) /4k_{t}^{\left( 1\right) }}. \end{matrix}[/math]

We can combine the two terms of equation (18) and express [math]\mathbf{Q}_{1}[/math] as

[math] \mathbf{Q}_{1}\tilde{\phi}\left( z\right) =\sum_{m=0}^{\infty }k_{m}^{\left( 1\right) }\left\langle \tilde{\phi}\left( z\right) ,\tau _{m}^{\left( 1\right) }\left( z\right) \right\rangle \,\tau _{m}^{\left( 1\right) }\left( z\right) (19) [/math]

where [math]k_{0}^{\left( 1\right) }=ik_{t}^{\left( 1\right) }[/math] and

[math] \tau _{0}^{\left( 1\right) }\left( z\right) =\left( \frac{1}{2}+\frac{\sin \left( 2k_{0}^{\left( 1\right) }\right) }{4k_{0}^{\left( 1\right) }}\right) ^{-\frac{1}{2}}\cos \left( k_{0}^{\left( 1\right) }\left( z+1\right) \right) . [/math]

As well as providing a more compact notation, equation (19) will be useful in the numerical calculation of [math]\mathbf{Q}_{1}.[/math]

Similarly, we now consider the potential in the region [math]\Omega ^{+}[/math] which satisfies

[math] \left. \begin{matrix} \nabla ^{2}\phi =0,\;\;\mathbf{x}\in \Omega ^{+}, \\ \phi _{n}-\nu \phi =0,\;\;z=0, \\ \phi _{n}=0, \;\;z=-H, \\ \phi =\tilde{\phi}\left( z\right) ,\;\;x=\,l, \\ \lim\limits_{x\rightarrow -\infty }\phi \left( x,z\right) =T\cosh \left( k_{t}^{\left( 2\right) }\left( z+H_{2}\right) \right) e^{ik_{t}^{\left( 2\right) }x}. \end{matrix} \right\} (20) [/math]

Solving equation (20) by separation of variables as before we obtain

[math] \left. \phi _{n}\right| _{x=l}=\mathbf{Q}_{2}\tilde{\phi}\left( z\right) , [/math]

where the outward normal is with respect to the [math]\Omega [/math] domain. The integral operator [math]\mathbf{Q}_{2}[/math] is given by

[math] \mathbf{Q}_{2}\tilde{\phi}\left( z\right) =\sum_{m=0}^{\infty }k_{m}^{\left( 2\right) }\left\langle \tilde{\phi}\left( z\right) ,\tau _{m}^{\left( 2\right) }\left( z\right) \right\rangle \tau _{m}^{\left( 2\right) }\left( z\right) . (21) [/math]

The orthonormal modes [math]\tau _{m}^{\left( 2\right) }[/math] are given by

[math] \tau _{m}^{\left( 2\right) }\left( z\right) =\left( \frac{H}{2}+\frac{\sin \left( 2k_{m}^{\left( 2\right) }\,H\right) }{4k_{m}^{\left( 2\right) }} \right) ^{-\frac{1}{2}}\cos \left( k_{m}^{\left( 2\right) }\left( z+H\right) \right) , [/math]

The eigenvalues [math]k_{m}^{\left( 2\right) }[/math] are the positive real solutions [math] \left( m\geq 1\right) [/math] and positive imaginary solutions [math]\left( m=0\right) [/math] of the dispersion equation

[math] -k_{m}^{\left( 2\right) }\,\tan \left( k_{m}^{\left( 2\right) }H\right) =\nu. [/math]

The inner product is the same as that given by equation (16) except that the integration is from [math]z=-H[/math] to [math]z=0.[/math] The transmission coefficient, [math]T,[/math] is given by

[math] T=\frac{\left\langle \tilde{\phi}\left( z\right) ,\cosh \left( k_{t}^{\left( 2\right) }\left( z+1\right) \right) \right\rangle }{\frac{1}{2 }+\sinh \left( 2k_{t}^{\left( 2\right) }H\right) /4k_{t}^{\left( 2\right) }} e^{-ik_{t}^{\left( 2\right) }l}. (22) [/math]

The equation for the finite domain

We now consider the finite domain [math]\Omega .[/math] In this domain, Laplace's equation is subject to the boundary conditions given by equations (\ref {integral_plate}), (18) and (21) as well as the free surface and sea floor boundary conditions. This gives the following equation for the potential in [math]\Omega ,[/math]

[math] \left. \begin{matrix} \nabla ^{2}\phi =0,\;\;\;\mathbf{x}\in \Omega , \\ \phi _{n}=\mathbf{Q}_{1}\phi -2ik_{t}^{\left( 1\right) }\cosh \left( k_{t}^{\left( 1\right) }\left( z+1\right) \right) \,e^{-ik_{t}^{\left( 1\right) }l},\;\;\;\mathbf{x}\in \partial \Omega _{1}, \\ \phi _{n}=\mathbf{Q}_{2}\phi ,\;\;\;\mathbf{x}\in \partial \Omega _{2}, \\ \phi _{n}=\nu \phi ,\;\;\;\mathbf{x}\in \partial \Omega _{3}\cup \partial \Omega _{5}, \\ \phi _{n}=\mathbf{g}\phi ,\;\;\;\mathbf{x}\in \partial \Omega _{4}, \\ \phi _{n}=0,\;\;\;\mathbf{x}\in \partial \Omega _{6}. \end{matrix} \right\} (23) [/math]

This boundary value problem is shown in Figure \ref{fig_sigma}. The boundary of [math]\Omega [/math] ([math]\partial \Omega )[/math] has been divided into the six boundary regions [math]\partial \Omega _{i}[/math] shown in the figure. They are respectively, the boundary of [math]\Omega ^{-}[/math] and [math]\Omega [/math] ([math]\partial \Omega _{1}),[/math] the boundary of [math]\Omega ^{+}[/math] and [math]\Omega [/math] ([math]\partial \Omega _{2}),[/math] the free water surface to the left ([math]\partial \Omega _{3}[/math]) and right of the plate ([math] \partial \Omega _{5}[/math]), the plate ([math]\partial \Omega _{4}[/math]), and the sea floor ([math]\partial \Omega _{6}[/math]). Equation (23) is the boundary value problem which we will solve numerically.

Numerical Solution Method

We have reduced the problem to Laplace's equation in a finite domain subject to certain boundary conditions (23). These boundary conditions give the outward normal derivative of the potential as a function of the potential but this is not always a point-wise condition; on some boundaries it is given by an integral equation. We must solve both Laplace's equation and the integral equations numerically. We will solve Laplace's equation by the Boundary Element Method and the integral equations by numerical integration. However, the same discretisation of the boundary will be used for both numerical solutions.

We begin by applying the Boundary Element Method. This gives us the following equation relating the potential and its outward normal derivative on the boundary [math]\partial \Omega [/math]

[math] \frac{1}{2}\phi \left( \mathbf{x}\right) =\int_{\partial \Omega }\left( G_{n}\left( \mathbf{x},\mathbf{x}^{\prime }\right) \phi \left( \mathbf{x} ^{\prime }\right) -G\left( \mathbf{x},\mathbf{x}^{\prime }\right) \phi _{n}\left( \mathbf{x}^{\prime }\right) \right) d\mathbf{x}^{\prime },\;\; \mathbf{x}\in \partial \Omega . (24) [/math]

In equation (24) [math]G\left( \mathbf{x},\mathbf{x}^{\prime }\right) [/math] is the free space Green function given by

[math] G\left( \mathbf{x},\mathbf{x}^{\prime }\right) =\frac{1}{2\pi }\ln \,\left| \mathbf{x}-\mathbf{x}^{\prime }\right| , (25) [/math]

and [math]G_{n}\left( \mathbf{x},\mathbf{x}^{\prime }\right) [/math] is the outward normal derivative of [math]G[/math] (with respect to the [math]\mathbf{x}^{\prime }[/math] coordinate).

We solve equation (24) by a modified constant panel method which reduces it to the following matrix equation

[math] \frac{1}{2}\vec{\phi}=\mathbf{G}_{n}\vec{\phi}-\mathbf{G}\vec{\phi}_{n}. (26) [/math]

In equation (26) [math]\vec{\phi}\mathcal{\ }[/math]and [math]\vec{\phi} _{n}[/math] are vectors which approximate the potential and its normal derivative around the boundary [math]\partial \Omega [/math], and [math]\mathbf{G}[/math] and [math]\mathbf{G}_{n}[/math] are matrices corresponding to the Green function and the outward normal derivative of the Green function respectively. The method used to calculate the elements of the matrices [math]\mathbf{G}[/math] and [math]\mathbf{G}_{n}[/math] will be discussed in section \ref{Green}.

The outward normal derivative of the potential, [math]\vec{\phi}_{n},[/math] and the potential, [math]\vec{\phi},[/math] are related by the conditions on the boundary [math] \partial \Omega [/math] in equation (23). This can be expressed as

[math] \vec{\phi}_{n}=\mathbf{A}\,\vec{\phi}-\vec{f}, (27) [/math]

where [math]\mathbf{A}[/math] is the block diagonal matrix given by

[math] \mathbf{A}\mathbb{=}\left[ \begin{matrix} \mathbf{Q}_{1} & & & & & \\ & \mathbf{Q}_{2} & & & & \\ & & \nu \,\mathbf{I} & & & \\ & & & \mathbf{g} & & \\ & & & & \nu \,\mathbf{I} & \\ & & & & & 0 \end{matrix} \right] , (28) [/math]

[math]\mathbf{Q}_{1}[/math], [math]\mathbf{Q}_{2}[/math], and [math]\mathbf{g}[/math] are matrix approximations of the integral operators of the same name and [math]\vec{f}[/math] is the vector

[math] \vec{f}=\left[ \begin{matrix} 2ik_{t}^{\left( 1\right) }\cosh \left( k_{t}^{\left( 1\right) }\left( z+1\right) \right) \,e^{ik_{t}^{\left( 1\right) }l} \\ 0 \\ \vdots \\ 0 \end{matrix} \right] . (29) [/math]

The methods used to construct the matrices [math]\mathbf{Q}_{1},\mathbf{Q}_{2}[/math] , and [math]\mathbf{g}[/math] will be described in sections 31 and \ref {numericalg} respectively.

Substituting equation (27) into equation (\ref {panelEqn_boundary}) we obtain the following matrix equation for the potential

[math] \left( \frac{1}{2}-\mathbf{G}_{n}+\mathbf{GA}\right) \vec{\phi}=\mathbf{G} \vec{f} [/math]

which can be solved straightforwardly. The reflection and transmission coefficients are calculated from [math]\vec{\phi}[/math] using equations (\ref {reflection}) and (22) respectively.

Numerical Calculation of [math]\mathbf{G}[/math] and [math]\mathbf{G}_{n}[/math]

The boundary element equation (24) is solved numerically by a modified constant panel method. In this method, the boundary is divided into panels over which the potential, [math]\phi ,[/math] or its outward normal derivative, [math]\phi _{n},[/math] are assumed to be constant. The free-space Green's function, [math]G,[/math] and its normal derivative, [math]G_{n}[/math] are more rapidly varying and have a singularity at [math]\mathbf{x}=\mathbf{x}^{\prime }[/math]. For this reason, over each panel, while [math]\phi [/math] and [math]\phi _{n}[/math] are assumed constant, [math]G[/math] and [math]G_{n}[/math] are integrated exactly. For example, we use the following approximation to calculate the integral of [math]G[/math] and [math]\phi [/math] over a single panel

[math] \int_{\mathbf{x}_{i}-h/2}^{\mathbf{x}_{i}+h/2}G\left( \mathbf{x},\mathbf{x} ^{\prime }\right) \phi \left( \mathbf{x}^{\prime }\right) d\mathbf{x} ^{\prime }\approx \phi \left( \mathbf{x}_{i}\right) \int_{\mathbf{x} _{i}-h/2}^{\mathbf{x}_{i}+h/2}G\left( \mathbf{x},\mathbf{x}^{\prime }\right) d\mathbf{x}^{\prime }, (30) [/math]

where [math]\mathbf{x}_{i}[/math] is the midpoint of the panel and [math]h[/math] is the panel length. The integral on the right hand side of equation (30), because of the simple structure of [math]G[/math], can be calculated exactly.

Numerical Calculation of [math]\mathbf{Q}_{1}[/math] and [math]\mathbf{Q}_{2} [/math]

We will discuss the numerical approximation of the operator [math]\mathbf{Q}_{1}[/math]. The operator [math]\mathbf{Q}_{2}[/math] is approximated in a similar fashion. We begin with equation (19) truncated to a finite number ([math]N[/math]) of evanescent modes,

[math] \mathbf{Q}_{1}\phi =\sum_{m=0}^{N}k_{m}^{\left( 1\right) }\left\langle \phi \left( z\right) ,\tau _{m}^{\left( 1\right) }\left( z\right) \right\rangle \tau _{m}^{\left( 1\right) }\left( z\right) . [/math]

We calculate the inner product

[math] \left\langle \phi \left( z\right) ,\tau _{m}^{\left( 1\right) }\left( z\right) \right\rangle =\int\nolimits_{-1}^{0}\phi \left( z\right) \,\tau _{m}^{\left( 1\right) }\left( z\right) \,dz [/math]

with the same panels as we used to approximate the integrals of the Green function and its normal derivative in subsection \ref{Green}. Similarly, we assume that [math]\phi [/math] is constant over each panel and integrate [math]\tau _{m}^{\left( 1\right) }\left( z\right) [/math] exactly. This gives us the following matrix factorisation of [math]\mathbf{Q}_{1},[/math]

[math] \mathbf{Q}_{1}\,\vec{\phi}=\mathbf{S}\,\mathbf{R}\,\vec{\phi}. [/math]

The components of the matrices [math]\mathbf{S}[/math] and [math]\mathbf{R}[/math] are

[math]\begin{matrix} s_{im} &=&\tau _{m}^{\left( 1\right) }\left( z_{i}\right) , \\ r_{mj} &=&k_{m}^{\left( 1\right) }\int_{z_{j}-h/2}^{z_{j}+h/2}\tau _{m}^{\left( 1\right) }\left( s\right) ds \end{matrix}[/math]

where [math]z_{j}[/math] is the value of the [math]z[/math] coordinate in the centre of the panel and [math]h[/math] is the panel length. The integral operator [math]\mathbf{Q}_{2}[/math] is approximated by a matrix in exactly the same manner.

Numerical Calculation of [math]\mathbf{g}[/math]

The method used to approximate [math]\mathbf{g}[/math] by a matrix is similar to the methods used for [math]\mathbf{Q}_{1}[/math] and [math]\mathbf{Q}_{2}[/math] and follows \cite {jgrfloe1d} with some modification. The Green function for the plate can be expressed as

[math] g\left( x,\xi \right) =\left\{ \begin{matrix}{c} A_{1}e^{i\alpha x}+B_{1}e^{-i\alpha x}+C_{1}e^{\alpha x}+D_{1}e^{-\alpha x},\qquad x\lt \xi , \\ A_{2}e^{i\alpha x}+B_{2}e^{-i\alpha x}+C_{2}e^{\alpha x}+D_{2}e^{-\alpha x},\qquad x\gt \xi , \end{matrix} \right. [/math]

where the coefficients are determined by solving a linear system \cite {jgrfloe1d}. Again the same panels are used to approximate the integral operator by a matrix as were used for the boundary integral equation (\ref {integral_eqn}). Over each panel we assume that the potential is constant and integrate the Green function [math]g[/math] exactly.

Results

We will now present some results, concentrating on comparing the reflection coefficient for constant and variable depth profiles. This will allow us to determine when the variable depth profile has a significant effect. To reduce the number of figures we restrict ourselves to four values of the stiffness [math]\beta [/math] and two variable depth profiles.

Profiles for the variable depth

We will consider two different profiles for the variable depth. The first will be the profile which was used by Staziker, Porter and Stirling 1996. This corresponds to a rise from a uniform depth to a maximum height of half the uniform depth at [math]x=0[/math]. The formula for this profile is the following

[math] d\left( x\right) =\left\{ \begin{matrix}{c} -1,\;\;x\lt -l, \\ -\left( \frac{1}{2}\left( \frac{x+l}{l}\right) ^{2}-\frac{x+l}{l} +1\right) ,\;\;-l\lt x\lt l, \\ -1,\;\;x\gt l. \end{matrix} \right. (33) [/math]

Following Staziker, Porter and Stirling 1996 we will refer to this variable depth profile as the hump.

In the second profile the depth rises linearly. The depth in the right hand region is half the depth in the left hand region. The formula for this profile is

[math] d\left( x\right) =\left\{ \begin{matrix}{c} -1,\;\;x\lt -l, \\ \frac{x+l}{4l}-1,\;\;-l\lt x\lt l, \\ -\frac{1}{2},\;\;x\gt l. \end{matrix} \right. (34) [/math]

We will refer to this variable depth profile as the simple slope .

Convergence study

We now present a convergence study. Since we have two parameters, the panel size used to discretise the boundary and the number of evanescent modes ([math]N[/math] ) used to approximate the integral equations [math]\mathbf{Q}_{1}[/math] and [math]\mathbf{Q} _{2},[/math] we must present two convergence studies considering each parameter separately. We begin by considering the panel size used to discretise the boundary. We expect that the panel size should be proportional to the wavelength and therefore inversely proportional to the wavenumber [math] k_{t}^{\left( 1\right) }[/math] (assuming that the water depth at either end is of a similar size). We therefore use the following formula for the panel size

[math] \mathrm{panel size} =\frac{1}{\kappa \,k_{t}^{\left( 1\right) }}, [/math]

where [math]\kappa [/math] is a constant of proportionality which will be determined from the convergence study. Table 35 shows the absolute value of the reflection coefficient for a plate of length [math]L=2.5,[/math] stiffness [math] \beta =1,[/math] and mass [math]\gamma =0[/math] for [math]\nu =1,[/math] 2, and 3 for a constant depth and for the hump with [math]l=2.5.[/math] The number of evanescent modes, [math]N,[/math] was fixed to be 5. Three values of [math]\kappa [/math] were considered, [math]\kappa =10,[/math] 20, and 40. The results in Table 35 show that good convergence is achieved when [math]\kappa =20.[/math]

Table 37 shows a similar convergence study for the number of evanescent modes. We have considered [math]0,[/math] 5 and [math]10[/math] evanescent modes ([math]N[/math]) and fixed [math]\kappa [/math] to be [math]\kappa =20[/math]. All other parameters as the same as those in Table 35. The results in Table 37 show that good convergence is achieved when the number of modes is 5. For all subsequent calculations the panel size will be determined by setting [math]\kappa =20[/math] and [math]N=5.[/math]

Comparison with existing results

Before presenting our results for a plate on water of variable depth we will make comparisons with two results from the literature. This is to establish that our method gives the correct solution for the simpler cases of either variable depth but no plate, or a plate floating on constant depth. We begin by comparing our results with Staziker, Porter and Stirling 1996 who solved for wave scattering by variable depth only. One problem which they solved was to determine the absolute value of the reflection coefficient for a hump depth profile with fixed frequency [math]\nu =1[/math] and variable hump length [math]l.[/math] The solution to this problem by our method is shown in Figure (\ref{staziker1.ps} ). This figure is identical to Figure 2 in Staziker, Porter and Stirling 1996 (p. 290) which establishes that our method gives the correct solution for a variable depth in the absence of the plate.

The second comparison is with Meylan and Squire 1994 in which the problem of a thin plate on water of constant depth was solved (to model an ice floe). One problem which they solved was the absolute value of the reflection coefficient as a function of plate length for fixed [math]\nu .[/math] The dimensional parameters which they used were, density [math]\rho ^{\prime }=922.5\,\mathrm{kgm}^{-3},[/math] thickness [math]h=1[/math]m, and bending rigidity [math]D=[/math]5.4945[math]\times 10^{8}\mathrm{kgm}^{2}\mathrm{s}^{-2}.[/math] The water density was 1025[math]\mathrm{kgm} ^{-3}[/math] and the incoming wave was chosen to have wavelength [math]100[/math]m. The solution to this problem by our method is shown in Figure (\ref{fig_mike_refl}). This figure is identical to Figure 3 in Meylan and Squire 1994 (p. 895) which establishes that our method gives the correct result for a plate on water of constant depth.

Reflection

We will consider the absolute value of the reflection coefficient as a function of wavenumber [math]\nu [/math] for various values of the parameters. Figure ( \ref{plot4hump}) shows the absolute value of the reflection coefficient as a function of [math]\nu [/math] with [math]\beta =0.01,[/math] [math]0.1,[/math] 1 and 10 and [math]\gamma =0[/math] for the hump depth profile. Both the length of the hump and the length of the plate were fixed to be [math]l=L=2.5[/math]. The solution for the plate and hump (solid line), plate with constant depth (dashed line) and hump only (dotted line) are shown. The two simpler solutions are drawn so that the full solution may be compared to these simpler cases. Figure (\ref{plot4hump}) shows the existence of two asymptotic regimes. When [math]\nu [/math] is small (low frequency or large wavelength)\ the reflection is dominated by the hump and the plate is transparent. For large [math]\nu [/math] the reflection is dominated by the plate and the hump is transparent. As the value of the stiffness [math]\beta [/math] is increased the hump dominated region becomes smaller and the plate dominated region becomes larger. This not unexpected because increasing the value of [math]\beta [/math] increases the influence of the plate. It is apparent that, especially for smaller values of stiffness, there is a large region where the hump and plate solution is significantly different from the plate only solution, even though the hump only reflection is practically zero. This is because the wavelength under the plate is larger than the open water wavelength so the depth variation is felt more strongly when the plate is present.

Figure (\ref{plot4plane}) is equivalent to Figure (\ref{plot4hump}) except that the depth profile is the simple slope and the results are also very similar[math].[/math] Figure (\ref{plot4hump_gamma}) is equivalent to Figure (\ref {plot4hump}) except that the value of stiffness is fixed ([math]\beta =0.1)[/math] and the value of [math]\gamma [/math] is varied. This figure shows that for realistic (small) values of [math]\gamma [/math] this parameter is not significant. This explains why [math]\gamma [/math] is often neglected (e.g. Namba and Ohkusu 1999) and why we have chosen [math]\gamma =0[/math] for Figures (\ref{plot4hump}) and (\ref{plot4plane}).

Figure (\ref{plot4ahump}) is equivalent to figure (\ref{plot4hump}) except that the hump has been moved by [math]L[/math] to the left so that the minimum depth is directly underneath the left (incoming) end of the plate. Comparing figure ( \ref{plot4ahump}) with figure (\ref{plot4ahump}) we see that moving the hump has made a significant change to the reflection coefficient for low frequencies, especially as the stiffness [math]\beta [/math] is increased.

Displacements

Finally we investigate the displacement of the plate for some of the regimes we have considered. We present the displacement for the variable and constant depth profiles so that we may compare the effect of the variable depth. We divide the displacement by [math]i\omega /k_{t}^{\left( 1\right) }\sinh \left( k_{t}^{\left( 1\right) }\right) [/math] so that the incoming wave now has unit amplitude in displacement at the water surface. Figure (\ref {deflbeta4hump_nu05}) shows the displacement of the plate for [math]\nu =0.5[/math] and [math]\beta =0.01,[/math] [math]0.1,[/math] 1 and 10 and [math]\gamma =0.[/math] The plate length is [math]L=2.5[/math]. The solid line is the real part of the displacement and the dotted line is the imaginary part of the displacement. The thicker line is the solution with the hump depth profile ([math]l=2.5)[/math] and the thinner line is the solution with a constant depth. It is apparent from this figure, that the bending of the plate is increased by the presence of the hump, but that this effect is not very strong. Figures (\ref{deflbeta4hump_nu1}) and (\ref {deflbeta4hump_nu2}) are identical to Figure (\ref{deflbeta4hump_nu05}) except that [math]v=1[/math] and [math]\nu =2[/math] respectively. These figures also show only a slight increase in the bending of the plate due to the hump. It appears that, while the variable depth does have a significant effect on the reflection coefficient, the effect on the plate displacement is not necessarily as strong.


Tables

\begin{table}[h] \centering

(35) \begin{tabular}{lll} [math]\nu =1[/math] & & \\ \hline \multicolumn{1}{|c}{[math]\kappa [/math]} & \multicolumn{1}{|c}{plate only} & \multicolumn{1}{|c|}{plate and hump} \\ \hline \multicolumn{1}{|c}{10} & \multicolumn{1}{|c}{0.2983746166} & \multicolumn{1}{|c|}{0.2491046427} \\ \hline \multicolumn{1}{|c}{20} & \multicolumn{1}{|c}{0.2957175612} & \multicolumn{1}{|c|}{0.2470511349} \\ \hline \multicolumn{1}{|c}{40} & \multicolumn{1}{|c}{0.2948480461} & \multicolumn{1}{|c|}{0.2465182527} \\ \hline [math]\nu =2[/math] & & \\ \hline \multicolumn{1}{|c}{[math]\kappa [/math]} & \multicolumn{1}{|c}{plate only} & \multicolumn{1}{|c|}{plate and hump} \\ \hline \multicolumn{1}{|c}{10} & \multicolumn{1}{|c}{0.3457203891} & \multicolumn{1}{|c|}{0.1934227806} \\ \hline \multicolumn{1}{|c}{20} & \multicolumn{1}{|c}{0.3461627544} & \multicolumn{1}{|c|}{0.1947144005} \\ \hline \multicolumn{1}{|c}{40} & \multicolumn{1}{|c}{0.3463681703} & \multicolumn{1}{|c|}{0.1952205755} \\ \hline [math]\nu =3[/math] & & \\ \hline \multicolumn{1}{|c}{[math]\kappa [/math]} & \multicolumn{1}{|c}{plate only} & \multicolumn{1}{|c|}{plate and hump} \\ \hline \multicolumn{1}{|c}{10} & \multicolumn{1}{|c}{0.0277944414} & \multicolumn{1}{|c|}{0.2605833203} \\ \hline \multicolumn{1}{|c}{20} & \multicolumn{1}{|c}{0.0249319083} & \multicolumn{1}{|c|}{0.2568361963} \\ \hline \multicolumn{1}{|c}{40} & \multicolumn{1}{|c}{0.0236686063} & \multicolumn{1}{|c|}{0.2550926942} \\ \hline \end{tabular} \caption{[math]\left|R\right|[/math] for [math]\kappa[/math] = 10, 20, and 40 and

[math]\nu [/math] = 1,2, and 3 for the plate only and the plate and hump.
[math]\beta=1[/math], [math]\gamma =0[/math] and 5 evanscent modes are 

used. (36)}


\end{table}

\pagebreak


\begin{table}[t] \centering

(37) \begin{tabular}{ccc} [math]\nu =1[/math] & & \\ \hline \multicolumn{1}{|c}{[math]N[/math]} & \multicolumn{1}{|c}{plate only} & \multicolumn{1}{|c|}{plate and hump} \\ \hline \multicolumn{1}{|c}{0} & \multicolumn{1}{|c}{0.2934754434} & \multicolumn{1}{|c|}{0.2448842047} \\ \hline \multicolumn{1}{|c}{5} & \multicolumn{1}{|c}{0.2957175612} & \multicolumn{1}{|c|}{0.2470511349} \\ \hline \multicolumn{1}{|c}{10} & \multicolumn{1}{|c}{0.2957697025} & \multicolumn{1}{|c|}{0.2470814786} \\ \hline [math]\nu =2[/math] & & \\ \hline \multicolumn{1}{|c}{[math]N[/math]} & \multicolumn{1}{|c}{plate only} & \multicolumn{1}{|c|}{plate and hump} \\ \hline \multicolumn{1}{|c}{0} & \multicolumn{1}{|c}{0.3392954189} & \multicolumn{1}{|c|}{0.2021691049} \\ \hline \multicolumn{1}{|c}{5} & \multicolumn{1}{|c}{0.3461627544} & \multicolumn{1}{|c|}{0.1947144005} \\ \hline \multicolumn{1}{|c}{10} & \multicolumn{1}{|c}{0.3459741914} & \multicolumn{1}{|c|}{0.1944525347} \\ \hline [math]\nu =3[/math] & & \\ \hline \multicolumn{1}{|c}{[math]N[/math]} & \multicolumn{1}{|c}{plate only} & \multicolumn{1}{|c|}{plate and hump} \\ \hline \multicolumn{1}{|c}{0} & \multicolumn{1}{|c}{0.0262065602} & \multicolumn{1}{|c|}{0.2296877448} \\ \hline \multicolumn{1}{|c}{5} & \multicolumn{1}{|c}{0.0249319083} & \multicolumn{1}{|c|}{0.2568361963} \\ \hline \multicolumn{1}{|c}{10} & \multicolumn{1}{|c}{0.0259483771} & \multicolumn{1}{|c|}{0.2578144528} \\ \hline \end{tabular} \caption{[math]\left|R\right|[/math] for 0, 5, and 10 evanescent modes ([math]N[/math]) and [math]\nu [/math] = 1,2, and 3 for the plate only and the plate and hump.

[math]\beta=1[/math], [math]\gamma =0[/math] and 

[math]\kappa = 20[/math].(38)}


\end{table}