Abstract
In this paper, the nonlinear interaction of regular water waves propagating over a fixed and submerged circular cylinder is numerically studied. At the structure’s lee side, the free surface profile experiences strong nonlinear deformation where the superharmonic free wave generated can be significant and is superposed on the transmitted wave. The wave profile then becomes asymmetric and skewed and may eventually reach the point of physical wave breaking. The governing equation and boundary conditions of this wave–structure interaction problem are formulated using both the fully nonlinear and the weak-scatterer theory. The corresponding boundary value problem is numerically solved by the immersed-boundary adaptive harmonic polynomial cell solver. In this study, a pragmatic wave-breaking suppression model is incorporated into the original solver. Both the harmonic free wave amplitudes at the structure’s lee side and the harmonic vertical forces on the cylinder are studied. The simulated harmonic wave amplitudes are compared to other published experiments and theoretical data. In general, good agreement is achieved. The effects of the incorporated wave-breaking suppression model on the simulated results are discussed. In our study, the incorporation of the pragmatic wave-breaking suppression model successfully extends the capabilities of the original fully nonlinear immersed-boundary adaptive harmonic polynomial cell solver.
1 Introduction
The generation and evolution of ocean waves as well as their interaction with realistic bathymetry and engineering structures involve complex kinematics and dynamics of the waves that are subjected to nonlinear, dispersive, and dissipative effects. Governed by these effects, even waves with mild steepness may undergo significant deformation when propagating over underwater obstacles, such as shallow reefs, ridges, or submerged marine structures. A set of free superharmonic waves is then generated at the lee side of the obstacles due to nonlinear free surface effects, accompanied by a transfer of wave energy from lower to higher order components. These free waves are of shorter wavelengths than the transmitted first-order wave, due to the dispersion relation being quadratic in wave frequency. The deformed waves may eventually break in the form of, e.g., spilling or plunging breakers. Once the wave breaking occurs, it contributes the most to wave energy dissipation [1].
For simulating wave processes in which viscous effects are dominant and wave-breaking details are desired, models based on Navier–Stokes (N–S) equations can be utilized [2,3]. While the N–S models involve more physics, simulations using these models are still computationally intensive, which makes them less suitable to be extensively applied in the early design stages of marine structures. When viscous effects become of second importance, linear potential-flow models such as WAMIT [4] and ANSYS-AQUA [5] are less computationally demanding and widely used in industry to provide a good approximation. Models that directly solve the fully nonlinear potential-flow (FNPF) equations are further developed to simulate problems where strong nonlinearities are involved [6–9]. Among these methods, the immersed-boundary adaptive harmonic polynomial cell (IB-AHPC) method proposed by Tong et al. [9] has been validated and verified to accurately and efficiently solve FNPF wave–structure problems in two dimensions. Generally, FNPF models are only valid up to wave breaking. A consequence of using these models is the possibility that the simulated wave steepens to the point of physical wave breaking and causes the numerical program to be unstable and eventually terminate.
However, when wave breaking is either local or not of primary concern, it is of practical interest to extend an existing FNPF model to study waves that are close to breaking without actually modeling the overturning wave. The aforementioned numerical instability can be suppressed by specifying artificial energy dissipation caused by wave breaking. Various breaking onset criteria have been proposed to predict the commencement of the breaking process and to instruct when and where the breaking suppression should be applied. These breaking criteria can be broadly classified into three categories geometric, kinematic, and dynamic ones [10]. More recent works suggested that there appears to be a universal breaking onset criterion, see Refs. [11–13]. Once the breaking onset is detected, a proper energy dissipation should be introduced, either in a physical manner using an absorbing surface pressure in the dynamic free surface boundary condition [13] or in a pragmatic way where the incipient breaking wave profile is smoothed by a special local filter, e.g., as implemented in the open-source code OceanWave3D [14,15]. Though less sophisticated than a physical wave-breaking suppression model, the implementation of a breaking filter is much more straightforward and effective. However, there still are uncertainties about the selection of values of a specified breaking criterion. Besides, the amount of energy dissipation introduced by the breaking filter is less controllable. The influence of these uncertainties on the FNPF wave–structure interaction problems is not fully explored and understood. This gap emphasizes the need for further research, leading us to examine these effects more comprehensively.
In this paper, we numerically simulate the interaction of regular water waves propagating over a fixed and submerged circular cylinder, which has been experimentally studied by Grue [16]. In this case, the Keulegan–Carpenter (KC) number [17,18] is small (), indicating that inertial forces dominate over viscous forces, making viscous effects insignificant. Due to nonlinear effects, significant deformation of waves was observed at the lee side of the cylinder. Therefore, we use the FNPF IB-AHPC method, which makes no assumptions about nonlinearity and is of accuracy and effectiveness, to perform the numerical simulations. The corresponding boundary value problems (BVPs) and the numerical schemes are described in Secs. 2.1–2.4. In Sec. 2.5, we incorporate a pragmatic wave-breaking suppression model into the original FN IB-AHPC solver, aiming to extend its capabilities that enable more simulation runs after physical wave breaking occurs. This breaking suppression model is the same as the one implemented in the open-source code OceanWave3D [14,15]. The simulation results and the influence of values of a specified breaking criterion factor on these results are presented and discussed in Sec. 3. Finally, conclusions are drawn in Sec. 4.
2 Wave–Structure Interaction Model
2.1 Fully Nonlinear Potential-Flow Formulation.
2.2 Weak-Scatterer Approximation.
Alternative to the FNPF model, the weak-scatterer (WS) approximation [19] is also able to solve higher than second-order wave–structure interaction problems [20–22]. In this method, the scattered and radiated waves are assumed to be small compared to the incident waves. The FN free surface boundary conditions (FSBCs) can then be linearized with respect to the incident wave. This linearization generally leads to more stable and robust numerical simulations. In this work, the WS IB-AHPC solver developed by Tong et al. [22] is also used to solve a wave–structure interaction problem with the studied cylinder submerged in water. The results obtained from the WS approximation can be used as reference values for cases where large wave amplitudes are involved and cause an early termination of the original FN IB-AHPC solver. Here, we only summarize the principles of this theory.
The no-flux boundary conditions on the surface of the fixed body and the tank bottom are .
To avoid unwanted reflection waves from the end walls of the numerical tank, a damping zone is applied in both FN and WS formulations to simulate tank beach.
Following both the FN and WS theories, two mixed Dirichlet–Neumann BVPs are established, respectively, to be solved for the velocity potential at each time-step. In this paper, the fourth-order explicit Runge–Kutta scheme is used for the time marching of the free surface boundary conditions. The selective filter [24] and the classical 13-point 10th-order Savitzky–Golay (SG) filter [25] are applied near the side walls of the numerical tank and the rest of the free surface, respectively, to reduce the numerical saw-tooth-type instability.
2.3 The Two-Dimensional Immersed-Boundary Adaptive Harmonic Polynomial Cell Method.
The harmonic polynomial cell (HPC) method is an efficient and accurate numerical method originally proposed by Shao and Faltinsen in both 2D [26] and 3D [27] to solve potential-flow boundary value problems governed by the Laplace equation. In this work, the newly developed IB-AHPC method, validated and verified for fully nonlinear problems [9,22], is adopted. Based on the original HPC method, the IB-AHPC is implemented with the immersed-boundary method (IBM) in adaptive cells. Within the computational domain, the adaptive cells allow local grid refinement only near areas of interest, such as areas near the free surface and body structures. This adaptive refinement enables the computing procedure to gain high efficiency without compromising accuracy. The application of IBM makes it straightforward and easy to model moving structures with complex boundaries. A brief summary of the principles of the IB-AHPC method is given in the following sections.
2.3.1 The Two-Dimensional Harmonic Polynomial Cell Method.
As a field solver, this method first discretizes the computational domain into overlapping cells. Each cell consists of four neighboring quadrilateral elements with eight grid points (with indices ) on its boundary and one (with indices ) located at the center of the cell, as illustrated in Fig. 2.

