Flexural-gravity wave interaction with multiple vertical cylinders of arbitrary cross section frozen in an ice sheet 

The interaction problem of flexural-gravity wave with multiple vertical cylinders frozen in an ice sheet on the surface of water with finite water depth is considered. The linearized velocity potential theory is adopted for fluid flow, and the thin elastic plate model is applied for ice sheet deflection. Each cylinder is bottom-mounted, and the shape of its cross section can be arbitrary while remaining constant in the vertical direction. The velocity potential is expanded into an eigenfunction series in the vertical direction, which satisfies the boundary condition on the ice sheet automatically. The horizontal modes, which satisfy the Helmholtz equations, are then transformed into a series of boundary integral equations along the ice sheet edges or the intersection of the ice sheet with the cylinders. The problem is then solved numerically by imposing the ice sheet edge condition together with the impermeable condition on the cylinders. The solution is exact in the sense that the error is only due to numerical discretization and truncation. Computations are first carried out for single and multiple vertical circular cylinders, and good agreements are obtained with the semi-analytical solution. To resolve the difficulty of excessive computation at a large number of cylinders, the effect of the evanescent wave of a cylinder on those at large distance is ignored. This allows for the case of a large number of cylinders in different arrangements to be simulated. Extensive results


I. INTRODUCTION
The interaction of flexural-gravity wave with single and multiple structures has received great attention.In the cold environment, the water surface may become frozen to form an ice sheet from time to time.In this case, water wave motion may interact with the sheet.When the wave encounters a structure in its path, the nature of wave diffraction and the wave force on the structure will be different from those in open water.Therefore, a better understanding of their features is of practical importance for designing a safe and reliable structure as well as for protecting the environment.
Extensive research works have been carried out on free surface gravity wave/structure interactions.In general, the flow is governed by the Navier-Stokes equations with the fully nonlinear boundary conditions on deforming free surface.When flow is dominated by gravity wave and the viscosity effect is small, the velocity potential theory may be used.The nonlinear boundary conditions on the free surface may be linearized when the amplitude of wave is relatively small compared with its length and the characteristic length of the body.This theory has been widely adopted in the body/wave interaction and made significant impact to the engineering practices.A particular type of structure is the vertical cylinder, for which extensive work has been done because of its wide arrange of applications in offshore platforms.For a cylinder with circular cross section, Havelock 1 solved the problem in infinite water depth.MacCamy and Fuchs 2 considered the case of finite water depth, using the Fourier series in the circumferential direction, vertical mode expansion along the depth, and Hankel function in the radial direction.Explicit equations for wave force and moment on a single circular cylinder were derived.Using the first order solution of the linearized theory and based on the procedure of Lighthill 3 and Molin, 4 Taylor and Hung 5 obtained the second order force on the cylinder, which did not require the solution of the second order potential itself.By adopting the eigenfunction expansions and boundary element method (BEM), Chau and Taylor 6 obtained full second order solution, which could provide many more other detailed results, such as pressure distribution and free surface elevation.
In some structures, the cross section of a cylindrical component may not be always circular.The problem of an elliptic cylinder cross section was considered by Chen and Mei. 7The expansion in the vertical direction is the same as that for a circular cylinder.In the horizontal plane, the elliptic cylindrical coordinate system was used, and the velocity potential was expanded into a series of Mathieu functions together with Fourier series.Later, for the same problem, Williams 8 used two different methods.The first one expanded the Mathieu functions into Bessel functions based on the elliptic eccentricity parameter, and the other is the boundary element method using the Green function.Mansour et al. 9 also used numerical method for wave diffraction of waves by a vertical cylinder.Liu, Guo, and Li 10 adopted a method by expanding the radial functions of cylinder surface and Bessel functions into Fourier series with unknown coefficients.The coefficients were determined in terms of cylinder boundary condition using the Galerkin method.Dişib€ uy€ uk, Korobkin, and Yilmaz 11 adopted a method by introducing a small parameter based on the variation of the radius from a constant.
For multiple vertical cylinders, notable works include those by Wu and Price, 12 Taylor and Chau, 13 and Evans and Porter. 14Most of them are based on either boundary element method or semi-analytical method.Different from a single cylinder, strong interactions may happen between cylinders.In particular, Maniar and Newman 15 considered wave diffraction by a finite array of identical cylinders.As the number of cylinders increases, near-resonant modes became more obvious at some wavenumbers, around which large forces were observed."Near trapping" was also observed by Evans and Porter 16 for the cylinders being arranged in circular arrangements.
When the upper water surface is covered by an ice sheet, the wave motion is affected by the properties of both fluid and ice sheet.The coupling of gravity wave in fluid and flexural wave in ice sheet results in the flexural-gravity wave in the icy water. 17][23][24] The three dimensional work includes that on wave scattering in an ice sheet covered harbor, 25 a floating body near a semi-infinite ice sheet, 26 a body in polynya. 27In some circumstances, the free surface may be fully frozen.Then, the structures will be connected directly to the edges of the ice sheets.For a circular cylinder, Brocklehurst, Korobkin, and P ar au 28 obtained a semi-analytical solution through Weber transform, and later, the same problem was solved by Korobkin, Malenica, and Khabakhpasheva 29 with the vertical mode expansion method.In the work, the velocity potential and ice deflection were represented by Hankel functions in the radial direction.Based on the eigenfunction expansion and the Green's second identity, Ren, Wu, and Ji 30 derived semi-analytical solutions for diffraction of flexural-gravity wave by multiple circular cylinders.Cases of various arrangements of circular cylinders were considered.For non-circular one, Dişib€ uy€ uk, Korobkin, and Yılmaz 31 adopted a procedure similar to that used in free surface and cylinder interaction. 11n this paper, we shall consider the problem of flexural-gravity wave interaction with multi vertical cylinders of arbitrary shapes piercing through an ice sheet.The vertical mode expansion method is used for potential.The three dimensional problem is reduced to an infinite number of two-dimensional problems, which is truncated at a number where convergence has been achieved.Each two-dimensional problem is solved using the boundary element method (BEM).The method is extended from that developed by the authors previously. 27In theory, such a method can be used for any number of cylinders.However, when the number of cylinders increases, the computation becomes very inefficient.On the other hand, as the number of cylinders increases, the distance between many of them will be much larger than their typical dimension.Thus, the large spacing approximation is adopted, which has been widely used in applications on problems such as wave scattering in multiple wide polynya 32 and problem of wave energy devices in free surface. 33In this work, each cylinder will be first considered on its own.Its effect on other cylinders will be taken into account through the traveling wave only, while the effect of the evanescent waves will be ignored or be included up to only a few modes, as they decay exponentially.Through this method, the computational requirement will not increase significantly, even when there is a very large number of cylinders.
The main feature of this work is, therefore, that we use twodimensional BEM combined with vertical mode expansion to consider multiple cylinders of arbitrary cross section, together with an accurate approximation method for a large number of cylinders.The rest of this paper is organized as follows.The mathematical model based on the linear velocity potential theory for fluid flow and thin elastic plate theory for ice sheet is provided in Sec.II.The numerical solution procedure is given in Sec.III, where the eigenfunction expansion is used in the vertical direction, and the boundary integral equation is applied in the horizontal plane.The convergence and accuracy of the method are demonstrated through comparison with semi-analytical solution for a single and multiple vertical circular cylinders in Sec.IV.In Sec.V, numerical results for vertical cylinders of non-circular cross sections and different arrangements are presented, and the characteristics of the hydrodynamic force and ice deflection are discussed in detail.A large number of cylinders are considered by adopting an approximation method.Finally, conclusions are drawn in Sec.VI.

