Diffraction of cnoidal waves by vertical cylinders in shallow water

Diffraction of nonlinear waves by single or multiple in-line vertical cylinders in shallow water is studied by use of different nonlinear, shallow-water wave theories. The fixed, in-line, vertical circular cylinders extend from the free surface to the seafloor and are located in a row parallel to the incident wave direction. The wave–structure interaction problem is studied by use of the nonlinear generalized Boussinesq equations, the Green–Naghdi shallow-water wave equations, and the linearized version of the shallow-water wave equations. The wave-induced force and moment of the Green–Naghdi and the Boussinesq equations are presented when the incoming waves are cnoidal, and the forces are compared with the experimental data when available. Results of the linearized equations are compared with the nonlinear results. It is observed that nonlinearity is very important in the calculation of the wave loads on circular cylinders in shallow water. The variation of wave loads with wave height, wavelength and the spacing between cylinders is studied. Effect of the neighboring cylinders, and the shielding effect of upwave cylinders on the wave-induced loads on downwave cylinders are discussed.

Havelock [29]. Scattering of waves for very long wavelength (solitary wave) by a cylindrical object (island) was first solved by Omer and Hall [54].
Isaacson [35] used the integral equation method to investigate the interaction of long waves with isolated cylinders. Abul-Azm and Williams [1] studied some of these effects by an approximate method. Chau and Taylor [8] used a second-order analysis and demonstrated that the second-order effects caused the hydrodynamic pressure to penetrate deeper into the water column. They also showed that the second-order diffraction effects in shallow water are substantially different from that in deep water. Malenica and Molin [42] used a third-order analysis to investigate the wave-cylinder interaction and showed that the third-order effects could alter the triple frequency component of the diffraction potential. Kriebel [39,40] studied the second-order diffraction of waves by a vertical cylinder and compared the calculations with experimental data, and Newman [51,52] investigated the diffraction of waves by a vertical cylinder up to third order through perturbation analysis.
Nonlinear wave loads on a slender vertical cylinder is studied by Faltinsen et al. [22] by assuming an inviscid and incompressible fluid and irrotational flow and solving nonlinear velocity potentials. Ferrant et al. [24] investigated the nonlinear wave scattering and diffraction forces on three-dimensional bodies, including a vertical cylinder, by use of the boundary-element method. Bai and Taylor [3] used the domain decomposition method to study the fully nonlinear wave diffraction of linear waves by a vertical cylinder in monochromatic waves as well as focused, group waves. Kim et al. [36,37] used the finite-element method to study the timedomain propagation of Stokes waves and their interaction with vertical cylinders and compared the run-up with the stream function theory of Ferrant [23]. Nonlinear loads on vertical cylinders and floating bodies in irregular waves are studied by Sclavounos [61,62] through the time derivative of the fluid impulse.
Majority of the studies on wave interaction with vertical cylinders have focused on deep water applications. In shallow water, and for large waves where nonlinearity and dispersion effects are significant, the use of shallow-water wave equations is inevitable, see, e.g., Kim et al. [37], Swan and Sheikh [66], for the effects of nonlinearity and wave steepness on wave loads on a vertical cylinder. Studies on solitary and cnoidal wave loads on vertical cylinders are very limited, see Apelt and Piorewicz [2], Wang et al. [71], Yang and Ertekin [80], Wang and Jiang [70], Yates and Wang [81], Zhong and Wang [86]. Recently, with the development of marine renewable energy devices, further attention is paid to nonlinear wave diffraction in shallow water by vertical cylinders, see, e.g., Paulsen et al. [55], who solved the Navier-Stokes equations (using OpenFOAM) to study the wave forces on vertical cylinders in water of finite depth. Similar study was performed by Mo et al. [45,46] by ignoring the effect of viscosity and solving Euler's equations.
Yet another approach to study the problem of interaction of nonlinear waves with vertical cylinders in shallow water is by use of nonlinear, shallow-water wave theories, for example, the Green-Naghdi (GN) equations (Green and Naghdi [25]) and the Boussinesq equations (Wu [77]). The GN equations and the Boussinesq equations both have been applied to many hydrodynamic problems concerning diffraction loads on structures in shallow water. The Boussinesq equations were used by, for example, Wang et al. [71] to study solitary wave scattering by a vertical cylinder, by Nadiga et al. [47] to study wave diffraction due to an uneven bottom, and by Schmitt et al. [60] to study wave loads on a bottom hinged oscillating wave converter. The Level I Green-Naghdi equations are used, for example, to determine the hydroelastic response of VLFS to solitary and cnoidal waves by Ertekin and Kim [14], Ertekin et al. [20], Xia et al. [78,79], Ertekin and Xia [16], and by Demirbilek and Webster [10], Ertekin and Sundararaghavan [15] to study refraction and diffraction of nonlinear waves in coastal waters and to study solitary and cnoidal wave loads on horizontal decks by Hayatdavoodi and Ertekin [30][31][32], Hayatdavoodi et al. [33]. Further discussion on the recent developments on the Green-Naghdi equations, with comparisons between high-level equations and laboratory experiments, can be found in Kim et al. [38], Ertekin et al. [21], Zhao et al. [83][84][85], among others. Theoretical models based on the Green-Naghdi and the Boussinesq equations were recently developed by Neill et al. [50] to study the interaction of solitary wave with multiple in-line, vertical cylinders. Good agreement was observed between the Green-Naghdi and the Boussinesq equations solitary wave results and laboratory measurements. Cnoidal wave run-up on a single vertical cylinder was studied by Zhang and Teng [82] by solving the Boussinesq equations by use of finite-element method.
Our objective in this work is to study the problem of diffraction of cnoidal waves by single or multiple inline vertical cylinders of finite diameter in shallow water. We will use the Level I GN equations, the generalized Boussinesq equations, and the linearized version of these shallow-water wave equations, to numerically solve the initial-boundary-value problem of cnoidal waves interaction with vertical cylinders and to determine the horizontal force and overturning moment on the cylinders in shallow water. Different cnoidal waves can be generated by use of different shallow-water wave equations, and hence, wave loads on the structures may be different. The effect of nonlinearity on the wave-induced loads on the vertical cylinders is discussed by comparing the nonlinear results with the linear solution.
The GN and Boussinesq nonlinear, shallow-water wave theories, along with the linearized version of the equations, are introduced next. This is followed by a discussion of the numerical modeling of the problem. Results of cnoidal wave diffraction by a single and multiple in-line vertical cylinders are discussed in Sect. 5. In presenting the results, and comparisons with other predictions and experiments when available, we first start with the single, isolated cylinder case. Results for multiple cylinders, namely two and three in-line cylinders, and the effect of neighboring cylinders on the wave loads are discussed after. The paper is closed with some concluding remarks.