Sketch of global overlapping quadrilateral cells and a standalone cell with its local indices and coordinate system
Solving this linear equation system and relabeling the suffices and , we obtain the superposition coefficient as follows:
Having the continuity constraints and the boundary conditions been established by Eqs. (17), (16), and (18), a global linear algebraic equation system can be obtained. The resulting global algebraic system is a sparse matrix, which has at most nine nonzero entries in each row. This linear equation system can be solved efficiently by choosing a proper preconditioner and an iterative solver, e.g., using Incomplete LU Decomposition with Threshold (ILUT) factorization for the preconditioning before solving the equation system with a Generalized Minimum Residual (GMRES) solver.
2.3.2 The Immersed-Boundary Method.
The IBM was first introduced to the HPC method by Hanssen et al. [28] for simulating a moving body in an infinite fluid domain. Later, this method was enhanced to model the free surface in a numerical wave tank (NWT), with various wave propagation problems examined to demonstrate its performance [29]. With the promising ability of the IBM, Hanssen [30] used this method to successfully study the potential-flow fully nonlinear wave-body interaction. Then, Tong et al. combine IBM with a quadtree-based adaptive grid refinement strategy, which strikes a balance between computational accuracy and efficiency. This IB-AHPC method has been successfully applied to fully nonlinear wave-body interaction problems involving moving boundaries and to weakly nonlinear wave-body interaction problems formulated by a WS approximation [9,22].
In our work, we adopted the IB-AHPC methods to solve the nonlinear wave-body interaction problems. The free surface and the structure are modeled as an immersed boundary as illustrated in Figs. 3 and 4, respectively. The free surface is represented by sets of tracking markers (circles with black stroke as in Fig. 3), which are used for imposing the Dirichlet boundary condition through interpolation by Eq. (16). The marker coordinates and the velocity potential at this location are obtained by the time integration of the FSBCs in Eqs. (3) and (4). This interpolation achieves its best accuracy along the middle line of harmonic polynomial cells as suggested by Ma et al. [31]. Hence, the position of a free surface marker is chosen to be located on the vertical grid lines of the computational domain. Those grid lines coincide with the vertical central line of interpolation cells to which free surface markers correspond. Moreover, those markers are distributed on the upper half of their interpolation cells which are defined as a free surface cell. In this work, the description of motions of free surface markers is thus within the semi-Lagrangian framework, i.e., the prescribed velocity in the material derivative operator is taken as , see Eq. (2).

