Optimization Study of Marine Energy Harvesting from Vortex-Induced Vibration Using a Response-Surface Method

Abstract

Vortex-induced vibration (VIV) of bluff bodies is one type of flow-induced vibration phenomenon, and the possibility of using it to harvest hydrokinetic energy from marine currents has recently been revealed. To develop an optimal harvester, various parameters such as mass ratio, structural stiffness, and inflow velocity need to be explored, resulting in a large number of test cases. This study primarily aims to examine the validity of a parameter optimization approach to maximize the energy capture efficiency using VIV. The Box–Behnken design response-surface method (RSM-BBD) applied in the present study, for an optimization purpose, allows for us to efficiently explore these parameters, consequently reducing the number of experiments. The proper combinations of these operating variables were then identified in this regard. Within this algorithm, the spring stiffness, the reduced velocity, and the vibrator diameter are set as level factors. Correspondingly, the energy conversion efficiency was taken as the observed value of the target. The predicted results were validated by comparing the optimized parameters to values collected from the literature, as well as to our simulations using a computational-fluid dynamics (CFD) model. Generally, the optimal operating conditions predicted using the response-surface method agreed well with those simulated using our CFD model. The number of experiments was successfully reduced somewhat, and the operating conditions that lead to the highest efficiency of energy harvesting using VIV were determined.

1. Introduction

In recent years, a strong demand has developed for transforming energy supply from fossil fuels towards renewable energy sources, such as solar, wind, and wave. In this context, marine energy, such as ocean currents, offers a significant potential for emissions-free and self-sufficient energy production. Nevertheless, the path to the utilization of currents has been very difficult, as the majority of ocean and river currents are slow, making harvesting of their energy by turbines or watermills challenging [1]. On the other hand, engineered structures with bluff cross-sections, such as circular cylinders and slender structures, may generate alternating or oscillatory lift in uniform flows. As is well known, vortex-induced vibration (VIV) generally occurs to a bluff body immersed in a cross flow, creating shedding vortices in the downstream, with oscillatory unsteady forces, and resulting in an oscillatory motion of the bluff body [2]. Lock-in is an important phenomenon associated with VIV of a bluff body, characterized by a high amplitude of body vibrations. To study the typical determining factors of a single-cylinder VIV, and characterize its lock-in regime, numerous experimental and numerical works have been carried out. For instance, Feng [3] performed wind tunnel experiments to investigate the issue of vortex vibration of a single cylindrical structure and first uncovered the phenomenon of vortex vibration. Khalak and Williamson [4,5] designed an experimental apparatus to study the effect of mass ratio, and the characteristics of the initial branch, the upper branch, and the lower branch, were illustrated using an amplitude graph of the vortex-excited vibrations of a cylindrical structure at different mass ratios. Alternatively, numerical computations are sometimes performed, based on two-dimensional, or more often three-dimensional, assumptions. For instance, Govardhan and Williamson studied VIV and vortex shedding modes using two- and three-dimensional numerical simulations [6,7]. Similarly, Shi et al. [8] studied the influence of corner rounding on the near-wake pattern and vortex structure of a square cylinder flow, by using two- and three-dimensional large eddy simulations. Cravero et al. [9] presented a methodology for the prediction of the recirculation length of the unsteady wake from a bluff body, and the quantitative correlation of higher recirculation length with lower shedding frequency was confirmed for considered applications at different operating conditions. For further works related to synchronization in VIV, see the detailed reviews of [10,11,12]. It is worth noting, that in practical offshore applications, provided that a structure is flexible or flexibly mounted, VIV may cause vibrations, leading to stresses and fatigue damage.

