Boussinesq Model and CFD Simulations of Non-Linear Wave Diffraction by a Floating Vertical Cylinder

A mathematical model for the problem of wave diffraction by a floating fixed truncated vertical cylinder is formulated based on Boussinesq equations (BEs). Using Bessel functions in the velocity potentials, the mathematical problem is solved for second-order wave amplitudes by applying a perturbation technique and matching conditions. On the other hand, computational fluid dynamics (CFD) simulation results of normalized free surface elevations and wave heights are compared against experimental fluid data (EFD) and numerical data available in the literature. In order to check the fidelity and accuracy of the Boussinesq model (BM), the results of the second-order super-harmonic wave amplitude around the vertical cylinder are compared with CFD results. The comparison shows a good level of agreement between Boussinesq, CFD, EFD, and numerical data. In addition, wave forces and moments acting on the cylinder and the pressure distribution around the vertical cylinder are analyzed from CFD simulations. Based on analytical solutions, the effects of radius, wave number, water depth, and depth parameters at specific elevations on the second-order sub-harmonic wave amplitudes are analyzed.


Introduction
Over the years, there has been considerable progress in the development of various design techniques for the study of floating structures for various purposes in different water depths, such as floating platforms, wave energy devices, and breakwaters. Due to the world's ever-increasing energy requirement, continuous depletion of fossil fuels and increasing environmental awareness, the energy industry has shifted focus towards renewable energy sources. Apart from solar and wind resources, waves have attracted high interest lately because of the promising prospect of wave energy converters. In recent years, several research groups have taken a keen interest in the development of efficient wave energy converters. In most cases, wave load analysis on the floating structures is performed using various linearized analytical models. For a better understanding the behavior of the wave energy devices under the action of waves, for the assessment of their motions, loads, and for power production, the wave energy industry puts trust in the utilization of these models. However, many physical problems associated with wave-structure interactions are essentially non-linear in nature and require non-linear models to properly realize the impact of wave load on structures.
The generic wave energy converter (WEC) concepts include point absorbers, attenuators, oscillating wave surge converters, oscillating water columns, overtopping devices and submerged pressure differentials. When it comes to the study of WECs, a large part of the topic consists of modeling the wave-structure interaction. Over recent decades, along with experimental studies, several linear, locations is performed between BM and CFD simulation results to ascertain their accuracy and to analyze the fidelity level of the developed BM. In addition, the effects of the radius of the cylinder, wave number, water depth, and depth parameters on sub-harmonic second-order wave amplitude at a specific location around the cylinder is analyzed. Finally, significant conclusions are drawn from the comparative study and future scope of the present model.

Mathematical Model
The mathematical formulation of the problem is considered in three-dimensional Cartesian coordinate system (x*, y*, z*), x*-y* being the horizontal plane and z*-being pointing upward starting at the undisturbed free surface. The fluid occupies the region -h* < z* < 0, -∞ < x*< ∞ and a cylinder of radius r = a is floating at z* = 0 with draft d and water depth at z = −h*. * indicates a dimensional variable. Hence, the fluid region is divided into two regions: an interior which corresponds to the area below the cylinder and exterior region constitute the rest of the region (Figure 1). Therefore, the velocity potential and the free surface displacement in the exterior region are denoted by Φ* and η*, respectively. Furthermore, the fluid is assumed to be inviscid, incompressible, and the motion is irrotational, and simple harmonic in time with angular frequency ω.
where P is the pressure, g = gravitational constant, and ρ = fluid density. From Equation (3), the pressure continuity at the free surface ** z   by considering atmospheric pressure zero leads to The velocity potential Φ* satisfy the Laplace equation as: where ∇ ≡ (∂/∂x * , ∂/∂y * ) in Cartesian co-ordinate system, the subscript 'z * z * ' denotes the double partial derivatives w.r.t to z * . As the bottom is considered as rigid, the bottom boundary condition yields.
The equations of motion for the exterior region can be expressed as:

of 27
where P is the pressure, g = gravitational constant, and ρ = fluid density. From Equation (3), the pressure continuity at the free surface z * = η * by considering atmospheric pressure zero leads to ∂Φ * ∂t * + 1 2 (∇Φ * ) 2 + gη * = 0 (4) Now, the problem is governed by one-linear equations and three boundary conditions which are non-linear in the potential variables Φ * and η * . Therefore, the problem itself is non-linear as η * is not known a priori.
The next subsection derives the nonlinear BEs for exterior region associated with free surface motions in the presence of floating fixed vertical cylinder.