Free surface modeled as an immersed boundary in Cartesian grids. The circles without outlines are fluid nodes in the computational domain. The circles with outlines are wave markers. The squares represent the free surface ghost node (including additional ghost nodes).

Body surface modeled as an immersed boundary. The squares are body ghost nodes. The circles on the body surface are discretized body points. The other types of nodes are the same as in Fig. 3. The black arrow represents the body’s normal vector pointing into the fluid.

Body surface modeled as an immersed boundary. The squares are body ghost nodes. The circles on the body surface are discretized body points. The other types of nodes are the same as in Fig. 3. The black arrow represents the body’s normal vector pointing into the fluid.
As the free surface is modeled as an immersed boundary, ghost nodes need to be found for the velocity potential interpolation at free surface markers. For a marker and its free surface cell, the node with local index 7 is denoted as its ghost node, referring to Fig. 2. However, there is another type of ghost node that cannot be found by index 7, see node-1, -2 and cell-1, -2 in Fig. 3. These nodes are defined as additional ghost nodes and they are used to make sure that the number of unknowns and equations are the same. For an additional ghost node upon which the Dirichlet boundary condition is imposed, it uses different interpolation cells from the interpolation cell its neighboring ghost node used. This neighboring ghost node shares the same horizontal coordinate as the additional ghost node. Referring to Fig. 3 as an example, additional ghost node-1 uses cell-1 as its interpolation cell, while the ghost node just below node-1 uses the cell at the right-below position of cell-1 as its interpolation cell. Though the Dirichlet boundary condition is imposed twice for the same free surface marker below node-1, due to different interpolation cells being used, the singularity in the matrix of the final algebraic equation system can thus be avoided.
Similar to the free surface, body ghost nodes and corresponding body ghost cells are to be found around the structure surface, see Fig. 4. Neumann boundary conditions are enforced on discrete body points through Eq. (18). The number of the body ghost nodes is equal to the boundary equations established, which are ensured by finding each body ghost node a closest body point including its coordinates and unit normal vector. This matching procedure for ghost nodes is accomplished by means of a cubic B-spline interpolation on the body surface.
For fluid nodes, i.e., the nodes located within the computational domain but not in its boundaries, the connectivity equation (17) will be used to ensure the continuity of the velocity potential.
2.3.3 Adaptive Cartesian Grid System.
An adaptive Cartesian grid algorithm is implemented in this work following Refs. [9,22] to refine the local discretization of the computational areas where the gradients of velocity potential are of interest. These areas, near the free surface and the body surface, are modeled as immersed boundaries; several layers of the Cartesian grid system are established as required, with cell sizes being gradually coarser away from the interested boundaries. The global number of unknowns is reduced by this dynamical cell refinement procedure, which strikes a balance between computational efficiency and accuracy. An illustration of the adaptive Cartesian grid system including the free surface and the submerged body boundary is depicted in Fig. 5.