Apart from its drawbacks, on the other hand, recent developments have revealed the possibility of extracting hydrokinetic energy using the VIV phenomena, such as the vortex-induced vibration aquatic clean energy (VIVACE) converter [13]. This converter maximizes the use of VIV energy by enhancing, rather than suppressing, the self-excited vibrations of a cylinder due to its vortex shedding. This concept has been studied extensively since then, consisting of model tests [14], computational fluid dynamics (CFD) [15], and mathematical models [16]. More specifically, a comprehensive physical model study was performed by Narendran et al. [17], to understand the underlying physics behind VIV. Their experimental set-up consisted of a linear generator, which had low mechanical losses, leading to a higher mechanical efficiency. Additionally, they also carried out a detailed parametric investigation, over a wide range of system parameters, such as mass ratio, damping ratio, and Reynolds number, using a new analytical model. Ding et al. [18] numerically studied flow-induced motion and energy harvesting of bluff bodies with different cross-sections, where time-dependent viscous flow solutions were obtained, by combining the incompressible unsteady Reynolds-averaged Navier–Stokes equations, with the Spalart–Allmaras turbulence model. Similarly, energy harvesting via VIV of bluff bodies with different cross-sections in tandem arrangement, was numerically studied by Zhang et al. [19]. Furthermore, Wu et al. [20] studied the VIV of a cylinder at a low Reynolds number, based on CFD simulations, cylinders with eight different gap ratios were considered. The responses of the cylinder system were discussed, following the classification of the flow patterns at different gap ratios. Zhang et al. [21] numerically investigated the effects of the damping-ratio and mass-ratio, on energy harnessing of a square cylinder in flow-induced motion. Their results showed that, as the velocity increased, the number of vortices shed per cycle increased, and the harnessed power increased, without an upper limit. The energy conversion efficiency increased up to the highest value, until the VIV upper branch. In summary, marine current energy may be harvested using one [22] or multiple cylindrical oscillators [23], where smooth cylinders or passive turbulence control could be adopted [24].

To develop an optimal harvester, various parameters such as mass ratio, structural stiffness, and inflow velocity need to be explored, resulting in a large number of test cases. In such contexts, the design-of-experiment (DOE) approach can assist, to perform a systematical and effective analysis [25]. By arranging and interpreting some of the main cases, a DOE approach derives optimal results through a minimum number of analysis cases. As a type of DOE approach, response-surface methodology (RSM) is an optimization method that integrates experimental design and mathematical modeling, to obtain an optimal level of each factor, by fitting the function between the factors and the results in the global range, through experiments at representative local points [26]. Combining optimization methodologies such as an RSM, with numerical models, avoids expensive and time-consuming experiments, while allowing for obtaining optimal conditions using different input combinations. A growing number of studies have adopted an RSM combined with CFD calculations, to identify the dominating factors affecting specific parameters and/or responses, that are vital in field operations. For instance, Seok et al. improved the bow shape of a tanker hull, where CFD simulations were minimized by applying an RSM optimization algorithm [27]. By modeling a dairy building through RSM models and CFD calculations, Yi et al. [28] examined the effects of opening ratio, building length to width ratio, and wind speed, on the air exchange rate. Typically, as an efficient model, with three-level experimental design, Box–Behnken design (BBD) is a practical design method within RSM, that can evaluate the nonlinear relationship between indicators and factors [29]. In comparison to other practical design methods with the same number of elements, BBD has the benefit that it does not necessitate several consecutive trials, which can significantly minimize the number of tests. This optimization method has been successfully applied in various engineering branches [30,31,32,33,34,35]. Nevertheless, its application in marine engineering has been limited. Specifically, to our best knowledge, an optimization study of marine current energy harvesting from VIV using such an RSM is still missing from the literature.

This study aims to demonstrate the applicability of an RSM optimization working scheme, that helps to maximize the energy capture efficiency using VIV. The adopted BBD optimization algorithm allows for us to efficiently explore various design parameters, as well as their influences on the energy capture efficiency. First, a three-dimensional VIV model is constructed, using the finite-volume method. Its numerical setup and model accuracy are then verified by performing grid and time-step sensitivity studies. The optimization algorithm is then introduced, to build models between the energy capture efficiency and its dominant design parameters, i.e., the spring stiffness, mass ratio, and incoming flow velocity, in the present study. The validity of the adopted BBD algorithm is then examined, by comparing its predictions with results obtained from the literature, as well as from our CFD simulations. Finally, the influences of the considered design parameters are systematically and effectively analyzed. This work lays a meaningful case foundation for the combination of a response-surface methodology and computational-fluid dynamics approach to optimize the efficiency of marine energy harvesting from VIV, demonstrating the effectiveness of our coupled design.

