**1. Introduction**

Wave energy conversion has been attracting more research attention due to the rising demand for renewable energy and its advantages including high power density [1], large energy potential, and consistency. Yet, wave energy conversion is still in the precommercial phase. In the last decade, the dynamics and control of WECs were significantly investigated. Starting from linear models, which are widely applied in developing optimal controls [2–5], to nonlinear hydrodynamics (especially nonlinear FK force) [6,7], to highfidelity CNWT simulations of wave energy converters [8,9]. In recent years, significant research effort has been made to improve the economic index of wave energy conversion. Yet, some challenges still exist; one of which is the need for reactive power to optimize energy production. Some conventional optimal control methods [4] aim to achieve resonance between the device and ocean waves to maximize energy production; yet a complex PTO unit would be required to produce reverse power flow, which represents a high expenditure and the generated power quality is poor (large fluctuation) [10,11]. On the other hand, using damping control does not require reactive power, but the performance of the WEC is degraded in terms of power production.

The concept of VSB WECs was recently proposed to address this challenge. The VSB WEC has a'soft' buoy, allowing relatively rapid shape-changing in response to pressure variations. This is unlike the conventional rigid WECs and unlike the concept of variable

**Citation:** Zou, S.; Abdelkhalik, O. A Numerical Simulation of a Variable-Shape Buoy Wave Energy Converter. *J. Mar. Sci. Eng.* **2021**, *9*, 625. https://doi.org/10.3390/ jmse9060625

Academic Editor: Eva Loukogeorgaki

Received: 10 May 2021 Accepted: 31 May 2021 Published: 4 June 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

geometry WECs [12–15] in which a rigid WEC changes shape occasionally and slowly in response to changing wave conditions. The VSB WEC involves nonlinear interaction between the device and the waves (continuous deformation in response to the waves). This interaction can be leveraged to produce more power without the need for reactive power. The motion of the device can be excited due to the deformation as pointed out in reference [16], which employs a physics-based low-fidelity model.

Since the proposed device involves complex and nonlinear hydrodynamics, a tool for high-fidelity simulation is developed in this paper to investigate the behavior of VSB WEC. Due to the nature of this problem, a FSI tool that can handle variable-shape (morphing) structure is needed in the CNWT. This is unlike most existing high-fidelity WEC simulation tools [8], which can handle only rigid body FSI. The FSI capability employed in this paper implements 2-way FSI; the data is transferred between the solid domain (finite element analysis (FEA) model) and the fluid domain (CFD model). As far as the simulation software is concerned, ANSYS® and OpenFOAM are both widely applied in developing CNWTs and solving FSI [17]. The ANSYS 2-way FSI (system coupling simulation) is used in this paper. As detailed in this paper, this type of problem requires large mesh motions; this is accommodated in this paper using the diffusion-based mesh smoothing method.

The main contribution of this paper is to develop a high-fidelity simulation tool that enables the investigation of the recently proposed VSB WECs, which are characterized by their nonlinear hydrodynamics. The tool presented in this paper will enable the validation of other low fidelity models for VSB WECs that are usually used for WEC control design and analysis. The paper is organized as follows: numerical modeling is introduced in Section 2 and mesh generation is introduced in Section 3. Sections 4 and 5 present the simulation results and discussion, respectively. Finally, the conclusion is drawn in Section 6 with the plan for future work.

#### **2. Numerical Modeling**

The details of the numerical simulation approach applied in this study are introduced in this section. A schematic view of the proposed device and the computational domain is presented in Figure 1. As shown in the figure, the dimensions of the computational domain are 80 × 60 × 60 m. The height of the Free Surface Level (FSL) is 40 m measured from the global coordinate system, which is also shown in the figure. Since the proposed device is allowed to change the shape in response to the waves, a significant interaction between the device and waves is expected. By applying ANSYS 2-way FSI, the data transfers between solid and fluid are created on the interface between solid and surrounding fluids (denoted as FSI 1) and the interface between solid and internal gas (denoted as FSI 2). The device initially has a shape of a hollow sphere with a radius of 2 m and a thickness of 0.3 m.

**Figure 1.** Computational domain layout.

#### *2.1. Governing Equations*

