Principal Parameters Analysis of the Double-Elastic-Constrained Flapping Hydrofoil for Tidal Current Energy Extraction

Abstract

Taking the rigid NACA0012 airfoil as the object, the key structural parameters of the spring–mass system that govern the dynamics of the double-elastic-constrained flapping hydrofoil are numerically studied in this paper. A two-dimensional numerical model, based on the CFD software FINE/Marine, is established to investigate the influence of the spring stiffness coefficient, frequency ratio, and damping coefficient on the motion and performance of the flapping hydrofoil. This study demonstrates that when the structural parameters are adequately adjusted, the power factor exceeding 1.0 has been achieved, and the corresponding efficiency is up to 37.8%. Moreover, this system can start and work within a wide range of damping coefficients. However, the hydraulic efficiency and power coefficient are sensitive to the change in damping coefficient, so it is very necessary to design an appropriate power output. Lastly, the most obvious parameter affecting the energy acquisition performance is the spring stiffness coefficients. Frequency ratios in the two directions have little influence on the peak value of the power coefficient, but they will cause the change of damping coefficients of the peak point. The key structural parameters studied in this paper provide a useful guideline for an optimized design of this interesting system through searching for the best performance.

1. Introduction

Tidal stream energy is an important marine renewable energy, which is widely distributed offshore and easy to extract. In the past decade, tidal stream energy converters have become a major focus for renewable energy research, with a number of turbine farms now in planning and development stages. At present, the majority of existing designs for tidal energy devices include three main types, namely horizontal-axis turbine-based energy converters [1,2,3], vertical-axis turbine-based energy converters [4,5], and flapping-hydrofoil energy converters. Compared with the first two conventional turbines, the flapping hydrofoil has several prominent features of relatively low tip speed, no centrifugal stress related to the rotating blade, and rectangular cross section of flow [6]. Subsequently, it has the advantages of low noise, friendliness to the marine biological environment, robust structure, and easy arrangement in shallow water.

Flapping-hydrofoil energy-harvesting devices have been developed over several decades. The application of flapping hydrofoils to extract energy from uniform flows was first proposed by McKinney and DeLaurier in 1981 [7]. They realized the heave and pitch motions of the flapping hydrofoil through a set of mechanisms and controlled the motion amplitude in two degrees of freedom to obtain the performance. Although these turbines were first identified to extract energy from the wind in McKinney’s work, they have mostly been developed as hydrokinetic turbines during the last decade, which extract energy from rivers or tidal currents. It is interesting to note that the flapping-foil turbine concept has been proven to be competitive with the horizontal-axis and vertical-axis turbine technologies, with efficiencies exceeding 40% [8,9]. Assuming that two degrees of freedom of pitch and heave are allowed, flapping-hydrofoil flow-energy harvesters are generally classified by their activation mode into three categories [9,10]: (i) a fully constrained system with motions that are fully prescribed in both heave and pitch motions, (ii) a semi-passive system with forced pitch and induced heave motions, and (iii) a fully passive system with free heave and pitch motions, fully determined by the fluid–structure interaction.

In the first mode, there are roughly two ways to realize the fully constrained flapping-hydrofoil system. One is the form of McKinney and DeLaurier’ work [7], that is, a linkage mechanism was used to realize the prescribed heave and pitch motions. The heave and pitch motions were coupled with each other. Kinsey et al. [11] also adopted this form in 2011, and the pitch motion of the flapping hydrofoil was coupled to the cyclic heave motion through four link mechanisms. However, Kinsey’s work effectively yielded a one-degree-of-freedom system, while flapping motions in McKinney’s work were both controlled by the execution unit. The other way is to use two separate actuators to realize heave and pitch motion control, respectively, such as the work by Simpson et al. [12] and Kim et al. [13]. Simpson et al. [12] conducted a comprehensive parameters study of this kind of foil system in a small water tank, including motion amplitude, St number, aspect ratio, etc.

With its simplicity in modeling, a foil undergoing fully constrained motions, which mostly are prescribed as sinusoidal, has been extensively studied. Many numerical studies have shown that energy-extraction performance is strongly related to pivot location, oscillation frequency and heave amplitude, and the relative phase of pitch and heave motions [8,14,15,16]. With the research development and deepening on fully constrained flapping hydrofoils, many scholars have found that flapping-hydrofoil devices essentially harvest the flow energy through the heave motion [17,18]. This means that the pitch motion produces or incurs only modest inputs/outputs of energy in the mean [19], which suggests the concept of a further simplified, semi-active mode.

In the second semi-active mode, the pitch motion of the foil is controlled by the actuator unit, while the induced heave motion is controlled by elastic constraints. These systems require controlling/actuating the pitch motion. When the pitch angle changes periodically between positive and negative, the foil system can move and work continuously in the heave motion. Thereby, energy input is needed to activate the pitch motion, whereas energy harvesting is achieved through the resulting heave motion generated by fluid dynamic lifting forces. Deng et al. [20] called this kind of foil system a semi-passive flapping hydrofoil, while Teng et al. [21] called it a semi-active flapping hydrofoil. From its structural characteristics, it could also be called a single-elastic-constrained flapping hydrofoil. In 2009, Hisanori et al. [22] completed the experimental study of this single-elastic-constrained flapping hydrofoil earlier, and the system’s energy-extraction efficiency reached 37%. In 2015, Deng et al. [20] discussed the influence of mass inertial effects on energy-extraction efficiency by a numerical method. In addition, they proposed that since the heave amplitude of this kind of flapping hydrofoil was uncertain, dimensionless coefficients of power should be used to measure its performance rather than efficiency. Then, Teng et al. [21] discussed the effects of non-sinusoidal pitch motion on the energy-extraction performance of this kind of foil in 2016. In the same year, Griffith et al. [23] analyzed the characteristics of a single-elastic-constrained flapping hydrofoil with different operating parameters by simulating a range of geometries, from the circular cylinder with an aspect ratio of 1.0 to elliptical cylinders up to an aspect ratio of 6.0. In addition to the studies on kinematic parameters, the research on fully constrained and semi-active foil systems tends to consider the influence of the realistic environment, where the flapping foils are deployed [24].

In the third type, the foil is connected with a torsional spring and a linear spring, passively oscillating in two degrees of freedom. The pitch motion and heave motion of this kind of foil are controlled by elastic constraints, which is a coupled-mode flutter phenomenon in its work. Compared with systems operating in the first two modes, fully passive systems do not involve any actuator, which simplifies the structure and reduces maintenance costs. Boudreau et al. [25] called it a fully passive flapping hydrofoil, while Young et al. [8] also defined the first kind of fully constrained system as fully passive. In order to avoid ambiguity, this paper considers its constraint characteristics and calls it a double-elastic-constrained flapping hydrofoil. Although this coupled-mode flutter phenomenon of airfoil has been thoroughly studied in the aeroelasticity field’s community, mainly with the purpose of determining the instability onset and ways to mitigate it [26,27], few studies have analyzed the large-amplitude motions that can result from it. In the use of tidal current energy, we prefer to have a more intense and greater amplitude of flutter, and higher efficiency will be obtained.

