# Nonlinear Hydroelastic Waves beneath a Floating Ice Sheet in a Fluid of Finite Depth

Authors: Ping Wang; Zunshui Cheng
Journal: Abstract and Applied Analysis
 (2013)

---

## Abstract

The nonlinear hydroelastic waves propagating beneath an infinite ice sheet floating on an inviscid fluid of finite depth are investigated analytically. The approximate series solutions for the velocity potential and the wave surface elevation are derived, respectively, by an analytic approximation technique named homotopy analysis method (HAM) and are presented for the second-order components. Also, homotopy squared residual technique is employed to guarantee the convergence of the series solutions. The present formulas, different from the perturbation solutions, are highly accurate and uniformly valid without assuming that these nonlinear partial differential equations (PDEs) have small parameters necessarily. It is noted that the effects of water depth, the ice sheet thickness, and Young’s modulus are analytically expressed in detail. We find that, in different water depths, the hydroelastic waves traveling beneath the thickest ice sheet always contain the largest wave energy. While with an increasing thickness of the sheet, the wave elevation tends to be smoothened at the crest and be sharpened at the trough. The larger Young’s modulus of the sheet also causes analogous effects. The results obtained show that the thickness and Young’s modulus of the floating ice sheet all greatly affect the wave energy and wave profile in different water depths.


---

## Body

## 1. Introduction

In recent decades, the ice cover in the polar region has attracted more and more attention in the field of ocean engineering and polar engineering in view of their practical importance and theoretical investigations. The motivations for the research work are to study damage to offshore constructions by floating ice sheets, the transportation systems in the cold region where the ice cover can be considered as roads and aircraft runways and air-cushioned vehicles are used to break the ice, for example. One of the important problems in this field would appear to be the accurate measurement of the characteristics of waves traveling beneath a floating ice sheet. And such wave may have been generated in the ice cover itself by the wind, or it may have originated by a moving load on the ice sheets. Considerable work has been done since the first theoretical model of wave propagation in sea ice was proposed by Greenhill [1] in 1887. A comprehensive summary on mathematical method and modeling for the problem can be found in some review articles such as Squire et al. [2, 3]. In addition to ice sheets, this work can apply to very large floating structures (VLFSs) such as floating airports, mobile offshore bases, offshore port facilities, offshore storage and waste disposal provisions, energy islands including some wave power configurations, and ultralarge ships, where there is an extensive complementary literature [4–6].Most theoretical works on the problem are still in the scope of linear theory based on the assumption that the wave amplitudes generated are very small in comparison with the wave lengths. So such models are not appropriate to describe waves of arbitrary amplitude considered here. According to hydrodynamics and elasticity, we can construct the nonlinear partial differential equations (PDEs) (1)–(5) to describe nonlinear hydroelastic waves of arbitrary amplitude traveling through water covered by an ice sheet in finite water depth. Unfortunately, it is very difficult to solve analytically the coupled nonlinear PDEs mathematically. Further, most of the most works literature on the nonlinear theory of sea waves ice sheet interaction are necessarily in the context of weakly nonlinear analysis due to the limitation of present mathematical tools. Now the main analytical study on such complex nonlinear PDEs still follows the well-known perturbation technique. For example, Forbes [7] derived nonlinear PDEs to describe two-dimensional periodic waves beneath an elastic sheet floating on the surface of an infinitely deep fluid. The periodic solutions are sought using the Fourier series and perturbation expansions for the Fourier coefficients. And it is found that the solutions have certain features in common with capillary-gravity waves. Following the framework in [7], Forbes continued his study of finite-amplitude surface waves beneath a floating elastic sheet in infinitely deep water [8], and optimized their previous perturbation technology directly by developing the Fourier coefficients as expansions in the wave height. Waves of extremely large amplitude are found to exist, and results are presented for waves belonging to several different nonlinear solution branches. Recently, Vanden-Broeck and Părău [9] further extended the results of Forbes for periodic waves to the arbitrary-amplitude waves. It is noted that perturbation and asymptotic approximations of nonlinear PDEs often break down as nonlinearity becomes strong. So the weakly nonlinear solutions of small-amplitude waves are derived by the perturbation approach, while fully nonlinear solutions of large-amplitude waves have to be calculated numerically by means of the numerical series truncation method in Vanden-Broeck’s study.Furthermore, perturbation and asymptotic techniques depend extremely on the small/large parameters in general, while our nonlinear PDEs have no any small/large parameters. Thus the perturbation techniques are not applicable to the nonlinear problem under consideration. In this paper, we apply a new analytic approximation method known as the homotopy analysis method (HAM) to effectively solve the nonlinear PDEs presented here. Based on the concept of homotopy in algebraic topology, the HAM was proposed by Liao [10] in 1992. Unlike the perturbation method, the HAM is entirely independent of any small/large parameters. Moreover, it provides us with extremely large freedom to choose base functions and initial approximations (16) and (17) of solutions and auxiliary linear operators (21)–(23) only under some basic rules [11, 12]. More importantly, in contrast to all other previous analytic techniques, the HAM provides us a convenient way to control and adjust the convergence of the approximate series solutions by means of introducing an auxiliary parameter c0. The method has been systematically described by Liao [11, 12]. Recently the HAM has been successfully applied to the study of a number of classical nonlinear differential equations including nonlinear equations arising in fluid mechanics [13–18], heat transfer [19, 20], solitons and integrable models [21–24], and finance [25, 26]. These aforementioned studies show the validity and generality of the HAM for some highly nonlinear PDEs with multiple solutions, singularity, and unknown boundary conditions.The objective of the present work is to analytically study the nonlinear hydroelastic waves under an ice sheet lying over an incompressible inviscid fluid of finite uniform depth by means of the HAM. According to the potential theory in hydrodynamics and elasticity, the nonlinear partial differential equations (PDEs) (1)–(5) are composed of the Laplace equation taken as the governing equation for inviscid flows, the kinematic and dynamic boundary conditions on the unknown ice sheet-water interface with a zero draft, a simple linear model for the thin sheet that includes the effects of flexural rigidity and vertical inertia, and a bottom boundary condition. The convergent homotopy-series solutions for the velocity potential and the wave surface elevation are formally derived by applying the HAM with the consideration of the minimum of the squared residual, respectively. It should be mentioned that we study the effects of the water depth and two important physical parameters including Young’s modulus and the thickness of the ice sheet on the wave energy and its elevation in detail. Discussion and conclusions are made in Sections 4 and 5, respectively. All of results obtained will help enrich our understanding of nonlinear hydroelastic waves propagating under a floating ice sheet on a fluid of finite depth.