Three stratified phases including the environmental air, environmental water and chamber gas are applied in this simulation. Although the environmental air and water are assumed to be incompressible, the chamber gas is assumed to be compressible since the chamber has a large volumetric change due to the interaction between VSB WEC and waves. The conservation equations for unsteady compressible fluid simulation can be written as:

$$\begin{aligned} \frac{\partial \rho}{\partial t} + \nabla \cdot (\rho \vec{v}) &= 0 \\ \frac{\partial}{\partial t} (\rho \vec{v}) + \nabla \cdot (\rho \vec{v} \vec{v}) &= -\nabla p + \nabla \cdot (\overleftarrow{\nabla}) + \rho \vec{g} + \vec{F} \\ \frac{\partial}{\partial t} (\rho E) + \nabla \cdot (\vec{v} (\rho E + p)) &= \nabla \cdot (k\_e \nabla T) \end{aligned} \tag{1}$$

where *ρ* is the fluid density, *v* is the flow velocity vector in the continuity equation. In the momentum equation, *τ* is the stress tensor, and *ρg* represents the gravitational body force. Since in this study, the wave absorption zone is included near the pressure outlet, the source term - *F* is added in the momentum conservation due to the applied numerical beach treatment. Further, in the energy equation, *E* is the energy, and *ke* is the effective thermal conductivity. A transient Reynolds Averaged Navier Stokes (RANS) solver is applied with the governing equations discretized over the fluid domain by using finite volume method. The fluid simulation (CFD) is further coupled with transient structural analysis (FEA).

#### *2.2. Physics Modeling*

The turbulence model applied in this study is the realizable *k* − model, which is considered robust and computational efficient [17]. Compared to the widely applied standard *k* − , the applied model provides more robustness, which is more suitable in this study since the interaction between the device and the waves is highly nonlinear. The transport equations for turbulent kinetic energy (*k*) and dissipation rate () will be solved. Details about the mathematical modeling of the k-epsilon turbulent model can be found in [18,19]. Furthermore, since the chamber gas has a negligible motion with respect to the device, the internal fluid domain is considered to be laminar.

The Volume of Fluid (VOF) method is applied to simulate the multi-phase flow. Environmental air is assumed to be the primary phase, and the environmental water and chamber gas are both assumed to be the secondary phase. To track the interface (e.g., free surface) between different phases, the volume fraction equation is applied:

$$\frac{\partial(\alpha\_q \rho\_q)}{\partial t} + \nabla \cdot (\alpha\_q \rho\_q \vec{v}\_q) = 0 \tag{2}$$

where *ρ<sup>q</sup>* and *α<sup>q</sup>* are the density and volume fraction of the *q*th fluid (phase *q*). The volume fraction *α<sup>q</sup>* varies between 0 and 1 in a cell. For each cell, it is possible that:

$$\begin{aligned} \kappa\_q &= 0 \text{ if the cell is empty of phase } q\\ 0 &< \kappa\_q < 1 \text{ if the cell contains interface} \\ \text{between phase } q \text{ and other phases} \\ \kappa\_q &= 1 \text{ if the cell is full of phase } q \end{aligned} \tag{3}$$

Furthermore, the summation of the volume fraction of all phases should be unity:

$$\sum\_{q=1}^{n} \alpha\_q = 1\tag{4}$$

To further sharpen the free surface between environmental air and water, the interfacial anti-diffusion is enabled by adding a negative diffusion source term in the volume fraction equation [20]:

$$\frac{\partial \mathfrak{a}\_{\eta}}{\partial t} + \nabla \cdot (\vec{v}\_{\eta} \mathfrak{a}\_{\eta}) = -\nabla \cdot (\vec{v}\_{\mathfrak{c}} \mathfrak{a}\_{\eta} (1 - \mathfrak{a}\_{\eta})) \tag{5}$$

$$\vec{v}\_{\mathfrak{c}} = \gamma |\vec{v}\_{\mathfrak{q}}| \frac{\nabla \mathfrak{a}\_{\mathfrak{q}}}{|\nabla \mathfrak{a}\_{\mathfrak{q}}|}$$

where *vc* and *vq* represent the compression velocity normal to the interface and cell velocity, respectively, and *γ* is the compression factor.

