1. Introduction
Hydroelasticity, which is concerned with the interaction of fluids with elastic bodies, is an important topic with numerous applications, including those in the cryosphere (ocean wave scattering by sea ice and ice shelves) and those related to the engineering of very large floating structures. A comprehensive account of the historical research can be found in the review articles of Squire (Reference Squire2007, Reference Squire2020). We give here a brief account of some of the hydroelastic models. We note that this extensive research has focussed almost entirely on isotropic plates, with a notable exception being the recent work of Thery et al. (Reference Thery, Peter, Bennetts and Guenneau2026), who used a space transformation approach to determine the parameters of an anisotropic plate mounted over shallow water for cloaking applications.
The two-dimensional water wave problem (one horizontal and one vertical) was solved by an eigenfunction expansion method by Fox & Squire (Reference Fox and Squire1994) and extended by Sahoo, Yip & Chwang (Reference Sahoo, Yip and Chwang2001), who analysed a semi-infinite elastic plate floating on the surface of water of finite depth. The Wiener–Hopf technique is particularly suited to semifinite problems and highly analytic solutions have been derived in this case (Balmforth & Craster Reference Balmforth and Craster1999; Chung & Fox Reference Chung and Fox2002; Tkacheva Reference Tkacheva2003; Williams & Meylan Reference Williams and Meylan2012; Smith et al. Reference Smith, Peter, Abrahams and Meylan2020). The finite elastic plate was solved by Meylan & Squire (Reference Meylan and Squire1994) using the Green’s function technique, which led to a Fredholm integral equation and was simultaneously solved using a modal expansion by Newman (Reference Newman1994). The three-dimensional water wave problem involving circular elastic discs was solved by Meylan & Squire (Reference Meylan and Squire1996) and Peter, Meylan & Chung (Reference Peter, Meylan and Chung2004). The solution for plates of arbitrary shapes was given in Meylan (Reference Meylan2001, Reference Meylan2002). Similar methods were developed by Hermans (Reference Hermans2004)and Andrianov & Hermans (Reference Andrianov and Hermans2006).
Recently, extensive studies on hydroelasticity have been driven by its potential applications in wave energy conversion, with flexible wave energy converters (WECs), continuing to be a focal point of current research (Babarit et al. Reference Babarit, Singh, Mélis, Wattez and Jean2017; Philen et al. Reference Philen, Squibb, Groo and Hagerman2018; Tay & Wei Reference Tay and Wei2020; Collins et al. Reference Collins, Hossain, Dettmer and Masters2021; Renzi et al. Reference Renzi, Michele, Zheng, Jin and Greaves2021; Michele, Zheng & Greaves Reference Michele, Zheng and Greaves2022; Valappil & Koley Reference Valappil and Koley2022; Zheng et al. Reference Zheng, Michele, Liang, Meylan and Greaves2022; Guo, Mohapatra & Soares Reference Guo, Mohapatra and Soares2023; Singh & Gayen Reference Singh and Gayen2023; Cheng et al. Reference Cheng, Du, Dai, Yuan and Incecik2024; Meylan et al. Reference Meylan, Challis, Thamwattana, Wegert and Wilks2025). A promising mechanism for flexible wave energy conversion is the piezoelectric effect, in which materials become electrically polarised in response to elastic deformation. First comprehensively modelled by Renzi (Reference Renzi2016), a piezoelectric WEC (usually modelled as a plate) would become polarised in response to deformation by ocean waves and thus induce a current in an external circuit. An important aspect of piezoelectric materials that has been hitherto ignored in the piezoelectric WEC literature is their anisotropy – while piezoceramics are usually transversely isotropic (Yang Reference Yang2018), the polymer Polyvinylidene fluoride, which is often proposed for piezoelectric WECs, is anisotropic (Vinogradov & Holloway Reference Vinogradov and Holloway1999). To date, authors have bypassed considerations of anisotropy by focussing on the one-dimensional WEC/two-dimensional fluid context, with submerged horizontal plates (Renzi Reference Renzi2016; Valappil & Koley Reference Valappil and Koley2022), surface-mounted plates (Meylan et al. Reference Meylan, Challis, Thamwattana, Wegert and Wilks2025) and vertical plates (Singh & Gayen Reference Singh and Gayen2023) being among the cases considered. The present paper seeks to investigate wave scattering by a two-dimensional anisotropic elastic plate in a three-dimensional fluid. In doing so, this paper will facilitate future extensions of earlier models of piezoelectric WECs (Renzi Reference Renzi2016; Meylan et al. Reference Meylan, Challis, Thamwattana, Wegert and Wilks2025) to three dimensions.
The outline of this paper is as follows. In § 2 we describe a model of a rectangular anisotropic elastic plate mounted on the water surface. Our method is an extension of that presented by Meylan et al. (Reference Meylan, Challis, Thamwattana, Wegert and Wilks2025) for the two-dimensional case, in which the solution is expanded as a sum of a diffraction potential and a superposition of radiation potentials. The diffraction potential, which assumes that the plate is fixed, is solved in § 3. The modes of vibration of the plate and the corresponding radiation potentials are found in § 4 and § 5, respectively. The coupled fluid/plate problem is solved in § 6. An energy balance identity used to verify our computations, which is based on the optical theorem, is derived in § 7. A large number of results are presented in § 8, which serve to illustrate the capabilities of the model and act as benchmark solutions for future studies. A conclusion is given in § 9.
2. Problem outline
2.1. Preliminaries
We consider the scattering of water waves by a rectangular anisotropic elastic plate situated at the surface of the water and having negligible draft, meaning that its draft is assumed to be much smaller than the plate’s in-plane dimensions and the wavelength. We adopt a three-dimensional Cartesian coordinate system
$(x,y,z)$
for the fluid domain, where
$x$
and
$y$
are horizontal coordinates parallel to the plate edges, and the
$z$
coordinate points vertically upward (i.e. it is normal to the plate at rest). At equilibrium, the fluid occupies the region
$\varOmega = {\mathbb{R}}^2\times (-H,0)$
, where
$H$
is the depth of the fluid, and the plate occupies the region
${\varGamma }\times \{0\}$
, where
${\varGamma }=(0,a)\times (0,b)$
. A schematic of the geometry is shown in figure 1.
(a) Side and (b) plan views of the scattering problem. The rectangular plate, labelled
$\varGamma$
, has side lengths
$a$
and
$b$
and the fluid is of depth
$H$
. The incident wave
$\phi ^{{inc}}$
excites the plate into motion, generating scattered waves
$\phi ^{{sc}}$
.
The fluid is assumed to be incompressible, inviscid and undergoing irrotational flow, meaning that the time-dependent velocity potential
$\varPhi$
satisfies
\begin{equation} \bigtriangleup \varPhi =0\quad \forall (x,y,z)\in \varOmega , \end{equation}
where
$\varPhi =\varPhi (x,y,z,t)$
and
$\bigtriangleup =\partial _x^2+\partial _y^2+\partial _z^2$
. The velocity of the fluid is given by
$\boldsymbol{v}=\boldsymbol{\nabla }\varPhi$
. There is no flow normal to the seabed at
$z=-H$
, which implies
\begin{equation} \partial _z\varPhi =0\quad z=-H, (x,y)\in {\mathbb{R}}^2. \end{equation}
The free surface of the fluid occupies the region
$(x,y)\in \mathbb{R}^2\setminus {\varGamma }$
. The free surface of the fluid is governed by the kinematic and dynamic boundary conditions, which are
\begin{align} \partial _t\zeta &= \partial _z\varPhi {\biggr \rvert }_{z=0}, \end{align}
\begin{align} \zeta &=-\frac {1}{g}\partial _t\varPhi {\biggr \rvert }_{z=0}, \end{align}
respectively, where
$\zeta =\zeta (x,y,t)$
is the displacement of the free surface of the fluid from equilibrium and
$g$
is the acceleration due to gravity.
The plate, which has uniform thickness
$h$
, is assumed to be thin, which means
$a,b\gg h$
. The motion of the plate is modelled using Kirchhoff–Love theory. We assume that in-plane deformations can be neglected because these will be negligible compared with bending. The Lagrangian of the plate is
\begin{equation} \mathcal{L}\{W\}=T\{W\}-U_{\textit{strain}}\{W\}-U_{\textit{pressure}}\{W\} \end{equation},
where
$W(x,y,t)$
is the vertical displacement of the plate. The kinetic energy, strain energy due to bending and potential energy due to pressure are
\begin{align} T\{W\}&=\frac {1}{2}\rho h\iint _\varGamma (\partial _t W)^2\mathrm{d} x\,\mathrm{d} y, \end{align}
\begin{align} U_{\textit{strain}}\{W\}&=\frac {1}{2}\iint _\varGamma \mathcal{B}(W,W)\mathrm{d} x\,\mathrm{d} y, \end{align}
\begin{align} U_{\textit{pressure}}\{W\}&=-\iint _\varGamma PW\mathrm{d} x\,\mathrm{d} y, \end{align}
respectively, where
$\rho$
is the uniform density of the plate and
$P(x,y,t)$
is the fluid pressure under the plate (Grossi & Nallim Reference Grossi and Nallim2003; Reddy Reference Reddy2006; Grossi & Nallim Reference Grossi and Nallim2008). The symmetric bilinear operator
$\mathcal{B}$
is given by
\begin{equation} \mathcal{B}(W_1,W_2)=\begin{bmatrix} \partial _x^2 W_1&\partial _y^2 W_1&\partial _x\partial _y W_1 \end{bmatrix}\begin{bmatrix} D_{11}&D_{12}&2D_{16}\\ D_{12}&D_{22}&2D_{26}\\ 2D_{16}&2D_{26}&4D_{66} \end{bmatrix}\begin{bmatrix} \partial _x^2W_2\\\partial _y^2 W_2\\\partial _x\partial _y W_2 \end{bmatrix}, \end{equation}
where
$D_{ij}$
are the rigidity coefficients which describe the material’s resistance to bending and twisting (Reddy Reference Reddy2006). The action over an arbitrary time interval
$(t_0,t_1)$
and its first variation are
\begin{align} \mathcal{A} &=\int _{t_0}^{t_1}\mathcal{L}\{W\}\mathrm{d} t, \end{align}
\begin{align} \delta \mathcal{A}&=\int _{t_0}^{t_1}\delta \mathcal{L}\{\delta W,W\}\mathrm{d} t \end{align}
\begin{align} &=\int _{t_0}^{t_1}\int _\varGamma \left \{\rho h\delta W \partial _t^2 W+\mathcal{B}(\delta W,W)-P\delta W \right \}\mathrm{d} x\,\mathrm{d} y\,\mathrm{d} t, \end{align}
respectively, where
$\delta W$
is an admissible variation of
$W$
. The principle of least action implies that
$\delta \mathcal{A}=0$
for arbitrary admissible variations
$\delta W$
over arbitrary time intervals
$(t_0,t_1)$
, which yields
\begin{equation} \iint _\varGamma \left \{\rho h\delta W\partial _t^2 W+\mathcal{B}(\delta W,W)-P\delta W\right \}\,\mathrm{d} x\,\mathrm{d} y=0. \end{equation}
The strong form of (2.10) is
\begin{equation} \rho h\partial _t^2{W}+\mathcal{D}W=P, \end{equation}
where the fourth-order differential operator
$\mathcal{D}$
is given by
\begin{equation} \mathcal{D}=D_{11}\partial _{x}^4+4D_{16}\partial _{x}^3\partial _{y}+2(D_{12}+2D_{66})\partial _{x}^2\partial _{y}^2+4D_{26}\partial _{x}\partial _{y}^3+D_{22}\partial _{y}^4. \end{equation}
Note that, in the isotropic case, we have
\begin{align} D_{11}=D_{22}&=D, \end{align}
\begin{align} D_{12}&=\nu D, \end{align}
\begin{align} D_{66}&=\dfrac {1}{2}(1-\nu )D, \end{align}
\begin{align} D_{16}=D_{26}&=0, \end{align}
\begin{align} \text{where }D&=\frac {Eh^3}{12(1-\nu ^2)}, \end{align}
in which
$E$
is Poisson’s ratio and
$\nu$
is the Young’s modulus. In this case, the operator
$\mathcal{D}$
is proportional to the biharmonic operator.
From the linearised Bernoulli equation, the fluid pressure is the sum of hydrostatic and hydrodynamic pressure, namely
$P=P_{\textit{hydrostatic}}+P_{\textit{hydrodynamic}}$
, where
\begin{equation} P_{\textit{hydrostatic}}=-\rho _wgz\quad \text{and}\quad P_{\textit{hydrodynamic}}=-\rho _w \partial _t\varPhi , \end{equation}
respectively, in which
$\rho _w$
is the fluid density. By assuming small motions, we linearise the hydrodynamic pressure about the equilibrium position of the plate at
$z=0$
, which leads to the following dynamic condition on the plate:
\begin{equation} \rho h \partial _t^2{W}+\mathcal{D}W=-\rho _wgW-\rho _w\partial _t\varPhi {\biggr \rvert }_{z=0}. \end{equation}
The coupling between the plate and the fluid is completed using the following linearised kinematic condition:
\begin{equation} \partial _z\varPhi {\biggr \rvert }_{z=0}=\partial _t{W}, \end{equation}
for
$(x,y)\in {\varGamma }$
. By equating the fluid velocity at the fluid–plate interface with the plate velocity, this condition assumes that the plate remains fully in contact with the fluid at all times. We note that linear hydroelasticity is well validated experimentally for sufficiently low ratios of amplitude to wavelength (Montiel et al. Reference Montiel, Bennetts, Squire, Bonnefoy and Ferrant2013a
,
Reference Montiel, Bonnefoy, Ferrant, Bennetts, Squire and Marsaultb
).
In this paper, the plate edges are assumed to be either clamped, implying
\begin{align} W=\partial _x W &= 0&x=0\text{ or }x=a, \end{align}
\begin{align} W=\partial _y W&=0&y=0\text{ or }y=b, \end{align}
or free, implying
\begin{align} D_{11}\partial _x^2 W +2D_{16}\partial _x\partial _y W+D_{12}\partial _y^2 W&=0&x=0\text{ or }x=a, \end{align}
\begin{align} D_{11}\partial _x^3 W+4 D_{16}\partial _x^2\partial _y W+(D_{12}+4D_{66})\partial _x\partial _y^2W+2D_{26}\partial _y^3 W&=0&x=0\text{ or }x=a, \end{align}
\begin{align} D_{22}\partial _y^2 W+2D_{26}\partial _x\partial _y W+D_{12}\partial _x^2 W&=0&y=0\text{ or }y=b, \end{align}
\begin{align} D_{22}\partial _y^3 W+4D_{26}\partial _y^2\partial _x W+(D_{12}+4D_{66})\partial _y\partial _x^2W+2D_{16}\partial _y^3 W&=0&y=0\text{ or }y=b, \end{align}
\begin{align} D_{16}\partial _x^2W+D_{26}\partial _y^2W+2D_{66}\partial _{x}\partial _yW&=0&(x,y)\in \varLambda , \end{align}
where
$\varLambda =\{(0,0),(a,0),(0,b),(a,b)\}$
are the corners of the plate (Grossi & Nallim Reference Grossi and Nallim2003), or simply supported, implying (§ 6.1.2 of Reddy Reference Reddy2006)
\begin{align} W=D_{11}\partial _x^2 W+D_{12}\partial _y^2 W &= 0&x=0\text{ or }x=a, \end{align}
\begin{align} W=D_{12}\partial _x^2 W+D_{22}\partial _y^2 W&=0&y=0\text{ or }y=b. \end{align}
We note that, while our framework can be extended to other boundary condition types (e.g. moored edges) or have different boundary conditions on different edges, we do not explore this here.
2.2. Time-harmonic formulation
We assume a time-harmonic motion with angular frequency
$\omega$
, which implies that we can write
\begin{align} \varPhi (x,y,z,t)&={\mathrm{Re}}\{\phi (x,y,z){\mathrm{e}}^{-\mathrm{i}\omega t}\}&(x,y,z)\in \varOmega , \end{align}
\begin{align} \zeta (x,y,t)&={\mathrm{Re}}\{\eta (x,y){\mathrm{e}}^{-\mathrm{i}\omega t}\}&(x,y)\in {\mathbb{R}}^2\setminus {\varGamma }, \end{align}
\begin{align} W(x,y,t)&={\mathrm{Re}}\{{w}(x,y){\mathrm{e}}^{-\mathrm{i}\omega t}\}&(x,y)\in {\varGamma }. \\[9pt] \nonumber \end{align}
This leads to the following boundary value problem for
$\phi$
:
\begin{align} \bigtriangleup \phi &=0&(x,y,z)\in \varOmega , \end{align}
\begin{align} \partial _z\phi &=0&z=-H,\quad (x,y)\in \mathbb{R}^2, \end{align}
\begin{align} \partial _z\phi &=\alpha \phi &z=0,\quad (x,y)\in \mathbb{R}^2\setminus {\varGamma }, \end{align}
\begin{align} \partial _z\phi &=-\mathrm{i}\omega {w}&z=0,\quad (x,y)\in {\varGamma }, \end{align}
\begin{align} \sqrt {r}(\partial _r-\mathrm{i} k)(\phi -\phi ^{{inc}})&\to 0&\text{as }r\to \infty , \end{align}
with
$\alpha = \omega ^2/g$
, which must be solved in tandem with the plate equation
\begin{equation} {\mathcal{D}}{w}-\omega ^2\rho h{w}+\rho _w g{w}=\mathrm{i}\omega \rho _w\phi {\biggr \rvert }_{z=0}, \end{equation}
for
$(x,y)\in {\varGamma }$
, where the plate satisfies either clamped (2.17), free (2.18) or simply supported (2.19) boundary conditions, as prescribed. In the Sommerfeld condition (2.21e),
$r=\sqrt {x^2+y^2}$
,
$k$
is the positive real root of the dispersion relation
\begin{equation} k\tanh (k H)=\alpha , \end{equation}
and
$\phi ^{{inc}}$
is the prescribed incident wave.
3. Diffraction by a rigid plate
We first consider the diffraction potential
$\phi ^{{di}}$
, which describes the scattering of a prescribed incident wave
$\phi ^{{inc}}$
by a plate held fixed in a flat configuration at the fluid surface. This potential satisfies the boundary value problem consisting of (2.21a
), (2.21b
), (2.21c
), (2.21e
), together with the no flow condition on the plate
\begin{equation} \partial _z[{\phi ^{{inc}}}+\phi ^{{di}}]=0\quad \text{for } z=0,\quad (x,y)\in {\varGamma }. \end{equation}
The unknown diffraction potential is obtained by solving a boundary integral equation.
3.1. Green’s function
The Green’s function for a three-dimensional fluid of constant depth with a point source on the surface of the fluid satisfies
\begin{align} \bigtriangleup G &={0}&(x,y,z)\in \varOmega , \end{align}
\begin{align} {\partial _z G- \alpha G} &={-\delta (x)\delta (y)}&z=0, \end{align}
\begin{align} \partial _zG&=0&{z=-H}, \end{align}
\begin{align} \sqrt {r}(\partial _rG-\mathrm{i} kG)&\to 0&\text{as }{r}\to \infty , \end{align}
where
$\delta (x)$
is the Dirac delta. For
$z=0$
, the Green’s function has the following representation:
\begin{equation} -4\pi G(r)=\frac {1}{r}+\frac {1}{\sqrt {r^2+4H^2}}+\mathop {{\smile }}\!\!\!\!\!\!\mathop {\int }_0^\infty \frac {2(\mu +\alpha ){\mathrm{e}}^{-\mu H}\cosh ^2\mu H}{\mu \sinh \mu H-\alpha \cosh \mu H} \mathrm{J}_0(\mu r)\mathrm{d}\mu , \end{equation}
where
$\mathrm{J}_0$
is the zeroth-order Bessel function of the first kind and
$\mathop {\smile }\!\!\!\!\!\mathop {\int }$
denotes complex contour integration beneath any poles on the real line. We note that (3.3) is equivalent to the submerged source Green’s function given by Linton & McIver (Reference Linton and McIver2001, equation B.86), in the limit as the point source approaches the surface. A derivation of the Green’s function using a Hankel transform is given in Appendix A. To evaluate the above integral numerically, we first write
$\mathrm{J}_0(x)= {1}/{2}(\mathrm{H}_0^{(1)}(x)+\mathrm{H}_0^{(2)}(x))$
(where
$\mathrm{H}_0^{(1)}$
and
$\mathrm{H}_0^{(2)}$
are the zeroth-order Hankel functions of the first and second kinds, respectively) and deform the contour of integration using residue calculus. We obtain the following rapidly converging representation of the Green’s function:
\begin{align} -4\pi G(r)&=\frac {1}{r}+\frac {1}{\sqrt {r^2+4H^2}}+2{\mathrm{Re}}\left \{\int _C\!\frac {(\mu +\alpha ){\mathrm{e}}^{-\mu H}\cosh ^2\mu H}{\mu \sinh \mu H-\alpha \cosh \mu H} \mathrm{H}_0^{(1)}(\mu r)\mathrm{d}\mu \right \}\nonumber \\ &\quad +{\frac {2\pi \mathrm{i}}{H}\left (1+\frac {\sinh 2kH}{2kH}\right )^{-1}}\cosh ^2(kH) \mathrm{H}_0^{(1)}(kr), \end{align}
where
$C$
is parametrised by
${\mathrm{e}}^{{i} \pi /4}t$
for
$0\leqslant t\lt \infty$
. Note that, in order to consider infinite depth, we need only replace (3.3) with the corresponding infinite depth Green’s function. The integral in (3.4) is calculated using the MATLAB routine integral.
3.2. Integral equation for diffraction problem
We assume
$\phi ^{{di}}$
can be written as a distribution of sources of the form
\begin{equation} \phi ^{{di}}(\boldsymbol{x})=\iint _{{\varGamma }\times \{0\}}\sigma (\boldsymbol{x}^\prime )G(\boldsymbol{x},\boldsymbol{x}^\prime )\mathrm{d} S_{\boldsymbol{x}^\prime }, \end{equation}
where, in abuse of notation,
$G(\boldsymbol{x},\boldsymbol{x}^\prime )$
denotes the value of the Green’s function at the point
$\boldsymbol{x}=(x,y,z)$
for a source located at
$\boldsymbol{x}^\prime =(x^\prime ,y^\prime ,z^\prime )$
, whenever its arguments are bold. Then for
$\boldsymbol{x}\in {\varGamma }\times \{0\}$
we compute
\begin{align} (\alpha -\partial _z)\phi ^{{di}}(\boldsymbol{x})&=\iint _{{\varGamma }\times \{0\}}\sigma (\boldsymbol{x}^\prime )(\alpha -\partial _z)G(\boldsymbol{x},\boldsymbol{x}^\prime )\mathrm{d} S_{\boldsymbol{x}^\prime }\nonumber \\ &=-\iint _{{\varGamma }\times \{0\}}\sigma (\boldsymbol{x}^\prime )\delta (x-x^\prime )\delta (y-y^\prime )\mathrm{d} S_{\boldsymbol{x}^\prime }\nonumber \\ &=-\sigma (\boldsymbol{x}). \end{align}
Thus
\begin{equation} \phi ^{{di}}(\boldsymbol{x})=\iint _{{\varGamma }\times \{0\}}(\alpha \phi ^{{di}}(\boldsymbol{x^\prime })-\partial _{z^\prime }\phi ^{{di}}(\boldsymbol{x^\prime }))G(\boldsymbol{x},\boldsymbol{x^\prime })\mathrm{d} S_{\boldsymbol{x^\prime }}. \end{equation}
An integral equation for
$\phi ^{{di}}$
is obtained from the above by restricting to the plate and using (3.1), namely
\begin{align} \phi ^{{di}}(\boldsymbol{x})=\alpha \iint _{{\varGamma }\times \{0\}}G(\boldsymbol{x}^\prime ,\boldsymbol{x})(\phi ^{{ inc}}(\boldsymbol{x}^\prime )+ \phi ^{{di}}(\boldsymbol{x}^\prime ))\mathrm{d} S_{\boldsymbol{x}^\prime }\quad \text{for all}\quad \boldsymbol{x}\in {\varGamma }\times \{0\}. \end{align}
We note that the same integral equation may be obtained by using the Green’s function for a submerged source (as given by Linton & McIver Reference Linton and McIver2001, equation B.86) and applying Green’s third identity to
$\phi ^{{di}}$
. This involves a subtle point which arises when the source point in on the free surface and we replace the normal derivative of the Green function using the free-surface condition. This point is discussed in Meylan & Squire (Reference Meylan and Squire1996) and Meylan (Reference Meylan2002), where a similar Green’s function method is used, and also in Buchner (Reference Buchner1993).
3.3. Numerical solution using a constant panel method
The boundary integral equation above is solved using a constant panel method. To do so, we discretise the plate surface
${\varGamma }\times \{0\}$
into a grid of
$N_x\times N_y$
squares
$\square _j$
. Each square has dimensions
$\Delta x\times \Delta x$
and midpoints
$\boldsymbol{x}_j=(x_j,y_j,0)^\intercal$
. Evaluated at nodes
$\boldsymbol{x}_i$
, (3.8) becomes
\begin{align} \phi ^{{di}}({\boldsymbol{x}_i})&=\alpha \sum _{j=1}^{N_x\times N_y}{\iint _{\square _j}} G({\boldsymbol{x}^\prime },{\boldsymbol{x}_i})\left ({\phi ^{{ inc}}}({\boldsymbol{x}^\prime })+\phi ^{{di}}({\boldsymbol{x}^\prime })\right ){\mathrm{d} S_{\boldsymbol{x}^\prime }} \end{align}
\begin{align} &\approx \alpha \sum _{j=1}^{N_x\times N_y}\left ({\phi ^{{ inc}}}({\boldsymbol{x}_j})+\phi ^{{di}}({\boldsymbol{x}_j})\right ){\iint _{\square _j}} G({\boldsymbol{x}^\prime },{\boldsymbol{x}_i}){\mathrm{d} S_{\boldsymbol{x}^\prime }}, \end{align}
where (3.9b) is the constant panel assumption, valid when
$\phi ^{{inc}}$
and
$\phi ^{{di}}$
are slowly varying relative to the grid spacing. Letting
$(\boldsymbol{\phi }^{{di}})_j=\phi ^{{di}}(x_j,y_j,0)$
and
$(\boldsymbol{\phi }^{{in}})_j={\phi ^{{inc}}}(x_j,y_j,0)$
, (3.9b) may be recast in matrix form as
\begin{equation} (I-\alpha K)(\boldsymbol{\phi }^{{in}}+\boldsymbol{\phi }^{{di}})=\boldsymbol{\phi }^{{in}}, \end{equation}
to be solved for
$\boldsymbol{\phi }^{{di}}$
, where
$I$
is the
$(N_x\times N_y)$
-dimensional identity matrix. The entries in the kernel matrix
$K$
are chosen to approximate the integrals over the squares
$\square _j$
appearing in (3.9b
), namely
\begin{equation} (K)_{ij}\approx {\iint _{\square _j} G(\boldsymbol{x}^\prime ,\boldsymbol{x}_i)\mathrm{d} S_{\boldsymbol{x}^\prime }}. \end{equation}
We approximate the off diagonal entries of
$K$
as
\begin{equation} {(K)_{ij} = G(\boldsymbol{x_j},\boldsymbol{x}_i)\Delta x^2,} \end{equation}
for
$i\neq j$
. This is equivalent to midpoint quadrature on the grid. Greater care is applied to the integrals associated with the diagonal entries of
$K$
, namely
\begin{align} {(K)_{ii}}&{=\iint _{\square _j} G(\boldsymbol{x}^\prime ,\boldsymbol{x}_i)\mathrm{d} S_{\boldsymbol{x}^\prime }} \end{align}
\begin{align} &=\int _{-\Delta x/2}^{\Delta x/2}\int _{-\Delta x/2}^{\Delta x/2} G\left (\sqrt {x^2+y^2},0\right )\mathrm{d} x\, \mathrm{d} y, \end{align}
for all
$i$
(meaning this quantity only needs to be calculated once), where the second line follows from translation. The double integral is then calculated by (i) directly integrating the term
$1/r$
in (3.4), and (ii) using the MATLAB routine integral2 to numerically integrate the remaining part. We remark that in preliminary tests of our method, the integral representation of the Green’s function (3.4) converged faster than methods based on an equivalent series representation (i.e. one analogous to equation B.91 in Linton & McIver (Reference Linton and McIver2001)), presumably due to direct integration of the
$1/r$
term that dominates in the near field. Our numerical code for the diffraction problem was validated for the case of a circular plate, using a code developed to solve the problem of wave scattering by a partially submerged truncated circular cylinder, which is valid in the limit of negligible submergence (Wilks et al. Reference Wilks, Meylan, Montiel, Balasooriya, Jauhar, Wheeler and Chalup2025). The aforementioned code uses separation of variables in cylindrical coordinates to solve the problem both beneath and exterior to the plate. Subsequently, continuity conditions between the interior and exterior regions reduce to integral equations that are solved numerically using a singularity respecting Galerkin method.
4. Dry modes of vibration of the plate
The modes of vibration of the plate are the eigenfunctions
${w}_j$
of the operator
$\mathcal{D}$
, which satisfy
\begin{equation} \mathcal{D}{w}_j=\rho h\omega _j^2 {w}_j, \end{equation}
and the boundary conditions (either (2.17), (2.18) or (2.19) for a clamped, free or simply supported plate, respectively), where
$\omega _j$
are the angular frequencies of the vibrational modes. They are called dry modes because they solve the dynamic equation of the plate (2.11) in the absence of any forcing from the fluid (i.e. with
$P=0$
). The corresponding weak form is
\begin{equation} \iint _{{\varGamma }} \mathcal{B}({w}_j,\delta {w})\mathrm{d} x\, \mathrm{d} y=\rho h\omega _j^2\iint _{{\varGamma }} {w}\delta {w}\mathrm{d} x \mathrm{d} y, \end{equation}
where
$\delta {w}$
is an admissible variation of
${w}_j$
. We seek a numerical solution of (4.2) using the Rayleigh–Ritz method, which is classical but included here for completeness (see Reddy Reference Reddy2006 for further discussion). To begin, we introduce a set of functions
$\{\psi _1,\ldots ,\psi _{N_{RR}}\}$
, which satisfy the boundary conditions (2.17). We approximate
${{w}}_j$
in the form of a superposition over this set, namely
\begin{equation} {w}_j\approx a_1\psi _1+\ldots +a_{N_{RR}}\psi _{N_{RR}}, \end{equation}
where
$a_1,\ldots ,a_{N_{RR}}$
are unknown coefficients. Lastly, we substitute (4.3) into (4.2) for all test functions
$\psi _1,\ldots ,\psi _{N_{RR}}$
, yielding the following system of equations for the unknown coefficients
$a_\theta$
:
\begin{equation} \sum _{\theta =1}^{N_{RR}} a_\theta \iint _{{\varGamma }}\mathcal{B}(\psi _\theta ,\psi _\tau )\mathrm{d} x\,\mathrm{d} y =\rho h \omega _j^2\sum _{\theta =1}^{N_{RR}}a_\theta \iint _{{\varGamma }}\psi _\theta \psi _\tau \,\mathrm{d} x\,\mathrm{d} y. \end{equation}
To proceed, a suitable set of functions must be selected for the Rayleigh–Ritz method. In the case that the plate is clamped, we construct such functions from the modes of vibration of a clamped-clamped Euler–Bernoulli beam of unit length, which, when normalised, satisfy
\begin{align} \partial _x^4 u_m&=\kappa _m^4u_m, \end{align}
\begin{align} u_m(0)&=u_m(1)=0, \end{align}
\begin{align} u_m^\prime (0)&=u_m^\prime (1)=0, \end{align}
\begin{align} \int _0^1 u_m(x)u_n(x)\mathrm{d} x&=\delta _{mn}, \end{align}
where
$\kappa _m$
are the beam eigenvalues and
$\delta _{mn}$
is the Kronecker delta. Further details about these eigenmodes are given in Appendix B. As considered by Young (Reference Young1950), the basis for the Rayleigh–Ritz method is then taken to be
\begin{align} \psi _\theta (x,y)&=u_{m(\theta )}\left (\frac {x}{a}\right ) u_{n(\theta )}\left (\frac {y}{b}\right ), \end{align}
\begin{align} a_\theta &=A_{m(\theta ) n(\theta )}^{(j)}, \end{align}
where
$\theta$
indexes all pairs
$(m,n)\in \mathbb{N}^2$
such that
$1\leqslant m\leqslant N_x$
and
$1\leqslant n\leqslant N_y$
, in which
$N_x$
and
$N_y$
are the number of beam modes used in the
$x$
and
$y$
dimensions, respectively. Thus
$N_{RR}=N_xN_y$
. We note that this choice of basis functions is advantageous since they are mutually orthogonal and satisfy the plate boundary conditions (2.17), although other choices of basis functions could be made (e.g. polynomials, see Reddy Reference Reddy2006, Chapter 7). Arguments of
$\theta$
in indices
$m$
and
$n$
will be dropped in what follows.
With the substitutions given in (4.6), (4.4) eventually becomes
\begin{align} &\sum _{m=1}^{N_x}\sum _{n=1}^{N_y}\left [\frac {D_{12}}{a^2b^2}\left ( \mathcal{U}_{mp}^{[0,2]}\mathcal{U}_{nq}^{[2,0]}+\mathcal{U}_{mp}^{[2,0]}\mathcal{U}_{nq}^{[0,2]}\right )\right .+\frac {2D_{16}}{a^3b}\left (\mathcal{U}_{mp}^{[2,1]}\mathcal{U}_{nq}^{[0,1]}+\mathcal{U}_{mp}^{[1,2]}\mathcal{U}_{nq}^{[1,0]}\right )\nonumber \\&\quad +\frac {2D_{26}}{ab^3}\left (\mathcal{U}_{mp}^{[0,1]}\mathcal{U}_{nq}^{[2,1]}+\mathcal{U}_{mp}^{[1,0]}\mathcal{U}_{nq}^{[1,2]}\right )+\frac {4D_{66}}{a^2b^2}\mathcal{U}_{mp}^{[1,1]}\mathcal{U}_{nq}^{[1,1]}\quad +\frac {D_{11}}{a^4}\mathcal{U}_{mp}^{[2,2]}\mathcal{U}_{nq}^{[0,0]}\nonumber \\&\left .\quad +\frac {D_{22}}{b^4}\mathcal{U}_{mp}^{[0,0]}\mathcal{U}_{nq}^{[2,2]}\right ]A_{mn}^{(j)}=\omega ^2\rho h A_{pq}^{(j)}, \end{align}
for
$1\leqslant p\leqslant N_x$
and
$1\leqslant q\leqslant N_y$
, where we have defined
\begin{equation} \mathcal{U}_{mp}^{[j,l]} = \int _0^1 u_m^{(j)}(x)u_p^{(l)}(x)\mathrm{d}x, \end{equation}
in which
$u_m^{(l)}$
denotes the
$l$
th derivative of
$u_m$
. Equation (4.7) is an eigenvalue problem of dimension
$N_xN_y$
for the coefficients
$A^{(j)}_{mn}$
, which we solve numerically. The eigenmodes
${w}_j(x,y)$
, which are recovered from the coefficients
$A^{(j)}_{mn}$
, are normalised so that
\begin{equation} \iint _{\varGamma } {{w}}_m(x,y){{w}}_n(x,y)^*\,\mathrm{d} x\mathrm{d} y= \delta _{mn}. \end{equation}
In the case that the edges of the plate are free, we construct instead the functions used in the Rayleigh–Ritz method from the modes of vibration of a free–free Euler–Bernoulli beam of unit length, which, when normalised, satisfy
\begin{align} \partial _x^4 u_m&=\kappa _m^4u_m, \end{align}
\begin{align} u_m^{(2)}(0)&=u_m^{(2)}(1)=0, \end{align}
\begin{align} u_m^{(3)}(0)&=u_m^{(3)}(1)=0, \end{align}
\begin{align} \int _0^1 u_m(x)u_n(x)\mathrm{d} x&=\delta _{mn}. \end{align}
Use of these functions instead of the corresponding modes of a clamped–clamped beam yields the modes of the free-edge plate, because free edges are the natural boundary conditions.
Finally, in the case that the edges of the plate are simply supported, we construct the functions used in the Rayleigh–Ritz method from the modes of vibration of a simply supported–simply supported Euler–Bernoulli beam of unit length, which, when normalised, satisfy
\begin{align} \partial _x^4 u_m&=\kappa _m^4u_m, \end{align}
\begin{align} u_m(0)&=u_m(1)=0, \end{align}
\begin{align} u_m^{(2)}(0)&=u_m^{(2)}(1)=0, \end{align}
\begin{align} \int _0^1 u_m(x)u_n(x)\mathrm{d} x&=\delta _{mn}. \end{align}
We take
$N_x/a=N_y/b = 20$
throughout this paper. We have validated our method for computing the dry modes using results provided by Leissa (Reference Leissa1973, table C1) in the case where the plate is isotropic, and results provided in An et al. (Reference An, Ni, Xu and Li2021) in the case where the plate is anisotropic, with excellent agreement – see the first and third columns of figure 2.
First four mode shapes of a square (a,d,g,j) isotropic plate, (b,e,h,k) orthotropic plate and (c,f,i,l) anisotropic plate. The rigidity coefficients associated with these terms are given in table 1. The remaining parameters are
$a=b=1$
m and
$\rho h = 1$
kg m
$^{-1}$
. Note the mode shapes have been rescaled as
${w}_j(x,y)/p_j$
, with
$p_j$
chosen for presentation.
5. Radiation potentials
For each vibration mode of the plate
${w}_j$
, we compute the radiation potential
$\phi ^{{ra}}_j$
, namely the velocity potential of the fluid induced by the plate motion in the form
${w}={\mathrm{Re}}\{{w}_j(x,y){\mathrm{e}}^{-\mathrm{i}\omega t}\}$
. The radiation potential
$\phi ^{{ra}}_j$
satisfies the boundary value problem consisting of (2.21a), (2.21b), (2.21c), (2.21e) in the case
$\phi ^{{inc}}=0$
, and the kinematic condition
\begin{equation} \partial _z\phi ^{{ra}}_j=-\mathrm{i}\omega {w}_j\quad {\text{for }(x,y)\in \varGamma ,\quad z=0}. \end{equation}
Following the steps outlined in § 3.2, we use the Green’s function (introduced in § 3.1) to derive an integral equation for the radiation potential. We obtain
\begin{equation} \phi ^{{ra}}_j(\boldsymbol{x}) = {\iint _{\varGamma \times \{0\}}}G(\boldsymbol{x}^\prime ,\boldsymbol{x})(\alpha \phi ^{{ra}}_j(\boldsymbol{x}^\prime )+\mathrm{i}\omega {w}_j(\boldsymbol{x}^\prime ))\mathrm{d} S_{\boldsymbol{x}^\prime }, \end{equation}
where, in abuse of notation, we have written
${{w}}_j(\boldsymbol{x^\prime })={{w}}_j([x^\prime ,y^\prime ,0]^\intercal )={{w}}_j(x^\prime ,y^\prime )$
. Applying the constant panel discretisation method used in (3.3) then yields the matrix equation
\begin{equation} (I-\alpha K)\left (\boldsymbol{\phi }^{{ra}}_j+\frac {\mathrm{i} g}{\omega }\boldsymbol{w}_j\right ) = \frac {\mathrm{i} g}{\omega }\boldsymbol{w}_j, \end{equation}
to be solved for
$\boldsymbol{\phi }^{{ra}}_j$
, where
$(\boldsymbol{\phi }^{{ra}}_j)_i=\phi ^{{ra}}_j(x_i,y_i)$
and
$(\boldsymbol{w}_j)_i = {w}_j(x_i,y_i)$
.
6. Scattering by an elastic plate
The solution to the coupled fluid–plate system of (2.21) is obtained by expanding in the dry modes of vibration of the plate, following the theory of Newman (Reference Newman1994). That is, the total potential is taken to be a superposition of the incident potential, the diffraction potential plus a weighted sum of the radiation potentials
\begin{equation} \phi = {\phi ^{{inc}}}+\phi ^{{di}}+\sum _{j=1}^\infty c_j \phi ^{{ra}}_j, \end{equation}
where the coefficients
$c_j$
are unknown. By the kinematic condition (2.21d
), the deformation of the plate is of the form
\begin{equation} {w} = \sum _{j=1}^\infty c_j {w}_j. \end{equation}
The dynamic condition (2.21f
), when multiplied by
${w}_m^*$
and integrated over the plate, then yields
\begin{equation} \rho h\omega _m^2c_m-\rho h\omega ^2c_m+\rho _wgc_m = \mathrm{i}\omega \rho _w{\iint _{\varGamma \times \{0\}} }\phi {w}_m^*\mathrm{d} S, \end{equation}
where we have used (4.1). Substituting (6.1) into the above then yields
\begin{align} &\rho h\omega _m^2 c_m -\rho h\omega ^2 c_m +\rho _wg c_m- \mathrm{i}\omega \rho _w\sum _{j=1}^\infty \left ({\iint _{\varGamma \times \{0\}}}\phi ^{{ra}}_j {w}_m^*\mathrm{d} S\right )c_j \nonumber\\&\quad = \mathrm{i}\omega \rho _w{\iint _{\varGamma \times \{0\}}}({\phi ^{{inc}}}+\phi ^{{di}}){w}_m^*\mathrm{d} S. \end{align}
Numerical evaluation of the integrals in the above, and truncation of the series after
$N$
terms, yields an
$N\times N$
system for the unknown coefficients
$c_1,\ldots ,c_N$
. (6.4) becomes
\begin{equation} ({K_{\textit{stiff}}}-\omega ^2M+C-\omega ^2A-\mathrm{i}\omega B)\boldsymbol{c}=\boldsymbol{f}, \end{equation}
where
$K_{\textit{stiff}}$
,
$M$
and
$C$
are the stiffness, mass and hydrostatic restoration matrices of the system, respectively, given by
\begin{align} {K_{\textit{stiff}}} &= \rho h\ {diag}(\omega _m^2), \end{align}
\begin{align} M &= -\rho h I, \end{align}
\begin{align} C &= \rho _w g I, \end{align}
where
$I$
is the
$N\times N$
identity matrix. Moreover, the real-valued added mass and damping matrices
$A$
and
$B$
, respectively, are such that
\begin{equation} -\omega ^2A_{mj}-\mathrm{i}\omega B_{mj}=-\mathrm{i}\omega \rho _w{\iint _{\varGamma \times \{0\}}}\phi ^{{ra}}_j {w}_m^*\mathrm{d} S, \end{equation}
the vector of unknown coefficients is
$\boldsymbol{c}=[c_1,\ldots ,c_N]^\intercal$
and the forcing vector is
\begin{equation} (\boldsymbol{f})_m=\mathrm{i}\omega \rho _w{\iint _{\varGamma \times \{0\}}}({\phi ^{{inc}}}+\phi ^{{di}}){w}_m^*\mathrm{d} S. \end{equation}
Note that in this paper
$N=30$
plate modes were considered, which was considered sufficient based on the energy balance validation in the following section. Also note that our implementation of this method has been validated using code developed by Meylan (Reference Meylan2002) for the case of a square isotropic plate with free edges.
7. Energy balance calculations
To verify our calculations, an energy balance identity is derived for the case in which the plate is excited by an incident plane wave of amplitude
$A$
and incident angle
$\chi$
, namely
\begin{equation} {\phi ^{{inc}}}(x,y,z)=\frac {-\mathrm{i} gA}{\omega }\frac {\cosh (k(z+H))}{\cosh (kH)}{\mathrm{e}}^{{i} k(x\cos \chi +y\sin \chi )}, \end{equation}
such that the corresponding free-surface elevation is
$\eta ^{{in}}(x,y)={\mathrm{e}}^{{i} kx}$
. Let
$(r,\theta ,z)$
be the system of cylindrical coordinates with
$x=r\cos \theta$
and
$y=r\sin \theta$
. For large
$r$
, the potential can be approximated as (Mei, Stiassnie & Yue Reference Mei, Stiassnie and Yue2005)
\begin{equation} \phi \sim \frac {-\mathrm{i} gA}{\omega }\frac {\cosh (k(z+H))}{\cosh (kH)}\left ({\mathrm{e}}^{{i} kr\cos (\theta -\chi )}+\sqrt {\frac {2}{\pi kr}}f(\theta ){\mathrm{e}}^{{i}(k r-\pi /4)}\right ). \end{equation}
In the above, the first term in the brackets describes the incident wave and the second term describes the scattered wave. Thus, the far-field pattern
$f$
must be derived from the scattered field. We first note that the total scattered field is of the form
\begin{equation} \phi -{\phi ^{{inc}}}={\iint _{\varGamma \times \{0\}}}G(\boldsymbol{x}^\prime ,\boldsymbol{x})u(\boldsymbol{x}^\prime )\mathrm{d} S_{\boldsymbol{x}^\prime }, \end{equation}
where, with reference to (3.7) and (5.2),
\begin{equation} u(x,y)=\alpha \left({\phi ^{{inc}}}(x,y,0)+\phi ^{{di}}(x,y,0)\right)+\sum _{j=1}^\infty c_j\left(\alpha \phi _j^{{ra}}(x,y,0)+\mathrm{i}\omega {w}_j(x,y)\right)\!. \end{equation}
Asymptotically, the Green’s function is
\begin{equation} G(r,\theta ,z;0,0)\sim -\frac {\mathrm{i}}{2}\frac {k\cosh (k(z+H))}{k H{\mathrm{sech}}(kH)+\sinh (kH)}\sqrt {\frac {2}{\pi kr}}{\mathrm{e}}^{{i}(kr-\pi /4)}. \end{equation}
Suppose the source point of the Green’s function is shifted to
$(x^\prime ,y^\prime )$
. Let
$d(x^\prime ,y^\prime )$
and
$\alpha (x^\prime ,y^\prime )$
, respectively, be the radial and angular coordinates of the point
$(x^\prime ,y^\prime )$
with respect to the global origin
$(0,0)$
. Then as
$r\to \infty$
, the Green’s function is asymptotically given by (Falnes Reference Falnes2002)
\begin{align} G(r,\theta ,z;x^\prime ,y^\prime )&\sim -\frac {\mathrm{i}}{2}\frac {k\cosh (k(z+H))}{k H{\mathrm{sech}}(kH)+\sinh (kH)}\sqrt {\frac {2}{\pi kr}}\nonumber \\ &\quad \times \exp (\mathrm{i}(kr-kd(x^\prime ,y^\prime )\cos (\alpha (x^\prime ,y^\prime )-\theta )-\pi /4)). \\[9pt] \nonumber \end{align}
With reference to (7.2), substitution of the above expression into (7.3) eventually yields
\begin{equation} f(\theta )=\frac {k\omega }{2g[kH{\mathrm{sech}}^2(kH)+\tanh (kH)]}\iint _{{\varGamma }}u(x^\prime ,y^\prime ){\mathrm{e}}^{-\mathrm{i} kd(x^\prime ,y^\prime )\cos (\alpha (x^\prime ,y^\prime )-\theta )}\mathrm{d} y^\prime \mathrm{d} x^\prime . \end{equation}
Conservation of energy implies that the optical theorem holds, which is (Mei et al. Reference Mei, Stiassnie and Yue2005)
\begin{equation} \frac {1}{2\pi }\int _{-\pi }^\pi |f(\theta )|^2\mathrm{d}\theta =-{{Re}}\{f(\chi )\}. \end{equation}
We have used this identity to verify our computations. Numerical parameters are chosen in order to obtain five figure agreement between the left- and right-hand sides of (7.8). This was obtained for
$\Delta x=(1/120)$
m.
8. Results
Throughout this results section, we take
$H=20$
m and
$\rho h=1$
kg m
$^{-1}$
. We consider three cases of the rigidity coefficients, as outlined in table 1. We take the incident wave to be an incident plane wave of the form (7.1) with
$A=1$
m and
$\chi =0$
(i.e. it is travelling in the positive
$x$
direction) unless otherwise stated. The time averaged kinetic energy of the plate is
\begin{equation} \bar {E}_{{kin}}=\frac {\rho h \omega ^2}{4}\int _0^a\int _0^b|{w}(x,y)|^2\mathrm{d} y \mathrm{d} x =\frac {\rho h\omega ^2}{4}\sum _{j=1}^\infty |c_j|^2. \end{equation}
This quantity is used to identify resonant frequencies of the plate.
Rigidity coefficients used to generate the results in this paper. The Poisson’s ratio associated with the isotropic plate rigidity coefficients is
$0.3$
. The rigidity coefficients of the orthotropic plate correspond to a material under plane-stress conditions for which
$E_y/E_x = 3/4$
,
$G_{xy}/E_x=1/2$
(where
$E_x$
and
$E_y$
are the Young’s moduli in the
$x$
and
$y$
directions, respectively, and
$G_{xy}$
is the shear modulus) and Poisson’s ratio is
$0.3$
. The coefficients for the anisotropic plate are as in An et al. (Reference An, Ni, Xu and Li2021). All units are in
${Pa}\,{m}^{{3}}$
.
Figure 3 shows the kinetic energy spectrum for a square isotropic plate. The two peaks correspond to resonances of the plate, which we investigate in figures 4 and 5. We observe that the first resonance is predominantly associated with excitation of the first mode of vibration of the plate, whereas the second resonance is predominantly associated with excitation of the second and third modes in equal amounts. This is unsurprising because the second and third modes are degenerate, as shown in figure 2. We note that, in both cases, the fourth mode of the plate is not excited (i.e.
$c_4=0$
) because it is antisymmetric with respect to the line
$y=b/2$
, while the incident wave is symmetric with respect to this line. We remark that the far-field patterns (figures 4
b and 4
d) are symmetric with respect to
$\theta =0$
, and the surface elevations (figure 5) are symmetric about
$y=b/2$
.
Kinetic energy spectrum for a square isotropic plate with clamped edges, with
$a=b=1$
m.
Dry mode expansion coefficient magnitudes (left column) and far-field patterns (i.e. polar plots of
$f(\theta )$
, right column) at resonant frequencies of the square isotropic plate with clamped edges and with
$a=b=1$
m. In particular, the resonant frequencies are (a,b)
$\omega =6.42$
s
$^{-1}$
and (c,d)
$\omega =10.08$
s
$^{-1}$
.
Surface elevation of the excited square isotropic plate with clamped edges (
$a=b=1$
m) at the resonant frequencies considered in figure 4, namely (a)
$\omega =6.42$
s
$^{-1}$
and (b)
$\omega =10.08$
s
$^{-1}$
. The surface elevation refers to
$|w|$
over the plate, and
$|\eta |$
otherwise. The absolute value of the surface elevation is indicated by the colour scale, with white isophasic contours indicating where the surface elevation is real. The boundary of the plate is marked with a thick white line.
Kinetic energy spectrum for a square orthotropic plate with clamped edges, with
$a=b=1$
m.
Dry mode expansion coefficient magnitudes (left column) and far-field patterns (right column) at resonant frequencies of the square orthotropic plate with clamped edges and with
$a=b=1$
m. The resonant frequencies are (a,b)
$\omega =6.38$
s
$^{-1}$
and (c,d)
$\omega =10.08$
s
$^{-1}$
.
Next, we consider a square orthotropic plate. The kinetic energy spectrum in figure 6 again has two peaks, which we investigate in figures 7 and 8. The first resonance is predominantly associated with excitation of the first plate mode, whereas the second resonance is predominantly associated with excitation of the third plate mode. In contrast with the square isotropic plate considered earlier, the second and third modes of the orthotropic plate are not degenerate because symmetry has been broken. In particular, the second mode of the orthotropic plate is antisymmetric with respect to the line
$y=b/2$
and hence cannot be excited by the symmetric incident plane wave. This also remains the case with the fourth mode. Thus we observe
$c_2=c_4=0$
in figure 7. We remark that the far-field patterns (figures 7
b and 7
d) are symmetric with respect to
$\theta =0$
, and the surface elevations (figure 8) are symmetric about
$y=b/2$
.
Surface elevation of the excited square orthotropic plate with clamped edges (
$a=b=1$
m) at the resonant frequencies considered in figure 7, namely (a)
$\omega =6.38$
s
$^{-1}$
and (b)
$\omega =10.08$
s
$^{-1}$
.
The next case we consider is a square anisotropic plate, as parametrised by An et al. (Reference An, Ni, Xu and Li2021). These parameters, which are reproduced in table 1, describe a graphite/epoxy laminated plate in which the tensor describing the material is rotated by 45
$^\circ$
in the
$xy$
-plane. As such, the modes of vibration in figure 2 are symmetric or antisymmetric with respect to the lines
$y=x$
and
$x+y=a$
. The response of the surface-mounted square anisotropic plate to a plane wave with incident angle
$\chi =0$
is considered in figures 9–11. We observe three resonant peaks in the kinetic energy spectrum (figure 9), which are predominantly associated with excitation of the first three modes of vibration in order of increasing frequency. The second and third modes of vibration are not degenerate, and neither are antisymmetric with respect to the line
$y=b/2$
, so both can be excited. We remark that the far-field patterns (figures 10
b, 10
d and 10
f) and the surface elevations (figure 11) do not possess any symmetry. Observe that while figures 10(b) and 11(a) may appear at first glance to suggest a symmetric response, this is due to the first mode shape (figure 2
i) being close to symmetric – there are in fact subtle asymmetries in figure 11.
Figures 12–14 also explore the square anisotropic plate, although this time the angle of the incident wave is modified to
$\chi =\pi /4$
. This incident wave is symmetric with respect to the line
$y=x$
, so modes that are antisymmetric with respect to this line should not resonate. Indeed, we observe that
$c_2=0$
, i.e. the second mode of vibration is not excited because it does have this antisymmetry (see figure 2). Moreover, there is no peak around
$\omega =9.98$
s
$^{-1}$
in the kinetic energy spectrum – which was associated with excitation of the second mode in the
$\chi =0$
case. We remark that the far-field patterns (figures 13
b and 13d) are symmetric with respect to
$\theta =\pi /4$
, and the surface elevations (figure 14) are symmetric about the line
$y=x$
.
Kinetic energy spectrum for a square anisotropic plate with clamped edges, with
$a=b=1$
m.
Dry mode expansion coefficient magnitudes (left column) and far-field patterns (right column) at resonant frequencies of the square anisotropic plate with clamped edges and with
$a=b=1$
m. The resonant frequencies are (a,b)
$\omega =6.51$
s
$^{-1}$
, (c,d)
$\omega =9.98$
s
$^{-1}$
and (e,f)
$\omega =10.73$
s
$^{-1}$
.
Surface elevation of the excited square anisotropic plate with clamped edges (
$a=b=1$
m) at the resonant frequencies considered in figure 10, namely (a)
$\omega =6.51$
s
$^{-1}$
, (b)
$\omega =9.98$
s
$^{-1}$
and (c)
$\omega =10.73$
s
$^{-1}$
.
Kinetic energy spectrum for a square anisotropic plate with clamped edges, with
$a=b=1$
m (as in figure 9) and with the direction of the incident wave being
$\pi /4$
.
Dry mode expansion coefficient magnitudes (left column) and far-field patterns (right column) at resonant frequencies of the square anisotropic plate with clamped edges and with
$a=b=1$
m, with the incident angle
$\pi /4$
. The resonant frequencies are (a)
$\omega =6.59$
s
$^{-1}$
and (b)
$\omega =10.70$
s
$^{-1}$
.
Surface elevation of the excited square anisotropic plate with clamped edges (
$a=b=1$
m) with incident angle
$\pi /4$
at the resonant frequencies considered in figure 13, namely (a)
$\omega =6.59$
s
$^{-1}$
and (b)
$\omega =10.70$
s
$^{-1}$
.
Next, we present results for rectangular clamped plates, with
$a=2$
m and
$b=1$
m, considering both the orthotropic and anisotropic cases. The first four mode shapes of both cases are given in figure 15. Their kinetic energy spectra are given in figures 16 and 19, each having three resonant peaks in the frequency interval considered. The relative excitations of the various plate modes at these resonant frequencies are given in figure 17 for the orthotropic case, and figure 20 for the anisotropic case. In both cases, we do not observe a strong resonant peak associated with the first plate mode. Moreover, in the orthotropic case, the fourth resonant mode is never excited because it is antisymmetric with respect to
$y=b/2$
, whereas the incident wave is symmetric with respect to this line. Finally, surface elevations of the first three resonances of orthotropic and anisotropic plates are given in figures 18 and 21, respectively.
First four mode shapes of an (a,c,e,g) orthotropic and (b,d,f,h) anisotropic rectangular plate, each with clamped edges and with
$a=2b=2$
m.
Kinetic energy spectrum for a rectangular orthotropic plate with clamped edges, with
$a=2$
m and
$b=1$
m.
Dry mode expansion coefficient magnitudes (left column) and far-field patterns (right column) at resonant frequencies of the rectangular orthotropic plate with clamped edges and with
$a=2b=2$
m. The resonant frequencies are (a,b)
$\omega =6.28$
s
$^{-1}$
, (c,d)
$\omega =7.88$
s
$^{-1}$
and (e,f)
$\omega =9.67$
s
$^{-1}$
.
Surface elevation of the excited rectangular orthotropic plate with clamped edges (
$a=2b=2$
m) at the resonant frequencies considered in figure 17, namely (a)
$\omega =6.28$
s
$^{-1}$
, (b)
$\omega =7.88$
s
$^{-1}$
and (c)
$\omega =9.67$
s
$^{-1}$
.
Kinetic energy spectrum for a rectangular anisotropic plate with clamped edges, with
$a=2b=2$
m.
Dry mode expansion coefficient magnitudes (left column) and far-field patterns (right column) at resonant frequencies of the rectangular anisotropic plate with clamped edges and with
$a=2\,m$
and
$b=1$
m. The resonant frequencies are (a,b)
$\omega =6.42$
s
$^{-1}$
, (c,d)
$\omega =8.09$
s
$^{-1}$
and (e,f)
$\omega =9.82$
s
$^{-1}$
.
Surface elevation of the excited rectangular anisotropic plate with clamped edges (
$a=2$
m and
$b=1$
m) at the resonant frequencies considered in figure 20, namely (a)
$\omega =6.42$
s
$^{-1}$
, (b)
$\omega =8.09$
s
$^{-1}$
and (c)
$\omega =9.82$
s
$^{-1}$
.
Next, we consider anisotropic plates with free edges, in both square (
$a=b=1$
m) and rectangular (
$a=2b=2$
m) geometries. Modes of vibration are illustrated in figure 22, which also includes the case of a square isotropic plate for the purpose of validation against Leissa (Reference Leissa1973, table C15). Kinetic energy spectra for square and rectangular anisotropic plates with free edges are given in figures 23 and 26, respectively. Each of these has only one resonant peak in the frequency interval considered, which is noticeably broader than peaks in the corresponding spectra for clamped plates (figures 9 and 19, respectively). Figures 24 and 27 show the dry mode coefficient amplitudes and far-field patterns at the resonant frequency. We observe that significant excitation occurs over many modes. In contrast, corresponding figures for clamped plates (figures 10 and 20) show that the resonances are predominantly due to excitation of a single mode. Figures 25 and 28 show the surface elevations for square and rectangular anisotropic plates with free edges, respectively.
First four non-rigid mode shapes and frequencies of vibration of plates with free edges. Panels (a,d,g,j) show the modes of a square isotropic plate (
$a=b=1$
m), (b,e,h,k) show the modes of a square anisotropic plate (
$a=b=1$
m) and (c,f,i,l) show the modes of a rectangular anisotropic plate (
$a=2b=2$
m). We note that, in each case, the first three modes, which are rigid body motions with frequency of vibration
$0$
s
$^{-1}$
, are not shown.
Kinetic energy spectrum for a square anisotropic plate with free edges, with
$a=b=1$
m.
Dry mode expansion coefficient magnitudes (a) and far-field patterns (b) at the resonant frequency
$\omega =8.18$
s
$^{-1}$
of the square anisotropic plate with free edges, with
$a=b=1\,$
m.
Surface elevation of the excited square anisotropic plate with free edges (
$a=b=1$
m) at the resonant frequency
$\omega =8.18$
s
$^{-1}$
. In contrast with earlier figures, the free boundary is marked with a black dashed line.
Kinetic energy spectrum for a rectangular anisotropic plate with free edges, with
$a=2b=2$
m.
Dry mode expansion coefficient magnitudes (a) and far-field patterns (b) at the resonant frequency
$\omega =7.22$
s
$^{-1}$
of the rectangular anisotropic plate with free edges, with
$a=2b=2\,$
m.
Finally, in figure 29 we present kinetic energy spectra for plates with simply supported edges. While this is primarily included to serve as a benchmark calculation, we note that the resonant frequencies are lower than those of clamped plates (e.g. compare figures 3 and 29 for isotropic plates), owing to their lower frequencies of vibration (Leissa Reference Leissa1973).
Surface elevation of the excited rectangular anisotropic plate with free edges (
$a=2b=2$
m) at the resonant frequency
$\omega =7.22$
s
$^{-1}$
.
9. Conclusion
This paper has considered the problem of water wave scattering by a surface-mounted rectangular anisotropic elastic plate, with either clamped or free boundaries. The problem was solved by expanding in the modes of vibration of the plate, with the required diffraction and radiation problems being solved using a boundary integral equation/constant panel method. Results for isotropic, orthotropic and anisotropic plates of either square or rectangular geometries were presented that illustrate the capabilities of the numerical method. Overall, our results show that the response by the plate can depend greatly on the rigidity coefficients
$D_{ij}$
. Our investigation has emphasised that the excitation of certain modes can be forbidden if they are symmetrically opposed to the incident wave.
Kinetic energy spectra for isotropic, orthotropic and anisotropic square plates (
$a=b=1$
m) with simply supported edges. Only the resonances of the isotropic case are marked.
The main motivation of this work was an application to the modelling of piezoelectric wave energy converters (PWECs), and this remains as an area of future work. The method presented here can be used (essentially verbatim) for a surface-mounted rectangular piezoelectric plate, with the rigidity coefficients
$D_{ij}$
becoming complex. The only necessary adjustment resulting from this would be to the dry mode expansion method presented in § 6, because the set of plate eigenmodes
$w_j$
may no longer be orthogonal in general, which would cause the hydrodynamic matrices
$K_{\textit{stiff}}$
,
$M$
and
$C$
to have non-zero off diagonal entries proportional to
$\iint _{\varGamma \times \{0\}} w_jw_m^*\,\mathrm{d} S$
. With this detail taken into account, we stress that our approach will allow the anisotropy of the piezoelectric material within the WEC to be modelled appropriately in three dimensions. We note that PWECs have not been considered in this paper because the determination of the complex rigidity coefficients from the poling direction, circuiting conditions and layering of the piezoelectric material is a complicated modelling problem (Renzi Reference Renzi2016). Another possible extension of this work is to arbitrary shaped plates, as considered by Meylan (Reference Meylan2002) for isotropic plates with free edges. Although the modes of vibration would need to be calculated using a different method to the one presented in § 4 (e.g. the finite element method could be used), the hydrodynamic component of the solution could remain unchanged. Lastly, we mention the problem of extending this study to submerged horizontal plates. This case is more complicated from a mathematical perspective because it leads to hypersingular integral equations.
Funding
This study was supported by the Australian Research Council’s Discovery Projects funding scheme (DP240102104).
Declaration of interests
The authors report no conflicts of interest.
Data availability
Data sharing does not apply to this article, as no new data were created or analysed in this study. The computer code used to generate the solutions is available from the corresponding author upon reasonable request.
Appendix A. Derivation of the Green’s function
In this appendix, we give a derivation of the Green’s function introduced in § 3.1. The equation satisfied by the Green’s function (3.2), when written in radially symmetric cylindrical coordinates, becomes
\begin{align} \partial _r^2 G(r,z)+\frac {1}{r}\partial _rG(r,z)+\partial _z^2 G(r,z)&=0&r\gt 0,\,-H\lt z\lt 0, \end{align}
\begin{align} \partial _zG(r,z)&=0&z=-H, \end{align}
\begin{align} \partial _zG(r,z)-\alpha G(r,z)&=-\frac {1}{2\pi r}\delta (r)&z=0, \end{align}
\begin{align} \sqrt {r}(\partial _r-\mathrm{i} k)G(r,z)&\to 0&r\to \infty . \end{align}
Let
$\hat {G}$
be the zeroth Hankel transform of
$G$
, i.e.
\begin{equation} \hat {G}(\mu ,z)=\int _0^\infty G(r,z) \mathrm{J}_0(\mu r)r\mathrm{d} r, \end{equation}
where
$\mathrm{J}_0$
is the Bessel function of the first kind of order zero. The inverse transform is given by
\begin{equation} G(r,z)=\int _0^\infty \hat {G}(\mu ,z) \mathrm{J}_0(\mu r)\mu \mathrm{d}\mu . \end{equation}
After transforming (A1a) and (A1b), we obtain
\begin{align} -\mu ^2\hat {G}(\mu ,z)+\partial _z^2\hat {G}(\mu ,z)&=0&\mu \gt 0,\,-H\lt z\lt 0, \end{align}
\begin{align} \partial _z\hat {G}(\mu ,z)&=0&z=-H, \\[9pt] \nonumber \end{align}
whereby
\begin{equation} \hat {G}(\mu ,z)=A(\mu )\cosh \mu (z+H). \end{equation}
The function
$A(\mu )$
must be determined from the boundary condition at
$z=0$
(A1c). The Hankel transform of this equation is
\begin{equation} \partial _z\hat {G}-\alpha \hat {G}=-\frac {1}{2\pi }, \end{equation}
for
$z=0$
< so we compute
\begin{align} A(\mu )&=-\frac {1}{2\pi }\frac {1}{\mu \sinh \mu H-\alpha \cosh \mu H}\nonumber \\ &=-\frac {{\mathrm{e}}^{-\mu H}}{2\pi \mu }\left (1+\frac {\mu +\alpha }{\mu \sinh \mu H-\alpha \cosh \mu H}\right ), \end{align}
where the second line follows from a substantial sequence of algebraic manipulations. Taking the inverse Hankel transform of
$\hat {G}$
(A3), with consideration of (A5) and (A7), yields
\begin{align} G(r,z)&=-\int _0^\infty \frac {{\mathrm{e}}^{-\mu H}}{2\pi \mu }\cosh \mu (z+H)\mathrm{J}_0(\mu r)\mu \mathrm{d}\mu \nonumber \\ &\qquad -\mathop {{\smile }}\!\!\!\!\!\!\mathop {\int }_0^\infty \frac {{\mathrm{e}}^{-\mu H}}{2\pi \mu }\frac {(\mu +\alpha )\cosh \mu H}{\mu \sinh \mu H-\alpha \cosh \mu H}\cosh \mu (z+H)\mathrm{J}_0(\mu r)\mu \mathrm{d}\mu . \end{align}
Note that we must integrate beneath the pole of the second integrand in order to satisfy the Sommerfeld radiation condition (A1d). After expressing
$\cosh \mu (z+H)= {1}/{2}({\mathrm{e}}^{\mu (z+H)}+{\mathrm{e}}^{-\mu (z+H)})$
and using the fact that (Bracewell Reference Bracewell2000)
\begin{equation} \int _0^{\infty }\frac {1}{\sqrt {r^2+\zeta ^2}}\mathrm{J}_0(\mu r) r\mathrm{d} r=\frac {{\mathrm{e}}^{-\mu |\zeta |}}{\mu }, \end{equation}
the first integral can be evaluated, yielding (3.3).
Appendix B. Beam eigenfunctions
The solutions of the eigenfunction problem associated with a beam clamped at both ends (4.5) may be expressed as
\begin{equation} u(x)=a_1\sin (\kappa x)+a_2\cos (\kappa x)+a_3{\mathrm{e}}^{-\kappa x}+a_4{\mathrm{e}}^{\kappa (x-1)}. \end{equation}
In order for the boundary conditions to be satisfied, we must find
$\kappa$
for which the following null-space problem has a non-trivial solution
\begin{equation} \begin{bmatrix} 0&1&1&{\mathrm{e}}^{-\kappa }\\ \sin (\kappa )&\cos (\kappa )&{\mathrm{e}}^{-\kappa }&1\\ 1&0&-1&{\mathrm{e}}^{-\kappa }\\ \cos (\kappa )&-\sin (\kappa )&-{\mathrm{e}}^{-\kappa }&1 \end{bmatrix}\begin{bmatrix} a_1\\a_2\\a_3\\a_4 \end{bmatrix}=\begin{bmatrix} 0\\0\\0\\0 \end{bmatrix}. \end{equation}
This is true when the determinant of the matrix is zero, which can be shown to occur when
\begin{equation} {\mathrm{sech}}(\kappa )=\cos (\kappa ). \end{equation}
We note that the clamped-clamped beam does not support modes for
$\kappa =0$
.
In the case that the beam is free at both ends, as in (4.10), the ansatz (B1), upon formation of a matrix analogous to that in (B2), gives rise to the same transcendental equation for the beam eigenvalues
$\kappa _m$
as the clamped case, namely
\begin{equation} {\mathrm{sech}}(\kappa )=\cos (\kappa ). \end{equation}
For non-zero roots of the above, the corresponding eigenmodes are of the form (B1). Note that the free–free beam also supports two modes for
$\kappa _1=\kappa _2=0$
, namely
\begin{align} u_1(x)&=1, \end{align}
\begin{align} u_2(x)&=2\sqrt {3}(x-\frac {1}{2}), \end{align}
which correspond to rigid body motions. The remaining modes
$u_j$
for
$j\geqslant 3$
are of the form (B1).
The eigenvalues and eigenmodes for the simply supported case (4.11) have the following simple form:
\begin{align} \kappa _n&=n\pi , \end{align}
\begin{align} u_n(x)&=\sqrt {2}\sin (n\pi ). \end{align}
We remark that the tasks of normalising of the eigenfunctions
$u_m$
in accordance with (4.5d), and computing the quantities
$\mathcal{U}_{mp}^{[j,l]}$
introduced in (4.8), require the computation of integrals that would be tedious to calculate by hand. Instead, symbolic computation was used to obtain suitable expressions, which were then verified by numerical integration.
Open access