Compared to fully constrained or semi-active systems, there has been much less work that investigates the double-elastic-constrained flapping-hydrofoil system. The reader could be referred to the recent reviews on this topic by Xiao and Zhu [9] and Young et al. [10]. This is partly because such a complex system involving coupled fluid–structure interaction and nonlinear coupling in two-degrees-of-freedom flapping motions has many parameters, and the experimental or numerical simulation is very challenging. In 2009, by using a Navier–Stokes model, Peng and Zhu [28] investigated the dynamics of a foil mounted on a rotational spring and a heave damper in a uniform incoming flow. In addition, they suggested that the pivot location is an important factor to the energy-harvesting capacity. Later, Zhu [29] suggested that double-elastic constraint systems may be much more adaptable to various real applications than the other two systems. Then, Veilleux and Dumas [30], Jiang et al. [31,32], and Wang et al. [33] conducted a series of numerical studies to investigate the influences of natural frequency, pivot location, spring stiffness, damping, and 3D effects on the system’s energy-extraction performance. Those promising numerical results were later verified by Boudreau [34], who performed the very first experimental study on the subject. They obtained a power coefficient of 0.86 and an efficiency of 31%. In addition, they also discussed the influence of multiple design parameters on performance and found that good performance can be attained over a wide range of flow and structural parameters, which is important from the practical point of view. Using a NACA0015 foil Re ≈ 6 × 10⁴, Duarte [35,36] experimentally studied the same type of system and found that the pivot axis should be located at least 0.29 chord length from the leading edge to achieve continuous operation. They also revealed that the system can perform much better by increasing the heaving natural frequency. More recently, in 2020, a follow-up numerical study was then conducted by them on this system, which obtained a much higher efficiency of 53.8% [26]. In 2022, R. Fernandez-Feria innovatively introduced the flexible foil into the double-elastic constraint system and studied its flutter instability [37].

In summary, despite the high hydraulic efficiencies of about 40% achieved by the first mode fully constrained flapping foil, it has the disadvantage of requiring very complex and costly constraining mechanisms. The alternative solution of a semi-active mode would have the advantage of being less complex from a technological point of view, while proving to be as efficient as the first type of fully constrained flapping foil devices. However, reports from the first full-scale prototype of a semi-passive flapping-foil turbine also highlighted prohibitively high maintenance costs related to this semi-active system. Compared with systems operating in the first two modes, double-elastic-constrained flapping-hydrofoil systems do not involve any actuator, which simplifies the structure and reduces maintenance costs. Therefore, more and more attention has been achieved both on numerical and experimental research of the double-elastic constrained flapping-hydrofoil turbine. The pitch motion and heave motion of this kind of foil both result from the interaction of the elastically supported foil and the fluid flow, thereby leaving the designer of such a turbine with only indirect control over the foil motions by adjusting the structural parameters that govern the foil dynamics. So, as discussed, although the research on the double-elastic constrained flapping hydrofoil has a period of time, its physical mechanism is still not clear in some respects, and its working characteristics and optimal working parameters still need a lot of data to confirm.

For this purpose, on the one hand, many of the above studies on foil mechanism parameters are completed under some specific working conditions, and the coverage of working condition parameters is not comprehensive, so it is difficult to obtain quantitative and universal laws. The systematic and parametric analysis will be conducted to correlate the influencing factors to the energy-conversion efficiency in this paper; on the other hand, either too small or too large damping is unfavorable for the system to obtain energy. Thereby, as an important parameter, the damping coefficient will be optimized to achieve the highest energy-conversion efficiency in this work. In addition, the access to the velocity and vorticity fields—motion and force characteristics as functions of space and time—could allow us to identify the underlying thrust production mechanisms more easily. In such a context, the spring–mass system parameters of the double-elastic-constrained flapping hydrofoil are studied by the numerical method in this paper. As a follow-up work of Boudreau and Dumas [25], the present investigation is an extension work in a more comprehensive range of parameters and achieves the highest energy-conversion efficiency by configuring reasonable damping. The double-elastic constrained flapping-foil turbine concept and the metrics that characterize its performance are presented in Section 2 and Section 3. Then, the numerical methodology and validation are presented in Section 4. The impact of various structural parameters affecting the turbine blade dynamics is analyzed in Section 5. In the range of parameters studied in this paper, the power coefficient of this turbine exceeding 1.0 has been achieved, and the corresponding efficiency is up to 37.8%. In addition, the variation of power coefficient with the damping coefficient is also obtained. The present study aims to provide a useful guideline on the design and operation of this interesting system through searching for the system’s best performance.

2. Motion Description

The double-elastic constrained flapping-hydrofoil turbine concept considered in this study is similar to the one described in the work of Veilleux and Dumas [30] and Matthieu Boudreau [25]. As shown in Figure 1, it consists of a rigid foil elastically supported by springs in heave and pitch, as well as by two dampers. The heave damping (D_h) models the energy conversion into electricity by an electric generator. The total pitch damping coefficient (D_θ) is corresponding to the sum of some undesired viscous pitch damping, as well as some coulomb friction (dry friction) in pitch, due to the inherent presence of friction in a real experimental setup. The foil is free to pitch (z direction, right-handed helix rule) about a pitch axis located at a distance c_EA from the leading edge and to heave in the y direction as schematically shown. The motion is not possible in any other direction nor about any other axis. In addition, gravity acts in the span direction, or z direction, hence playing no role in the blade dynamics.

The tidal current along the negative x direction flows to the foil system at the speed V_∞. In case there are small force disturbances, such as uneven incoming flow or pulsating force caused by turbulence, the foil will leave the equilibrium position and oscillate in the two directions of heave and pitch. This oscillation is caused by the flutter instability, i.e., the coalescence flutter phenomenon. The energy to stimulate the coalescence flutter comes from the kinetic energy of the tidal current, so the tidal current energy can be obtained through this phenomenon. Generally, we obtain energy in the heave direction, because according to the existing research, the energy harvested by the pitch motion is negligible compared with the heaving motion [36].

When the foil deflects under the action of disturbing force, it is subjected to the action of lift (i.e., y direction force) and head rolling moment at the same time and begins to move. According to the dynamic relationship between the components of the system, the elastic-constrained foil could be considered as a linear mass–spring system, and its heave and pitch motions are governed by the second-order oscillator equations as follows.

$$\left(\begin{array}{cc}{m}_h& S\cos \theta \\ S\cos \theta & I_{\theta }\end{array}\right)\left(\begin{array}{c}\ddot{y}\\ \ddot{\theta }\end{array}\right)+\left(\begin{array}{cc}{D}_h& 0\\ 0& D_{\theta }\end{array}\right)\left(\begin{array}{c}\dot{y}\\ \dot{\theta }\end{array}\right)+\left(\begin{array}{cc}{k}_h& 0\\ 0& k_{\theta }\end{array}\right)\left(\begin{array}{c}y\\ \theta \end{array}\right)=\left(\begin{array}{c}{F}_y-S{\dot{\theta }}^2\sin \theta \\ M_z\end{array}\right)$$

(1)

Here, the superscript (·) denotes differentiation with respect to time. In addition, referring to the definition method of Boudreau [25], the definitions of all the parameters involved in Equation (1) are provided in

Table 1. List of the parameters involved in the equations of motion.