Moreover, the large mesh motion due to the highly nonlinear interaction is accounted by diffusion-based smoothing dynamic mesh method. The pressure–velocity coupling problem is solved by the coupled algorithm to achieve aggressive convergence in transient simulation of the proposed challenging model. The density, momentum, turbulence, and energy are spatially discretized based on the second-order upwind scheme. The applied transient formulation is first-order implicit.

#### *2.3. Domain and Boundary Conditions*

The boundary conditions applied in this study are summarized in Figure 2. The waves are generated from the velocity inlet, and the downstream boundary is specified as a pressure outlet. Symmetry boundary conditions are applied for the top, bottom, and sides boundary conditions. The purpose of setting the symmetry boundary conditions is to construct an open sea condition for the simulation where the deep water and infinite air assumptions are applied. As addressed in [21], velocity inlet and slip wall are also considered as appropriate selections for the top, bottom, and sides boundary conditions if the boundaries are located far away from the device (which is the case in this study). The interface between the surrounding fluids and the solid and between the chamber gas and the solid are both defined as adiabatic walls. A mesh will be morphed around the solid to account for the motion and deformation of the device.

**Figure 2.** Boundary conditions applied in the simulation.

Since open sea conditions are applied in this study, the generated waves are assumed to be deep water waves with regard to a water depth 60 m. To prevent the drifting of the device, the flow current velocity is assumed to be zero. The Person–Monskowitz (PM) wave spectrum is applied as the irregular wave spectrum:

$$S\_{PM}(\omega) = \frac{5}{16} \frac{H\_s^2 \omega\_p^4}{\omega^5} e^{-\frac{5\omega\_p^4}{4\omega^4}}\tag{6}$$

where *Hs* is the significant wave height and *ω<sup>p</sup>* = <sup>2</sup>*<sup>π</sup> Tp* is the peak frequency. The irregular wave that has a significant height of *Hs* = 1 m and a peak period of *Tp* = 6 s is applied in this paper given the fact that typical peak period of ocean waves varies from 6 s to 12 s. The frequency range of the wave is from 0.2 rad·s−<sup>1</sup> to 4 rad·s−<sup>1</sup> and the number of frequencies is 20. Moreover, deep water waves are normally multidirectional waves, thus the directional spreading function also needs to be specified. The widely applied cos-2*s* directional spreading function [22–24] is applied in this study:

$$G(\theta) = G(s) \cos^{2s}(\frac{\pi}{2\theta\_s}(\theta - \theta\_{\mathbb{P}})) \tag{7}$$

where *s* is the frequency independent cosine exponent, which is a positive integer. *θ<sup>p</sup>* is the mean wave heading angle, and *θ<sup>s</sup>* is angle spreading from *θp*. In this paper, *s* is selected as 30 for deep water waves as suggested in [25]. The mean wave heading angle is 0◦ and the spreading angle is 90◦, which is found to be a good assumption for directional spreading when *s* > 5 [25]. Additionally, the number of angular components is chosen as 10. The resulting wave spectrum is the multiplication of the frequency spreading and directional spreading:

$$S(\omega, \theta) = S\_{PM}(\omega) G(\theta) \tag{8}$$

The final wave field will be generated based on the wave spectrum by taking superposition of different wave components:

$$\log(r,t) = \sum\_{i=1}^{N\_f} \sum\_{j=1}^{N\_d} a\_{ij} \cos(k\_{x,ij}x + k\_{y,ij}y - \omega\_i t + \phi\_{ij}) \tag{9}$$

where *a* = 2 *<sup>θ</sup>*<sup>+</sup> *δθ* 2 *<sup>θ</sup>*<sup>−</sup> *δθ* 2 *<sup>ω</sup>*<sup>+</sup> *δω* 2 *<sup>ω</sup>*<sup>−</sup> *δω* 2 *S*(*ω*, *θ*)*dωdθ* is the amplitude of each wave component. Further, *kx*,*ij* = *ki* cos(*θj*) and *ky*,*ij* = *ki* sin(*θj*), where *ki* is the wave number of the *i*th wave component (which has a frequency of *ωi*) and *θ<sup>j</sup>* is the *j*th wave heading. Random phase shift *φij* is used in generating the wave field, which is uniformly distributed in [0, 2*π*].

