Boussinesq-type equations of hydroelastic waves in shallow water

Abstract

Accurate computation of hydroelastic waves in shallow water is critical because many hydroelastic wave applications are nearshores, such as sea-ice and floating infrastructures. In this paper, Boussinesq assumptions for shallow water are employed to derive nonlinear Boussinesq-type equations of hydroelastic waves, in which non-uniform distribution of structural stiffness and varying water depth are considered rigorously. Application of Boussinesq assumptions enables complicated three-dimensional problems to be reduced and formulated on the two-dimensional horizontal plane, therefore the proposed Boussinesq-type models are straightforward and versatile for a wide range of hydroelastic wave applications. Two configurations, a floating plate and a submerged plate, are studied. The first-order linear governing equations are solved analytically with periodic conditions assuming constant depth and uniform stiffness, and the linear dispersion relations are subsequently derived for both configurations. For flexural-gravity waves of a floating plate, unique behaviours of flexural-gravity waves different from shallow-water waves are discussed, and a generalized solitary wave solution is investigated. A nonlinear numerical solver is developed, and nonlinear flexural-gravity waves are found to have smaller wavelength and celerity than their linear counterparts. For hydroelastic waves of a submerged plate, dual-mode analytical solutions are discovered for the first time. Numerical computation has demonstrated that a plate with decreasing submerged depth is able to transfer wave energy from the deeper water to the surface layer.

1. Introduction

Hydroelastic waves propagate along the interface between elastic structures and hosting fluid. Flexural-gravity (FG) waves are a specific type of hydroelastic waves that occur in structures floating on the fluid free surface, dominated by both flexural and gravitational effects. Submerged hydroelastic (SH) waves are another type of hydroelastic waves that appear in structures fully submerged in the fluid. Most of existing studies on hydroelastic waves are based on potential flow theory, and thin plate models are commonly applied as wave amplitude and structure thickness are often small compared to hydroelastic wavelength.

FG waves have attracted significant research attentions due to their relevance to practical applications, such as analysis of very large floating structures (VLFS) (Ohmatsu 2005; Wang et al. 2010; Lamas-Pardo, Iglesias & Carral 2015; Bispo, Mohapatra & Soares 2021) and hydroelastic response of large sea-ice sheets in cold regions (Squire et al. 1995; Squire 2007, 2011, 2020). Flexural stiffness can hardly distribute uniformly because of variable thickness and structural irregularities. Different types of irregularities in floating plates and the resulting effects on FG waves have been investigated, such as cracks (Marchenko 1999; Squire & Dixon 2000; Evans & Porter 2003; Porter & Evans 2006, 2007; Ren, Wu & Li 2020; Barman et al. 2021), leads (Chuang & Linton 2005; Williams & Squire 2006; Shi, Li & Wu 2019; Zeng et al. 2021), ridges (Marchenko 1996; Marchenko & Voliak 1997; Williams & Squire 2004; Vaughan, Williams & Squire 2007; Squire & Williams 2008) and variable thickness (Porter & Porter 2004; Bennetts, Biggs & Porter 2007; Vaughan & Squire 2007; Smith & Meylan 2011). Apart from the above-mentioned investigations on periodic FG waves, solitary FG waves attracted renewed interest in the last two decades due to their relevance to the effects of a moving load and tsunami impact. Solitary FG wave solutions were derived theoretically using a fifth-order Korteweg–de Vries (KdV) equation (Hărăguş-Courcelle & Il'ichev 1998) and a cubic nonlinear Schorödinger equation (Guyenne & Părău 2014a,b), and computed numerically using Hamiltonian truncation methods (Vanden-Broeck & Părău 2011; Guyenne & Părău 2012; Wang, Milewski & Vanden-Broeck 2014), conformal mapping methods (Milewski, Vanden-Broeck & Wang 2011; Wang, Vanden-Broeck & Milewski 2013; Gao, Wang & Vanden-Broeck 2016), a high-order spectral model (Guyenne & Părău 2017), and a boundary integral technique (Trichtchenko et al. 2018). The existing studies on solitary FG waves are largely limited to finite-amplitude waves in deep water and small-amplitude waves in shallow water.

Investigations on SH waves of submerged plates are limited compared with FG waves of floating plates. Hassan, Meylan & Peter (2005) first developed the governing equations for semi-infinite and finite submerged elastic plates, and solved the two-dimensional (2-D) velocity potentials using the eigenfunction-matching method. Williams & Meylan (2012) considered the linear wave scattering by a semi-infinite submerged plate, and Smith et al. (2020) extended the Williams and Meylan solution to a submerged porous plate. Mohapatra, Sahoo & Soares (2018) and Mohapatra & Soares (2019) studied the interaction between surface gravity waves and an infinite submerged plate, and the interaction between FG and SH waves, respectively. Submerged plate structures have great potential for a wide range of engineering applications, such as submerged breakwaters, wave energy converters, wave cloaking and wave trapping. Our understanding on the hydroelastic response of submerged structures is still rather limited.

Most hydroelastic waves can be considered as long waves in the sense that (1) many hosting structures (e.g. ice covers and VLFS) are located in shallow-water regions, and (2) hydroelastic waves have larger wavelengths than gravity waves with the same frequency. Marchenko (1988) first derived a linear long wave model for FG waves in shallow water of constant depth, and considered inhomogeneous ice sheets in his latter studies (Marchenko 1996; Marchenko & Voliak 1997). Xia & Shen (2002) represented weakly nonlinear FG waves in shallow water of constant depth, and the resulting model was employed to study collision between two solitary FG waves by Bhatti & Lu (2018). Varying bathymetry has significant effects on propagation of hydroelastic waves in shallow water. Sturova (2008, 2009) investigated hydroelastic responses of homogeneous and heterogeneous elastic plates in shallow water of variable depth in the 2-D vertical plane. The wave responses of VLFS on water of variable depth were solved numerically in three-dimensional (3-D) domains using the finite element method (Kyoung et al. 2005), and in the 2-D vertical plane using various methods (Wang & Meylan 2002; Belibassakis & Athanassoulis 2006; Ding et al. 2017; Wu et al. 2018; Yang et al. 2019; Liu, Wang & Xu 2020). Recently, Ding et al. (2017), Wu et al. (2018) and Yang et al. (2019) developed a hybrid method in which the Boussinesq equations were applied to solve the shallow-water wave field of the outer region near islands and reefs, while the Rankine source method was used to solve the hydroelastic responses of VLFS in the inner region.

The aforementioned studies all fall into the framework of potential flow theory, in which hydroelastic and bathymetry effects are incorporated in the boundary conditions for Laplace's equation of velocity potential. When variable plate stiffness or varying bathymetry of shallow water are considered, the hydroelastic model becomes 3-D, and complicated numerical methods are inevitable to compute the 3-D solutions of velocity potential. In contrast, Boussinesq wave theory, one of the most successful shallow-water wave theories, applies the Boussinesq approximations and vertical integral to shallow-water waves for dimension reduction, and formulates the 3-D wave problems in the 2-D horizontal plane, enabling more efficient computation compared to 3-D potential flow theory. Boussinesq approximations assume that the square of the ratio between water depth and wavelength is of the same order as the ratio between wave amplitude and water depth, and both are small values for long waves in shallow water. The classic Boussinesq equations were derived from the Euler equations to model weakly dispersive and weakly nonlinear surface gravity waves in shallow water of variable depth (Peregrine 1967). Many researchers have contributed to developing various extended Boussinesq models to improve the dispersive properties (Madsen, Murray & Sørensen 1991; Madsen & Sørensen 1992; Nwogu 1993; Schäffer & Madsen 1995; Beji & Nadaoka 1996) and to enhance the nonlinearity (Wei et al. 1995; Madsen & Schäffer 1998; Zou 1999, 2000; Gobbi, Kirby & Wei 2000). However, to the best of the authors’ knowledge, the Boussinesq-type model has not yet been developed for hydroelastic waves. As the fundamental is the same (inviscid, incompressible, irrotational and continuous flow), velocity potential can be introduced in the derivation of Boussinesq models (Wei et al. 1995), yet we follow the original derivation in Peregrine (1967) without introducing velocity potential, which allows a straightforward non-dimensionalization process with explicit physical implications. Moreover, the small-displacement assumption is applied to the plate because wave amplitude is small compared to both wavelength and water depth in Boussinesq approximations. The thin-plate assumption is also applicable as the thickness is considered to be smaller than water depth and wavelength. Therefore in this study, the linear Kirchhoff–Love thin-plate theory is applied with Boussinesq approximations for hydroelastic waves in shallow water. Nonlinear elastic theories have been applied to study hydroelastic response to moving loads (Părău & Dias 2002; Bonnefoy, Meylan & Ferrant 2009), hydroelastic solitary waves (Milewski et al. 2011; Guyenne & Părău 2012, 2014b) and FG waves (Xia & Shen 2002; Toland 2008; Milewski & Wang 2013; Wang et al. 2013). However, we show in Appendix A that a nonlinear plate theory may not be necessary for typical Boussinesq-type models as the nonlinear elastic terms are higher-order terms that can be neglected under Boussinesq approximations, while the inertial term cannot be neglected.

Motivated to investigate the effects of variable plate stiffness and varying water depth on hydroelastic waves in shallow water, in this study, two new theoretical models that combine shallow-water Boussinesq approximations and Kirchhoff–Love plate theory are derived respectively for FG waves of an infinite floating plate and SH waves of an infinite submerged plate. Variable stiffness of the plates, varying submerged depth of the submerged plate, and varying bathymetry of shallow water are all considered rigorously. The resulting nonlinear governing equations are formulated in the 2-D horizontal plane, which allows easy application to nearshore problems and enables efficient numerical computations. The periodic analytical solutions to the first-order linear equations are derived, as well as the linear dispersion relations. A solitary solution to the nonlinear FG waves in shallow water is derived via analysis of an equivalent dynamical system, which is found to be small-amplitude as Wang et al. (2014) and Trichtchenko et al. (2018). Some unique behaviours of FG waves that are different from shallow-water waves are discussed, demonstrating the combining effects of gravitation and elastic flexion. We also reveal the dual-mode solutions to the gravity–hydroelastic coupled waves with a submerged plate, which are similar to the fast and slow modes of internal waves in the two-layer stratified fluids (Lamb 1932). A submerged plate with varying submerged depth is found to display the ability to allocate wave energy in the space domain, which is appealing to various engineering applications like wave cloaking and wave energy focusing. This paper is organized as follows. Derivations of the two new Boussinesq-type models are presented in § 2, with mathematical details given in Appendix A. The analyses of FG and SH waves are presented in §§ 3 and 4, respectively, including analytical solutions, numerical computations and new findings. The conclusions are given in § 5.

2. Boussinesq-type equations

Considering inviscid and incompressible fluid, the equations of motion of the fluid are presented by the Euler equations

(2.1)$$\begin{gather} \frac{{\partial} \boldsymbol{u}^*}{{\partial} t^*}+(\boldsymbol{u}^* \boldsymbol{\cdot} \boldsymbol{\nabla}^*) \boldsymbol{u}^* +w^*\left(\frac{{\partial} \boldsymbol{u}^*}{{\partial} z^*}\right) + \frac{1}{\rho}\,\boldsymbol{\nabla}^*P^* = 0, \end{gather}$$

(2.2)$$\begin{gather}\frac{{\partial} w^*}{{\partial} t^*}+(\boldsymbol{u}^* \boldsymbol{\cdot} \boldsymbol{\nabla}^*)w^* + w^*\left(\frac{{\partial} w^*}{{\partial} z^*}\right) + \frac{1}{\rho}\,\frac{{\partial} P^*}{{\partial} z^*} +g = 0, \end{gather}$$

in which superscript * indicates dimensional variables, the positive $z^*-$direction is vertically upwards (figure 1), $\boldsymbol {u}^*=(u^*, v^*)$ is a horizontal velocity vector, $w^*$ is the vertical velocity component, $\rho$ is fluid density, $P^*$ is fluid pressure, $g$ is the gravitational acceleration, and $\boldsymbol {\nabla }^*=({ {\partial }}/{ {\partial } x^*}, { {\partial }}/{ {\partial } y^*})$.

Figure 1. Floating elastic plate over varying bathymetry.

The continuity equation is

(2.3)\begin{equation} \boldsymbol{\nabla}^* \boldsymbol{\cdot} \boldsymbol{u}^* + \frac{{\partial} w^*}{{\partial} z^*} = 0.\end{equation}

Fluid motions are assumed to be irrotational such that

(2.4)\begin{equation} \frac{{\partial} \boldsymbol{u}^*}{{\partial} z^*} = \boldsymbol{\nabla}^* w^*. \end{equation}

The bottom boundary condition is given as

(2.5)\begin{equation} w^* + (\boldsymbol{u}^* \boldsymbol{\cdot} \boldsymbol{\nabla}^*)H^* = 0,\quad {\rm at}\ z^*=-H^*(x^*, y^*), \end{equation}

where $H^*(x^*, y^*)$ is the total depth representing a varying bathymetry.

The dimensional variables are non-dimensionalized using three different length scales: the representative wavelength $L_0$, typical water depth $H_0$, and wave amplitude $A_0$. The horizontal coordinates are scaled by $L_0$: $x=x^*/L_0$, $y=y^*/L_0$. The vertical coordinate and water depth are scaled by $H_0$: $z=z^*/H_0$, $H=H^*/H_0$. The floating plate displacement $\eta ^*(x^*, y^*, t^*)$ and submerged plate displacement $\zeta ^*(x^*, y^*, t^*)$ are scaled by $A_0$: $\eta =\eta ^*/A_0$, $\zeta =\zeta ^*/A_0$. Time is scaled by wave period, $t=t^*\sqrt {gH_0}/L_0$, and the non-dimensional fluid pressure is $P=P^*/(\rho g A_0)$. Two non-dimensional characteristic ratios, $\epsilon =A_0/H_0$ and $\mu =H_0/L_0$, are introduced. Then the horizontal velocity components are scaled as $u=u^*/(\epsilon \sqrt {gH_0})$, $v=v^*/(\epsilon \sqrt {gH_0})$, and the vertical velocity is scaled as $w=w^*/(\epsilon \mu \sqrt {gH_0})$. Mathematically speaking, the selection of scaling factors is arbitrary, yet an ideal non-dimensionalization scheme would enable non-dimensional parameters to have an order close to $O(1)$, allowing physical implications to be represented by the coefficient of each term in the non-dimensional equations. Therefore, typical realistic values are used as scaling factors in the above non-dimensionalization scheme. The resulting non-dimensional equations of motion and the bottom boundary condition are then obtained:

(2.6)$$\begin{gather} \frac{{\partial} \boldsymbol{u}}{{\partial} t} + \epsilon(\boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla})\boldsymbol{u} + \epsilon w\left(\frac{{\partial} \boldsymbol{u}}{{\partial} z}\right) + \boldsymbol{\nabla} P = 0, \end{gather}$$