Non-Linear Boussinesq Equations (BEs) in the Exterior Region
The following non-dimensional variables are used for the derivation of nonlinear BEs with a floating truncated vertical cylinder (as in [25,26]): where h= length representative of the water depth scaled by h 0 , t * = time, p * = water pressure, λ is the wave length, I 0 is the wave amplitude, η * = free surface elevation, and r * = radius of the cylinder. In order to characterizing long waves, the following two non-dimensional parameters ε and µ are used as: where ε and µ in Equation (6) are the nonlinearity and frequency dispersion and are assumed to be small. The non-dimensional three-dimensional governing equation in terms of velocity potential Φ is expressed as (hereafter all the non-dimensional symbols are neglected): where ∇ 2 = ∂ 2 ∂r 2 + 1 r ∂ ∂r + 1 r 2 ∂ ∂θ 2 , with r is the radial distance from the z-axis and θ is the angle about the x-axis.
The generalized Bernoulli equation read as dynamic boundary condition: where the subscripts 't' and z denote the partial derivative w. r. t. t and z, respectively. U r = ε/µ 2 is known as the Ursell number. The kinematic condition at z = εη in non-dimensional form is given by: No flow condition at the bottom boundary yields: Integrating Equation (7) w. r. t. z and using boundary conditions (9-10), results in: The kinematic boundary condition on the cylinder is given by: The velocity potential satisfies the far field radiation condition as The next section presents the non-linear BEs in terms of velocity potential with depth parameters at specific elevations in the exterior region based on perturbation technique and second-order wave amplitudes around the floating fixed truncated vertical cylinder in the application range of weakly dispersive waves.

BEs for Exterior Region Based on Perturbation Technique
The velocity potential Φ in polynomial form can be expanded as: where x = r cos θ and y = r sin θ. Substituting Equation (5) into Equations (7) and (10) equating the powers of µ in terms of the specific elevation at z = z α (−h ≤ z α ≤ 0) and related velocity potential Φ α , one can obtain (as in [8,27]).
Again, substituting Equation (15) into Equation (11), and on simplification, the BE in terms of Φ α and η in the exterior region is obtained as: Substituting Φ from Equation (15) into Equation (8), the nonlinear free surface condition is derived as: where the higher-order non-linear terms are neglected.
The next subsection presents the second-order free surface elevations and velocity potentials in the exterior and interior regions using suitable Bessel's and Hankel's functions associated with wave diffraction by a floating fixed truncated vertical cylinder.

Second-Order Free Surface Wave and Velocity Potentials for Exterior and Interior Region
The second-order free surface displacement in exterior region r ≥ a can be expressed as: where A ± mn are the unknown second-order wave amplitude to be determined. Further, the + sign and -sign denote the super-and sub-harmonic parts generated by the interactions of the linear waves (m-wave and n-wave), respectively, J l (k ± mn r) is the Bessel function of first kind of order l (as in [3]), k ± mn = k m ± k n , ω ± mn = ω m ± ω n , and ϕ ± mn = ϕ m ± ϕ n with k m and k n are the wave numbers of the m and n-waves associated with the wave frequency ω m and ω n , respectively. Further, ϕ m and ϕ n are the phases of m-and n-wave, respectively.
The incident velocity potential Φ inc αmn and the diffracted velocity potential Φ di f αmn can be expanded as: where H l (k ± mn r) is the Hankel function of first kind of order l. B ± mn and C ± mn are the unknown coefficients of super-and sub-harmonic waves.
Hence, using Equations (19) and (20), the total potential Φ ext αmn can be expressed as: The second-order velocity potential Φ int αmn in the interior region can be written as: where D ± mn is the unknown coefficient and I l (k ± mn r) are the modified Bessel functions of first kind of order l.
The next subsection determines the unknown coefficients in terms of second-order superand sub-harmonic amplitudes associated with free surface elevations and velocity potentials using matching conditions.

Determination of Second-Order Free Surface Amplitude and Unknown Coefficients
The unknowns A ± mn , B ± mn , C ± mn and D ± mn associated with Equations (18)- (20) and Equation (22) can be determined by applying the structural boundary conditions and matching the pressure and normal velocities across the boundary at r = a. The boundary conditions are given by: ∂φ ext The wave amplitudes A ± mn defined in Equation (18) can be determined by applying the structural boundary condition (23) and using Equation (18) into Equation (15). However, hereafter the super-harmonic wave amplitude is considered for comparing with CFD model simulations and the effect of sub-harmonic second-order wave amplitudes on different parametric values in subsequent sections.
Furthermore, the velocity continuity condition (25) gives Further, applying the structural boundary condition (23) and using Equation (18) into Equation (17), one can obtain the second-order wave amplitude A ± mn as Now solving Equations (26)- (28), the second-order wave amplitude D ± mn is obtained as where where H l , H l and H l denote to the 1st, 2nd, and 3rd derivatives of Hankel's function, respectively, and the same as with Bessel's function J l . The next section describes the numerical CFD model based on OpenFOAM solver in terms of CFD, mathematical model equations, computational mesh, the case set up, and computational resource associated with the aforementioned problem.