II. MATHEMATICAL MODEL
The problem of flexure-gravity wave interaction with multi vertical cylinders of arbitrary horizontal cross section is considered.As sketched in Fig. 1, the upper water surface is covered by an ice sheet of infinite extent.Each cylinder is bottom-mounted with its side surface being vertical.The water depth H is finite and constant.To describe the problem, a Cartesian coordinate system is introduced with its origin o being on the mean water surface.x axis and y axis of the system are in the horizontal plane, and z axis is vertically upward.
The fluid with constant density q w is assumed to be inviscid and homogeneous, and its flow to be irrotational.With the assumption that the wave amplitude is much smaller than its length and the typical body dimension, the linearized velocity potential theory can be used.For sinusoidal wave motion in time with radian frequency x, the total velocity potential can be written as Uðx; y; z; tÞ ¼ Re ixg/ðx; y; zÞe ixt Â Ã where g is the complex amplitude of incident flexural-gravity wave, i ¼ ffiffiffiffiffiffi À1 p , and / I and / D denote the incident velocity potential and diffracted velocity potential, respectively.
The conservation of mass requires that the velocity potential should satisfy the Laplace equation throughout the fluid, where The ice sheet extending to infinity is modeled as a thin elastic plate with uniform properties.The governing equation of the ice sheet deflection is where is the effective flexural rigidity with Young's modulus E, Poisson's ratio , and thickness h I , m I ¼ q I h I is the mass per unit area with density q I , and g denotes the acceleration due to gravity.The right hand side of Eq. ( 4) is, in fact, the hydrodynamic pressure.The impermeable kinematic condition on the ice sheet yields Similar to the velocity potential, the deflection of the ice sheet can be written as Combing Eqs. ( 4) and ( 5) together with Eqs. ( 1) and ( 6), the boundary condition on the ice sheet can be written as Without losing generality, the ice sheet edge is assumed to be clamped to each cylinder, which means that the defection and slope of the ice sheet at its edge are equal to zero, i.e., or where n is the two-dimensional unit normal vector of the ice sheet edge C. The impermeable boundary condition on the cylinders and seabed z ¼ ÀH reads The radiation condition requires the flexure-gravity wave propagating outward lim in which r 2 h ¼ x 2 þ y 2 , and j 0 denotes the purely positive real root of the dispersion equation

