Comparison of wave–structure interaction dynamics of a submerged cylindrical point absorber with three degrees of freedom using potential flow and computational fluid dynamics models

Power production using wave energy gained momentum in the 1970s during the oil crisis. This field is regaining a renewed interest in the marine hydrokinetic research community that is aiming to reduce the current carbon footprint of power production. Despite the abundantly available wave power in the oceans and seas worldwide¹ and research efforts dating back since the 1970s,² no commercial-scale wave power production operations exist today. Consequently, various wave energy conversion (WEC) concepts have been proposed and implemented, yet no single device architecture has been recognized as the ultimate solution.

Point absorber (PA) is a type of WEC system, for which the device characteristic dimensions are relatively small compared to the wavelength of the site.^3,4 Depending on the wave energy extraction mechanism and the power take-off (PTO) system employed, PAs can be further categorized into different subtypes. To name a few, Inertial Sea Wave Energy Converter (ISWEC) developed by the Polytechnic University of Turin is a floating point absorber (FPA) that converts the pitching motion of the hull to an electrical output using a gyroscopic PTO system.^5–7 PowerBuoy is a two-body FPA developed by Ocean Power Technologies that uses the heave mode to extract energy from the waves.^8–10 Although FPAs have the advantage of receiving a dense concentration of wave energy from the ocean or sea surface, they are also prone to extreme waves and other severe weather conditions that can limit their operability and long-term survivability. Fully submerged point absorbers (SPAs) have been designed to overcome these issues. CETO is a cylindrical shaped SPA developed by Carnegie Wave Energy that is able to absorb wave energy using multiple degrees of freedom (DOF).^11,12 An added advantage of SPAs is their zero visual impact on the ocean or sea shorelines.¹³ Currently, efforts are underway that are testing point absorber devices at various locations around the world, including but not limited to, the Pacific Ocean,¹⁴ the Atlantic Ocean,¹⁵ and the Mediterranean Sea.^16,17 Some of us are also directly involved in testing and improving WEC devices at various sea locations.^5–7,17

Numerical models based on frequency- or time-domain methods are commonly used to study the performance of point absorbers.^18–20 The hydrodynamic loads in these methods are calculated using the boundary element method (BEM) approach, which is based on a linear potential flow (LPF) formulation. The BEM approach to WEC modeling solves a radiation and a diffraction problem of the oscillating converter separately. The pressure solutions from the radiation and diffraction problems and the pressure field of the undisturbed incident wave are superimposed to obtain the net hydrodynamic load on the wetted surface of the converter. Frequency- or time-domain methods ignore the viscous phenomena and the nonlinear convective terms from the equations of motion. As a result, these methods cannot capture highly nonlinear phenomena such as wave-breaking and wave-overtopping. Moreover, these methods overpredict the dynamics and the wave absorption efficiency of the WEC systems.¹⁰ An improvement over LPF based models is the fully nonlinear potential flow (FNPF) formulation, which permits large-amplitude displacements of the WECs and modeling of nonlinear free-surface.²¹ Furthermore, FNPF models impose body boundary conditions based on the instantaneous location of the WEC in the computational domain, rather than assuming the free-surface and WECs at their equilibrium positions.

A considerable amount of accuracy in WEC modeling is achieved by solving the nonlinear incompressible Navier–Stokes (INS) equations of motion,^10,22–25 albeit at a higher computational cost compared to the BEM technique. Several approaches to fully resolved wave–structure interaction (WSI) modeling have been adopted in the literature. The two main categories are (i) the overset or the Chimera grid-based methods and (ii) the fictitious domain-based methods. The overset method employs an unstructured mesh for the solid structure and a background fluid grid, which is generally taken as block structured.^26–29 The fictitious domain (FD) approach to fluid–structure interaction (FSI) modeling is a Cartesian grid-based method in which the fluid equations are extended into the solid domain, and a common set of equations are solved for the two domains. Fictitious domain methods have been used to simulate FSI of porous structures,³⁰ elastic boundaries,³¹ and rigid bodies of complex shapes.^32,33 FD methods can be implemented in several ways, for example, by using the distributed Lagrange multiplier (DLM) technique^34,35 or employing the Brinkman penalization (BP) approach.^36–38 Ghasemi et al.²⁴ and Anbarsooz et al.²⁵ have studied the performance of a submerged cylindrical shaped PA with two degrees of freedom using the FD/DLM approach. The cylindrical buoy was constrained to move in two orthogonal directions (heave and surge) in these studies. One of the limitations of the FD/DLM method is that any external force acting on the immersed object (e.g., via tethered PTO system) needs to be expressed as a distributed body force density in the INS momentum equation. This is not a versatile approach, but it can work for simple scenarios.⁸⁸ However, applying external torque on the structure is not straightforward for the FD/DLM method because the velocity, and not vorticity, is generally solved for in the INS equations. In contrast, it is straightforward to include both external forces and torques on the immersed object using the FD/BP methodology. Moreover, the FD/BP method is a fully Eulerian approach to FSI. This makes the parallel implementation of the technique relatively easier compared to the FD/DLM method, which is typically implemented using two (Eulerian and Lagrangian) grids. Since we are interested in studying the dynamics of a three degrees of freedom (three-DOF) buoy under the action of external forces and torques, we employ the more versatile FD/BP approach in our computational fluid dynamic (CFD) model. To our knowledge, this works presents the first application of the FD/BP method for simulating WEC devices.

Using the FD/BP framework, we simulate the wave–structure interactions of a cylindrical buoy in one, two, and three degrees of freedom. The CFD solution is compared against the potential flow model. We find that the Cummins model overpredicts the heave and surge amplitude and does not capture the slow drift in the surge dynamics. Moreover, the potential flow model results in an insignificant pitch of an axisymmetric converter. Finally, we study the wave absorption efficiency of a three-DOF cylindrical buoy under varying PTO coefficients, mass density of the buoy, and incoming wave heights using the CFD method. The resolved simulations provide useful insights into an efficient design procedure for a simple WEC.