Symbol	Units	Definition
V∞	m/s	Freestream velocity
cEA	m	Distance between the elastic axis and the leading edge of the airfoil
y	m	Instantaneous displacement in heaving motion (upward positive in Figure 1)
θ	rad	Pitch angle (clockwise negative in Figure 1)
$m_h$	Kg	Mass of all the components undergoing the heaving motion (including foil)
$D_h$	Ns/m	Total heave damping
$k_h$	N/m	Heave spring stiffness
$I_{\theta }$	Kg m2	Moment of inertia about the pitching axis
$D_{\theta }$	Nm s/rad	Total pitch damping
$k_{\theta }$	Nm/rad	Pitch spring stiffness
$S$	Kg m	mθxθ Static moment (mass of the components only undergoing the pitch motion times xθ)
$m_{\theta }$	Kg	Mass of all the components undergoing the pitch motion (including foil)
$x_{\theta }$	m	Distance between the pitch axis and the center of pitching mass (defined positive when the pitch axis is upstream of the center of mass)
$F_y$	N	Hydrodynamic force component in the heave (y) direction
$M_z$	N m	Hydrodynamic moment about the pitch axis
b	m	Blade span length

.

3. Methodology

As mentioned in the previous section, since the coupled-mode flutter of the system results from the interaction of the double-elastic-constrained foil and the fluid flow, the designer of such a turbine could only indirectly control the foil motions by adjusting the structural parameters that govern the system’s dynamics. Therefore, in order to obtain the energy efficiently, on the one hand, the system parameters in the motion equation should be reasonably configured to make the hydrofoil flutter continuously; on the other hand, the appropriate damping coefficient (representing power output) should also be discussed to achieve the highest energy-conversion efficiency. Here, on the basis of a series of research results conducted by Boudreau et al. in 2018 [24] and 2020 [25], we selected the dimensionless parameters recognized by most researchers and further analyzed and simplified the key parameters according to the existing research results.

Since this paper mainly focuses on the energy acquisition efficiency of the system, for a given incoming flow, this paper redefined 8 dimensionless structural parameters which explicitly play a role in the turbine dynamics. Their definitions and significance are shown in

Table 2. Non-dimensional parameters of double-elastic-constrained foil system.

Non-Dimensional Parameters	Definition
$k_h^{*}=\frac{k_h}{\rho V_{\infty }^{2}b}$	Heave spring stiffness coefficient, represents the magnitude of heave amplitude
$k_{\theta }^{*}=\frac{k_{\theta }}{\rho V_{\infty }^{2}b{c}^2}$	Pitch spring stiffness coefficient, represents the magnitude of pitch amplitude
${\omega }_h^{*}=\frac{c}{V_{\infty }}\sqrt{\frac{k_h}{m_h}}$	Frequency ratio in the heave direction, represents the natural frequency in the heave direction
${\omega }_{\theta }^{*}=\frac{c}{V_{\infty }}\sqrt{\frac{k_{\theta }}{I_{\theta }}}$	Frequency ratio in the pitch direction, represents the natural frequency in the pitch direction
${\zeta }_h=\frac{D_h}{2\sqrt{m_h{k}_h}}$	Total linear heave damping coefficient
${\zeta }_{\theta }=\frac{D_{\theta }}{2\sqrt{I_{\theta }{k}_{\theta }}}$	Total linear pitch damping coefficient
$S^{*}=\frac{S}{m_{h}c}$	Mathematically static imbalance arm length
$c_{EA}^{*}=\frac{c_{EA}}{c}$	Position of pitch axis

. In this table, b is the span of the flapping hydrofoil, c is the chord length of the foil, and ρ represents the density of the fluid. In particular, S is the static moment, and it is equal to the product of the mass of the components undergoing the pitch motion with the distance between the center of mass and the pitch axis, expressed as m_θx_θ. However, from Equation (1), there is no necessary connection between S and m_θ. In order not to introduce new non-dimensional parameters, m_h is used instead of m_θ in the dimensionless parameter construction of S. S* is expressed as the mathematically static imbalance arm length. As the key parameters affecting the energy acquisition efficiency of this system are very rich, including not only the parameters listed in Table 2 but also airfoil, three-dimensional effect, Re number, non-uniform incoming flow, oblique flow, etc., a comprehensive parametric study is thus beyond the scope of the present work. Considering the focus of this paper, the parameter setting of this problem is simplified as follows.

In this paper, a rigid and symmetrical NACA0012 airfoil is preliminarily considered for analysis. Because NACA0012 airfoil has a broader research basis, and it was also chosen to allow a direct comparison of the current results with previous studies conducted by the authors on a semi-passive flapping-foil turbine. At the same time, because the flutter of the hydrofoil system is always accompanied by large separation, which is similar to the separation of the slender cylinder, its three-dimensional (3D) effect is inevitable. The 2D analysis method is adopted in this paper because of the following reasons. On the one hand, in the experimental cases referred to in this paper (Boudreau [25]) and in most real engineering projects, the flapping foil is equipped with end plates at both ends or is close to the wall, which reduces the 3D effect and is similar to the 2D flow. Therefore, 2D analysis in this paper could be able to provide useful insights for real (3D) cases; on the other hand, this paper focuses on the influence of system parameters on energy efficiency. The 3D effect is not the focus of this paper. In addition, a large number of parameter analyses will be carried out in this study, and 3D analysis is strictly out of reach. However, the study of the 3D effect should be and will be further studied in the follow-up work.

The turbine operates at a Reynolds number (Re = V_∞c/γ) of 2 × 10⁵ based on the freestream velocity and the chord length. Such a large Reynolds number also ensures that the boundary layers are turbulent, hence making the use of a turbulence model in the fully turbulent mode adequate, as stated in Section 4.1. According to the previous research, the tidal current energy obtained by the flapping hydrofoil mainly comes from the heave motion. Therefore, in order to simplify the analysis, this paper does not consider the influence of pitch damping and sets the pitch damping to zero. Meanwhile, from the research of Boudreau et al. [25] and Goyaniuk et al. [38], s* is very small, so this paper sets s* as zero. That is, the mass center of the foil coincides with the position of the rotating axis. The position of the rotating axis of the head heave motion is selected at 35% chord length from the leading edge [35]. Since the turbine dynamics is analyzed via 2D numerical simulations, forces per unit span are obtained, and the span length is therefore considered to be equal to one (b = 1). Through the above simplification, there are only five main parameters affecting the energy acquisition efficiency of the double-elastic-constrained flapping hydrofoil, including k_h*, k_θ*, ω_h*, ω_θ*, and ζ_h. The final selected baseline case and the parameters describing this baseline are shown in

Table 3. Baseline case description.

Airfoil	Dimension	${\mathit{k}}_{\mathit{h}}^{*}$	${\mathit{k}}_{\mathit{\theta }}^{*}$	${\mathit{\omega }}_{\mathit{h}}^{*}$	${\mathit{\omega }}_{\mathit{\theta }}^{*}$	${\mathit{\zeta }}_{\mathit{h}}$	${\mathit{\zeta }}_{\mathit{\theta }}$	${\mathit{S}}^{*}$	${\mathit{c}}_{\mathit{E}\mathit{A}}^{*}$	b
NACA0012	2D	2.0	0.08	0.707	0.894	0.0~6.0	0.0	0.0	0.35	1

