Fully Nonlinear Interaction of Water Waves and a Submerged Cylinder With Wave Breaking Detection and Suppression

Wu, Qihao; Chen, Yujing; Zhang, Min; Shao, Yanlin

doi:10.1115/1.4067502

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 (⁠ $≪ 1$ ⁠), 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.

In this work, the wave–structure interaction problems are formulated in 2D. The computational domain is assumed to be filled with a homogeneous, inviscid, and incompressible fluid, and the flow is irrotational. The fluid particle velocity can then be described by the gradient of the resulting velocity potential

φ (x, z, t)

⁠. A Cartesian coordinate system

O x z

is defined with its

O z

axis oriented positively upward and its

O x

axis oriented positively rightward coinciding with the undisturbed free surface. The origin is located on the horizontal center plane of the numerical tank, as depicted in Fig. 1. The Laplace equation is applied as the governing equation, that is

\nabla^{2} φ (x, z, t) = 0

(1)

A semi-Lagrangian method is adopted to track markers representing the instantaneous free surface using a fully nonlinear (FN) manner. These markers on the free surface follow the material time derivative defined as

\frac{δ}{δ t} = \frac{\partial}{\partial t} + v \cdot \nabla

(2)

where

v

is a prescribed velocity. Thus, the fully nonlinear free surface boundary conditions in the inertial reference frame are given as