(2.7)$$\begin{gather}\mu^2\,\frac{{\partial} w}{{\partial} t} + \epsilon\mu^2(\boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla})w + \epsilon\mu^2w \left(\frac{{\partial} w}{{\partial} z}\right) + \frac{{\partial} P}{{\partial} z} + \frac{1}{\epsilon} = 0, \end{gather}$$

(2.8)$$\begin{gather}\boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{u} + \frac{{\partial} w}{{\partial} z} = 0, \end{gather}$$

(2.9)$$\begin{gather}\frac{{\partial} \boldsymbol{u}}{{\partial} z} - \mu^2\,\boldsymbol{\nabla} w = 0, \end{gather}$$

(2.10)$$\begin{gather}w + (\boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla})H = 0,\quad {\rm at}\ z=-H(x, y). \end{gather}$$

2.1. Floating plate in shallow water

We first consider an infinite elastic flexible plate floating on surface of the fluid (figure 1). Assuming that no separation occurs, surface elevation is the same as displacement of the plate. The kinematic and dynamic surface boundary conditions resulting from the Kirchhoff–Love plate theory are given as

(2.11)$$\begin{gather} \frac{{\partial} \eta^*}{{\partial} t^*} + (\boldsymbol{u}^* \boldsymbol{\cdot} \boldsymbol{\nabla}^*)\eta^* = w^*,\quad {\rm at}\ z^*=\eta^*(x^*, y^*, t^*), \end{gather}$$

(2.12)$$\begin{gather}\alpha^*\,\frac{{\partial}^2 \eta^*}{{\partial} {t^*}^2} + {\nabla^*}^2(\beta^*\,{\nabla^*}^2\eta^*) + P^*_a - P^* = 0,\quad {\rm at}\ z^*=\eta^*(x^*, y^*, t^*), \end{gather}$$

in which ${\nabla ^*}^2={ {\partial }^2}/{ {\partial } {x^*}^2}+{ {\partial }^2}/{ {\partial } {y^*}^2}$ is the Laplacian operator, $P^*_a$ is the ambient air pressure on top of the plate, assumed to be zero, $\alpha ^*(x^*,y^*)=\rho _{plt}\, d^*(x^*,y^*)$ is the mass of the plate per unit area, $\beta ^*(x^*,y^*)={E{d^*}^3}/{12(1-\nu ^2)}$ is the flexural stiffness of the plate, $d^*$ is the thickness of the plate, assumed to be negligible, and $\rho _{plt}$, $E$ and $\nu$ are plate density, Young's modulus and Poisson's ratio, respectively. Introducing the additional non-dimensional variables

(2.13a,b)\begin{equation} \alpha=\alpha^*/(\rho H_0),\quad \beta=\beta^*/(\rho g{L_0}^4), \end{equation}

(2.11) and (2.12) can be non-dimensionalized as

(2.14)$$\begin{gather} \frac{{\partial} \eta}{{\partial} t} + \epsilon(\boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla})\eta -w = 0,\quad {\rm at}\ z=\epsilon\,\eta(x, y, t), \end{gather}$$

(2.15)$$\begin{gather}\mu^2\alpha\,\frac{{\partial}^2 \eta}{{\partial} {t}^2} + {\nabla}^2(\beta\,{\nabla}^2\eta) -P = 0,\quad {\rm at}\ z=\epsilon\,\eta(x, y, t). \end{gather}$$

An expansion method is employed to expand non-dimensional variables with respect to $\epsilon$ for the derivation of Boussinesq-type equations, which is similar to Peregrine's work (Peregrine 1967) but with major modifications to include the elastic and inertial effects of the plate. Details of the mathematical derivation are given in Appendix A. The resulting governing equations are presented in the following subsections.

2.1.1. The zeroth-order solutions and first-order equations

The zeroth-order solutions result from the still-water condition with no displacement:

(2.16a–c)\begin{equation} \boldsymbol{u}_0=\boldsymbol{Q}_0=\boldsymbol{0},\quad w_0=0\quad {\rm and}\quad P_0=-z/\epsilon, \end{equation}

in which $\boldsymbol {Q}$ is the horizontal flux of the flow, and the subscript indicates the associated order with respect to $\epsilon$ in the expansion series.

The first-order linear governing equation that describes FG waves propagating in a floating plate with variable stiffness over a varying bathymetry are derived as

(2.17)\begin{equation} \frac{{\partial}^2 \eta_1}{{\partial} t^2} - \boldsymbol{\nabla} \boldsymbol{\cdot}(H\,\boldsymbol{\nabla}\eta_1) - \boldsymbol{\nabla} \boldsymbol{\cdot}\{H\,\boldsymbol{\nabla}[\nabla^2(\beta\,\nabla^2\eta_1)]\} = 0. \end{equation}

Equation (2.17) is applicable to mild variations of water depth and plate stiffness, and can be regarded as the ‘mild-variation’ equation of FG waves, which is similar to the time-dependent mild-slope equation of shallow-water waves in Bin (1994). Equation (2.17) reduces to the standard mild-slope equation when $\beta =0$.

Once (2.17) is solved for $\eta _1$, other first-order variables can be obtained: $P_1$ can be solved from (A21); $\boldsymbol {u}_1$ and $\boldsymbol {Q}_1$ can be computed using (A22) and (A24); and $w_1$ can be calculated using the following equation that results from the integration of (A15),

(2.18)\begin{equation} w_1 =-z\,\boldsymbol{\nabla} \boldsymbol{\cdot}\boldsymbol{u}_1 - \boldsymbol{\nabla} \boldsymbol{\cdot}(H\boldsymbol{u}_1), \end{equation}

which satisfies the boundary conditions of (A18) and (A20).

2.1.2. Boussinesq-type equations