III. SOLUTION PROCEDURES
The aforementioned boundary value problem may be converted into a boundary integral equation over the body surface with a vertical orthogonal product. 34Here, it is solved through using the eigenfunction expansion in the vertical direction and a boundary integral equation in the horizontal plane.In the vertical direction, the diffracted velocity potential is expanded into eigenfunction series, i.e., where in which j m is the m th root of the dispersion equation for ice sheet in Eq. (12).It should be noted that the two complex roots with negative imaginary parts j À2 and j À1 are symmetric about the imaginary axis.j 0 is the purely positive real root, and j m (m ¼ 1; …; 1) are an infinite number of purely negative imaginary roots.Substituting Eq. ( 13) into Eq.( 2), we have By applying Green's identity over the boundary of the horizontal plane, we can convert this two-dimensional Helmholtz equation into an integral equation, where aðpÞ is the solid angle at field point p.Gðp; qÞ in Eq. ( 16) is Green's function, which can be written in terms of the zeroth order Hankel function of second kind where R is the horizontal distance between the field point pðx; y; zÞ and the source point qðn; g; fÞ.It should be noticed that G m ðp; qÞ in Eq. ( 17) satisfies the radiation condition in Eq. (11).In fact, apart from when m ¼ 0, G m ðp; qÞ decays exponentially at far field, while at m ¼ 0, the term in the brackets becomes zero because of the radiation condition.As a result, only the integral along the ice sheet edge C is kept in Eq. ( 16).
To satisfy the ice sheet edge condition ( 9) and the impermeable condition (10) on the cylinder, the orthogonal inner product for the eigenfunction w m may be used, 35 hw m ; w m i ¼ This gives hw m ; w m i ¼ 0 for m 6 ¼ m and hw m ; Applying the inner product to @/=@n and w m along the vertical surface of a cylinder, we have According to Eqs. ( 10) and ( 9), both the first and second terms on the right hand side of Eq. (20) equal to zero, which provides where is an unknown function and is to be determined.Following Li, Shi, and Wu, 27 the incident potential can be written as where Àij0ðx cos bþy sin bÞ ; (24)   and b is the incident wave angle.Substituting Eqs. ( 13) and ( 23) into (21), we have Invoking Eq. ( 24), Eq. ( 16) can be rewritten as Here, C includes the ice edge of all cylinders where } is the number of cylinders.To obtain numerical solution, C k of cylinder k is discretized into N k straight line segments, and all together, there will be N tot ¼ P N k (k ¼ 1; 2; …}) segments.The variables on each segment are assumed to be constant.Then, the boundary integral Eq. ( 26) can be discretized into in which the field point p i is the center of element i, and the source point q j of the Green function varies with l within element j.Equation ( 27) can be also written in the matrix form or where Through Eq. ( 28), we have where Substituting Eq. ( 33) into (13), truncating the series at m ¼ M À 3, or keeping only the first M terms, and invoking the clamped edge condition in Eq. ( 9), we have or After ½! is obtained through the solution of Eq. ( 36), Eq. ( 25) provides Equations ( 28), (36), and (37 Noticing there are M Â N tot unknowns u m , M Â N tot unknown @u m =@n, and N tot unknown ! at the centers of the segments, the total number of the unknown is the same as that of the equations.
After the diffraction potential is found, the horizontal force F ðkÞ on the k th cylinder can be obtained through Based on the kinematic condition ( 5), the deflection of ice sheet can be written as Together with Eqs. ( 2), (22), and (40), the vertical shear stress Q k can be derived as In such way, the shear force Q ðkÞ t along the cylinder edge can be easily obtained by

IV. VALIDATION OF THE METHODOLOGY
We first consider the case of a single circular cylinder with radius a.The semi-analytical solution of this case has been obtained by Brocklehurst, Korobkin, and P ar au 28 through the Weber transform and by Ren, Wu, and Ji 30 through the matched eigenfunction expansion.To verify the methodology, the same parameters as those in Ren, Wu, and Ji 30 are adopted.Radius of the circular cylinder is chosen as a ¼ 2 m, together with the ice density q I ¼ 917 kg m À3 , Young's modulus E ¼ 4:2 GPa, ice thickness h I ¼ 1:6 m, Poisson's ratio v ¼ 0:33, water depth H ¼ 350 m, and water density q w ¼ 1026 kgm À3 .Two sets of discretization for the case are taken, namely, (1) N ¼ 47 and M ¼ 20, (2) N ¼ 94 and M ¼ 30.In Fig. 2, the dimensional horizontal force F x on a single circular cylinder is shown against the j 0 H, where the semi-analytical solutions in Ren, Wu, and Ji 30 are also provided for comparison.It can be observed from Fig. 2(a) that there is no visible difference between the results from two sets of discretization, and the present results are in good agreement with the semi-analytical solutions of Ren, Wu, and Ji. 30 Furthermore, comparison for the dimensional vertical shear force on a circular cylinder in an ice sheet with a ¼ 3:5 m and h I ¼ 1:6 m is made with that Ren, Wu, and Ji. 30 The result is shown in Fig. 2(b), and a good agreement can be seen.
We then consider the case of four identical circular cylinders of equal radius arranged in square of length 2d.The centers of cylinders marked 1-4 are located at À ffiffi ffi . The radius a is taken as the characteristic length.The ice sheet thickness h I ¼ 0:1, the length of the square d ¼ 2, water depth H ¼ 10, fluid density q w ¼ 1025 kg=m 3 , Young's modulus E ¼ 5:0 GPa, Poisson's ratio v ¼ 0:3, and ice density q I ¼ 922:5 kg=m 3 , which gives L ¼ 4:5582 and m I ¼ 0:09.Figure 3 depicts the horizontal force F ðkÞ x , nondimensionalized by q w ga 3 Á ðg=aÞ on the k -th cylinder against j 0 a.From the figure, a good agreement with the semi-analytical solution of Ren, Wu, and Ji 30 can be seen, which further confirms the accuracy of the present numerical methodology for multiple cylinders in frozen sea.

V. NON-CIRCULAR CYLINDERS
After verification, the flexural-gravity wave interaction with vertical cylinders of more general shape is considered.Similar to previous works, the parameters of ice sheet and fluid are taken as E ¼ 5 GPa; ¼ 0:3; q I ¼ 922:5 kgm À3 ; q w ¼ 1025 kgm À3 ; g ¼ 9:80 ms À2 : Unless it is specified, the following numerical results are given in the nondimensional form based on three basic parameters, the characteristic length of the cylinder, density of water q w , and the acceleration due to gravity g.

A. Single cylinder
We first consider the case of a single cylinder.The geometry of the cylinder we taken is sketched in Fig. 4, which is a rectangle rounded at each corner by a quarter of a circle.2a and 2b denote the length and the width of the rectangle, respectively, and four quarter circles at the four corners have the same radius r. a is taken to be the characteristic length scale l.The other parameters are set to be b ¼ 1, h I ¼ 0:1, H ¼ 10, and b ¼ 0.
Figure 5(a) provides the horizontal force F x nondimensionalized by q w ga 3 Á ðg=aÞ on the single cylinder against j 0 a, and the corresponding results in open water are plotted in Fig. 5(b) as comparisons.It may be noticed that for open water with h I ¼ 0, j 0 represents the wavenumber for free surface gravity wave.Five sets of the corner radiuses are taken, i.e., r ¼ 0, 0.2, 0:4, 0:8, and 1:0.It is noted that r ¼ 0 refers to cylinder with rectangular cross section.Since the cylinder is symmetric with respect to the o À xz plane, the horizontal force F y equals zero when the incident wave angle b ¼ 0.
In Fig. 5, when j 0 a !0, the wavelength tends to infinity, and the variation of the wave over the cylinder is negligible and F x !0. As j 0 a increases, the dynamic effect becomes more obvious as can be seen in Fig. 5(b).For smaller r, or when the shape is closer to a rectangle, the force is larger.It is smallest for a circular shape or when r ¼ 1.The results at different r all reach a peak around a similar k 0 a, or j 0 a ¼ 0:18, and then decay at a similar rate within the range of j 0 a calculated.For the free surface in Fig. 5(b), at smaller j 0 a, smaller r gives the larger force.However, the curves cross each other around j 0 a ¼ 1:6, from where the curves are also close to each other within the range calculated.The results with ice cover and in the open water of course cannot be expected to be the same.In particular, it should be noted at same j 0 , and the frequencies in these two cases are different.
The shear force may sometimes be strong enough, which may pose a risk to integrity of the cylinder.Figure 6 displays the shear force Q t nondimensionalized by q w ga 3 Á ðg=aÞ, against the j 0 a.As in Fig. 5, five different r are considered.The ice sheet thickness is taken as h I ¼ 0:1.From the figure, we may find that when j 0 a tends to zero, the shearing force tends to a non-zero value.For each cylinder with different r, the value is not the same.In fact at j 0 !0, the ice sheet deflection due to the incident potential is w I ¼ @/ I =@z ! 1.For deflection due to the diffraction potential w d ¼ @/ d =@z, Eq. ( 7) gives where k À 4 ¼ q w g=L.The boundary conditions of w d give w d ¼ 0 at infinity, and @w d =@n ¼ 0 and w d ¼ À1 on the cylinder.We may consider a single circular cylinder in the polar system (q; h).w d is then a function of r only.This gives a solution in the form 36 where kei 0 , ker 0 , bei 0 , and ber 0 are the Kelvin functions of the zeroth order. 37Because of the condition at infinity, only the first two terms are kept.From the edge conditions at q ¼ a, we have (46) where v ¼ kei 1 ðk À aÞ þ ker 1 ðk À aÞ ½ kei 0 ðk À aÞ þ ker 1 ðk À aÞ À kei 1 ðk À aÞ ½ ker 0 ðk À aÞ: (47) The nondimensionalized shear force due to diffraction can be obtained through where the derivatives of Kevin functions have been obtained from Abramowitz and Stegun. 37For the case r ¼ 1:0 in Fig. 6, Eq. (48) gives jQ t;d j ¼ 27:7793, which agrees well with the result in the figure.In Fig. 6, the shear forces on cylinders with different r are very close to each other.When wavenumber increases from zero, Q t first decreases until it reaches a local trough at about j 0 a ¼ 0:36.It then increases mildly before it continues to decrease.

B. Cylinders arranged in square
Then, we shall consider flexural-gravity wave interactions with multiple cylinders.Four identical cylinders of the shape in Fig. 4 with a ¼ b, arranged in square shown in Fig. 7, are investigated in this section.The distances between the centers of cylinders 1 and 2 and between cylinders 2 and 3 are both 2d. a is taken as the characteristic length scale.
Figure 8 provides the horizontal force F ðkÞ x and F ðkÞ y on the k -th cylinder against j 0 a, with r ¼ 0, r ¼ 0:2, r ¼ 0:4, r ¼ 0:8, and r ¼ 1:0, nondimensionalized by q w ga 3 =ðg=aÞ.As a comparison, the results with r ¼ 0:2 in open water are also provided.The incident  When j 0 a is small, the forces on four cylinder are very close to each other and to those on the cylinder in open water, which is similar to the ice channel problem solved by Li, Wu, and Ren. 24Overall, the curves show an oscillatory behavior, and the reason for which has been extensively discussed by Li, Shi, and Wu 23 through detailed mathematical analysis.Principally, this is caused by the wave created by one cylinder, reached other cylinders, and reflected back.Peaks of the hydrodynamic force can be observed at j 0 a ¼ 1:573, 1.562, 1.498, 1.384, and 1.339 for r ¼ 0, 0.2, 0.4, 0.8, and 1.0, respectively.Away from the peaks, the hydrodynamic forces at different r follow the same trend.The result corresponding to the ice sheet can be seen very much different from that in open water.It should be noted that as j 0 increases, x 2 increases linearly at large j 0 in open water, while it increases at a rate j 5 0 in the case of ice sheet.1.5 j 0 a ¼ 2:904 j 0 a ¼ 2:904 j 0 a ¼ 2:904 j 0 a ¼ 2:904 2.0 j 0 a ¼ 1:562 j 0 a ¼ 1:528 j 0 a ¼ 1:528 j 0 a ¼ 1:528 3.0 j 0 a ¼ 0:778 j 0 a ¼ 0:753 j 0 a ¼ 0:778 j 0 a ¼ 0:778 It can be seen in Fig. 8 that the wave force on each cylinder in icy water becomes more oscillatory, compared with that in open water.Further than that magnitude of the force peak in icy water is generally higher than that in open water.occur around j 0 a ¼ 1:53, which is close to that in Fig. 8.However, the amplitudes of peaks are different in the two figures.
The effect of the distance between cylinders is then considered, with r ¼ 0:2 and b ¼ 0. Figure 10 provides the horizontal force F ðkÞ x and F ðkÞ y on the k -th cylinder against j 0 a, with three different distances to dimension ratio, namely, 2d=2a ¼ 1:5, 2d=2a ¼ 2, and 2d=2a ¼ 3.In Fig. 10(a), jF ð1Þ x j for a single cylinder in icy water is also shown, as a comparison.It can be observed from the figure that when the distance 2d between cylinders becomes larger, the hydrodynamic force becomes more oscillatory, and the position of main peak changes, as given in Table I.As the oscillation is very much due to the reflection by other cylinders, its period can be expected to be very much related to j 0 d, which is similar to the problem of a body in a polynya. 23In fact, as can be seen in Fig. 11, the gaps between local peaks are similar, although the peaks happen at different wavenumbers.
C. An approximate solution method for a large number of cylinders For multiple cylinders, with the number increases, direct solution based on the aforementioned procedure will become very timeconsuming or impractical.Here, procedure based on an approximation is proposed.For } cylinders, k ¼ 1; …; }, in Eq. ( 28), when both the field point p and the source point q in the Green function are located on the same cylinder k, we keep all the modes m ¼ À2; À1; ….We notice that m ¼ 0 corresponds to a traveling wave, while other modes correspond to evanescent waves.Thus, when the source point q is located on the other cylinders, especially when the distance between cylinders is large, we can keep only the mode of m ¼ 0 or a few more modes, which have not substantially decayed.These are moved to the right hand of equations for each cylinder.The solution of Eqs. ( 28), (36), and (37) can be obtained through iteration.In such a case, the whole solution of multiple cylinders in Sec.III can be simplified as the solutions of individual cylinders.In the previous solution procedure without approximation, the size of coefficient matrix of Eqs. ( 28), (36), and ( 37 With the approximate solution method, the matrix size is reduced into For identical cylinders, all the matrixes will be the same.The computing time will be dramatically reduced. We first consider cylinders arranged in an array and two arrays, sketched in Fig. 12.For convenience, all the cylinders are identical, 2a k ¼ 2a and r ¼ 0:2.The distances between two adjacent cylinders are 2d kÀ1 ¼ 2d.The ice sheet thickness is taken as h I ¼ 0:1, and the incident wave angle as b ¼ 0. Figures 13 and 14, respectively, present the horizontal force on the four cylinders in an array and two arrays against nondimensional wavenumber j 0 d=p with 2d=2a ¼ 10.In the figures, both approximate solution with only m ¼ 0 kept and exact numerical solution based on Eqs. ( 27), (36), and (37) are given.It can be seen that there are no distinguishable discrepancies between the results from these two solutions, which indicates the accuracy and effectiveness of the approximate method when 2d=2a becomes larger.
The approximate solution method is then applied to two arrays of cylinders arranged in Fig. 12(b), with large number } ¼ 2 Â 25.Due to symmetry at b ¼ 0, the forces on the corresponding cylinders in the two arrays are the same.Figure 15 provides the horizontal forces on the k th cylinders against j 0 d=p, where k ¼ 1, 7, 13, 19, and 25.The corresponding shear forces are exhibited in Fig. 16.It can be seen in Fig. 15(a) that there are some local peaks occurred around j 0 d=p ¼ n=2, n ¼ 1; 2; 3; ….The distribution of the force on each cylinder at three of these frequencies j 0 d=p ¼ n=2 (n ¼ 1; 2; 4) is shown in Fig. 17 for } ¼ 2 Â 25 and } ¼ 2 Â 26.It can be seen from Figs. 15-17 that the forces acting on the first cylinder are much larger than those on other cylinders, in general.However, around j 0 d=p ¼ 1:53 where a sharp peak exists, the variation of the force on the middle cylinder, or k ¼ 13, is largest.The behavior of the forces on a large array of cylinders in the free surface wave has been extensively discussed by Newman, 15 and large force can be observed at certain frequencies related to resonance and trapped modes.The result here exhibits some similar behavior.