The rest of the paper is organized as follows. We first describe the potential flow formulation and the time-domain Cummins model in Sec. II. Next, we describe the continuous and discrete equations for the multiphase fluid–structure system in Secs. III A and III D, respectively. Validation cases for the FD/BP framework are presented in Sec. III F. The tests also highlight the solver stability in the presence of large density contrasts. Section IV compares the dynamics of the cylindrical buoy using potential flow and CFD models. The performance of a submerged cylindrical buoy using the CFD model under various scenarios is presented in Sec. V.

Using the potential flow model, the velocity potential Φ of an inviscid and incompressible fluid under the assumption of irrotational flow is obtained by solving the Laplace equation in the water domain,

using suitable kinematic and dynamic boundary conditions.^39,40 The fluid velocity is expressed as the gradient of the velocity potential, u = (u, $v$⁠) = ∇Φ. Once the solution to the Laplace equation is found, the fluid pressure p is obtained from the linearized Bernoulli equation ∂Φ/∂t + p/ρ_w + gy = 0 as

in which g is the acceleration due to gravity, ρ_w is the density of water, and y = 0 represents the undisturbed free-water surface. The hydrodynamic force on a submerged body is obtained by integrating pressure forces on the wetted surface of the body.

In the linear wave theory, the wave amplitude and the resulting body motions are assumed to be small compared to the wavelength of the incident wave. Under this assumption, the flow potential can be divided into three distinct parts,⁴¹ 

in which Φ_I is the undisturbed (assuming no body in the domain) wave potential of the incident wave, Φ_D is the diffraction potential of the incident wave about the stationary body, and Φ_R is the radiation potential due to an oscillatory motion of the body in still water. If the motion of an oscillating body such as a wave energy converter is not affected by nonlinearities in the system (like those arising from nonlinear power take-off units), frequency-domain models are generally used to obtain the solution of motion X due to a monochromatic harmonic wave excitation of angular frequency ω,

Here, M is the mass matrix of the buoy, A(ω) and B(ω) are the frequency-dependent added mass and damping matrix of the buoy, respectively, K_hydro is the linear hydrostatic stiffness arising from the buoyancy force for a floating buoy, and F_e(ω) is the wave excitation force (Froude–Krylov and diffraction). The displacement, velocity, and acceleration of the body are denoted by X, $Ẋ$⁠, and $Ẍ$⁠, respectively. The dimensions of the matrices and vectors depend on the degrees of freedom.⁴² The frequency-domain models are typically resolved using the boundary element method (BEM), and this approach has been widely adopted in commercial codes such as ANSYS AQWA⁴³ or WAMIT.⁴⁴ 

When nonlinear effects such as viscous forces or nonlinear PTO interactions with the buoy are considered, the linearity assumption of frequency-domain models is no longer valid. A common approach to overcome this limitation for many seakeeping applications is to use the time-domain model based on the Cummins equation,^45,46 in which the nonlinearities are modeled as time-varying coefficients of a system of ordinary differential equations.^47,48 The Cummins equation is a vector integro-differential equation, which involves a convolution of radiation impulse response function (RIRF) with the velocity of the body and reads as

in which A_∞ is the added mass at infinite frequency, given by

and K_r is the radiation impulse response function (RIRF). The radiation convolution function is also referred to as memory function because it represents a fluid memory effect due to the radiation forces emanated by the oscillating body in the past. In the Cummins equation, all nonlinear effects are lumped into the F_ext term, which represents external forces applied to the system. The nonlinear external forces could arise, for example, due to viscous drag or PTO/mooring forces. One of the computational challenges to the solution of the time-domain Cummins equation is the convolution integral involving the memory function. The time-varying RIRF can be evaluated by a number of numerical methods; see, for example, works of Yu and Falnes,⁴⁹ Jefferys,⁵⁰ Damaren,⁵¹ McCabe et al.,⁵² and Clément.⁵³ In this work, we follow the state-space representation approach of Perez and Fossen^54,55 to approximate the radiation convolution integral by

in which $Ar$⁠, $Br$⁠, and $Cr$ are the state-space matrices for carrying out the time-domain analysis. It is also possible to evaluate the RIRF in the frequency domain and subsequently transform it back to the time domain.^50,56,57

Figure 1 shows the schematic representation of a submerged point absorber tethered to the sea floor by the PTO system, as modeled in this work. In our model, the PTO provides both damping and stiffness loads on the submerged buoy. For the PA system considered here, the external forces arise from the nonlinear viscous drag F_drag and the PTO stiffness and damping loads, which are denoted by F_m and F_PTO, respectively. Since the point absorber is taken to be completely submerged under water, the hydrostatic stiffness arising from the buoyancy force is K_hydro = 0. Instead, the buoy experiences a permanent hydrostatic force in the upward direction. Accounting for all the external forces acting on the buoy, the Cummins equation for a submerged two-dimensional buoy of density ρ_s, diameter D, and volume V_buoy = πD²/4 reads as

in which k_PTO and b_PTO are the stiffness and damping constants of the PTO, Δl is the elongation of the PTO from a reference length, Δl₀ = |F_hydrostatic|/k_PTO is the permanent extension of the PTO to balance the hydrostatic force F_hydrostatic, and $l^$ is a unit vector along PTO in the current configuration. The nonlinear viscous drag on the two-dimensional cylindrical buoy is modeled following the resistive drag model of Ding et al.,¹⁹ 