Brief Description of CFD Model
For CFD simulations in the study, the open source CFD tooklit OpenFOAM (Open Source Field Operation and Manipulation) was used. A relatively old version of OpenFOAM ws used (version 2.4.0), coupled with waves2Foam utility for wave generation. A detailed description on the CFD solver can be found at works by Jasak [28].
OpenFOAM is governed by the Navier-Stokes and continuity equation, and models an incompressible laminar flow of a Newtonian fluid. The Navier-Stokes equation is given by: where U is the velocity, p is the pressure, µ is the dynamic viscosity, g is acceleration due to gravity, and ∇ 2 is the Laplace operator. Furthermore, the continuity equation is of the form: The governing equations are discretized using the Finite Volume Method (FVM). Further special and temporal discretization are performed using second and first order methods. Pressure velocity coupling is done using the PISO-SIMPLE (PIMPLE) algorithm. The solver follows the Cartesian coordinate system. The volume of fluid (VoF) method together with the artificial compression technique tracks the free surface elevation, using the volume fraction technique. The equation for the volume fraction is obtained as where U is the velocity field, α is the volume fraction of water in the cell and varies from 0 to 1, representing full of air to full of water, respectively. OpenFOAM incorporates several two-equation turbulence models for turbulence modelling and also has incorporated a wall function model. For this study, kOmegaSST turbulence model was used, and the waves were generated following Stokes second-order wave theory.

Computational Mesh
For mesh generation for the case, the OpenFOAM mesh generation utility was used. Depending on the settings, the utility generates hybrid mesh that incorporates both structures and unstructured mesh for an efficient solution. Initially, blockMesh utility is used to generate the rectangular simulation domain using hexahedral mesh. In the mesh, y ' + ve' represents vertical positive direction, and the horizontal free surface is in the x-z plane. The domain was further refined using toposet and refinmesh commands in specific locations, near the geometry and free surface, for better flow prediction around the geometry and better propagation of waves. Finally, the cylinder geometry was integrated to the domain using the snappyHexMesh utility, no refinement or mesh layer was applied in snappyHexMesh.
Two different mesh resolutions were used for simulating the cases to ensure that the mesh satisfies the minimum number of cell requirements per wave height and wave length, while consuming minimum resources. The domain lengths and depths were also varied for the various cases to ensure enough space for the wave relaxation zone and propagation zone, and to maintain a relatively shallow water case for comparing the results with the Boussinesq model. Mesh resolution was maintained such that for the smallest wave length and wave height case, there were 77 cells per wave length and 8 cells per wave height. However, this was the configuration for the shortest wave length and wave height, and the later waves had more cells per wave length and height. For the first five cases, a cell spacing of 0.01 in the X and Z direction was used, and for Y it was 0.005. For the rest of the cases, the cell spacing was 0.05 for the X and Z direction and 0.025 for Y. An average mesh resolution of 6 million cells (mostly hexahedral) was used for the cases. The mesh setup for the first cases is shown in Figure 2. As for the later cases (with longer wave lengths), the refinement was confined to the near cylinder area, where wave probes have been placed, instead of the entire wave propagation area. Since the meshes were generated following standard guidelines for wave simulation, an independent grid dependency study was not performed. Nevertheless, the authors' previous experience from wave simulations indicate that the mesh resolution and distribution used is sufficient to produce reliable results.

Case Setup
The CFD simulations were performed following the case experimental setup by Mercier and Niedzwecki [25]. Few changes were added to the setup, and all simulations depths were adjusted to maintain a water depth to wavelength ratio (h/L) of 0.5 (water depth/wavelength). This was done to ensure comparability with weakly dispersive Boussinesq model results. Furthermore, considering the relatively small diameter of the cylinder (0.23 m), the domain in lateral direction was reduced from 2.8 m to 1.5 m. The domain length was also adjusted depending on applied wave lengths to ensure enough propagation and dissipation area (at the end of domain) for the incoming waves. On average, the domain length was 12 m.
Following the experimental setup, five wave elevation probes were placed in the simulation domain, three near the cylinder wall and two at a distance, to capture the free surface elevation. The wave probes are denoted by A, B, C, D and E. The locations of the probes are provided in Table 1. The locations of the probes/gauges near the cylinder were slightly modified to place them at the first cell after the cylinder in order to avoid numerical exception. The displacements were at the fraction of millimeters, so it should not impact the simulation outcome. A general arrangement of the experimental setup and simulation domain are shown in Figures 3 and 4, respectively.  Following the experimental setup, five wave elevation probes were placed in the simulation domain, three near the cylinder wall and two at a distance, to capture the free surface elevation. The wave probes are denoted by A, B, C, D and E. The locations of the probes are provided in Table 1. The locations of the probes/gauges near the cylinder were slightly modified to place them at the first cell after the cylinder in order to avoid numerical exception. The displacements were at the fraction of millimeters, so it should not impact the simulation outcome. A general arrangement of the experimental setup and simulation domain are shown in Figures 3 and 4, respectively.  All CFD simulations were performed in a desktop machine with Intel ® core i7-7700 CPU with 3.60 GHz processor and 32 GB of RAM. The data was written on an solid state drive (SSD)disk. On average, each case was run for roughly 100 s (simulation time) to ensure around 10 wave encounters and the physical time took to run each case was roughly 75 h.
In the context of the present paper, the 3D numerical simulations were performed for the comparison of results with other methods like experimental data, immersed boundary method (IBM), numerical wave tank (NS) and the present developed BM. The comparison should provide a relative idea about the efficiency and limitations of the methods and models. Therefore, the results obtained from the present BMs and CFD simulations and different method results available in the literature are presented in the subsequent subsections.  All CFD simulations were performed in a desktop machine with Intel ® core i7-7700 CPU with 3.60 GHz processor and 32 GB of RAM. The data was written on an solid state drive (SSD)disk. On average, each case was run for roughly 100 s (simulation time) to ensure around 10 wave encounters and the physical time took to run each case was roughly 75 h.
In the context of the present paper, the 3D numerical simulations were performed for the comparison of results with other methods like experimental data, immersed boundary method (IBM), numerical wave tank (NS) and the present developed BM. The comparison should provide a relative idea about the efficiency and limitations of the methods and models. Therefore, the results obtained from the present BMs and CFD simulations and different method results available in the literature are presented in the subsequent subsections.