. Considering that the damping ratio directly affects the amount of power obtained, there should be an optimal damping. So, it is necessary to analyze the influence of the damping ratio (0.0~6.0) for each parameter configuration.

Because the flutter of the double-elastic constraint system is formed by the interaction of the elastically supported foil and the fluid flow, there will be multiple fluctuations in the amplitude of heave or pitch in a cycle. Therefore, in order to measure the motion characteristics and performance of the flapping hydrofoil, this paper takes the difference between the peak value and valley value in one cycle as the heave or pitch displacement in one stroke and defines its maximum heave amplitude and pitch amplitude as half of the average value in N strokes. The maximum heave amplitude ΔY_max and pitch amplitude Δθ_max are shown in Equation (2). Taking heaving motion as an example, a stroke refers to the process that the heave coordinate starts from the zero point and then passes through the zero point twice; heave amplitude in one stroke Δ Y_i refers to the distance (in y direction) from the lowest edge to the uppermost edge in a single movement.

$$\begin{gathered} \Delta Y_{\max }=\frac{1}{2N}{\displaystyle \sum _i\Delta Y_i} \\ \Delta {\theta }_{\max }=\frac{1}{2N}{\displaystyle \sum _i\Delta {\theta }_i} \end{gathered}$$

(2)

Here, ΔY_i is the heave amplitude of the pitch axis location in one stroke, m; Δθ_i is the pitch amplitude in one stroke, radian; N represents the number of cycles and is a set of positive integers.

The average mechanical power provided to produce electricity by a double-elastic-constrained turbine is defined through F_yv_y in Equation (3). Therefore, two metrics characterizing the energy extraction, namely the hydraulic efficiency (η) and the average power coefficient (c_P), both based on the cycle-averaged power dissipated in the eddy-current brake which models the energy extraction, are also defined in Equations (4) and (5). The hydraulic efficiency (η) of the turbine expresses the ratio between the total power harvested and the hydraulic power available in the cross section of the flow swept by the foil, whereas the average power coefficient (c_P) of total power harvested by the turbine is normalized in terms of the projected surface of the foil instead, where the unknown maximum heave amplitude (ΔY_max) is replaced by the known chord length (c).

$$\overline{Power}=\frac{1}{t_1-t_0}{\displaystyle {\int }_{t_0}^{t_1}{F}_y{v}_{y}dt}$$

(3)

$$\eta =\frac{\overline{Power}}{\frac{1}{2}\rho V_{\infty }^{3}A}$$

(4)

$$c_P=\frac{\overline{Power}}{\frac{1}{2}\rho v_{\infty }^{3}bc}=\frac{2\Delta Y_{\max }}{c}\eta$$

(5)

Here, t₀ and t₁ are the start time and end time of total N cycles, respectively; A is defined as swept area ($A=2b\Delta Y_{\max }$), the product between the span length b and the maximum heave amplitude ΔY_max during its motion; V_∞ is freestream velocity.

4. Numerical Method

A 2D numerical model for simulating the double-elastic-constrained flapping-hydrofoil turbine has been implemented on the CFD software FINE/Marine (a software package of NUMECA company, EURANUS solver is developed by the European Space Agency). This solver adopts internal implicit iteration within a time step iteration to ensure strong flow/motion coupling. The main features of the model—namely the fluid-solid coupling strategy, mesh and boundary conditions, turbulence modeling, and validation—will be discussed in this section.

4.1. Governing Equations

Fluid–structure interaction is a multiphysics problem that involves a two-way interaction between the fluid and the structure. Therefore, the hydrodynamic performance of this double-elastic-constrained flapping-hydrofoil turbine can be obtained by solving the fluid hydrodynamic equation and the system motion equations simultaneously. The equation of motion of the foil (Equation (1)) was introduced above. The integral incompressible viscous fluid dynamics equation is adopted, and the gravity factor is not considered. The continuity equation and momentum equation can be written as Equations (6) and (7) separately. In addition, the Menter’s k-ω shear-stress transport turbulent model was used [39].

$$\frac{\partial }{\partial t}{\displaystyle \underset{\Omega }{\int }\rho d\Omega }+{\displaystyle \underset{S}{\oint }\rho \overrightarrow{v}\cdot d\overrightarrow{S}}=0$$

(6)

$$\frac{\partial }{\partial t}{\displaystyle \underset{\Omega }{\int }\rho v_{i}d\Omega }+{\displaystyle \underset{S}{\oint }\rho v_i\overrightarrow{v}\cdot d\overrightarrow{S}}={\displaystyle \underset{S}{\oint }{\tau }_{ij}d{S}_j}-{\displaystyle \underset{S}{\oint }pd\overrightarrow{S}}$$

(7)

Here, Ω is the element volume; $\overrightarrow{v}$ is the flow velocity; $v_i$ is velocity component; ${\tau }_{ij}$ is the sum of viscous stress and Reynolds stress and p is the pressure; $\rho$ is the density of water; and $\overrightarrow{S}$ is the area vector of the element’s surface.

4.2. Mesh and Method

The computational domain and mesh structures are shown in Figure 2. A hexahedron/polyhedron hybrid mesh with local mesh refinement was chosen for the computational domain. In order to accurately simulate the separation on the foil’s surface and the propagation of the wake vortex, further mesh refinement is carried out on the foil’s surface and in the wake region, respectively. The enlarged detail of the near-wall grid could also be shown in Figure 2, in which the maximum mesh size is less than 1/2⁸ of the chord length. In addition, the refined near-wall mesh could ensure the y plus value is about 1.0.

The setting method of boundary conditions is shown in Figure 3. The simulation is carried out in a computational domain of 20c × 40c, with a circular subdomain and a hollow rectangular subdomain. The circular subdomain around the foil is 6c in diameter and can rotate to adapt to the foil’s head-pitching motion. In addition, when the inner circular subdomain was rotational, the outer rectangular one was relatively stationary. The sliding mesh model was coupled at the interface between the two subdomains. When the foil heaves in the y direction, the two subdomains heave synchronously. In the calculation domain, the velocity conditions are given at the right side (inflow), and the pressure conditions are given at the left side (outflow). The velocity and pressure gradient conditions are given on upper and lower boundaries, which are set as the far-field boundary, as shown in Figure 3.

4.3. Validation

In order to validate the numerical method used in the current work, two parts of validation are carried out in this part: the verification of grid independence and simulation verification of selected benchmark conditions in Boudreau’s experiments [25].

Firstly, a two-dimensional NACA0012 airfoil with a chord length of 0.2m is taken as the object to verify the grid independence. The results show that the grid size and area near the airfoil have a great impact on the performance simulation of the foil. Combined with our previous study, four sets of grids with different mesh refinement schemes were designed under an appropriate working condition. According to different grid growth rates and refinement areas, four mesh refinement schemes are shown in Figure 4. The number of grid layer refinements each time is 2, 4, 8, and 12 in different refinement schemes, respectively. In the selected working condition, k_h* = 2.0, ω_h* = 0.707, k_θ* = 0.08, and ω_θ* = 0.894, and the comparison results of foil performance with different meshes are shown in Figure 5. The incoming flow velocity is defined as 1 m/s, and the Re number is about 2 × 10⁵. The time step Δt = 0.001 s is selected in the simulation. The position of the rotation axis coincides with the center of gravity and is set at 0.35c.