in which C_x, C_y, and C_θ are the drag coefficients in the surge (⁠$x^$⁠), heave (⁠$y^$⁠), and pitch (⁠$z^$⁠) directions, respectively, and S_x = S_y = D is the planar cross-sectional area of the disk.

To solve the time-domain Cummins equation (8) for a cylindrical shaped submerged point absorber, we use ANSYS AQWA to extract the frequency-dependent quantities such as: (1) added mass matrix A(ω), (2) wave excitation force F_e(ω), and (3) radiation damping matrix B(ω). Further post-processing is done to transform the AQWA data to time-domain amenable quantities such as A_∞, $Ar$⁠, $Br$⁠, and $Cr$⁠. The approach described in Perez and Fossen^54,55 is implemented using custom MATLAB scripts for this purpose. Finally, the time-domain equations (8)–(11) are solved using an in-house Simulink-based code.

We remark that ANSYS AQWA provides frequency-dependent data for three-dimensional geometries. To obtain the two-dimensional data for a submerged disk, we simulate a three-dimensional cylinder of length L and diameter D and normalize the extracted quantities by L. In order to ensure that the extracted quantities are not affected by finite length truncation effects, we perform three BEM simulations in AQWA by taking L = 8D, 15D, and 30D. Figure 2 shows the time history of normalized surge and heave force on a three-dimensional cylinder with different lengths. Since the length-normalized forces are almost identical for the three chosen lengths, the finite length truncation effects are negligible for a cylinder of length 8D and beyond. Similar trends are also obtained for normalized added mass A(ω) and radiation damping B(ω) matrices in Sec. V A.

We use the fictitious domain Brinkman penalization (FD/BP) method to perform fully resolved wave–structure interaction simulations. The FD/BP method is a fully Eulerian approach to FSI modeling. Contrary to the body-conforming mesh techniques, the FD/BP method extends the fluid domain equations into the solid domain, and a common set of equations are written for the two domains. This modeling approach makes the fictitious domain methods computationally more efficient for moving bodies compared to the body-conforming mesh techniques.

The WSI framework is implemented within IBAMR,⁵⁸ which is an open-source C++ library providing support for immersed boundary methods with adaptive mesh refinement.⁵⁹ IBAMR relies on SAMRAI^60,61 for Cartesian grid management and the AMR framework and PETSc^62–64 for linear solver support.

We begin by first describing the continuous equations of motion for the multiphase FD/BP method and thereafter detail its spatiotemporal discretization.

We state the governing equations for a coupled multiphase fluid–structure system in a fixed region of space $Ω⊂Rd$⁠, for d = 2 spatial dimensions. A fixed Eulerian coordinate system x = (x, y) ∈ Ω is used to describe the momentum and incompressibility of fluid and solid domains. The spatial location of the immersed body Ω_b(t) ⊂ Ω is tracked using an indicator function χ(x, t), which is nonzero in the solid domain and zero in the fluid domain Ω_f(t) ⊂ Ω. The time-varying solid and fluid domains are non-overlapping, and their union occupies the entire domain Ω = Ω_f(t) ∪ Ω_b(t). We employ spatially and temporally varying density ρ(x, t) and viscosity μ(x, t) fields to describe the coupled three phase system. The equations of motion for the multiphase fluid–structure interaction system read as

Equation (15) is the conservative form of the momentum equation, whereas Eq. (16) describes the incompressibility of the system. The quantity u(x, t) expresses the velocity, p(x, t) represents the mechanical pressure, and f_c(x, t) represents the Eulerian constraint force density, which is nonzero only in the solid domain. The acceleration due to gravity is taken to be in the negative y direction g = (0, −g). In the FD/BP method, the rigidity-enforcing constraint force f_c(x, t) is defined as a penalization force that enforces a rigid body velocity u_b(x, t) in Ω_b(t). By treating the immersed structure as a porous region of vanishing permeability κ ≪ 1, the Brinkman penalized constraint force is formulated as

Section III E 3 describes the fluid–structure coupling algorithm and the rigid body velocity u_b(x, t) calculation.

We use scalar level set function ϕ(x, t) to identify liquid and gas regions, Ω_l ⊂ Ω and Ω_g ⊂ Ω, respectively, in the computational domain. The combined liquid and gas region denotes the fluid region described in Sec. III A, i.e., Ω_l ∪ Ω_g = Ω_f. The zero-contour of ϕ function implicitly defines the liquid–gas interface Γ(t) = Ω_l ∩ Ω_g. Similarly, the Eulerian indicator function χ(x, t) of the immersed body is expressed in terms of the level set function ψ(x, t),⁸⁹ whose zero-contour implicitly defines the surface of the body S_b(t) = ∂V_b(t); see Fig. 3. Without loss of generality, ϕ level set (signed distance) values are taken to be negative in the liquid phase and positive in the air phase. Similarly, ψ level set values are taken to be negative inside the solid body, whereas they are taken to be positive outside the solid region. Using the signed distance property of the level set functions, the material properties such as density and viscosity can be conveniently expressed as a function of ϕ(x, t) and ψ(x, t) fields,

As the simulation advances in time, the phase transport is described by the advection of level set fields by the local fluid velocity, which in the conservative form reads as

The signed distance property of ϕ and ψ is generally disrupted under linear advection, Eqs. (20) and (21). Let $ϕ̃n+1$ denote the flow level set function following an advective transport after time stepping through the interval $tn,tn+1$⁠. The flow level set is reinitialized to obtain a signed distance field ϕⁿ⁺¹ by computing a steady-state solution to the Hamilton–Jacobi equation,

which yields a solution to the Eikonal equation ∥∇ϕ∥= 1 at the end of each time step. We refer the readers to Ref. 65 for more details on the specific discretization of Eqs. (22) and (23), which employs second-order ENO finite differences combined with a subcell-fix method described by Min⁶⁶ and an immobile interface condition described by Son.⁶⁷ Both the subcell-fix method and the immobile interface condition have been shown to be effective in conserving mass of the flowing phases in the literature.