A portion of the immersed boundary in adaptive cells, including still water surface and a submerged cylinder. The maximum quadtree refinement level is 2 here. The meanings of the different nodes are the same as those mentioned in Figs. 3 and 4. In addition, there are fluid nodes represented by triangles (defined as type-1 nodes), circles in between two adjacent triangles (type-2 nodes). These two kinds of nodes are used for information interchange between different quadtree layers.

A portion of the immersed boundary in adaptive cells, including still water surface and a submerged cylinder. The maximum quadtree refinement level is 2 here. The meanings of the different nodes are the same as those mentioned in Figs. 3 and 4. In addition, there are fluid nodes represented by triangles (defined as type-1 nodes), circles in between two adjacent triangles (type-2 nodes). These two kinds of nodes are used for information interchange between different quadtree layers.
2.4 The Lagrangian Acceleration Potential Method.
In our work, is obtained using a Lagrangian acceleration potential following Ref. [23]. For a fixed body, . To obtain , we first define a BVP for . On the free surface, the Dirichlet boundary condition . On the body surface, the impermeable condition is . Once the above BVP defined for is solved, the time derivative of the velocity potential is also obtained.
2.5 Wave-Breaking Suppression Model.
While the occurrence of wave breaking is apparent when the free surface profile becomes disturbed or collapsed, the mechanism leading to the breaking is elusive. Hence, as mentioned in Sec. 1, various breaking onset criteria [10–13] were proposed to help detect the commencement of a breaking process. Once the breaking onset is detected according to the specified criterion, different methods can be utilized to introduce the breaking dissipation [13,14]. In our work, we directly follow Schløer [14] and the open-source code OceanWave3D [15] for selecting the breaking onset criterion and the breaking suppression method.
The exact amount of energy dissipated by this filter cannot be controlled. However, more dissipation can be introduced as it is desired to simulate cases with very steep waves. This can be done, for example, by decreasing the time-marching step to increase the number of filtering operations per unit time.
3 A Fixed Submerged Circular Cylinder in Regular Waves
3.1 Numerical Setup.
In this section, the nonlinear interaction of regular waves propagating through a fixed circular cylinder located under the free surface is simulated in a NWT, see Fig. 6. The coordinate system is defined as described in Sec. 2.1. This numerical wave tank is 12 m long and filled with water to a depth of . The studied cylinder is placed at the longitudinal center of the NWT, which coincides with the z-axis of the coordinate system. The radius of the cylinder is . The distance between the uppermost point of the cylinder and the mean free surface, i.e., the submergence , is either 50 mm or 100 mm, depending on the case.

Schematic representation of the numerical wave tank with a fixed submerged circular cylinder under waves
The incident wave is described by the stream function theory [33] with amplitude varying between and . There are two wave frequencies selected as either or , with corresponding wavelengths being and , and the maximum wave steepness for the two waves with different frequencies are approximately and , respectively. The simulated wave profile is described by sets of free surface markers. The movement of these markers is tracked in a semi-Lagrangian manner, with the prescribed velocity given as , see Eq. (2). Thus, these free surface markers are only allowed to move vertically, which coincides with the vertical grid lines of the deepest quadtree grid layer, as in Fig. 5.
3.2 Numerical Convergence Study.
Numerical tests were carried out to assess the convergence of spatial and temporal discretization. The cases with an incident wave frequency , wave amplitude , and cylinder submergence were considered. Figure 7(a) shows different snapshots of free surface elevation along the numerical wave tank taken at an instant time , which were obtained using a coarse mesh (20 elements per cylinder diameter , i.e., , in the deepest quadtree layer), an intermediate mesh (), a fine mesh (), respectively. A small time-step of () is used. It can be seen that the three obtained results clearly converge, and the results from the intermediate and fine meshes are almost identical. For the vertical force on the cylinder, the convergence of results can also be observed, as depicted in Fig. 7(b). Hence, the intermediate mesh is used in the following study.

Comparison of different space-step results, with incident wave amplitude mm and wave frequency Hz. The submergence of the cylinder is mm. The time-step . (a) Snapshot of free surface elevation along the NWT at and (b) time history of vertical wave force on the cylinder.
Likewise, tests were conducted using the intermediate spatial mesh and different time-steps (). Figure 8(a) shows the obtained time histories of free surface elevation recorded at a probe location , while Fig. 8(b) shows the obtained time histories of the vertical wave force on the cylinder. Both sets of results clearly indicate that convergence is achieved for the tested time-steps. For the following study, at least an intermediate time-step of is used.