VI. CONCLUSIONS
A method has been proposed for the problem of flexural-gravity wave interaction with vertical cylinders of arbitrary shape frozen in an ice sheet.It combines the eigenfunction expansions in the vertical direction and the boundary element method in the horizontal plane.The verification of the method is conducted through single and multiple circular cylinders, and results are an excellent agreement with the semi-analytical solution.Extensive calculations are then made for multiple non-circular cylinders.When the number of the cylinder increases, an approximation is introduced, which is efficient and accurate.Based on these, the following conclusions can be drawn: (1) For cylinder in non-circular shape with a rectangle rounded at each corner r by a quarter of a circle, the forces are larger with smaller r.It is smallest for r ¼ 1, i.e., a circular shape.The horizontal forces at different r all reach a peak around a similar k 0 a and then decrease.On the other hand, the shear forces on cylinders with different r are very close to each other.When r ¼ 1, an analytical solution can be obtained and shows that the shear force tends to a non-zero value when k 0 a !0. (2) When the four non-circular cylinders arranged in square, wave force on each cylinder in icy water becomes more oscillatory, compared to that in open water.Peaks of the forces can be observed to decay when r increases.Away from the peaks, the hydrodynamic forces at different r follow the same trend.(3) When the number of the body becomes large, an approximate solution method introduced is found accurate and efficient.
From the obtained results, large force variations with frequency and with the location of the cylinder have been observed, similar to those observed in the free surface problem.

FIG. 1 .
FIG. 1. Sketch of the problem.(a) Side view and (b) bird view.
wave angle is chosen as b ¼ 0. Because of symmetry, we have jF ð3Þ Thus, the forces with k ¼ 3 and k ¼ 4 are omitted in figure.

Figure 9
Figure9provides the horizontal force F ðkÞ x and F ðkÞ y on the k -th (j ¼ 1; …; 4) cylinder against j 0 a, at quarter sea of b ¼ 0:25p.The radius of the cylinder corner is taken as r ¼ 0:2, and the ice sheet thickness is set to be h I ¼ 0:1 and 2d=2a ¼ 2. Due to the arrangement Moreover, as shown in Fig. 9, the positions of main peak of the forces F ðkÞ x and F ðkÞ y

TABLE I .
The wavenumber of the main peak of horizontal forces for cylinders with different distances.