Numerical Study on the Energy Extraction Performance of Coupled Tandem Flapping Hydrofoils

Tidal current energy is a promising renewable energy source for future electricity supply. The flapping hydrofoil is regarded as a useful tool to extract the tidal current energy in shallow water. A concept of coupled tandem flapping hydrofoils under semi-activated mode was proposed in the present study. A two-dimensional numerical model, based on the computational fluid dynamics software ANSYS-Fluent, was established to investigate the power extraction performance of the coupled tandem flapping hydrofoils. The effects of the reduced frequency, pitching amplitude, moment of inertia, damping coefficient, and longitudinal distance between hydrofoils were studied. The vortices, pressure distribution, and kinetic characteristics of hydrofoils under various conditions were analyzed to reveal the interaction between the shedding vortex and hydrofoils. The energy extraction mechanism and hydrodynamic performance were analyzed. The positive interactions for energy harvesting were identified for improvements of the further performance. The peak values of efficiency and power coefficient were achieved at 0.69 and 2.13, respectively.


Introduction
Tidal current energy converter can be classified into three major types, horizontal-axis turbine, vertical-axis turbine and flapping hydrofoil (Xiao and Zhu, 2014;Zhang et al., 2019). Originally, the flapping hydrofoils were used to provide thrust for underwater vehicles (Su et al., 2012;Tian et al., 2012). Flapping hydrofoil is believed to have some advantages over rotational turbines as the tidal current energy converter. Low tip speed ratio contributes to reducing the negative impact on the environment; Flapping hydrofoil can operate in shallow water due to the rectangular swept window (Young et al., 2014;Rostami and Armandei, 2017). Extensive studies on the flapping hydrofoil for energy harvesting continued.
Considering the high power extraction efficiency requirement and the special structure of flapping hydrofoil, the tandem flapping hydrofoils are proposed. The tandem flapping hydrofoils adopt series arrangement along the flow direction, therefore, the same flow window can be shared by the two hydrofoils, which allows the hydrofoils to reach higher efficiencies. The concept of tandem flapping hydrofoils was tested by the prototype experiment (Kinsey et al., 2011). The primary pitching motion can be transmitted to the electric generator through a series of link mechanisms. The experimental results showed that the hydrodynamic power extraction efficiency can reach up to 40%. A series of numerical studies were carried out by Kinsey and Dumas (2012) for optimizing the configuration of tandem flapping hydrofoils. The positive interactions between the downstream hydrofoil and wake vortices generated by the upstream hydrofoil contribute to a higher efficiency which is up to 64%. It is reported that the leading-edge vortex (LEV) plays an important role in the enhancement of hydrodynamics force during flapping motions, and the energy extraction performance of hydrofoil is closely related to the evolution process of the LEV (Shyy and Liu, 2007;Siala et al., 2017).
The interactions were reported to be determined by the relative position and motion phase shift of the two hydrofoils, and changes in the position and phase shift can be used to change the force production and efficiency (Broering and Lian, 2012;Kim et al., 2015). The effects of the global phase shift determined by the longitudinal distance and the phase shift of the motions of the two hydrofoils on energy extraction performance were studied (Xu and Xu, 2017;Xu et al., 2019), and the highest efficiency is found when St = 0.24 and the global phase shift is about 0.3. It is reported that the highest efficiency is obtained when the distance between the two hydrofoils is equal to 5.4 times the chord length (Xu et al., 2016). Furthermore, the energy extraction performance of tandem flapping hydrofoils would decrease with lower water depth due to the interaction between the seabed and hydrofoils (Liu, 2015(Liu, , 2017Pourmahdavi et al., 2018).
The power extraction performance of multiple hydrofoils in tandem formation was studied by Karbasian et al. (2015). It is found that the overall power efficiency would increase with the increasing number of hydrofoils, and the overall power efficiency would be independent of hydrofoil numbers after a certain number. The peak of the overall power efficiency reported in the paper is up to 95% within 15 stages of hydrofoils. 3D numerical simulations were conducted to investigate the vortex structure and interaction between the upstream and downstream hydrofoils (Broering and Lian, 2015;Pourmostafa and Ghadimi, 2018). The 3D simulation results showed that a weaker leading-edge vortex was generated due to the 3D effect, and the interaction between upstream and downstream hydrofoils was weaker than that predicted by the 2D simulation. Besides hydrofoils in tandem, another type of multi-hydrofoils in parallel was studied (Isogai and Abiru, 2012;Ma et al., 2018aMa et al., , 2018b. The parallel hydrofoils were arranged perpendicularly to the flow direction which is different from the tandem hydrofoils. Most previous studies of tandem flapping hydrofoils concentrated on the power extraction performance optimization and the interaction between the upstream and downstream hydrofoils via experiments and numerical simulations. The tandem hydrofoils of the prototype experiment (Kinsey and Dumas, 2012) were actuated by the flow without external drive train, while the motions of hydrofoils especially the amplitudes of primary pitching and pitching motions are limited by the link mechanism. Meanwhile the upstream and downstream hydrofoils are connected by one rod, which means the motions of the two hydrofoils would be interacted with each other. Most numerical studies applied the full-activated controlling strategy to prescribe the typical sinusoidal primary pitching and secondary pitching motions. Thus, the motion interactions of the two tandem hydrofoils are ignored, which is not realistic for the practical application. Therefore, the semi-activated con-trolling strategy is adopted in the present study. It is reported that the moment of inertia and damping coefficient are critical parameters for semi-activated hydrofoils (Deng et al., 2015;Zhu and Peng, 2009). The effects of multiple parameters on the performance of coupled tandem hydrofoils would be studied.
In this paper, the numerical simulation based on a semiactivated controlling strategy was performed to investigate the coupled tandem flapping hydrofoils. The governing motion equations of semi-activated coupled tandem flapping hydrofoils and the numerical model based on the CFD software ANSYS-Fluent were described in Section 2. In Section 3, variations in flow fields, hydrodynamic torques, forces, and power efficiencies under various reduced frequency, secondary pitching amplitude, moment of inertia, damping coefficient, and distance between hydrofoils were analyzed. The main conclusions derived from the results were drawn in Section 4.

Semi-activated coupled tandem flapping hydrofoils
The 2D tandem configuration of two coupled hydrofoils is shown in Fig. 1. The configuration consists of an upstream hydrofoil (UF) and a downstream hydrofoil (DF) interconnected by a rod. UF and DF pitch about pivot P U and P D , respectively. And the pitching motions about the pivot are defined as secondary pitching motion which is activated by external devices. The rod is restricted to pitch only around the shaft O, and the pitching motion about the shaft O is defined as primary pitching motion which is passively driven by the periodically changing forces of water flow. Considering the motions of hydrofoils consist of primary and secondary pitching motion, the device in this paper is defined as coupled-pitching type. Most hydrofoils undergo a combined heave-pitch motion about pivot which is defined as heaving-pitching type. And the heaving-pitching type was applied in the tandem flapping hydrofoils (Xu et al., 2019;Broering and Lian, 2015) and parallel flapping hydrofoils (Isogai and Abiru, 2012;Ma et al., 2018aMa et al., , 2018b. And the guide rail is adopted to support the heaving motion. The incident velocity U is fixed at 2 m/s. The profile of hydrofoils is NACA0015, and the chord length of hydrofoils is defined as c, which is fixed at 0.24 m. The span of the hydrofoil is defined as S, and it is set as 1 m in the reference value panel of Fluent. The pivot is located at c/3 away from the leading edge. The length of the rod is defined as L and the radius of primary pitching motion R = L/2. Besides, the positive directions of the moment and angular velocity are defined to be counterclockwise. In theory, the primary pitching angle α(t) should be within the range of −90° and 90° during the stable stage. And α(t) is effected by environmental dynamics, kinematic parameters, and structural parameters. In practical application, by setting appropriate para-meters (such as reduced frequency and secondary pitching amplitude), flip-over behavior should be avoided to realize limit-cycle oscillation.

Motion and performance evaluation equations
As mentioned that the hydrofoils secondary pitching motion is imposed, and the motion equations can be defined as (Xu and Xu, 2017): (1) where θ U (t) and θ D (t) are instantaneous pitching angles of UF and DF; θ 0 is the secondary pitching amplitude; f is the frequency. Furthermore, the reduced frequency can be written as k = fc/U (Xu et al., 2016). And t is time.
The primary pitching motion of the hydrofoils under the restriction of the rod is determined by Newton's second law and realized by using the user-defined function (UDF) in the computation platform. The fundamentals and the simulation procedures are as follows. First, the pressure distribution all over the hydrofoil can be derived after initializing the calculation. The resultant torques about the shaft and the pivot are calculated by integrating the pressure and shear stress. According to Newton's second law, the primary pitching angular acceleration and the angular velocity at this time step can be calculated. The angular positions of the hydrofoils are updated, and the new solid boundary to the water is reformed. At the next time step, the renewed flow field is subsequently employed in the next-round calculation. By repeating the above procedures, the numerical simulation of the flow field and the motion of hydrofoils under the semiactivated mode is realized. The flow chart of simulation procedures is illustrated in Fig. 2.
According to the simulation procedures, the governing motion equations of hydrofoil are established. The primary pitching motion is governed by Newton's second law, and the specific equation is defined as: I dω dt where I is the moment of inertia about the shaft. The dimensionless moment of inertia is defined as m fluid is the mass of fluid with the same volume as the hydrofoil. ω is the primary pitching angular velocity. T L denotes the resistive load torque from the PTO damper which is proportional to the primary pitching angular velocity.
ρ represents the water density. T U and T D are the hydrodynamic torque about the shaft integrated from the upstream and downstream hydrofoils, respectively. T U and T D can be defined as: where i and j are node index numbers of UF and DF, respectively. x i and y i are the coordinates of node i in the x and y directions. x j and y j are the coordinates of node j in the x and y directions. x 0 and y 0 are the coordinates of shaft O in the x and y directions. p i , τ i , p j and τ j are the pressure and shear stress on node i and j, respectively. A xi , A yi , A xj and A yj are the project areas in the X and Y directions of node i and j, respectively. The dimensionless hydrodynamic torque can be written as C T U = T U /(0.5ρU 2 Sc 2 ) and C T D = T D / (0.5ρU 2 Sc 2 ).
The primary pitching angular velocity at time t+Δt can be calculated as: Furthermore, the instantaneous primary pitching angular position at the same time t +Δt can be calculated as: The dimensionless cycle-averaged conversion efficiency η is defined for evaluating the hydrofoil performance of power extraction (Xu et al., 2016): where P = T L ω + Mω P , M is the pitching hydrodynamic torque act on the hydrofoil about the pitching axis, and ω P is the pitching angular velocity. Mω P is the secondary pitch- ing hydrodynamic power, and it is equal to the external power required to actuate the pitching motion in a complete period without considering the mechanical and electric losses. d is the maximum translational displacement of the leading or trailing edges in the Y direction. The dimensionless cycle-averaged power coefficient C P is defined as (Xu et al., 2016): The dimensionless instantaneous hydrodynamic power coefficients for the primary pitching motions of UF and DF are defined as: (10) It is mentioned that the contribution to energy extraction from primary pitching motion plays a dominant role compared with the contribution from secondary pitching motion (Wang et al., 2017). Thus, the power coefficients for the primary pitching motions are analyzed in the time domain for brevity. And the power from both the primary and secondary pitching motions are considered in the analysis of η and C P .

Numerical model setup and validation
The water is considered as the incompressible fluid and employed in this study. The hydrodynamics are governed by the continuity equation ∇ · U = 0 (12) and the Navier−Stokes equations where u and v are the velocity components in the x and y directions, respectively; p represents the pressure force; f x and f y are the body forces in the x and y directions, respectively; ν is the kinematic viscosity coefficient. The continuity equation and the Navier−Stokes equation were solved by using the finite volume method through the commercial CFD software ANSYS-Fluent 16.0. The computational domain and mesh structures are shown in Fig. 3. The straight line, located at 6c away from the upstream hydrofoil pivot at the left side, was set as velocity inlet boundary. The straight line, located at 9c away from the downstream hydrofoil pivot at the right side, was set as the pressure outlet boundary. The upper and lower boundaries were both 5c away from the hydrofoil and set as the sym-metry boundaries, respectively.
The sliding and dynamic mesh technologies were employed for the mesh treatment to realize the coupled pitching motions of hydrofoils. Thus, the computational domain was divided into three parts. The two circular domains with a radius of 2.5c move with the hydrofoil which is governed by the UDF and are connected to the outer domain through the non-conformal interfaces for an accurate flux exchange. The large deformation of outer domain meshes induced by the movement inner domains requires unstructured mess. The structured meshes were used for the boundary layer around the hydrofoil surface. And y+ was controlled to be below 1. The total number of meshes was approximately 3.8×10 5 . The Reynolds number is 4.8×10 5 . The standard k−ω was chosen to deal with turbulent problems. The enhanced wall treatment was applied for the near wall function. The pressure-velocity coupling employed the SIMPLE (Semi-implicit method for pressure-linked equations) scheme. The derivatives were calculated based on the Green-Gauss mode. The second-order upwind scheme was applied for the spatial discretization of momentum, turbulent kinetic energy and specific dissipation rate. The temporal term was discretized by using a second-order implicit scheme. The absolute convergence criteria were chosen as the scale of residual of 10 −6 .
As previously mentioned, the coupled tandem hydrofoils undergo the coupled pitching motion under the semiactivated mode. Studies focused on this type of system are very limited. Due to the lack of relevant experimental data, the numerical model was validated by using an experimental study of an isolated semi-activated hydrofoil in a water tunnel (Huxham et al., 2012). The experimental model device consisted a coupled-pitching hydrofoil with a crosssection profile of NACA0012. The chord length is 0.1 m, and the span length is 0.34 m. The rod is 0.3 m in length. The incident velocity is 0.5 m/s, and the reduced frequency is 0.1. The secondary pitching amplitude ranges from 0° to 60°.
Studies of grid independence and time-step sensitivity were performed. Three types of grids were chosen with coarse, medium, and fine meshes. The number of grid nodes on the boundaries of the outer domain and hydrofoil sur- face were 100 and 208 for the coarse meshes, 200 and 523 for the medium meshes, and 400 and 1029 for the fine meshes. The total numbers of three types of meshes were approximately 2.7×10 4 , 8.7×10 4 , and 12.3×10 4 , respectively. Effects of the mesh number and time step on the calculation accuracy are shown in Fig. 4. The average difference is 8.3% between the results of coarse mesh and fine mesh, and the average difference is 3.2% between the res-ults of medium mesh and fine mesh. It can be seen from Fig. 4a that the torque variations for medium mesh agreed well with that for fine mesh. Similar to the mesh number, as shown in Fig. 4b, the larger time step shows a larger difference in the torque variations than the other two time steps. Considering the computational cost and the calculation accuracy, the medium mesh and 0.0025T time step were chosen for subsequent simulations. The experimental and numerical efficiencies under various θ 0 were compared in Fig. 5. The average error for various θ 0 was 8.6%. The overestimation mainly results from the neglecting of mechanical loss, 3D hydrodynamic losses, and the effects of the rods and shaft. The validation results demonstrated that the capability and potential of the present model on numerical simulations of the coupled tandem hydrofoils and can be applied in further investigations. Given that the numerical model satisfies the requirement of accuracy, the rod and shaft were not considered to simplify the 2D numerical model.  Fig. 6c consists of two individual hydrofoils. The two hydrofoils can independently pitch around O in the upstream and downstream region, respectively. Meanwhile, the simulation procedures and governing motion equations for coupled tandem hydrofoils can be adopted for the above configurations under semi-activated mode. The reduced frequency k is 0.12. The value of θ 0 varies from 30° to 90° with an interval of 15° and the primary pitching radius is R =2.7c. The dimensionless moment of inertia of single hydrofoil is I * = 10, and the dimensionless damping coefficient of single hydrofoil is C d * = 10.6.
The instantaneous torque coefficients, angular velocity, and power coefficients of primary pitching motion for UF and DF are shown in Fig. 7. As shown in Fig. 3a, C TU and angular velocity of isolated hydrofoil and uncoupled tandem hydrofoil are almost similar. It means that the uncoupled downstream hydrofoil has a negligible influence on the performance of uncoupled upstream hydrofoil, but C TU and angular velocity of coupled tandem hydrofoils show slight difference from other configurations, which mainly comes from the effect of DF. As shown in Fig. 3b, C TD and angular velocity of three configurations show distinct differences. C TD of coupled tandem hydrofoils is higher than those of other configurations. The differences between isolated hydrofoil and uncoupled tandem hydrofoils are mainly caused by the wake generated by upstream hydrofoil, and the differences of coupled tandem hydrofoils are also caused by the effect of coupled UF. As shown in Fig. 7c, the peak values of C PU and C PD for isolated hydrofoil are 1.66 and 1.62, respectively. The peak values of C PU and C PD for uncoupled tandem hydrofoils are 1.65 and 0.93, respectively. And the peak values of C PU and C PD for coupled tandem hydrofoils are 1.56 and 1.36, respectively. The average values of C PU and C PD for isolated hydrofoil are 0.80 and 0.76, respectively. The average values of C PU and C PD for uncoupled hydrofoil are 0.79 and 0.51, respectively. The average values of C PU and C PD for uncoupled hydrofoil are 0.74 and 0.65, respectively. According to the statistic, the ratios of averaged power coefficient of secondary pitching motion to primary pitching motion for isolated upstream hydrofoil, isolated downstream hydrofoil, uncoupled upstream hydro-foil, uncoupled downstream hydrofoil, UF, and DF are 10.1%, 7.8%, 10.7%, 15.7%, 10.3%, and 14.5%. It can be concluded that the contribution from secondary pitching motion is much lower than that from primary pitching motion. The power coefficient of isolated upstream hydrofoil is slightly higher than that of isolated downstream hydrofoil. The power coefficient of uncoupled upstream hydrofoil is similar to that of isolated upstream hydrofoil, while the power coefficient of uncoupled downstream hydrofoil is the lowest for the negative interaction with upstream wake. The power coefficients of UF and DF of coupled tandem hydrofoils are moderate. The energy extraction performance of tandem hydrofoils consists of contributions of upstream and downstream hydrofoils, and therefore the energy extraction performance of tandem hydrofoils is higher than that of isolated hydrofoil. Meanwhile, the energy extraction performance of coupled tandem hydrofoils is higher than that of uncoupled tandem hydrofoils.
The average power coefficients of various configurations are compared in Fig. 8. It can be seen that the overall power coefficient of coupled tandem hydrofoils is higher than those of other configurations. It should be noted that the hydrofoil in upstream would flip over when the secondary pitching amplitude is large. So the results of isolated hydrofoil in upstream and uncoupled tandem hydrofoils when θ 0 is 75° and 90° are missing. The motions of UF and DF of coupled tandem foils are connected, so UF would not flip over. The coupled tandem hydrofoils can operate stably and efficiently. Considering the high power extraction performance and wide operation range, the coupled tandem hydrofoils are systematically studied.   QU Heng-liang, LIU Zhen China Ocean Eng., 2022, Vol. 36, No. 1, P. 38-49 43 3.2 Effects of reduced frequency and pitching amplitude The effects of reduced frequency k and secondary pitching amplitude θ 0 on hydrodynamic and power extraction performance are studied in this subsection. The reduced frequency k ranges from 0.06 to 0.15 with an interval of 0.03. The pitching amplitude θ 0 ranges from 30° to 90° with an interval of 15°. The primary pitching radius R = 0.648 m (2.7c). The dimensionless moment of inertia I* is set as 28, and the damping coefficient C d * is set as 21.2.
The vorticity contours around the hydrofoil in a stable cycle for k = 0.06 and 0.12 are illustrated in Fig. 9. The flow structures at t =3T/4 and t = T for each value of k were observed symmetrical to those at t = T/4 and t = T/2, respectively. Therefore, only the vortex structures in the first half cycle were analyzed. As can be seen, the reduced frequency has a significant effect on the timing of the generation, evolution, and shedding of vortices from the hydrofoil surface. At t = T/4, significant LEV is observed at the middle chord of the hydrofoil and the trailing edge for both k = 0.06 and 0.12. And the differences in vortex size and strength mainly cause the slight difference of C TU for k = 0.06 and 0.12. A pair of vortexes shed off from UF and travel downstream for k = 0.06. DF has a wake interaction with upstream wake when k = 0.06, and the upstream vortex just brushes past the lower surface of DF. When k = 0.12, DF just encounters the upstream vortex, and the vortex around DF is different to that when k = 0.06. The LEV remains attached to the lower surface, and TEV is more intense. The strong interaction causes a high torque coefficient. At t = T/2, an intense LEV develops and moves to the middle chord of the lower surface when k = 0.06, and a tiny LEV grows in the front of the lower surface when k = 0.12. When k = 0.06, vortexes that shed off from the surface of UF have traveled downstream from the upward side of the DF, and the interaction is relatively weak. And an intense LEV develops at the upper surface DF. The upstream vortex moves to the middle position of tandem hydrofoils, so the DF does not interact with vortex when k = 0.12. It can be concluded that, due to the interaction between vortices of UF and DF, the effect of reduced frequency on the flow structure around the tandem hydrofoil configuration would be more significant than that around the single hydrofoil. The reduced frequency would influence the flow structure around UF. And the wake would be influenced accordingly. Different wakes caused by UF would influence the hydrodynamic performance of DF. The mechanism of the overall process is complicated which needs specific analysis.
The time histories of C TU , C TD and primary pitching angular velocity for k = 0.06 and 0.12 are shown in Fig. 10. It can be seen that the amplitude of C TU for k = 0.12 is similar to that for k = 0.06, while the amplitude of C TD for k = 0.12 is larger than that for k = 0.06. There is a more significant phase difference between C TU and C TD for k = 0.12 than that for k = 0.06. Furthermore, when k = 0.12, the amplitude of C TD is larger than that of C TU , and the difference between C TU and C TD for k = 0.06 is smaller than that for k = 0.12. The angular velocity of passive primary pitch motion is determined by the hydrodynamic torques of UF and DF. The amplitude of ω(t) for k = 0.12 is larger than that for k = 0.06. It should be noted that the minimum values of C TU and C TD for various k are not obtained exactly at t = T/4 and t = 3T/4. As shown in Fig. 7, the pitching angles of UF and DF are 0°, and the projected area of hydrofoils is minimum. Nonetheless, the pressure differences due to the asymmetric vortex of upper and lower surfaces can generate vertical forces, and the vertical forces generate the torque around shaft O.
The time histories of primary pitching angles for various k are shown in Fig. 11. It can be seen that the coupled tandem hydrofoils undergo highly periodic passive primary pitching motion. The primary pitching motion is stable and symmetric about 0°. The amplitudes of α(t) for k = 0.06, 0.09 ,0.12, and 0.15 are 26.4°, 19.7°, 17.7°, and 15.5°, respectively. It can be concluded that the amplitude of α(t) is influenced by k. And the amplitude of α(t) for k = 0.09 decreases with the increasing k. The high smoothness is shown in the primary pitching motion.
The vorticity contours around the hydrofoil in a stable cycle for k = 0.12 and θ 0 = 45°are illustrated in Fig. 12. It can be seen that the timing of the vortex generation for θ 0 = 45° is similar to that for θ 0 = 75° seen in Fig. 9b. However, the vortex strength and scale for θ 0 = 45° is significantly lower than that for θ 0 = 75°. Therefore, the interaction between DF and upstream vortex is relatively weak. At t = T/4, the LEV is generated on the upper surface of UF, and the LEV reattaches in the middle of the hydrofoil. A wake vortex is behind the trailing edge. The clockwise and anticlockwise vortexes are attached to the upper and lower surfaces of DF, respectively. The DF encounters the upstream vortex. At t = T/2, the clockwise and anticlockwise attached vortexes are generated on the upper surface and the lower surface, respectively. The wake vortex shedding from the trailing edge moves to the middle position between UF and DF. The instantaneous torque coefficients C TU and C TD for different θ 0 are shown in Fig. 13. It can be seen that the secondary pitching amplitude has significant effects on the hydrodynamic performance of coupled tandem hydrofoils. C TU and C TD for θ 0 = 45° are relatively lower than those for θ 0 = 75°. And the amplitudes of C TD for both θ 0 = 45° and 75°a re higher than those of C TU . It means that intense FSI contributes to the enhancement of the hydrodynamic force. Meanwhile, the phase differences between C TU and C TD for θ 0 = 45° and 75° are similar. The power extraction efficiency and power coefficient variations with secondary pitching amplitude θ 0 at various reduced frequencies are shown in Fig. 14. It can be seen that η increases first and then decreases with the increasing θ 0 . The peak values of 0.32, 0.53, and 0.64 can be achieved at θ 0 = 75° for k = 0.06, 0.09 and 0.12, respectively. While the peak value of 0.52 for k = 0.15 is obtained at θ 0 = 60°, and η decreases dramatically to 0.23 at θ 0 = 90°. And it should be noted that the motion for k = 0.15 and θ 0 = 90° would not enter a stable stage. The trend of power coefficient variation is similar to the efficiency. As can be seen, C P increases first and then decreases as pitching amplitude increases except for the case of k = 0.09. The peak value of C P is 1.59 for k = 0.12 at θ 0 = 75°. Considering optimal power extraction performance, k = 0.12 and θ 0 = 75° would be adopted in subsequent study.

Effect of moment of inertia
The vorticity contours around the hydrofoil in a stable cycle for I * = 20 and θ 0 = 75°are illustrated in Fig. 15. Compared with Fig. 9b, it can be seen that the moment of inertia I * has a slight influence on the flow structure. The distribution, scale, and intensity of vortex under different I * are similar.
The effects of the moment of inertia on power extraction performance are studied. The reduced frequency k is set as 0.12. The pitching amplitude θ 0 ranges from 30° to 90°w ith an interval of 15°. The primary pitching radius R = 0.648 m (2.7c). The dimensionless moment of inertia I * ranges from 16 to 28 with an interval of 4, and the damping coefficient C d * is set as 21.2. The power extraction efficiency and power coefficient variations with secondary pitching amplitude θ 0 at various moments of inertia are shown in Fig. 16. The maximum difference of efficiency under various conditions is below 12.1%, and the maximum difference of power coefficient under various condi-  QU Heng-liang, LIU Zhen China Ocean Eng., 2022, Vol. 36, No. 1, P. 38-49 tions is below 12.9%. It can be seen that the moment of inertia has a slight influence on power extraction performance. However, the efficiency for I * = 20 is slightly higher than that of other conditions. And I * = 20 would be used in the subsequent study.

Effect of damping coefficient
The effects of the damping coefficient on power extraction performance are studied. The reduced frequency k is set as 0.12. The pitching amplitude θ 0 ranges from 30° to 90°w ith an interval of 15°. The primary pitching radius R = 0.648 m (2.7c). The dimensionless moment of inertia I * is 20, and the damping coefficient C d * ranges from 10.6 to 42.4.
The vorticity contours around the hydrofoil in a stable cycle for C d * = 21.2 and C d * = 42.4 are illustrated in Fig. 17.
It is shown that the vortex for C d * = 42.4 is relatively more intense than that for C d * = 21.2. At t = T/4, the clock-wise LEV is formed on the upper surface of UF for both C d * . The obvious attached anticlockwise vortex is observed near the trailing edge on the upper surface of UF. And the clockwise vortex shed from the trailing edge of UF. The scale and intensity of shedding vortex increase as C d * increases. The DF encounters upstream vortex for both C d * . And the upstream vortex for C d * = 42.4 is more intense than that for C d * = 21.2. The anticlockwise LEV is generated on the lower surface of DF. The clockwise attached vortex is observed on both sides of DF. The anticlockwise vortex shed from trailing edge. The scale and intensity of shedding vortex for C d * = 42.4 are larger than those for C d * = 21.2. At t = T/2, LEV is generated near the leading edge on the lower surface, and the attached vortex is generated in the rear section of the lower surface. The scale and intensity of LEV for C d * = 42.4 are larger than those for C d * = 21.2. The clockwise vortex shed from the trailing edge and moves to the downstream, and the scale and intensity of shedding vortex     Fig. 18. It can be seen that the amplitude of C TU for C d * = 21.2 is smaller than that for C d * = 42.4. The amplitudes of C TD for different C d * are similar. The cycle averaged C TU for C d * = 21.2 is 3% lower than that for C d * = 42.4. However, the cycle averaged C TD for C d * = 21.2 is 18% larger than that for C d * = 42.4. It means that the DF contributes to more power than the UF does when C d * = 21.2. The overall power output for C d * = 21.2 is higher than that for C d * = 42.4. The power extraction efficiency and power coefficient variations with secondary pitching amplitude θ 0 at various damping coefficients are shown in Fig. 19. The value of η for C d * = 10.6 is lower than that for other C d * . Under a given value of C d * , the highest values of η and C P are achieved when θ 0 = 75°. The values of η first increase and then decrease with the increase of C d * , while the values of C P decrease monotonously with the increase of C d * . The peak value of η is achieved at 0.66 for C d * = 21.2, and the peak value of C P is achieved at 1.95 for C d * = 10.6. C d * = 10.6 would be used in the subsequent study.

Effect of distance between UF and DF
The effects of the distance between UF and DF on power extraction performance are studied. The reduced frequency k is set as 0.12. The pitching amplitude θ 0 ranges from 30° to 90° with an interval of 15°. The dimensionless moment of inertia I* is 20, and the damping coefficient C d * is set as 21.2. The distance L between UF and DF ranges from 3.6c to 9.0c with an interval of 1.8c.
The vorticity contours around the hydrofoil in a stable cycle for L = 7.2c are illustrated in Fig. 20. Compared with Fig. 12a, for each instant, the flow structure of UF for L = 7.2c is similar to that for L = 5.4c. However, L plays an important role in the flow structure of DF. It can be seen that the DF encounters the upstream vortex at t = T/4 for L = 5.4c. The upstream vortex is just in the front of DF without interaction for L = 7.2c. The upstream shedding vortex has not moved to the DF due to the increasing L. At t = T/2, the scale and intensity of LEV of UF for L = 7.2c are smaller  QU Heng-liang, LIU Zhen China Ocean Eng., 2022, Vol. 36, No. 1, P. 38-49 47 than those for L = 5.4c. The upstream shedding vortex moves to the middle position between UF and DF. The power extraction efficiency and power coefficient variations with secondary pitching amplitude θ 0 at various L are shown in Fig. 21. It can be seen that under a given θ 0 , η first increases and then decreases with the increase of L. The peak values of η and C P are achieved at 0.69 and 2.13 for L =7.2c.

Conclusions and future works
A 2D numerical model was set up based on ANSYS-Fluent to investigate the power extraction performance of coupled tandem flapping hydrofoils. The effects of reduced frequency, pitching amplitude, moment of inertia, damping coefficient, and distance on the power coefficient, efficiency, and flow field were studied. The reduced frequency and pitching amplitude have a significant effect on the timing of the generation, evolution, and shedding of vortices from the hydrofoil surface. The interaction between the vortices and the downstream hydrofoil plays an important role in the performance of power extraction. The peak values of η and C P are achieved when θ 0 = 75° and k = 0.12. The moment of inertia has a slight influence on the flow structure and energy extraction performance. The optimal damping coefficient and distance are 21.2 and 7.2c, respectively. The peak value of η and C P are 0.69 and 2.13, respectively.
In future studies, more parameters will be employed to reveal the mechanism of energy extraction of tandem coupled hydrofoils, and the power extraction performance will be improved. A fully-passive pattern will be induced in the 2D numerical simulation. A 3D numerical model will be established to understand the detailed flow structures and pressure distribution around the coupled tandem hydrofoils.