The shallow-water wave theories
A right-handed Cartesian coordinate system, whose origin is at the still-water level (SWL) and on the entrance boundary where the numerical wavemaker is located, is used. The incident waves propagate in the positive x-direction, along the line of symmetry. The y-direction is parallel to the entrance boundary, and the z-direction is vertical, with positive z pointing to the opposite of the gravitational acceleration direction. A schematic of the theoretical wave tank is shown in Fig 1. The fixed vertical cylinders are of constant circular cross section and extend from the free surface to the flat and stationary sea floor. The problem is symmetric with respect to the line x-direction that passes through the cylinder center, and hence, only one half of the physical region is considered for the computations. The symmetric numerical tank will be discussed in Sect. 3.4. In this study, the fluid is initially quiescent and at time t = 0, the wavemaker starts generating cnoidal waves. The initial-boundary-value problem of wave interaction with the vertical cylinders will be solved by both the GN and the gB equations, or their linearized version. A summary of the GN and the gB equations are presented here. Further details of these models can be found in Neill et al. [50], who studied the problem of solitary wave interaction with multiple, in-line vertical cylinders by use of the GN and the gB equations. Unless otherwise mentioned, all variables are dimensionless using (ρ, g, h) as a dimensionally independent set, where ρ is the mass density of the fluid, g is the gravitational acceleration and h is the water depth.

The Green-Naghdi equations
The nonlinear Green-Naghdi (GN) equations for water waves are first given by Green and Naghdi [26] based on the theory of directed fluid sheets and by Green and Naghdi [25,27], through the exact equations for propagation of waves in a homogeneous, incompressible, inviscid fluid. In these equations, irrotationality of the flow is not required, although the equations can also be derived for irrotational flows. Using a direct approach, Green and Naghdi [26] postulated the integral balance equations of conservation of mass, momentum, moment of momentum (or director momentum), and energy to a deformable fluid sheet. The resulting GN equations satisfy the nonlinear boundary conditions exactly, identically conserve mass, and exactly satisfy the integrated momentum equations in an average sense. Unlike the KdV equations, the GN equations are invariant under Galilean transformation.
The GN equations are classified based on the level of the polynomial functions used to prescribe the distribution of the vertical velocity along the water column. The Level I GN equations, also known as the restricted theory, are a specialized version of the more general results given by Green and Naghdi [28] and assume a linear distribution of the vertical velocity along the water column. This assumption, which is the only one made about the kinematics of the fluid sheet, results in the horizontal velocities being invariant in the vertical direction for an incompressible fluid. The resultant Level I GN equations, hence, are mostly applicable to the propagation of fairly long waves in shallow water. High-level GN equations, applicable to wave propagation in any water depth, are obtained by assuming higher-order polynomials for the distribution of the vertical velocity along the water column. Among others, see Shields and Webster [64], Webster and Kim [72], Demirbilek and Webster [9], and more recently Zhao et al. [83][84][85], for discussion on high-level GN equations.
The Level I GN equations are given here in the same form first presented by Ertekin [11] and Ertekin et al. [17]. In this form of the equations, the fluid is inviscid, and the flow is incompressible, but irrotationality is not required. The free surface, ζ(x, y, t), is measured from the SWL. The Level I GN equations as used here are given by the conservation of mass and momentum in dimensionless form as where α(x, y, t) is the elevation of the bottom of the fluid sheet, and superposed dot is the material time derivative. A double superposed dot denotes the second material time derivative, and ∇ is the gradient vector operator. The subscripts x, y and t indicate partial differentiation with respect to the variables. P I is the integrated (over the water column) pressure,p is the pressure on the bottom curve (α), andp(x, y, t) is the pressure on the top curve of the fluid sheet. The velocity field is defined by V = ue 1 + ve 2 + we 3 , where e 1 , e 2 and e 3 are the unit base vectors in the x, y and z directions, respectively. The horizontal velocity components, u and v, do not depend on the vertical z coordinate in the Level I GN equations.
In the absence of any perturbation expansion and scaling parameter, the order of error of the GN equations cannot be defined. Applicability and accuracy of the GN equations for various wave conditions must be determined by comparison with laboratory experiments. See, for example, Ertekin [12] and Webster and Wehausen [73] for further discussion on this.
In this work, the GN equations are applied to the problem of wave interaction with vertical cylinders following the same approach discussed in Neill et al. [50].

The generalized Boussinesq equations
We also use the Boussinesq equations to study the problem of interaction of nonlinear, shallow-water waves with vertical cylinders. For a homogeneous, incompressible and inviscid fluid, and unlike the GN equations, the Boussinesq equations are derived by assuming that the flow is irrotational. Hence, the velocity field is described by the velocity potential φ. A fundamental difference between the Boussinesq equations, of any form, and the GN equations, is the use of perturbation theory in deriving the equations. In the Boussinesq equations, the unknown velocity potential φ and surface elevation ζ are expanded over the linear solutions by use of scaling parameters, see Boussinesq [4]. Hence, the boundary conditions are satisfied approximately by the Boussinesq equations. Similarly, the conservation of momentum is satisfied only approximately.
Various versions of Boussinesq equations are discussed in the literature. Here, we use the generalized Boussinesq (gB) equations in the form first derived by Wu [77] for propagation of waves in constant water depth. See, for example, Nwogu [53] for an alternative form of the Boussinesq equations. The gB equations are given by the conservation of mass and momentum in dimensionless form by where φ is the layer-mean velocity potential. In this form of the equations, the atmospheric pressure is assumed zero. In deriving the above gB equations, it is assumed that σ = H/ h and = h/λ, the nonlinearity and dispersion terms, respectively, are small, where H is the wave height and λ the wavelength. Moreover, Hence, the gB equations are mostly applicable when the Ursell parameter, U r = σ/μ 2 , is of O(1). In the gB equations, the momentum equation (Eq. (8)) is obtained using a perturbation method, and hence, the conservation of momentum is satisfied only approximately. The error is of order (σ 4 , σ 2 2 ), see Wu [77] for details.
Further details about the gB equations as applied to the problem of nonlinear wave interaction with vertical cylinders can be found in Neill et al. [50].