## 2. Mathematical Description

The problem under consideration is a train of nonlinear hydroelastic waves propagating beneath a two-dimensional infinite elastic plate floating on a fluid of finite depthh and a zero draft. A Cartesian coordinate oxz is used in which the z-axis points vertically upward, while z=0 represents the undisturbed surface. We follow Greenhill in [1] assuming that this problem is capable of modeling ocean waves in the presence of sea ice when the fluid is inviscid and incompressible and the flow is irrotational, and the ice sheet is mathematically idealized as a thin elastic plate. Then the governing equations for a velocity potential ϕ(x,z,t) can be written as
(1)∂2ϕ∂x2+∂2ϕ∂z2=0,(-h≤z≤ζ(x,t)),
where ζ(x,t) is the wave surface elevation. The bottom boundary condition reads
(2)∂ϕ∂z=0,(z=-h).The motion of the fluid and the plate is coupled through the dynamic free-surface condition. We also assume that any particle which is once between the elastic plate and the water surface remains on it. So the kinematic and dynamic boundary conditions on the unknown surfacez=ζ(x,t) are, respectively, modeled as
(3)∂ζ∂t+∂ϕ∂x∂ζ∂x-∂ϕ∂z=0,(4)∂ϕ∂t+12|∇ϕ|2+peρ+gζ=0,
where pe is the water-plate interface pressure, ρ is the fluid density, and g is the gravitational acceleration, for a thin homogeneous elastic plate with uniform mass density ρe and constant thickness d.Since we are considering long waves here, the linear Kirchhoff (Euler-Bernoulli) beam theory is applied to the floating elastic plate as follows:(5)pe=D∂4ζ∂x4+me(∂2ζ∂t2+g),
where me=ρed, D=Ed3/[12(1-ν2)] is the flexural rigidity of the plate, E is the effective Young’s modulus of the plate, and ν Poisson’s ratio. We substitute (5) into (4) to derive a new form of the dynamic boundary condition as follows:
(6)∂ϕ∂t+12|∇ϕ|2+gζ+1ρ[D∂4ζ∂x4+me(∂2ζ∂t2+g)]=0.Here, we consider a train of nonlinear waves traveling beneath an elastic plate with constant wave numberk and constant angular frequency ω of the incident wave. For a general case it should be emphasized that, by means of the traveling-wave method directly, the progressive waves are transferred from the temporal differentiation into the spatial one, which is very different from the mathematical model obtained by simply eliminating the time-dependent terms from the kinematic and dynamic boundary conditions on the unknown free surface [7–9]. Namely, we introduce an independent variable transformation
(7)X=kx-ωt,
where the angular frequency ω and the wave number k are given. Thus, we can express the potential function ϕ(x,z,t)=ϕ(X,z) and the traveling wave profile ζ(x,t)=ζ(X).Then the governing equation and the bottom boundary condition for the velocity potential are transformed, respectively, by(8)k2∂2ϕ∂X2+∂2ϕ∂z2=0,(-h≤z≤ζ(X)),(9)∂ϕ∂z=0,(z=-h).
With the transformation (7), (3), and (6) on z=ζ(X) are given by
(10)-ωdζdX+k2∂ϕ∂XdζdX-∂ϕ∂z=0,(11)-ω∂ϕ∂X+f+gζ+1ρ[Dk4d4ζdX4+me(ω2d2ζdX2+g)]=0,
respectively, where
(12)f=12[k2(∂ϕ∂X)2+(∂ϕ∂z)2].
We combine partially (10) and (11) to gain the boundary conditions on z=ζ(X) as follows:
(13)ω2∂2ϕ∂X2+g∂ϕ∂z-ω∂f∂X-ωρ(Dk4d5ζdX5+meω2d3ζdX3)-k2g∂ϕ∂XdζdX=0.
Now the corresponding unknown potential function ϕ(X,z) and the wave surface elevation ζ(X) are governed by (8), (9), (11), and (13).


## 3. Analytic Approach Based on the Homotopy Analysis Method

### 3.1. Solution Expression and Initial Approximation