Since we consider a simple geometry (cylinder) in this work, the solid level set ψⁿ⁺¹ is analytically calculated and reinitialized by using the new location of the center of mass at tⁿ⁺¹. For more complex structures, computational geometry techniques can be employed to compute the signed distance function. (In our code, we use surface triangulation of complex geometries to compute the signed distance function.)

The continuous equations of motion [Eqs. (15) and (16)] are discretized on a staggered Cartesian grid, as shown in Fig. 3. The discretized domain is made up of N_x × N_y rectangular cells that cover the physical domain Ω. The mesh spacing in the two directions is denoted by Δx and Δy. Taking the lower left corner of the rectangular domain to be the origin (0, 0) of the coordinate system, each cell center of the grid has position $xi,j=(i+12)Δx,(j+12)Δy$ for i = 0, …, N_x − 1 and j = 0, …, N_y − 1. The physical location of the vertical cell face is half a grid space away from x_i,j in the x direction and is given by $xi−12,j=iΔx,(j+12)Δy$⁠. Similarly, $xi,j−12=(i+12)Δx,jΔy$ is the physical location of the horizontal cell face that is half a grid cell away from x_i,j in the y direction. The level set fields, pressure degrees of freedom, and the material properties are all approximated at cell centers and are denoted by $ϕi,jn≈ϕxi,j,tn$⁠, $ψi,jn≈ψxi,j,tn$⁠, $pi,jn≈pxi,j,tn$⁠, $ρi,jn≈ρxi,j,tn$⁠, and $μi,jn≈μxi,j,tn$⁠, respectively; some of these quantities are interpolated onto the required degrees of freedom as needed (see Ref. 65 for further details). Here, tⁿ denotes the time at time step n. The velocity degrees of freedom are defined on the cell faces and are denoted by $ui−12,jn≈uxi−12,j,tn$ and $vi,j−12n≈vxi,j−12,tn$⁠. Additional body forces on the right-hand side of the momentum equation are approximated on the cell faces of the staggered grid.

We use second-order finite differences to discretize all spatial derivative operators. The discretized version of the spatial operator is denoted with a h subscript, i.e., ∇ ≈ ∇_h. Further details on the spatial discretization can be obtained in prior studies.^65,68–70

Discretely, we use smoothed Heaviside functions to transition from one material phase to the other. The transition zone occurs within n_cells grid cells on either side of water–air interface Γ or fluid–solid interface S_b. Correspondingly, two numerical Heaviside functions are defined,

The number of transition cells across the Γ or S_b interface is assumed to be the same. This is not a strict requirement of our numerical method, but it holds true for all the WSI cases simulated in this work. A two-step process (see Fig. 4) is used to prescribe a given material property $I$ (such as ρ or μ) in the whole domain Ω:

• First, the material property in the “flowing” phase is set via the liquid–gas level set function,

• Next, $Iflow$ is corrected by accounting for the structural material property to obtain $Ii,jfull$ throughout the computational domain,

We employ a fixed-point iteration time stepping scheme with n_cycles = 2 cycles per time step to advance quantities of interest at time level tⁿ to time level tⁿ⁺¹ = tⁿ + Δt. If k superscript denotes the cycle number of the fixed-point iteration, then at the beginning of each time step, we set k = 0, u^n+1,0 = uⁿ, $pn+12,0=pn−12$⁠, ϕ^n+1,0 = ϕⁿ, and ψ^n+1,0 = ψⁿ. For n = 0 initial time level, all the physical quantities have a prescribed initial condition.

The level set functions are time-marched using an explicit advection scheme,

in which Q(·, ·) represents an explicit piecewise parabolic method (xsPPM7-limited) approximation to the linear advection terms on cell centers.^69,71

The spatiotemporal discretization of the conservative form of incompressible Navier–Stokes equations [(15) and (16)] reads as

in which the newest approximation to density $ρ̆n+1,k+1$ and the discretization of the convective term C^n+1,k are computed such that they satisfy consistent mass/momentum transport, which is required to maintain numerical stability for a high value of the air–water density ratio. We refer the reader to the work of Nangia et al.^65,58 for more details on the consistent mass/momentum transport scheme. Note that $Lμun+12,k+1=12Lμun+1,k+1+Lμun$ is a semi-implicit approximation to the viscous strain rate with $Lμn+1=∇h⋅μn+1∇hu+∇huTn+1$⁠. The newest approximation to viscosity μ^n+1,k+1 is obtained via the two-stage process described in Eqs. (26) and (27).

In Eq. (30), ρⁿ is the face-centered value of the density field ρ^full [refer to Eq. (27) for the definition of ρ^full], and the density field $ρ̆n+1,k+1$ is obtained by advecting ρⁿ discretely. The specific value of the density field ℘ used to compute the gravitational body force ℘g is explained next in the context of the FD/BP fluid–structure coupling algorithm.

The Brinkman penalization term that imposes the rigidity constraint in the solid region is treated implicitly in Eq. (30) and is expressed as

in which $χ̃=1−H̃body$⁠, $H̃body$ represents the regularized structure Heaviside function [Eq. (25)], and $κ∼O(10−8)$⁠; this permeability value has been found sufficiently small to enforce the rigidity constraint effectively in the prior studies.^72,58 In Eq. (32), u_b is the solid body velocity, which can be expressed as a sum of translational U_r and rotational W_r velocities,

The center of the mass point $XCOMn+1,k+1$ is updated using the rigid body translational velocity as

The rigid body velocity is computed by integrating Newton’s second law of motion,