Numerical Results
In this section, several numerical results associated with the comparison among CFD, experimental, numerical data, and analytical BM results, wave forces, moment, pressure distribution   All CFD simulations were performed in a desktop machine with Intel ® core i7-7700 CPU with 3.60 GHz processor and 32 GB of RAM. The data was written on an solid state drive (SSD)disk. On average, each case was run for roughly 100 s (simulation time) to ensure around 10 wave encounters and the physical time took to run each case was roughly 75 h.
In the context of the present paper, the 3D numerical simulations were performed for the comparison of results with other methods like experimental data, immersed boundary method (IBM), numerical wave tank (NS) and the present developed BM. The comparison should provide a relative idea about the efficiency and limitations of the methods and models. Therefore, the results obtained from the present BMs and CFD simulations and different method results available in the literature are presented in the subsequent subsections.

Numerical Results
In this section, several numerical results associated with the comparison among CFD, experimental, numerical data, and analytical BM results, wave forces, moment, pressure distribution

Numerical Results
In this section, several numerical results associated with the comparison among CFD, experimental, numerical data, and analytical BM results, wave forces, moment, pressure distribution around the cylinder, effects of radius, wave number, water depth, and depth parameters at specific elevations are analyzed. First, comparisons are made of normalized free surface elevations/wave heights among CFD, EFD, and numerical data available in the literature in the probe locations mentioned in the subsequent subsection. Then, the simulation results of wave forces and moments acting on the cylinder and pressure distribution around the vertical cylinder are presented. To check the accuracy of the second-order super-harmonic wave amplitude Boussinesq results, comparisons between the Boussinesq and CFD model simulations on the second-order super-harmonic wave amplitudes for the same probe locations are presented. Finally, to analyze the effects of radius, wave number, water depth, and depth parameters at specific elevations on the second-order sub-harmonic wave amplitudes have been studied.