A multiple-directional numerical beach (as shown in Figure 1) is implemented in the vicinity of the downstream and the sides to suppress the wave reflection. A damping sink term is added in the momentum equation for this damping zone:

$$S = -[\mathbb{C}\_1 \rho V + \frac{1}{2}\mathbb{C}\_2 \rho \mid V \mid V] f(z) f(\alpha) \tag{10}$$

where *C*<sup>1</sup> = 116.24631/s and *C*<sup>2</sup> = 3192.6571/m are the linear and quadratic damping resistance, respectively. *V* is the vertical velocity, *z* is the distance measured from FSL, and *x* is the distance along the flow direction. The length of the multi-directional beach in the downstream (beach direction [1, 0, 0]) is specified as twice the wave length (*D* = 2*λ*) [26,27] and the length of the damping zone in the sides (beach direction [0, 1, 0] and [0, −1, 0]) is defined as one wave length (*D* = *λ*).

#### **3. Mesh Generation**

Mesh quality is critical for free surface flows and the proposed challenging model, which involves highly nonlinear interaction between VSB WEC and the waves. Although a structured mesh provides better efficiency and accuracy in the simulation, it has the difficulty in tracking curved or complex geometries. On the other hand, the unstructured mesh is suitable to track complex geometries and change the mesh size locally, while a significant number of cells are required to ensure accuracy [28]. Therefore a hybrid mesh generation is performed in the fluid domain using ANSYS FLUENT meshing application. The structured mesh is applied in the background with a larger element size and the unstructured mesh is applied near the hull of the device and the internal gas domain with

smaller element size. Though hybrid mesh has the advantages of computationally efficient and tracking features accurately, the resulting mesh is non-conformal due to different mesh topologies applied in different zones. Thus, the meshes are matched in the faces that are meshed with different topologies in the calculation.

The mesh of the overall computational domain is presented in Figure 3. As shown in the figure, the unstructured mesh is applied in a refined region, which has a box shape and internal gas domain with a smaller mesh size. The mesh is also locally refined in the regions that require extra accuracy. It is clearly visible in Figure 3 that the mesh in the expected free surface is refined. In addition, as shown in Figure 4, the mesh is also refined in the interface between surrounding fluids and solid and the interface between the internal fluid and the solid to capture the geometry and motion of the device. Moreover, inflation layers are placed in the boundary layer along the surface (with respect to the surrounding fluids) of VSB WEC to capture the rapid velocity change. As shown in Figure 5, the unstructured mesh is also applied in the solid domain to capture the geometry. Since the applied geometry is a simple hollow sphere, no local mesh refinement is used in the solid domain.

**Figure 3.** Mesh generation of the overall computational domain. Faces selected for presentation are bottom, outlet, and a cross-sectional face located at *y* = 30 m.

**Figure 4.** Mesh generation around the body surface with the mesh refinement at the interfaces and the boundary layer.

**Figure 5.** Mesh generation for the solid domain.

To select an appropriate grid resolution, a mesh sensitivity analysis is performed with three different levels of mesh resolutions (coarse, medium, and fine). The device is assumed to be rigid body (fixed-shape) in this analysis, and the details of each grid resolution are summarized in Table 1. The structured mesh is applied in the background, which occupies around 99.8% of the fluid domain, while only 18.4% of the cells of the fluid domain are formed in the background when the medium mesh is applied. On the other hand, nearly 81.6% of the cells (unstructured mesh) are concentrated in less than 0.2% of the fluid volume. The presented mesh arrangement allows a significantly reduced number of cells in the background to accelerate the calculation and keeps the accuracy of the simulation near the wave and structure interaction. The three-dimensional motions of the center of gravity (CoG) of the device (assumed rigid body) versus different mesh resolutions are plotted in Figures 6–8. As shown in the figures, the medium mesh is considered to be an appropriate mesh resolution since it provides close accuracy as the fine mesh while significantly saves the computational cost. Notably, the motion of the device in the y-direction is predicted inaccurately by applying the coarse mesh. The corresponding mesh size applied in different domains of the medium mesh is summarized in the following. The grid resolution applied on the body surfaces is 0.06 m (most refined) and applied in the background is 1.7 m. Moreover, the mesh is refined in the Refined Region (unstructured) with a mesh size of 0.42 m and the free surface is locally refined with a grid resolution of 0.17 m. The resulting minimum cell volume and maximum cell volume of the fluid domain are 2.12 × <sup>10</sup>−<sup>6</sup> <sup>m</sup><sup>3</sup> and 3.25 m3, respectively. The similar mesh size (medium mesh) will be applied for the following simulations of VSB WEC (without assuming fixed-shape).