The linear equations
We also use the linearized version of the GN and the gB equations in this study to assess the nonlinear effects. The linearized equations are obtained by neglecting the nonlinear terms in the equations, after fully expanding the material derivative terms. The linearized versions of the GN equations and the gB equations are identical as shown by Ertekin [11]. Naghdi [48] earlier showed that the GN equations also reduce to a number of other equations, including the one obtained originally by Boussinesq [4], when the equations are linearized.
In dimensionless form, after removing the nonlinear terms, the linearized version of the GN momentum equation, Eq. (2), can be written as Similarly, substituting u = ∂φ/∂x and v = ∂φ/∂y in Eq. (8) and removing the nonlinear terms gives as the order of mixed derivatives is interchangeable. Similar arguments can be made for the linear momentum equation in the y direction. Therefore, the linearized forms of the gB and GN equations are the same, and in dimensionless form, are given by

Numerical modeling
The nonlinear GN and gB equations, and the linearized version of these equations, are solved numerically to study the wave-structure interaction problem. In this section, we shall discuss the numerical modeling and intorduce the numerical wave tank used to study the problem.

Initial and boundary conditions
The fluid is initially quiescent and at time t = 0, generation of the cnoidal waves at the upwave boundary begins. The waves propagate in the positive x direction toward the open boundary. The seafloor, at the vicinity of the cylinders, is flat and stationary. The top surface is free, with uniform atmospheric pressure distribution. The seafloor no-flux boundary condition and the free surface kinematic and dynamic boundary conditions are embedded in the derivations of the equations, and hence, they do not need to be considered here. The free surface kinematic and dynamic boundary conditions are satisfied exactly by the GN equations and satisfied approximately by the gB equations and the linear equations. Details on how these conditions are embedded in the derivation of the equations can be found in Green and Naghdi [26] and Wu [77].
There are different types of cnoidal waves that are appropriate for nonlinear wave propagation in shallow waters, see, for example, Wiegel [75] and Whitham [74]. In comparison with linear waves, cnoidal waves have flatter troughs and peaked crests as a result of the nonlinearity. In some cases, semi-analytic and in other cases numerical solutions of cnoidal equations can be obtained by writing the governing shallow-water wave equation in Galilean coordinates and then solve the resulting ordinary differential equation by also imposing the periodicity condition. In the infinite wavelength limit, a cnoidal wave becomes a solitary wave.
In the following sections, the wavemaker solutions used in solving the three different governing equations employed in this work is introduced. In all periodic wave cases, however, it is necessary to modulate the waves at the wavemaker and around the initial time so that unrealistic pressure loads are not imposed on the fluid. This is because at time t = 0, the quiescent water condition is imposed throughout the domain and it is necessary to modulate (or ramp) the waves initially. This is accomplished by defining the modulated surface elevation by where e is the modulation parameter and L and T are the wavelength and wave period, respectively. Smaller value of e corresponds to longer and stronger modulation of the wave front. In this work, e = 0.5 is used. Velocity is determined by use of the exact conservation of mass equation, and hence, it is modulated directly through the modulation of the surface elevation ζ.

GN cnoidal wavemaker
The cnoidal wave solution used as the wavemaker in the case of GN equations is provided by Sun [65] and Ertekin and Becker [13] in a semi-closed form when the water depth is constant at the vicinity of the wavemaker. Let us define, x = x − ct, where c is the phase speed of waves, i.e., the constant speed of the right-moving coordinate. A periodic solution (dimensionless) to the GN equations for a wave of height H can be written in the moving coordinate as where the constant phase speed c is given by and where k 2 = H/(ζ 3 − ζ 1 ), and Cn is the Jacobian elliptic cosine function, and K and E are the complete elliptic integrals of the first and second kind, respectively. The wavelength L can be calculated using the following dispersion relation: Equations (15) and (18) are solved iteratively (see Ertekin and Becker [13]). We first assume a value for k and then iterate Eq. (18) until we find its root at the nth iteration, k (n) . The elliptic integrals are calculated by a polynomial approximation as described by Ertekin et al. [19]. Within the accuracy of the numerical calculations, the resultant waves satisfy the Level I GN equations exactly. Once the surface elevation, ζ , is determined, the velocity can be calculated from the conservation of mass equation in the moving coordinates:

gB cnoidal wavemaker
Following Qian [56,57], the cnoidal wave solution of the gB equations is obtained numerically in the x − z plane only (see also Teng and Wu [67]). By use of the definition of velocity potential, the layer-mean velocity potential φ can be deduced from the gB equations, Eqs. (7) and (8). In dimensionless form, the resultant gB equations read (see Ertekin et al. [18]) where is the 2-D Laplacian on the horizontal plane. By combining the gB equations Eqs. (19), (20) into a single nonlinear, ordinary differential equation in the moving (or wave) coordinate, x = x − ct, where c is the phase speed, we obtain where Fr = c/ √ gh is the depth Froude number, and C and D are integration constants. The solution must also satisfy the periodicity conditions: where A C is the crest amplitude and A T is the trough amplitude of the cnoidal wave, and therefore, the wave height is H = A C + A T , and L is the wavelength. Equation (22) requires an iterative solution: Given Fr = L/T, where T is the wave period, and wave height H, A k C is estimated, where k is the iteration counter, and A k T and C and D are calculated. Equation (22) is then solved by the fourth-order Runge Kutta method and A k+1 The iterations continue until the solution is converged within a specified tolerance, and the wavelength is then determined by L = 2x L .

Linear wavemaker
The long-wave solution used for the linear wave problem is the sinusoidal wave solution, where k = 2π/L is the wave number and ω = 2π/T is the wave frequency. Unlike small-amplitude linear waves, the shallow-water linear waves, as used here, do not follow the ω 2 = ghk 2 relation, and a specific dispersion relation must be defined.
To obtain the dispersion relation governing the 2-D linearized equations, Eqs. (11), (12), we substitute the dimensional forms of Eq. (11) in (12) and set v = 0 (velocity in y direction) to obtain Substituting Eq. (25) in (26), results in The wave speed, c, therefore, is in contrast with the linear wave phase speed for long waves, c = √ gh, and even with the two-term Taylor series expansion of it given by Mei [44], p. 30, as Finally, substituting Eq. (25) (for v = 0) in Eq. (11) results in the horizontal component of the particle velocity generated by the linear wavemaker as

Open boundary
An absorbtion boundary is used on the downwave side of the numerical wave tank to eliminate, as much as possible, the reflection of waves back into the domain. Potentially, this will alow minimizing the length of the computational domain. For this, the values of surface elevation and velocities must be known at the open boundary, the downwave end of the tank. However, due to the diffraction of the waves by the cylinders, there is no a priori knowledge of these variables at the open boundary. Previous works of Wu and Wu [76] and Ertekin [11] showed that the relatively simple Orlanski's condition with constant phase speed c = ± √ gh can be used to approximate the value of the variables and prevents significant reflections from the open boundary. This open-boundary condition, used in this work, is given by where Ω may be ζ(x, t) or u(x, t) at the downwave boundary. The wave surface elevation is recorded and monitored by numerical wave gauges throughout the wave tank. It is observed that the open-boundary condition utilized here works well with negligible reflections.