2. Optimization Methodology

This section covers the optimization methodology adopted in the present study. A brief introduction of its basic concept is first given. Its validity is then roughly examined, by comparing its prediction with results collected from the literature [36,37,38]. A systematic verification of the adopted optimization methodology is performed in Section 5, where its predictions are validated against our CFD simulations.

2.1. Concept of Response-Surface Method

Response-surface methodology (RSM) is a computational and scientific technique for modeling and analyzing situations, that considers various factors impacting the desired response and attempts to maximize the result [39]:

$$y=f\left(x_1,x_2\right)+ϵ$$

(1)

Here, variables $x_1$ and $x_2$ maximize the yield of a process y. $ϵ$ is the error variable. Its efficiency depends on the accuracy of y at various points throughout the response surface. Box–Behnken design (BBD) is a three-level RSM design, fitting a second-order response-surface model [40,41]:

$$y={\beta }_0+\sum _{i=1}^k{\beta }_i{x}_i+\sum _{i=1}^k{\beta }_{ii}{x}_i^2+\sum _{i<j}\sum {\beta }_{ij}{x}_i{x}_j+ϵ$$

(2)

where i and j are the linear and quadratic coefficients, respectively. $x_i$ and $x_j$ are the uncoded independent variables, and k is the number of factors studied and optimized in the experiment. Its parameters $\beta$, are estimated by using the least-square approach, where a surface-response analysis is performed on the linked surface.

A BBD requires only the low, middle, and high levels for each design variable. As shown in Figure 1, the levels of factors are at the midpoints of the edges (blue dots), and in the center is the center point (orange dot). It does not account for corner points (i.e., points with design factors at their highest levels), because it tends to avoid extreme settings for experiments, which can potentially lead to higher prediction variance near the vertices [42]. The three-dimensional coordinates of each design point in the figure represent the three test levels of each test point, and the test itself requires that the three levels are equally distributed throughout the domain. $N=2k\times \left(k-1\right)+C_0$, where k represents the number of factors and $C_0$ represents the number of central test point, is used to estimate the test error [22].

2.2. Feasibility Verification of Response-Surface Method

To verify the feasibility of the adopted optimization algorithm, before starting our simulations, various operational variables were collected from the literature [38], and their proper combinations are validated against our predictions using BBD in this subsection. The considered verification case consisted of double cylinders with various diameter ratios and spacing ratios. The associated energy capture efficiency was simulated using CFD, for various reduced velocities. In this regard, the cylinders’ diameters, inlet velocities, and spacing ratios between cylinders, were selected as level factors. The associated energy capture efficiencies were treated as target observation, as shown in Figure 2. Herein, a three-dimensional (3D) surface graph allows the correlation between a response variable and two operating variables to be viewed. Furthermore, to determine the influence degrees of the two operating variables on the response value, one may observe the slope of a response surface. The steeper the slope, the more significant the corresponding operating variable is, making it intuitive to find an optimal range of operating variables.

We see that the chosen operating variables (i.e., spacing ratio, diameter ratio, and reduced velocity) all have a significant influence on the response variable (i.e., efficiency). There is a relationship between each of the two operating variables, and their optimal combinations are identified, leading to the maximum response value. Specifically, as shown in Figure 2a, at small diameter ratios, the efficiency increases with the increase in spacing ratio. However, this is not the case for large diameter ratios, where the efficiency decreases with the increase in spacing ratio. The influence of diameter ratio on the efficiency, has more or less the same trend of that of the spacing ratio. The interactions of spacing ratio and diameter ratio with reduced velocity, are given in Figure 2b and Figure 2c, respectively. We see that the efficiencies increase first, and then decrease, as the inlet velocity increases. This is because there was a locking interval of incoming velocity on the vibration energy of the oscillator, and within the locking interval, the motion system reached the maximum energy conversion efficiency.

Table 1. Energy capture efficiency collected from [38] under various operating conditions, together with the optimal condition (OP) predicted using BBD (the data was from [38]).