Comparison of different time-step results, with incident wave amplitude mm and wave frequency Hz. The submergence of the cylinder is mm. The space-step . (a) Time history of free surface elevation recorded at and (b) time history of vertical wave force on the cylinder.
3.3 Harmonic Wave Amplitudes Around the Cylinder.
In Grue’s work, it was suggested that the free harmonic amplitudes obtained with this method are almost independent of the values of and . The exact position where the wave gauges were placed in his experiment was not given. In our numerical simulation, the surface elevation was obtained by the wave markers tracking the free surface profile. For runs with a relatively small incident wave amplitude, we observed similar situations where the wave amplitudes obtained using this method were almost independent of the values of as in Fig. 9(a). However, for runs with a larger incident wave amplitude, fluctuations of the obtained amplitudes were observed with the increasing values of , i.e., going downstream away from the cylinder’s lee side, as shown in Fig. 9(b). Hence, additional treatment of averaging the obtained surface harmonic free wave amplitudes along a given distance was introduced in our simulation. This distance was given as in length and began at one wavelength downstream away from the center of the cylinder as illustrated in the shaded area in Figs. 9(a) and 9(b). This averaging method is used in the following studies. The influence of the value of was relatively small as observed also in our work. We set for the following studies, where is the incident wavelength.

Normalized amplitudes of the first three harmonic components of the free surface elevation. Here, these results are for the case with the incident wave frequency Hz, and the submergence mm. The horizontal coordinate represents the distance from the cylinder center at its lee side, normalized by incident wavelength . (a) mm and (b) mm.

Normalized amplitudes of the first three harmonic components of the free surface elevation. Here, these results are for the case with the incident wave frequency Hz, and the submergence mm. The horizontal coordinate represents the distance from the cylinder center at its lee side, normalized by incident wavelength . (a) mm and (b) mm.
The first and second harmonic free waves obtained by our FN and WS methods are presented in Fig. 10, for the circular cylinder with submergence mm, i.e., the distance between the free surface and its uppermost point. The incident wave amplitudes vary from to mm. Two wave frequencies are simulated, namely Hz and Hz, which are shown in Figs. 10(a) and 10(b), respectively. For comparison, the experimental results of Grue [16] and the second-order results of Friis [35] are also shown. In general, the FN and WS simulations of our work agree well with those reference results, despite small discrepancies occurring when the incident wave amplitude increases.
For the first harmonic free wave amplitudes, the FN, WS, and experimental results compare well with each other for the runs for which the incident amplitude is relatively small. For runs with larger incident amplitude, the FN results deviate from the WS, but still are in good agreement with the experimental results. The results obtained from the experiment and the FN solution are slightly smaller than unity in these runs. It can be observed in Fig. 10(b) that approximately after the wave amplitude exceeds mm, the first-order amplitude gradually decays. Since there is no wave breaking observed in these cases, the reason for the decay can be partly inferred as the energy transfer to the second harmonic free wave and the energy dissipation due to viscosity.
For the second harmonic free wave amplitudes, similar results can be found, where those results compare well with each other in runs with small . The results obtained from the experiment, the FN and WS solutions are increased linearly for the smaller incidents wave cases and are much the same as the results obtained from Friis’s [35] second-order theory. As the incident wave amplitude gets larger, the three results clearly deviate from the second-order results, while they still reach a good agreement with each other.
Then, we further consider the cases where the distance between the free surface and the uppermost point of the cylinder, i.e., the submergence , is reduced to mm. The comparisons of the first and second harmonic free wave amplitudes and are shown in Fig. 11, for the incident wave frequency Hz, and in Fig. 12, for the incident wave frequency Hz. For comparison, the experimental results of Grue [16] and the second-order results of Friis [35] are included as benchmarks. Spilling and plunging breakers were observed in Grue’s experiment study, and we use arrows annotated with the corresponding initials (S) and (P) to denote these limits, as in Figs. 11 and 12.

As Fig. 10(a), with submergence mm, Hz. The arrows denote the spilling (S) and plunging (P) limits.

As Fig. 10(a), with submergence mm, Hz. The arrows denote the spilling (S) and plunging (P) limits.

As Fig. 10(b), with submergence mm, Hz. The arrows denote the spilling (S) and plunging (P) limits.