Grid generation
To facilitate the use of finite-difference methods in solving the governing equations in a Cartesian coordinate system, the physical domain, including the curved body boundaries, is mapped into a regular rectangular computational domain. The transformation is carried out by use of an elliptic grid generation method, where there is a one-to-one mapping between the physical and the computational domains; see, for example, Thompson et al. [68]. Although there is one or more bodies (cylinders) in the physical region, the domain can be reduced to a single, simply connected region by use of the symmetry of the problem. The axis of symmetry and cylinder surface of the physical problem map to a single side of the transformed domain. An example of the physical and transformed domains is shown in Fig. 2 for a single cylinder. Shown in this figure, the upwave boundary is where the numerical wavemaker is located and the downwave boundary is the open, or absorbing, boundary to prevent reflections as much as possible. On the symmetry, far wall and the cylinder boundaries, the normal component of the fluid velocities must vanish but we allow the tangential component as the fluid is assumed to be inviscid.
We note that with this grid generation approach, N number of vertical cylinders lined in a row parallel to the incident wave can be considered. For offset distribution of vertical cylinders, a modification of the grid generation approach is required, while the problem formulation remains the same. We only consider a row of in-line vertical cylinders in this study.

Boundary-fitted coordinates and transformed equations
An elliptical grid generation technique is used to relate the physical coordinates x(ξ, η) and y(ξ, η) to the curvilinear coordinates ξ(x, y) and η(x, y). For each node, x(i, j) and y(i, j), in the physical domain, there is a corresponding point, ξ(i, j) and η(i, j), in the computational domain. The governing equations for the elliptical generation system are found through the minimization of the Euler integral: This leads to the elliptical equations These equations are the governing equations for the elliptical grid transformation and ensure a one-to-one transformation between the physical and the computational domains since the determinant of the Jacobian is required not to vanish. The curvilinear coordinates ξ(x, y) and η(x, y) must be determined inside the physical region.
To solve Eq. (32) on the transformed domain, the derivatives in terms of x and y must be transformed into derivatives in terms of ξ and η. This is accomplished by the repeated application of the chain rule of differentiation, i.e., and similarly for ∂ f /∂ y, where f is a generic function. This results in the following equations: where For details of the derivation of Eqs. (33) and (34), see, for example, Thompson et al. [69]. Equations (33) and (34) provide the location of each node in the physical region corresponding to a node in the computational region. The transformed Eqs. (33) and (34), are solved iteratively using a finite-difference technique.
The derivatives in terms of x and y are then expressed in terms of ξ and η. Some of these relationships are listed below: where f is a generic variable and J is the Jacobian of the transformation, i.e., The transformation of the governing equations is explained in greater detail in Qian [56].

Grid orthogonality
The accuracy of a numerical model is significantly influenced by the orthogonality of the grid lines to the boundary. As the grid lines become less orthogonal, generally the numerical solution becomes less accurate. This inaccuracy can also lead to instability. To minimize this problem, we use grid attraction to maximize orthogonality by rather using the Poisson equation (see also, Thompson et al. [69]): where P and Q are control functions which determine the amount of grid attraction at each point and is the 2-D Laplacian. For positive values of P and Q, the lines of constant ξ and η are attracted; negative values repel these lines. Following the recommendation of Thompson et al. [68], P and Q are chosen to be exponential functions of the following form: whereā andc are the line attraction coefficients, b and d are the decay coefficients, and sgn is the sign of the difference in parentheses. This allows either line attraction or repulsion in either the ξ or η direction. After transforming the derivatives from the physical (x, y) plane to the computational (ξ, η) plane by the repeated application of the chain rule of differentiation, Eq. (40) become where J = x ξ y η − x η y ξ is the Jacobian of the transformation that should not vanish in the domain, and α, β and γ are defined by Eq. (35).

Stagnation points
The stagnation points are located where the line of symmetry bisects the cylinders. These locations are the most upwave and downwave locations on each cylinder; the points B and C in Fig. 2. The use of numerical grid generation requires that the grid lines be approximately orthogonal along the boundaries. Since there are 90 • bends at these points, the incident angle is undefined and the grid lines cannot be orthogonal. If the same transformation techniques used throughout the region are applied here, erroneous results would be produced. By definition, the velocities are zero at the stagnation points. For the problem to be symmetric, the flow through the line of symmetry, in the y direction, must be zero. Therefore, the y-component of velocity, v, at the stagnation points must be zero. On the surface of each fixed cylinder, the no-flux boundary condition requires the normal velocity to be zero also. At the stagnation points, this requires that the x-component, u, be equal to zero. By setting both velocities equal to zero at the stagnation points, the analysis can be significantly simplified. The surface elevation, however, must still be analyzed at these stagnation points.
At the stagnation point locations, the governing equations can be solved for the surface elevation without transforming into the computational domain. This boundary corresponds to the j = 1 line or constant η = 0 line. At point B, the line initially lies along the x axis and then, after the 90 • turn, lies approximately along the y axis. Along this boundary line, the grid points are evenly spaced. This line is used to create a Cartesian coordinate system, with an origin at the stagnation point. At each stagnation point, the x − y physical coordinate system is used. Conventional finite-difference formulas are applied without the transformation of the derivatives to the computational plane. The stagnation points are analyzed in the physical plane, while the rest of the region is analyzed in the computational plan. This method is also recommended by Thompson et al. [69].

Numerical solution
Second-order central difference schemes are used to solve the governing equations in the spatial domain. Fictitious points are used at the domain boundaries to carry out the spatial derivation, see Roddier [58] and Roddier and Ertekin [59] for details. The second-order, modified Euler method is used for time marching. This two-step method has been used successfully by Ertekin [11] and Hayatdavoodi and Ertekin [32], among others, in similar problems.
The successive over-relaxation method (SOR) is used to iteratively solve Eq. (40). The rate of convergence of the iterations is determined through the acceleration parameter following the method described by Burden and Faires [5]. Solving Eq. (40) by the SOR method yields a solution for each point in the physical domain x(i, j) and y(i, j), where i, j are the nodal points, with a corresponding point in the computational domain, ξ(i, j) and η(i, j). Therefore, initial guesses for x(i, j) and y(i, j) are required. We used a weighted average of the boundary points for the initial guesses as suggested by Thompson et al. [68] and used by Qian [56] and Ertekin et al. [18].
The Courant condition is satisfied in the numerical solution of the equations considered here as Von Neumann stability analysis done on the linear equations dictate that t < x or y.
It is observed that the results may experience saw-teeth noise, mainly due to the use of the central difference scheme. This is not physical and also may cause instability in the numerical solution. Hence, numerical filtering of the result is necessary. Shown by Ertekin et al. [17] for similar problems as that studied here, a combination of the five-(second-order) and seven-point (third-order) linear filtering schemes, developed by Shapiro [63], would ensure stability, while not altering the shape of the incident wave or the results. This includes the use of a third-order filtering in the direction normal to the prevailing wave crests, the ξ direction, and a second-order filtering parallel to the wave crests, the η direction. These filtering schemes read as where f is a generic variable that can represent ζ, u or v.