As can be seen from Figure 5, when the mesh density near the airfoil is lower and the mesh refinement area is small (e.g., Mesh 0), the foil hydraulic efficiency and power coefficient decrease significantly. With the increase of the mesh refinement, the impact of the mesh on the foil performance becomes smaller. With Mesh 2 and Mesh 3, the foil performance and motion curves in Figure 5 almost coincide. Therefore, in order to give consideration to both computational efficiency and mesh refinement, Mesh 2’s scheme is selected for subsequent simulation in this paper.

To provide further validation of the numerical method, a study of the double-elastic-constrained flapping hydrofoil was conducted to directly compare with the previous work of Boudreau [24]. The airfoil used in the experiment is NACA0015, the chord length is 50 mm, and the experiment is completed in a water tank with an incoming flow velocity of 0.38 m/s. Because Boudreau’s experiment adds end plates at both ends of the blade to reduce the 3D effect, the 2D numerical method is used to verify this example in this paper. The efficiency (η), heave amplitude (ΔY_max), pitch angle (Δθ_max), and reduced frequency (f*) are compared, respectively, as shown in Figure 6. It can be seen that although there are still some errors between the two results, the overall trend is consistent, and the coincidence degree of some working conditions is relatively high. It is speculated that there are two main reasons for this problem. On the one hand, the flutter of this system is always accompanied by a large angle of attack and complex separation, and the accurate simulation of the separation has always been the difficulty of this kind of problem. On the other hand, the attitude of the oscillating foil changes rapidly and violently at any time in the experiment, and it is difficult to accurately evaluate the influence of the water attached to the oscillating foil in the numerical method. In addition, this is obviously not a simple problem. It requires a lot of research, which is not what this paper can carry out. Follow-up research can be carried out in the future.

In conclusion, the numerical method used in this paper is feasible and reliable for the hydrodynamic performance analysis of the double-elastic-constrained flapping hydrofoil.

5. Results and Analysis

In this section, the impact of various structural parameters affecting the flapping-hydrofoil turbine dynamics is analyzed. The effects of each structural parameter are investigated individually by keeping all the other structural parameters constant with their baseline value. Firstly, the working characteristics of the flapping hydrofoil under a wide range of damping coefficient ζ_h are analyzed. Secondly, the flow field and force at the working point with high efficiency and power coefficient are analyzed; then, the effects of the dimensionless spring stiffness coefficient (k_h*, k_θ*) and the frequency ratio (ω_h*, ω_θ*) on the performance of the flapping hydrofoil are analyzed.

5.1. Influence of Damping Coefficient on Energy Acquisition Efficiency

In the double-elastic-constrained flapping-hydrofoil system for tidal current energy extraction, damping is used to control the load. Obviously, when the damping is zero, that is, there is no load, the flapping-hydrofoil system cannot extract energy from the tidal current; when the load is very large and the system cannot bear it, it cannot work normally, so it cannot extract energy either. Therefore, there should be an optimal load for tidal current energy extraction. This part will analyze the performance of the flapping hydrofoil under different damping coefficients, establish the performance curve, and find the working characteristics of the flapping hydrofoil.

Firstly, in the baseline case presented in Table 3, dynamic behaviors of the hydrofoil are studied under a wide range of damping coefficients from 0 to 0.55. The characteristics (including hydraulic efficiency η, power coefficient c_p, heave amplitude ΔY_max, and pitch amplitude Δθ_max) with different damping coefficients are given in Figure 7. As can be seen from Figure 7, on the whole, this flapping-hydrofoil system can start and work regularly within a wide range of damping coefficients, but the hydraulic efficiency and power coefficient are sensitive to the change of the damping coefficient. With the increase in the damping coefficient, the hydraulic efficiency increases rapidly at first and then remains almost unchanged. The highest efficiency can reach about 50%, which is very close to the highest efficiency of 53.8% achieved by Matthieu Boudreau in 2020 [26]; a noticeable feature of the power coefficient is that it does not change monotonically with the change of the damping coefficient, but there is a peak. Under the condition of a small damping coefficient, the power coefficient increases rapidly with the increase of the damping coefficient, but with the further increase of the damping coefficient, the heave amplitude decreases suddenly, and the power coefficient decreases slowly. Under the baseline case working condition, the maximum power coefficient of 0.75 is reached with the damping coefficient of 0.15.

Evaluating of energy extraction has been discussed in some papers [25,29,36], and the difference between hydraulic efficiency η and power coefficient c_p is due to the different definitions of the expression. That is, η is the total harvested power normalized by the power available in the flow window harvested by the device, while c_p is the total harvested power normalized by a characteristic power based on the chord length of the foil. The hydraulic efficiency η defined in Equation (4) in this paper takes the heave amplitude as the characteristic length. Because the heave amplitude decreases with the increase of the damping coefficient, the efficiency does not show a decreasing trend under a large damping coefficient but continues to increase, which is inconsistent with the actual energy obtained by the flapping hydrofoil; as shown in Equation (5), the power coefficient c_p takes the fixed chord length c as the characteristic length, which represents the actual energy-extraction performance. When the damping coefficient is too small or too large, its energy-extraction performance is not good. Therefore, it is more reasonable to use the power coefficient c_p to evaluate the energy-extraction performance of this system. For example, a very high value of hydraulic efficiency η combined with small amplitudes of motion would result in a low value of c_p. In this scenario, the device would be very efficient, but very little power would be harvested from the flow. However, considering that the hydraulic efficiency η is also of certain significance to the working characteristics analysis, and both of the two metrics need to be optimized in order to obtain an interesting turbine, the hydraulic efficiency η will still be listed and discussed in this paper.

In Figure 7, the pitch amplitude of the foil keeps a large value under all damping coefficients, about 70~90°, which is very close to the experimental results (about 50~90°) and simulation results (about 70~90°) of Boudreau [25,26]. The large pitch amplitude will inevitably make the foil work at a large angle of attack. The wake vortex diagrams of the flapping hydrofoil with four damping coefficients under the baseline case are observed and shown in Figure 8. It shows that the flow field of the double-elastic-constrained flapping hydrofoil is mostly unsteady separation. As shown, except for the case without damping (Figure 8a), there are obvious vortex street phenomena in the other three damping coefficient cases, and the shape of the separated vortex is also very similar, showing a standard 2S wake pattern (comprising two opposite-sign single vortices shed per cycle). This is consistent with the high-efficiency 2S wake discussed by Zhuo Wang et al. [40], and the efficiency of all the three working conditions is more than 30%. The difference between the latter three cases is just the width of the wake. With the increase of the damping coefficient, the width of the wake gradually decreases. Obviously, this is caused by the decrease of the heave amplitude, which is consistent with the variation of ΔY_max in Figure 7.

In order to further analyze the wake dispersion when ζ = 0.0, Figure 9 further present the time history curves of heave motion under two damping coefficients (ζ = 0.0, ζ = 0.177). We find that without damping (ζ = 0.0), the nonperiodic motion law of the double-elastic-constrained flapping hydrofoil is very obvious, but under a large damping coefficient (ζ = 0.177), the foil movement is close to periodicity. We speculate that this is the reason why the wake cannot show a regular vortex street when the damping coefficient is 0 in Figure 8.