The higher-order nonlinear governing equations are derived in Appendix A, and the resulting governing equations are presented as

(2.19)\begin{gather} \begin{aligned} & \frac{{\partial} \bar{\boldsymbol{u}}}{{\partial} t} + \epsilon(\bar{\boldsymbol{u}}\boldsymbol{\cdot}\boldsymbol{\nabla}) \bar{\boldsymbol{u}} + \boldsymbol{\nabla}\eta + \boldsymbol{\nabla}[\nabla^2(\beta\,\nabla^2\eta)] + \mu^2\,\boldsymbol{\nabla}\left(\alpha\,\frac{{\partial}^2\eta}{{\partial} t^2}\right) \nonumber\\ &\quad = \mu^2\,\frac{H}{2}\,\frac{{\partial}}{{\partial} t} \boldsymbol{\nabla}[\boldsymbol{\nabla}\boldsymbol{\cdot}(H\bar{\boldsymbol{u}})] - \mu^2\,\frac{H^2}{6}\,\frac{{\partial}}{{\partial} t}\boldsymbol{\nabla}(\boldsymbol{\nabla}\boldsymbol{\cdot}\bar{\boldsymbol{u}}), \end{aligned} \end{gather}

(2.20)\begin{gather} \boldsymbol{\nabla}\boldsymbol{\cdot}[(H + \epsilon\eta)\bar{\boldsymbol{u}}] + \frac{{\partial} \eta}{{\partial} t}=0, \end{gather}

in which $\bar {\boldsymbol {u}}$ is the mean horizontal velocity defined by (A39).

Equations (2.19) and (2.20) are the Boussinesq-type equations for FG waves of a floating plate in shallow water, in which the elastic and inertial effects of the plate are represented respectively by the terms $\boldsymbol {\nabla }[\nabla ^2(\beta \,\nabla ^2\eta )]$ and $\mu ^2\,\boldsymbol {\nabla }(\alpha ({ {\partial }^2\eta }/{ {\partial } t^2}))$. These two coupled equations are used to solve for $\bar {\boldsymbol {u}}(x,y,t)$ and $\eta (x,y,t)$. They are reduced to the classic Boussinesq equations for shallow-water waves in Peregrine (1967) when $\alpha =\beta =0$. The Boussinesq-type equations can be represented by (A43) and (A44) in a dimensional form as well.

2.2. Submerged plate in shallow water

Consider an infinite elastic flexible plate submerged in the fluid (figure 2). The free-surface boundary conditions are applicable, whose non-dimensional forms are presented as

(2.21)$$\begin{gather} \frac{{\partial} \hat{\eta}}{{\partial} t} + \epsilon(\boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla}) \hat{\eta} - w = 0,\quad {\rm at}\ z=\epsilon\,\hat{\eta}(x, y, t), \end{gather}$$

(2.22)$$\begin{gather}P = P_a = 0,\quad {\rm at}\ z=\epsilon\,\hat{\eta}(x, y, t), \end{gather}$$

where $\hat {\eta }=\hat {\eta }^*/A_0$ denotes non-dimensional free-surface elevation, which is different from the floating plate displacement $\eta$ as they satisfy different dynamic boundary conditions.

Figure 2. Submerged elastic plate over varying bathymetry.

Due to the presence of the submerged plate, an additional set of boundary conditions is introduced to represent the kinematics and dynamics of the plate:

(2.23)$$\begin{gather} w^+- (\boldsymbol{u}^+ \boldsymbol{\cdot} \boldsymbol{\nabla})(\epsilon\zeta-H_I) = w^\unicode{x2013}- (\boldsymbol{u}^- \boldsymbol{\cdot}\boldsymbol{\nabla})(\epsilon\zeta-H_I) = \frac{{\partial} \zeta}{{\partial} t},\quad {\rm at}\ z=\epsilon\zeta-H_I, \end{gather}$$

(2.24)$$\begin{gather}\mu^2\alpha\,\frac{{\partial}^2 \zeta}{{\partial} {t}^2} + {\nabla}^2(\beta\,{\nabla}^2\zeta) + P^+- P^-= 0,\quad {\rm at}\ z=\epsilon\zeta-H_I, \end{gather}$$

in which $H_I=H_I^*/H_0$ is non-dimensional submerged depth, $w^+ = w|_{z=-H_I+\epsilon \zeta ^+}$, $w^- = w|_{z=-H_I+\epsilon \zeta ^-}$, $\boldsymbol {u}^+ = \boldsymbol {u}|_{z=-H_I+\epsilon \zeta ^+}$ and $\boldsymbol {u}^- = \boldsymbol {u}|_{z=-H_I+\epsilon \zeta ^-}$ are fluid velocities at the upper and lower surfaces of the plate, and $P^+ = P|_{z=-H_I+\epsilon \zeta ^+}$ and $P^- = P|_{z=-H_I+\epsilon \zeta ^-}$ are fluid pressures on the upper and lower plate surfaces.

The flow field is separated by the submerged plate into two layers, each of which has its own boundary conditions but conforms to the same equations of motion as (2.6)–(2.9). We define the horizontal flux of the upper flow as $\boldsymbol {Q}_I=\int _{-H_I+\epsilon \zeta ^+}^{\epsilon \hat {\eta }} \boldsymbol {u}\,{\rm d} z$, and the lower horizontal flux as $\boldsymbol {Q}_{II}=\int _{-H}^{-H_I+\epsilon \zeta ^-} \boldsymbol {u}\,{\rm d} z$. Therefore, integrating (2.8) from $z=-H_I+\epsilon \zeta ^+$ to $z=\epsilon \hat {\eta }$, and from $z=-H$ to $z=-H_I+\epsilon \zeta ^-$, respectively, we have

(2.25)$$\begin{gather} \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{Q}_I + \frac{{\partial} \hat{\eta}}{{\partial} t} - \frac{{\partial} \zeta}{{\partial} t} = 0, \end{gather}$$

(2.26)$$\begin{gather}\boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{Q}_{II} + \frac{{\partial} \zeta}{{\partial} t} = 0, \end{gather}$$

in which the boundary conditions of $w$ from (2.10), (2.21) and (2.23) have been applied.

2.2.1. The zeroth-order solutions and first-order equations

The free surface elevation $\hat {\eta }$, the submerged plate displacement $\zeta$, and the horizontal fluxes $\boldsymbol {Q}_I$ and $\boldsymbol {Q}_{II}$ are similarly expanded with respect to $\epsilon$ as in Appendix A. The system of governing equations is then obtained by substituting the expanded variables into the equations of motion and boundary conditions, including (2.6), (2.7), (2.9), (2.10) and (2.21)–(2.26).

The zeroth-order solution of these governing equations is the still-water condition:

(2.27a–d)\begin{align} \boldsymbol{u}_0=\boldsymbol{Q}_{I0}=\boldsymbol{Q}_{II0}=\boldsymbol{0},\quad w_0=0,\quad P_0=-z/\epsilon\quad {\rm and}\quad P_0^+=P_0^-=-H_I/\epsilon. \end{align}

Following steps similar to those in Appendix A, we can derive the first-order linear equations describing propagation of the coupled surface gravity wave and SH wave in shallow water. The resulting equations are

(2.28)$$\begin{gather} \frac{{\partial}^2 \hat{\eta}_1}{{\partial} t^2} - \frac{{\partial}^2 \zeta_1}{{\partial} t^2} - \boldsymbol{\nabla} \boldsymbol{\cdot}(H_I\,\boldsymbol{\nabla}\hat{\eta}_1) = 0, \end{gather}$$

(2.29)$$\begin{gather}\frac{{\partial}^2 \zeta_1}{{\partial} t^2} - \boldsymbol{\nabla}\boldsymbol{\cdot}(H_{II}\,\boldsymbol{\nabla}\hat{\eta}_1) - \boldsymbol{\nabla}\boldsymbol{\cdot}\{H_{II}\,\boldsymbol{\nabla}[\nabla^2(\beta\,\nabla^2\zeta_1)]\}= 0, \end{gather}$$

in which $H_{II}=H-H_I$ is the non-dimensional distance from the submerged plate to the bottom. Equation (2.29) reduces to (2.17) when $H_I$ goes to zero and $\zeta _1$ becomes $\eta _1$. Similar to (2.17), then (2.28) and (2.29) are the ‘mild-variation’ equations of SH waves applicable to mild variations of $H_I$, $H_{II}$ and $\beta$. The first-order surface elevation and submerged plate displacement are solved using (2.28) and (2.29), and other first-order variables can be computed subsequently.

2.2.2. Boussinesq-type equations

The higher-order equations can be derived similarly to the derivation in Appendix A. Because the submerged plate separates the flow into two layers, the number of unknown variables is doubled. We define the mean horizontal velocities of the upper and lower flows as

(2.30a,b)\begin{equation} \bar{\boldsymbol{u}}_I = \frac{\boldsymbol{Q}_I}{H_I+\epsilon\hat{\eta}-\epsilon\zeta}\quad {\rm and}\quad \bar{\boldsymbol{u}}_{II} = \frac{\boldsymbol{Q}_{II}}{H_{II} + \epsilon\zeta}. \end{equation}

The Boussinesq-type equations for the gravity–hydroelastic coupled wave system that consists of surface gravity waves and SH waves are derived in terms of mean horizontal velocities and surface displacements:

(2.31)\begin{gather} \frac{{\partial} \bar{\boldsymbol{u}}_I}{{\partial} t} + \epsilon(\bar{\boldsymbol{u}}_I\boldsymbol{\cdot}\boldsymbol{\nabla}) \bar{\boldsymbol{u}}_I + \boldsymbol{\nabla}\hat{\eta} = \mu^2\,\frac{H_I}{2}\,\frac{{\partial}}{{\partial} t} \boldsymbol{\nabla}[\boldsymbol{\nabla}\boldsymbol{\cdot}(H\bar{\boldsymbol{u}}_I)] - \mu^2\,\frac{H_I^2}{6}\,\frac{{\partial}}{{\partial} t} \boldsymbol{\nabla}(\boldsymbol{\nabla}\boldsymbol{\cdot}\bar{\boldsymbol{u}}_I), \end{gather}