Error monitoring
The conservation of mass is monitored in the domain to ensure that the computations, and the open-boundary condition does not violate the law. The change in total mass, δ M, inside the physical domain is monitored by: where A C denotes the dimensional area of the computational domain in the horizontal plane and U W and DW represent the open boundaries on the upwave and downwave sides. The highest relative percentage change in mass determined in this study was less than 0.13%. Example error calculations and further discussion on the error monitoring of the numerical model used in this study can be found in Neill [49] and Neill et al. [50].

The Green-Naghdi equations
Closed-form equations for the pressure on the bottom surface of the fluid sheet (p, Eq. (6)) and the integrated pressure (P I , Eq. (5)) of the GN equations were first presented by Ertekin [11]. The wave-induced force on each cylinder can be determined by integrating the instantaneous integrated pressure, Eq. (5), around the cylinder.
In deriving the GN equations, however, the conservation of momentum and director momentum are satisfied in an averaged sense, and information about pressure distribution at every point of space is lost. In this study, the assumption made by Neill et al. [50] for pressure distribution over the water column is followed, i.e., that the pressure varies linearly from the free surface to the seafloor. Hence, is the expression for pressure at any point in the space.
The overturning moment (about the y axis) on the cylinder can then be obtained by multiplying Eq. (46) by the moment arm (with respect to the seafloor) and integrating over the circumference of the cylinder. The GN moment on a vertical cylinder can be obtained by The error associated with the assumption of linear pressure distribution in the vertical direction can be measured; see Neill et al. [50] for details. In this study, the maximum error, for cnoidal wave diffraction by vertical cylinders, is 0.29%. This is a small and reasonable error. In Sect. 5, further discussion on the pressure distribution over the water column by the GN equations will be provided.

The gB equations
The gB pressure distribution over the vertical direction in dimensionless form is given by (Wu [77]) The force on the vertical cylinder can be found by integrating Eq. (48) in the vertical direction, i.e., Similarly, the moment about the y axis can be found by multiplying Eq. (48) by the moment arm (with respect to the seafloor) and integrating it over the water column. The gB moment on a vertical cylinder can be obtained by

Linear equations
The depth-integrated pressure equation can be found by linearizing either the depth-integrated gB pressure equation, Eq. (49), or the GN integrated pressure equation, Eq. (5). Note thatζ when linearized, and hence, the linear pressure distribution reads which is identical for the GN and the gB equations.
Since the moment equation used for the GN equations is an approximation, the linearized gB moment equation is used to determine the overturning moment, The integrated pressure and moment, P I and M I , are numerically integrated around the circumference of the cylinder to determine the horizontal force and the overturning moment.

Numerical wave tank setup
Unless otherwise stated, the cylinder diameter is D = 4.0h in this study. This cylinder diameter is large enough to produce significant diffraction, while allows for reasonable numerical computations. Smaller cylinders would require finer grids for the same accuracy and viscous forces may become important. A larger diameter cylinder would require a larger domain. Clearly, the latter two factors would increase the computational time significantly. While the primary wave height is H = 0.2h, various wavelength and wave heights are considered to study the effect of nonlinearity and dispersion on cylinders in shallow waters.
The domain used includes a 20h distance from the upwave boundary to the first cylinder surface, a 20h distance from the last cylinder surface to the downwave boundary and a 20h distance from the far wall to the symmetry axis. It will be shown later that this domain is large enough to avoid problems of wave interactions at the boundaries that affect the resulting forces and moments on the cylinders.
Results are presented in dimensionless form. Unless stated otherwise, the following dimensionless variables are used in this study by selecting (ρ, g, h) as a dimensionally independent set: where the bars represent the dimensionless quantities, ρ is the mass density, F is the horizontal force on the cylinder, M is the overturning moment with respect to the sea floor, P is the pressure, and R is the cylinder radius. Any quantity whose dimension is length is scaled by water depth h. When comparing the results with the available laboratory measurements, a different dimensionless form of the force and moment is utilized; this will be discussed later. The nominal (dimensionless) grid sizes used in this domain are x = 0.25 and y = 0.33. These sizes are small enough to adequately model the surface displacements and large enough to not require excessive CPU (central processing unit) time. To insure stability, the time step must be smaller than the grid size as discussed before. Therefore, the time step is chosen as t = 0.2. Results of a grid convergence study, examining three different grids, namely Grid I, Grid II and Grid III configurations, are shown in Fig. 3. Table 1 shows the dimensionless spatial and time discretization used in these three configurations.

Results and discussion
In this section, results of the GN equations, gB equations and the linearized equations for cnoidal wave interaction with vertical cylinders are presented. First, cnoidal wave diffraction by a single cylinder is considered, and results are compared with other theories and available laboratory measurements. This is followed by cnoidal wave diffraction by multiple, in-line vertical cylinders, where comparison between the different equations is also presented. Namely, results are presented for two and three cylinders lined up in a row. All vertical   cylinders have the same diameters. The spacing between the cylinders, S, is the same and defined in terms of some product of the diameter of the cylinders, S = α D, as the distance between the two closest points on the cylinders. In some other works, see, e.g., Wang and Jiang [70], the spacing is defined as the distance between the axes (center to center) of the cylinders, and therefore, our definition of spacing needs to be modified by adding the cylinder diameter D = R/2 to it when one compares with results by others.