Using the homotopy analysis method, we should first of all start from a set of base functions and solution expression which are very important to approximate the unknown solutions of the nonlinear boundary problem under consideration. Mathematically, it seems impossible to guess the expression forms of the unknown potential function and the wave vertical displacement. Fortunately, considering the physical background of our problem, we may gain proper solution expressions of it. From viewpoints of the physical considerations here, our problem is composed of a train of progressive waves cause by a load moving on the ice sheet, an infinite elastic plate acting as an ice sheet floating on an fluid of finite depth. As is well known, in case of the pure water waves, the progressive wave elevation can be expressed as(14)ζ(X)=∑n=0+∞‍βncos(nX)
by a set of base functions {cos(nX),n≥0}, where βn are unknown coefficients. In the case of plate-covered surface, since we assume that there is no gap between the bottom surface of the thin elastic plate and the top surface of the fluid layer and a zero draft, the vertical displacement of the thin plate is still periodic in the X direction. Therefore, we clearly know that ζ(X) can be expressed in the above form (14) too.Besides, according to the linear wave theory, we can find the solutions of the Laplace equation (8) by the separation of variables method. To acquire those solutions, we have to use kinematic and dynamic boundary conditions of the free surface and the boundary condition in finite water depth, and we consider the solution derived here as the solution expression of potential function
(15)ϕ(X,z)=∑n=1+∞‍αncosh[nk(z+h)]cosh(nkh)sin(nX)
by a set of base functions {cosh[nk(z+h)]/cosh(nkh)sin(nX),n≥0}, where αn are unknown coefficients. Note that the potential function ϕ(X,z) defined by (15) automatically satisfies the governing equation (8) and the bottom boundary condition (9). The above expressions (14) and (15) are called the solution expressions of ϕ(X,z) and ζ(X), respectively, which play important roles in the method of homotopy analysis.According to the solution expression (15) and the boundary condition (9), we construct the initial approximation of the potential function:
(16)ϕ0(X,z)=α0,1cosh[k(z+h)]cosh(kh)sin(X),
where α0,1 is an unknown coefficient. We choose
(17)ζ0(X)=0.
as the initial approximation of wave profile ζ(X) to simplify the subsequent solution procedure [18, 20]. It should be emphasized that higher order terms can hold the corrections of the analytic series solutions due to the nonlinearity inherent in (11) and (13) although the initial guess ζ0(X) is zero.


### 3.2. Continuous Variation

The HAM is based on a kind of continuous mapping of an initial approximation to the exact solution through a series of deformation equations. For simplicity, based on the nonlinear boundary condition (13) and (11), we define the two following nonlinear operators 𝒩1 and 𝒩2 as follows
(18)𝒩1[Φ(X,z;q),η(X;q)]=ω2∂2Φ(X,z;q)∂X2+g∂Φ(X,z;q)∂z-ω∂F∂X-ωρ(Dk4∂5η(X;q)∂X5+ω2me∂3η(X;q)∂X3)-k2g∂Φ(X,z;q)∂X∂η(X;q)∂X,(19)𝒩2[η(X;q),Φ(X,z;q)]=-ω∂Φ(X,z;q)∂X+F+gη(X;q)+1ρ[Dk4∂4η(X;q)∂X4+me(ω2∂2η(X;q)∂X2+g)],
where
(20)F=12[k2(∂Φ∂X)2+(∂Φ∂z)2]
and q∈[0,1] is the embedding parameter of the HAM.Here, it should be emphasized that, as mentioned by Liao and Cheung and Tao et al. [14, 15], the HAM provides us with extremely large freedom to choose the auxiliary linear operators and the initial guess. Note that both linear terms of Φ(X,z;q) and linear terms of η(X;q) are all contained in (18). If we choose all linear terms, the subsequent iterative procedure will become very complex. Fortunately, based on the HAM, we can completely forget the linear terms in (13) and choose proper auxiliary linear operator of Φ(X,z;q) by means of the solution expression (15) which is obtained under the physical considerations as
(21)ℒ-1[Φ(X,z;q)]=ω2∂2Φ(X,z;q)∂X2+g∂Φ(X,z;q)∂z.In particular, if the angular frequencyω is given, we can choose such an approximation based on the linear wave theory to simplify the subsequent resolution of the nonlinear PDEs as follows:
(22)ω≈gktanh(kh).
So we simplify the auxiliary linear operator in (21) as follows:
(23)ℒ1[Φ(X,z;q)]=gktanh(kh)∂2Φ(X,z;q)∂X2+g∂Φ(X,z;q)∂z,
where ℒ1[0]=0. Note that, due to the weakly nonlinear effects, the actual frequency ω is often slightly different from the linear dispersion relation ω0=gktanh(kh). In Section 4, ω/ω0=1.01 is chosen so that the perturbation theory is valid and corresponding results are highly accurate, and then we can compare our results with those obtained by the perturbation method.Based on the linear operator of the wave profile functionη(X;q) in the nonlinear operator 𝒩2, for simplicity, we may choose another auxiliary linear operator:
(24)ℒ2[η(X;q)]=∂4η(X;q)∂X4+∂2η(X;q)∂X2+η(X;q),
where ℒ2[0]=0.We letc0 be an nonzero convergence-control parameter. It is noted that both c0 and q in the HAM are auxiliary parameters without any physical meaning. Instead of the nonlinear PDEs (8), (9), (11), and (13), we reconstruct the so-called zeroth-order deformation equations as follows:
(25)k2∂2Φ(X,z;q)∂X2+∂2Φ(X,z;q)∂z2=0,(-h≤z≤η(X;q)),(26)∂Φ(X,z;q)∂z=0,(z=-h),(27)(1-q)ℒ1[Φ(X,z;q)-ϕ0(X,z)]=qc0𝒩1[Φ(X,z;q),η(X;q)],(z=η(X;q)),(28)(1-q)ℒ2[η(X;q)-ζ0(X)]=qc0𝒩2[η(X;q),Φ(X,z;q)],(z=η(X;q)).
Then, from (27) and (28), two mapping functions Φ(X,z;q) and η(X;q) vary respectively continuously from their initial approximation ϕ0(X,z) and ζ0(X) to the exact solutions ϕ(X,z) and ζ(X) of the original problem. The Taylor series of Φ(X,z;q) and η(X;q) at q=0 are
(29)Φ(X,z;q)=ϕ0(X,z)+∑m=1+∞‍ϕm(X,z)qm,(30)η(X;q)=ζ0(X)+∑m=1+∞‍ζm(X)qm,
where
(31){ϕm(X,z),ζm(X)}=1m!∂m∂qm{Φ(X,z;q),η(X;q)}|q=0.Assume thatc0 is so properly chosen that the series in (29) and (30) converge at q=1; then we have the so-called homotopy-series solutions as follows:
(32)ϕ(X,z)=Φ(X,z;1)=ϕ0(X,z)+∑m=1+∞‍ϕm(X,z),ζ(X)=η(X;1)=ζ0(X)+∑m=1+∞‍ζm(X).At thenth-order of approximations, we have
(33)ϕ(X,z)≈ϕ0(X,z)+∑m=1+n‍ϕm(X,z),ζ(X)≈ζ0(X)+∑m=1+n‍ζm(X).As shown later in the following section, the unknown termsϕm(X,z) and ζm(X) are governed by the linear PDEs (34)–(36).