(2.32)\begin{gather} \begin{aligned} & \frac{{\partial} \bar{\boldsymbol{u}}_{II}}{{\partial} t} + \epsilon(\bar{\boldsymbol{u}}_{II}\boldsymbol{\cdot}\boldsymbol{\nabla})\bar{\boldsymbol{u}}_{II} + \boldsymbol{\nabla}\hat{\eta} + \boldsymbol{\nabla}[\nabla^2(\beta\,\nabla^2\zeta)] + \mu^2\,\boldsymbol{\nabla}\left(\alpha\,\frac{{\partial}^2\zeta}{{\partial} t^2}\right) \nonumber\\ &\quad =\mu^2\,\frac{(H^2-H_I^2)}{2H_{II}}\,\frac{{\partial}}{{\partial} t} \boldsymbol{\nabla}[\boldsymbol{\nabla}\boldsymbol{\cdot}(H\bar{\boldsymbol{u}}_{II})] - \mu^2\, \frac{(H^3-H_I^3)}{6H_{II}}\,\frac{{\partial}}{{\partial} t}\boldsymbol{\nabla} (\boldsymbol{\nabla}\boldsymbol{\cdot}\bar{\boldsymbol{u}}_{II}), \end{aligned}\end{gather}

(2.33)$$\begin{gather} \boldsymbol{\nabla}\boldsymbol{\cdot}[(H_I+\epsilon\hat{\eta}-\epsilon\zeta)\bar{\boldsymbol{u}}_I] + \frac{{\partial} \hat{\eta}}{{\partial} t} - \frac{{\partial} \zeta}{{\partial} t} = 0, \end{gather}$$

(2.34)$$\begin{gather}\boldsymbol{\nabla}\boldsymbol{\cdot}[(H_{II}+\epsilon\zeta)\bar{\boldsymbol{u}}_{II}] + \frac{{\partial} \zeta}{{\partial} t} = 0. \end{gather}$$

3. Analyses of FG waves in shallow water

The solutions to FG waves in shallow water are studied in this section, including the linear dispersion relation and nonlinear effects. The existence of a solitary wave solution is also investigated, and the solitary waves propagating through varying stiffness and depth are demonstrated.

3.1. First-order analytical solution

Analytical solutions are highly desired because they enable straightforward computation and provide meaningful insights into the fluid–structure coupled dynamics. However, it is difficult to analytically solve the nonlinear Boussinesq-type equations. Instead, we derive the analytical solution to the first-order, linear governing equation, while numerically solving the nonlinear Boussinesq-type equations of FG waves.

The first-order, linear, non-dimensional governing equation of shallow-water FG waves is given by (2.17). Assuming that water depth $H$ and plate stiffness $\beta$ are constants that do not vary in space, (2.17) is reduced to

(3.1)\begin{equation} \frac{{\partial}^2 \eta_1}{{\partial} t^2} - H\,\frac{{\partial}^2 \eta_1}{{\partial} x^2} - H\beta\,\frac{{\partial}^6 \eta_1}{{\partial} x^6} = 0, \end{equation}

in which $x$ is the direction of wave propagation. Assume that the solution to (3.1) is periodic: $\eta _1=\eta _a \exp ({{\rm i}(kx\pm \sigma t)})$, in which $k$ and $\sigma$ are non-dimensional wavenumber and angular frequency, respectively, and $\eta _a$ is the amplitude determined by boundary and initial conditions. Substituting the periodic solution into (3.1), we can derive the first-order dispersion relation for FG waves in shallow water:

(3.2)\begin{equation} \sigma^2=Hk^2+H\beta k^6. \end{equation}

Equation (3.2) reduces to the dispersion relation of linear shallow-water waves when flexural stiffness $\beta$ goes to zero. Substituting $\gamma =k^2$ into (3.2), and solving the resulting cubic equation for the real root of $\gamma$ using Cardano's formula, we can compute the wavenumber for the given wave frequency, water depth and plate stiffness:

(3.3)\begin{equation} {k}^2=\sqrt[3]{\frac{{\sigma}^2}{2H\beta}-\sqrt{\left(\frac{{\sigma}^2}{2H\beta}\right)^2+ \left(\frac{1}{3\beta}\right)^3}}+\sqrt[3]{\frac{{\sigma}^2}{2H\beta}+ \sqrt{{\left(\frac{{\sigma}^2}{2H\beta}\right)}^2+{\left(\frac{1}{3\beta}\right)}^3}}. \end{equation}

Additionally, assume that the periodic horizontal velocity is $u_1=u_a \exp ({{\rm i}(kx\pm \sigma t)})$, where $u_a$ is the amplitude of the velocity. Substituting the periodic solutions of $\eta _1$ and $u_1$ into (A16), (A22) and (A24), and combining the resulting equations, we have

(3.4)\begin{equation} u_a=\frac{k+\beta k^5}{\sigma}\,\eta_a=\frac{\sigma}{Hk}\,\eta_a. \end{equation}

Equations (3.3) and (3.4) offer the complete periodic solutions to the first-order linear system of FG waves in shallow water.

3.2. Linear dispersion relations

Equation (3.2) presents the dispersion relation derived from the first-order equation, in which all the terms associated with $\epsilon$ and $\mu ^2$ in the Boussinesq-type equations are neglected. Now we derive the linear dispersion relation from the reduced equations of (2.19) and (2.20), in which only the nonlinear terms associated with $\epsilon$ are neglected. Similarly assuming constant depth and uniform plate, the reduced equations are

(3.5)$$\begin{gather} \frac{{\partial} \bar{u}_x}{{\partial} t} + \frac{{\partial} \eta}{{\partial} x} + \beta\,\frac{{\partial}^5 \eta}{{\partial} x^5} + \mu^2\alpha\,\frac{{\partial}^3\eta}{{\partial} x\,{\partial} t^2} - \mu^2\,\frac{H^2}{3}\,\frac{{\partial}^3 \bar{u}_x}{{\partial} x^2\,{\partial} t} = 0, \end{gather}$$

(3.6)$$\begin{gather}H\,\frac{{\partial} \bar{u}_x}{{\partial} x} + \frac{{\partial} \eta}{{\partial} t} = 0, \end{gather}$$

in which $\bar {u}_x$ is the $x$-component of $\bar {\boldsymbol {u}}$ in the propagating direction.

The linear dispersion relation is obtained by substituting harmonically periodic $\bar {u}_x$ and $\eta$ into (3.5) and (3.6):

(3.7)\begin{equation} \sigma^2=\frac{Hk^2+H\beta k^6}{1+\mu^2H\alpha k^2+\tfrac{1}{3}\mu^2k^2H^2}. \end{equation}

Similar to § 3.1, the linear analytic solutions of $\bar {u}_x$ and $\eta$ can be obtained by solving the linear dispersion relation (3.7).

The linear dispersion relation of FG waves for finite water depth was derived by Fox & Squire (1994) and Sahoo, Yip & Chwang (1994), and it can be non-dimensionalized as

(3.8)\begin{equation} \sigma^2=(Hk^2-\mu^2H\alpha k^2\sigma^2+H\beta k^6)\,\frac{\tanh(\mu kH)}{\mu kH}. \end{equation}

When $\mu$ is small such that the [$1/2$]-order Padé approximation of $\tanh (\mu kH)$ gives ${\mu k H}/{(1+(\mu kH)^2/3)}$, (3.8) is reduced to

(3.9)\begin{equation} \sigma^2=(Hk^2-\mu^2H\alpha k^2\sigma^2+H\beta k^6)\,\frac{1}{1+\tfrac{1}{3}\mu^2k^2H^2}, \end{equation}

which is equivalent to (3.7) if $\sigma ^2$ is moved to the left-hand side of the equation. Furthermore, (3.7) and(3.8) both converge to (3.2) if the higher-order terms associated with $\mu ^2$ are neglected.

The phase and group velocities of FG waves can be derived using different linear dispersion relations for shallow water and finite depth, respectively. For the first-order shallow-water model, its phase and group velocities are derived from (3.2) as

(3.10)$$\begin{gather} C_{p\_first}^2 = \frac{\sigma^2}{k^2}=H+H\beta k^4, \end{gather}$$

(3.11)$$\begin{gather}C_{g\_first} = \frac{{\partial} \sigma}{{\partial} k}=C_{p\_first}\,\frac{1+3\beta k^4}{1+\beta k^4}. \end{gather}$$

For the linear shallow-water model, its phase and group velocities are derived from (3.7):

(3.12)$$\begin{gather} C_{p\_linear}^2 = \frac{H+H\beta k^4}{1+\mu^2H\alpha k^2+\tfrac{1}{3}\mu^2k^2H^2}, \end{gather}$$

(3.13)$$\begin{gather}C_{g\_linear} = C_{p\_linear}\Big(\frac{2\beta k^4}{1+\beta k^4}+\frac{1}{1+\mu^2 H\alpha k^2+ \tfrac{1}{3}\mu^2 k^2H^2}\Big)\!. \end{gather}$$

For the linear finite depth model, its phase and group velocities are derived from (3.8):

(3.14)$$\begin{gather} C_{p\_finite}^2 = \frac{(H+H\beta k^4)\,\dfrac{\tanh(\mu kH)}{\mu kH}}{1+\mu^2H\alpha k^2\, \dfrac{\tanh(\mu kH)}{\mu kH}}, \end{gather}$$

(3.15)$$\begin{gather}C_{g\_finite} = C_{p\_finite}\Big[\frac{2\beta k^4}{1+\beta k^4}+\frac{\mu k H +\tanh(\mu k H)- \mu k H\tanh^2(\mu k H)}{2\tanh(\mu k H)+2\mu\alpha k\tanh^2(\mu k H)}\Big]\!. \end{gather}$$

To evaluate the effectiveness of the proposed shallow-water models for various relative depth $\mu$, we compute the phase and group velocities of FG waves propagating in sea-ice (table 1) using the first-order, linear shallow-water and linear finite-depth dispersion relations. The phase and group velocities of linear water waves (Airy waves) are included for comparison as well. Figure 3 shows the changes of phase and group velocities against $\mu$ with different ice thicknesses when a typical wave period $T^*=10$ s is considered. The phase and group velocities predicted by the first-order dispersion relation generally agree well with those of the finite-depth model up to $\mu =0.05$ (1.7 % error in phase velocity for $d^*=1.0$ m), while the linear shallow-water model can be applied up to $\mu =0.2$ (2.0 % error in phase velocity for $d^*=1.0$ m). It is noticed that the elastic effects are prominent in shallow water when we compare the water waves with FG waves: the predictions by Airy wave theory are closer to FG waves in finite depth for the larger $\mu$.