As Fig. 10(b), with submergence mm, Hz. The arrows denote the spilling (S) and plunging (P) limits.
For runs with small incident wave amplitudes, the results obtained by the second-order theory [35], the FN and WS solutions are generally in good agreement with the experiment data, which can be observed in Figs. 11 and 12. It can be seen that in the two cases, the values of normalized first harmonic free waves obtained from the experiment decay more significantly with increasing incoming amplitude than cases shown in Figs. 10(a) and 10(b), as nonlinear effects become more significant due to the smaller submergence. Hence, in addition to energy dissipation caused by wave breaking and viscous effects, energy transfer from first harmonic to second harmonic are accounted important for this decay.
For the normalized second harmonic wave, obtained from the FN solution agrees well with the experiment results, as can be seen in Fig. 11. Since the second-order theory cannot accurately describe the energy transfer from first to second harmonic free wave, as reported by Grue’s work [16], deviations between the results of second-order theory and the FN results can be observed, when the incident wave amplitude gets larger, which is shown in Figs. 11 and 12. We further notice that in the case with a small incident wave frequency as in Fig. 11, the FN results agree well with the second-order theory results and the experiment data, while in the case with a larger incident wave frequency as in Fig. 12, small discrepancies can be observed between the FN and the experiment results. In Fig. 12, the FN and second-order results still agree for runs with smaller incident amplitudes, but both results deviate from the experiment data. This may partly result from the first-order results that are much smaller than unity even when the incident wave amplitudes are quite small. However, the explanation for the mechanisms of this phenomenon may require further studies.
It is worth noting that, even with the simulated wave amplitudes exceeding the spilling limits, the FN results are consistent with the experiment one, for example, in Fig. 11. For small incident wave amplitude, the second-order results obtained is a monotonically increasing function. The curve based on second-order theory will keep increasing even for larger steepness, no matter whether the waves break or not. At the observed spilling limit, the normalized second harmonic wave reaches its maximum value, which is approximately . Then, for larger incoming waves, turns to be a monotonically decaying function.
It seems that overall agreement with the experiment data is achieved by the WS when the incident wave amplitude is relatively small, as in Figs. 11 and 12. However, with incident wave amplitude exceeding the spilling limit, results inconsistent with the experiment data are obtained using this method. A closer look at the last several data points of WS in Fig. 12 shows that increases and even exceeds unity. A similar phenomenon is observed for the first harmonic vertical force results in our work, see Fig. 12, and in Ref. [22], which are more apparent. From their work, it is suggested that this may result from the increasing of the scattered wave elevation which makes itself no more a small value compared to the incident wave elevation.
For the WS results, large discrepancies from the experiment data can be observed in Figs. 11 and 12. It is therefore reasonable to suppose that as reduced models which neglect nonlinearity to some extent, both Friis’s second-order theory [35] and the WS formulation cannot completely describe the energy transfer from lower harmonic free waves to higher ones due to nonlinear effects. It appears that the values of obtained by WS will eventually peak at the case whose incident wave amplitude exceeds the breaking limit, though the corresponding case is much delayed as compared to the observation in its FN and experiment counterparts; then these values of cease to increase, i.e., to maintain or decrease.
Interestingly, this second-order inconsistency seems eventually to have no significant influence on the dynamic forces exerted on the cylinder, as can be seen in Figs. 15 and 16. However, a nonphysical increase in the values of both the first-order wave amplitude and vertical forces is observed as the incident wave amplitude becomes larger. Nonetheless, it is suggested that the basic assumptions of this theory must not be violated, if trustworthy results are desired [22].
Generally, a fully nonlinear potential-flow model remains valid up to physical wave breaking. However, if implemented with no other treatments, the computation procedure usually experiences numerical instabilities, due to very steep waves, and terminates when actual wave breaking occurs. In this work, simulations using our original FN IB-AHPC solver without implementation of a wave-breaking suppression model will terminate soon, as expected, after the amplitude of waves at the cylinder’s lee side exceeds the spilling limit. This termination can be observed in Figs. 11 and 12. It should be noted that for dealing with the saw-tooth numerical instability, a classical 13 point 10th-order SG filter is utilized, as described in Sec. 2.1. It is supposed that implementing a wave-breaking suppression model which somehow introduces more dissipation, either being physical or nonphysical, the computation capability of a fully nonlinear potential-flow model may then be extended.
By implementing the pragmatic wave-breaking suppression model, as described in Sec. 2.5, the computation capability of the presented FN IB-AHPC solver is extended. To demonstrate, results obtained by the FN model with/without implementing the wave-breaking suppression model are compared with the experiment data as shown in Figs. 13 and 14. For the simulations using the wave-breaking suppression model, two wave-breaking criteria are used, namely and . The selection of a smaller breaking criterion means an early commencement of the breaking suppression procedure, which results in more dissipation. The use of such a wave-breaking suppression model clearly enables the previously terminated simulations; and a smaller breaking criterion enables more runs than those with a higher criterion, as expected.