Comparison and CFD Simulations
Validation studies for CFD simulation results were undertaken for cases performed in [5,29]. Nine cases were simulated, for different wave lengths and wave heights. The test cases are shown in Table 2. For validation, the simulation results were compared with EFD reported in [28], IBM data reported [5], and numerical wave tank (NWT) data generated using Navier-Stokes (NS) and marker and cell (MAC) methods reported in [30]. For normalization of the free surface elevation data from CFD, the total free surface displacement was divided by the incoming wave height, H/H 0 . Calculations were also made based on amplitude since the analytical Boussinesq method calculated the free surface displacement based on amplitude, A/A 0 . The target of the simulations was to analyze the diffraction of the incoming waves due to interaction with a floating fixed vertical truncated cylinder. Thus, three wave probes were placed at locations close to the cylinder body, one at 0.48 m before the cylinder and one at 0.48 m at the side ( Figure 3). Figure 5 shows the normalized free surface elevation just before the cylinder at location A. The results show that almost for all cases, the free surface shows higher elevation than the incoming incident waves. Especially for the low wave period cases, the waves surge up to 2.4 times of initial wave height on the cylinder wall. The CFD results here show reasonable agreement with experimental data, although the values do not exactly agree, the trend is well captured. For almost all cases, CFD under predicts the results compared to EFD results. Slight deviation might be explained to some extent by the location of the probe in the CFD model. In CFD, the probe was set 5 mm away from the cylinder, whereas in experiments it was exactly on the cylinder. This could not be done in the CFD model since probes in simulation can only capture data from a cell, not at voids. Furthermore, the numerical dissipation also contributes to the under prediction. The presented CFD results also agree better with the NS model results, since both follow a similar mathematical model, compared to the immersed boundary method results. Figure 6 shows surface elevation at the side of the cylinder, at location B. The results show that only for short wave length cases, the incident wave becomes distorted at the side; however, for relatively longer wave lengths, the incident wave height remains almost unchanged. The CFD results here capture the resulting trend, however, they show oscillatory predictions around the EFD data, that is, some cases are over predicted and others are under predicted. Nevertheless, the overall deviation is minimum. Also, for these cases, CFD shows better agreement than the NS model.
here capture the resulting trend, however, they show oscillatory predictions around the EFD data, that is, some cases are over predicted and others are under predicted. Nevertheless, the overall deviation is minimum. Also, for these cases, CFD shows better agreement than the NS model.   here capture the resulting trend, however, they show oscillatory predictions around the EFD data, that is, some cases are over predicted and others are under predicted. Nevertheless, the overall deviation is minimum. Also, for these cases, CFD shows better agreement than the NS model.    Figure 7 gives free surface elevation just behind the cylinder (location C). It is observed that the diffraction effect is minimum behind the cylinder and the simulations results here capture the trend well, with relatively notable deviation for three cases. For short wavelengths, the waves lose amplitude while interacting with the cylinder and thus show a lower steepness behind the cylinder. As the wave period and height increases, the loss in steepness decreases and the waves end up gaining slight amplitude due to the interaction with the cylinder. However, none of the amplitude gains actually represent a gain in energy, rather just a sudden transfer of energy.  Figure 7 gives free surface elevation just behind the cylinder (location C). It is observed that the diffraction effect is minimum behind the cylinder and the simulations results here capture the trend well, with relatively notable deviation for three cases. For short wavelengths, the waves lose amplitude while interacting with the cylinder and thus show a lower steepness behind the cylinder. As the wave period and height increases, the loss in steepness decreases and the waves end up gaining slight amplitude due to the interaction with the cylinder. However, none of the amplitude gains actually represent a gain in energy, rather just a sudden transfer of energy. In Figures 8 and 9, the comparison between CFD, experimental, and numerical data in location D and E are shown, which represent free surface elevation at 0.48 m before and at the side of the cylinder. The results show that the effect of the wave and structure interaction can also be observed at roughly one radial distance from the cylinder. As before, the CFD results capture the trend well, with deviation in the values for several cases. In Figures 8 and 9, the comparison between CFD, experimental, and numerical data in location D and E are shown, which represent free surface elevation at 0.48 m before and at the side of the cylinder. The results show that the effect of the wave and structure interaction can also be observed at roughly one radial distance from the cylinder. As before, the CFD results capture the trend well, with deviation in the values for several cases.  Figure 7 gives free surface elevation just behind the cylinder (location C). It is observed that the diffraction effect is minimum behind the cylinder and the simulations results here capture the trend well, with relatively notable deviation for three cases. For short wavelengths, the waves lose amplitude while interacting with the cylinder and thus show a lower steepness behind the cylinder. As the wave period and height increases, the loss in steepness decreases and the waves end up gaining slight amplitude due to the interaction with the cylinder. However, none of the amplitude gains actually represent a gain in energy, rather just a sudden transfer of energy. In Figures 8 and 9, the comparison between CFD, experimental, and numerical data in location D and E are shown, which represent free surface elevation at 0.48 m before and at the side of the cylinder. The results show that the effect of the wave and structure interaction can also be observed at roughly one radial distance from the cylinder. As before, the CFD results capture the trend well, with deviation in the values for several cases. Especially, in case of Figure 9 (position E), the CFD values show relatively high deviation from experimental values, except for three cases. The reason for this over prediction may partly be explained by the relatively short width of the domain in the sides. Because of the slight reflection of the waves by the symmetry planes, the elevations might have been over predicted. Increasing the domain size in the lateral direction might have improved the prediction, however, considering the added computational expense, it was avoided. Next, horizontal and vertical forces and pitch moment acting on the cylinder were computed from the simulation results. For non-dimensionalization of the results, the following equations were used. ' ' where Fx, Fy and Mz represent horizontal and vertical force and pitch moment, respectively,  is the fluid density, g is gravitational acceleration, is the cylinder radius, T is the cylinder draft and H is the incident wave height. In Figure 10, the horizontal force, vertical force, and moment acting on the cylinder versus wave period are plotted. It is observed that both the horizontal force and pitch moment show similar response with peak amplitude being observed close to the 1.70 wave period. The force and moment decrease for relatively low and high frequency range. Whereas, for the heaving force, it keeps on increasing with the increasing wave period and wave height. This might be because the structure is fixed and is not free to heave with incoming incident waves. With heave-free motion, the results might be slightly different. Especially, in case of Figure 9 (position E), the CFD values show relatively high deviation from experimental values, except for three cases. The reason for this over prediction may partly be explained by the relatively short width of the domain in the sides. Because of the slight reflection of the waves by the symmetry planes, the elevations might have been over predicted. Increasing the domain size in the lateral direction might have improved the prediction, however, considering the added computational expense, it was avoided.
Next, horizontal and vertical forces and pitch moment acting on the cylinder were computed from the simulation results. For non-dimensionalization of the results, the following equations were used.
where F x , F y and M z represent horizontal and vertical force and pitch moment, respectively, ρ is the fluid density, g is gravitational acceleration, r is the cylinder radius, T is the cylinder draft and H is the incident wave height. In Figure 10, the horizontal force, vertical force, and moment acting on the cylinder versus wave period are plotted. It is observed that both the horizontal force and pitch moment show similar response with peak amplitude being observed close to the 1.70 wave period. The force and moment decrease for relatively low and high frequency range. Whereas, for the heaving force, it keeps on increasing with the increasing wave period and wave height. This might be because the structure is fixed and is not free to heave with incoming incident waves. With heave-free motion, the results might be slightly different.  To further illustrate the CFD results, the pressure distribution in the simulation domain is shown in Figures 11 and 12. The figures show the change of pressure in the domain due to the interaction between the cylinder and the incoming waves. They also show the disturbance in the incoming incident waves due to interaction with the diffracted waves by the cylinder. Figure 11 shows the interaction of the wave crest with the vertical cylinder. The cylinder encounters high horizontal pressure during the interaction, as can be seen from the deep red color near the free surface. The horizontal pressure gradually fades away with increasing water depth. The contours on the cylinder also show the vertical (upward) force encountered by the cylinder during the wave interaction. Figure 12 shows both the pressure distribution on the free surface during the wave encounter by the cylinder and also reference axis for observing the free-surface elevation. The figure shows how the free surface deforms around the cylinder during wave propagation. The cylinder faces low pressure at the wave trough and the free surface is close to −0.23. Next, the wave passes and when the cylinder encounters the wave peak, there is a jump of the free surface at the front of the cylinder and the wave reaches a height of roughly 0.3. During the wave encounter, the elevation of the free surface at the side and behind the cylinder is lower than that at the front.