Single vertical cylinder
Several difficulties exist in using the existing experimental data to verify the cnoidal wave forces. Most periodic wave experiments use small-diameter cylinders and short waves that are sinusoidal. This is generally done to obtain the inertia and drag coefficients to be used in a Morison-type analysis of wave forces. This study is primarily interested in large diameter cylinders, long waves and nonlinear cnoidal waves. Because of lack of other data for the conditions considered here, we choose the experimental studies of Chakrabarti and Tam [7] and Chakrabarti [6] for comparison. In these studies, experiments were conducted on a 1.72h diameter cylinder. Although this cylinder is less than one half the diameter used as the standard in our study (4.0h), it is still much larger than the cylinders used in most studies, and it is large enough to provide a qualitative comparison.
In the experiments of Chakrabarti and Tam [7] and Chakrabarti [6], sinusoidal waves were used. Since sinusoidal waves are symmetric about the still-water level, the maximum wave force, F MAX , and minimum wave force, F MIN , produced by sinusoidal waves are nearly equal in magnitude in experiments. Cnoidal waves are not symmetric about the still-water level and do not, in general, produce maximum and minimum wave forces of nearly equal magnitude unless the wave height is very small. Hence, in this study we define the equivalent wave force, F EQ , so that the cnoidal wave forces can be compared to the maximum force amplitude, F AMP , used in Chakrabarti and Tam [7] and Chakrabarti [6]. This procedure removes most of the effect that the non-symmetry of the cnoidal waves has on the resulting forces. The equivalent force is defined here as the average of the magnitudes of the maximum and minimum cnoidal wave forces, i.e., The maximum force amplitude, F AMP , is then taken as the greater of the two values, the magnitude of maximum wave force, F MAX , or the magnitude of the minimum force, F MIN , i.e., F AMP = |F| MAX . Although the equivalent force of the cnoidal wave solutions is compared to the maximum amplitude solution for the sinusoidal wave, the comparison of the forces directly would not have significantly altered the outcome. All the statements made above, in reference to the wave forces, also apply equally well to the wave-induced moments. Five different wavelengths, L = 5.42h, 6.28h, 7.22h, 10.83h, and 15.48h, are used in this comparison of the numerical solutions with the experimental data. The distances from the upwave boundary to the cylinder, from the cylinder to the downwave boundary and from the symmetry boundary to the far wall were all 20h. Solitary waves also are used, and they correspond to kh → 0. See Neill et al. [50] for details on solution of solitary wave impact on vertical cylinders by the GN and gB equations.
The wavelengths were chosen to determine the applicability of the wavelength range for both the gB and GN equations. The 6.28h wavelength corresponds to μ 2 = 1, where μ = kh = 2π , and = h/L, defined in Sect. 2.2. Recall that the gB equations are valid for O(μ 2 ) < 1. Although the GN equations do not have an explicit limit, they must, nevertheless, have similar implicit limitations. Any such limitations of the GN equations must be judged by comparison with experiments. Shields and Webster [64] investigated the applicable range for modeling a two-dimensional wave using the Level I equations used here. They determined that the Level I GN equations are appropriate for wavelengths greater than approximately 7.0h, see also Zhao et al. [84,85]. The wave height of 0.2h is chosen because this wave height was used by Chakrabarti [6]. It is believed to be large enough to produce some nonlinearity.
The present force and moment results numerically determined for this configuration are shown in Fig. 4 for the gB equations and in Fig. 5 for the GN equations for the cylinder diameter of D = 2R = 1.72h. These figures also show the solitary wave limit (L → ∞) of the cnoidal wave diffraction. The forces are compared to the experimental data of Chakrabarti and Tam [7] and to the linear solutions in Fig. 6. These linear solutions are based on the diffraction theory of MacCamy and Fuchs [41], which is for a finite water depth but without the shallow-water approximation.
In Fig. 6, F MAX refers to either F AMP or F EQ . In this figure, only, the force and moment are dimensionless in the form given by Chakrabarti and Tam [7] and Chakrabarti [6], which is different from the normalizations used throughout this study. The normalizations used in Fig. 6 arē where F and M are the dimensional force and overturning moment, respectively. The forces and moments are normalized by the wave height and therefore, one may compare waves of several different wave heights since the experimental forces varied linearly with wave height. However, note that nonlinear forces/moments do not necessarily vary linearly with wave height, or even with its square in some cases; we use this normalization for the purpose of comparison with experiments. Chakrabarti and Tam [7] used waves that varied in heights and μ 2 = (kh) 2 , of experiments of Chakrabarti [6] and Chakrabarti and Tam [7] vs computations of the GN and the Boussinesq equations, and the linear theory solution of MacCamy and Fuchs [41], for the single-cylinder case, D = 1.72h, H = 0.2h. Note that forces and moments are nondimensionalized using Eq. (55) from 0.08h to 0.21h. In a later study, Chakrabarti [6] used waves only with wave heights near 0.20h. To be consistent with these works, the wave forces are plotted versus k R, the wave number multiplied by the cylinder radius. As before, all the statements referring to the forces are also applicable to the moments. Figure 6 shows a very good correlation for the longer wavelength periodic waves between the GN solution and the existing experimental data for both the force and moment. For the shorter wavelengths, the linear solution accurately predicts the force and moment, see Huseby and Grue [34] for further discussion on the agreement of the linear solution with laboratory experiments. The linear solution predicts that the force approaches zero as the wavelength approaches infinity (k R → 0), whereas nonlinear solutions predict the soliton force which is nonzero. (Both the GN and gB solitary wave force and moment limits are shown in Fig. 6.) For the shorter wavelengths, between μ 2 = 0.336 and 1.0, the gB equations appear to not accurately predict either the wave force or moment. An exception is when the gB predictions cross the linear curve and the solutions must agree. For these shorter wavelengths, the gB equations underpredict the dynamic force and moment. This causes the gB equations to predict forces and moments that are less than the linear solutions for these shorter wavelengths. For longer wavelengths, μ 2 < 0.336, the gB equations are in general agreement with the GN equations. Next, we shall study the problem of cnoidal wave interaction with a single cylinder by use of the different equations developed here. Note that hereafter, variables are given in dimensionless form as shown by Eq. (53). Also, the standard numerical setup, discussed in Sect. 4, is used.
For cnoidal waves, it is necessary to verify that the domain size is adequate. After the waves diffract from the cylinder, some of the waves reflect off of the far wall, the upwave boundary and the downwave boundary, and then move back toward the cylinder. This can potentially cause significant errors in the analysis. The larger the domain size is, the longer it takes for the reflected waves to return to the cylinder, and therefore, the less effect these reflections would have on the results. To verify that the standard domain size is sufficiently large, a longer distance (40h) from the wavemaker to the cylinder is also considered in addition to the standard distance (20h). Both the 9.0h and the 11.0h wavelength cases are considered. As shown in Fig. 7, except for the time phase shift (which is expected as it takes longer for the waves to arrive at the cylinder), the forces are nearly identical for the two different domain lengths and for both equations. A similar result was obtained for the longer wavelength case of 11.0h.
For the linear case, the standard 4.0h diameter cylinder and 0.2h wave height is used. Both the 9.0h and the 11.0h wavelength sinusoidal waves are investigated. For the 9.0h wavelength case, the wave heights of 0.1h and 0.3h are also included. For the 9.0h wavelength and 0.2h wave height case, shown in Fig. 8, substantial increase in the forces and moments for both the GN and gB equations is observed. The maximum force (0.42) from the gB case is 14% larger than the linear case (0.37). The maximum force (0.48) for the GN case is 30% larger than the linear case (0.37). The maximum moments for the gB case (0.24) and the GN case (0.25) are 26 and 32% larger, respectively, than the linear case (0.19).
The 11.0h wavelength also showed substantially larger forces and moments for both the GN and gB equations, compared with the linear case. These forces and moments are shown in Fig. 9. The maximum force (0.46) from the gB case is 18% larger than the linear case (0.39). The maximum force (0.54) for the GN case is 38% larger than the linear case (0.39). The maximum moments for the gB case (0.23) and the GN case (0.27) are 21 and 42% larger than the linear case (0.19), respectively.
The nonlinear wave solutions produce larger forces and moments primarily because of the greater run-up of the surface elevation on the side of the cylinder. This is partially caused by the lack of symmetry of the incident surface elevation about the still-water level for the nonlinear waves. Unlike the linear waves, the nonlinear waves have maximum surface elevation greater than one half the wave height. The greater run-up is also caused by the increase in the interaction between the waves and the cylinder. To demonstrate the effect of nonlinearity on the force and the moment, the wave heights of 0.1h, 0.3h, 0.4h, and 0.5h are also included for the 9.0h wavelength case. These forces and moments, shown in Figs. 10, 11, 12, and 13, respectively, demonstrate that as the wave height is increased, the difference in the results from the linear and nonlinear equations increase. This is partially a result of the increase in non-symmetry of the surface elevation. Figure 14 shows the sample surface elevations for the gB and GN equations for the large cylinder case.
The pressure as a function of the vertical coordinate, z, for both the GN equations and the gB equations is shown in Fig. 15 for comparison purposes. This is the total pressure under the wave crest, for a freely       propagating wave, unaffected by diffraction. The wavelength and height used are 9.0h and 0.2h, respectively. The pressure distribution for the gB equations is determined through Eq. (48), and the pressure distribution for the GN equations is determined through Eq. (46). Figure 15 demonstrates that the pressure distribution predicted by the gB and GN equations is nearly linear with depth. This supports the assumption of a linear pressure distribution used for the GN equations. Note, however, that the total pressure is not hydrostatic as expected.
As we have seen, nonlinearity can significantly increase the forces and moments, relative to the linear case. Figure 16 demonstrates this for the case of 9.0h wavelength and D = 4.0h as a function of the dimensionless wave height H. The GN equations predict the highest wave forces among the three models, and unless the wave height is about 0.05h or less, the linear theory should not be used at this wavelength or longer.