in which M_b is the mass, I_b is the moment of the inertia matrix, $F$ and $M$ are the net hydrodynamic force and torque, respectively, M_bg is the net gravitational force, F_m and F_PTO are the PTO stiffness and damping force, respectively [see Eqs. (10) and (11)], and M_m and M_PTO are the corresponding PTO stiffness and damping torque about the center of mass point X_COM, respectively. The net hydrodynamic force $F$ and torque $M$ on the immersed body are computed by summing the contributions from pressure and viscous forces acting on the body,

The hydrodynamic traction in Eqs. (37) and (38) is evaluated on Cartesian grid faces that define a stair-step representation of the body on the Eulerian mesh.⁵⁸ In Eqs. (37) and (38), n_face denotes the unit normal to the face (pointing away from the solid and into the fluid), and ΔA_face is the area of the face.

Since the net gravitational force on the body is already included in Eq. (35), we exclude it from the body region in the INS momentum equation (30). Therefore, we set ℘g = ρ^flowg in Eq. (30), which also avoids spurious currents due to the large density variation near the fluid–solid interface.⁷³ More specifically, ρ^flow is obtained using Eq. (26) and not through Eq. (27).

To ensure that the fully resolved CFD model performs reliably for the submerged point absorber problem and remains stable under high density contrast between different phases, two validation cases are considered. For the first validation case, damped oscillations of a dense and a light cylinder in various damping regimes are considered, and the simulated results are compared against analytical solutions.⁷⁴ For the second validation case, we perform grid convergence tests for a submerged point absorber that oscillates in the heave direction under wave excitation loads. For validation of other aspects of FD/BP method and numerical wave tank (NWT) implementation, we refer readers to our prior works.^58,73

For all the cases considered in this work, water and air densities are taken to be 1025 kg/m³ and 1.2 kg/m³, respectively, and their respective viscosities are taken to be 10⁻³ Pa s and 1.8 × 10⁻⁵ Pa s. Surface tension effects are neglected as they do not affect the wave and converter dynamics at the scale of the problems considered in this work.

We begin by considering a simplified version of the point absorber problem. A planar disk of diameter D and mass density ρ_s (or mass M = ρ_sπD²/4) is attached to a spring of stiffness value k_spring and a mechanical damper of damping coefficient c_damper/c_critical = ζ. Here, $ccritical=2kspringM$ is the critical damping coefficient. The value of ζ determines the behavior of the system: ζ < 1 leads to under-damping, ζ > 1 leads to over-damping, and ζ = 1 results in a critically damped system.

The computational domain is taken to be Ω: [0, 5D] × [0, 10D], and the cylinder’s initial center of the mass point is released from (X₀, Y₀) = (2.5D, 8D). The rest length of the spring is taken as 6.5D, which gives the initial extension of the spring to be 1.5D. The spring and the damper are connected to the bottom of the domain, as shown in Fig. 5(a). The cylinder is surrounded by the air phase. The density ratio between the solid and fluid phases is denoted by m^*.

Next, we consider the damped oscillation of the cylinder at different density ratios m^* = 100, 2, and 0.8. We fix D = 0.2 m and k_spring = 500 N/m in these simulations. Four test cases corresponding to ζ = 0.1, 0.5, 1.0, and 1.5 are simulated with each density ratio.

We discretize the domain by a 100 × 200 uniform grid and use a constant time step size of Δt = 2.5 × 10⁻⁴ s to simulate the four damping ratio cases. Figure 5(b) compares the center of mass vertical position as a function of time with the analytical solutions in different damping regimes. The analytical solutions are derived by neglecting the fluid forces on the mass. With a large density difference between the solid and air phases, the hydrodynamic forces are small compared to the inertia of the system, and the fluid–structure interaction solution matches the analytical solution quite well, as observed in Fig. 5(b). Moreover, the results also show that the FD/BP methodology remains stable for large density ratios.

To ensure that the aforementioned spatiotemporal resolution produces a converged FSI solution, we perform a convergence study with coarse, medium, and fine grids for the ζ = 0.5 case. Table I tabulates the spatiotemporal resolutions used for these grids. Figure 6 shows the grid convergence plots for the cylinder displacement and the hydrodynamic force acting in the vertical direction. As observed in Fig. 6, a medium grid resolution of 100 × 200 is sufficient to resolve the FSI dynamics. Moreover, the spurious oscillations (as a function of time) in the hydrodynamic forces are significantly reduced at higher (medium and fine in this case) spatiotemporal resolutions, which is an expected trend for fictitious domain/immersed techniques.⁷⁵ Figure 7 shows the pressure and flow distribution around the cylinder at t = 0.05 s.

If the density of the solid is comparable to the surrounding fluid, the fluid forces are not negligible, and they affect the motion of the solid, particularly at low damping ratios. We simulate the previous case with low density ratios of m^* = 2 and 0.8. Figure 8 compares the cylinder displacement with the analytical solution. As observed in Figs. 8(a) and 8(b), the CFD solution deviates significantly from the analytical solution for ζ = 0.1 and 0.5 with m^* = 2 and for ζ = 0.1, 0.5, and 1 for m^* = 0.8, respectively. This deviation is expected physically for these cases, however.

The results in this section show that the fully resolved FD/BP model can accurately capture the rigid body dynamics of a system in which an external forcing is provided with mechanical springs and dampers. The method remains stable for a large density ratio between different phases, and it does not suffer from the added mass effects that become prominent when m^* ⪅ 1.^76–80 We also refer the readers to our prior work on the FD/BP method⁵⁸ that considers additional low density ratio cases and validates them against numerical and experimental results from the literature.

As a next validation case of our fully resolved WSI model, we consider the heave motion of a submerged cylindrical buoy induced by the incoming water waves. This case is close to the actual problem of interest and serves to provide an estimate of mesh resolution needed to resolve the WSI dynamics adequately. First, we describe the problem setup.