Comparison of Analytical Boussinesq Model (BM) Results against CFD Simulations
A MATLAB code is developed to compute the Bessel and Hankel functions in the analytical solutions using built-in sub-routine Matlab to present numerical results of Boussinesq solutions for locations A, B, C, D, and E (as discussed in Section 3 and labeled in Figure 3). Figures 13-17 show the comparisons of results of Boussinesq second-order super-harmonic wave amplitude A + mn against CFD model simulations in the exterior region for locations A, B, D, C, and E, respectively, versus time period T(s). From Figures 13 and 14, it is observed that the peak amplitude of the Boussinesq result agreed well and the wave amplitude trend is similar with CFD model results. However, the wave amplitude in location A is higher than that of location B.
In Figure 15, the results of second-order super-harmonic wave amplitude A + mn are in good agreement between the Boussinesq and CFD model simulations for a smaller wave period in location C (behind the cylinder). However, for higher wave period the results between the two models did not agree well although the trend is similar. This may be due to the numerical errors arising from time and spatial discretization that are more evident with the disagreement of the wave amplitudes in between two models. Furthermore, unlike the CFD simulations, the BM is based on the frequency domain and the numerical model requires some time for the forcing wave to be applied and to stabilize behind the cylinder (location C). This may be the reason to lead the discrepancy in comparison of the Boussinesq results and those of CFD simulations. To further illustrate the CFD results, the pressure distribution in the simulation domain is shown in Figures 11 and 12. The figures show the change of pressure in the domain due to the interaction between the cylinder and the incoming waves. They also show the disturbance in the incoming incident waves due to interaction with the diffracted waves by the cylinder. Figure 11 shows the interaction of the wave crest with the vertical cylinder. The cylinder encounters high horizontal pressure during the interaction, as can be seen from the deep red color near the free surface. The horizontal pressure gradually fades away with increasing water depth. The contours on the cylinder also show the vertical (upward) force encountered by the cylinder during the wave interaction.   Figure 12 shows both the pressure distribution on the free surface during the wave encounter by the cylinder and also reference axis for observing the free-surface elevation. The figure shows how the free surface deforms around the cylinder during wave propagation. The cylinder faces low pressure at the wave trough and the free surface is close to −0.23. Next, the wave passes and when the cylinder encounters the wave peak, there is a jump of the free surface at the front of the cylinder and the wave reaches a height of roughly 0.3. During the wave encounter, the elevation of the free surface at the side and behind the cylinder is lower than that at the front.