Run	1	2	3	4	5	6	7	8	9	10	11	12	13	14	15	16	17	OP
$\eta$ [%]	29.87	18.33	35.41	21.66	29.87	27.08	14.16	29.87	25.25	16.25	22.54	9.58	28.33	29.87	16.25	29.87	26.66	36.10

lists the energy capture efficiencies collected from the literature under various operating conditions. Meanwhile, the predicted maximum efficiency using the optimal operating conditions is also included, where the spacing ratio is 5.78, diameter ratio is 0.99, and the reduced velocity is 5.23. In such an optimal condition, the corresponding efficiency is 36.10%.

Overall, for all the considered operating variables, the surface curves are steep, indicating that the chosen operating variables within this validation study all have a pronounced influence on the target response variable. From which, one may possibly identify the optimum setting of each operating parameter, leading to the maximum energy capture efficiency. We conclude that the adopted BBD algorithm is feasible in the VIV research area, evidenced by the fact that these peaks in the response surfaces, represent the most significant response values among the factors and levels under study. The maximum efficiency was successfully predicted using the optimal conditions. It must be emphasized that the predicted maximum efficiency was only roughly evaluated against the used input data rather than a systematic comparison against CFD simulations. To verify the correct and reliable application of the present method for VIV, the latter will be addressed in Section 5.

3. Numerical Methods of Modeling Fluid–Structure Interaction

To model the VIV of a single cylinder, the numerical methods used for modeling fluid–structure interaction are addressed in this section. The fluid dynamics is solved using the finite-volume method, while the equation of motion of a rigid body is treated as a classical spring-damper system. They are solved within Star-CCM+, in an implicitly coupled manner.

3.1. Computational Fluid Dynamics

For an incompressible and viscous flow, i.e., water in the present study, it is governed by the continuity and the momentum conservation laws:

$$\nabla ·\mathbf{v}=0$$

(3)

$$\frac{\partial \mathbf{v}}{\partial t}+\left(\mathbf{v}·\nabla \right)\mathbf{v}=\nu {\nabla }^2\mathbf{v}-\frac{1}{\rho }\nabla p+\mathbf{g}$$

(4)

where $\mathbf{v}$ is the fluid velocity filed vector, $\nu$ and $\rho$ are the kinematic viscosity and the density field, respectively, p is the pressure, and $\mathbf{g}$ is the gravity vector. The force acting on the body is given as follows [43]:

$$\mathbf{F}={\mathbf{F}}_p+{\mathbf{F}}_v={\int }_S\mathbf{n}pdS+{\int }_S{\rho }^2{\nu }^2\mathbf{n}·\left[\nabla \mathbf{v}+{\left(\nabla \mathbf{v}\right)}^T\right]dS$$

(5)

where subscripts p and v denote the pressure and viscous component of the fluid-induced force, respectively. $\mathbf{n}$ stands for surface normal vector of the body surface.

To solve turbulent flows, the two-equation model K-epsilon [44], was used in this study. This model has been applied to vortex shedding in previous work [45] Consequently, the velocity and pressure may be decomposed into time-averaged values $\overline{\mathbf{v}}$, $\overline{p}$, and turbulent fluctuations ${\mathbf{v}}^{\prime }$, $p^{\prime }$, the time average of which is zero:

$$\begin{array}{c}\hfill \mathbf{v}=\overline{\mathbf{v}}+{\mathbf{v}}^{\prime }\\ \hfill p=\overline{p}+p^{\prime }\end{array}$$

(6)

Inserting Equation (6) into Equations (3) and (4), results in the Reynolds Averaged Navier–Stokes (RANS) equations for incompressible flow. The solution domain is discretized into a finite number of control volumes, which may be of arbitrary shape. Second-order temporal and spatial discretization schemes were adopted. The convergence criterion, which is the residual of velocity components and pressure in the momentum equations, was set to ${10}^{-4}$.

3.2. Rigid Body Motions

In the Star-CCM+ software, the motion of a rigid body was solved using the dynamic fluid body interaction (DFBI) tool. Considering the transverse motion of an elastically mounted cylinder only, its equation of motion can be simplified as the classical mass-spring-damper oscillator model in one degree of freedom (DOF):

