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

TL;DR

This study compares potential flow and CFD models for a submerged cylindrical wave energy converter, revealing differences in predicted motions and efficiencies, and exploring effects of PTO tuning, mass density, and wave height on performance.

Contribution

It provides a detailed comparison between potential flow and CFD models for a 3-DOF submerged wave energy converter, highlighting the limitations and advantages of each approach.

Findings

Potential flow over-predicts heave and surge amplitudes.

CFD captures slow drift and detailed dynamics.

Higher PTO damping improves wave absorption efficiency.

Abstract

In this paper we compare the heave, surge, and pitch dynamics of a submerged cylindrical point absorber, simulated using potential flow and fully-resolved computational fluid dynamics (CFD) models. The potential flow model is based on the time-domain Cummins equation, whereas the CFD model uses the fictitious domain Brinkman penalization (FD/BP) technique. The submerged cylinder is tethered to the seabed using a power take-off (PTO) unit which restrains the heave, surge, and pitch motions of the converter, and absorbs energy from all three modes. It is demonstrated that the potential theory over-predicts the amplitudes of heave and surge motions, whereas it results in an insignificant pitch for a fully-submerged axisymmetric converter. It also under-estimates the slow drift of the buoy, which the CFD model is able to capture reliably. Further, we use fully-resolved CFD simulations to study the performance of a three degrees of freedom (DOF) cylindrical buoy 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 than the predicted value of the PTO damping is selected. Simulations with different mass densities of the buoy show that converters with low mass densities have an increased tension 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.

Tables

Table 1. Table 1: Grid resolutions to simulate the forced damped oscillation of a cylinder.

Resolution	Coarse	Medium	Fine
$N_x$	50	100	400
$N_y$	100	200	800
$\Delta t$ (s)	$2.5\times {10}^{-4}$	$2.5\times {10}^{-4}$	${10}^{-4}$

Table 2. Table 2: Grid resolution to simulate the heave dynamics of a submerged point absorber.

Resolution	Coarse	Medium	Fine
$N_y$	125	250	400
$N_x$	750	1500	2400
$\Delta t$ (s)	${10}^{-3}$	$5\times {10}^{-4}$	$2.5\times {10}^{-4}$

Table 3. Table 3: Wave characteristics.

Wave period (s)	Wave length (m)	Wave steepness
0.625	0.6099	0.0328
0.666	0.6925	0.0288
0.714	0.7959	0.0251
0.7692	0.9235	0.0216
0.833	1.0822	0.0184
0.909	1.2856	0.0155
1	1.5456	0.0129
1.11	1.875	0.0107

Table 4. Table 4: Characteristics of the submerged buoy with different relative densities.

Relative density	Mass (kg)	Optimal wave period range (s)	Permanent tension (N)
0.5	10.30	[0.74-0.8]	101.08
0.7	14.42	[0.8-0.85]	60.65
0.9	18.54	[0.85-0.9]	20.21

Equations

$${\nabla }^2\Phi =0,$$

$$p(\mathbf{x},t)$$

$$=-{\rho }_{\text{w}}\frac{\partial \Phi }{\partial t}(\mathbf{x},t)-{\rho }_{\text{w}}gy,$$

$$\Phi ={\Phi }_{\text{I}}+{\Phi }_{\text{D}}+{\Phi }_{\text{R}},$$

$$\left(\mathbf{M}+\mathbf{A}(\omega )\right)\ddot{\mathbf{X}}+\mathbf{B}(\omega )\dot{\mathbf{X}}+{\mathbf{K}}_{\text{hydro}}\mathbf{X}={\mathbf{F}}_{\text{e}}(\omega ).$$

$$\left(\mathbf{M}+{\mathbf{A}}_{\infty }\right)\ddot{\mathbf{X}}(t)+{\int }_0^t{\mathbf{K}}_{\text{r}}(t-\tau )\dot{\mathbf{X}}(\tau )\text{d}\tau +{\mathbf{K}}_{\text{hydro}}\mathbf{X}(t)={\mathbf{F}}_{\text{e}}(t)+{\mathbf{F}}_{\text{ext}}(\mathbf{X}(t),\dot{\mathbf{X}}(t),t),$$

$${\mathbf{A}}_{\infty }=\text{scalebox}0.8\underset{\omega \to \infty }{\lim }\mathbf{A}(\omega ),$$