Comparison of Analytical Boussinesq Model (BM) Results against CFD Simulations
A MATLAB code is developed to compute the Bessel and Hankel functions in the analytical solutions using built-in sub-routine Matlab to present numerical results of Boussinesq solutions for locations A, B, C, D, and E (as discussed in Section 3 and labeled in Figure 3).  From Figures 13-14, it is observed that the peak amplitude of the Boussinesq result agreed well and the wave amplitude trend is similar with CFD model results. However, the wave amplitude in location A is higher than that of location B.    In Figure 15, the results of second-order super-harmonic wave amplitude mn A  are in good agreement between the Boussinesq and CFD model simulations for a smaller wave period in location C (behind the cylinder). However, for higher wave period the results between the two models did not agree well although the trend is similar. This may be due to the numerical errors arising from time and spatial discretization that are more evident with the disagreement of the wave amplitudes in between two models. Furthermore, unlike the CFD simulations, the BM is based on the frequency domain and the numerical model requires some time for the forcing wave to be applied and to   Comparison of Boussinesq, second-order super-harmonic wave amplitude A + mn in Equation (29) and CFD models of free surface elevation for location C. Figure 16 compares the results of second-order super-harmonic wave amplitude A + mn between Boussinesq and CFD model for location D. It is observed that the peak amplitude and trend between the Boussinesq and CFD models are fully agreed. It suggests that the interaction of second-order wave amplitude and vertical cylinder exhibits similar behavior in their trends for both models. In Figure 17, it is found that the observations are similar to those of Figures 13 and 14. However, for the higher wave period, the CFD data points are a little far and they only agree in the trend, not their values. The reason is similar to those of Figure 15.
The next subsection analyses the effect of second-order sub-harmonic waves on the different parameters associated with wave diffraction by a floating fixed truncated vertical cylinder.    Boussinesq and CFD model for location D. It is observed that the peak amplitude and trend between the Boussinesq and CFD models are fully agreed. It suggests that the interaction of second-order wave amplitude and vertical cylinder exhibits similar behavior in their trends for both models. In Figure 17, it is found that the observations are similar to those of Figures 13-14. However, for the higher wave period, the CFD data points are a little far and they only agree in the trend, not their values. The reason is similar to those of Figure 15. The next subsection analyses the effect of second-order sub-harmonic waves on the different parameters associated with wave diffraction by a floating fixed truncated vertical cylinder.

Effect of Design Parameters on Second-Order Sub-Harmonic Wave Amplitudes
In Figure 18, the effects of radius r on second-order sub-harmonic wave amplitude. A  at Figure 17. Comparison of Boussinesq, second-order super-harmonic wave amplitude A + mn in Equation (29) and CFD models of free surface elevation for location E.

Effect of Design Parameters on Second-Order Sub-Harmonic Wave Amplitudes
In Figure 18, the effects of radius r on second-order sub-harmonic wave amplitude. A − mn at location B for h = 4.0 m, k m = 0.84 and z α = 0.5 versus wave period T(s) is plotted. It is observed that the second-order sub-harmonic wave amplitude A − mn increases with an increase in the radius of the cylinder. The observations are similar to [31] in the case of wave diffraction by a cylinder in two-layer fluids. This may be due to fact that as the cylinder occupied more water plane area, the large amount of incident wave energy concentrates around the structure and on the other hand, physically, this result may be from an increase in fluid momentum as the flow negotiates the cylinder.   Overall, it is observed that the sub-harmonic second-order wave amplitude decreases with the increase in water depths for higher values of the wave period. This is due to an increase in hydrostatic pressure, the force per unit is exerted by water on the floating cylinder. In Figure 21, the effects of depth parameters on secondorder sub-harmonic wave amplitude   Figure 19 shows the effects of wavenumber k m on second-order sub-harmonic wave amplitude A − mn at location B for h = 4.0 m, r = 0.35 m and z α = 0.5 versus wave period T(s). It is observed that the second-order sub-harmonic wave amplitude A − mn decreases with increase in wave number values. Moreover, the observation is similar to [8]. Furthermore, it is seen that for higher wave period the variations of second-order sub-harmonic wave amplitude becomes negligible for fixed values of radius and for different wave number.   Overall, it is observed that the sub-harmonic second-order wave amplitude decreases with the increase in water depths for higher values of the wave period. This is due to an increase in hydrostatic pressure, the force per unit is exerted by water on the floating cylinder. In Figure 21,   Overall, it is observed that the sub-harmonic second-order wave amplitude decreases with the increase in water depths for higher values of the wave period. This is due to an increase in hydrostatic pressure, the force per unit is exerted by water on the floating cylinder. In Figure 21, the effects of depth parameters on second-order sub-harmonic wave amplitude A − mn at location B for r = 0.35 m, k m = 0.84 and h = 4.0 m, and versus wave period T(s) is plotted. It is observed that the second-order sub-harmonic wave amplitude A − mn increases with increase in depth parameters whilst the effect is negligible for a smaller wave period.   All numerical computations associated with analytical expressions were performed in a desktop machine with Intel ® core i7-4790 CPU with 3.60GHz processor and 32 GB of ram memory. On an   All numerical computations associated with analytical expressions were performed in a desktop machine with Intel ® core i7-4790 CPU with 3.60GHz processor and 32 GB of ram memory. On an average, each case took roughly 10-15 min to finish. All numerical computations associated with analytical expressions were performed in a desktop machine with Intel ® core i7-4790 CPU with 3.60GHz processor and 32 GB of ram memory. On an average, each case took roughly 10-15 min to finish.