$$m\ddot{z}+c\dot{z}+kz=F_z$$

(7)

where z is the position of the rigid body in the transverse direction, and $F_z$ is the corresponding fluid force in that direction. m, c, and k are the structural system mass, damping, and stiffness, respectively. After a dimensionless analysis, the following parameters are introduced:

$$\begin{array}{cccc}\hfill m^{*}& =\frac{m}{m_f}=\frac{4m}{\rho \pi d^{2}l}\hfill & \hfill \zeta & =\frac{c}{2\sqrt{km}}\hfill \\ \hfill f_n& =\frac{1}{2\pi }\sqrt{\frac{k}{m}}\hfill & \hfill v_r& =\frac{v}{d{f}_n}\hfill \end{array}$$

(8)

where $m^{*}$ is the mass ratio and $m_f$ is the mass of fluid occupied by the structure. d and l are the body’s diameter and length, respectively. $\zeta$ is the damping ratio and $f_n$ is the natural frequency, and $v_r$ is the reduced velocity. The drag and lift coefficients are defined as:

$$c_d=\frac{2F_x}{\rho v^{2}dl}{c}_l=\frac{2F_z}{\rho v^{2}dl}$$

(9)

where v is the free stream velocity, and $F_x$ and $F_z$ are the drag and lift forces exerted by fluid flow loads. It is worth noting that a global earth fixed coordinate system and body local coordinate system were used to predict the rigid body motion. The fluid force (Equation (5)) acting on the rigid body was calculated on the earth fixed coordinate system, while the rigid body equation was solved in the local coordinate system and again converted back to the earth coordinate system. The new position of the body was updated using mesh deformation, with respect to the global earth fixed axes.

4. Test Case Description and Solution Verification

The computational domain is sketched in Figure 3a, where its dimensions were chosen to exclude the blockage effect [46]. Grids near the cylinder and its downstream were refined, to better capture the shedding vertices. In our simulations, the incoming velocity direction was parallel to the x-axis, only vibration responses in the transverse direction were considered (i.e., parallel to the z-axis). The diameter of the cylinder is $d=0.08$ m, and its length is $l=0.44$ m, which is identical to the width of the computational domain. The velocity inlet is at 15d away from the center of cylinder at the left end, and the pressure outlet is at 30d away at the right end. The top and bottom sides were set as slip walls, which are 15d away from the cylinder. A non-slip wall condition was applied for the cylinder. Prism layers were used around the cylinder, to make sure the viscous sublayer was directly resolved (i.e., $y^{+}<1$). A grading ratio of 1.1 was used and a total of 30 prism layers were generated.

On three successively finer grids, using a refinement factor of 1.26, we performed a grid-spacing and a time-step convergence study to estimate numerical errors. For the three grids, with number of cells of $G_1$ = 413,742, $G_2$ = 784,448, and $G_3$ = 1,510,836, three time-step sizes were selected, so that the Courant–Friedrichs–Lewy number satisfies that $CFL=v\Delta t/\Delta x\le 1$. Defined as functions of the inlet velocity, our time-step sizes were $\Delta t_1={\left(455v\right)}^{-1}$, $\Delta t_2={\left(910v\right)}^{-1}$, and $\Delta t_1={\left(1820v\right)}^{-1}$, respectively. In the solution verification study, the mass ratio was $m^{*}=1.36$, and the effect of structural damping was excluded, i.e., $\zeta =0.0$, to maximize the vibrations amplitude. The considered structural stiffness was of $k=400$ N/m and the reduced velocity was $v_r=4.92$, resulting in a Reynolds number of $Re=49,674$ and a Strouhal number of $St=0.20$. Five simulations were performed using the designed meshes and time steps, and assessments of numerical errors were performed for numerically predicted amplitude ratios $A^{*}=A/d$. The convergence ratio R is defined as:

$$R=\frac{{\varphi }_3-{\varphi }_2}{{\varphi }_2-{\varphi }_1}$$

(10)