Two vertical cylinders
The two-cylinder cnoidal wave case uses the same D = 4.0h diameter cylinder used before. Five different spacings are used: 0.5D, 0.75D, 1.00D, 2.00D, and 3.00D. For the 9.0h-long wave, the spacing of 1.25D = 5.0h is also included to study the case of two in-line cylinders with the cylinder center-to-center distance equal to the wavelength. The cnoidal wavelengths are the same 9.0h and 11.0h as used before, as well as the wave height of 0.2h.   Fig. 17. For the 9.0h wavelength, the gB equations show significant interaction between the cylinders, see Fig. 18. They demonstrate that the forces on the first cylinder could either be increased (maximum 26% increase at S = 1.0D) or decreased (maximum 12% decrease at S = 0.5D) by the presence of the second cylinder, in comparison with the singlecylinder case. The force on the second cylinder is decreased (maximum 28% decrease at S = 0.5D, minimum decrease 3.8% at S = 1.0D) by the presence of the first cylinder.
The GN equations show similar interaction between the cylinders for the 9.0h wavelength case, see Fig.  20. Forces on the first cylinder could either increase (max 28% increase at S = 1.0D) or decrease (max 15% decrease at S = 0.5D). The force on the second cylinder is decreased (maximum 35% decrease at S = 0.5D, minimum decrease 4% at S = 2.0D) by the presence of the first cylinder. For the GN case, the maximum force is predicted at the S = 1.0D spacing. For the gB case the maximum force is predicted at a spacing of S = 1.25D (not shown here but will be discussed later). This difference is caused by a difference in dispersion. The velocity of the waves is affected differently by the wave height for the two cases. Therefore, the incident and the reflected waves combine differently and produce maximum forces at different spacings. In both the gB and GN cases, and for the 9.0h-long wave, the force on the first cylinder either increases or decreases because of the presence of the second cylinder, see Figs. 18 and 20. The force on the second cylinder in general decreases, relative to the force on the first cylinder, by the presence of the first cylinder. In the 11.0h wavelength case, the gB equations show significant interaction between the cylinders, see Fig.  22. The forces on the first cylinder could either increase (max 2% increase at S = 3.0D) or decrease (max 19% decrease at S = 0.75D). The force on the second cylinder could also either increase (max 7% increase at S = 3.0D) or decrease (max decrease 28% at S = 0.5D). This is the only case, among the two-cylinder cases, where the force on the second cylinder is increased by the presence of the first cylinder, relative to the single-cylinder case.
The GN equations also show significant interaction similar to the gB equations for the 11.0h wavelength case, see Fig. 10. The forces on the first cylinder could either be increased (max 14% increase at S = 1.0D) or decreased (max 20% decrease at S = 0.5D) by the presence of the second cylinder. Unlike in the case of the gB equations, however, the force on the second cylinder is only decreased (maximum 24% decrease at S = 0.5D, minimum decrease 4% at S = 3.0D) by the presence of the first cylinder. This difference between the GN and the gB cases is also caused by the greater run-up from the gB case which causes greater wave reflections and therefore stronger interference effects. For the 11.0h-long wave, the force on the first cylinder can either be increased or decreased by the presence of the second cylinder, see Figs. 22 and 23. The two equations do not agree on whether the force on the second cylinder could be increased by the presence of the first cylinder.  The effect of a second cylinder on the resulting moment is very similar to the effect on the forces, see, e.g., Figs. 19 and 21. All the interactions described above also apply to the moment. The values, however, are slightly different.