Discussion
In the theoretical perspective, the contribution of the present boundary value problem is the formulation associated with wave diffraction by a floating fixed truncated vertical cylinder based on Boussinesq-type equations using Bessel's and Hankel's functions, which lead to further formulation of the second-order free surface elevation displacement and velocity potentials in the exterior domain. Therefore, importantly, the comparisons of super-harmonic wave amplitude for different probe locations between BM, numerical CFD simulations, experimental, and other numerical data, and some other effects are also discussed in the succeeding paragraphs.
Initially, in order to show efficiency and level of exactness of the CFD model simulations relative to other methods, the results of wave amplitudes and wave heights from CFD are compared against EFD and numerical data available in the literature for the different probe locations. It is found that the CFD results here show reasonable agreement with experimental data; although the values do not exactly agree, the trend is well captured. For almost all cases, CFD under predicts the results compared to EFD results. The slight deviation might be explained to some extent by the location of the probe in CFD model. In CFD, the probe was set 5 mm away from the cylinder, whereas, in experiments, it was exactly on the cylinder. This could not be done in the CFD model since probes in simulation can only capture data from a cell, not at voids. The numerical dissipation of incoming waves and instability of the turbulence model at the free surface also contributes to the under prediction. In general, the CFD results agree better with the NS model results, compared to the immersed boundary method results, since the present CFD model also follows the RANS model. Overall, the CFD results show reasonable agreement with the experimental and other numerical data, as presented in Figures 5-9.
Then, the numerical results of second-order super-harmonic wave amplitudes obtained from BM are compared with CFD results in the same probe locations. It is observed that the peak amplitudes are in good agreement between the two models and the trend of the wave profiles are similar in nature for all cases, as shown in Figures 13-17. However, the CFD data points are little far from the BM results. These differences can be justified by the full non-linearity solution and viscous effects of the CFD model, while the solution of the analytical model is only of second-order. However, the comparison showed that the developed BM is supported by the numerical CFD model which indicates that the BM may be reliable for concept validation and further enhancement of the model to study motion characteristics of the cylinder to design a floating point absorber wave energy converter.
Regarding the force and pressure distribution analysis, it may be mentioned that the value of the horizontal force is in between the vertical force and moment acting on the vertical cylinder. The peak amplitude of moment is higher than those of horizontal and vertical force. Regarding the pressure distribution around the cylinder, it is found that during the wave encounter, the free surface elevation at the side and behind the cylinder is lower than that at the front probe location.
Finally, the effects of radius, wave number, water depth, and depth parameters on the second-order sub-harmonic wave amplitudes from the Boussinesq results are analyzed, as presented in Figures 18-21. It is observed that the second-order subharmonic waves around the cylinder are negative for all cases whilst, in the same range the super-harmonic waves are positive as presented in Figures 13-17. Furthermore, it is seen that the effects of the second-order on super-harmonic waves for free surface wave profiles around the cylinder are stronger than those of sub-harmonics as in [8].

Conclusions
The difference between the paper [2] and the present one is the presence of detailed Boussinesq formulation and the analytical expressions for the second-order super-and sub-harmonic wave amplitudes of the BM and more wave period cases in CFD simulation under the consideration of turbulent flow. Furthermore, one more comparison result for the wave gauge behind the cylinder between Boussinesq and CFD simulations is studied, the CFD results for wave forces, moment, and pressure distribution around the cylinder are presented and discussed, and finally, a new section is added to investigate the effect of second-order sub-harmonic waves on the cylinder for different design parameters from BM.
The innovation of this study is the development of BM based on BEs and numerical CFD simulations associated with the vertical truncated cylinder and providing benchmark results by comparing various methodologies of non-linear/second-order super-and sub-harmonic incoming wave amplitudes/wave heights interactions, as these phenomena have central importance during the design of WECs. On the other hand, the horizontal and vertical forces and pitch moment acting on the vertical truncated cylinder are newly added to this study. It has been observed that: 1.
The CFD results generated using OpenFOAM for wave heights/wave elevations showed reasonable agreement with experimental and other numerical models for all the cases. CFD simulations also predicted results for horizontal and vertical forces and moments encountered by the cylinder.

2.
The results of second-order super-harmonic wave amplitude were compared with CFD results and it was found that they are in good agreement and the trends are also well captured in all cases.

3.
The effects of second-order super-harmonic waves are stronger than those of sub-harmonic waves around the cylinder, and the complicated transformation process in the wave amplitude distribution is found to be due to the phase interaction between the flows around and below the cylinder.

4.
The CFD pressure distribution indicates that the wave elevation at the side and behind of the cylinder is lower than that at the front. From the computational point of view, the present study indicates that the CFD is far more expensive compared to BM and perhaps is best suited for final design validation for WEC devices, rather than the design iteration stage.

5.
As a result, the present BM can be used with confidence to generalize for heaving motion characteristics of the floating vertical truncated cylinder to design point absorber floating WEC. 6.
Furthermore, the wave force results of the CFD modeling may help to use in the future BM model to analyze the level of fidelity by comparing them.