#### 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.

#### 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 . 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 depth and a zero draft. A Cartesian coordinate is used in which the -axis points vertically upward, while 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 , is the flexural rigidity of the plate, 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:

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–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).

#### 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 by a set of base functions , 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 , 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 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 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 and as follows where and 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 . Note that, due to the weakly nonlinear effects, the actual frequency is often slightly different from the linear dispersion relation . In Section 4, 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 in the nonlinear operator , for simplicity, we may choose another auxiliary linear operator: where .

We let be an nonzero convergence-control parameter. It is noted that both and 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: Then, from (27) and (28), two mapping functions and vary respectively continuously from their initial approximation and to the exact solutions and of the original problem. The Taylor series of and at are where

Assume that is so properly chosen that the series in (29) and (30) converge at ; then we have the so-called homotopy-series solutions as follows:

At the th-order of approximations, we have

As shown later in the following section, the unknown terms and are governed by the linear PDEs (34)–(36).

##### 3.3. 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 .

Note that, at the unknown surface may be expressed in terms of the Taylor expansion at instead of . The detailed derivation of the expansion of at the unknown surface is given in Appendices –. Upon the substitution of appropriate series and (30) into the boundary conditions (27) and (28), we have two linear boundary conditions on 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 . Furthermore, the original nonlinear DPEs (1)–(5) are transferred into an infinite number of linear decoupled high-order deformation equations (34)–(36). Namely, given and , 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.

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

Substituting initial approximations (16) and (17) into (36), we can get using the inverse linear operator in (36) as follows:

But now the coefficient in the initial approximation of 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 .

Further, in the analogous manner as for the first-order approximation, by using the inverse linear operator in (35), it is easy to get the solution of , especially by means of the symbolic computation software such as Mathematica:

We find the common solution has one unknown coefficient which can be determined by avoiding the “secular” term in . 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 and 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 .

##### 3.5. Optimal Convergence-Control Parameter

If we fix all model parameters in our approximate series solutions, there is still an unknown convergence control parameter in them, which is used to guarantee the convergence of our approximation solutions. According to Liao [12], it is the convergence control parameter that essentially differs the HAM from all other analytic methods. And the optimal value of is determined by the minimum of the total squared-residual 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 .

Theorem 2.1 given by Liao in [12] can guarantee the rationality of (42). So we obtain the optimal convergence control parameter by the minimum of the squared-residual , generally corresponding to .

#### 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 of , m, , , , m, , and 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 is shown in Figure 1. It is found that at every order has the smallest values which corresponds to the optimal . For example, as , the optimal , and the smallest value of . So, let the optimal convergence-control parameter , the total residual square error decreases quickly as the order increases, as shown in Table 1. It is also found that is down to 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 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.

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–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 , , and . 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 to . In Figures 7 and 8, we show several displacements with , , and , respectively. It indicates that the results are very similar to the effects due to different Young’s moduli of the plate.

#### 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 . 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.

#### Appendices

#### A. The Detailed Derivation of (35) and (36) and the Expressions for and

Let For any , we have a Maclaurin series as follows: For , it follows from (A.1) and (A.2) that where Thus we have, for , where

Substituting the series expansions (A.1) and (A.5) into the boundary conditions (27) and (28) and then equating the like-power of the embedding parameter , we have two linear boundary conditions (35) and (36), respectively. And the explicit expressions for , , , and in these two conditions are given by where

#### B. Expressions of the Coefficients

#### Conflict of Interests

There is no conflict of interests in the paper. The authors themselves used the program of the symbolic computation software named Mathematica independently to gain the approximate analytical solutions of the PDEs considered here.

#### Acknowledgments

This work was supported in part by China Postdoctoral Science Foundation funded Project 20100481088 and in part by the Natural Science Foundation of Shandong Province of China under Grant ZR2010FL016. The authors would like to thank the reviewer for his constructive comments.