5.2. Flow Field, Motion, and Force Analysis

According to the efficiency analysis under different damping coefficients in Section 5.1, the working point of ζ = 0.177 is selected to further analyze the flow field characteristics, motion, and force of the flapping hydrofoil. Figure 10 shows the time evolution of wake vortex diagrams being formed during one cycle with ζ =0.177 and C_p =0.7 in the baseline case. In accordance with the result in Figure 9, there could be a great difference between the motion of the double-elastic-constrained system and the sinusoidal motion of the conventional fully constrained flapping hydrofoil. In an effort to display the foil wake vortex diagrams of different attitudes and positions, the times selected in Figure 10 are not divided equally but according to the time when different pitching angles (θ = 0, θ = θ_max/2, θ = θ_max) are reached. It can be seen from Figure 10 that, during every half cycle, a large leading-edge vortex is produced (Figure 10c,g), developed backward (Figure 10d,h), and shed off at the lowest pitch position (Figure 10a,e). It shows an obvious well-defined 2S wake pattern, which has been discussed by Zhuo Wang et al. [40]. Furthermore, it should be noted that when the flapping hydrofoil approaches the equilibrium position, that is, when the pitch angle is close to 0°, such as in Figure 10a,e, a topical separation vortex will be generated on the leading edge of the foil. The topical separation vortex is marked by the arrow in Figure 10a,b,e,f. However, due to the pitching motion of the hydrofoil, the separated vortex turns to the upstream side of the hydrofoil. Then, the topical separation vortex could be impacted by the incoming flow, and it does not have the conditions for continuous development. This separation vortex gradually disappears in the subsequent development, which could be confirmed in Figure 10c,g. In addition, it seems to not affect the subsequent vortex’s street development. According to the time evolution of vortex development in Figure 10, it is speculated that this topical separation vortex is caused by the excessively high pitching speed of the foil when passing through the 0° position.

The movement trajectory and phase between pitch and heave motions have crucial impacts on the energy-extraction efficiency of a flapping hydrofoil. Many previous studies of fully constrained [18] and semi-active flapping hydrofoils have been conducted on the basis that high efficiency is achieved when the pitch–heave phase is near 90°. However, in the case of the double-elastic-constrained system, the pitch–heave phase is self-sustained, which is determined by the fluid–structure interaction of the system. Here, in Figure 11, we give the time traces of the double-elastic-constrained foil motions of the highest power coefficient cases (ζ = 0.177). Firstly, the results show that there is great symmetry between the downstroke and the upstroke as a whole. This has been verified by the results of optimized cases in a lot of research [27,30,40]. Another noticeable feature that can be observed is the heave motion is very close to being sinusoidal. Nevertheless, this is not the case for the pitch motion; the pitch motion curve is relatively plump and full. In other words, the pure heaving phases last for a shorter time, which leads to longer stroke reversal phases. Therefore, the pitching motion could be divided into three distinct phases during the cycle, including fast pitching regime (region I), slow pitching regime (region II), and transitional pitching regime (region III). The transitional pitching regime happens during both transitions between the upstroke and the downstroke. That is, when the flapping hydrofoil moves to the position near the zero-degree pitch angle, the pitch rate takes a fairly high value, resulting in the topical separation on the downstream side of the hydrofoil; at the position near the maximum pitch angle, the pitch rate is much lower than previously, which gives time for the full development of the leading-edge vortex. All of these are consistent with the results displayed in Figure 8 and Figure 10. Lastly, according to the data analysis in Figure 11, the phase difference between pitch and heave motion is about 70°, which is also very close to the phase difference obtained in the research of Mumtaz Qadri [27] and Zhuo Wang [40] of optimized cases. All of these phase differences achieved are lower than those optimized in fully constrained flapping motions. While the reasons for this discrepancy are not presently clear, there are two possible reasons. First, the phase associated with high efficiency is a function of the pivot location, which has been discussed in our former research [19]; second, it might not be possible to achieve such a large pitch–heave phase of 90° in a self-sustained system. Further investigations about this are necessary in the future.

According to Equation (3), the average mechanical power provided by the turbine in order to produce electricity could be defined through the F_yv_y curve. Figure 12 shows the time traces of F_y, v_y, and F_yv_y with damping coefficient ζ = 0.177 in the baseline case. As known from the v_y curve, the motion of the flapping hydrofoil in the y direction is close to sinusoidal, which is similar to that of many fully constrained foil systems, so we could make a further comparison in future work. Nevertheless, this is not the case for F_y. From the F_y curve in Figure 12, there is a sudden jump near the moment of v_y = 0, which corresponds to the stage of rapid pitch movement of the flapping hydrofoil. Combined with the wake vortex diagrams in Figure 8 and Figure 10, it can be analyzed that the rapid pitching process of the flapping hydrofoil is corresponding to the shedding of the large separation vortex on the leading edge of the foil. Due to a large shedding vortex and a small separation vortex successively forming on the foil, which have been described in the analysis of Figure 10, F_y has two peaks, one large and one small. The large shedding vortex has completely fallen off at the end of the rapid pitching stage; hence, the hydrodynamic effect is stronger, and the F_y shows a significantly larger peak. However, at this time, there is a certain phase difference between F_y and v_y, so the average mechanical power at this stage is negative on the F_yv_y curve, which has an adverse impact on energy acquisition. As such, the power F_yv_y is positive in the pure heaving phases but negative in the stroke reversal phases. However, the F_yv_y curve shows that the positive work is still large in the whole stroke. Zhuo Wang [40] has also given the time trace curve of C_fy in the research of the double-elastic flapping system but has not discussed the F_yv_y curve, and their C_fy curve characteristics are quite similar to the F_y curve in this paper.

5.3. Effects of Pitch Motions

In the second section of this paper, the main parameters affecting the energy acquisition efficiency of the double-elastic-constrained system have been simplified to the following five ones: k_h*, k_θ*, ω_h*, ω_θ*, and ζ_h. Based on the previous discussion of damping coefficient ζ_h, in this part, we further study the influence of pitch spring stiffness coefficient k_θ* and frequency ratio in the pitch direction ω_θ* on the performance of the flapping hydrofoil.

In the previous analysis from Figure 7, we found that the pitch amplitude of the double-elastic-constrained foil keeps a large value under all damping coefficients, about 70~90°. In this case, the flapping hydrofoil works at a large angle of attack, which can be verified by the obvious separation in the flow field shown in Figure 8 and Figure 11. The force state of the flapping hydrofoil in the case of a large angle of attack is obviously not ideal, which has been discussed in a lot of research. Considering that increasing the stiffness of the torsion spring may reduce the pitch angle amplitude, this part first analyzes the performance variation of this system with different damping ratios under the different pitch spring stiffness coefficient k_θ*. In order to ensure that the frequency ratio of the pitching motion ω_θ* remains unchanged, while changing the spring stiffness ratio, the moment of inertia of the flapping-hydrofoil system is also adjusted, and other settings are consistent with the baseline case in Table 3. The performances as a function of the damping coefficient ζ for different pitch spring stiffness coefficient k_θ* are prescribed in Figure 13.