where ${\varphi }_1$, ${\varphi }_2$, and ${\varphi }_3$ are solutions obtained on the coarse, medium, and fine grids, respectively. For results that converge monotonically, i.e., $0<R<1$, the Richardson extrapolation was then applied, to derive the estimated numerical error, ${\delta }_{RE}$, and the observed order of accuracy, $p_{RE}$.

$$p_{RE}=\frac{ln\left(1/R\right)}{ln\kappa },{\delta }_{RE}=\frac{{\varphi }_3-{\varphi }_2}{1-{\kappa }^{p_{RE}}}$$

(11)

Following [47], a ratio of convergence P, was applied, as a measure to define the deviation of solutions from the asymptotic range:

$$P=\frac{p_{RE}}{p_{th}}$$

(12)

where $p_{th}$ is the theoretical order of convergence limited by the discretization schemes, and $p_{th}=2$ was adopted here for a second-order accuracy method. The numerical error, ${\delta }_D$, and the numerical benchmark result, ${\varphi }_{\infty }$, were obtained as follows:

$${\delta }_D=P{\delta }_{RE},{\varphi }_{\infty }={\varphi }_3-{\delta }_D$$

(13)

Following [48], our observed order of accuracy was limited to $p_{RE}\ge 0.5$ and correspondingly, a P value of close to 1 indicates grids in the asymptotic region. The uncertainty, $U_D$, was then calculated as follows: [47]:

$$U_D=\left\{\begin{array}{cc}\left(2.45-0.85P\right)\left|{\delta }_{RE}\right|,& 0<P\le 1\\ \left(16.4P-14.8\right)\left|{\delta }_{RE}\right|,& P>1\end{array}\right.$$

(14)

Table 2. Estimated errors and uncertainties of response ratio ($A^{*}$) based on various time-step and grid-spacing sizes.

Property	${\mathit{A}}_1^{*}$ [-]	${\mathit{A}}_2^{*}$ [-]	${\mathit{A}}_3^{*}$ [-]	R [-]	${\mathit{A}}_{\mathit{\infty }}^{*}$ [-]	${\mathit{ϵ}}_1$ [%]	${\mathit{ϵ}}_2$ [%]	${\mathit{ϵ}}_3$ [%]	U [%]
$\Delta t$	0.5625	0.6043	0.6075	0.077	0.6080	−7.48	−0.61	−0.08	0.41
$\Delta x$	0.5950	0.6075	0.6087	0.096	0.6089	−2.29	−0.76	−0.04	0.16

lists the estimated errors and uncertainties based on $A^{*}$. The time-step study was performed using three different sizes of $\Delta t_1$, $\Delta t_2$, and $\Delta t_3$, based on the medium grid $G_2$. Large discrepancies were observed between the results obtained from $\Delta t_1$ and those obtained from $\Delta t_2$ and $\Delta t_3$ (see Figure 4). Then, the grid-spacing study was performed, using the selected time-step $\Delta t_3$ for three grid sizes $G_1$, $G_2$, and $G_3$. We see that decreasing the grid-spacing size from $G_2$ to $G_3$ did not improve the results significantly. Both grids can provide favorable results, whose numerical errors are below 1%. We chose the medium grid size $G_2$, with the fine time step $\Delta t_3$, for the subsequent simulations.

5. Simulation Design Using a Response-Surface Method

This section presents simulations designed using the BBD algorithm. It starts with an introduction to our simulation design. Following this is a statistical analysis of the associated results. It ends with the verification of the optimum operating parameters, using CFD simulations.

5.1. Simulation Design and Statistical Analysis

To design the simulation matrix, the associated design factors and levels needed to be identified. Based on historical data of previous studies, the energy capture efficiency of VIV was dominated by its incoming flow velocity, structural stiffness, and oscillator mass, which were selected as operating variables in the present study. A three-factor, three-level response-surface test was then conducted using these operating variables, and, correspondingly, the energy capture efficiency, $\eta$, of VIV was treated as the response value.

Table 3. Level of operating variables selected for the Box–Behnken design.

Variables	Symbols	Level −1	Level 0	Level 1
Velocity [m/s]	A	0.55	0.65	0.75
Stiffness [N/m]	B	300	400	500
Mass [kg]	C	2.60	3.00	3.40

lists these independent variables, as well as their levels for the Box–Behnken design. It is noted that for these operating variables, good energy capture efficiencies were observed for the selected operating ranges [49,50,51].

For the adopted BBD algorithm with three factors and three levels, a total of 13 sets of simulations were proposed.

Table 4. Designed test cases, together with the estimated results using our CFD simulations.

Run	Factor A Velocity [m/s]	Factor B Stiffness [N/m]	Factor C Mass [kg]	Response $\mathit{\eta }$ Efficiency [-]
1	0.55	300	3.0	0.140
2	0.75	300	3.0	0.127
3	0.55	500	3.0	0.166
4	0.75	500	3.0	0.113
5	0.55	400	2.6	0.139
6	0.75	400	2.6	0.095
7	0.55	400	3.4	0.128
8	0.75	400	3.4	0.112
9	0.65	300	2.6	0.120
10	0.65	500	2.6	0.129
11	0.65	300	3.4	0.124
12	0.65	500	3.4	0.116
13	0.65	400	3.0	0.124

lists the sets of operating conditions used in each test and the simulated results using our CFD model. It must be noted, that as the structural damping was excluded, the conventional approach [52] of energy conversion efficiency is not valid for the present study. Inspired by the field of wave energy utilization [53], alternatively, the energy capture efficiency can be represented, to some extent, using the oscillator’s amplitude and the flow velocity:

$$\eta =\frac{A}{Tv}$$

(15)

where T is the structure’s oscillatory period. As the kinetic energy of the structure is a function of its oscillatory amplitude and period, and the energy of fluid over the area swept by the structure is dominated by its velocity.

To assess the significance and suitability of the coupled RSM prediction with CFD simulations, a statistical analysis is carried out in this section. According to the adopted design model, an empirical relationship between the CFD simulations and the BBD input variables was fitted, resulting in a second-order polynomial equation, with the interaction terms (i.e., Equation (2)):

$$\begin{array}{cc}\hfill \eta =& 0.124-0.01575\times \mathrm{A}+0.001625\times \mathrm{B}-0.000375\times \mathrm{C}-0.01\times \mathrm{AB}+0.007\times \mathrm{AC}\hfill \\ \hfill & -0.00425\times \mathrm{BC}+0.004375{\times \mathrm{A}}^2+0.008125{\times \mathrm{A}}^2-0.009875\times {\mathrm{C}}^2\hfill \end{array}$$

(16)

In terms of these coefficients, the magnitudes of A, B, and C reflect the magnitude and direction of their influences on the energy conversion efficiency of VIV. The estimated results, together with the analysis of variance (ANOVA), are summarized in

Table 5. Estimated results, based on the BBD algorithm, together with the analysis of variance.

Source	Coefficient	Sum of Squares	DOF	Mean Square	f-Value	p-Value
Model	0.1240	0.0034	9	0.0004	55.50	<0.0001
A	−0.0158	0.0020	1	0.0020	290.9	<0.0001
B	0.0016	0.0000	1	0.0000	3.100	0.1218
C	−0.0004	1.1 × 10${}^{-6}$	1	1.1 × 10${}^{-6}$	0.165	0.6968
AB	−0.0100	0.0004	1	0.0004	58.64	0.0001
AC	0.0070	0.0002	1	0.0002	28.73	0.0011
BC	−0.0042	0.0001	1	0.0001	10.59	0.0140
A${}^2$	0.0044	0.0001	1	0.0001	11.81	0.0109
B${}^2$	0.0081	0.0003	1	0.0003	40.75	0.0004
C${}^2$	−0.0099	0.0004	1	0.0004	60.19	0.0001
Residual	-	0.0000	7	6.8 × 10${}^{-6}$	-	-
Adj. R${}^2$	0.9684	-	-	-	-	-
Pre. R${}^2$	0.7789	-	-	-	-	-

.

Herein, the model row shows how much variation in the response is explained by the model, along with the overall model test. In total, nine individual terms were tested independently in this study. The residual indicates how much variation in the response is still unexplained. Sum of squares represents the difference between the overall average and the amount of variation. The total degrees of freedom (DOF) are the amount of estimated parameters. Additionally, the sum of squares divided by the DOF is the mean square. Specifically, the f-value is used for comparing the source’s mean square to the residual mean square, and the large f-value of 55.50 in Table 5 implies that the adopted model is valid. The p-value denotes the probability of seeing the observed f-value if the null hypothesis is true, i.e., there are no factor effects. A small probability value calls for rejection of the null hypothesis. From Table 5, we see that the extremely small p-value gives confidence that the probability of occurrence of such an f-value due to noise is less than 0.01%.

Generally, a p-value less than 0.05 gives confidence that the associated model term is statistically significant, and a value greater than 0.10 indicates the model term is not significant. In this study, A, AB, AC, BC, A${}^2$, B${}^2$ and C${}^2$ are significant model terms. We see that for the primary term, only factor A is significant. The influence of each independent variable on the composite score became successively more pronounced for the factors C, B, and A. Furthermore, we see that all the secondary terms were significant for the considered case, indicating that the interaction terms of incoming velocity, structural stiffness, and mass, all had significant effects on the target response. The adjusted coefficient of determination (adj. R${}^2$) and the predicted coefficient of determination (pre. R${}^2$), were estimated as 0.9684 and 0.7789, respectively. Correspondingly, a quadratic model (Equation (16)) was suitable to describe the relationship between the factors and responses for the present study. We may conclude that the accuracy of the obtained model is high enough and it provides a good fit between our simulations and the predicted values.

5.2. Verification of Optimum Parameters Using CFD Simulations

Based on the designed test cases, an optimal operating condition was given, using statistical analysis in Section 5.1. The optimum operating variables were identified as $v=0.55$ m/s, $k=500$ N/m, and $m=2.798$ kg, resulting in an energy capture ratio of $\eta =0.166$. To validate the optimum ratio predicted by the BBD algorithm, we performed a CFD simulation using these optimum operating variables. Figure 5 plots the amplitude ratio, A${}^{*}$, together with the lift coefficient, C${}_l$, in the time and frequency domains. The associated distributions of velocity and vorticity for the considered cylinder at its optimal operating condition are given in Figure 6. The simulated energy capture ratio is $\eta =0.169$, which agrees well with the BBD predicted one, with a discrepancy of 1.77%. Although the adopted optimization algorithm slightly underpredicts the energy capture ratio, the combination of a response-surface method with computational fluid dynamics shows its capability to identity the optimum operating conditions for marine energy harvesting from vortex-induced vibrations.

In addition, Figure 7 plots a comparison of the predicted energy capture ratio by the BBD algorithm, versus the actual energy capture ratio simulated using CFD. We observe that all the points of the predicted and simulated response values are very close to the ${45}^{\circ }$ line, indicating that the developed model successfully captures the correlation between our CFD simulations and the BBD-predicted values.

6. Concluding Remarks

We coupled a response-surface method with computational fluid dynamics, to perform an optimization study of marine energy harvesting from vortex-induced vibration. The validity of the adopted response-surface method for optimization purposes was first examined using the inputs collected from literature, and a reasonable agreement was achieved. To numerically simulate the cylindrical vortex-induced vibration, a numerical model, based on the finite-volume method, was then built, and its solutions were verified through time-step and grid-spacing studies. Afterwards, a set of CFD simulations were carried out, whose operating variables were recommended by the design tool. The number of simulation tests, under the influence of multiple operating variables, were significantly reduced.

The final optimal operating variables were then estimated using the second-order polynomial equation, which represents the empirical relationship between the CFD simulations and the optimization input variables. To obtain the highest energy capture ratio, the optimum conditions of the independent operating variables, i.e., inlet velocity, stiffness, and mass, were predicted. The predicted ratio from the optimization algorithm was eventually validated against that calculated using our CFD model. Our study demonstrated the great potential of the combination of a response-surface methodology and computational-fluid dynamics, to optimize the efficiency of marine energy harvesting from VIV, proving the effectiveness of the coupled design.