Table 1. Physical parameters of floating sea-ice.

Figure 3. Comparison of phase and group velocities computed using different wave models with various $\mu$ ($T^*=10$ s).

The effects of wave frequency are also investigated. Figure 4 shows the phase and group velocities of water waves on surface and FG waves in sea-ice for different wave periods. The wavelength of FG waves increases with wave period but its velocities do not, and the group velocity of FG waves is greater than the phase velocity for higher frequencies. These characteristics are completely different from water waves. For pure bending waves, the phase velocity is $\sqrt [4]{\sigma ^2\beta /\alpha }$, which decreases with wave period, and the group velocity is twice as large as the phase velocity. The FG waves result from the combining effects of gravitation and elastic flexion, and it is noticed that the elastic effects are prominent for higher frequencies.

Figure 4. Comparison of phase and group velocities computed for different wave periods: (a,b) $d^*=0$ (water wave); (c,d) $d^*=0.3$ m (FG wave); (e,f) $d^*=1.0$ m (FG wave).

3.3. Nonlinear numerical solutions

The nonlinear Boussinesq-type equations of FG waves in shallow water are solved numerically. Equations (2.19) and (2.20) are computed using an implicit, time-marching scheme:

(3.16)\begin{gather} \eta_{s+1^*} = \eta_{s} - \Delta t\,\boldsymbol{\nabla} \boldsymbol{\cdot} [(\epsilon\eta_{s}+H)\bar{\boldsymbol{u}}_{s}], \end{gather}

(3.17)\begin{gather} \begin{aligned}& \frac{\bar{\boldsymbol{u}}_{s+1}-\bar{\boldsymbol{u}}_{s}}{\Delta t} +\epsilon(\bar{\boldsymbol{u}}_s\boldsymbol{\cdot}\boldsymbol{\nabla})\bar{\boldsymbol{u}}_s -\frac{\mu^2}{2}\,H\,\boldsymbol{\nabla}\left[\boldsymbol{\nabla}\boldsymbol{\cdot}\left(H\,\frac{\bar{\boldsymbol{u}}_{s+1}- \bar{\boldsymbol{u}}_s}{\Delta t}\right)\right] \nonumber\\ &\qquad +\frac{\mu^2}{6}\,H^2\,\boldsymbol{\nabla} \left[\boldsymbol{\nabla}\boldsymbol{\cdot}\left(\frac{\bar{\boldsymbol{u}}_{s+1}-\bar{\boldsymbol{u}}_s}{\Delta t}\right)\right] \nonumber\\ &\quad =-\boldsymbol{\nabla}\left(\frac{\eta_{s}+\eta_{s+1^*}}{2}\right) - \boldsymbol{\nabla}\left\{\nabla^2\left[\beta\,\nabla^2\left(\frac{\eta_{s}+\eta_{s+1^*}}{2}\right)\right]\right\} \nonumber\\ &\qquad -\mu^2\,\boldsymbol{\nabla}\left(\alpha\,\frac{\eta_{s+1^*}-2\eta_{s}+\eta_{s-1}}{\Delta t^2}\right)\!, \end{aligned}\end{gather}

(3.18)\begin{gather} \eta_{s+1} = \eta_{s} - \Delta t\,\boldsymbol{\nabla} \boldsymbol{\cdot} [(\epsilon\eta_{s+1^*}+H)\bar{\boldsymbol{u}}_{s+1}], \end{gather}

in which $\Delta t$ is the time step size, subscript $s$ indicates the $s$th time step, and $\eta _{s+1^*}$ is the provisional value of $\eta _{s+1}$. Equation (3.17) is solved implicitly to update $\bar {\boldsymbol {u}}_{s+1}$. Application of the provisional step allows conditionally stable computation of (2.19) and (2.20) using the central difference schemes for spatial derivatives.

To demonstrate the nonlinear effects, FG waves propagating in sea-ice in shallow water are computed using the first-order linear analytical solution and the proposed nonlinear numerical scheme, respectively. The parameters used in these computations are presented in tables 1 and 2, and thickness 0.1 m is employed. By changing water depth and the amplitude of boundary excitation, we compare the linear analytical solutions and the nonlinear Boussinesq solutions based on different values of $\epsilon$ and $\mu$ that generally conform to the Boussinesq approximations: $O(\epsilon ) \approx O(\mu ^2)\ll O(1)$. Figure 5 shows that the nonlinear FG waves have smaller wavelength than their first-order linear counterparts.

Table 2. Simulation cases.

Figure 5. Comparison between the linear and nonlinear solutions of FG waves with various $\epsilon$ and $\mu$. (a) $\epsilon = 0.0010,\ \mu = 0.0309$; (b) $\epsilon = 0.0051,\ \mu = 0.0713$; (c) $\epsilon = 0.0100,\ \mu = 0.1000$; (d) $\epsilon = 0.0204, \mu = 0.1428$.

The shoaling effect of FG waves is illustrated by the nonlinear numerical solutions. Figure 6 presents the snapshots of nonlinear FG waves in sea-ice with different thicknesses propagating over decreasing water depth. The wavelength reduces in all three cases during shoaling, regardless of the different ice thickness, as predicted by the linear dispersion relation. However, the amplitude of FG waves does not necessarily increase during shoaling. For flexural deflection, its elastic strain energy is associated with the curvature of deflection, which increases with decreasing wavelength if wave amplitude remains. When the elastic effect is dominant in an FG wave system (with relatively large plate stiffness), its amplitude must decrease with decreasing wavelength during shoaling for the conversation of energy. Recalling the relation between wave velocities and wave frequency shown in figure 4, this is another unique behaviour of shallow FG waves that is different from shallow-water waves, demonstrating the combining effects of gravitation and elastic flexion.

Figure 6. Shoaling effect of nonlinear FG waves. (a) $d^{\ast} = 0.1\ {\rm m}$; (b) $d^{\ast} = 0.3\ {\rm m}$; (c) $d^{\ast} = 0.5\ {\rm m}$; (d) Varying depth.

3.4. Solitary wave solution

A solitary wave solution is derived for nonlinear FG waves in shallow water. Similarly to the derivation of linear analytical solutions, we assume that water depth $H$, flexural stiffness $\beta$, and unit mass $\alpha$ are positive constants that do not change with respect to $x$ and $y$, such that (2.19) and (2.20) are reduced to

(3.19)$$\begin{gather} \frac{{\partial} \bar{u}_x}{{\partial} t} + \epsilon\bar{u}_x\,\frac{{\partial} \bar{u}_x}{{\partial} x} + \frac{{\partial} \eta}{{\partial} x} + \beta\,\frac{{\partial}^5 \eta}{{\partial} x^5} + \mu^2\alpha\,\frac{{\partial}^3\eta}{{\partial} x\,{\partial} t^2} - \frac{\mu^2 H^2}{3}\,\frac{{\partial}^3 \bar{u}_x}{{\partial} x^2\,{\partial} t} = 0, \end{gather}$$

(3.20)$$\begin{gather}\epsilon\bar{u}_x\,\frac{{\partial} \eta}{{\partial} x} + \epsilon\eta\,\frac{{\partial} \bar{u}_x}{{\partial} x} + H\,\frac{{\partial} \bar{u}_x}{{\partial} x} + \frac{{\partial} \eta}{{\partial} t} = 0. \end{gather}$$

Introducing a moving coordinate $\xi =x-ct$, we have ${ {\partial }}/{ {\partial } x}={ {\partial }}/{ {\partial } \xi }$ and ${ {\partial }}/{ {\partial } t}=-c({ {\partial }}/{ {\partial } \xi })$, such that (3.19) and (3.20) are rewritten as

(3.21)$$\begin{gather} -c\bar{u}_x' + \epsilon\bar{u}_x\bar{u}_x' + \eta' + \beta\eta^{(5)} + \mu^2c^2\alpha\eta''' + \frac{\mu^2 c H^2}{3}\,\bar{u}_x''' = 0, \end{gather}$$

(3.22)$$\begin{gather}\epsilon\bar{u}_x\eta' + \epsilon\eta\bar{u}_x' + H\bar{u}_x' - c\eta' = 0, \end{gather}$$

in which the superscripts indicate the order of derivative with respect to $\xi$.

Equation (3.22) is linearized by assuming $O(\epsilon )\ll O(1)$, which gives $H\bar {u}_x'-c\eta '=0$. Because solitary solutions approach zero for $\xi \rightarrow \pm \infty$, integral constants must be zero and we have $\bar {u}_x=({c}/{H})\eta$, $\bar {u}_x'=({c}/{H})\eta '$ and $\bar {u}_x'''=({c}/{H})\eta '''$. Then (3.21) can be rewritten as a fifth-order nonlinear ordinary differential equation (ODE) with respect to $\eta$:

(3.23)\begin{equation} \left(1-\frac{c^2}{H}\right)\eta' + \epsilon\,\frac{c^2}{H^2}\,\eta\eta' + \mu^2 c^2\left(\alpha+\frac{H}{3}\right)\eta''' + \beta\eta^{(5)} = 0. \end{equation}

It is noted that (3.23) is also equivalent to a generalized fifth-order KdV equation:

(3.24)\begin{equation} \left(\frac{c}{H}-\frac{1}{c}\right)\frac{{\partial}\eta}{{\partial} t} + \epsilon\,\frac{c^2}{H^2}\,\eta\,\frac{{\partial}\eta}{{\partial} x} + \mu^2 c^2 \left(\alpha+\frac{H}{3}\right)\frac{{\partial}^3\eta}{{\partial} x^3} + \beta\frac{{\partial}^5\eta}{{\partial} x^5} = 0. \end{equation}

Substituting a general harmonic solution $\eta =A\exp ({{\rm i}(kx+\sigma t)})$ into the linear part of (3.24), we can derive the approximate dispersion relation of the system: $({(c^2-H)}/{Hc})\sigma =\mu ^2 c^2(\alpha +{H}/{3})k^3+\beta k^5$. The system bifurcates when ${ {\partial } c}/{ {\partial } k}={ {\partial }^2 \sigma }/{ {\partial } k^2}=0$, such that a solitary solution exists at $k\rightarrow 0$ as $c^2\rightarrow H$. Substituting $c^2=H$ into (3.23) and integrating once yields

(3.25)\begin{equation} \eta^{(4)} =-\left(\frac{\mu^2\alpha H}{\beta}+\frac{\mu^2H^2}{3\beta}\right)\eta'' - \frac{\epsilon}{2\beta H}\,\eta^2. \end{equation}