On the whole, the variation of the stiffness coefficient of the pitch spring will not affect the variation trends of these parameters with the damping coefficient ζ. Hence, the variation of these four parameters with the damping coefficient, shown in Figure 13, is basically consistent with Figure 7. As known from Figure 13a,b, in the process of increasing the pitch spring stiffness coefficient from 0.02 to 0.16, the efficiency η and power coefficient c_p of the flapping-hydrofoil system are significantly improved as a whole. The power coefficient c_p reaches as high as 0.8, the efficiency η reaches 38% when the spring stiffness is 0.16, and the damping coefficient is about 0.177. However, the change of the pitch spring stiffness coefficient has little effects on the heave and pitch amplitudes. The heave amplitude only increases slightly with the increase of the spring stiffness coefficient, and the pitch amplitude has little change, which is between 75~90°. Moreover, in addition to the four torsion spring stiffness ratios listed in Figure 13, we also analyzed another two cases where the spring stiffness ratio is 0.25 and 0.5. Unfortunately, with these two settings, the double-elastic-constrained foil cannot obtain sustained motions, which is consistent with the conclusion that the system cannot be stable after the spring stiffness ratio is increased to a certain range obtained by Duarte in 2019 [35]. For the consideration of reliability and energy efficiency of the system, pitch spring stiffness coefficient k_θ* = 0.08 is selected in the later analysis in this paper.

From the dynamic equation of pitch motion, another key parameter affecting pitching motion is the frequency ratio in the pitch direction ω_θ*. With the baseline case, we keep the stiffness ratio of the pitch spring at 0.08 and increase the frequency ratio in the pitch direction ω_θ* from 0.632 to 1.789 by changing the inertia moment I_θ of the flapping hydrofoil. Figure 14 lists the performance curves of the flapping hydrofoil with four different ω_θ*, which show that the frequency ratio in the pitch direction also has a great impact on the system performance. Noticeable, at a small frequency ratio, such as 0.632, the efficiency and power coefficient are both particularly small. In addition, in this working condition, the heave amplitude is very small, and the pitch angle is more than 100°. This phenomenon can be explained; due to the moment of inertia I_θ being large, the pitch motion of the flapping hydrofoil has poor responsiveness to the flow, resulting in a large deviation from the usual working state. Hence, the energy acquisition performance of system is poor. With the other three frequency ratios (0.894, 1.265, and 1.789), the performance variation trend shown in Figure 14 is basically consistent. Firstly, from Figure 14a, when the frequency ratio ω_θ* increased from 0.894 to 1.789, the overall trend of the efficiency η is almost the same. The efficiency with these three frequency ratios reaches more than 50% when the damping ratio is larger than 0.5, which basically reaches the highest efficiency achieved by Boudreau in 2020 [26], whereas the peak point of power coefficient c_p shown in Figure 14b is very sensitive to the frequency ratio ω_θ*. On one hand, the value of the peak point varies from 0.6 to 0.8 with the increase of the frequency ratio ω_θ*, but the law of change is uncertain; on the other hand, the damping ratio of the power coefficient peak point is different. With the increase of frequency ratio ω_θ*, the damping ratio of the power coefficient peak increases gradually, which can be seen from the three dotted lines representing the peak position in Figure 14b. Then, as seen from Figure 14c, the heave amplitude ΔY_max/c of these three frequency ratios ω_θ* does not change obviously and only increases slightly; lastly, in Figure 14d, the pitch amplitude Δθ_max decreases significantly with the increase of frequency ratio ω_θ^*, which decreases from 80°~90° to 60°~70° as a whole.

In conclusion, in the range of parameters studied in this paper, compared with the frequency ratio in the pitch direction ω_θ*, pitch spring stiffness coefficient k_θ* has a more obvious impact on the performance of the system, such as the efficiency and power coefficient.

5.4. Effects of Heave Motions

At the end of this section, the influence of the heave spring stiffness coefficient k_h* and frequency ratio in the heave direction ω_h* on the performance of the flapping hydrofoil is discussed. The double-elastic-constrained flapping-hydrofoil system obtains energy in the heaving direction, that is, the damping is applied to the heaving motion. Thus, the damping coefficients ζ will be greatly affected when the spring stiffness and the mass of the heave system are adjusted and then affect the performance of the whole system.

Firstly, increase the spring stiffness and mass in the same proportion to ensure that the spring frequency ratio in the heave direction ω_h* remains unchanged. Then, the spring stiffness ratio is increased from 0.2 to 2.0, and a series of performance curves are obtained, as shown in Figure 15. Firstly, from the efficiency curve shown in Figure 15a, it can be seen that the system with large spring stiffness ratio coefficient k_h* has higher efficiency η with the same damping coefficient, and the maximum efficiency even reaches 67.7%, which exceeds the maximum value of 59.26% by Bates’s theory. This is because the definition of hydraulic efficiency η in Equation (4) used maximum heave amplitude at the pitching axis position ΔY_max as the sweep width. When the pitch angle of the flapping hydrofoil is large, the actual swept width of the leading edge or trailing edge of the foil is much larger than that of the pitch axis position. In order to compare these two kinds of efficiencies, with k_h* = 2.0, the hydraulic efficiency η₂ calculated by the larger value of the maximum heaving amplitude in the leading-edge point or trailing-edge point is also listed in Figure 15e. The curve of hydraulic efficiency η₂ with the damping coefficient shows a characteristic of a single peak, which is very close to that of the power coefficient c_p. In addition, the peak value of efficiency η₂ is about 30%. It can be further concluded that as metrics characterizing the energy extraction of the double-elastic-constrained flapping-hydrofoil system, power coefficient c_p should be more reasonable than hydrodynamic efficiency η. As shown in Figure 15b, the power coefficient c_p obviously decreases with the increase of the heave spring stiffness coefficient k_h*, which is opposite to that of efficiency η. When the heave spring stiffness ratio k_h* is 0.2, the system does not form stable oscillation under a small damping coefficient. Hence, the highest power coefficient c occurs when the heave spring stiffness ratio k_h* is 0.5, reaching 1.05, and its corresponding efficiency η is 37.8%. Although it is possible to further optimize the power coefficient c_p in the range of k_h* = 0.2 ~0.5, considering the risk that the flapping hydrofoil cannot work normally, the working condition of k_h* = 0.5 is selected in the later analysis. In Figure 15c,d, the variation trend of the heaving amplitude and pitching amplitude with the heave spring stiffness ratio k_h* is consistent, and both of them increase significantly with the decrease of the spring stiffness ratio.