Three vertical cylinders
For the three-cylinder case, the resulting forces and moments on the first, second, and third cylinders for the 9.0h wavelength gB case are shown in Figs. 25, 26, and 27. A sample surface elevation snapshot for this case is shown in Fig. 24. The resulting forces and moments for the 9.0h wavelength GN case are shown in Figs. 28, 29, and 30.
The effects of the third cylinder on the first cylinder are moderate, see Figs. 25 and 28, and this is expected. The effect of the third cylinder on the first cylinder is to moderately amplify the effects of the two-cylinder case. The lowest force values are slightly lower, and the highest force values are slightly higher on the first cylinder in the three-cylinder case than on the first cylinder in the two-cylinder case. The forces on the first cylinder, in the three-cylinder case, behave similarly to the forces on the first cylinder in the two-cylinder case.
The effect of the third cylinder on the second cylinder is also to amplify the effects of the two-cylinder cases, see Figs. 26 and 29. This amplification is significantly greater for the second cylinder than for the first cylinder. For both equations and both wavelengths, the force on the second cylinder in the three-cylinder case has greater extremes than the force on the second cylinder in the two-cylinder case. The trend of the forces on the third cylinder, in the three-cylinder case, is similar to the second cylinder in the two-cylinder case, see Figs. 27 and 30. The forces, however, are generally less.      The effect of a third cylinder on the resulting moment is very similar to the effect on the forces, see Figs. 25, 26, 27, 28, 29, and 30. Similar to the two-cylinder case, all the interactions described above also apply to the moment. The values, however, are slightly different.

Further discussion on the wave forces
The cnoidal wave forces for both the two-cylinder and three-cylinder cases are shown in Figs. 31, 32, and 33. The interaction from cnoidal waves on multiple cylinders is much more complex than the interactions from the solitary waves. The forces on each cylinder can be increased or decreased by the influence of neighboring cylinders.
For the 9.0h-wavelength case, the forces on the first and second cylinder increases for the spacings near 1.25D (5.0h), see Figs. 31, 32, and 33. This spacing is nearly one half the wavelength (9.0h). The cylinder center-to-center distance for this configuration, however, is equal to the wavelength. Apparently, the reflections from the downwave cylinder meet the upwave side of the upwave cylinder at the same time as the next wave. The waves interfere constructively, and the force on the first cylinder is increased. The wave reflected from the downwave cylinder also reflects off of the downwave side of the upwave cylinder and back toward the downwave cylinder. The reflected wave and the next incoming wave combine and encounter the downwave cylinder at the same instant. This also increases the force on the downwave cylinder. This effect, described above, occurs in the results of both equations and in both the two-and three-cylinder cases. These results are in agreement with the study of Maniar and Newman [43], who found that near-resonant modes can occur between adjacent cylinders at critical wavelengths, resulting in relatively larger wave loads. The 11.0h wavelength case does not show the same constructive interference near the 1.25D spacing as the 9.0h wavelength case, see Figs. 31, 32, and 33. It does, however, show a similar effect for the spacing of 3.0D. This spacing of 3.0D (12.0h) is close in length to the wavelength (11.0h). For this case, the trough of the wave, reflected off of the downwave cylinder, meets the downwave side of the upwave cylinder, at the same instant as the new wave impacts the upwave side. This effect is also found in the results of both equations and for both the two-and three-cylinder cases.
The effect of the third cylinder is to amplify the effects relative to the two-cylinder case. Both the constructive and destructive interferences are increased. For example, as described above for the 9.0h wavelength and 1.0D spacing, the interactions with the waves with the first and second cylinders increase the force on both cylinders as seen in Figs. 31 and 32. A similar effect is found between the second and third cylinders. The result is a significant increase in the force on the second cylinder. These increases are 25% for the gB equations and 14% for the GN equations.
The forces on the third cylinder behave similar to the forces on the second cylinder in the two-cylinder case. The last cylinder in each case exhibited similar behavior. The third cylinder, in general, has more shielding and smaller forces than the second cylinder in the two-cylinder case.
For cnoidal waves, the combined effect of both nonlinearity and multiple cylinder constructive interference can significantly increase the resulting forces on a vertical cylinder. In the 0.2h wave height case, for both the gB and the GN equations, the combined effects are shown to increase the forces up to approximately 45% over the linear single-cylinder case. A larger wave height would produce an even greater difference.
All of the above statements in regard to the forces also apply to the moments.

Summary and conclusions
The problem of interaction of nonlinear waves of cnoidal type with vertical cylinders in shallow water is studied by use of the Green-Naghdi equations, the Boussinesq equations, and linearized form of the two equations. The governing equations are solved in boundary-fitted curvilinear coordinate computational domain by use of a mapping technique. Results are presented in the form of the wave-induced force and moment on the vertical cylinders. Results of the GN, the gB, and the linearized equations are compared with each other in all cases and with the available data for the single-cylinder case. Effect of neighboring cylinders on the wave diffraction is discussed quantitatively. Overall, close agreement is observed between the results of the Green-Naghdi equations and the Boussinesq equations with laboratory measurements and existing theoretical solutions. The agreement is significantly better for the range of long waves. The performance of the Green-Naghdi equations is found to be generally better than the Boussinesq equations. They produce values for the forces and the moments that are in slightly closer agreement with both the experimental data and other predictions. This is true for both the solitary wave and cnoidal wave cases.
It is shown that nonlinearity can significantly increase the forces and moments, relative to the linear case. For the cnoidal waves, the Green-Naghdi equations show a greater increase than the Boussinesq equations.
As the wave height increases, the nonlinearity increases, which causes a corresponding increase in the force and moment. Larger wave heights produce greater difference between the linear and nonlinear equations as expected. Similarly, it is found that the diffraction theory of MacCamy and Fuchs [41] significantly underestimates the loads on the (single) cylinder for long waves.
It is found that the presence of the second and third cylinders on the wave loads on all cylinders is significant in general. In a number of cases studied here, the resultant loads on the first cylinder have increased due to the second and third cylinders, particularly in the case of cnoidal waves. Such effect is found to be a function of the distance between the cylinders. This is in qualitative agreement with the results obtained for wave interaction with an array of vertical cylinders in deep water.
The shielding effect of the leading cylinder on the downwave cylinders varies nonlinearly with the wavelength and the cylinder spacing. In some cases, the loads on the second cylinder are larger than those on the first cylinder. Almost in all cases, however, the force on the third cylinder is smaller than that on the first cylinder. The impact of neighboring cylinders on the cnoidal wave loads is much more complex than solitary wave loads.
The Level I GN equations do not present any relation for the pressure distribution over the water column. It is shown that the (total) pressure varies almost linearly from the free surface of cnoidal waves to the seafloor. The error associated with this assumption is shown to be very small, even for large waves. Further information about the fluid flow properties over the water column can be obtained by use of the high-level GN equations.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http:// creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.