An equivalent dynamical system of (3.25) is derived as

(3.26)\begin{equation} \begin{Bmatrix} \eta \\ \eta' \\ \eta'' \\ \eta''' \end{Bmatrix}'= \begin{bmatrix} 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \\ 0 & 0 & -\left(\dfrac{\mu^2\alpha H}{\beta}+\dfrac{\mu^2 H^2}{3\beta}\right) & 0 \end{bmatrix} \begin{Bmatrix} \eta \\ \eta' \\ \eta'' \\ \eta''' \end{Bmatrix}+ \begin{Bmatrix} 0 \\ 0 \\ 0 \\ -\dfrac{\epsilon\eta^2}{2\beta H} \end{Bmatrix}. \end{equation}

The eigenvalues of the matrix in (3.26) are 0, 0 and $\pm ik_o$, given that $\alpha$, $\beta$, $\epsilon$ and $H$ are positive real numbers. Because the dynamical system has a pair of imaginary eigenvalues, a generalized solitary solution appears, which consists of a ‘${\rm sech}^2$’ core soliton that decays at infinity, and a small-amplitude oscillatory tail with wavenumber $k_o=\sqrt {{\mu ^2\alpha H}/{\beta }+{\mu ^2 H^2}/{3\beta }}$ (Il'ichev 1997, 2000). The generalized solitary solution for FG waves in shallow water is derived as

(3.27)\begin{equation} \eta_{sol} = \eta_c+\eta_o = 3H\,\frac{c^2-H}{\epsilon c^2}\,{\rm sech}^2 \left[\frac{\tau}{2}\,(x-ct)\right] + A_o\exp\left({-\frac{\rm \pi}{\tau}\,k_o}\right)\exp({{\rm i}k_o(x-ct)}), \end{equation}

in which $\tau =\sqrt {{3(c^2-H)}/{\mu ^2 c^2 H(3\alpha +H)}}$, and $A_o\exp ({-({{\rm \pi} }/{\tau })k_o})$ is the amplitude of the oscillatory tail that is exponentially small with respect to $\sqrt {\beta }$. We present the derivation of the solitary solution $\eta _{sol}$ in Appendix B.

The nonlinear solver developed in § 3.3 is employed to compute the solitary wave given by (3.27), in which the physical parameters from table 1 and water depth 10 m are used. The core soliton $\eta _c$ is used as the initial condition for numerical computations to evaluate the evolution of the solitary solution. The propagating speed $c$ plays an important role in the solitary solution. The solitary solution exists only if $c^2\rightarrow H$. The computation results show that the initial soliton hardly changes when $c$ is close to $\sqrt {H}$ (figure 7a), but breaks down into two major peaks when $c$ gets larger (figure 7b): one peak moves at a lower speed 1.08$\sqrt {H}$ and maintains its shape, while the other propagates at a higher speed 1.32$\sqrt {H}$ and contentiously decays. The amplitude of the core soliton is also dependent on $c$, therefore only small-amplitude solitary waves are permitted.

Figure 7. Solitary waves on sea-ice of constant thickness $d^*=0.1$ m with different wave celerity: (a) $c=1.001\sqrt {H}$, (b) $c=1.1\sqrt {H}$.

Elastic effects are introduced by the fifth-order term in the equivalent KdV equation that generates the oscillatory tail. As demonstrated in Benilov, Grimshaw & Kuznetsova (1993) and Grimshaw & Joshi (1995), the oscillatory tail of a generalized solitary solution is one-sided in a realistic wave system, otherwise a pair of source and sink must be required for energy conservation. The oscillatory tail propagates either in front of or behind the core soliton, depending on the sign of the coefficients in the equivalent KdV equation, and radiates energy from the core. Moreover, a small value of $\beta$ leads to a large wavenumber $k_o$ and results in an exponentially small amplitude of the oscillation with order $\exp [({{-{\rm \pi} c\mu ^2 H(3\alpha +H)}/{3})\sqrt {{1}/{\beta (c^2-H)}}}]$. To evaluate the effects of flexural stiffness $\beta$, solitary waves are simulated on sea-ice with different thicknesses. The numerical results conform to the theoretical predictions. Figure 8 shows that the oscillatory tail appears only in front of the core soliton and has larger wavelength and amplitude for a higher flexural stiffness, resulting in a quicker decay of the core soliton. The simulation of a solitary wave over sea-ice with varying thickness shows similar results (figure 9).

Figure 8. Solitary waves on sea-ices with different thickness (oscillatory tails at $t^*=900$ s are shown in the insets), for $c=1.001\sqrt {H}$ and (a) $d^*=0.6$ m, (b) $d^*=1.2$ m, (c) $d^*=1.8$ m.

Figure 9. Solitary wave propagating over sea-ice with varying thickness (0.1 m to 3 m).

To demonstrate the shoaling effect, a solitary wave propagating in uniform sea-ice ($d^*=0.1$ m) over varying bathymetry is simulated. Figure 10 shows the evolution of the solitary wave when water depth reduces from 10 to 2 m in the range from $x^*=3000$ to $x^*=6000$ m. The speed of the solitary wave reduces, but its amplitude and steepness increase with decreasing water depth. It is expected that the solitary wave would break if it propagates to even shallower water.

Figure 10. Solitary wave on uniform sea-ice over varying water depth ($c=1.001\sqrt {H}$, $d^*=0.1$ m, and $H^*$ changes from $-$10 to $-$2 m).

It should be noted that the solitary solution given by (3.27) is derived assuming that the equivalent KdV equation is singularly perturbed on the fifth-order term, and the KdV equation itself is a reduced model of the proposed Boussinesq equations. Therefore, (3.27) is an approximate solitary solution for FG waves in shallow water with limitations. Exact quantitative analysis of solitary solutions for FG waves in shallow water requires further investigation on the Boussinesq-type equations.

4. Analyses of SH waves in shallow water

Analytical solutions to the first-order gravity–hydroelastic coupled waves with a submerged plate in shallow water are derived, resulting in newly discovered dual-mode solutions. The resulting dispersion relation for SH waves is evaluated. Co-existence of the dual-mode solutions is demonstrated in numerical simulations.

4.1. Dual-mode analytical solutions

The first-order linear governing equations for a submerged plate in shallow water are derived for constant water depth and uniform plate stiffness from (2.28) and (2.29):

(4.1)$$\begin{gather} \frac{{\partial}^2 \hat{\eta}_1}{{\partial} t^2} - \frac{{\partial}^2 \zeta_1}{{\partial} t^2} - H_I\, \frac{{\partial}^2 \hat{\eta}_1}{{\partial} x^2} = 0, \end{gather}$$

(4.2)$$\begin{gather}\frac{{\partial}^2 \zeta_1}{{\partial} t^2} - H_{II}\,\frac{{\partial}^2 \hat{\eta}_1}{{\partial} x^2} - H_{II}\beta\, \frac{{\partial}^6 \zeta_1}{{\partial} x^6} = 0. \end{gather}$$

We solve these equations analytically using separation of variables and periodic conditions by assuming that surface gravity wave $\hat {\eta }_1$ and SH wave $\zeta _1$ have the same periodic characteristics in both space and time, such that $\zeta _1(x,t)=C_3\,\hat {\eta }_1(x,t)=C_3\,\hat {A}(x)\,\hat {B}(t)$. Then (4.1) and (4.2) become

(4.3)$$\begin{gather} \frac{1}{\hat{B}(t)}\,\frac{{\rm d}^2 \hat{B}(t)}{{\rm d} t^2} -\frac{C_3}{\hat{B}(t)}\,\frac{{\rm d}^2 \hat{B}(t)}{{\rm d} t^2} = \frac{H_I}{\hat{A}(x)}\,\frac{{\rm d}^2 \hat{A}(x)}{{{\rm d}x}^2} = C_1, \end{gather}$$

(4.4)$$\begin{gather}\frac{1}{\hat{B}(t)}\,\frac{{\rm d}^2 \hat{B}(t)}{{\rm d} t^2} = \frac{H_{II}}{C_3\hat{A}(x)}\,\frac{{\rm d}^2 \hat{A}(x)}{{{\rm d}x}^2} + \frac{H_{II}\beta}{\hat{A}(x)}\,\frac{{\rm d}^6 \hat{A}(x)}{{{\rm d}x}^6} = C_2, \end{gather}$$

where $C_1$, $C_2$ and $C_3$ are real constants.

Substituting (4.4) into (4.3), we derive a sixth-order homogeneous ODE with respect to $\hat {A}(x)$, whose characteristic equation can be reduced to a cubic equation. Here, $\hat {A}(x)$ is solved by computing the cubic characteristic equation, and the resulting solution for $\hat {A}(x)$ is substituted back into (4.4) to solve $\hat {B}(t)$. Periodic solutions exist in only two scenarios: ${H_{II}}/{(H_I+H_{II})}< C_3<1$ and $C_3<0$. The general periodic solution of $\hat {\eta }_1$ is derived as

(4.5)\begin{equation} \hat{\eta}_1(x,t)=C_4\exp({{\rm i}(\hat{k}x \pm \hat{\sigma} t)}), \end{equation}

in which $\hat {\sigma }$ and $\hat {k}$ are the non-dimensional angular frequency and wavenumber, respectively, and $C_4$ is a real constant determined by boundary and initial conditions.

The SH wave along the submerged plate shares the same frequency and wavenumber as the coupled surface wave as $\zeta _1(x,t)=C_3\,\hat {\eta }_1(x,t)$. Because $C_3$ can be either positive or negative, the coupled hydroelastic wave system has dual-mode analytical periodic solutions: when $C_3$ is positive and always less than $1$, the SH wave has the same phase but a smaller amplitude compared to the coupled gravity wave on the free surface; when $C_3$ is negative, the SH wave is 180$^{\circ }$ out of phase with the surface gravity wave.

In the derivation of analytical solutions, we have

(4.6)\begin{equation} C_3=1-\frac{\hat{k}^2}{\hat{\sigma}^2}\,H_I, \end{equation}

and the dispersion relation in an implicit form is

(4.7)\begin{equation} \hat{k}^2=\sqrt{\frac{\dfrac{\hat{\sigma}^2}{\hat{k}^2}-H_I-H_{II}}{H_{II} \beta-\dfrac{\hat{k}^2}{\hat{\sigma}^2}\,H_IH_{II}\beta}}. \end{equation}

Equation (4.7) can be further reformulated into a quartic equation with respect to $\lambda =\hat {k}^2$:

(4.8)\begin{equation} \frac{H_I H_{II} \beta}{\hat{\sigma}^2}\,\lambda^4 - H_{II}\beta\lambda^3 - (H_I+H_{II})\lambda + \hat{\sigma}^2=0. \end{equation}

Solving the quartic equation with Ferrari's method, (4.8) has two positive real roots and a pair of complex conjugate roots, given that $H_I$, $H_{II}$, $\beta$ and $\hat {\sigma }$ are positive real constants. The wavenumbers of the periodic solutions are the square roots of the two positive real roots of $\lambda$, corresponding to ${H_{II}}/{(H_I+H_{II})}< C_3<1$ and $C_3<0$, respectively. Considering (4.6), the smaller root of $\lambda$ gives the wavenumber of the in-phase mode with larger wavelength and celerity $({L}/{T})^2>H_I+H_{II}$, and the greater root of $\lambda$ yields the wavenumber of the out-of-phase mode with smaller wavelength and celerity $({L}/{T})^2< H_I$. Hereafter, we employ superscripts to denote the wavenumber and $C_3$: $\hat {k}^+$ and $C_3^+$ are associated with the in-phase mode, and $\hat {k}^-$ and $C_3^-$ are associated with the out-of-phase mode.

To demonstrate the dual-mode solutions and the effects of submerged depth, we compute the gravity–hydroelastic coupled waves on a submerged polyvinyl chloride (PVC) plastic plate, using the analytical solutions. Table 3 presents the physical parameters used for the computations. Computing and substituting the corresponding dimensionless parameters ($H_I$, $H_{II}$, $\beta$ and $\hat {\sigma }$) into (4.8) to solve for $\lambda$ and $\hat {k}$, the dual-mode solutions of $\hat {\eta }_1$ and $\zeta _1$ are obtained using (4.5) and (4.6). Figures 11 and 12 show the in-phase and out-of-phase coupled waves computed for various submerged depths, respectively.

Table 3. Physical parameters for submerged PVC plate computation.

Figure 11. Analytical solutions (in-phase mode) with a PVC plate submerged at different depths (wave amplitudes are exaggerated by a factor of 10). Submerged depths are (a) $1$ m, (b) $2$ m, (c) $3$ m.

Figure 12. Analytical solutions (out-of-phase mode) with a PVC plate submerged at different depths (wave amplitudes are exaggerated by a factor of 10). Submerged depths are (a) $1$ m, (b) $2$ m, (c) $3$ m.

For the in-phase mode (figure 11), the amplitude of the SH wave decreases with increasing submerged depth, while wavelength and celerity of the coupled waves decrease with increasing submerged depth. For the out-of-phase mode (figure 12), the wavelength, wave celerity and amplitude of the SH wave all increase with increasing submerged depth. Amplitude of the SH wave can be even greater than that of the surface wave in the out-of-phase mode when $C_3<-1$ (see figure 12c).

4.2. First-order linear dispersion relation

To further investigate the characteristics of the dual-mode solutions, we rewrite (4.7) and (4.8) for the first-order linear dispersion relation:

(4.9)\begin{equation} \frac{H_I H_{II} \beta}{\hat{\sigma}^2}\,\hat{k}^8 - H_{II}\beta \hat{k}^6 - (H_I+H_{II})\hat{k}^2 + \hat{\sigma}^2=0. \end{equation}

Equation (4.9) reduces to (3.2), the first-order dispersion relation of FG waves, when the submerged depth $H_I$ goes to zero. Equation (4.9) also reduces to the first-order dispersion relation of shallow-water waves when $H_{II}$ goes to zero ($H_I\rightarrow H$) or $\beta$ becomes zero. The linear dispersion relation for SH waves in water of finite depth was derived by Hassan et al. (2005), which can be non-dimensionalized as

(4.10)\begin{align} & \mu \hat{\sigma}^2 \beta \hat{k}^5 - \beta \hat{k}^6\tanh(\mu \hat{k}H_I) - \mu^3 \alpha \hat{\sigma}^4\hat{k} + \mu^2\alpha\hat{\sigma}^2\hat{k}^2\tanh(\mu \hat{k}H_I) \nonumber\\ &\quad =\mu^2\hat{\sigma}^4\tanh(\mu \hat{k}H_I)-\mu\hat{\sigma}^2\hat{k} - \frac{\mu\hat{\sigma}^2\hat{k}\tanh(\mu \hat{k}H_I)-\mu^2\hat{\sigma}^4}{\tanh(\mu \hat{k}H_{II})}. \end{align}

Equation (4.10) reduces to (4.9) for shallow water, when $\mu$ is small ($\tanh \mu \hat {k}H\rightarrow \mu \hat {k}H$) and its higher-order terms are neglected.

The effects of submerged depth on SH waves are investigated for various wave frequencies and total depths. The dispersion relation (4.9) is solved for $\hat {k}^+$ and $\hat {k}^-$ for different submerged depths $H_I$, and the resulting $C_3^+$ and $C_3^-$ are computed as well. Figures 13 and 14 show that the wavelength and $C_3^+$ decrease with increasing submerged depth for the in-phase mode. For the out-of-phase mode, on the other hand, its wavelength generally increases with increasing submerged depth. Such results can be explained by combining (2.28) and (2.29), which yields

(4.11)\begin{equation} \frac{{\partial}^2 \hat{\eta}_1}{{\partial} t^2} - \boldsymbol{\nabla}\boldsymbol{\cdot}(H\,\boldsymbol{\nabla}\hat{\eta}_1) - \boldsymbol{\nabla}\boldsymbol{\cdot}\{H_{II}\,\boldsymbol{\nabla}[\nabla^2(\beta\,\nabla^2\zeta_1)]\} =0. \end{equation}

The first two terms in (4.11) are similar to those in the linear FG wave equation (2.17), representing the gravitational effects of surface waves, and the third term represents the elastic effects of the submerged plate. Comparing (4.11) with (2.17), it is noted that the effective depth of the elastic term becomes the depth below the submerged plate ($H_{II}$), such that the wavelength decreases with decreasing $H_{II}$ (increasing $H_I$) for the in-phase mode. Substituting the out-of-phase mode $\zeta _1=-|C_3|\,\hat {\eta }_1$ into the third term of (4.11), the sign of the elastic term flips comparing to the in-phase mode, resulting in the opposite influence of submerged depth on elastic effects.

Figure 13. Dispersion relations of SH waves in shallow water with various wave periods ($L_0=50$ m, $H_0 =4\ {\rm m}$, $\beta =3.155\times 10^{-4}$, PVC thickness 0.4 m).

Figure 14. Dispersion relations of SH waves in shallow water with various water depths ($L_0=50$ m, $T^{*}=5\ {\rm s}$, $\beta =3.155\times 10^{-4}$, PVC thickness 0.4 m).

The effects of flexural stiffness on the first-order dispersion relation are investigated similarly. Figure 15 shows that increasing plate stiffness results in increase of wavelength for both in-phase and out-of-phase modes. When $\beta \rightarrow 0$, the elastic effects fade away and the in-phase mode converges to the normal shallow-water wave, such that $\hat {k}^+$ becomes independent of submerged depth, and $C_3^+$ presents motion magnitude of a submerged water particle relative to the surface. When $\beta \rightarrow \infty$, the submerged plate plays a role similar to the rigid bottom boundary, and the out-of-phase mode converges to the normal shallow-water wave over the rigid plate, such that $\hat {k}^-$ approaches the wavenumber of the shallow-water wave, and $C_3^-$ reaches zero.

Figure 15. Dispersion relations of SH waves in shallow water with various flexural stiffnesses ($L_0=50\ {\rm m}$, $H_0=4$ m, $T^*=5$ s).

4.3. Dual-mode numerical solutions

We develop a numerical solver to compute the first-order gravity–hydroelastic coupled wave system, which allows for demonstration of the dual-mode solutions with varying flexural stiffness and varying submerged depth.

The equations preceding (2.28) and (2.29), other than themselves, are solved to reduce the order of time derivatives:

(4.12)$$\begin{gather} \frac{{\partial} \hat{\eta}_1}{{\partial} t} - \frac{{\partial} \zeta_1}{{\partial} t} =- \boldsymbol{\nabla}\boldsymbol{\cdot}(H_I\bar{\boldsymbol{u}}_{I1}), \end{gather}$$

(4.13)$$\begin{gather}\frac{{\partial} \zeta_1}{{\partial} t} =- \boldsymbol{\nabla}\boldsymbol{\cdot}(H_{II}\bar{\boldsymbol{u}}_{II1}), \end{gather}$$

(4.14)$$\begin{gather}\frac{{\partial} \bar{\boldsymbol{u}}_{I1}}{{\partial} t} =- \boldsymbol{\nabla}\hat{\eta}_1, \end{gather}$$

(4.15)$$\begin{gather}\frac{{\partial} \bar{\boldsymbol{u}}_{II1}}{{\partial} t} =- \boldsymbol{\nabla}\hat{\eta}_1 - \boldsymbol{\nabla}[\nabla^2(\beta\,\nabla^2\zeta_1)]\}, \end{gather}$$

in which $\bar {\boldsymbol {u}}_{I1}$ and $\bar {\boldsymbol {u}}_{II1}$ are the first-order horizontal velocities above and below the submerged plate, respectively. The use of provisional values allows conditionally stable computation using the central difference schemes for spatial derivatives. Therefore, provisional steps are used in an explicit time-marching scheme to solve (4.12)–(4.15):

(4.16)$$\begin{gather} \zeta_{1,s+1^*} = \zeta_{1,s} - \Delta t\,\boldsymbol{\nabla}\boldsymbol{\cdot}(H_{II}\bar{\boldsymbol{u}}_{II1,s}), \end{gather}$$

(4.17)$$\begin{gather}\hat{\eta}_{1,s+1^*} = \hat{\eta}_{1,s} - \Delta t\,\boldsymbol{\nabla}\boldsymbol{\cdot}(H_{II}\bar{\boldsymbol{u}}_{II1,s}) - \Delta t\,\boldsymbol{\nabla}\boldsymbol{\cdot}(H_{I}\bar{\boldsymbol{u}}_{I1,s}), \end{gather}$$

(4.18)$$\begin{gather}\bar{\boldsymbol{u}}_{I1,s+1} = \bar{\boldsymbol{u}}_{I1,s} - \Delta t\,\boldsymbol{\nabla}\left(\frac{\hat{\eta}_{1,s}+\hat{\eta}_{1,s+1^*}}{2}\right)\!, \end{gather}$$

(4.19)$$\begin{gather}\bar{\boldsymbol{u}}_{II1,s+1} = \bar{\boldsymbol{u}}_{II1,s} - \Delta t\,\boldsymbol{\nabla}\left(\frac{\hat{\eta}_{1,s}+\hat{\eta}_{1,s+1^*}}{2}\right) - \Delta t\,\boldsymbol{\nabla}\left\{\nabla^2\left[\beta\,\nabla^2\left(\frac{\zeta_{1,s}+\zeta_{1,s+1^*}}{2}\right)\right]\right\}\!, \end{gather}$$

(4.20)$$\begin{gather}\zeta_{1,s+1} = \zeta_{1,s} - \Delta t\,\boldsymbol{\nabla}\boldsymbol{\cdot}(H_{II}\bar{\boldsymbol{u}}_{II1,s+1}), \end{gather}$$

(4.21)$$\begin{gather}\hat{\eta}_{1,s+1} = \hat{\eta}_{1,s} - \Delta t\,\boldsymbol{\nabla}\boldsymbol{\cdot}(H_{II}\bar{\boldsymbol{u}}_{II1,s+1}) - \Delta t\,\boldsymbol{\nabla}\boldsymbol{\cdot}(H_{I}\bar{\boldsymbol{u}}_{I1,s+1}), \end{gather}$$

in which $\hat {\eta }_{1,s+1^*}$ and $\zeta _{1,s+1^*}$ are the provisional values of $\hat {\eta }_{1,s+1}$ and $\zeta _{1,s+1}$, respectively. A linear numerical solver is then developed based on (4.16)–(4.21).

The linear gravity–hydroelastic coupled wave system is computed numerically to demonstrate the influence of varying flexural stiffness. A harmonically oscillating boundary excitation is applied to the left end on the free surface to induce a regular incident wave. Figure 16 presents the responses of the free surface and the PVC plates with different stiffnesses, in which the parameters in table 3 are employed for computations, except the thickness. Figure 17 shows the coupled wave responses with plate thickness varying from $d^*=0.1$ to $d^*=0.6$ m. It is clear that both the in-phase and out-of-phase modes co-exist. The surface wave does not change noticeably due to the varying stiffness, while the SH wave is significantly affected by the plate stiffness: the out-of-phase mode is more prominent in the plate with lower stiffness, which conforms to the observation in figure 15(d).

Figure 16. Linear gravity–hydroelastic coupled waves with different flexural stiffnesses ($H_0=4$ m, $H_I=0.5$, wave amplitudes are exaggerated by a factor of 10). (a) $\beta = 1.204 \times 10^{-5}\ (d^{\ast} = 0.1\ {\rm m})$; (b) $\beta = 3.250 \times 10^{-4} (d^{\ast} = 0.3\ {\rm m})$; (c) $\beta = 2.600 \times 10^{-3}\ (d^{\ast} = 0.6\ {\rm m})$.

Figure 17. Linear gravity–hydroelastic coupled waves with varying flexural stiffnesses ($H_0=4$ m, $H_I=0.5$, wave amplitudes are exaggerated by a factor of 10). (a) $\beta{:}\ 1.204 \times 10^{-5}\unicode{x2013}2.600 \times 10^{-3}$; (b) $\beta{:}\ 2.600 \times 10^{-3}\unicode{x2013}1.204 \times 10^{-5}$.

The influence of varying submerged depth on the coupled wave system is more complicated than that of varying stiffness. Figure 18 shows the responses of the free surface and the PVC plate submerged at different depths subject to the harmonic boundary excitation. Both the surface elevation and plate motion are sensitive to submerged depth: the dominant mode of the coupled wave system changes from the in-phase mode to the out-of-phase mode with increasing submerged depth. Figure 19 presents the coupled wave responses with submerged depth varying from 0.5 to 3.5 m in total water depth of 4 m. When submerged depth gradually increases from $H_I=0.125$ to $H_I=0.875$ (figure 19a), the in-phase mode is dominant: the amplitude of the SH wave decreases dramatically with increasing submerged depth, while the amplitude of the surface wave sightly increases, which conforms to figures 13(b) and 14(b). When submerged depth gradually decreases from $H_I=0.875$ to $H_I=0.125$ (figure 19b), the out-of-phase mode component in the surface wave significantly increases, while the in-phase mode component in the SH wave increases.

Figure 18. Linear gravity–hydroelastic coupled waves with various submerged depths ($d^*=0.2$ m, $\beta =9.629\times 10^{-5}$, wave amplitudes are exaggerated by a factor of 10). Submerged depths are (a) $1$ m, (b) $2$ m, (c) $3$ m.

Figure 19. Linear gravity-hydroelastic-coupled waves with varying submerged depth ($d^*=0.4$ m, $\beta =3.155\times 10^{-4}$, wave amplitudes are exaggerated by a factor of 10). (a) $H_I : 0.125\unicode{x2013}0.875$; (b) $H_I : 0.875\unicode{x2013}0.125$.

Energy flows into the coupled wave system from the boundary excitation and is distributed between the in-phase and out-of-phase modes. For each mode, the modal energy is further distributed between the free surface and the submerged plate based on the values of $C_3^+$ and $C_3^-$. The coupled wave responses are the combined results of these energy distributions. The result in figure 19(b) is appealing because we notice that the surface wave energy density remarkably increases due to the plate with decreasing submerged depth. The energy input from the oscillating boundary is allocated to the kinetic energies of the upper and lower flows, the gravitational potential at the free surface, and the elastic and kinetic energies of the submerged plate. The kinetic energy associated with the plate is neglected in the linear model, while the elastic potential is included. The free surface wave is associated with the gravitation potential and the kinetic energy of the upper flow, and similarly the SH wave is an indicator of the elastic potential of the plate and the kinetic energy of the lower flow. Figure 20 shows the horizontal velocities of the upper and lower flows. The changes in flow velocities follow the coupled waves in figure 19(b): the amplitude of the upper flow velocity increases, and its dominant mode becomes the out-of-phase mode as the submerged depth decreases. The average energy per unit surface area over a wave period can be computed based on the formula

(4.22)\begin{align} P(x^*,t^*) &= \frac{1}{2T^*}\,\rho \int_{t^*}^{t^*+T^*} \int_{-H_I^*+\zeta_1^*}^{\hat{\eta}_1^*}\bar{u}_{I1}^{*2} \,{\rm d} z^*\,{\rm d} t^* + \frac{1}{T^*}\,\rho \int_{t^*}^{t^*+T^*}\int_{0}^{\hat{\eta}_1^*} gz^*\,{\rm d} z^*\,{\rm d} t^* \nonumber\\ &\approx \frac{1}{2T^*}\,\rho H_I^* \int_{t^*}^{t^*+T^*} \bar{u}_{I1}^{*2} \,{\rm d} t^* + \frac{1}{2T^*}\,\rho g \int_{t^*}^{t^*+T^*} \hat{\eta}_1^{*2} \,{\rm d} t^*. \end{align}

Both the in-phase and out-of-phase modes share the same period as the boundary excitation ($T^*=10$ s in these computations). Substituting the time histories of $\bar {u}_{I1}^*$ and $\hat {\eta }_1^*$ into (4.22), we compute that the surface wave energy density at $x^*=600$ m is nearly nine times greater than that at $x^*=100$ m (7.86 J m$^{-2}$ versus 0.79 J m$^{-2}$). Therefore, the plate with decreasing submerged depth effectively extracts wave energy from the deeper water and transfers it to the surface layer. Such effect is enabled by energy redistributions between the in-phase and out-of-phase modes, and between the surface and SH waves. We have investigated the energy distribution between the free surface and SH waves based on the values of $C_3^+$ and $C_3^-$ for various submerged depth and flexural stiffness. However, the energy distribution between the in-phase and out-of-phase modes is not trivial. That depends on total water depth, submerged depth, properties of the submerged plate, and the means of energy input, such as frequency and location of the excitation. Analysis on the dual-mode solutions of SH waves will be within the scope of our future study.

Figure 20. Averaged horizontal velocities of the upper and lower flows with varying submerged depths (3.5–0.5 m).

5. Conclusions

Boussinesq approximations and Kirchhoff–Love plate theory are applied to study hydroelastic waves in shallow water, in which variable stiffness, varying submerged depth, and varying bathymetry are considered rigorously. The first-order linear and the nonlinear Boussinesq-type governing equations are derived for both FG waves of a floating plate and gravity–hydroelastic coupled waves with a submerged plate. Analytical periodic solutions to the linear equations are derived, and the linear dispersion relations are obtained based on the analytical solutions. These developments in this work provide two new theoretical models for hydroelastic waves in shallow water. Both linear and nonlinear numerical solvers are developed using the finite-difference spatial discretization and provisional time-marching scheme, to investigate the behaviours of hydroelastic waves.

For periodic FG waves of a floating plate, we found that the wave velocities do not necessarily increase with increasing wave period, and the wave amplitude does not necessarily increase during shoaling. The wavelength and celerity of nonlinear FG waves have been demonstrated to be smaller than those of linear counterparts. Variable stiffness and varying bathymetry have significant effects on wavelength and celerity, which should be considered in the design of VLFS. A small-amplitude generalized solitary wave solution is derived and computed numerically, which consists of a ‘$\textrm {sech}^2$’ core soliton and a one-sided oscillatory tail in front of the soliton. The oscillatory tail is found to be sensitive to flexural stiffness, while the core soliton is greatly affected by water depth.

For periodic gravity–hydroelastic coupled waves with a submerged plate, dual-mode analytical solutions are derived, and the numerical solutions have demonstrated that the in-phase and out-of-phase modes co-exist. Changing flexural stiffness or varying submerged depth not only significantly alter the wavelength and celerity of both modes, but also change the energy distributions between the two modes as well as between the free surface and the submerged plate. A submerged plate with decreasing submerged depth has been shown to enable wave energy transmission from the deeper water to the surface layer, which can be employed for wave energy extraction. These findings provide valuable insights to future development of offshore structures for wave focusing, trapping and cloaking.