\begin{aligned} {\begin{matrix} \frac{δ x}{δ t} = \frac{\partial φ}{\partial x} + (v - \nabla φ) \cdot \nabla x \\ \frac{δ η}{δ t} = \frac{\partial φ}{\partial z} + (v - \nabla φ) \cdot \nabla η - ν (x) η \end{matrix} on z = η \end{aligned}

(3)

\frac{δ φ}{δ t} = - \frac{1}{2} {| \nabla φ |}^{2} - g η + v \cdot \nabla φ - ν (x) φ on z = η

(4)

where

g

is the gravitational acceleration,

\nabla x = (1, 0)

⁠,

\nabla η =

(\partial η / \partial x, 0)

⁠, and

v = (0, \partial η / \partial t)

⁠. Equation is derived from Bernoulli’s equation assuming that the free surface pressure is equal to the constant atmospheric pressure.

Fig. 1

View large Download slide

Sketch of the wave-body interaction model in a numerical wave tank

To solve the governing Laplace equation, the computational domain is either enclosed by Dirichlet boundaries

Γ_{D}

or Neumann boundaries

Γ_{N}

⁠, where corresponding boundary conditions are imposed. The two kinds of boundary conditions are generally defined as follows, respectively:

φ (x, z) = φ_{D} on Γ_{D}

(5)

\frac{\partial φ (x, z)}{\partial n} = V_{N} on Γ_{N}

(6)

where

V_{N}

is the normal component of the rigid-body velocity. Specifically, for a fixed rigid boundary, e.g., tank bottom, body surface, etc., the no-penetration or no-flux condition shall be imposed so that the right-hand side of Eq. on those boundaries equals zero. For the free surface, a Dirichlet boundary, the right-hand side of Eq. is known from the time integration of its dynamic boundary condition (Eq. ).

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.

It is first assumed that the total velocity potential and the free surface displacement

(φ, r)

can be split into the incident component

(φ_{I}, r_{I})

and the scattered component

(φ_{S}, r_{S})

⁠:

{\begin{matrix} φ = φ_{I} + φ_{S} \\ r = r_{I} + r_{S} \end{matrix}

(7)

where

r_{I} = (x_{I}, η_{I})

and

r_{S} = (x_{S}, η_{S})

⁠. Then, a small parameter

ε

is introduced to denote the smallness of the quantities associated with the scattered waves, i.e.,

φ_{S} / φ_{I} \in O (ε)

and

r_{i}^{S} / r_{i}^{I} \in O (ε)

⁠,

i = 1, 2

⁠, where

r_{i}^{I}

and

r_{i}^{S}

are components of the position vectors

r_{I}

and

r_{S}

⁠. With inserting Eq. into Eq. and taking the Taylor expansion of the FSBC at the position of the incident wave and dropping small terms of higher order than

O (ε)

⁠, the free surface conditions expressed by the WS are

\begin{aligned} \frac{δ η_{S}}{δ t} & = \frac{\partial φ_{S}}{\partial z} - (\nabla φ_{I} - v) \cdot \nabla η_{S} - \nabla φ_{S} \cdot \nabla η_{I} - ν (r) η_{S} \\ - \nabla (\frac{\partial η_{I}}{\partial t} - \frac{\partial φ_{I}}{\partial z} + \nabla φ_{I} \cdot \nabla η_{I}) \cdot r_{S} at r = r_{I} \end{aligned}

(8)

\begin{aligned} \frac{δ φ_{S}}{δ t} & = - g η_{S} - (\nabla φ_{I} - v) \cdot \nabla φ_{S} - ν (r) φ_{S} \\ - \nabla (\frac{\partial φ_{I}}{\partial t} + g η_{I} + \frac{1}{2} \nabla φ_{I} \cdot φ_{I}) \cdot r_{S} at r = r_{I} \end{aligned}

(9)

The no-flux boundary conditions on the surface of the fixed body and the tank bottom are $\partial φ_{S} / \partial n = 0$ ⁠.

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.

The damping coefficient

ν (r)

in the boundary conditions (Eqs. and ) is only active in the wave-damping zone and is zero elsewhere. In this study,

ν (r)

is defined as in Ref. [23]:

ν (r) = {\begin{matrix} 0, & r \leq r_{0} \\ ν_{max} (- 2 {\tilde{r}}^{3} + 3 {\tilde{r}}^{2}), & r_{0} < r \leq r_{1} (downstream) \\ ν_{max} (1 + 2 {\tilde{r}}^{3} - 3 {\tilde{r}}^{2}), & r_{0} < r \leq r_{1} (upstream) \end{matrix}

(10)

where coefficient

ν_{max}

is empirically taken as

ν_{max} = 2.0

in this study;

r_{0}

and

r_{1}

are the start and end position of the damping zone.

\tilde{r}

is an auxiliary normalized coordinate defined as

\tilde{r} = \frac{r - r_{0}}{L_{damp}}

(11)

where

L_{damp} = x_{1} - x_{0}

is the length of the wave-damping zone.

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 $i = 1, \dots, 8$ ⁠) on its boundary and one (with indices $i = 9$ ⁠) located at the center of the cell, as illustrated in Fig. 2.

Fig. 2

View large Download slide

Sketch of global overlapping quadrilateral cells and a standalone cell with its local indices and coordinate system

A local

O ξ ζ

coordinate system is introduced in each cell with its origin located at the center node. Then, the velocity potential

φ (ξ, ζ)

within each cell is represented by a weighted linear superposition of the first eight harmonic polynomials, which satisfy the Laplace equation everywhere in space. Hence, the interpolated velocity potential in each cell automatically satisfies the Laplace equation. In two dimensions, those harmonic polynomials are given by the real and imaginary parts of the complex polynomial

(ξ + i ζ)^{n} = u_{n} (ξ, ζ) + i v_{n} (ξ, ζ)

(12)

where n is an integer number indicating the order of the polynomial and

i = \sqrt{- 1}

is the imaginary unit. Table 1 summarizes the harmonic polynomial for

n \leq 4

⁠.

Table 1

The harmonic polynomials up to fourth order

$z^{n}$	$u_{n} (ξ, ζ)$	$v_{n} (ξ, ζ)$
$(u_{n} + v_{n})^{0}$	1	0
$(u_{n} + v_{n})^{1}$	$ξ$	$ζ$
$(u_{n} + v_{n})^{2}$	$ξ^{2} - ζ^{2}$	$2 ξ ζ$
$(u_{n} + v_{n})^{3}$	$ξ^{3} - 3 ξ ζ^{2}$	$3 ξ^{2} ζ - ζ^{3}$
$(u_{n} + v_{n})^{4}$	$ξ^{4} - 6 ξ^{2} ζ^{2} + ζ^{4}$	$4 ξ^{3} ζ - 4 ξ ζ^{3}$

The local interpolation of velocity potential within a cell can be given as the following form, which makes use of the first eight harmonic polynomials in Table 1:

φ (ξ, ζ) = \sum_{j = 1}^{8} b_{j} f_{j} (ξ, ζ)

(13)

where

φ

is the velocity potential. The function

f_{j} (j = 1, \dots, 8)

are the first eight nonzero harmonic polynomials listed in Table 1, whereas

b_{j} (j = 1, \dots, 8)

are unknown superposition coefficients to be determined. Note that the term

ξ^{4} - 6 ξ^{2} ζ^{2} + ζ^{4}

is exclusively included between the two fourth-order harmonic polynomials in order to reduce the free surface wave dispersion errors in time domain analyses [26]. In other words, the fourth-order harmonic polynomial term

(ξ^{3} ζ - ξ ζ^{3})

is neglected.

Inserting the local coordinates

(ξ_{i}, ζ_{i})

(i = 1, 2, \dots, 8)

of the border nodes of a given cell into Eq. yields

φ_{i} = \sum_{j = 1}^{8} b_{j} f_{j} (ξ_{i}, ζ_{i})

(14)

Solving this linear equation system and relabeling the suffices $i$ and $j$ ⁠, we obtain the superposition coefficient $b_{i}$ as follows:

b_{i} = \sum_{j = 1}^{8} c_{i, j} φ_{j} (i = 1, \dots, 8)

(15)

where

c_{i, j}

is the element of the inverse of the matrix whose elements are

f_{j} (ξ_{i}, ζ_{i})

⁠. Here,

i

and

j

mean the row index and column index of the matrix, respectively. Substituting these coefficients obtained by Eq. back into the interpolation expression yields the interpolation expression using harmonic polynomials:

φ (ξ, ζ) = \sum_{i = 1}^{8} [\sum_{j = 1}^{8} c_{j, i} f_{j} (ξ, ζ)] φ_{i}

(16)

Note that the center node in a cell is also the boundary node of a neighboring cell of this cell, as in Fig. 2. We can then utilize this to establish the connectivity equations for the velocity potential in the fluid by plugging the local coordinates of the center node (with local index

i = 9

⁠) and the velocity potential

φ_{i}

(i = 1, \dots, 8)

on the cell boundary nodes into Eq. which yields

φ_{9} = φ (ξ = 0, ζ = 0) = \sum_{i = 1}^{8} c_{1, i} φ_{i}

(17)

which is used to assemble the global linear algebraic equation system of this potential solver.

On the boundaries of the computational domain, the Dirichlet boundary conditions can be enforced through Eq. , if the boundary is immersed in the grid system, or, directly specified as the known velocity potential on corresponding nodes, if these nodes themselves form the boundary. The Neumann boundaries can readily be enforced by analytically taking the normal derivative of the same equation, i.e.,

\frac{\partial φ}{\partial n} (ξ, ζ) = \sum_{i = 1}^{8} [\sum_{j = 1}^{8} c_{j, i} \nabla f_{j} (ξ, ζ) \cdot n (ξ, ζ)] φ_{i}

(18)

where

n

is the unit normal vector of a Neumann boundary node.

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 $(ξ_{m}, ζ_{m})$ and the velocity potential at this location $φ (ξ_{m}, ζ_{m})$ 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 $v = (0, \partial η / \partial t)$ ⁠, see Eq. (2).

Fig. 3

View large Download slide

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).

Fig. 4

View large Download slide

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 $n$ 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.

Since this work mainly considers the submerged structures under waves, the modeling procedure for waterlines on surface-piercing structures is therefore not presented, we refer readers the original work [9,22] for more details.

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.

Fig. 5

View large Download slide

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.

Special fluids are identified to enforce the connectivity conditions between two neighboring layers of the Cartesian grid system. Those fluid nodes are labeled as type-1 nodes (triangles in Fig. 5) and as type-2 nodes (circles in Fig. 5), respectively. The type-1 nodes are at the center node of parent cells in parent levels. For this type of node, the connectivity equation is enforced, with the coefficients

c_{1, i} (1, \dots, 8)

from the parent cells. The type-2 nodes lie at the horizontal or vertical middle lines of the parent cells. For this type of nodes, with local coordinates

(ξ, ζ)

in parent cells, the velocity potential

φ_{type2}

can be interpolated by

φ_{type 2} = \sum_{i = 1}^{8} [\sum_{j = i}^{8} c_{j, i} f_{j} (ξ, ζ)] φ_{i}

(19)

where

c_{j, i} (i, j = 1, \dots, 8)

are from the parent cells. It is noted that the left-hand side value of this equation is taken as an unknown to be solved in the global linear algebraic equation system, which is different from the Dirichlet boundary conditions using Eq. , where the right-hand side value is treated as a known value.

2.4 The Lagrangian Acceleration Potential Method.

Once the governing equations of the nonlinear wave-body interaction problems are solved, the forces and moments acting on the body due to the dynamic pressure of the fluid can be calculated by integrating the pressure over the body’s wet surface

Γ_{S}

as follows:

F_{hydro} = - \int_{Γ_{S}} P n d s

(20)

M_{hydro} = - \int_{Γ_{S}} P r \times n d x

(21)

where the pressure

P

can be obtained from the Bernoulli’s equation

P = - ρ (φ_{t} + \frac{1}{2} {| \nabla φ |}^{2} + g z)

(22)

where

φ_{t}

is the time derivative of the velocity potential

φ

⁠.

In our work, $φ_{t}$ is obtained using a Lagrangian acceleration potential $Ψ$ following Ref. [23]. For a fixed body, $Ψ = φ_{t}$ ⁠. To obtain $φ_{t}$ ⁠, we first define a BVP for $Ψ$ ⁠. On the free surface, the Dirichlet boundary condition $Ψ = - 1 / 2 {| \nabla φ |}^{2} - g η$ ⁠. On the body surface, the impermeable condition is $\partial Ψ / \partial n = 0$ ⁠. Once the above BVP defined for $Ψ$ is solved, the time derivative of the velocity potential $φ_{t}$ 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 wave-breaking suppression model incorporated into the original FN IB-AHPC solver is mainly composed of two steps. First, the wave-breaking onset is detected by a threshold value, namely the breaking criterion, which is a dynamic one related to the Lagrangian downward vertical acceleration of the free surface wave particle, as in

a_{downward} > γ g

(23)

where

γ

is a specified factor. This criterion instructs when and where the breaking suppression should be applied to the surface wave profile; lower values of

γ

indicate an earlier commencement of the smoothing procedure. For a monochromatic wave, the value of the factor is generally assumed as

γ = 0.5

⁠; for real seas, it is assumed that

γ = 1

⁠. Some field observations further indicated that the value of

γ

could be even lower than 0.4 [32]. The downward acceleration

a_{downward}

can be calculated as

a_{downward} = \frac{d w}{d t} = \frac{\partial w}{\partial t} + \frac{\partial w}{\partial x} \frac{d x}{d t} + \frac{\partial w}{\partial z} \frac{d z}{d t}

(24)

where

w

is the vertical velocity of a free surface fluid particle and

\partial w / \partial t

is calculated using a three-point backward difference method. In our work, two values of the factor

γ = 1.0

and

γ = 0.425

are studied.

Second, after the breaking onset is detected at a free surface marker with horizontal coordinate

x_{i}

⁠, a local three-point filter is applied in a 10-point region to smooth the surface profile in each time-step. The filter is with coefficients

m_{k} = (0.25, 0.5, 0.25)

⁠,

k = 0, 1, 2

⁠. Within the smoothing region, each neighbor marker with coordinate

x_{j}

⁠, where

j \in (j = i - 5, i + 5)

⁠, need to be smoothed. The smoothed data

η^{*} (x_{j})

are expressed as a weighted average of the original data

η (x_{j})

as

η^{*} (x_{j}) = \sum_{n = j - 1}^{j + 1} m_{(n - j + 1)} f (x_{n})

(25)

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 $O x z$ 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 $h = 0.46 m$ ⁠. 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 $R = 100 mm$ ⁠. The distance between the uppermost point of the cylinder and the mean free surface, i.e., the submergence $Z_{c}$ ⁠, is either 50 mm or 100 mm, depending on the case.

Fig. 6

View large Download slide

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 $a$ varying between $2$ and $28 mm$ ⁠. There are two wave frequencies selected as either $ω = 2 π \times 1.05 Hz$ or $2 π \times 1.22 Hz$ ⁠, with corresponding wavelengths $L_{w}$ being $1.38$ and $1.05$ ⁠, and the maximum wave steepness $k a$ for the two waves with different frequencies are approximately $0.13$ and $0.17$ ⁠, 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 $v = (0, \partial η / \partial t)$ ⁠, 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.

These waves are generated in the NWT by means of a wave-making zone, where a widely used and straightforward relaxation technique [34] is performed after each time integration of the free surface condition (Eqs. and ). After relaxation, the free surface elevation and velocity potential are given as

\begin{matrix} η_{new} (x, t) = γ_{r} (x) η_{ref} (x, t) + [1 - γ_{r} (x)] η_{num} (x, t) \\ φ_{new} (x, t) = γ_{r} (x) φ_{ref} (x, t) [1 - γ_{r} (x)] φ_{num} (x, t) \end{matrix}

(26)

where the subscripts {•}_num and {•}_ref denote the numerical results obtained from time integration and the reference solution from the analytical stream function theory, respectively. The relaxation functions

γ_{r} (x)

used in this work are proposed by Jacobsen et al. [34] as

γ_{r} (x) = {\begin{matrix} \frac{\exp ((1 - x_{r})^{3.5}) - 1}{\exp (1) - 1}, & x - x_{0} \leq L_{WMZ} \\ 0, & x - x_{0} \geq L_{WMZ} \end{matrix}

(27)

where

x_{r} = (x - x_{0}) / L_{WMZ}

is the local coordinate within the wave-making zone; the parameter

x_{0}

is the horizontal coordinate of the left boundary. In this work, a wave-making zone is

L_{WMZ} = 2 L_{w}

and a wave-damping zone is

L_{WDZ} = L_{w}

in length, where the

L_{w}

is the incident wavelength. The wave-damping zone is described in Eqs. and .

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 $ω = 2 π \times$ $1.05 Hz$ ⁠, wave amplitude $a = 10 mm$ ⁠, and cylinder submergence $Z_{c} = 100 mm$ were considered. Figure 7(a) shows different snapshots of free surface elevation along the numerical wave tank taken at an instant time $t = 19 T$ ⁠, which were obtained using a coarse mesh (20 elements per cylinder diameter $D$ ⁠, i.e., $Δ x = D / 20$ ⁠, in the deepest quadtree layer), an intermediate mesh (⁠ $Δ x = D / 30$ ⁠), a fine mesh (⁠ $Δ x = D / 40$ ⁠), respectively. A small time-step of (⁠ $Δ t = T / 80$ ⁠) 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.

Fig. 7

Comparison of different space-step Δx results, with incident wave amplitude a=10 mm and wave frequency ω=2π×1.05 Hz. The submergence of the cylinder is Zc=100 mm. The time-step Δt=T/80. (a) Snapshot of free surface elevation along the NWT at t=19T and (b) time history of vertical wave force on the cylinder.

View large Download slide

Comparison of different space-step $Δ x$ results, with incident wave amplitude $a = 10$ mm and wave frequency $ω = 2 π \times 1.05$ Hz. The submergence of the cylinder is $Z_{c} = 100$ mm. The time-step $Δ t = T / 80$ ⁠. (a) Snapshot of free surface elevation along the NWT at $t = 19 T$ and (b) time history of vertical wave force on the cylinder.

Likewise, tests were conducted using the intermediate spatial mesh and different time-steps (⁠ $Δ t = T / 30, T / 40, T / 80$ ⁠). Figure 8(a) shows the obtained time histories of free surface elevation recorded at a probe location $x = 0.75 L_{w}$ ⁠, 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 $Δ t = T / 40$ is used.

Fig. 8

Comparison of different time-step Δt results, with incident wave amplitude a=10 mm and wave frequency ω=2π×1.05 Hz. The submergence of the cylinder is Zc=100 mm. The space-step Δx=D/30. (a) Time history of free surface elevation recorded at x=0.75Lw and (b) time history of vertical wave force on the cylinder.

View large Download slide

Comparison of different time-step $Δ t$ results, with incident wave amplitude $a = 10$ mm and wave frequency $ω = 2 π \times 1.05$ Hz. The submergence of the cylinder is $Z_{c} = 100$ mm. The space-step $Δ x = D / 30$ ⁠. (a) Time history of free surface elevation recorded at $x = 0.75 L_{w}$ and (b) time history of vertical wave force on the cylinder.

3.3 Harmonic Wave Amplitudes Around the Cylinder.

Grue [16] conducted this experiment in a wave channel and particularly analyzed the first and second harmonic free wave amplitudes at the lee side of the circular cylinder against incoming wave amplitude. At the lee side of the cylinder, the waves are assumed given by

\begin{aligned} η_{+} (x, t) & = a_{+}^{(1)} \cos (k x - ω t + δ_{+}^{(1)}) + \sum_{n > 1} a_{l +}^{(n)} \cos n (k x - ω t + δ_{+}^{(1)}) \\ + \sum_{n > 1} a_{+}^{(n)} \cos (k_{n} x - n ω t + δ_{+}^{(n)}) \end{aligned}

(28)

where

a_{+}^{(n)}

⁠,

n = 1, 2, \dots,

are the

n

th harmonic transmitted free wave amplitudes, and

a_{l +}^{(n)}

(⁠

n = 2, 3, \dots

⁠) are the amplitudes of the locked

n

th harmonic wave components connected to the first harmonic transmitted wave;

δ_{+}^{(n)}

(⁠

n = 1, 2, \dots

⁠) are phase angles. For a given frequency

ω

and water depth

h

⁠, the wave number

k

is approximated by the linear dispersion relation

ω^{2} = g k \tanh k h

(29)

The wave numbers

k_{n}

(⁠

n = 2, 3, \dots

⁠) for free waves can be obtained by

(n ω)^{2} = g k_{n} \tanh k_{n} h, n = 2, 3, \dots

(30)

Following Grue [16], the free harmonic wave amplitudes were extracted by first taking the Fourier transform of the surface elevation

η (x, t)

obtained at the cylinder’s lee side as

{\hat{η}}^{(n)} (x) = \frac{ω}{2 π} \int_{0}^{2 π / ω} η (x, t) \exp (- i n ω t) d t, n = 1, 2, \dots

(31)

Then, the first harmonic transmitted wave and the higher harmonic free wave amplitudes at the lee side can be obtained by

a_{+}^{(1)} = \frac{1}{| \sin (k Δ x) |} | {\hat{η}}^{(1)} (x_{1}) - {\hat{η}}^{(1)} (x_{1} + Δ x) \exp (- i k Δ x) |

(32)

a_{+}^{(n)} = \frac{| {\hat{η}}^{(n)} (x_{1}) - {\hat{η}}^{(n)} (x_{1} + Δ x) \exp (i n k Δ x) |}{| \sin ((1 / 2) (k_{n} - n k) Δ x) |}, n = 2, 3, \dots

(33)

where

Δ x

is the distance between two wave gauges selected, which record the surface elevation at the cylinder’s lee side.

In Grue’s work, it was suggested that the free harmonic amplitudes obtained with this method are almost independent of the values of $Δ x$ and $x_{1}$ ⁠. 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 $x_{1}$ 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 $x_{1}$ ⁠, 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 $0.5 L_{W}$ 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 $Δ x$ was relatively small as observed also in our work. We set $Δ x \approx 1 / 5 L_{W}$ for the following studies, where $L_{W}$ is the incident wavelength.

Fig. 9

View large Download slide

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 $ω = 2 π \times 1.05$ Hz, and the submergence $Z_{c} = 100$ mm. The horizontal coordinate represents the distance from the cylinder center at its lee side, normalized by incident wavelength $L_{W}$ ⁠. (a) $a = 11$ mm and (b) $a = 19$ 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 $Z_{c} = R = 100$ mm, i.e., the distance between the free surface and its uppermost point. The incident wave amplitudes vary from $2$ to $28$ mm. Two wave frequencies are simulated, namely $ω = 2 π \times 1.05$ Hz and $ω = 2 π \times 1.22$ 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.

Fig. 10

View large Download slide

Normalized amplitudes of first and second harmonic free waves at the lee side of the circular with submergence $Z_{c} = 100$ mm. The FN and WS results are obtained using the presented IB-AHPC method, which are compared to Grue’s experiments [16] and the second-order theory of Friis [35]. (a) $ω = 2 π \times 1.05$ Hz and (b) $ω = 2 π \times 1.22$ Hz.

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 $a$ exceeds $20$ 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 $a$ ⁠. 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 $a$ 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 $Z_{c}$ ⁠, is reduced to $1 / 2 R = 50$ mm. The comparisons of the first and second harmonic free wave amplitudes $a_{+}^{(1)}$ and $a_{+}^{(2)}$ are shown in Fig. 11, for the incident wave frequency $ω = 2 π \times 1.05$ Hz, and in Fig. 12, for the incident wave frequency $ω = 2 π \times 1.22$ 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.

Fig. 11

As Fig. 10(a), with submergence Zc=50 mm, ω=2π×1.05 Hz. The arrows denote the spilling (S) and plunging (P) limits.

View large Download slide

As Fig. 10(a), with submergence $Z_{c} = 50$ mm, $ω = 2 π \times 1.05$ Hz. The arrows denote the spilling (S) and plunging (P) limits.

Fig. 12

As Fig. 10(b), with submergence Zc=50 mm, ω=2π×1.22 Hz. The arrows denote the spilling (S) and plunging (P) limits.

View large Download slide

As Fig. 10(b), with submergence $Z_{c} = 50$ mm, $ω = 2 π \times 1.22$ 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 $a_{+}^{(1)} / a$ 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, $a_{+}^{(2)} / a$ 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 $a_{+}^{(2)} / a$ 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 $a_{+}^{(2)} / a$ reaches its maximum value, which is approximately $0.4$ ⁠. Then, for larger incoming waves, $a_{+}^{(2)} / a$ 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 $a_{+}^{(1)} / a$ 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 $a_{+}^{(2)} / a$ 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 $a_{+}^{(2)} / a$ 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 $a_{+}^{(2)} / a$ 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 $γ = 0.425$ and $γ = 1.0$ ⁠. 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.

Fig. 13

As Fig. 11 and compared with results from using different breaking criterion γ. Submergence Zc=50 mm, ω=2π×1.05 Hz. The arrows denote spilling (S) and plunging (P) limits.

View large Download slide

As Fig. 11 and compared with results from using different breaking criterion $γ$ ⁠. Submergence $Z_{c} = 50$ mm, $ω = 2 π \times 1.05$ Hz. The arrows denote spilling (S) and plunging (P) limits.

Fig. 14

As Fig. 12 and compared with the results from using the wave-breaking suppression model. Submergence Zc=50 mm, ω=2π×1.22 Hz. The arrows denote the spilling (S) and plunging (P) limits.

View large Download slide

As Fig. 12 and compared with the results from using the wave-breaking suppression model. Submergence $Z_{c} = 50$ mm, $ω = 2 π \times 1.22$ 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 $D = 0.102$ m, fixed under the still water surface with submergence $Z_{c} = D$ in a wave channel. The wave channel was filled in water with a depth of $0.85$ m. The period of the incident wave is $1$ s, with wave height varying from $0.0018$ m to $0.0720$ 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 vertical force exerted on the cylinder can be represented by the Fourier series

F_{z} (t) = F_{z}^{(0 ω)} + \sum_{n = 1}^{N} F_{z}^{(n ω)} \sin (n ω t + δ^{(n ω)})

(34)

In our work, the harmonic force amplitudes are obtained using a least-square fitting technique. After a steady state has been achieved, at least three periods of the wave force time series are applied to the fitting. The obtained vertical forces are normalized by the KC number defined as

KC = \frac{2 π a}{D} e^{- k Z_{c}}

(35)

where

a

is the incident wave amplitude and

k

denotes the wave number.

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].

Fig. 15

View large Download slide

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 $γ_{0} = 0.425$ and $γ_{1} = 1.0$ ⁠. The incident wave frequency $ω = 2 π \times 1.05$ Hz, and the submergence of the cylinder $Z_{c} = 50$ mm.

Fig. 16

View large Download slide

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 $γ_{0} = 0.425$ and $γ_{1} = 1.0$ ⁠. The incident wave frequency $ω = 2 π \times 1.22$ Hz, and the submergence of the cylinder $Z_{c} = 50$ 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 $γ_{0} = 0.425$ and $γ_{1} = 1.0$ ⁠. 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

$g$ =: gravitational acceleration
$k$ =: wave number
$n$ =: normal vector
$r$ =: displacement vector
$v$ =: prescribed velocity
$D$ =: diameter of the cylinder
$R$ =: radius of the cylinder
$T$ =: incident wave period
$F$ =: body forces due to dynamic pressure of the fluid
$M$ =: body moments due to dynamic pressure of the fluid
$a_{downward}$ =: Lagrangian downward vertical acceleration
$L_{w}$ =: wavelength
$L_{WMZ}$ =: length of the wave-making zone
$L_{WDZ}$ =: length of the wave-damping zone
$V_{N}$ =: normal components of the rigid-body velocity
$Z_{c}$ =: submergence of the cylinder
$a_{+}^{(n)}$ =: harmonic free wave amplitudes
$x, z$ =: global coordinates
$ξ, ζ$ =: local coordinates in the harmonic polynomial cell

Greek Symbols

$γ$ =: specified factor related to breaking criterion
$η$ =: free surface elevation
$ν$ =: damping coefficient
$φ$ =: velocity potential
$φ_{t}$ =: time derivative of the velocity potential
$Ψ$ =: Lagrangian acceleration potential
$ω$ =: wave frequency

Dimensionless Groups

KC =: Keulegan–Carpenter number, $(2 π a / D) e^{- k Z_{c}}$

Abbreviations

FN =: fully nonlinear
FNPF =: fully nonlinear potential flow
HPC =: harmonic polynomial method
IB-AHPC =: immersed-boundary adaptive harmonic polynomial cell
NWT =: numerical wave tank
WS =: weak scatterer

References

1.

Babanin

,

A.

,

2011

,

Breaking and Dissipation of Ocean Surface Waves

, 1st ed.,

Cambridge University Press

,

Cambridge

.

Google Scholar

Crossref

2.

Guignard

,

S.

,

Marcer

,

R.

,

Rey

,

V.

,

Kharif

,

C.

, and

Fraunié

,

P.

,

2001

, “

Solitary Wave Breaking on Sloping Beaches: 2-D Two Phase Flow Numerical Simulation by SL-VOF Method

,”

Eur. J. Mech. B/Fluids

,

20

(

1

), pp.

57

–

74

.

Google Scholar

Crossref

3.

Abadie

,

S.

,

Morichon

,

D.

,

Grilli

,

S.

, and

Glockner

,

S.

,

2010

, “

Numerical Simulation of Waves Generated by Landslides Using a Multiple-Fluid Navier–Stokes Model

,”

Coastal Eng.

,

57

(

9

), pp.

779

–

794

.

Google Scholar

Crossref

4.

Lee

,

C.-H.

,

1995

,

WAMIT Theory Manual

,

Massachusetts Institute of Technology, Department of Ocean Engineering

,

Cambridge

.

5.

Ansys

,

A.

,

2013

,

AQWA Theory Manual.

,

ANSYS

,

Canonsburg, PA

.

6.

Bingham

,

H. B.

, and

Zhang

,

H.

,

2007

, “

On the Accuracy of Finite-Difference Solutions for Nonlinear Water Waves

,”

J. Eng. Math.

,

58

(

1

), pp.

211

–

228

.

Google Scholar

Crossref

7.

Yates

,

M. L.

, and

Benoit

,

M.

,

2015

, “

Accuracy and Efficiency of Two Numerical Methods of Solving the Potential Flow Problem for Highly Nonlinear and Dispersive Water Waves

,”

Int. J. Numer. Methods Fluids

,

77

(

10

), pp.

616

–

640

.

Google Scholar

Crossref

8.

Ducrozet

,

G.

,

Bonnefoy

,

F.

, and

Perignon

,

Y.

,

2017

, “

Applicability and Limitations of Highly Non-linear Potential Flow Solvers in the Context of Water Waves

,”

Ocean Eng.

,

142

, pp.

233

–

244

.

Google Scholar

Crossref

9.

Tong

,

C.

,

Shao

,

Y.

,

Bingham

,

H. B.

, and

Hanssen

,

F. W.

,

2021

, “

An Adaptive Harmonic Polynomial Cell Method With Immersed Boundaries: Accuracy, Stability, and Applications

,”

Int. J. Numer. Methods Eng.

,

122

(

12

), pp.

2945

–

2980

.

Google Scholar

Crossref

10.

Wu

,

C. H.

, and

Nepf

,

H. M.

,

2002

, “

Breaking Criteria and Energy Losses for Three-Dimensional Wave Breaking

,”

J. Geophys. Res. Oceans

,

107

(

C10

), pp.

1

–

18

.

Google Scholar

Crossref

11.

Barthelemy

,

X.

,

Banner

,

M. L.

,

Peirson

,

W. L.

,

Fedele

,

F.

,

Allis

,

M.

, and

Dias

,

F.

,

2018

, “

On a Unified Breaking Onset Threshold for Gravity Waves in Deep and Intermediate Depth Water

,”

J. Fluid Mech.

,

841

, pp.

463

–

488

.

Google Scholar

Crossref

12.

Derakhti

,

M.

,

Kirby

,

J. T.

,

Banner

,

M. L.

,

Grilli

,

S. T.

, and

Thomson

,

J.

,

2020

, “

A Unified Breaking Onset Criterion for Surface Gravity Water Waves in Arbitrary Depth

,”

J. Geophys. Res. Oceans

,

125

(

7

), p.

e2019JC015886

.

Google Scholar

Crossref

13.

Mohanlal

,

S.

,

Harris

,

J. C.

,

Yates

,

M. L.

, and

Grilli

,

S. T.

,

2023

, “

Unified Depth-Limited Wave Breaking Detection and Dissipation in Fully Nonlinear Potential Flow Models

,”

Coastal Eng.

,

183

, p.

104316

.

Google Scholar

Crossref

14.

Schløer

,

S.

,

2013

, “

Fatigue and Extreme Wave Loads on Bottom Fixed Offshore Wind Turbines. Effects From Fully Nonlinear Wave Forcing on the Structural Dynamics

,” Ph.D. thesis,

Technical University of Denmark

,

Lyngby

.

15.

Engsig-Karup

,

A.

,

Bingham

,

H.

, and

Lindberg

,

O.

,

2009

, “

An Efficient Flexible-Order Model for 3D Nonlinear Water Waves

,”

J. Comput. Phys.

,

228

(

6

), pp.

2100

–

2118

.

Google Scholar

Crossref

16.

Grue

,

J.

,

1992

, “

Nonlinear Water Waves at a Submerged Obstacle or Bottom Topography

,”

J. Fluid Mech.

,

244

(

1

), p.

455

.

Google Scholar

Crossref

17.

Keulegan

,

G.

, and

Carpenter

,

L.

,

1958

, “

Forces on Cylinders and Plates in an Oscillating Fluid

,”

J. Res. Natl. Bur. Stand.

,

60

(

5

), p.

423

.

Google Scholar

Crossref

18.

Chaplin

,

J. R.

,

1984

, “

Nonlinear Forces on a Horizontal Cylinder Beneath Waves

,”

J. Fluid Mech.

,

147

(

1

), p.

449

.

Google Scholar

Crossref

19.

Pawlowski

,

J.

,

1994

, “

A Nonlinear Theory of Ship Motion in Waves

,”

19th Symposium on Naval Hydrodynamics

,

Seoul, South Korea

,

Aug. 23, 1992

, pp.

33

–

58

.

20.

Letournel

,

L.

,

Chauvigné

,

C.

,

Gelly

,

B.

,

Babarit

,

A.

,

Ducrozet

,

G.

, and

Ferrant

,

P.

,

2018

, “

Weakly Nonlinear Modeling of Submerged Wave Energy Converters

,”

Appl. Ocean Res.

,

75

, pp.

201

–

222

.

Google Scholar

Crossref

21.

Zhang

,

Y.

, and

Teng

,

B.

,

2021

, “

A Nonlinear Potential Flow Model for Higher-Harmonic Wave Loads and Ringing Response of a Monopile

,”

Ocean Eng.

,

222

, p.

108574

.

Google Scholar

Crossref

22.

Tong

,

C.

,

Shao

,

Y.

,

Bingham

,

H. B.

, and

Hanssen

,

F.-C. W.

,

2022

, “

A Generalized Weak-Scatterer Approximation for Nonlinear Wave–Structure Interaction in Marine Hydrodynamics

,”

Marine Struct.

,

86

, p.

103292

.

Google Scholar

Crossref

23.

Greco

,

M.

,

2001

, “

A Two-Dimensional Study of Green-Water Loading

,” Ph.D.thesis,

Norwegian University of Science and Technology

,

Trondheim, Norway

.

24.

Berland

,

J.

,

Bogey

,

C.

,

Marsden

,

O.

, and

Bailly

,

C.

,

2007

, “

High-Order, Low Dispersive and Low Dissipative Explicit Schemes for Multiple-Scale and Boundary Problems

,”

J. Comput. Phys.

,

224

(

2

), pp.

637

–

662

.

Google Scholar

Crossref

25.

Savitzky

,

A.

, and

Golay

,

M. J. E.

,

1964

, “

Smoothing and Differentiation of Data by Simplified Least Squares Procedures

,”

Anal. Chem.

,

36

(

8

), pp.

1627

–

1639

.

Google Scholar

Crossref

26.

Shao

,

Y.-L.

, and

Faltinsen

,

O. M.

,

2012

, “

Towards Efficient Fully-Nonlinear Potential-Flow Solvers in Marine Hydrodynamics

,”

31st International Conference on Ocean, Offshore and Arctic Engineering

,

Rio de Janeiro, Brazil

,

July 1

, pp.

369

–

380

.

Google Scholar

Crossref

27.

Shao

,

Y.-L.

, and

Faltinsen

,

O. M.

,

2014

, “

A Harmonic Polynomial Cell (HPC) Method for 3D Laplace Equation With Application in Marine Hydrodynamics

,”

J. Comput. Phys.

,

274

, pp.

312

–

332

.

Google Scholar

Crossref

28.

Hanssen

,

F.-C. W.

,

Greco

,

M.

, and

Shao

,

Y.

,

2015

, “

The Harmonic Polynomial Cell Method for Moving Bodies Immersed in a Cartesian Background Grid

,”

34th International Conference on Ocean, Offshore and Arctic Engineering

,

St. John’s, Newfoundland, Canada

,

May 31

.

Google Scholar

Crossref

29.

Hanssen

,

F.-C.

,

Bardazzi

,

A.

,

Lugni

,

C.

, and

Greco

,

M.

,

2018

, “

Free-Surface Tracking in 2D With the Harmonic Polynomial Cell Method: Two Alternative Strategies

,”

Int. J. Numer. Methods Eng.

,

113

(

2

), pp.

311

–

351

.

Google Scholar

Crossref

30.

Hanssen

,

F.-C. W.

,

2019

, “

Non-Linear Wave-Body Interaction in Severe Waves

,” Ph.D. thesis,

Norwegian University of Science and Technology

,

Trondheim, Norway

.

31.

Ma

,

S.

,

Hanssen

,

F.-C. W.

,

Siddiqui

,

M. A.

,

Greco

,

M.

, and

Faltinsen

,

O. M.

,

2018

, “

Local and Global Properties of the Harmonic Polynomial Cell Method: In-Depth Analysis in Two Dimensions

,”

Int. J. Numer. Methods Eng.

,

113

(

4

), pp.

681

–

718

.

Google Scholar

Crossref

32.

Liu

,

P. C.

, and

Babanin

,

A. V.

,

2004

, “

Using Wavelet Spectrum Analysis to Resolve Breaking Events in the Wind Wave Time Series

,”

Ann. Geophys.

,

22

(

10

), pp.

3335

–

3345

.

Google Scholar

Crossref

33.

Fenton

,

J.

,

1988

, “

The Numerical Solution of Steady Water Wave Problems

,”

Comput. Geosci.

,

14

(

3

), pp.

357

–

368

.

Google Scholar

Crossref

34.

Jacobsen

,

N. G.

,

Fuhrman

,

D. R.

, and

Fredsøe

,

J.

,

2012

, “

A Wave Generation Toolbox for the Open-Source CFD Library: OpenFoam^®

,”

Numer. Methods Fluids

,

70

(

9

), pp.

1073

–

1088

.

Google Scholar

Crossref

35.

Friis

,

A.

,

Grue

,

J.

, and

Palm

,

E.

,

1991

, “Application of Fourier Transform to the Second Order 2D Wave Diffraction Problem,”

M. P. Tulin’s Festschrift: Mathematical Approaches in Hydrodynamics

,

T.

Miloh

, ed.,

SIAM

,

Philadelphia, PA

.

Fully Nonlinear Interaction of Water Waves and a Submerged Cylinder With Wave Breaking Detection and Suppression

Abstract

1 Introduction

2 Wave–Structure Interaction Model

2.1 Fully Nonlinear Potential-Flow Formulation.

2.2 Weak-Scatterer Approximation.

2.3 The Two-Dimensional Immersed-Boundary Adaptive Harmonic Polynomial Cell Method.

2.3.1 The Two-Dimensional Harmonic Polynomial Cell Method.

2.3.2 The Immersed-Boundary Method.

2.3.3 Adaptive Cartesian Grid System.

2.4 The Lagrangian Acceleration Potential Method.

2.5 Wave-Breaking Suppression Model.

3 A Fixed Submerged Circular Cylinder in Regular Waves

3.1 Numerical Setup.

3.2 Numerical Convergence Study.

3.3 Harmonic Wave Amplitudes Around the Cylinder.

3.4 Harmonic Wave Forces on the Cylinder.

4 Conclusion

Acknowledgment

Conflict of Interest

Data Availability Statement

Nomenclature

Greek Symbols

Dimensionless Groups

Abbreviations

References

Get Email Alerts

Cited By

ASME Journals

ASME Conference Proceedings

ASME eBooks

Resources

Opportunities

Fully Nonlinear Interaction of Water Waves and a Submerged Cylinder With Wave Breaking Detection and Suppression

Abstract

1 Introduction

2 Wave–Structure Interaction Model

2.1 Fully Nonlinear Potential-Flow Formulation.

2.2 Weak-Scatterer Approximation.

2.3 The Two-Dimensional Immersed-Boundary Adaptive Harmonic Polynomial Cell Method.

2.3.1 The Two-Dimensional Harmonic Polynomial Cell Method.

2.3.2 The Immersed-Boundary Method.

2.3.3 Adaptive Cartesian Grid System.

2.4 The Lagrangian Acceleration Potential Method.

2.5 Wave-Breaking Suppression Model.

3 A Fixed Submerged Circular Cylinder in Regular Waves

3.1 Numerical Setup.

3.2 Numerical Convergence Study.

3.3 Harmonic Wave Amplitudes Around the Cylinder.

3.4 Harmonic Wave Forces on the Cylinder.

4 Conclusion

Acknowledgment

Conflict of Interest

Data Availability Statement

Nomenclature

Greek Symbols

Dimensionless Groups

Abbreviations

References

Get Email Alerts

Cited By

Related Articles

Related Proceedings Papers

Related Chapters

ASME Journals

ASME Conference Proceedings

ASME eBooks

Resources

Opportunities

This Feature Is Available To Subscribers Only