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

and Applied Analysis 3 where p e 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:


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][5][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  0 .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][14][15][16][17][18], heat transfer [19,20], solitons and integrable models [21][22][23][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.

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 depth ℎ and a zero draft.A Cartesian coordinate  is used in which the -axis points vertically upward, while  = 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 (, , ) can be written as where (, ) is the wave surface elevation.The bottom boundary condition reads 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 surface  = (, ) are, respectively, modeled as where   is the water-plate interface pressure,  is the fluid density, and  is the gravitational acceleration, for a thin homogeneous elastic plate with uniform mass density   and constant thickness .
Since we are considering long waves here, the linear Kirchhoff (Euler-Bernoulli) beam theory is applied to the floating elastic plate as follows: where   =   ,  =  Here, we consider a train of nonlinear waves traveling beneath an elastic plate with constant wave number  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][8][9].Namely, we introduce an independent variable transformation where the angular frequency  and the wave number  are given.Thus, we can express the potential function (, , ) = (, ) and the traveling wave profile (, ) = ().
Then the governing equation and the bottom boundary condition for the velocity potential are transformed, respectively, by With the transformation ( 7), ( 3), and ( 6) on  = () are given by respectively, where We combine partially (10) and (11) to gain the boundary conditions on  = () as follows: Now the corresponding unknown potential function (, ) and the wave surface elevation () are governed by ( 8), ( 9), (11), and (13).

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 by a set of base functions {cos(),  ≥ 0}, where   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  direction.Therefore, we clearly know that () 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 by a set of base functions {cosh[(+ℎ)]/cosh(ℎ) sin(),  ≥ 0}, where   are unknown coefficients.Note that the potential function (, ) 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 (, ) and (), 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: where  0,1 is an unknown coefficient.We choose as the initial approximation of wave profile () 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 () is zero.

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 N 1 and N 2 as follows where and  ∈ [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 Φ(, ; ) and linear terms of (; ) 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 Φ(, ; ) by means of the solution expression (15) which is obtained under the physical considerations as 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: So we simplify the auxiliary linear operator in (21) as follows: where L 1 [0] = 0. Note that, due to the weakly nonlinear effects, the actual frequency  is often slightly different from the linear dispersion relation  0 = √ tanh(ℎ).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.

High-Order Deformation Equations.
High-order deformation equations for the unknown   (, ),   () 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 , we have where  ≥ 1.
Note that, Φ(, ; ) at the unknown surface  = (; ) may be expressed in terms of the Taylor expansion at  = 0 instead of  = (;).The detailed derivation of the expansion of Φ(, ; ) 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  = 0 as follows: where The detailed derivation of the above equations and the expression for   and   are given in Appendix A. It should be noted that ( 27) and (28) holds on the unknown boundary  = (; ), while (35) and (36) hold on  = 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  −1 and  −1 ,   and   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   and   are presented to the second order in the coming subsection.

First-Order and
Second-Order Approximations.Substituting initial approximations ( 16) and ( 17) into (36), we can get  1 () using the inverse linear operator L 2 in (36) as follows: But now the coefficient  0,1 in the initial approximation of  0 (, z) in ( 16) is still unknown.So we introduce an additional equation to relate the solutions with the wave height: in which  is an even integer,  is an odd integer, and  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 L 1 in (35), it is easy to get the solution of  1 (, ), especially by means of the symbolic computation software such as Mathematica: We find the common solution  1 (, ) has one unknown coefficient  1,1 which can be determined by avoiding the "secular" term sin() in  2 (, ).We note that all subsequent functions occur recursively.Utilizing the linear equations ( 35) and (36) to continue with the first-order approximations we have where  , is the th unknown coefficient of   (, ) and  , is the th unknown coefficient of   ().The detailed expressions of these coefficients for  2 and  2 are given in Appendix B.
In order to obtain higher-order functions   (, ) and   (), 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  0 .

Optimal Convergence-Control Parameter.
If we fix all model parameters in our approximate series solutions, there is still an unknown convergence control parameter  0 in them, which is used to guarantee the convergence of our approximation solutions.According to Liao [12], it is the convergence control parameter  0 that essentially differs the HAM from all other analytic methods.And the optimal value of  0 is determined by the minimum of the total squaredresidual    of our nonlinear problem, defined by where where    and    are two residual square errors of boundary conditions ( 27) and (28), respectively. is the number of the discrete points, and Δ = /.In this paper, we choose  = 10.Theorem 2.1 given by Liao in [12] can guarantee the rationality of (42).So we obtain the optimal convergence control parameter  0 by the minimum of the squared-residual    , generally corresponding to    / 0 = 0.

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 of  = /5 m −1 ,  = 0.01 m,   = 900 kgm −3 , ] = 0.33,  = 10 10 Nm −2 , ℎ = 5 m,  = 0.1 m, and / 0 = 1.01 and take these data hereinafter for computation unless otherwise stated.The total residual square error    at several orders of approximation versus the convergence-control parameter  0 is shown in Figure 1.It is found that    at every order has the smallest values which corresponds to the optimal  0 .For example, as  = 7, the optimal  0 = −0.18,and the smallest value of   7 = 7.910 × 10 −6 .So, let the optimal convergence-control parameter  0 = −0.18, the total residual square error    decreases quickly as the order  increases, as shown in Table 1.It is also found that   15 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.
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 Figure 2. It should be noted that the perturbationseries 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.
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 the -axis [27].In terms of the wave surface elevation function, the energy density can be written as Different from all research objectives in [7][8][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  and different plate thicknesses ℎ in various water depths are as shown in Figures 3 and 4 and Tables 2 and 3, respectively.We find that, when water depth ℎ 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.
The effect of Young's modulus  of the plate on the wave elevation () under a floating elastic plate is studied.Figures 5 and 6 show the differences in () for  = 10 8 , 10 9 , and 10 10 .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 .Finally, we consider the impact the plate thickness  by increasing  from 0.001 to 0.01.In Figures 7 and 8, we show several displacements () with  = 0.001,  = 0.005, and  = 0.01, respectively.It  indicates that the results are very similar to the effects due to different Young's moduli  of the plate.

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][8][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 .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 () 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.the explicit expressions for Δ  −1 ,  −1 ,   , and Δ  −1 in these two conditions are given by where

Table 1 :
The total residual square error    for different approximation order  with  0 = −0.18.

Table 2 :
P.E. for (44) with different plate thicknesses and various water depths ℎ.

Table 3 :
P.E. for (44) with different values of Young's modulus of the plate  and various water depths ℎ.