$${\mathbf{F}}_{\text{r}}(t)={\int }_0^t{\mathbf{K}}_{\text{r}}(t-\tau )\dot{\mathbf{X}}(\tau )\text{d}\tau \simeq \{\begin{array}{l}{\dot{\mathbf{\zeta }}}_{\text{r}}(t)={\mathcal{A}}_{\text{r}}{\mathbf{\zeta }}_{\text{r}}(t)+{\mathcal{B}}_{\text{r}}\dot{\mathbf{X}}(t),\\ {\mathbf{F}}_{\text{r}}(t)={\mathcal{C}}_{\text{r}}{\mathbf{\zeta }}_{\text{r}}(t),\end{array}$$

$$\left(\mathbf{M}+{\mathbf{A}}_{\infty }\right)\ddot{\mathbf{X}}(t)+{\mathbf{F}}_{\text{r}}$$

$${\mathbf{F}}_{\text{hydrostatic}}$$

$${\mathbf{F}}_{\text{m}}$$

$${\mathbf{F}}_{\text{PTO}}$$

$${\mathbf{F}}_{\text{drag,x}}$$

$${\mathbf{F}}_{\text{drag,y}}$$

$${\mathbf{M}}_{\text{drag},\theta }$$

$$\frac{\partial \rho \mathbf{u}(\mathbf{x},t)}{\partial t}+\nabla \cdot \text{mbox}\rho \mathbf{u}(\mathbf{x},t)\mathbf{u}(\mathbf{x},t)$$

$$\nabla \cdot \text{mbox}\mathbf{u}(\mathbf{x},t)$$

$${\mathbf{f}}_{\text{c}}(\mathbf{x},t)$$

$$\rho (\mathbf{x},t)$$

$$\mu (\mathbf{x},t)$$

$$\frac{\partial \varphi }{\partial t}+\nabla \cdot \text{mbox}(\varphi \mathbf{u})$$

$$\frac{\partial \psi }{\partial t}+\nabla \cdot \text{mbox}(\psi \mathbf{u})$$

$$\frac{\partial \varphi }{\partial \tau }+\text{sgn}\left({\tilde{\varphi }}^{n+1}\right)\left(\Vert \nabla \varphi \Vert -1\right)=0,$$

$$\varphi (\mathbf{x},\tau =0)={\tilde{\varphi }}^{n+1}(\mathbf{x}),$$

$${\tilde{H}}_{i,j}^{\text{flow}}$$

$${\tilde{H}}_{i,j}^{\text{body}}$$

$${\Im }_{i,j}^{\text{flow}}={\Im }_{\text{l}}+({\Im }_{\text{g}}-{\Im }_{\text{l}}){\tilde{H}}_{i,j}^{\text{flow}}.$$

$${\Im }_{i,j}^{\text{full}}={\Im }_{\text{s}}+({\Im }_{i,j}^{\text{flow}}-{\Im }_{\text{s}}){\tilde{H}}_{i,j}^{\text{body}}.$$

$$\frac{{\varphi }^{n+1,k+1}-{\varphi }^n}{\Delta t}+Q\left({\mathbf{u}}^{n+\frac{1}{2},k},{\varphi }^{n+\frac{1}{2},k}\right)$$

$$\frac{{\psi }^{n+1,k+1}-{\psi }^n}{\Delta t}+Q\left({\mathbf{u}}^{n+\frac{1}{2},k},{\psi }^{n+\frac{1}{2},k}\right)$$

$$\frac{{\breve{\mathbf{\rho }}}^{n+1,k+1}{\mathbf{u}}^{n+1,k+1}-{\mathbf{\rho }}^n{\mathbf{u}}^n}{\Delta t}+{\mathbf{C}}^{n+1,k}$$

$${\nabla }_h\cdot {\mathbf{u}}^{n+1,k+1}$$

$${\mathbf{f}}_{\text{c}}^{n+1,k+1}=\frac{\tilde{\chi }}{\kappa }\left({\mathbf{u}}_{\text{b}}^{n+1,k+1}-{\mathbf{u}}^{n+1,k+1}\right),$$

$${\mathbf{u}}_{\text{b}}^{n+1,k+1}={\mathbf{U}}_{\text{r}}^{n+1,k+1}+{\mathbf{W}}_{\text{r}}^{n+1,k+1}\times \left(\mathbf{x}-{\mathbf{X}}_{\text{COM}}^{n+1,k+1}\right).$$

$${\mathbf{X}}_{\text{COM}}^{n+1,k+1}={\mathbf{X}}_{\text{COM}}^n+\Delta t{\mathbf{U}}_{\text{r}}^{n+1,k+1}.$$

$${\text{M}}_{\text{b}}\frac{{\mathbf{U}}_{\text{r}}^{n+1,k+1}-{\mathbf{U}}_{\text{r}}^n}{\Delta t}$$

$${\mathbf{I}}_{\text{b}}\frac{{\mathbf{W}}_{\text{r}}^{n+1,k+1}-{\mathbf{W}}_{\text{r}}^n}{\Delta t}$$

Full text

Keywords: incompressible Navier-Stokes equations, Brinkman penalization method, numerical wave tank, potential flow theory, Cummins equation

1 Introduction

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. In spite of the abundantly available wave power in the oceans and seas worldwide [1], and research efforts dating back since the seventies [2], 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 wave length of the site [3, 4]. Depending upon 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, 6, 7]. PowerBuoy is a two-body FPA developed by Ocean Power Technologies that uses heave mode to extract energy from the waves [8, 9, 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 (SPA) 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 [11, 12]. An added advantage of SPAs is their zero visual impact on the ocean or sea shorelines [13]. Currently, efforts are underway that are testing point absorber devices at various locations around the world, including but not limited to, the Pacific Ocean [14], the Atlantic Ocean [15], and the Mediterranean Sea [16, 17]. Some of us are also directly involved in testing and improving WEC devices at various sea locations [5, 6, 7, 17].

Numerical models based on frequency- or time-domain methods are commonly used to study the performance of point absorbers [18, 19, 20]. The hydrodynamic loads in these methods are calculated using the boundary element method (BEM) approach which is based on 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 like wave-breaking and wave-overtopping. Moreover, these methods over-predict the dynamics and the wave absorption efficiency of the WEC systems [10]. An improvement over LPF based models is the fully nonlinear potential flow (FNPF) formulation, that permits large-amplitude displacements of the WECs and modeling of nonlinear free-surface [21]. 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 [22, 23, 10, 24, 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, 27, 28, 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 [30], elastic boundaries [31], 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 by employing the Brinkman penalization (BP) approach [36, 37, 38]. Ghasemi et al. [24] and Anbarsooz et al. [25] 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 111For example in cases where linear and angular momentum due to additional external force and torque can be included during the momentum conservation stage of the FD/DLM method.. However, applying external torque on the structure is not straightforward for FD/DLM method because 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 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 (3-DOF) buoy under the action of external forces and torques, we employ the more versatile FD/BP approach in our CFD model. To our knowledge, this works presents the first application of 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 over-predicts 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 3-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. 2. Next, we describe the continuous and discrete equations for the multiphase fluid-structure system in Secs. 3.1 and 3.4, respectively. Validation cases for the FD/BP framework are presented in Sec. 3.6. The tests also highlight the solver stability in presence of large density contrasts. Sec. 4 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. 5.

2 Numerical model based on the potential flow theory

2.1 State-space fluid-structure interaction formulation

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

[TABLE]

using suitable kinematic and dynamic boundary conditions [39, 40]. The fluid velocity is expressed as gradient of the velocity potential, $\bm{u}=(u,v)=\bm{\nabla}\Phi$. Once the solution to the Laplace equation is found, the fluid pressure $p$ is obtained from the linearized Bernoulli equation $\partial\Phi/\partial t+p/\rho_{\text{w}}+gy=0$ as

[TABLE]

in which $g$ is the acceleration due to gravity, $\rho_{\text{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 [41]

[TABLE]

in which $\Phi_{\text{I}}$ is the undisturbed (assuming no body in the domain) wave potential of the incident wave, $\Phi_{\text{D}}$ is the diffraction potential of the incident wave about the stationary body, and $\Phi_{\text{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 $\bm{\bm{X}}$ due to a monochromatic harmonic wave excitation of angular frequency $\omega$

[TABLE]

Here $\bm{M}$ is the mass matrix of the buoy, $\bm{A}(\omega)$ and $\bm{B}(\omega)$ are the frequency-dependent added mass and damping matrix of the buoy, respectively, $\bm{K}_{\text{hydro}}$ is the linear hydrostatic stiffness arising from the buoyancy force for a floating buoy, and $\bm{F}_{\text{e}}(\omega)$ is the wave excitation force (Froude-Krylov and diffraction). The displacement, velocity, and acceleration of the body are denoted by $\bm{X}$, $\dot{\bm{X}}$, and $\ddot{\bm{X}}$, respectively. The dimensions of the matrices and vectors depend on the degrees of freedom [42]. The frequency-domain models are typically resolved using the boundary element method (BEM) and this approach has been widely adopted in commercial codes like ANSYS AQWA [43] or WAMIT [44].

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

[TABLE]

in which $\bm{A}_{\infty}$ is the added mass at infinite frequency, given by

[TABLE]

and $\bm{K}_{\text{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 $\bm{F}_{\text{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 [49], Jefferys [50], Damaren [51], McCabe et al. [52], and Clément [53]. In this work we follow the state-space representation approach of Fossen and Perez [54, 55] to approximate the radiation convolution integral by

[TABLE]

in which $\bm{\mathcal{A}}_{\text{r}}$, $\bm{\mathcal{B}}_{\text{r}}$, and $\bm{\mathcal{C}}_{\text{r}}$ are the state-space matrices for carrying out time-domain analysis. It is also possible to evaluate the RIRF in frequency domain and subsequently transform it back to the time domain [56, 57, 50].

Fig. 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 $\bm{F}_{\text{drag}}$, and the PTO stiffness and damping loads, which are denoted by $\bm{F}_{\text{m}}$ and $\bm{F}_{\text{PTO}}$, respectively. Since the point absorber is taken to be completely submerged under water, the hydrostatic stiffness arising from the buoyancy force is $\bm{K}_{\text{hydro}}=\bm{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 $\rho_{\text{s}}$, diameter $D$, and volume $V_{\text{buoy}}=\pi D^{2}/4$ reads as

[TABLE]

in which $k_{\text{PTO}}$ and $b_{\text{PTO}}$ are the stiffness and damping constants of the PTO, $\Delta l$ is the elongation of the PTO from a reference length, $\Delta l_{0}=|\bm{F}_{\text{hydrostatic}}|/k_{\text{PTO}}$ is the permanent extension of the PTO to balance the hydrostatic force $\bm{F}_{\text{hydrostatic}}$, and $\hat{\bm{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. [19]

[TABLE]

in which $C_{x}$, $C_{y}$, and $C_{\theta}$ are the drag coefficients in the surge ($\hat{\bm{x}}$), heave ($\hat{\bm{y}}$) and pitch ($\hat{\bm{z}}$) directions, respectively, and $S_{x}=S_{y}=D$ is the planar cross-section area of the disk.

2.2 Calculation of 2D coefficients

To solve the time-domain Cummins Eq. (8) for a cylindrical shaped submerged point absorber, we use ANSYS AQWA to extract the frequency-dependent quantities such as: 1) added mass matrix $\bm{A}(\omega)$; 2) wave excitation force $\bm{F}_{\text{e}}(\omega)$; 3) and radiation damping matrix $\bm{B}(\omega)$. Further post-processing is done to transform the AQWA data to time-domain amenable quantities like $\bm{A}_{\infty}$, $\bm{\mathcal{A}}_{\text{r}}$, $\bm{\mathcal{B}}_{\text{r}}$, and $\bm{\mathcal{C}}_{\text{r}}$. The approach described in Fossen and Perez [54, 55] is implemented using custom MATLAB scripts for this purpose. Finally, the time-domain Eqs. (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$. Fig. 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 $\bm{A}(\omega)$ and radiation damping $\bm{B}(\omega)$ matrices in Sec. 5.1.

3 Numerical model based on the incompressible Navier-Stokes equations

We use 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, 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 [58], which is an open-source C++ library providing support for immersed boundary methods with adaptive mesh refinement [59]. IBAMR relies on SAMRAI [60, 61] for Cartesian grid management and the AMR framework and on PETSc [62, 63, 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.

3.1 Continuous equations of motion

We state the governing equations for a coupled multiphase fluid-structure system in a fixed region of space $\Omega\subset\mathbb{R}^{d}$, for $d=2$ spatial dimensions. A fixed Eulerian coordinate system $\bm{x}=(x,y)\in\Omega$ is used to describe the momentum and incompressibility of fluid and solid domains. The spatial location of the immersed body $\Omega_{\rm b}(t)\subset\Omega$ is tracked using an indicator function $\chi(\bm{x},t)$ which is nonzero in the solid domain and zero in the fluid domain $\Omega_{\rm f}(t)\subset\Omega$. The time-varying solid and fluid domains are non-overlapping and their union occupies the entire domain $\Omega=\Omega_{\rm f}(t)\cup\Omega_{\rm b}(t)$. We employ spatially and temporally varying density $\rho(\bm{x},t)$ and viscosity $\mu(\bm{x},t)$ fields to describe the coupled three phase system. The equations of motion for the multiphase fluid-structure interaction system read as

[TABLE]

Eq. (15) is the conservative form of momentum equation, whereas Eq. (16) describes the incompressibility of the system. The quantity $\bm{u}(\bm{x},t)$ expresses the velocity, $p(\bm{x},t)$ the mechanical pressure, and $\bm{f}_{\text{c}}(\bm{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 negative $y$-direction $\bm{g}=(0,-g)$. In the FD/BP method, the rigidity-enforcing constraint force $\bm{f}_{\text{c}}(\bm{x},t)$ is defined as a penalization force that enforces a rigid body velocity $\bm{u}_{\text{b}}(\bm{x},t)$ in $\Omega_{b}(t)$. By treating the immersed structure as a porous region of vanishing permeability $\kappa\ll 1$, the Brinkman penalized constraint force is formulated as

[TABLE]

Sec. 3.5.3 describes the fluid-structure coupling algorithm and the rigid body velocity $\bm{u}_{\text{b}}(\bm{x},t)$ calculation.

3.2 Phase tracking

We use scalar level set function $\phi(\bm{x},t)$ to identify liquid and gas regions, $\Omega_{\text{l}}\subset\Omega$ and $\Omega_{\text{g}}\subset\Omega$, respectively, in the computational domain. The combined liquid and gas region denotes the fluid region described in the previous section, i.e., $\Omega_{\text{l}}\cup\Omega_{\text{g}}=\Omega_{\rm f}$. The zero-contour of $\phi$ function implicitly defines the liquid-gas interface $\Gamma(t)=\Omega_{\text{l}}\cap\Omega_{\text{g}}$. Similarly, the Eulerian indicator function $\chi(\bm{x},t)$ of the immersed body is expressed in terms of the level set function $\psi(\bm{x},t)$ 222$\chi(\bm{x},t)=1-H^{\rm body}$, in which $H^{\rm body}$ is the body Heaviside function defined in Eq. 25 using $\psi(\bm{x},t)$. , whose zero-contour implicitly defines the surface of the body $S_{\rm b}(t)=\partial V_{\rm b}(t)$; see Fig. 3. Without loss of generality, $\phi$ level set (signed distance) values are taken to be negative in the liquid phase and positive in the air phase. Similarly, $\psi$ 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 like density and viscosity can be conveniently expressed as a function of $\phi(\bm{x},t)$ and $\psi(\bm{x},t)$ fields

[TABLE]

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 conservative form reads as

[TABLE]

The signed distance property of $\phi$ and $\psi$ is generally disrupted under linear advection, Eqs. (20) and (21). Let $\widetilde{\phi}^{n+1}$ denote the flow level set function following an advective transport after time stepping through the interval $\left[t^{n},t^{n+1}\right]$. The flow level set is reinitialized to obtain a signed distance field $\phi^{n+1}$ by computing a steady-state solution to the Hamilton-Jacobi equation

[TABLE]

which yields a solution to the Eikonal equation $\|\nabla\phi\|=1$ at the end of each time step. We refer the readers to [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 [66], and an immobile interface condition described by Son [67]. Both 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 $\psi^{n+1}$ is analytically calculated and reinitialized by using the new location of center of mass at $t^{n+1}$. For more complex structures, computational geometry techniques can be employed to compute the signed distance function 333In our code we use surface triangulation of complex geometries to compute the signed distance function..

3.3 Spatial discretization

The continuous equations of motion Eqs. (15)-(16) are discretized on a staggered Cartesian grid as shown in Fig. 3. The discretized domain is made up of $N_{x}\times N_{y}$ rectangular cells that cover the physical domain $\Omega$. The mesh spacing in the two directions is denoted by $\Delta x$ and $\Delta 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 $\bm{x}_{i,j}=\left((i+\frac{1}{2})\Delta x,(j+\frac{1}{2})\Delta y\right)$ for $i=0,\ldots,N_{x}-1$ and $j=0,\ldots,N_{y}-1$. The physical location of the vertical cell face is half a grid space away from $\bm{x}_{i,j}$ in the $x$-direction and is given by $\bm{x}_{i-\frac{1}{2},j}=\left(i\Delta x,(j+\frac{1}{2})\Delta y\right)$. Similarly, $\bm{x}_{i,j-\frac{1}{2}}=\left((i+\frac{1}{2})\Delta x,j\Delta y\right)$ is the physical location of the horizontal cell face that is half a grid cell away from $\bm{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 $\phi_{i,j}^{n}\approx\phi\left(\bm{x}_{i,j},t^{n}\right)$, $\psi_{i,j}^{n}\approx\psi\left(\bm{x}_{i,j},t^{n}\right)$, $p_{i,j}^{n}\approx p\left(\bm{x}_{i,j},t^{n}\right)$, $\rho_{i,j}^{n}\approx\rho\left(\bm{x}_{i,j},t^{n}\right)$ and $\mu_{i,j}^{n}\approx\mu\left(\bm{x}_{i,j},t^{n}\right)$, respectively; some of these quantities are interpolated onto the required degrees of freedom as needed (see [65] for further details). Here, $t^{n}$ denotes the time at time step $n$. The velocity degrees of freedom are defined on the cell faces and are denoted by $u_{i-\frac{1}{2},j}^{n}\approx u\left(\bm{x}_{i-\frac{1}{2},j},t^{n}\right)$, and $v_{i,j-\frac{1}{2}}^{n}\approx v\left(\bm{x}_{i,j-\frac{1}{2}},t^{n}\right)$. 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., $\nabla\approx\nabla_{h}$. Further details on the spatial discretization can be obtained in prior studies [65, 68, 69, 70].

3.4 Solution methodology

3.4.1 Material property specification

Discretely, we use smoothed Heaviside functions to transition from one material phase to the other. The transition zone occurs within $n_{\text{cells}}$ grid cells on either side of water-air interface $\Gamma$ or fluid-solid interface $S_{\rm b}$. Correspondingly, two numerical Heaviside functions are defined:

[TABLE]

The number of transition cells across $\Gamma$ or $S_{\rm b}$ interface is assumed to be the same. This is not a strict requirement of our numerical method, but is 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 $\Im$ (such as $\rho$ or $\mu$) in the whole domain $\Omega$:

• •

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

[TABLE]

• •

Next, $\Im^{\text{flow}}$ is corrected by accounting for the structural material property to obtain $\Im_{i,j}^{\text{full}}$ throughout the computational domain

[TABLE]

3.5 Time stepping scheme

We employ a fixed-point iteration time stepping scheme with $n_{\text{cycles}}=2$ cycles per time step to advance quantities of interest at time level $t^{n}$ to time level $t^{n+1}=t^{n}+\Delta 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$, $\bm{u}^{n+1,0}=\bm{u}^{n}$, $p^{n+\frac{1}{2},0}=p^{n-\frac{1}{2}}$, $\phi^{n+1,0}=\phi^{n}$, and $\psi^{n+1,0}=\psi^{n}$. For $n=0$ initial time level, all of the physical quantities have a prescribed initial condition.

3.5.1 Level set advection

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

[TABLE]

in which $Q(\cdot,\cdot)$ represents an explicit piecewise parabolic method (xsPPM7-limited) approximation to the linear advection terms on cell centers [69, 71].

3.5.2 Incompressible Navier-Stokes solution

The spatiotemporal discretization of the conservative form of incompressible Navier-Stokes Eqs. (15)-(16) reads as

[TABLE]

in which the newest approximation to density $\breve{\bm{\rho}}^{n+1,k+1}$ and the discretization of the convective term $\bm{C}^{n+1,k}$ are computed such that they satisfy consistent mass/momentum transport, which is required to maintain numerical stability for high value of air-water density ratio. We refer the reader to Nangia et al. [65, 58] for more details on consistent mass/momentum transport scheme. Note that $\left(\bm{L}_{\mu}\bm{u}\right)^{n+\frac{1}{2},k+1}=\frac{1}{2}\left[\left(\bm{L}_{\mu}\bm{u}\right)^{n+1,k+1}+\left(\bm{L}_{\mu}\bm{u}\right)^{n}\right]$ is a semi-implicit approximation to the viscous strain rate with $\left(\bm{L}_{\mu}\right)^{n+1}=\nabla_{h}\cdot\left[\mu^{n+1}\left(\nabla_{h}\bm{u}+\nabla_{h}\bm{u}^{T}\right)^{n+1}\right]$. The newest approximation to viscosity $\mu^{n+1,k+1}$ is obtained via the two-stage process described in Eqs. (26) and (27).

In Eq. (30), ${\bm{\rho}}^{n}$ is the face-centered value of the density field $\rho^{\rm full}$ (refer Eq. (27) for the definition of $\rho^{\rm full}$), and the density field $\breve{\bm{\rho}}^{n+1,k+1}$ is obtained by advecting ${\bm{\rho}}^{n}$ discretely. The specific value of the density field $\bm{\wp}$ used to compute the gravitational body force $\bm{\wp}\bm{g}$ is explained next in the context of FD/BP fluid-structure coupling algorithm.

3.5.3 Fluid-structure coupling

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

[TABLE]

in which $\widetilde{\chi}=1-\widetilde{H}^{\text{body}}$, $\widetilde{H}^{\text{body}}$ represents the regularized structure Heaviside function (Eq. (25)) and $\kappa\sim\mathcal{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), $\bm{u}_{\text{b}}$ is the solid body velocity which can be expressed as a sum of translational $\bm{U}_{\text{r}}$ and rotational $\bm{W}_{\text{r}}$ velocities

[TABLE]

The center of mass point $\bm{X}_{\text{COM}}^{n+1,k+1}$ is updated using rigid body translational velocity as

[TABLE]

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

[TABLE]

in which $\text{M}_{\text{b}}$ is the mass, $\bm{I}_{\text{b}}$ is the moment of inertia matrix, $\bm{\mathcal{F}}$ and $\bm{\mathcal{M}}$ are the net hydrodynamic force and torque, respectively, $\text{M}_{\text{b}}\bm{g}$ is the net gravitational force, $\bm{F}_{\text{m}}$ and $\bm{F}_{\text{PTO}}$ are the PTO stiffness and damping force, respectively (see Eqs. (10)-(11)), and $\bm{M}_{\text{m}}$ and $\bm{M}_{\text{PTO}}$ are the corresponding PTO stiffness and damping torque about the center of mass point $\bm{X}_{\text{COM}}$, respectively. The net hydrodynamic force $\bm{\mathcal{F}}$ and torque $\bm{\mathcal{M}}$ on the immersed body is computed by summing the contributions from pressure and viscous forces acting on the body

[TABLE]

The hydrodynamic traction in these two equations is evaluated on Cartesian grid faces that define a stair-step representation of the body on the Eulerian mesh [58]. In Eqs. (37) and (38), $\bm{n}_{\rm face}$ denotes the unit normal to the face (pointing away from the solid and into the fluid), and $\Delta A_{\rm 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 Eq. (30). Therefore we set $\bm{\wp}\bm{g}=\bm{\rho}^{\text{flow}}\bm{g}$ in Eq. (30), which also avoids spurious currents due to large density variation near the fluid-solid interface [73]. More specifically, $\bm{\rho}^{\text{flow}}$ is obtained using Eq. (26) and not through Eq. (27).

3.6 FD/BP validation

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 [74]. 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 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/m3 and $1.2$ kg/m3, respectively, and their respective viscosities are taken to be $10^{-3}$ Pa$\cdot$s and $1.8\times 10^{-5}$ Pa$\cdot$s. Surface tensions effects are neglected as they do not affect the wave and converter dynamics at the scale of the problems considered in this work.

3.6.1 Damped oscillations of a cylinder

We begin by considering a simplified version of the point absorber problem. A planar disk of diameter $D$, mass density $\rho_{\text{s}}$ (or mass $M=\rho_{\text{s}}\pi D^{2}/4$), is attached to a spring of stiffness value $k_{\text{spring}}$, and a mechanical damper of damping coefficient $c_{\text{damper}}/c_{\text{critical}}=\zeta$. Here, $c_{\text{critical}}=2\sqrt{k_{\text{spring}}M}$ is the critical damping coefficient. The value of $\zeta$ determines the behavior of the system: $\zeta<1$ leads to under-damping, $\zeta>1$ leads to over-damping, and $\zeta=1$ results in a critically damped system.

The computational domain is taken to be $\Omega$: $[0,5D]\times[0,10D]$ and the cylinder’s initial center of mass point is released from $(X_{0},Y_{0})=(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 the schematic 5(a). The cylinder is surrounded by air phase. The density ratio between the solid and the fluid phase 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_{\rm spring}=500$ N/m in these simulations. Four test cases corresponding to $\zeta=$ 0.1, 0.5, 1.0, and 1.5 are simulated with each density ratio.

High density ratio $\bm{m^{*}=100}$:

We discretize the domain by a $100\times 200$ uniform grid and use a constant time step size of $\Delta t=2.5\times 10^{-4}$ s to simulate the four damping ratio cases. Fig. 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 phase, 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 $\zeta=0.5$ case. Table 1 tabulates the spatiotemporal resolutions used for these grids. Fig. 6 shows the grid convergence plots for the cylinder displacement and the hydrodynamic force acting in the vertical direction. As observed in the figure, a medium grid resolution of $100\times 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 [75]. Fig. 7 shows the pressure and flow distribution around the cylinder at $t=0.05$ s.

Low density ratios $\bm{m^{*}}$ = 2 and $\bm{m^{*}}$ = 0.8:

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. Fig. 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 $\zeta=$ 0.1 and 0.5 with $m^{*}=2$, and for $\zeta=$ 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 by mechanical springs and dampers. The method remains stable for large density ratio between different phases and it does not suffer from the added mass effects that become prominent when $m^{*}\lessapprox 1$ [76, 77, 78, 79, 80]. We also refer the readers to our prior work on FD/BP method [58] that considers additional low-density ratio cases and validates them against numerical and experimental results from the literature.

3.6.2 Motion under wave excitation

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.

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

[TABLE]

and are in deep water regime $d/\lambda>0.5$.

The diameter of the buoy with relative density $\rho_{\text{s}}/\rho_{\text{w}}=0.9$ is taken to be $D=0.16$ m. The PTO stiffness and damping coefficients are taken as $k_{\text{PTO}}=1995.2$ N/m and $b_{\text{PTO}}=80.64$ N$\cdot$s/m, respectively. The initial submergence depth of the disk is $d_{\text{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 2. Fig. 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 approximately 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.

4 Comparison between Cummins and CFD models

In this section we compare the simulated dynamics of a submerged cylindrical point absorber using Cummins equation-based Simulink model and fully-resolved CFD model implemented in IBAMR. The dynamics of the buoy and the generated PTO power at steady-state are compared using the two models. We take the same buoy setup, PTO coefficients, and wave characteristics of the previous section 3.6.2 to compare the dynamics. Since 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é [82], it can be verified that wave parameters of Sec. 3.6.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.

4.1 One degree of freedom

We first simulate one degree of freedom (1-DOF) buoy using the two models. Fig. 11 compares the heave dynamics of the submerged buoy and the generated PTO power during steady-state. As observed in Fig. 11(a), Cummins model over-predicts the heave amplitude due to the fact that linear potential theory over-estimates the Froude-Krylov forces or the wave excitation loads on the submerged buoy [23, 10, 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.

4.2 Two degrees of freedom

Next, we simulate a 2-DOF buoy that oscillates in heave and surge directions using the two models. Fig. 12 compares the heave and surge dynamics of the submerged buoy and the generated PTO power during steady-state. The amplitude of the heave motion remains almost the same as that of 1-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 1-DOF buoy as shown in Fig. 12(b). This is in agreement with Falnes [83], who also predicts an increase in power absorption efficiency of a point absorber with more than one degree of freedom using theoretical analyses.

Notice that for the 2-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 over-predicts the surge amplitude. The inability of the Cummins model to capture the slow drift phenomenon can be explained as follows. According to Chakrabarti [84], the drift forces estimated by the potential flow theory tend to zero for $ka<0.5$, in which $k=2\pi/\lambda$ is the wave number and $a=\mathcal{H}/2$ is the wave amplitude. In the present simulation, $ka=0.0258$, and consequently, the Cummins equation-based model under-predicts the drift force. We remark that estimating drift forces is an important criteria for designing a robust mooring system for a practical WEC [85]. Therefore, fully-resolved simulations are more reliable for such calculations.

4.3 Three degrees of freedom

Lastly, we analyze the 3-DOF buoy that can undergo heave, surge, and pitching motions. Fig. 13 shows the heave, surge, and pitch dynamics of the buoy. Compared to the 2-DOF buoy dynamics, the heave and surge dynamics (Fig. 13(a) and Fig. 13(b), respectively) show a larger oscillatory behavior when 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 2- and 3-DOF buoys. Fig. 14(a) shows the complete trajectory of the center of mass point, including the slow drift, traced during few steady-state periods of motion ($64-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 3-DOF buoy. The buoy undergoes an average rotation of $5^{\circ}$ at 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 [86]. Consequently, the resulting pitch velocity $\dot{\theta}$, and the viscous drag torque $\bm{M}_{\text{drag},\theta}$ are negligible. The simulated dynamics of a 3-DOF buoy using the potential flow model is therefore quite similar to the prior 2-DOF buoy case.

Fig. 15 shows the power production of the 3-DOF buoy and contrasts it with 2-DOF PTO power generation. As observed in Fig. 15, 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 like ISWEC, the pitch motion of the hull is the primary source of power production. Therefore, fully-resolved WSI framework can be reliably used to resolve all modes of motion of a PA.

4.4 Vortex shedding

The prior sections highlighted the differences in the dynamic response of the buoy in one, two and three degrees of freedom. Here, we analyze the differences in vortex shedding pattern for a buoy having different degrees of freedom. Fig. 16 shows the vorticity production and transport due WSI of 1- and 2-DOF buoys. The vortex production of a 3-DOF buoy is qualitatively similar to the 2-DOF buoy and is not shown.

As observed in Fig. 16, the 1-DOF buoy sheds vortices in an inclined direction, whereas the 2-DOF one sheds it in the vertical direction. This observation can be explained as follows. For the 1-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 2-DOF buoy, Fig. 17 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, as well as the surge velocity of the buoy using the CFD model are also shown in Fig. 17. It can be observed that the surge displacement and the horizontal flow velocity are in phase with each other. Also, 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 heave velocity of the buoy and $v$ component of the flow field is relatively large, as shown in Fig. 17. It can also be observed that the heave displacement and heave velocity of the buoy are out of phase ($\pi/2$ and $\pi$ radians, respectively) with the $v$ component of the flow field. Therefore, for the 2-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. For the 2-DOF (and also 3-DOF) buoy, the surge velocity of the buoy is out-of-phase (by approximately $\pi/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). The difference in the magnitude of surge and $u$ velocities becomes even smaller if $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 2-DOF buoy.

For these simulations, the PTO coefficients are taken to be $k_{\text{PTO}}=1995.2$ N/m and $b_{\text{PTO}}=80.64$ N$\cdot$s/m. Increasing or decreasing these coefficients by a factor of five changed the phase difference between surge and horizontal velocities only slightly, suggesting that the 2- or 3-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 444They 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.

5 Results and discussion

In this section we use the WSI framework to study the conversion efficiency of a 3-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.

• •

Wave height.

The wave absorption or the conversion efficiency of a wave energy converter is defined as the ratio of the mean absorbed power $\overline{P}_{\text{absorbed}}$ to the mean wave power $\overline{P}_{\text{wave}}$. For a regular wave, the wave absorption efficiency can be written as

[TABLE]

Here, $c_{\text{g}}$ is the wave group velocity and $P_{\text{absorbed}}(t)$ is the instantaneous power absorbed by the PTO damper

[TABLE]

5.1 PTO coefficients

Wave energy absorption efficiency of a converter is maximized if it resonates with the incoming waves [83]. 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

[TABLE]

in which $M$ is the mass of the converter, and $A_{33}(\omega)$ and $B_{33}(\omega)$ 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 $\mathcal{T}\in[0.625-1.1]$ s are obtained using AQWA. Fig. 18 shows the variation of frequency-dependent hydrodynamic coefficients in the selected $\mathcal{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: $k_{\text{PTO}}=\omega^{2}\left(M+A_{33}(\omega)\right)$ and $b_{\text{PTO}}=B_{33}(\omega)$.

• •

Optimal control coefficients: $k_{\text{PTO}}=\omega^{2}\left(M+A_{33}(\omega)\right)$ and $b_{\text{PTO}}=80.64$ N$\cdot$ m/s.

Fig. 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 $\mathcal{T}=0.9$ s. Therefore, a higher damping coefficient could be used to narrow down the optimal performance period range and to possibly enhance the efficiency of the converter. This is achieved by using the optimal control damping coefficient value of $b_{\text{PTO}}=80.64$ N$\cdot$ m/s, which is approximately four times larger than the maximum value of $B_{33}(\omega)$ predicted by the linear wave theory for a wave period around $\mathcal{T}=0.9$ s (see Fig. 18(b)). This optimal value of $b_{\text{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 $b_{\text{PTO}}$ value did not enhance the performance of the converter significantly (data not shown). Similar observations were made in Anbarsooz et al. [25]. Notice that extremely large values of $b_{\text{PTO}}$ would lead to over-damping of the system, and therefore should also be avoided. Table 3 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 $\mathcal{H}/\lambda>0.01$. Therefore, linear wave theory does not predict optimal PTO coefficients for steeper waves as it assumes low values of amplitude to wavelength ratio. The sub-optimal values of PTO coefficients predicted by linear wave theory is also reported in [25].

5.2 Buoy density

Next, we consider the effect of buoy density on its conversion performance. We fix the PTO coefficients to $k_{\text{PTO}}=1995.2$ N/m and $b_{\text{PTO}}=80.64$ N$\cdot$ m/s but vary the mass density for this study. Three relative densities of the buoy are considered: $\rho_{\text{s}}/\rho_{\text{w}}=$ 0.5, 0.7, and 0.9. Fig. 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]

Table 4 shows the optimal performing wave period range calculated using Eq. (44) for a fixed value of $k_{\text{PTO}}=1995.2$ N/m and by using lowest and highest values of $A_{33}(\omega)$ from Fig. 18(a). The simulated results and the analytically predicted range suggest that the optimal period range increases with increasing mass density of the buoy. Notice that since $k_{\text{PTO}}=1995.2$ N/m is an optimal reactive PTO stiffness for a relative density 0.9, the efficiency curve is higher for $\rho_{\text{s}}/\rho_{\text{w}}=0.9$ compared to other densities. Fig. 20(b) plots the reactive PTO stiffness for different mass densities of the buoy using Eq. (44): $k_{\rm PTO}=4\pi^{2}\left(M+A_{33}(\omega)\right)/\mathcal{T}^{2}$. Lower PTO stiffness coefficients are obtained for lower mass densities from this relationship. Therefore, lowering $k_{\text{PTO}}$ value from 1995.2 N/m (or in other words using a more optimal $k_{\text{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 4. 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.

5.3 Wave height

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 $\mathcal{H}=$ 0.01, 0.02 and 0.03 m having time periods in the range $\mathcal{T}\in[0.625-1.1]$ s. The absorption efficiency of the buoy for different wave heights is shown in Fig. 21. As observed in the figure, 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 [83, 87, 3]. 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 [3].

6 Conclusions

In this study, we compared the dynamics of a 3-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 FD/BP method over the FD/DLM method is its ability to incorporate external forces and torques in 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 show that the Cummins model over-predicts 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.

Acknowledgements

A.P.S.B. acknowledges support from NSF award OAC 1931368. NSF XSEDE and SDSU Fermi compute resources are particularly acknowledged.

Data availability statement

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