Figure 9 shows the schematic of the problem and the numerical wave tank (NWT) layout. Water waves of height $H$⁠, wavelength λ, and time period $T$ are generated at the left end of the domain using fifth-order Stokes wave theory.⁸¹ The waves propagate a distance of 8λ in the positive x direction before getting attenuated in the damping zone of length of 2λ located toward the right end of the domain. The submerged point absorber is placed midway at a distance of 4λ in the working zone of NWT. The depth of submergence of the point absorber is d_s, whereas the mean depth of water in the NWT is d. Regular water waves of $H=0.01$ m, λ = 1.216 m, $T=0.8838$ s, and d = 0.65 m are generated in the NWT. The simulated waves satisfy the dispersion relationship

and are in deep water regime d/λ > 0.5.

The diameter of the buoy with relative density ρ_s/ρ_w = 0.9 is taken to be D = 0.16 m. The PTO stiffness and damping coefficients are taken as k_PTO = 1995.2 N/m and b_PTO = 80.64 N s/m, respectively. The initial submergence depth of the disk is d_s = 0.25 m. Three grid resolutions corresponding to coarse, medium, and fine grids are used to simulate the heave dynamics of the buoy. The spatiotemporal resolutions for the three simulations are tabulated in Table II. Figure 10 compares the vertical center of mass displacement (heave motion) as a function of time for the three grid resolutions. The results suggest that the medium grid resolution with ∼20 cells per diameter of the buoy, 150 cells per wavelength, and 5 cells per wave height is sufficient to resolve the WSI dynamics adequately. We use this resolution to simulate the rest of the CFD cases considered in this paper.

In this section, we compare the simulated dynamics of a submerged cylindrical point absorber using the Cummins equation-based Simulink model and fully resolved CFD model implemented in IBAMR. The dynamics of the buoy and the generated PTO power at the steady state are compared using the two models. We take the same buoy setup, PTO coefficients, and wave characteristics of Sec. III F 2 to compare the dynamics. Since the Cummins equation is based on linear wave theory, the wave parameters should correspond to first-order Stokes/Airy wave theory. Using the wave classification phase space described by Le Méhauté,⁸² it can be verified that wave parameters of Sec. III F 2 closely correspond to the linear wave theory. Next, we compare the dynamic response of the buoy in various degrees of freedom using the two models to assess the performance of Cummins equation-based models for practical wave energy conversion applications.

We first simulate one degree of freedom (one-DOF) buoy using the two models. Figure 11 compares the heave dynamics of the submerged buoy and the generated PTO power during the steady state. As observed in Fig. 11(a), the Cummins model overpredicts the heave amplitude due to the fact that linear potential theory overestimates the Froude–Krylov forces or the wave excitation loads on the submerged buoy.^10,23,25 The generated PTO power is also higher using the Cummins model, as observed in Fig. 11(b). Nevertheless, the dynamics obtained using the two models are qualitatively the same.

Next, we simulate a two-DOF buoy that oscillates in heave and surge directions using the two models. Figure 12 compares the heave and surge dynamics of the submerged buoy and the generated PTO power during the steady state. The amplitude of the heave motion remains almost the same as that of the one-DOF buoy. The amplitude of surge is lower than heave, which can be observed in Fig. 12(d). However, both models estimate slightly higher PTO power (7%–10% more for this case) compared to the one-DOF buoy, as shown in Fig. 12(b). This is in agreement with the work of Falnes,⁸³ who also predicts an increase in the power absorption efficiency of a point absorber with more than one degree of freedom using theoretical analyses.

Note that for the two-DOF buoy, the fully resolved WSI model is able to capture the slow drift phenomenon induced by the background flow. On the contrary, the potential flow model is unable to capture the slow drift of the buoy, as evidenced in Fig. 12(c). To compare the amplitude of the high-frequency surge motion, Fig. 12(d) is plotted by eliminating the slow drift frequency from the surge dynamics obtained from CFD using a high-pass filter. Similar to the heave dynamics, the potential flow model overpredicts the surge amplitude. The inability of the Cummins model to capture the slow drift phenomenon can be explained as follows. According to Chakrabarti,⁸⁴ the drift forces estimated by the potential flow theory tend to zero for ka < 0.5, in which k = 2π/λ is the wave number and $a=H/2$ is the wave amplitude. In the present simulation, ka = 0.0258, and consequently, the Cummins equation-based model underpredicts the drift force. We remark that estimating drift forces is an important criterion for designing a robust mooring system for a practical WEC.⁸⁵ Therefore, fully resolved simulations are more reliable for such calculations.

Finally, we analyze the three-DOF buoy that can undergo heave, surge, and pitching motions. Figure 13 shows the heave, surge, and pitch dynamics of the buoy. Compared to the two-DOF buoy dynamics, the heave and surge dynamics [Figs. 13(a) and 13(b), respectively] show a larger oscillatory behavior when the pitch is included in the WSI model. This can be further confirmed from Fig. 14, which shows the trajectory traced by the center of mass points of two- and three-DOF buoys. Figure 14(a) shows the complete trajectory of the center of the mass point, including the slow drift, traced during few steady-state periods of motion (64 s–70 s). After eliminating the slow drift from the complete trajectory using a high-pass filter, Fig. 14(b) shows a larger trajectory for the three-DOF buoy. The buoy undergoes an average rotation of 5° at the steady state, which can be observed in Fig. 13(c). In contrast, the potential flow model results in an insignificant pitch. This can be explained as follows. For a fully submerged axisymmetric body, the wave excitation forces due to the induced fluid pressure do not produce any rotational torque as the pressure forces pass through the center of mass of the buoy. Moreover, the hydrodynamic coupling between various degrees of freedom is also weak when linear potential flow equations are considered.⁸⁶ Consequently, the resulting pitch velocity $θ̇$ and the viscous drag torque M_drag,θ are negligible. The simulated dynamics of a three-DOF buoy using the potential flow model is, therefore, quite similar to the prior two-DOF buoy case.

