A Numerical Simulation of a Variable-Shape Buoy Wave Energy Converter

Wave energy converters (WECs) usually require reactive power for increased levels of energy conversion, resulting in the need for more complex power take-off (PTO) units, compared to WECs that do not require reactive power. A WEC without reactive power produces much less energy, though. The concept of Variable Shape Buoy Wave Energy Converters (VSB WECs) is proposed to allow continuous shape-change aiming at eliminating the need for reactive power, while converting power at a high level. The proposed concept involves complex and nonlinear interactions between the device and the waves. This paper presents a Computational Fluid Dynamics (CFD) tool that is set up to simulate VSB WECs, using the ANSYS 2-way fluid–structure interaction (FSI) tool. The dynamic behavior of a VSB WEC is simulated in this CFD-based Numerical Wave Tank (CNWT), in open sea conditions. The simulation results show that the tested device undergoes a significant deformation in response to the incoming waves, before it reaches a steady-state behavior. This is in agreement with a low-fidelity dynamic model developed in earlier work. The resulting motion is significantly different from the motion of a rigid body WEC. The difference in the motion can be leveraged for better energy capture without the need for reactive power.


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][3][4][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 geometry WECs [12][13][14][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.

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.

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: 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 k e 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).

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: where ρ q and α q are the density and volume fraction of the qth fluid (phase q). The volume fraction α q varies between 0 and 1 in a cell. For each cell, it is possible that: Furthermore, the summation of the volume fraction of all phases should be unity: 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]: where v c and v q 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.

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. 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: where H s is the significant wave height and ω p = 2π T p is the peak frequency. The irregular wave that has a significant height of H s = 1 m and a peak period of T p = 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 −1 to 4 rad·s −1 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-2s directional spreading function [22][23][24] is applied in this study: where s is the frequency independent cosine exponent, which is a positive integer. θ p is the mean wave heading angle, and θ s 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: The final wave field will be generated based on the wave spectrum by taking superposition of different wave components: where a = 2 S(ω, θ)dωdθ is the amplitude of each wave component. Further, k x,ij = k i cos(θ j ) and k y,ij = k i sin(θ j ), where k i is the wave number of the ith wave component (which has a frequency of ω i ) and θ j is the jth 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: where C 1 = 116.24631/s and C 2 = 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 = λ).

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.   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 × 10 −6 m 3 and 3.25 m 3 , respectively. The similar mesh size (medium mesh) will be applied for the following simulations of VSB WEC (without assuming fixed-shape).

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.

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

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 1L = 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 × 10 3 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 × 10 3 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.

VSB WEC Deformation
Since the proposed device allows significant deformation in response to the waves, it is interesting to study the geometry variation in free motion. As shown in Figure 17, although the geometry of the device varies significantly from 0 s to 1.5 s, the shape of the device at 1.5 s indicates a trend of recovering to the original shape (at 0 s). Thus, a periodic change of the shape is expected, which is also reasonable since the gas chamber is a restoring mechanism. This periodic phenomenon can be better explained in the following figures. Figures 18 and 19 show the vertical motion of two nodes where the U=pper Node denotes the node defined originally at (40, 30,42) m and the Lower Node denotes the node defined at (40, 30,38). Figure 18 shows the individual motion of two nodes with respect to their original location. For a conventional fixed-shape buoy WEC (FSB WEC), the motions of these two nodes are expected to be identical since rigid body motion is assumed. Although, in this study, their motion patterns are different. The relative vertical distance between these two nodes is plotted in Figure 19. The original distance between these two nodes is 4 m, which will be kept as a constant for a FSB WEC. While the vertical distance of these two nodes for the VSB WEC changes periodically and finally reaches a steady-state deformation (approximately when t > 15 s). The natural period of steady-state deformation is around 1.96 s, which is almost one third of the wave peak period. The relation between the natural period of the shape deformation and the wave period can be more rigorously studied in the future by simulating the dynamic behavior of VSB WEC with varied wave conditions, which may benefit future optimization of the proposed device.

VSB WEC Motion
The motion of VSB WEC is presented in detail in this section and it will be compared with FSB WEC motion under the same wave condition to highlight the difference in motions introduced by changing the shape. The motion difference can be later controlled and utilized to produce more wave power. Additionally, the mass and dimensions of FSB WEC are the same as VSB WEC (undeformed shape). The motions of the center of gravity of VSB WEC and FSB WEC in x, y, and z directions are shown from Figures 20-22, respectively. The motion of VSB WEC significantly differs from FSB WEC in three dimensions both in terms of phase and magnitude. The difference in vertical motion between the two devices is presented in Figure 23. The maximum difference is around 0.26 m, and the maximum z-directional motion of FSB WEC is around 0.64 m, which indicates a maximum 40.6% change in the heave motion caused by the introduction of geometry variation. This large difference represents a room for improving the performance of wave energy conversion by introducing appropriate control.

Discussion
In this paper, a 3D numerical simulation architecture is introduced to simulate the fluid-structure interaction of the proposed VSB WEC. The simulation results show the proposed device has a significant shape-changing due to the highly nonlinear interaction between waves and the device. The resulting motion of VSB WEC also significantly differs from the motion of a conventional FSB WEC with the same mass and dimensions. It is noted that the design of VSB WEC can be further optimized in terms of energy extraction since this paper focuses on presenting a framework of the high-fidelity simulation of this device. The applied design is simple in terms of geometry and mechanism. For instance, a possible design presented in [16] is a VSB WEC that has a rigid body part (cylindrical shape) and a shape-changing part (truncated conical shape). A control valve is included to control the pressure oscillation in the gas chamber. As found in [16], a restoring mechanism (gas chamber in this study) is required to keep a steady-state deformation of the device in order to harvest more energy. The presented framework is applicable to simulate the performance of other designs of VSB WEC.
Hyperelastic material is applied in this study since the proposed device is required to be 'soft'. The material applied in the simulation obeys Neo-Hookean hyperelasticity with an initial shear modulus of 0.1 Mpa and an initial bulk modulus of 2000 Mpa, and the density of the material is 1000 kg·m −3 . The corresponding strain-stress curve of the applied material is shown in Figure 24. Unlike vulcanized rubber [30] (initial shear modulus is 0.41 Mpa and initial bulk modulus is 414.5 Mpa), which only allows small deformation, the material used in this study is 'softer', which allows more deformation. In addition, unlike neoprene rubber (initial shear modulus is 0.027 Mpa and initial bulk modulus is 13.86 Mpa), which is too soft such that the deformation is too large to keep a stable simulation. More advanced dynamic mesh techniques need to be applied to address this challenging mesh motion. The study of the effect of different materials is interesting, though beyond the scope of this paper, therefore a material that has the hyperelastic properties between vulcanized rubber and neoprene rubber is applied such that the simulation is stable with a considerable deformation of the device.

Conclusions
A high-fidelity CFD-based numerical wave tank simulation for a VSB WEC is presented in this paper. This numerical tool is demonstrated to be able to simulate a VSB spherical WEC. This highly nonlinear interaction between the device and the waves is simulated using a 2-way FSI technique. Open sea conditions are applied in this study by assuming infinite air, deep water, and multi-directional waves. The numerical results in this paper capture the shape deformation in response to the varying surface pressure from the simulated ocean waves. It was shown that the VSB spherical WEC exhibits a transient response period before it reaches a steady-state motion and deformation. It was shown that this resulting motion of the VSB WEC significantly differs from that of a FSB WEC of the same shape as that of the non-deformed VSB WEC. This difference in response characteristics (difference in motion trajectories) is key for investigating the advantage of the VSB WEC power production over the FSB WEC.