As Fig. 11 and compared with results from using different breaking criterion . Submergence mm, Hz. The arrows denote spilling (S) and plunging (P) limits.

As Fig. 11 and compared with results from using different breaking criterion . Submergence mm, Hz. The arrows denote spilling (S) and plunging (P) limits.

As Fig. 12 and compared with the results from using the wave-breaking suppression model. Submergence mm, Hz. The arrows denote the spilling (S) and plunging (P) limits.

As Fig. 12 and compared with the results from using the wave-breaking suppression model. Submergence mm, Hz. The arrows denote the spilling (S) and plunging (P) limits.
More strikingly, as shown in Figs. 13 and 14, the dissipation introduced by the wave-breaking suppression model is reasonably small and almost independent of the specified breaking criterion factor so that the extended results are in line with the corresponding experiment data; hence, the free surface kinematics are not sabotaged by the additional dissipation introduced. To enable more simulation runs than those presented in this paper, which require more dissipation, one can increase the number of smoothing operations in a given time, for example, by decreasing the simulation time-step. This encouraging result motivates further examinations on the influence of this wave-breaking suppression upon the wave forces exerting on the cylinder, which will be presented in the next section.
3.4 Harmonic Wave Forces on the Cylinder.
Chaplin [18] conducted an experiment with a cylinder of diameter m, fixed under the still water surface with submergence in a wave channel. The wave channel was filled in water with a depth of m. The period of the incident wave is s, with wave height varying from m to m. Recently, using the FN and WS formulation with the IB-AHPC method, Tong et al. [22] numerically simulated this experiment. Detailed discussion and thorough comparisons of the dynamic harmonic vertical forces on the cylinder were made between their simulated results and the experiment data. In addition, results obtained by other theoretical results were also compared. The ability of the WS and FN IB-AHPC methods to accurately calculate the dynamic forces in the wave–structure interaction problems is thus validated and verified. For details, we refer readers to Sec. 5.1 in Ref. [22].
In this work, a similar but different case concerning the harmonic free wave amplitude, as described in the previous section, is simulated and compared to the experiment conducted by Grue [16], using the same FN and WS IB-AHPC solvers as in Ref. [22]. But we further incorporated a pragmatic wave-breaking suppression model into the original FN IB-AHPC solver to successfully enable more simulation runs in cases where large amplitudes could have led to an early termination of the original numerical procedure.
The results obtained in the previous section show that the incorporation of the pragmatic wave-breaking suppression model into the original FN IB-AHPC solver successfully extends the solver’s capability which enables the simulation even after the physical wave breaking occurs. It is shown in Figs. 13 and 14 that the incorporated wave-breaking suppression model does not sabotage the free surface kinematics in terms of harmonic free wave amplitude; the extended results are in line with experiment data. In this section, we continue to study the dynamic forces on the same setup as in Figs. 13 and 14.
The FN results of the vertical force harmonic components are presented in Figs. 15 and 16, corresponding to the same test cases in Figs. 13 and 14, respectively. For comparison, the WS results are shown as reference results, since, in our work, the cases whose wave breaking were observed in the experiment do not lead to an early termination of the WS numerical program. In general, the FN results are in good agreement with the WS results, especially for the zeroth-, first-, and second-order harmonic vertical force components. Though, for the third-order harmonics, small discrepancies are observed, it is still acceptable as the WS, being a weakly nonlinear theory, inheres less nonlinearity compared to its FN counterpart. Similar findings for the comparison between the FN and WS results are reported in Ref. [22].

Normalized amplitude of the various harmonic components of the vertical wave force on the circular cylinder. The results in the shaded area are obtained using the wave-breaking suppression model, with breaking criterion factor and . The incident wave frequency Hz, and the submergence of the cylinder mm.

Normalized amplitude of the various harmonic components of the vertical wave force on the circular cylinder. The results in the shaded area are obtained using the wave-breaking suppression model, with breaking criterion factor and . The incident wave frequency Hz, and the submergence of the cylinder mm.

Normalized amplitude of the various harmonic components of the vertical wave force on the circular cylinder. The results in the shaded area are obtained using the wave-breaking suppression model, with breaking criterion factor and . The incident wave frequency Hz, and the submergence of the cylinder mm.