Figure 15 shows the power production of the three-DOF buoy and contrasts it with two-DOF PTO power generation. As observed in Fig. 15(b), the pitching motion does not contribute substantially to the power generation compared to heave and surge motions for the simulated wave parameters. However, for higher wave heights and, hence, more energetic waves, power contribution from pitching motion could, in general, be significant. Moreover, for point absorbers such as ISWEC, the pitch motion of the hull is the primary source of power production. Therefore, the fully resolved WSI framework can be reliably used to resolve all modes of motion of a PA.

Sections IV A–IV C highlighted the differences in the dynamic response of the buoy in one, two, and three degrees of freedom. Here, we analyze the differences in the vortex shedding pattern for a buoy having different degrees of freedom. Figure 16 shows the vorticity production and transport due to WSI of one- and two-DOF buoys. The vortex production of a three-DOF buoy is qualitatively similar to the two-DOF buoy and is not shown.

As observed in Fig. 16, the one-DOF buoy sheds vortices in an inclined direction, whereas the two-DOF one sheds it in the vertical direction. This observation can be explained as follows. For the one-DOF buoy, vortex structures generated at the surface of the buoy is transported by wave induced horizontal flow velocity u, along with a vertical flow produced by the heaving buoy. Therefore, the net transport of vorticity is in an inclined direction. For the two-DOF buoy, Fig. 17(a) plots the undisturbed horizontal flow velocity u at the bottom most point of the buoy, when the buoy is at its initial equilibrium position. The heave and surge displacements and the surge velocity of the buoy using the CFD model are also shown in Fig. 17(a). It can be observed that the surge displacement and the horizontal flow velocity are in phase with each other. In addition, the difference in the magnitude of the surge velocity and the u component of the flow velocity is small. These two factors mitigate flow separation in the horizontal direction. In contrast, the difference in the magnitude of the heave velocity of the buoy and v component of the flow field is relatively large, as shown in Fig. 17(b). It can also be observed that the heave displacement and heave velocity of the buoy are out of phase (π/2 and π radians, respectively) with the $v$ component of the flow field. Therefore, for the two-DOF buoy, vortex shedding primarily happens due to the heave motion. Moreover, a vortex pair is shed at the end of a heave stroke (at the extremum of the heave curve) when the flow and heave velocities begin to increase from their zero values, but in opposite directions; see Fig. 17(b). For the two-DOF (and also three-DOF) buoy, the surge velocity of the buoy is out of phase (by ∼π/2 radians) with the horizontal flow velocity. At the time of vortex shedding (which will be at some but not all peaks of the heave curve in this case), the surge velocity of the buoy and the horizontal component of the flow velocity are of similar magnitude but have opposite signs [Fig. 17(a)]. The difference in the magnitude of surge and u velocities becomes even smaller if the u velocity is considered at the lowest y location attained by the buoy. Hence, the main component of the transport velocity is in the vertical direction when a vortex pair is released for the two-DOF buoy.

For these simulations, the PTO coefficients are taken to be k_PTO = 1995.2 N/m and b_PTO = 80.64 N s/m. Increasing or decreasing these coefficients by a factor of five changed the phase difference between the surge and horizontal velocities only slightly, suggesting that the two- or three-DOF buoy will shed vortices in the vertical direction irrespective of (physically reasonable) PTO coefficients.

In our simulations, we impose zero pressure boundary condition on the top boundary, which is taken rather close to the air–water interface to reduce the overall computational cost. Within this narrow space, the vortex structures in the air phase can reach the top boundary, especially when they are generated by waves of large amplitude. To reduce the boundary artifacts (they do not affect the wave and buoy dynamics though) observed in Fig. 16, the top boundary could be modeled further away, or a vorticity damping zone can be added near the boundary to dissipate the incoming vortex structures. At the bottom of the NWT, we use no-slip velocity boundary conditions. This produces a thin boundary layer region at the channel bottom, which can be observed in Fig. 16.

In this section, we use the WSI framework to study the conversion efficiency of a three-DOF submerged cylindrical buoy. More specifically, we analyze the performance of the buoy by varying

• PTO stiffness and damping coefficients,

• density of the submerged buoy, and

• wave height.

The wave absorption or the conversion efficiency of a wave energy converter is defined as the ratio of the mean absorbed power $P¯absorbed$ to the mean wave power $P¯wave$⁠. For a regular wave, the wave absorption efficiency can be written as

Here, c_g is the wave group velocity, and P_absorbed(t) is the instantaneous power absorbed with the PTO damper,

Wave energy absorption efficiency of a converter is maximized if it resonates with the incoming waves.⁸³ In this regard, PTO stiffness and damping coefficients are two tuning parameters that can be effectively controlled to synchronize the natural frequency of the mechanical oscillator with the incoming wave frequency. Accordingly, the linear wave theory suggests a reactive control-based strategy to select the PTO stiffness and damping coefficients as

in which M is the mass of the converter, and A₃₃(ω) and B₃₃(ω) are the frequency-dependent added mass and radiation damping coefficients in the heave direction, respectively. The length-normalized hydrodynamic coefficients of a submerged cylinder of length 8D for a wave period range of $T∈[0.625−1.1]$ s are obtained using AQWA. Figure 18 shows the variation of frequency-dependent hydrodynamic coefficients in the selected $T$ range.

To study the effect of PTO stiffness and damping on the converter efficiency, we perform fully resolved WSI simulations of the submerged buoy for two sets of PTO coefficients:

• Reactive control coefficients: $kPTO=ω2M+A33(ω)$ and b_PTO = B₃₃(ω).

• Optimal control coefficients: $kPTO=ω2M+A33(ω)$ and b_PTO = 80.64 N m/s.

Figure 19 shows the absorption efficiency of the buoy at various wave frequencies using reactive control and optimal control PTO coefficients. The curve obtained using reactive control indicates that the converter efficiency saturates around $T=0.9$ s. Therefore, a higher damping coefficient could be used to narrow down the optimal performance period range and possibly enhance the efficiency of the converter. This is achieved by using the optimal control damping coefficient value of b_PTO = 80.64 N m/s, which is approximately four times larger than the maximum value of B₃₃(ω) predicted by the linear wave theory for a wave period around $T=0.9$ s [see Fig. 18(b)]. This optimal value of b_PTO increases the absorption efficiency of the buoy for all wave frequencies, which suggests that the reactive control damping coefficients predicted by the linear wave theory do not lead to an optimal performance. Further increasing the b_PTO value did not enhance the performance of the converter significantly (data not shown). Similar observations were made in the work of Anbarsooz et al.²⁵ Note that extremely large values of b_PTO would lead to over-damping of the system and, therefore, should also be avoided. Table III shows the characteristics of the simulated waves including their steepness values. The reactive control curve of Fig. 19 shows a sub-optimal performance of the buoy for lower wave periods at which the wave steepness is high, i.e., the ratio $H/λ>0.01$⁠. Therefore, linear wave theory does not predict optimal PTO coefficients for steeper waves as it assumes low values of amplitude to the wavelength ratio. The sub-optimal values of PTO coefficients predicted by linear wave theory are also reported in Ref. 25.

Next, we consider the effect of buoy density on its conversion performance. We fix the PTO coefficients to k_PTO = 1995.2 N/m and b_PTO = 80.64 N m/s but vary the mass density for this study. Three relative densities of the buoy are considered: ρ_s/ρ_w = 0.5, 0.7, and 0.9. Figure 20(a) shows that each density curve results in an optimal performance of the converter in a certain range of wave periods. This optimal range can also be predicted from the reactive control theory,

Table IV shows the optimal performing wave period range calculated using Eq. (44) for a fixed value of k_PTO = 1995.2 N/m and lowest and highest values of A₃₃(ω) from Fig. 18(a). The simulated results and the analytically predicted range suggest that the optimal period range increases with an increase in mass density of the buoy. Note that since k_PTO = 1995.2 N/m is an optimal reactive PTO stiffness for a relative density of 0.9, the efficiency curve is higher for ρ_s/ρ_w = 0.9 compared to other densities. Figure 20(b) plots the reactive PTO stiffness for different mass densities of the buoy using Eq. (44): $kPTO=4π2M+A33(ω)/T2$⁠. Lower PTO stiffness coefficients are obtained for lower mass densities from this relationship. Therefore, lowering the k_PTO value from 1995.2 N/m (or in other words using a more optimal k_PTO value) is expected to increase the conversion efficiency for lighter buoys.

We remark that the relative density of the buoy directly influences the permanent tension in its mooring line. For lower mass density values, there is an increased upward hydrostatic force, which is written in the last column of Table IV. The submerged buoy having 50% water density has five times greater permanent load than a buoy having 90% water density. This factor should be considered while selecting the wave energy converter density.

As a last parametric study of converter performance, we simulate the dynamic response of the submerged point absorber in regular waves of different heights. The buoy is subject to incident waves of height $H=$ 0.01 m, 0.02 m, and 0.03 m having time periods in the range $T∈[0.625−1.1]$ s. The absorption efficiency of the buoy for different wave heights is shown in Fig. 21. As observed in Fig. 21, the absorption efficiency decreases with increased wave heights. Waves with greater heights are more energetic, and they result in an increased dynamic response of the submerged object. However, the point absorber fails to absorb a significant portion of the available wave energy. One possible way to achieve an optimal performance for more energetic waves is to increase the size of the power take-off unit and the WEC device.^3,83,87 However, bigger wave energy converters are more costly, and therefore, a balance between cost and efficiency should be considered during the design stage of WECs.³ 

In this study, we compared the dynamics of a three-DOF cylindrical buoy using potential flow and CFD models. The potential flow model is based on the time-domain Cummins equation, whereas the CFD model employs the FD/BP method—a fully Eulerian technique for modeling FSI problems. An advantage of the FD/BP method over the FD/DLM method is its ability to incorporate external forces and torques into the equations of motion, and it enabled us to solve the coupled translational and rotational degrees of freedom of the buoy.

The comparison of the dynamics shows that the Cummins model overpredicts the amplitude of heave and surge motions, whereas it results in an insignificant pitch of the buoy. Moreover, it does not capture the slow drift phenomena in the surge dynamics. The CFD model was, then, used to study the wave absorption efficiency of the converter under varying PTO coefficients, mass density of the buoy, and incoming wave heights. It is demonstrated that the PTO coefficients predicted by the linear potential theory are sub-optimal for waves of moderate and high steepness. The wave absorption efficiency improves significantly when higher values of PTO damping coefficients are used. Simulations with different mass densities of the buoy show that converters with low mass density have an increased permanent load in their PTO and mooring lines. Moreover, the mass density also influences the range of resonance periods of the device. Finally, simulations with different wave heights show that at higher heights, the wave absorption efficiency of the converter decreases, and a large portion of available wave power remains unabsorbed.

The data that support the findings of this study are openly available in IBAMR GitHub repository at https://github.com/IBAMR/IBAMR.

A.P.S.B. acknowledges support from the NSF, Award No. OAC 1931368. The NSF XSEDE and SDSU Fermi compute resources are particularly acknowledged.