Hydroelastic interaction of nonlinear waves with floating sheets

Hydroelastic responses of floating elastic surfaces to incident nonlinear waves of solitary and cnoidal type are studied. There are N number of the deformable surfaces, and these are represented by thin elastic plates of variable properties and different sizes and rigidity. The coupled motion of the elastic surfaces and the fluid are solved simultaneously within the framework of linear beam theory for the structures and the nonlinear Level I Green–Naghdi theory for the fluid. The water surface elevation, deformations of the elastic surfaces, velocity and pressure fields, wave reflection and transmission coefficients are calculated and presented. Results of the model are compared with existing laboratory measurements and other numerical solutions. In the absence of any restriction on the nonlinearity of the wave field, number of surfaces, their sizes and rigidities, a wide range of wave–structure conditions are considered. It is found that wave reflection from an elastic surface changes significantly with the rigidity, and the highest reflection is observed when the plate is rigid (not elastic). It is also found that due to the wave–structure interaction, local wave fields with different length and celerity are formed under the plates. In the case of multiple floating surfaces, it is observed that the spacing between plates has more significant effect on the wave field than their lengths. Also, the presence of relatively smaller floating plates upwave modifies remarkably the deformation and response of the downwave floating surface.


Introduction
Floating elastic surfaces have been the subject of extensive research because of their widespread applications. For example, floating ice sheets of different sizes are commonly observed in marginal ice zones [55]. Wave interaction with floating ice sheets is of interest due to the effect that the ice has on the wave field and the impact of waves on the ice floes, particularly when it results in breaking of the ice to smaller pieces [38]. An elastic thin plate on a fluid foundation flexing under the action of waves is commonly used as a model for a floating ice sheet. The alternative use of a thick-plate model has been shown to be unnecessary except for the case of very short waves or very thick ice, shown by Fox and Squire [14]. Another application of floating deformable surfaces is related to Very Large Floating Structures (VLFS), which have been extensively studied as a model for big ships, floating ferry piers, bridges and airports, see, e.g., Watanabe et al. [62]. As these structures have a large surface area and a relatively small thickness, they behave elastically under wave action and can be considered in the framework of shallow-water wave theory. Therefore, the study of wave action on elastic floating bodies has received considerable attention. The next step in investigation of VLFS is to consider several separate structures situated near each other to study how they interact with each other. In this sense, a set of floating elastic sheets can model, for example, the floating fuel storage base or an offshore field of solar panels.
An early solution to the problem involving thin elastic plates on water can be found in Evans and Davies [12], where the two-dimensional problem of wave impact on a floating semi-infinite thin elastic plate was solved using the Weiner-Hopf technique. Newman [47] developed a general theory for the interaction of water waves with elastic structures by defining the deformation of the body with an expansion of arbitrary modal shape functions, so that the response in each mode was obtained as rigid body motion.
The interaction of waves with flexible structures of finite length can hardly be analyzed analytically. One of the attempts was made by Ertekin and Kim who constructed explicit solution of the Level I Green-Naghdi (GN) equations by the use of the wide-plate approximation [8,30]. Korobkin and Khabakhpasheva [36] found the way to construct exact solutions of the problem of the flexural-gravitational oscillations of a floating elastic plate by an inverse method, but only due to the small-amplitude periodic external loads. Therefore, the general problem of two-dimensional water wave scattering by a single elastic plate of finite length has typically been treated numerically. There are series of numerical methods that can be applied to hydrodynamic analysis of floating elastic surfaces, such as the finite-element method (FEM) [59], boundary element method (BEM) [28], the modal expansion method [35,37,64] and the direct method [1]. Utsunomiya et al. [60] and Liu and Sakai [39] combined BEM for evaluating the fluid motion with FEM for calculating the response of the elastic structure. According to the modal expansion method, the deformation of the plate is decomposed into vibration modes that are chosen arbitrarily. In this regard, researchers have adopted different modal functions: eigenfunctions [64], bi-cubic B-spline functions [29], to name a few. For example, Ren et al. [51], in studying the hydroelastic waves in a channel covered by cracked ice with arbitrary edge conditions at the walls, expanded the solution into a series of cosine functions with unknown coefficients, which can be found analytically. Usually, the modal expansion method is limited to periodic waves in frequency domain. Meylan [42] applied the spectral theory to determine the time-dependent motion of a thin plate in waves. In the direct method, the deformation of the sheet is determined by solving the equations directly without expanding the plate motion into eigenmodes. Andrianov and Hermans [1] examined wave scattering by a finite elastic plate for infinite, finite and shallow water depths by deriving an integro-differential equation.
More complicated problem statements as the motion of a floating plate in water of variable depth have been studied by Wang and Meylan [61] with BEM for homogenous plate, by Porter and Porter [50] by variational approach for heterogeneous plate and by time-domain modal expansion method by Sturova [56]. Three-dimensional model of wave interaction with circular ice floes has been presented by Meylan and Squire [44] and later extended to ice floes of arbitrary shape and properties. Bennetts and Williams [3] proposed an efficient solution for a floating elastic plate of arbitrary shape by reducing the problem of water-wave scattering to the system of one-dimensional integro-differential equations. Meylan [43] used the eigenfunction-matching method to solve for the forced vibration of a plate of constant thickness and extended it to three dimensions with little modification of the underlying numerical method.
In this study, we consider a two-dimensional fluid of finite depth covered by a finite number of elastic surfaces of different properties. One way to deal with this problem is to consider the series of plates as a collection of cracks in an infinite or semi-infinite plate. Evans and Porter [13] provided an explicit solution for the scattering of an incident flexural-gravity wave by a narrow straight-line crack in an infinite elastic plate and later extended it to the case of multiple cracks [49]. Their solution expresses the potential in terms of a linear combination of the incident wave and certain source functions located at each of the cracks. This solution, however, is limited to plates of identical properties. Similarly, Guyenne and Pǎrǎu [21] approached the problem of wave attenuation by irregular array of ice floes by using mathematical formulation with elasticity function varying between zero in open waters and nonzero in the ice-covered region. Their solution was based on the full time-dependent equations for nonlinear potential flow, but was limited to a solitary wave case and the ice of negligible mass and thickness. Another way to solve the problem was proposed by Ogasawara and Sakai [48], where BEM and FEM of Liu and Sakai [39] were incorporated to numerically solve the problem for a set of arbitrary plates of finite length by the use of a time-domain solution. Kohout et al. [34] solved a twodimensional multiple floating elastic plate problem via an extension of the matched eigenfunction expansion method of Fox and Squire [15]. Within this approach, they considered regions of open water as arising from limiting cases of plates of vanishing thickness. In three-dimensional context, the multiple plate problem has been studied by Bennetts and Squire [2], who constructed the model of wave scattering by multiple rows of identical and equally spaced ice floes of circular shape, requiring periodical distribution of the floes on the water surface. Unlike most of the works described thus far, where the draft of the floating plates has been neglected, in the present analysis the draft is relevant.
Most of the methods discussed above operate within the limits of linear potential theory and/or are timeconsuming as they require large memory capacity on a computer. In a real sea state, when the wave steepness is large or when nonlinearity is not negligible, the linear solution to this problem will no longer be satisfactory. There are only a few numerical models that consider global hydroelastic response of a floating elastic surface to nonlinear waves. Recently, Dinvay et al. [4,5] with the use of Hamiltonian formulation derived the fully dispersive and weakly nonlinear equation of Whitham type. They demonstrated that small-amplitude solution to this equation corresponds well to the fully nonlinear solution of hydro-elastic Euler equations. This approach allows for a description of waves of any wavelength propagating in a thin elastic sheet lying on the surface of an inviscid fluid. Xia et al. [65] and Ertekin and Xia [9] used the nonlinear GN equations, following the theory of directed fluid sheets, to study the elastic response of a thin plate to nonlinear waves, and showed the importance of nonlinearity in the structural response of a floating elastic plate. They coupled the Level I GN equations with elastic beam theory through the kinematic and dynamic boundary conditions to obtain a new set of modified GN equations. Later Hayatdavoodi and Ertekin [23] applied successfully the same model to nonlinear problem of wave scattering by a submerged rigid plate. A comparative study of the results of a select group of hydroelastic approaches, including a GN-based model, is provided by Riggs et al. [52].
To investigate the interaction of nonlinear water waves with a set of elastic surfaces of arbitrary properties, we use the Level I GN theory coupled with the thin plate theory. For this purpose, the numerical scheme developed by Xia et al. [65] for a single plate is revisited and it is modified by considering N number of arbitrary deformable surfaces and by considering appropriate matching and jump conditions at the edges of the plates. The jump conditions derived by Xia et al. [65] through the use of the postulated conservation laws of mass, momentum and mechanical energy are substituted in the present study with the continuity of pressure along the bottom surface, shown originally by Hayatdavoodi and Ertekin [24].
The rest of the paper is organized as follows. The mathematical problem is formulated in the next section, which is then solved in Sect. 3, where the numerical scheme is explained. The results of the model are compared with available numerical and experimental data in Sect. 4. The velocity field and pressure distribution in the fluid domain are presented in Sect. 5. In Sect. 6, wave scattering by the deformable plates is studied by defining the reflection and transmission coefficients. The results for wave interaction with multiple surfaces are presented in Sect. 7, and this is followed by the conclusions of this study in Sect. 8.