Normalized amplitude of the various harmonic components of the vertical wave force on the circular cylinder. The results in the shaded area are obtained using the wave-breaking suppression model, with breaking criterion factor and . The incident wave frequency Hz, and the submergence of the cylinder mm.
Besides the results obtained from using the same numerical routines as in Ref. [22], additional results are obtained by utilizing a slightly different FN IB-AHPC solver, namely the solver incorporated a pragmatic wave-breaking suppression model as in Sec. 2.5. We first conducted simulations without implementation of the wave-breaking suppression model; the result is only available up to the KC number indicated by the left boundary of the shaded area in Figs. 15 and 16. For a larger KC number, this numerical procedure inevitably terminates. Then, simulations with the implementation of the wave-breaking suppression model are conducted. Two alternative values of breaking criterion factor are specified, namely and . A smaller factor represents an early commencement of the wave profile smoothing procedure, which brings more dissipation into the numerical simulation. More simulation runs are thus made available through this extended FN solver, as indicated by the shaded area in Figs. 15 and 16. It is encouraging that, in this extended area, the zeroth-, first-, and second- order components of the extended FN results conform well with the WS solutions and thus the extended line of the original FN results is obtained without implementing such a wave-breaking suppression model. The extended line is not drawn on the figures, as the FN and WS results before the extended area already agree well in these harmonic components. For the third-order components, the extended FN results are slightly lower than the extended line of the original FN results. However, it is still considered acceptable, as these values appear to gradually approach their WS counterparts.
It is noticed that a smaller value of breaking criterion allows more simulation runs with larger KC numbers. And the specific value of this criterion has limited influence on the magnitude of the higher order harmonic forces. To enable more runs with even larger KC numbers than the cases presented in this paper, more dissipation can be further introduced. For example, one can decrease the time-marching step to increase the number of smoothing operations per unit time.
4 Conclusion
In this paper, we have numerically studied the nonlinear interaction of regular water waves propagating over a fixed and submerged circular cylinder. The BVP is formulated within the framework of the potential-flow theory in both FN and WS manners. The established BVPs are solved by the IB-AHPC method, where the free surface and body boundaries are immersed in an adaptively refined Cartesian quadtree grid system. To further extend the computational capability of the original FN IB-AHPC method, we incorporate a pragmatic wave-breaking suppression model into the solver to stabilize the numerical program even after physical wave breaking occurs. The free surface is tracked using wave makers by a semi-Lagrangian approach. The explicit fourth-order Runge–Kutta time integration scheme is used to update the free surface boundary conditions.
Both the harmonic free wave amplitudes at the structure’s lee side and the harmonic vertical forces on the cylinder are studied. Generally, the FN and WS simulation results of our work agree well with other experimental and theoretical works. The extended FN IB-AHPC solver provides trustworthy results in terms of both wave kinematics and dynamics after wave breaking occurs. We found that different choices of the breaking criterion factor have limited influence on the extended results. In other words, the magnitudes of the extended results are almost independent of the breaking criterion factor. Therefore, we conclude that the incorporation of the pragmatic wave-breaking suppression model successfully extends the computational capabilities of the original FN IB-AHPC solver. More dissipation can further be introduced if more cases with larger wave amplitudes are desired. For example, one can decrease the time-marching step to increase the number of the smoothing operations involved in the numerical simulation.
Future work should extend the wave-breaking suppression model to handle three-dimensional wave–structure interactions, with a focus on making the amount of energy dissipation more controllable.
Acknowledgment
The authors gratefully acknowledge Dr. Chao Tong for providing the FN and WS IB-AHPC solvers. The authors Qihao Wu, Yujing Chen, and Min Zhang are supported by the Key R&D Program of Shandong Province, China (2021ZLGX04) and the 10.13039/501100001809 National Natural Science Foundation of China (Nos. 52088102 and 52171281).
Conflict of Interest
There are no conflicts of interest.
Data Availability Statement
The datasets generated and supporting the findings of this article are obtainable from the corresponding author upon reasonable request.
Nomenclature
- =
gravitational acceleration
- =
wave number
- =
normal vector
- =
displacement vector
- =
prescribed velocity
- =
diameter of the cylinder
- =
radius of the cylinder
- =
incident wave period
- =
body forces due to dynamic pressure of the fluid
- =
body moments due to dynamic pressure of the fluid
- =
Lagrangian downward vertical acceleration
- =
wavelength
- =
length of the wave-making zone
- =
length of the wave-damping zone
- =
normal components of the rigid-body velocity
- =
submergence of the cylinder
- =
harmonic free wave amplitudes
- =
global coordinates
- =
local coordinates in the harmonic polynomial cell