### 3.3. High-Order Deformation Equations

High-order deformation equations for the unknownϕm(X,z), ζm(X) can be derived directly from the zeroth-order deformation equations. Firstly, substituting the homotopy-Maclaurin series (29) and (30) into the governing equation (25) and the boundary condition in finite water depth (26) and then equating the like-power of the embedding parameter q, we have
(34)k2∂2ϕm(X,z)∂X2+∂2ϕm(X,z)∂z2=0,(-h≤z≤0),∂ϕm(X,z)∂z=0,(z=-h),
where m≥1.Note that,Φ(X,z;q) at the unknown surface z=η(X;q) may be expressed in terms of the Taylor expansion at z=0 instead of z=η(X;q). The detailed derivation of the expansion of Φ(X,z;q) at the unknown surface is given in Appendices (A.1)–(A.5). Upon the substitution of appropriate series (A.5) and (30) into the boundary conditions (27) and (28), we have two linear boundary conditions on z=0 as follows:
(35)ℒ1(ϕm)|z=0=c0Δm-1ϕ+χmSm-1-S-m,(36)ℒ2(ζm)=c0Δm-1ζ+χm(d4ζm-1dX4+d2ζm-1dX2+ζm-1),
where
(37)χm={0,m⩽11,m>1.The detailed derivation of the above equations and the expression forϕm and ζm are given in Appendix A. It should be noted that (27) and (28) holds on the unknown boundary z=η(X;q), while (35) and (36) hold on z=0. Furthermore, the original nonlinear DPEs (1)–(5) are transferred into an infinite number of linear decoupled high-order deformation equations (34)–(36). Namely, given ϕm-1 and ζm-1, ϕm and ζm can be obtained easily by means of the inverse operators of the right-hand sides of (35) and (36), respectively, and a computer algebra system such as Mathematica. The resulting expressions for ϕm and ζm are presented to the second order in the coming subsection.


### 3.4. First-Order and Second-Order Approximations

Substituting initial approximations (16) and (17) into (36), we can get ζ1(X) using the inverse linear operator ℒ2 in (36) as follows:
(38)ζ1(X)=14[4dgc0+c0a0,12+k2c0a0,12tanh2(hk)]-ωc0α0,1cos(X)+152[c0a0,12-k2c0a0,12tanh2(hk)]cos(2X).But now the coefficientα0,1 in the initial approximation of ϕ0(X,z) in (16) is still unknown. So we introduce an additional equation to relate the solutions with the wave height:
(39)ζ1(mπ)-ζ1(nπ)=H,
in which m is an even integer, n is an odd integer, and H is the wave height to the first order based on the HAM. The relation (39) for the wave height and its vertical displacement can determine the solution of α0,1.Further, in the analogous manner as for the first-order approximation, by using the inverse linear operatorℒ1 in (35), it is easy to get the solution of ϕ1(X,z), especially by means of the symbolic computation software such as Mathematica:
(40)α0,1=-H2ωc0,ϕ1(X,z)=α1,1cosh[k(h+z)]cosh(kh)sin(X)+-H2+H2k2tanh2(hk)16gkωc0[2tanh(hk)-tanh(2hk)]×cosh[2k(h+z)]cosh(2kh)sin(2X).We find the common solutionϕ1(X,z) has one unknown coefficient α1,1 which can be determined by avoiding the “secular” term sin(X) in ϕ2(X,z). We note that all subsequent functions occur recursively. Utilizing the linear equations (35) and (36) to continue with the first-order approximations we have
(41)ζ2(X)=β2,0+β2,1cos(X)+β2,2cos(2X)+β2,3cos(3X)+β2,4cos(4X),ϕ2(X,z)=α2,1cosh[k(h+z)]cosh(kh)sin(X)+α2,2cosh[2k(h+z)]cosh(2kh)sin(2X)+α2,3cosh[3k(h+z)]cosh(3kh)sin(3X)+α2,4cosh[4k(h+z)]cosh(4kh)sin(4X)+α2,5cosh[5k(h+z)]cosh(5kh)sin(5X),
where αi,j is the jth unknown coefficient of ϕi(X,z) and βi,j is the jth unknown coefficient of ζi(X). The detailed expressions of these coefficients for ϕ2 and ζ2 are given in Appendix B.In order to obtain higher-order functionsϕm(X,z) and ζm(X), we need only to continue this approach. In principle, we can acquire infinite-order solutions for our physical model. It is also worthwhile to mention that these solutions will retain model parameters and the convergence control parameter c0.


### 3.5. Optimal Convergence-Control Parameter

If we fix all model parameters in our approximate series solutions, there is still an unknown convergence control parameterc0 in them, which is used to guarantee the convergence of our approximation solutions. According to Liao [12], it is the convergence control parameter c0 that essentially differs the HAM from all other analytic methods. And the optimal value of c0 is determined by the minimum of the total squared-residual εmT of our nonlinear problem, defined by
(42)εmT=εmϕ+εmζ,
where
(43)εmϕ=11+M∑i=0M‍(𝒩1[ϕ(X,z),ζ(X)]|X=iΔX)2,εmζ=11+M∑i=0M‍(𝒩2[ϕ(X,z),ζ(X)]|X=iΔX)2,
where εmϕ and εmζ are two residual square errors of boundary conditions (27) and (28), respectively. M is the number of the discrete points, and ΔX=π/M. In this paper, we choose M=10.Theorem 2.1 given by Liao in [12] can guarantee the rationality of (42). So we obtain the optimal convergence control parameter c0 by the minimum of the squared-residual εmT, generally corresponding to dεmT/dc0=0.


## 3.1. Solution Expression and Initial Approximation

Using the homotopy analysis method, we should first of all start from a set of base functions and solution expression which are very important to approximate the unknown solutions of the nonlinear boundary problem under consideration. Mathematically, it seems impossible to guess the expression forms of the unknown potential function and the wave vertical displacement. Fortunately, considering the physical background of our problem, we may gain proper solution expressions of it. From viewpoints of the physical considerations here, our problem is composed of a train of progressive waves cause by a load moving on the ice sheet, an infinite elastic plate acting as an ice sheet floating on an fluid of finite depth. As is well known, in case of the pure water waves, the progressive wave elevation can be expressed as(14)ζ(X)=∑n=0+∞‍βncos(nX)
by a set of base functions {cos(nX),n≥0}, where βn are unknown coefficients. In the case of plate-covered surface, since we assume that there is no gap between the bottom surface of the thin elastic plate and the top surface of the fluid layer and a zero draft, the vertical displacement of the thin plate is still periodic in the X direction. Therefore, we clearly know that ζ(X) can be expressed in the above form (14) too.Besides, according to the linear wave theory, we can find the solutions of the Laplace equation (8) by the separation of variables method. To acquire those solutions, we have to use kinematic and dynamic boundary conditions of the free surface and the boundary condition in finite water depth, and we consider the solution derived here as the solution expression of potential function
(15)ϕ(X,z)=∑n=1+∞‍αncosh[nk(z+h)]cosh(nkh)sin(nX)
by a set of base functions {cosh[nk(z+h)]/cosh(nkh)sin(nX),n≥0}, where αn are unknown coefficients. Note that the potential function ϕ(X,z) defined by (15) automatically satisfies the governing equation (8) and the bottom boundary condition (9). The above expressions (14) and (15) are called the solution expressions of ϕ(X,z) and ζ(X), respectively, which play important roles in the method of homotopy analysis.According to the solution expression (15) and the boundary condition (9), we construct the initial approximation of the potential function:
(16)ϕ0(X,z)=α0,1cosh[k(z+h)]cosh(kh)sin(X),
where α0,1 is an unknown coefficient. We choose
(17)ζ0(X)=0.
as the initial approximation of wave profile ζ(X) to simplify the subsequent solution procedure [18, 20]. It should be emphasized that higher order terms can hold the corrections of the analytic series solutions due to the nonlinearity inherent in (11) and (13) although the initial guess ζ0(X) is zero.


## 3.2. Continuous Variation

The HAM is based on a kind of continuous mapping of an initial approximation to the exact solution through a series of deformation equations. For simplicity, based on the nonlinear boundary condition (13) and (11), we define the two following nonlinear operators 𝒩1 and 𝒩2 as follows
(18)𝒩1[Φ(X,z;q),η(X;q)]=ω2∂2Φ(X,z;q)∂X2+g∂Φ(X,z;q)∂z-ω∂F∂X-ωρ(Dk4∂5η(X;q)∂X5+ω2me∂3η(X;q)∂X3)-k2g∂Φ(X,z;q)∂X∂η(X;q)∂X,(19)𝒩2[η(X;q),Φ(X,z;q)]=-ω∂Φ(X,z;q)∂X+F+gη(X;q)+1ρ[Dk4∂4η(X;q)∂X4+me(ω2∂2η(X;q)∂X2+g)],
where
(20)F=12[k2(∂Φ∂X)2+(∂Φ∂z)2]
and q∈[0,1] is the embedding parameter of the HAM.Here, it should be emphasized that, as mentioned by Liao and Cheung and Tao et al. [14, 15], the HAM provides us with extremely large freedom to choose the auxiliary linear operators and the initial guess. Note that both linear terms of Φ(X,z;q) and linear terms of η(X;q) are all contained in (18). If we choose all linear terms, the subsequent iterative procedure will become very complex. Fortunately, based on the HAM, we can completely forget the linear terms in (13) and choose proper auxiliary linear operator of Φ(X,z;q) by means of the solution expression (15) which is obtained under the physical considerations as
(21)ℒ-1[Φ(X,z;q)]=ω2∂2Φ(X,z;q)∂X2+g∂Φ(X,z;q)∂z.In particular, if the angular frequencyω is given, we can choose such an approximation based on the linear wave theory to simplify the subsequent resolution of the nonlinear PDEs as follows:
(22)ω≈gktanh(kh).
So we simplify the auxiliary linear operator in (21) as follows:
(23)ℒ1[Φ(X,z;q)]=gktanh(kh)∂2Φ(X,z;q)∂X2+g∂Φ(X,z;q)∂z,
where ℒ1[0]=0. Note that, due to the weakly nonlinear effects, the actual frequency ω is often slightly different from the linear dispersion relation ω0=gktanh(kh). In Section 4, ω/ω0=1.01 is chosen so that the perturbation theory is valid and corresponding results are highly accurate, and then we can compare our results with those obtained by the perturbation method.Based on the linear operator of the wave profile functionη(X;q) in the nonlinear operator 𝒩2, for simplicity, we may choose another auxiliary linear operator:
(24)ℒ2[η(X;q)]=∂4η(X;q)∂X4+∂2η(X;q)∂X2+η(X;q),
where ℒ2[0]=0.We letc0 be an nonzero convergence-control parameter. It is noted that both c0 and q in the HAM are auxiliary parameters without any physical meaning. Instead of the nonlinear PDEs (8), (9), (11), and (13), we reconstruct the so-called zeroth-order deformation equations as follows:
(25)k2∂2Φ(X,z;q)∂X2+∂2Φ(X,z;q)∂z2=0,(-h≤z≤η(X;q)),(26)∂Φ(X,z;q)∂z=0,(z=-h),(27)(1-q)ℒ1[Φ(X,z;q)-ϕ0(X,z)]=qc0𝒩1[Φ(X,z;q),η(X;q)],(z=η(X;q)),(28)(1-q)ℒ2[η(X;q)-ζ0(X)]=qc0𝒩2[η(X;q),Φ(X,z;q)],(z=η(X;q)).
Then, from (27) and (28), two mapping functions Φ(X,z;q) and η(X;q) vary respectively continuously from their initial approximation ϕ0(X,z) and ζ0(X) to the exact solutions ϕ(X,z) and ζ(X) of the original problem. The Taylor series of Φ(X,z;q) and η(X;q) at q=0 are
(29)Φ(X,z;q)=ϕ0(X,z)+∑m=1+∞‍ϕm(X,z)qm,(30)η(X;q)=ζ0(X)+∑m=1+∞‍ζm(X)qm,
where
(31){ϕm(X,z),ζm(X)}=1m!∂m∂qm{Φ(X,z;q),η(X;q)}|q=0.Assume thatc0 is so properly chosen that the series in (29) and (30) converge at q=1; then we have the so-called homotopy-series solutions as follows:
(32)ϕ(X,z)=Φ(X,z;1)=ϕ0(X,z)+∑m=1+∞‍ϕm(X,z),ζ(X)=η(X;1)=ζ0(X)+∑m=1+∞‍ζm(X).At thenth-order of approximations, we have
(33)ϕ(X,z)≈ϕ0(X,z)+∑m=1+n‍ϕm(X,z),ζ(X)≈ζ0(X)+∑m=1+n‍ζm(X).As shown later in the following section, the unknown termsϕm(X,z) and ζm(X) are governed by the linear PDEs (34)–(36).


## 3.3. High-Order Deformation Equations

High-order deformation equations for the unknownϕm(X,z), ζm(X) can be derived directly from the zeroth-order deformation equations. Firstly, substituting the homotopy-Maclaurin series (29) and (30) into the governing equation (25) and the boundary condition in finite water depth (26) and then equating the like-power of the embedding parameter q, we have
(34)k2∂2ϕm(X,z)∂X2+∂2ϕm(X,z)∂z2=0,(-h≤z≤0),∂ϕm(X,z)∂z=0,(z=-h),
where m≥1.Note that,Φ(X,z;q) at the unknown surface z=η(X;q) may be expressed in terms of the Taylor expansion at z=0 instead of z=η(X;q). The detailed derivation of the expansion of Φ(X,z;q) at the unknown surface is given in Appendices (A.1)–(A.5). Upon the substitution of appropriate series (A.5) and (30) into the boundary conditions (27) and (28), we have two linear boundary conditions on z=0 as follows:
(35)ℒ1(ϕm)|z=0=c0Δm-1ϕ+χmSm-1-S-m,(36)ℒ2(ζm)=c0Δm-1ζ+χm(d4ζm-1dX4+d2ζm-1dX2+ζm-1),
where
(37)χm={0,m⩽11,m>1.The detailed derivation of the above equations and the expression forϕm and ζm are given in Appendix A. It should be noted that (27) and (28) holds on the unknown boundary z=η(X;q), while (35) and (36) hold on z=0. Furthermore, the original nonlinear DPEs (1)–(5) are transferred into an infinite number of linear decoupled high-order deformation equations (34)–(36). Namely, given ϕm-1 and ζm-1, ϕm and ζm can be obtained easily by means of the inverse operators of the right-hand sides of (35) and (36), respectively, and a computer algebra system such as Mathematica. The resulting expressions for ϕm and ζm are presented to the second order in the coming subsection.


## 3.4. First-Order and Second-Order Approximations

Substituting initial approximations (16) and (17) into (36), we can get ζ1(X) using the inverse linear operator ℒ2 in (36) as follows:
(38)ζ1(X)=14[4dgc0+c0a0,12+k2c0a0,12tanh2(hk)]-ωc0α0,1cos(X)+152[c0a0,12-k2c0a0,12tanh2(hk)]cos(2X).But now the coefficientα0,1 in the initial approximation of ϕ0(X,z) in (16) is still unknown. So we introduce an additional equation to relate the solutions with the wave height:
(39)ζ1(mπ)-ζ1(nπ)=H,
in which m is an even integer, n is an odd integer, and H is the wave height to the first order based on the HAM. The relation (39) for the wave height and its vertical displacement can determine the solution of α0,1.Further, in the analogous manner as for the first-order approximation, by using the inverse linear operatorℒ1 in (35), it is easy to get the solution of ϕ1(X,z), especially by means of the symbolic computation software such as Mathematica:
(40)α0,1=-H2ωc0,ϕ1(X,z)=α1,1cosh[k(h+z)]cosh(kh)sin(X)+-H2+H2k2tanh2(hk)16gkωc0[2tanh(hk)-tanh(2hk)]×cosh[2k(h+z)]cosh(2kh)sin(2X).We find the common solutionϕ1(X,z) has one unknown coefficient α1,1 which can be determined by avoiding the “secular” term sin(X) in ϕ2(X,z). We note that all subsequent functions occur recursively. Utilizing the linear equations (35) and (36) to continue with the first-order approximations we have
(41)ζ2(X)=β2,0+β2,1cos(X)+β2,2cos(2X)+β2,3cos(3X)+β2,4cos(4X),ϕ2(X,z)=α2,1cosh[k(h+z)]cosh(kh)sin(X)+α2,2cosh[2k(h+z)]cosh(2kh)sin(2X)+α2,3cosh[3k(h+z)]cosh(3kh)sin(3X)+α2,4cosh[4k(h+z)]cosh(4kh)sin(4X)+α2,5cosh[5k(h+z)]cosh(5kh)sin(5X),
where αi,j is the jth unknown coefficient of ϕi(X,z) and βi,j is the jth unknown coefficient of ζi(X). The detailed expressions of these coefficients for ϕ2 and ζ2 are given in Appendix B.In order to obtain higher-order functionsϕm(X,z) and ζm(X), we need only to continue this approach. In principle, we can acquire infinite-order solutions for our physical model. It is also worthwhile to mention that these solutions will retain model parameters and the convergence control parameter c0.


## 3.5. Optimal Convergence-Control Parameter

If we fix all model parameters in our approximate series solutions, there is still an unknown convergence control parameterc0 in them, which is used to guarantee the convergence of our approximation solutions. According to Liao [12], it is the convergence control parameter c0 that essentially differs the HAM from all other analytic methods. And the optimal value of c0 is determined by the minimum of the total squared-residual εmT of our nonlinear problem, defined by
(42)εmT=εmϕ+εmζ,
where
(43)εmϕ=11+M∑i=0M‍(𝒩1[ϕ(X,z),ζ(X)]|X=iΔX)2,εmζ=11+M∑i=0M‍(𝒩2[ϕ(X,z),ζ(X)]|X=iΔX)2,
where εmϕ and εmζ are two residual square errors of boundary conditions (27) and (28), respectively. M is the number of the discrete points, and ΔX=π/M. In this paper, we choose M=10.Theorem 2.1 given by Liao in [12] can guarantee the rationality of (42). So we obtain the optimal convergence control parameter c0 by the minimum of the squared-residual εmT, generally corresponding to dεmT/dc0=0.


## 4. Results and Analysis

In order to show the convergence of the analytical series solution to our problems by means of the HAM, we consider the cases ofk=π/5m-1, d=0.01m, ρe=900kgm-3, ν=0.33, E=1010Nm-2, h=5m, H=0.1m, and ω/ω0=1.01 and take these data hereinafter for computation unless otherwise stated. The total residual square error εmT at several orders of approximation versus the convergence-control parameter c0 is shown in Figure 1. It is found that εmT at every order has the smallest values which corresponds to the optimal c0. For example, as m=7, the optimal c0=-0.18, and the smallest value of ε7T=7.910×10-6. So, let the optimal convergence-control parameter c0=-0.18, the total residual square error εmT decreases quickly as the order m increases, as shown in Table 1. It is also found that ε15T is down to 5.382×10-11 at the 15th-order of approximation, which indicates the convergence of our series solutions. In this way, we ensure that all our solutions are highly accurate.Table 1
The total residual square errorεmT for different approximation order m with c0=-0.18.
m
ε
m
T
1
3.497
×
1
0
-
3
3
3.404
×
1
0
-
4
5
3.700
×
1
0
-
5
7
7.910
×
1
0
-
6
10
4.803
×
1
0
-
8
15
5.382
×
1
0
-
11Figure 1
Residual squares oflog10εmT versus c0. Solid line: first-order approximation; dashed line: third-order approximation; dash-dotted line: fifth-order approximation; dash-dot-dotted line: seventh-order approximation.Also, we compare our HAM solutions of waves propagating beneath an elastic plate floating on a fluid of finite depth with those results obtained by perturbation techniques, as shown in Figure2. It should be noted that the perturbation-series solution is derived by substituting the series expansions (4.5) and (4.6) in [9] into the nonlinear PDEs (8)–(12), and equating power of small parameter ϵ leads to a succession of linear PDEs, and then the linear PDEs can be solved by the separation of variables. In Figure 2. It is seen that our homotopy-series approximation of the surface elevation ζ agrees well with the perturbation-series approximation, and only slight derivations occur at the trough of the wave profile as in Figure 2, which further indicates the validity of our present theory about nonlinear hydroelastic waves beneath a floating ice sheet.Figure 2
Comparison of our present 3rd-order surface elevationζ with those obtained by the perturbation method. Solid line: perturbation-series solution; dashed line: homotopy-series solution.We define quantities which measure how much energy there is in the wave propagating beneath an infinite elastic plate. Let P.E. be the mean potential density per unit length in theX-axis [27]. In terms of the wave surface elevation function, the energy density can be written as
(44)P.E.=14π2∫02π‍ζ2(X)dX.Different from all research objectives in [7–9], we firstly consider in this paper the effect of water depth on nonlinear hydroelastic waves beneath a floating elastic plate in detail. The energy of hydroelastic waves for different Young’s moduli of the plate E and different plate thicknesses h in various water depths are as shown in Figures 3 and 4 and Tables 2 and 3, respectively. We find that, when water depth h is about more than 2, the hydroelastic waves traveling beneath the thickest plate always contain the largest wave energy in different water depths. And with an increasing Young’s modulus of the plate, the wave energy becomes large too.Table 2
P.E. for (44) with different plate thicknesses and various water depths h.
h
P.E.
P.E.
P.E.
(
d
=
0.001
)
(
d
=
0.005
)
(
d
=
0.01
)
1
0.00067245
0.00100093
0.00040836
3
0.00031356
0.00055509
0.00069567
5
0.00030650
0.00053375
0.00074396
10
0.00030590
0.00053158
0.00074694
15
0.00030592
0.00053159
0.00074694
20
0.00030595
0.00053159
0.00074696
30
0.00030600
0.00053159
0.00074696
∞
0.00030600
0.00053159
0.00074696Table 3
P.E. for (44) with different values of Young's modulus of the plate E and various water depths h.
h
P.E.
P.E.
P.E.
(E=108)
(E=109)
(E=1010)
1
0.00076484
0.0009946
0.00040836
3
0.00038884
0.00054698
0.00069567
5
0.00030138
0.00047932
0.00074396
10
0.00028886
0.00046970
0.00074694
15
0.00028884
0.00046969
0.00074694
20
0.00028884
0.00046969
0.00074696
30
0.00028884
0.00046969
0.00074696
∞
0.00028884
0.00046969
0.00074696Figure 3
P.E. for (44) versus the water depth h for different plate thicknesses d. Solid line: d=0.001; dashed line: d=0.005; dash-dot-dotted line: d=0.01.Figure 4
P.E. for (44) versus the water depth h for different Young’s moduli of the plate E. Solid line: E=108; dashed line: E=109; dash-dot-dotted line: E=1010.The effect of Young’s modulusE of the plate on the wave elevation ζ(X) under a floating elastic plate is studied. Figures 5 and 6 show the differences in ζ(X) for E=108, 109, and 1010. According to Figures 5 and 6, respectively, we can see that the nonlinear hydroelastic response of the waves becomes flatter at the crest and steeper at the trough due to the larger value of Young’s modulus E. Finally, we consider the impact the plate thickness d by increasing d from 0.001 to 0.01. In Figures 7 and 8, we show several displacements ζ(X) with d=0.001, d=0.005, and d=0.01, respectively. It indicates that the results are very similar to the effects due to different Young’s moduli E of the plate.Figure 5
Variation of the plate deflectionζ(X) near the crest versus X for different Young’s moduli of the plate E. Solid line: E=108; dashed line: E=109; dash-dot-dotted line: E=1010.Figure 6
Variation of the plate deflectionζ(X) near the trough versus X for different Young’s moduli of the plate E. Solid line: E=108; dashed line: E=109; dash-dot-dotted line: E=1010.Figure 7
Variation of the plate deflectionζ(X) near the crest versus X for different plate thicknesses d. Solid line: d=0.001; dashed line: d=0.005; dash-dot-dotted line: d=0.01.Figure 8
Variation of the plate deflectionζ(X) near the trough versus X for different plate thicknesses d. Solid line: d=0.001; dashed line: d=0.005; dash-dot-dotted line: d=0.01.


## 5. Conclusions

In this paper, the nonlinear hydroelastic waves propagating beneath a two-dimensional infinite elastic plate floating on a fluid of finite depth are investigated analytically by the HAM. Mathematically, for a train of nonlinear hydroelastic waves traveling at a constant velocity in a fluid of finite or infinite depth, the PDEs in [7–9] were obtained by simply eliminating the time-dependent terms from the kinematic and dynamic boundary conditions on the unknown free surface in the frame of reference moving with the wave. Here, for a general case it should be noted that we construct the PDEs by directly applying the traveling-wave method to transfer the temporal differentiation into the spatial one in a fixed Cartesian coordinate oxz. Furthermore, the convergent homotopy-series solutions for the PDES are derived by the HAM with the optimal convergence control parameter.Physically, we study the effect of the water depth on the nonlinear hydroelastic waves under an elastic plate in detail. It is found that, in different water depths, the wave energy density (P.E.) tends to become larger with an increasing thickness of the sheet. The same conclusions are obtained in various water depths by means of different values of Young’s modulus of the plate. Additionally, the influences of Young’s modulus and the thickness of the plate on the wave elevationζ(X) are investigated, respectively. As Young’s modulus of the plate increases, the wave elevation becomes lower. And the increasing thickness of the plate flattens the crest and sharpens the trough of the wave profile. The results obtained here demonstrate that Young’s modulus and the thickness of the sheet have important effects on the energy and the profile of nonlinear hydroelastic waves under an ice sheet floating on a fluid of finite depth.



---

*Source: 108026-2013-10-27.xml*