Then, on the basis of k_h* = 0.5, this part further discusses the influence of frequency ratio in the heave direction ω_h* on the performance of the flapping hydrofoil. As shown in Figure 16, this paper gives the foil performance and motion characteristic curves with four different frequency ratios ω_h* (0.577–1.414). In Figure 16a, the influence of frequency ratio ω_h* on hydraulic efficiency η is not very clear. On the whole, the efficiency curves first increase rapidly with the increase of the damping coefficient and then gradually tend to increase slowly. Only under different frequency ratio ω_h* conditions, the rate of hydraulic efficiency η increasing with the damping coefficient is slightly different. As seen in Figure 16b, the effect of frequency ratio ω_h* on power coefficient c_p is obvious. The peak values of the power coefficient are different with different frequency ratio ω_h*, and the damping coefficients ζ of the peak point are also different. On the whole, the damping coefficients of the peak points increase gradually with the increase of frequency ratio ω_h*. It is noticeable that, when the frequency ratio in the heave direction ω_h* is one, the peak of the power factor is the lowest. With the increase or decrease of frequency ratio ω_h*, the peak of the power factor increases for both. The power factor reaches the peak value of 1.05 when the frequency ratio ω_h* is 0.707 and 1.414. It is conjectured that when the frequency ratio in the heave direction ω_h* = 1, it is close to pitch motion frequency ratio ω_θ* = 0.894. Hence, the system may tend toward resonance, which may result in a decrease of the system power factor. The particularity of this operating point (ω_h* = 1) can also be seen in the other three figures in Figure 16. Compared with the other three working conditions of frequency ratio ω_h*, the performance curve and motion curve of this working condition (ω_h* = 1) fluctuate more irregularly. Therefore, the influence of frequency ratio ω_h* on the power factor of the flapping hydrofoil needs more and further research. The heave amplitude and the pitch amplitude, as shown in Figure 16c,d, are also obviously affected by frequency ratio ω_θ*. On the whole, with the increase of the frequency ratio in heave motion ω_θ*, the heave amplitude and the pitch amplitude both increase, except for some working points.

Since the reduced frequency f* is an important parameter in the engineering design of a flapping-hydrofoil turbine, Matthieu Boudreau obtained the variation law of the reduced frequency as a function of the inflow velocity and several dimensionless structural parameters in his research in 2018 [25]. Here, we further enrich the influence parameters. Five structural parameters have been independently varied around a baseline case. Through the statistics of all working conditions involved in this paper, in a large working condition range, the reduced frequency f* of a double-elastic-constrained flapping-hydrofoil system keeps in the range of 0.11 to 0.15. The range of reduced frequency f* is shown by the two blue dotted lines in Figure 17. From the previous results, it can be concluded that the working condition with the highest power coefficient c_p of the double-elastic-constrained system is usually near the damping ratio of 1.0–1.5. From the range in the red circle in Figure 17, the reduced frequency of 0.15 is a more moderate choice, which can provide a suggestion for the mechanical design of the double-elastic-constrained system.

6. Conclusions

The double-elastic-constrained flapping-foil turbine, for which the foil is supported by two springs both in heave and in pitch, has been investigated in this paper through a two-dimensional numerical method at a Reynolds number of 2 × 10⁵. Through the analysis performance of the turbine blade dynamics, flow field, motion, and force of the flapping hydrofoil, the following conclusions are obtained:

(1)

In the range of parameters studied in this paper, the power coefficient c_p of the double-elastic-constrained flapping-foil turbine exceeding 1.0 has been achieved, the corresponding efficiency is up to 37.8%, and it is even comparable with the most efficient modern turbines, which capture around 35~45% of available flow energy [26]. The optimal system dynamic parameters are k_h* = 0.5, ω_h* = 0.707, k_θ* = 0.08, ω_θ* = 1.265, and ζ_h = 0.53, respectively; in addition, compared with hydrodynamic efficiency η, power coefficient c_p is more reasonable to measure the energy-extraction ability of the double-elastic-constrained flapping-hydrofoil system.

(2)

On the whole, this double-elastic-constrained flapping-foil system can start and work regularly within a wide range of damping coefficients. However, the hydraulic efficiency and power coefficient are sensitive to the change of the damping coefficient. When other parameters remain unchanged, with the increase of the damping coefficient, the hydraulic efficiency increases rapidly at first and then gradually tends toward a maximum. The highest efficiency can reach about 50%, whereas the power coefficient does not change monotonically with the increase of the damping coefficient; there is a peak. Noticeably, for different working conditions, the optimal damping coefficient changes greatly. For example, the optimal damping coefficient of the basic case in this paper is 0.17, and the optimal damping coefficient obtained after optimizing the parameters is 0.53. So, it is very necessary to optimize the damping coefficient for different working conditions.

(3)

The motion parameters of the double-elastic-constrained flapping-foil system are also affected by the variation of the damping coefficient. The heave amplitudes of the foil always decrease gradually with the increase of the damping coefficient. Although the pitching amplitudes also show a downward trend, the pitch amplitude of the foil keeps a large value under all damping coefficients (about 70~90°). The large pitch amplitude will inevitably make the foil work at a large angle of attack. In addition, the damping coefficient will also affect the periodic characteristics of the flapping hydrofoil’s motion. At a higher damping coefficient, the motion of the foil presents obvious periodic characteristics. The heave displacement is close to the sinusoidal curve, while the pitch angle curve is much fuller; at a small damping coefficient, the periodicity of the foil motion is not obvious, and the randomness is enhanced. This randomness characteristic is more obvious when there is no damping.

(4)

In the study of the four main parameters of heave spring stiffness coefficient k_h*, pitch spring stiffness coefficient k_θ*, frequency ratio in the heave direction ω_h*, and frequency ratio in the pitch direction ω_θ*, the most obvious parameter affecting the energy acquisition performance of the system is the spring stiffness coefficients. With the increase of pitch spring stiffness coefficient k_θ*, the efficiency and power coefficient of the flapping hydrofoil increase significantly; the influence trend of heave spring stiffness coefficient k_h* on the efficiency and power factor is different. With the increase of the stiffness ratio of the heave spring stiffness coefficient k_h*, the efficiency increases obviously, while the power coefficient decreases. However, frequency ratios in the heave direction and in the pitch direction (ω_h*, ω_θ*) both have little influence on the peak value of the efficiency and power coefficient, but they will cause the change of damping coefficients of the peak point, and the law is not very clear, which needs further research.

(5)

Under all working conditions involved in this paper, the reduced frequency of the double-elastic-constrained system is roughly in the range of 0.11~0.16. If the range of maximum power factor is considered, the optimal reduced frequency is about 0.15, which can provide a suggestion for the mechanism design of the double-elastic-constrained system.

In summary, high-efficiency results are obtained for many combinations of the parameters, and the investigation could help us to identify the range of the parameters providing high energy-extraction efficiency. The numerical results agree well with previous experimental measurements and, in addition, allow access to the velocity and vorticity fields, as functions of space and time, and input power characteristics and the force characteristics, which in turn allows us to identify the underlying energy-extraction mechanisms more easily.

These findings are based on results obtained by varying one parameter while keeping the other parameters fixed. Hence, for definitive conclusions, a more exhaustive investigation of the parameters should be performed in the future. On the one hand, this paper found the frequency ratios in the heave direction and in the pitch direction (ω_h*, ω_θ*) both have little influence on the peak value of the efficiency and power coefficient, but they will cause the change of the damping coefficients of the peak point, and the law is not very clear. Therefore, further analysis of various parameters can be made in this regard; on the other hand, all the conclusions of this paper are based on the analysis results of numerical calculations, which can be verified by a series of experiments in future work, enrich the test data of the flapping foil, and further demonstrate the correctness of this paper. Furthermore, the 3D effects are inevitable when performing such simulations. This task can be imagined to be very arduous. However, it is of great significance for foil-based flow structure, aerodynamic/hydrodynamic performance, and energy collection efficiency.