**Figure 6.** x-directional motion of the center of gravity of the device with different mesh resolutions.

**Figure 7.** y-directional motion of the center of gravity of the device with different mesh resolutions.

**Figure 8.** z-directional motion of the center of gravity of the device with different mesh resolutions.


**Table 1.** The breakdown of different types of mesh applied in mesh sensitivity analysis.

#### **4. Numerical Results**

The system coupling simulations are carried out on Iowa State University's High-Performance Computer (HPC). The number of processors used in the parallel simulation is 36 and the corresponding elapsed time for the simulation of VSB WEC with medium mesh is approximately 96.6 hrs including post-processing. To guarantee the flow current number is close to 1, the time step is selected as 0.01 s and a total of 3000 steps will be solved (simulation end time is 30 s). Although, a more robust analysis (may be conducted in the future) can be conducted by significantly extending the simulation timeseries (e.g., 20 min simulation for varied wave conditions). Giving that the computational cost of this simulation is extremely high, only a short time period (30 s) will be simulated and analyzed. Additionally, the applied simulation end time covers part of the steady-state responses, which provides meaningful results for both transient and steady-state analysis.

#### *4.1. Free Decay Test*

To validate the proposed simulator, a free heave decay test is first performed. The device (assumed rigid body in the decay test) is initially placed 1.3 m away from the equilibrium point with no initial velocity. As shown in Figure 9, the natural decay period of the device is found to be around 2.6 s. Moreover, the energy given to the device can be successfully dissipated by the implemented multi-directional numerical beach, and the device stays at the equilibrium without being excited by reflected waves (absorbed).

**Figure 9.** Heave decay test.

#### *4.2. Irregular Wave Generation*

As described in Section 2.3, multi-directional irregular waves will be applied in this study to simulate open sea conditions. A wave probe is placed between the inlet and the device (at *y* = 30 m) to monitor the free surface elevation (as shown in Figure 1). As suggested in [29], the distance between the inlet and the wave probe is 1*L* = 11 m. The surface elevation measured at the free surface (*αair* = 0.5) by the wave probe is shown in Figure 10.

The wave pattern on the free surface is presented in Figures 11–14 at different time instants when the most significant interaction between waves and the device happens and at the end of the simulation. It is visible in the figure that the wave pattern (which is multidirectional) is significantly affected around the device due to the strong interaction. Although, a wavy water surface is assumed initially, the wave pattern near the outlet and the sides is flattened during the simulation (as shown in Figure 14). To better illustrate the influence of the wave forces on the motion and deformation of VSB WEC, Figures 15 and 16 show the dynamic pressure field around the device at two different time instants. A large dynamic pressure acts on the bottom of the device at *t* = 0.5 s with a peak pressure around 2.89 × 103 Pa, which not only causes a motion, but also a deformation of the hull (compress the device). At *t* = 1.5 s, there is also a significant deformation of the bottom half of the device (mainly at the sides) and peak dynamic pressure is around 0.696 × <sup>10</sup><sup>3</sup> Pa acts on the wet surface in upstream. By presenting these figures, we can claim the challenging interaction between the waves and the device is successfully captured by the proposed CFD model.

**Figure 10.** Free surface elevation measured at the wave probe.

**Figure 11.** Wave pattern on the free surface at *t* = 0.5 s.

**Figure 12.** Wave pattern on the free surface at *t* = 1 s.

**Figure 13.** Wave pattern on the free surface at *t* = 1.5 s.

**Figure 14.** Wave pattern on the free surface at *t* = 30 s.

**Figure 15.** Dynamic pressure field measured around the device at *t* = 0.5 s.

**Figure 16.** Dynamic pressure field measured around the device at *t* = 1.5 s.