Problem formulation
The plane flow of inviscid fluid is considered in the Cartesian coordinate system Ox y with the horizontal axis lying on the undisturbed free surface and vertical axis directed upwards. Incident waves propagate in the positive x-direction and excite the motion of a set of floating deformable sheets being initially at rest. The sketch of the problem is shown in Fig. 1. Each sheet has the width L i , draft d i , uniform mass per unit width m i and flexural rigidity D i . The subscript i denotes the i-th sheet in the sheet formation and in the case of a single sheet is not used. The depth of the fluid under the i-th sheet at rest is defined as h i = h − d i . The elastic sheets are restricted to vertical motion only and are always in contact with the fluid. The loss of energy due to the structural damping is neglected and over-washing is disregarded. The role of overwash in the wave-plate interaction was studied in the analysis of Nelli et al. [46]. The physical variables are given in dimensionless form by use of ρ, g and h, selected as a dimensionally independent set, where ρ is the mass density of the fluid, g is the acceleration due to gravity and h is the constant depth of the fluid, i.e., Due to the equilibrium condition, mass per unit width is expressed as m i = ρd i in relation (2), which means that dimensionless magnitudes of draftd i and mass per unit lengthm i are equal. For consistency, all equations, input parameters and resulting figures will be given in dimensionless form with the tildes dropped.
In this study, we will use the nonlinear Level I GN equations, originally developed by Green and Naghdi [16,17], to study the wave-structure interaction problem. The GN equations can be obtained through various approaches, conserve mass identically and satisfy the integrated momentum equations and nonlinear boundary conditions exactly. The GN equations are classified based on the functions used to prescribe the distribution of the vertical velocity along the fluid sheet. The Level I GN equations used here assume a linear distribution of the vertical velocity along the water column and are the original set of the GN equations given by Green et al. [19], and a specialized version of the more general equations given later by Green and Naghdi [18]. Prescription of the variation in the vertical velocity is the only assumption made about the kinematics of the fluid flow in the GN equations. Given the linear variation in the vertical velocity over the water column (and consequently, invariant horizontal velocity over the water column), the Level I GN equations are mostly applicable to the propagation of long waves in shallow water.
In the context of the Level I GN equations, the problem is best studied by dividing the fluid domain into two types of regions. Region I (R I ) is formed by a flat floor at the bottom and by a free surface at the top, where the fluid pressure is constant atmospheric pressure. Region II (R I I ) is formed by a flat floor at the bottom and by an elastic plate at the top, where, as opposed to Region I, the fluid pressure is variable. Solutions, obtained in each region, are connected through the proper matching and jump conditions at the interfaces.
The governing equations for the motion of the fluid in R I are provided by the Level I GN theory for a flat and stationary seafloor. They can be written in a compact form as [10]: where u(x, t) is the horizontal component of fluid velocity and η(x, t) is the free surface elevation, measured from the still-water level. In our notation, subscript after comma denotes the partial differentiation with respect to the given variable and superimposed dot specifies the total time derivative. Similarly, the governing equations for coupled motion of the fluid and elastic plate in R I I can be written as: where ζ(x, t) denotes the plate deformation, measured from its still position. The wave-induced pressure under the plate is given by the thin elastic plate theory [57]: where the flexural rigidity is defined by with δ i , E i , and ν i being the thickness, Young's modulus and Poisson's ratio of the i-th plate, respectively. System of Eqs. (3)-(7) can be integrated for the unknown functions η( Ertekin [6] provided explicit relations for the vertical velocity of the fluid in Level I GN equations, given as: and pressure on the bottom (y = −1), given as: The nonlinear kinematic and dynamic conditions on top and bottom boundaries of the flow domain are satisfied exactly by Eqs. (3)-(10).
In deriving the GN equations, information about pressure over the water column is limited to that of the integrated pressure, pressures at the bottom and top surfaces. Pressure distribution over the water column can be found from Euler's equations for the vertical component of fluid velocity, written in dimensionless form as: Substituting the relation for vertical velocity (9) into Eq. (11) and integrating with respect to y yield: In integration, the condition of zero pressure at the free surface y = η(x, t) and pressure distribution (7) at elastic surfaces y = ζ(x, t) should be used. Pressure (12) at the bottom surface y = −1 appears to be equal to the bottom pressure (10), predicted by the GN theory. This is of no surprise, as Green and Naghdi [16] showed that the differential equations of the direct theory can be obtained from the exact equations of an incompressible and inviscid fluid, using only a single approximation for the velocity field, equivalent to the Level I assumption in the direct approach. We note that High-level GN equations, applicable to wave propagation in any water depth, can be obtained by assuming higher-order polynomials (or alternatively exponential functions) for the vertical velocity distribution over the water column. As such, more information about the vertical distribution of the velocity and pressure fields over the water column can be described. Further discussion about the High-level GN equations can be found in Shields and Webster [54], Webster and Kim [63], and more recently in Zhao et al. [66][67][68], who discussed the High-level GN equations.
To formulate the problem fully, suitable boundary conditions at the fixed positions of the leading (x = x L i ) and trailing edges (x = x T i ) of each floating sheet must be specified. Since the elastic sheets are floating freely, the bending moment and the shear stress should vanish at the edges, i.e., Moreover, as there is no gap, by assumption, between the bottom surfaces of the sheets and the top surface of the fluid layer, the mass continuity Eq. (5) together with vanishing bending moment condition (13) at the edges implies that: In the approach discussed above, the floating elastic surfaces cause discontinuities of the unknown functions and their derivatives at the interfaces between regions. Appropriate jump conditions should be called for to provide the matching of the solution. The theory demands the conservation of mass across the discontinuity curves, which is given by where x L i ± 0 and x T i ± 0 denote the single-sided limiting values of x L i and x T i , respectively. The physics of the problem and the theory require the continuity of the bottom pressure, obtained from Eq. (10), across the interfaces between regions: On the left side of the domain, a numerical wavemaker is installed [7]. The wavemaker is capable of generating GN solitary waves of amplitude A and cnoidal waves of height H and length λ. On the right side of the domain, Orlanski's condition is used to minimize the wave reflection back to the tank: where c is the phase velocity of the incoming wave, see Ertekin [6] for more details.
The set of Eqs. (3)- (19) formulate fully the problem of interaction of surface waves with the multiple floating elastic sheets. This nonlinear system is solved numerically, and this is discussed in the next section.

Numerical solution
The governing equations in both types of regions, along with the boundary conditions at the discontinuity curves at the leading and trailing edges, are solved simultaneously for the unknowns. Following the same approach, discussed in Ertekin [6], we eliminate the time derivatives of ζ from equations (5)-(7) by applying the mass continuity Eq. (5) in the momentum Eq. (6). As a result, we obtain where functions Y I I f and Y p , accounting for the fluid motion in R I I and plate deformation, respectively, include only spatial derivatives of the unknown functions and can be written as Equation (20) where An explicit time stepping method is required to solve equations (20)-(24) numerically. In this study, the modified Euler method is applied to perform the time integration. It is completed in two steps by using the first-order Euler method twice. The spatial discretization of the derivatives is performed using the central difference formulas for the nodes inside the domain. Both time and space stepping methods are second-order accurate. The unknown functions η(x, t), ζ(x, t) and u(x, t) are approximated by discrete sets η n j , ζ n j and u n j , respectively, where j is the mesh point in spatial domain and n specifies a mesh point on the time axis. In this paper, the space grid x and time step t are kept constant. At each given time, such discretization results in an array of unknowns and in a banded matrix with nonzero elements below and above the main diagonal. The Gaussian elimination algorithm for a banded diagonal matrix is used to solve the system of equations. This method has been applied to GN equations by Xia et al. [65], Hayatdavoodi and Ertekin [23] and Ertekin [6] successfully. To enforce the free-end boundary conditions (13)- (15) and jump conditions (16)- (18), fictitious points before the leading and after the trailing edges of each plate are introduced so that the real nodes are reserved for equations (20) and (23). Moreover, the fictitious points compensate the missing nodes for calculation of the higher derivatives at the edges.
The computational procedure is explained in detail as follows. Time-stepping calculation starts from the still-water surface by prescribing water elevation and horizontal velocity on the wave maker to the left and at the wave absorber to the right. In the solitary wave case, the surface elevation and velocity are prescribed at the  (5) according to the region type. Then, surface elevation at fictitious nodes between the regions is calculated using the free-free beam end boundary condition (13). Next, coefficients in the momentum Eqs. (20) and (23) are evaluated as well as in the boundary and jump conditions (14)- (18). Finally, the set of linear equations on the internal nodes, resulting from momentum Eqs. (20) and (23), together with those at the fictitious points are solved by the Gaussian elimination method to obtain the velocity at subsequent time step.
Results of calculations may experience saw-teeth noise, mainly due to the use of the central difference scheme, especially in Eq. (20) containing the derivatives up to the fifth order. This is not physical and may cause instability in numerical solution. To make the numerical scheme stable without affecting the results, we adopt the five-point filtering formula, developed by Shapiro [53] and utilized by Ertekin et al. [10] for a similar problem. The scheme reads as where f can represent η, ζ or u. It should be noted that although saving the calculations from unexpected oscillations, frequent filtering may lead to numerical damping. The filtering time step n F , most appropriate for the given discretization, is determined small enough to prevent the solution from diverging, but big enough to avoid unwanted dissipation. Figure 2 shows the convergence of the model comparing the time series of deformation at the center of the plate obtained from the same simulation run using different grids, given in Table 1. To ensure numerical stability, the ratio x/ t must be no less than the wave propagation speed both on the free surface and under the structure. Tsubogo [58] has shown that the wave propagation speed under the elastic plate depends strongly on its flexural rigidity. With dispersion relation analysis, he proved that the hydroelastic waves propagate faster than the incident water waves and the group velocity of the hydroelastic waves increases according to the increase in the bending rigidity. Hence, the bigger the rigidity parameter D is, the smaller should the ratio t/ x be for convergence of the numerical scheme. For rigidities considered in this work, the inequality t/ x ≤ 0.1 should be satisfied. Figure 2 demonstrates that smaller ratio of t/ x provides better convergence of the solution, as expected. The grid specification, used in this work, corresponds to Grid 6 in Table 1. It is optimal to adequately model the problem and does not require excessive computational time. To be exact, for 2.7 GHz Intel Core i5 Central Processor Unit it takes 2.15 seconds to run 100 time units in a numerical tank of total length of 100 water depths.

Comparisons with experiments and numerical data
To confirm the validity of the model described in the previous sections, comparisons are made with the available numerical and experimental data. The experiments cited here were conducted in two-dimensional wave flumes filled with water of mass density 1000 kg/m 3 . The gravitational acceleration g is 9.8 m/s 2 . Few data are available on nonlinear periodic waves interacting with the floating elastic sheets in shallow waters. Most of the existing studies are focused on linear and regular waves. In linear theory, it is conventional to describe the dynamic response of the plate to regular waves by plotting the amplitude of deformation and bending moment at each point of the structure, scaled by the incident wave amplitude H/2 and by the product d L H/2, respectively. In this work, however, we use nonlinear cnoidal waves. Hence, the wave conditions of the existing data and the wave conditions of our study differ slightly, mainly depending on how large the incident wave height is. To compare the present theory with existing numerical results and experimental data, the absolute values of  Figures 3 and 4 show that the plate experiences less deformation and hence accumulates more stress in response to shorter waves. However, Wu et al. [64], having considered a wide range of the incoming wavelengths, took notice that the peak of the maximum bending moment does not occur at the lowest possible wavelength. They stated also that the maximum bending moment becomes smaller as the floating body becomes more flexible. Figures 3 and 4 demonstrate good correspondence between the results. In obtaining the results presented in Figs. 3a and 4a, a finer mesh size ( x = 0.1, t = 0.01) than that discussed in the previous section is used. This is because the wave is significantly shorter in this case than the rest of the cases considered in this study. Figures 5 and 6 bring together results of the present theory, the numerical and experimental data of Kohout et al. [34] for two and four elastic plates, respectively. The experiments were conducted in a wave flume of 8 m long at a water depth of 0.6 m with two 4 m sheets and four 2 m sheets of 20 mm thickness, mass density 914  [34]. The only divergence of both numerical solutions from the experimental data observed in the area around the fluid gap between two plates (Fig. 5) may be attributed to a potential defect in experimental measurements in this particular case. In fact, experimental and calculational data provided by Kohout et al. [34] for other wave cases with collections of two and four plates showed good correlation without any contradictions in the regions dividing the adjacent plates (e.g. Fig. 6).  Figure 7 shows excellent agreement between the results. Figure 7 shows that the floating structure introduces sufficient transformations into the initial form of the solitary wave. The reduction in the wave amplitude results in increase in its length. As soon as the wave enters the plate region, the small-amplitude leading waves appear before the main wave crest. These waves grow larger as the solitary wave propagates along the plate. The above comparisons show very good correspondence of the present theory with previous numerical methods and experimental data for both cnoidal and solitary wave cases, as well as single and multiple plates. From these cases, we conclude that the numerical solution developed here can be used to predict the hydroelastic response of the complex systems of floating elastic plates under various nonlinear wave conditions. Extensive comparisons of the GN results with laboratory experiments of solitary and cnoidal wave loads on horizontal, submerged rigid plates can also be found in Hayatdavoodi et al. [22,23] and more recently by Hayatdavoodi et al. [27], where similar conclusions about the agreement of the GN model with laboratory measurements were made.

Velocity field and pressure distribution
We recall that the only assumption made about the fluid motion in the Level I GN equations is that the vertical component v of the particle velocity varies linearly along the fluid depth. This assumption, along with the incompressibility condition (3) or (5), results in the horizontal component u of the particle velocity to be uniform over the fluid thickness. Since no assumption is made about the irrotational character of the flow, such a velocity field allows for rotational motion. Indeed, in Level I GN equations, the vorticity at any point is equal to v x , (given that u y = 0), which is a non-zero function. We note that GN equations can also be obtained for irrotational flows (known as Irrotational Green-Naghdi equations, IGN), see, e.g., Kim (23), while the vertical component v is determined analytically from Eq. (9). Greater wave attenuation is observed in longer plates with higher rigidity. Hence, to be able to distinguish the nonlinear wave motion in full details, a long plate of sufficiently large rigidity has been chosen. Figure 8 shows that due to the small-amplitude leading waves propagating with higher velocity along the elastic surfaces, the fluid under the trailing edge of the plate feels the wave motion long before the wave reaches it. Moreover, the solitary wave drives the perturbations along the floating surface mainly in one direction, so that the fluid and the structure behind the wave are almost at rest. Due to these elastic waves, the plate attenuates the main wave and consequently reduces its propagation speed. In the moment, the solitary wave reaches the trailing edge of the plate, its propagation speed is reduced by half, and the wave energy is transmitted to the fluid downwave.
It is interesting to observe in detail the process of a cnoidal wave train transformation at the edges of the plate. Figure 9 shows snapshots of the fluid flow close to the leading edge of the plate at nine successive time moments during one incoming wave period (the same observation can be carried out at the trailing edge of the plate). As expected, the flow is discontinuous at the interface between Regions R I and R I I . It appears as if the fluid is stored for some time before it collects the required momentum to break through the boundary. This happens when the wave crest reaches the plate and the leading edge starts to uplift. In this moment, the wave penetrates into the region under the plate. When the wave enters completely into Region R I I , the free end of the plate drops down and the fluid is pushed back under the plate, resulting in a reversed motion. The process is accompanied by formation of zones close to the edge with extremely high fluid velocity. In fact, the largest plate deformation is observed at the edges (see Fig. 3).
In Fig. 10, the difference between rigid and flexible plates and different effects they have on the wave field is studied by considering the velocity field under the plates of low, high and infinite rigidity. In case of a rigid plate, ζ ,x = ζ ,x x = 0 and the equations of motion (5)- (6) in Region R I I are reduced to: Hayatdavoodi and Ertekin [23] showed that fluid motion under the rigid plate is calculated by integrating Eq. (27) Figure 11 presents the pressure distribution p under the plate calculated by use of Eq. (12). Isolines of fluid pressure reproduce the free surface shape in Region R I , but diverge from the shape of the plate in Region R I I , because of nonzero pressure along the fluid-structure interface. Although being discontinuous at the plate edges, the fluid pressure remains continuous along the bottom surface, which is provided by the jump conditions (17)- (18). As shown in Fig. 12, the pressure p calculated by Eq. (12) is compared with the linear  Figure 12 reveals the linear distribution of pressure along the water column for both solitary and cnoidal waves. This is in line with the conclusions of Neil et al. [45] and Hayatdavoodi et al. [26] (who used the GN, Boussinesq, and linear equations) about the stratification of pressure during nonlinear wave interaction with the structure in shallow water.

Wave reflection and transmission
In this section, the problem of wave scattering by elastic plates on the surface is analyzed. When an incident wave interacts with floating elastic plates, one part of the wave energy is transmitted past the body, and another part is reflected back, while a small portion of the energy is locked between the leading and trailing edges of the plates [25]. To quantify this phenomenon, reflection and transmission coefficients are considered. These coefficients use the change in wave height upwave and downwave to represent the energy reflection and attenuation. In this study, the four-gauge (two gauges upwave and two gauges downwave) method introduced by Grue [20] is used. The same method has been applied successfully by Hayatdavoodi et al. [25] to determine wave reflection and transmission by a submerged horizontal plate who used the Level I GN equations. Upwave of the plates, the surface elevation is a superposition of the surface elevations of the incoming wave η I and the reflected wave η R . In this approach, the nonlinear wave is reconstructed by a superposition of linear waves. Let a I and a R be the amplitudes of the fundamental harmonics of the incident wave and the reflected wave, respectively. Then, these waves can be approximated in their linear form as where k = 2π/λ is the wave number, ω is the angular wave frequency, ε I and ε R are the phases of the incident and reflected waves, respectively. Generally, η I and η R are nonlinear and hence can be decomposed, apart from the fundamental frequency component, into free and bound (phase-locked) harmonics. Nevertheless, as it was shown in Hayatdavoodi et al. [25], the fundamental harmonics may be sufficient for constructing the reflection and transmission coefficients. Suppose that surface elevation is recorded at two adjacent gauges with coordinates x 1 and x 2 = x 1 + δx. The observed profiles can also be presented by the fundamental harmonics as The unknown amplitudes A 1 , B 1 , A 2 and B 2 are evaluated for the given time series of η(x 1 , t) and η(x 2 , t) by use of the Fourier transform for the fundamental frequency ω. This frequency is determined either as an input parameter or it is estimated for any wave number k from the dispersion relation of the linearized Level I GN equations as [19] Eqs. (29)-(32) can be solved for a I and a R as follows: The amplitude of the transmitted wave a T is calculated the same way as the amplitude of the incident wave a I by the use of the two wave gauges located downwave of the plate with the coordinates x 3 and x 4 = x 3 + δx . Once the amplitudes of the incident, reflected and transmitted waves are determined, the reflection and transmission coefficients are defined by In this study, reflection and transmission coefficients are considered only for the cnoidal wave case. See Hayatdavoodi et al. [25] for a discussion on solitary wave scattering. In the discussion below, the plots of the reflection and transmission coefficients are smoothed to eliminate noise using the MATLAB package function smoothdata.
Comparison of the reflection and transmission coefficients, due to cnoidal wave interaction with an elastic plate on the surface, calculated here by the GN equations and the results of Maiti and Mandal [40], is shown in Fig. 13. Maiti and Mandal [40] used the eigenexpansion method for the potential flow to solve the wave scattering problem for the elastic plate floating over a porous floor. As shown in Fig. 13, the reflection coefficient increases rapidly with a decrease in wavelength, i.e., as the ratio λ/L becomes smaller. However, maximum wave reflection is less than 10% of the incident wave, and most of the wave energy is transmitted through the plate. The transmission coefficient, as expected, increases nonlinearly with increasing wavelength. Overall, good agreement is observed between the results of the GN equations and by linear equations of Maiti and Mandal [40] for the rigid bottom as a special case of the porous bed, shown in Fig. 13. The minor differences are likely due to nonlinear effects on the transmitted wave, which are not considered in the definition of the reflection and transmission coefficients, discussed in this section. We also note that nonlinear cnoidal wave is used in the GN model which differs from the linear wave considered by Maiti and Mandal [40].
Variation in the reflection coefficient with the plate length is shown in Fig. 14 and compared with the results of Meylan and Squire [41], who solved the problem in linearized form by the use of the Green function. Figure 14 demonstrates an infinite set of resonances, where, as it was noticed by Meylan and Squire [41], nearly perfect transmission occurs. Minor discrepancy between the results is observed for the smaller plate lengths, which is likely due to the effect of the draft of the plate not considered in Meylan and Squire [41] and at the peaks of the plots, due to the reasons discussed previously. Overall, Fig. 14 shows good agreement between both approaches.
In Fig. 15, the effect of plate rigidity D and mass m on wave transformation is studied. In comparison with elastic plates of different rigidities, the rigid plate is characterized by the highest wave reflection. On the contrary, the plate of least rigidity remains almost invisible to the waves of any length. Reflection coefficient also increases with mass of the plate, as it is directly connected with its draft.
The relation c 2 R +c 2 T = 1, applicable to linear waves, is not always exactly satisfied for the waves considered in this study. Figure 16a shows that with increase in wavelength the relation c 2 R + c 2 T goes to unity, but it sags down for greater values of H . This is largely due to the wave nonlinearity, the use of linear dispersion relation (33) and the effect of the higher harmonics, not included in decompositions (29)- (32). By considering the next harmonic of the wave signal measured at the first gauge upwave the plate, Eq. (31) changes to η(x 1 , t) = a 1 1 cos(ωt + ϕ 1 ) + a 2 1 cos(2ωt + ϕ 2 ), where a 1 1 and a 2 1 are the amplitudes of the fundamental and the second-order harmonics, respectively. Figure 16b shows that contribution from the second harmonic for short incoming waves is small, when compared to the fundamental harmonic, and thus can be neglected. However, contribution of the second harmonic increases with longer wavelength. Therefore, to minimize the error, in the analysis of the wave transformation by multiple plates in the next section, we confine our attention to those incoming waves for which the scattering process can be described in terms of the fundamental harmonics only, i.e., waves with heights as small as H = 0.1. This approach in determining the reflection and transmission coefficients can be expanded to become more suitable to the incoming waves of large height and long length by considering the contribution of the higher harmonics in Fourier expansion of the signals given in Eqs. (31)-(32).

Wave interaction with multiple flexible surfaces
In this section, we investigate the flexural motion and scattering characteristics of floating formations consisting of more than one elastic surface under cnoidal wave forcing. In such formations, each plate affects the motion of others and vice versa. Generally, plates may have different lengths, mass and flexural rigidities, and they may be located at arbitrary distances from each other. Section 6 discusses that higher mass and rigidity of the plate cause greater wave attenuation. Thus, in the current analysis, all elastic sheets are made of the same material, i.e., have the same dimensionless mass m = 0.01 and flexural rigidity D = 1. The formation of elastic sheets is then characterized by the number of plates, their lengths and distances between them. The four-gauge method explained in the previous section is used to determine the reflection and transmission coefficients. From the long-wave perspective (λ/L > 1.4), scattering characteristics are more or less the same regardless of the case: c R is sufficiently small (c R < 0.1) and c T is very close to 1. Indeed, long waves may reach tens to hundreds of kilometers into the ice-covered ocean [55]. Very short waves (λ/L < 0.4) are out of scope of the current model that covers the shallow water case. Therefore, in subsequent analysis, we focus on the waves having the wavelength λ ∈ (2L/5, 7L/5). Figure 18 compares the variations in the reflection c R and transmission c T coefficients with the incident wavelength to plate length ratio λ/L (relative wavelength) in Case A [N ] for different values of integer parameter N . Figure 18 shows that coefficients c R and c T vary nonlinearly with λ/L having distinguished peaks and troughs at different positions for different number of plates N . An interesting observation can be made from Fig. 18 that a single elastic plate is characterized by lower wave reflection and higher transmission than the same plate broken into equal parts. In other words, cracks in the ice make it more resistant to the incident waves and reduce the wave transmission further into the ice field. That is accounted for by multiple reflection from the edges of the plates. Figure 18 shows that smaller waves are better reflected by the group of smaller sheets, which conforms with Squire et al. [55], who stated that ocean waves better interact with ice floes of comparable lengths. Nevertheless, it should be noted that for higher N , when the plate is decomposed into numerous arrays of tiny plates, the reflection coefficient drops to zero quickly, since the size of each plate in such a configuration is small even in comparison with the relatively short incoming waves.
In Fig. 19, coefficients c R and c T in Case B[l] are plotted as functions of λ/L for several values of the distance parameter l. Figure 19 shows that distance between two equally sized sheets has little effect on their scattering performance. Significant increase in relative position of two plates from L/10 to L/2 leads to insignificant changes in the location and magnitude of peaks in the plots of the wave reflection coefficient. At the same time, the transmission coefficient becomes smaller for larger l, since more wave energy can be trapped within wider space between the plates.
It is expected that with increasing distance between two floating sheets, the interaction between them becomes weaker. Figure 20 compares the amplitudes of plate deformation and bending moments at each point of the structures (see Sect. 4 for details) in one plate and in two plates divided by the fluid gap. Figure 20 shows that, as the distance between the plates increases, the amplitudes of deformation and bending moments in the second plate decrease gradually. The first plate oscillates with almost the same amplitude as a single plate, but experiences greater stress. In contrast with the second plate, the amplitudes of deformation and bending moments in the first plate are less affected by the distance parameter l.
In Figs. 21 and 22 , velocity field and pressure distribution between two plates in Case B[L/5] at two successive time moments are captured, respectively. During one time unit, the magnitude and direction of the flow in the region between the plates change significantly. As shown previously in Sect. 5, the size of the fluid pockets corresponds to the incoming wavelength (see Fig. 10). But when the distance between the plates is less than the wavelength, regular character of the fluid pocket formation and motion is disturbed. This has an effect on the fluid motion downwave and causes reduction in stress and deformation in the second plate, as shown in Fig. 20. In Case C[L 1 , l] we investigate the effect that the position of a crack in a plate has on resulting scattering characteristics. Figure 23 provides the plots of reflection and transmission coefficients for three configurations of two plates located at the same distance l from each other. In a wide range of incoming wavelengths, the non-symmetric configurations (L 1 < L/2) have bigger wave reflection than symmetric couple of equally sized plates (L 1 = L/2), while wave transmission remains the same in all these cases. Comparing the plots of the transmission coefficients in Figs. 18, 19 and 23 , it is observed that the key factors affecting the wave transmission are the number of plates and the distance between them, but not their total lengths. This correlates with observations of Kohout and Meylan [33] who found out that attenuation coefficient in the marginal ice zones is a function of number of floes and is independent of average floe length.
In Case D[L 1 , L 2 ], we investigate the question of, if wave scattering can be manipulated by installing relatively smaller plates (wave breakers) of different lengths upwave the floating structure. In Fig. 24, the variations in reflection and transmission coefficients with λ/L for the long plate preceded by two relatively smaller plates of different lengths are shown and compared with a single plate case. Figure 24 demonstrates that such plate formations are more resistant to the incident waves than a single plate. Moreover, non-uniform size distribution in an array of wave breakers should provide better wave reflection than an array of equally sized plates, since different plates have different scattering properties, and thus, a wider wave spectrum can be covered. The order in which the wave breakers are installed plays an important role in transforming the wave motion.  Figure 25 shows the contour plots of reflection and transmission coefficients for two-parametric families of Cases C[L 1 , l] and D[L 1 , L 2 ], so that interplay between parameters can be studied. In contrast with Figs. 23-24, Fig. 25 gives more information about the wave transformation by different plate configurations presented in Cases C and D. Figure 25 allows to observe clearly the optimum pairs of parameters (l, L 1 ) and (L 1 , L 2 ) in Cases C and D, respectively, for one specific wavelength. Figure 25a shows that compact zones of maximum wave reflection are distributed uniformly in two-dimensional space of parameters L 1 and l. The diameter of these zones depends on the incoming wavelength. The plots of c R as a function of parameter L 1 with l fixed, or as a function of l with L 1 fixed, have the basic structure, consisting of a series of arched segments, dropping rapidly to zero at their ends. The same structure of reflection coefficient as a function of ice floe diameter has been observed in Meylan and Squire [41] and has been explained by the resonance effect, when the wave suits perfectly to the length of the floe. Figure 25b shows that for particular incoming wave properties, there is a compact region of parameters L 1 and L 2 where Case D[L 1 , L 2 ] reflects the waves most effectively. Comparing Cases C[L 1 , l] and D[L 1 , L 2 ] in Fig. 25, it may be concluded that breaking the structure into two parts can provide similar wave reflection as installation of two additional plates upwave.
From this analysis, it can be concluded that for any given incoming wavelength λ and plate properties L i , m i , D i (i = 1, . . . , N ) it is possible to optimize wave scattering by choosing appropriate number of plates N and the distances between them. Given the extremely low computational cost of the GN equations, the GN model makes a superior tool to study such nonlinear problems.

Conclusions
In this paper, the problem of interaction of nonlinear water waves with a set of N number of floating elastic sheets of arbitrary properties is solved by use of the Level I GN theory. The final scheme is stable and demonstrates high efficiency, performing calculation with multiple sheets in a short time period. One of the advantages of using the GN equations, when compared to existing methods based on the computational fluid dynamics, is that it allows to construct the velocity field and pressure distribution in the fluid domain at low computational cost. The velocity field analysis reveals the irregular character of the flow near the plate edges. The accuracy of the computed solution is validated by comparing with experimental data and other numerical predictions obtained by the use of alternative methods. The effect of multiple sheets is studied by wave scattering analysis on several basic models. The following main effects are observed: (i) the set of plates is more resistant to incident waves than a single plate of the same total length; (ii) increase in open water space inside the floating sheet system leads to decrease in wave transmission; (iii) the last element in the floating sheet system experiences the least stress and less likely to break under the wave impact; (iv) addition of the relatively smaller sheets upwave the floating structure changes the diffraction image, so that the waves of specific length are better reflected. These provide an insight into the problem of wave interaction with floating ice sheets, and particularly the wave effect on the ice sheets (whether it is strong enough to crack them) and the effect of ice sheets on the wave field. Moreover, these observations can contribute to the design of a floating wave filtering devices, protecting offshore structures from strong wave conditions, or to optimize the performance of floating wave absorbers. Since the number of parameters, describing the model, increases with increase in the elements constituting the model, the interplay between individual effects of the sheet properties can be explored thoroughly by the GN model developed in this study.