Investigation of phase difference and separation distance effects in the design of a dual flapping hydrofoil turbine

In this study, a parametric analysis was conducted for the design of a suggested dual flapping hydrofoil turbine. A Navier–Stokes‐based computational fluid dynamics code was utilized for the analysis while varying the pitch angle, reduced frequency, separation distance, and phase difference. The smallest pitch angle of 60° and the reduced frequency of 0.1 among the given ranges of the parameters were fixed for the first parametric analysis. After an assessment, the 90° front‐lead with the longest distance of 6c and the 90° rear‐lead with the shortest distance of 2c were chosen for further parametric analysis based on a suggested performance indicator and a margin for improvement. It was derived from a subsequent parametric analysis that a 70° pitch angle, 0.12 reduced frequency, and a 4c distance for the 90° front‐lead was the final optimum with 59.48% efficiency and 55.44% fluctuation according to the performance indicator as well as the system length and power balance of dual hydrofoils.


| INTRODUCTION
Due to the climate change crisis, the importance of renewable energy as an alternative to fossil fuel energy, which has caused the greenhouse effect, continues to increase. [1][2][3] The tidal stream turbine, which extracts tidal energy as a type of renewable energy resource has the advantage of being environmentally friendly and predictable. Accordingly, many studies have focused on this type of turbine. 4,5 Tidal stream turbines come in a variety of types, such as horizontal axis turbines, vertical axis turbines, and a type known as the flapping hydrofoil turbine (FHT), inspired by the movement of aquatic animals.
Flapping-type turbines were initially developed for wind turbines by McKinney et al. and, as mentioned earlier, generate electricity through the periodic flapping of hydrofoils, similar to the flapping of the fins and tails of aquatic creatures and the wings of birds. 6 Unlike the conventional method that uses rotating blades, this type has the advantage of being able to be installed in shallow water. For this reason, many studies have attempted to apply this system for tidal stream power generation. 7 FHTs can be divided into the pitch-heave, rightswing, and left-swing types according to the moving trajectory of the hydrofoil. Sitorus and Ko compared the performance capabilities of these three trajectories and found that the left-swing trajectory has the best output. 8 Ma et al. conducted a study using computational fluid dynamics (CFD) to confirm the effect of the trajectory parameter on the output, finding that the swing radius and heaving amplitude have important effects on the output. 9 Nguyen Le Dang et al. studied the camber effect while focusing on the stability and efficiency of the rightswing trajectory. 10 As part of the effort to improve the output efficiency, studies of dual or tandem systems, referring to those that use of two or more hydrofoils, have been conducted. Lindsey designed a system in which two hydrofoils have a right-swing trajectory, a phase difference of 90°, and a separation distance between the two hydrofoils of 4.8c, and verified the feasibility of developing this type of dual system through a water tunnel test and CFD research. 11 Platzer et al. using CFD confirmed the potential to develop tidal stream turbines that use pitch-heave motion with single and dual hydrofoils. 12 Kinsey et al. designed a prototype that used a crank-rocker linkage with hydrofoils with pitch-heave motion in a seesaw shape and verified it experimentally. 13 Afterwards, Kinsey and Dumas of the same group conducted a study on the optimal design by setting the phase difference to range from 90°to 180°and the separation distance from 3.6 to 7.5c through a CFD analysis of a similar system 14 ; by comparing the output of a single system, the result of the dual system was shown to be superior. 15 Similarly, Ma et al. studied the movable characteristics while varying the separation distance between the two hydrofoils (3-6c) and the initial pitch angle to explore the effect of the wake of the front and rear hydrofoils of a dual system with a pitch-heave trajectory. 16 Kim et al. conducted a test in a towing tank of a small-scale prototype designed with two hydrofoils in a dual system configuration with a 90°phase difference, a right-swing trajectory for the front hydrofoil, and a left-swing trajectory for the rear hydrofoil. They were able to verify the performance of the prototype using 2D CFD simultaneously. 17 Sitorus et al. studied a dual system with a left-swing trajectory for the front hydrofoil and a right-swing trajectory for the rear hydrofoil through a CFD analysis. 18 In other studies of flapping mechanisms, Muscutt et al. studied the motion of two pairs of flippers in ancient marine creatures and found that the phase difference and separation distance between the two flippers have considerable effects on the propulsion efficiency. 19 Ji et al. suggested an optimization framework based on an active learning method to determine the conditions of tandem flapping wings for optimal performance in terms of thrust or efficiency. 20 In this study, to design a turbine system consisting of two right-swing hydrofoils mimicking a four-legged marine creature, the phase difference and the separation distance of the dual system are determined through a parametric study using CFD.

| Configuration of the targeted dual FHT
A dual FHT with two right-swing hydrofoils located downstream w.r.t. the swing axis is designed, inspired by the locomotion of fore-limb underwater fossils and living reptiles, such as the Plesiosauria and the sea turtle. 21 This configuration was also applied in an earlier FHT study. 11,12 In addition, it showed good stability when the hydrofoil is located downstream according to a more recent study. 10 In this dual system, the power of the two swing hydrofoils is transmitted to the final unidirectional rotating shaft by means of the crank-rocker linkage shown in Figure 1, a typical mechanism for converting reciprocating rotation into unidirectional rotation, as applied in previous experimental studies. 13,22 In this system, a phase difference of ±90°was selected between the front and rear hydrofoils according to earlier findings showing that this offers the advantage of smoothing power fluctuations. 17 To realize the actual F I G U R E 1 Dual flapping hydrofoil turbine conceptual drawings. mechanical configuration, the phase difference was imposed by setting the difference between the front and rear crank angles, as shown in Figure 1. The separation distance between the two hydrofoils can be imposed by adjusting the distance between the two cranks. As the separation distance (Lx) increases, the length of the entire system increases; thus, 6c was selected as the longest distance.
Once the design of this system is finished, the goal was to conduct functional and power generation experiments in a circulating water tank or towing tank in Korea. Accordingly, earlier dual FHT specifications and experimental conditions 17,23 were referred to and system information was then devised, as listed in Table 1. Because the NACA0015 type is frequently used as the cross-section, whereas the NACA0020 type, which is 20% thicker, is also known to have similar hydrodynamic properties, 7 the latter given its 20% greater thickness was selected for the purpose of guaranteeing the structural safety of the small-size model devised here.
First, the design parameters to be determined in this study were the phase difference angle and the separation distance, and the selected minimum separation distance was 2c. This value was also used as the separation distance in previous studies. 12,17 Next, the ranges of kinematic parameters that affect the main performance, in this case the power efficiency and fluctuation, are determined as follows: a high-efficiency frequency range in a tandem configuration similar to that of our system, 23 0.1-0.14, was chosen; and for the pitch angle, a range from 60°to 75°was used, as it includes 70°u sed in the aforementioned study 23 (Figure 2). Typically, a flapping hydrofoil for harnessing tidal stream energy consists of two simultaneous motions: the flapping motion of the arm ψ t ( ) and the pitching motion θ t ( ). Both pitch motions of the front and rear hydrofoils are described as follows: The flapping motions are expressed as The phase difference φ between the pitching and flapping motions has a constant value of 90°.
The instantaneous positions of the pitch axis are then given, as follows: The velocities in the x-and y-directions as well as the angular velocity become

| CFD code
An in-house code, in this case a parallelized multiblock-structural Navier-Stokes solver, is adopted; it was validated and has been utilized in several FHTrelated studies. 8,18 To deal with the relative motion between the body-fitted mesh and the domain mesh for the FHT, the multigrid method known as the chimera method is used. Details of the chimera grid system for dual hydrofoils are shown in Figure 3. The square background domain (1) has a size of 40c in the x-direction and 40c in the y-direction. The mesh at the center is finer than that at other locations. Subfigure (2) shows the location of the body-fitted mesh in the grid system. This composition consists of body-fitted C-type meshes covering the boundary of the hydrofoils and an H-type domain mesh. The body-fitted meshes have 1e − 4 as the first layer with orthogonality ensured along the chordwise direction of the hydrofoils. Subfigures (3) and (4) present zoomed-in views of the leading edge and the trailing edge, respectively.

| CFD code validation
To verify the accuracy of this code, the mesh density was initially divided into three meshes (mesh 1: 297 × 50 F I G U R E 2 Trajectory of the dual right-swing configuration with a phase difference. Referring to (A) in Figure 4, the error at the peak point (t/T = 0.5) was 0.82% in mesh 1 and 0.09% in mesh 2, showing a relatively small error value based on mesh 3. In (B), at t/T = 0.75, mesh 1 had an error of 8.43% and mesh 2 showed 1.20%, confirming that the relative errors of the result of mesh 2 compared with that of mesh 3 were very small for both the front and rear hydrofoils. At this time, the result of mesh 1 is relatively poor and a phase difference between the peak points is also present. Considering the error range and the time taken for the CFD simulation, the quality level of this mesh approximated that of mesh 2.
Additionally, the results of a mesh with quality similar to that of mesh 2 were compared with those of a previous study 14 using the NACA0015 foil in which Re = 500,000, f* = 0.14, Lx = 5.4c, x p = c/3, h 0 = 1c, θ = 70°, and ϕ 1−2 = −180°. As shown in Figure 5, the C L curves of the front and rear hydrofoils are very similar to those in the previous study.

| Performance measures
In this study, a mixed value representing the power efficiency and fluctuation is used as a performance indicator. First, the definition of the efficiency is as follows. The lift force or vertical force based on C y is expressed as The drag force or the horizontal force based on C x is described as The moment at the pitch axis based on C m is calculated by F I G U R E 5 C L comparison of front and rear hydrofoils.
Hence, the total power of the three components is as follows: The time-averaged value of the total power is obtained via The efficiency is defined as summed average power levels (P ) S of front and rear hydrofoils divided by the available power: The power available from the flow passing through the frontal area swept by the foil is defined as Here, d is the maximum swept distance of the hydrofoil heave motion, as shown in Figure 6.
Next, the calculation formula of the variability from a previous study 24 was adopted here to determine the fluctuation level, and for the FHT, both the torque and rotation speed change. Accordingly, the power was used here instead of the torque as used in earlier work.
In this equation, n is the number of data instances per period. When the fluctuation level is high, it is difficult to control the power generation output; thus, a mechanical element such as a flywheel 23 is used for smoothing. The greater the variability is, the larger is the flywheel mass required. Considering the difficulty of power generation control and the burden of additional components, the two measures were regarded as the main factors affecting the power generation performance, and the following value was used as a single performance indicator. In this case, because the variability is considered to have a negative effect, the term for the variability was multiplied by a negative coefficient (−f 2 ).

| CFD simulation plan for a parametric study
To determine the design parameters, a simulation procedure was devised, as follows. First, the pitch angle and reduced frequency are set to the corresponding minimum values within the available ranges established in Section 2.1, and the effects of the two design variables, in particular the separation distance and phase difference, were assessed. At this time, as a fixed parameter, the pitch angle is set to 60°and f* is fixed at 0.1. The values of these parameters serve as the low bound of the corresponding ranges in this study, and they are chosen because efficiency rates that exceed 40% are obtained with these values, in a previous parametric study. 14 As a variable parameter, the separation distances used here are 2c, 3c, 4c, 5c, and 6c and the phase differences are 90 and −90 to interpret 10 cases. Next, we fix the reduced frequency to the minimum value of 0.1 and increase the pitch angle, which is a key kinematic factor, to 65, 70, and 75 to check the effects on the two design variables. Also, we fix the pitch angle to the minimum value of 60°and increase the reduced frequency, which is another key kinematic factor, to 0.12, 0.14, and 0.16 to check the effects on the two design variables. In this case, the separation distances of the minimum 2c and maximum 6c are used. If the basic condition set is extracted in this way, the first selected sets can be organized as shown in Table 2.
The optimal solution is primarily obtained through a parametric analysis of the 30 condition sets above. Subsequently, we attempt to derive the final optimal solution by exploring the effects of additional adjustments based on the first optimal solution. The global phase difference Φ 1-2 was used as an additional measure F I G U R E 6 Swept area of the flapping foil.
to find the optimal point; it was reported from earlier work 14,23 that high efficiency could be obtained at the global phase difference between 15% and 70% of 360°, and the peak efficiency was observed in the vicinity of 25% of the global phase difference, with these ranges changing slightly depending on the Reynolds number and reduced frequency. The global phase difference Φ 1-2 is defined as T A B L E 2 Condition sets for the first parametric analysis and results. 3 | RESULTS AND DISCUSSION 3.1 | Effects of the separation distance and phase difference in the basic condition sets First, when f* and the pitch angle are set to the minimum values in the usable range, to check the effects of the phase difference and the separation distance between the two hydrofoils on the output, increases of 1c from 2c to 6c and in the front-lead (90°d ifference) and rear-lead (−90°difference) conditions were compared.
First, when we check the efficiency in Figure 7, the rear-lead showed a similar level of efficiency between 39.38% and 40.89% over all separation distances, whereas for the front-lead, as the separation distance increased, the efficiency increased from 39.16% at 2c to 47.82% at 6c. The average change rate of the efficien- is calculated by 1.38. For the frontlead, the global phase difference over 360°increases by 95%, 5%, 15%, 25%, and 35%, as listed in Table 2, and it is speculated that the efficiency is increased by moving from the low-efficiency band to the high-efficiency band, 20%-80%. On the other hand, for the rear-lead, the global phase differences are 45%, 55%, 65%, 85%, and 95%, and the overall efficiency value is low and tends to decrease due to the effect of moving from the high-efficiency band to the low-efficiency band. In addition, it was also found that an efficiency rate over 47% could be obtained at 25% and 35% in the highefficiency band.
The fluctuation, as shown in Figure 7, was between 31.34% and 61.92% in the rear-lead. On the other hand, for the front-lead, the minimum fluctuation value was 38.15%, whereas 58.63% was obtained as the maximum value. The average change rate of the fluctuation,  (13) was therefore calculated based on the fact that the ratio of f 1 to f 2 is set by 9-1. The results for all condition sets are listed in Table 2. When both the efficiency and fluctuation were considered, it was found that the front-lead in 5c (case 08) is a local optimum with the highest performance indicator, 38.66 due to the highest efficiency at 25% Φ 1-2 /360 as similarly as in previous studies. 14,23 Meanwhile, for further parametric analysis, the cases that have a margin for improvement by varying the global phase difference were required; thus, we chose the rear-lead (performance indicator value, 33.67) in the proximity condition of 2c and the front-lead (performance indicator value, 37.35) in the longest condition of 6c. Accordingly, when selecting the next condition sets, the intermediate distances were omitted; that is, only 2c and 6c were considered.

| Pitch effect
In the case of a single FHT, the efficiency increases and then decreases as the pitch angle increases, and the peak angle tends to change according to the reduced frequency. 8 To check the effects of the design variables when the pitch angle was increased for this DFHT, the efficiency and fluctuation values were checked by increasing the pitch angle by 5°from 60°to 75°with separation distances of 2c and 6c.
When comparing the efficiency rates, as depicted in Figure 8A(a), the efficiency of the rear-lead efficiency is shown to be higher than that of the front-lead in the 2c condition regardless of the change in the pitch angle. As shown in Figure 8A(b), at a 60°pitch angle, the efficiency of the front-lead was 7.57% better, but when the pitch angle was increased to 75°, the rear-lead was higher by 8.07% than the front-lead. It is showed from the rear-lead and the front-lead curves in Figure 8B that the fluctuation increases as the pitch angle increase in the 2c case. In particular, in the 6c case, the levels of fluctuation at a high pitch angle were quite high at 60.62% for the front-lead and 77.31% for the rear-lead.
When checking the performance indicator according to Table 2, increasing the pitch, which does not change the global phase difference, had a negative effect on the front-lead cases due to the decrease in the efficiency and increase in the fluctuation except in the case of 6c front-lead of which the pitch angle rises from 60°to 65°while increasing the pitch had little effect on the rear-lead cases when applying a low bound of 0.1 of the f* range.

| Frequency effect
When the reduced frequency, which affects the global phase difference, was varied, a quite high efficiency, 50.46% and the lowest fluctuation, 25.53% are acquired at 0.12f* of 6c front-lead and 60°pitch angle (case 30) as shown in Figure 9; thus, according to the performance indicators in Table 2, the case 30 is superior to all other cases and it also has a margin for improvement due to that its global phase difference, 47% is a little far from 25%.

| Additional conditions and final optimal solution
Among the 34 condition sets, the case 30 was extracted as the optimal condition in the first parametric analysis due to the highest performance indicator value. The power curves of the optimal case and, those of case 23 as its counterpart in terms of the phase difference and the separation distance are provided below and were subject to a further investigation.
As shown in Figure 10, the P x and P y curves for the rear in cases 23 and 30 show similar trends when considering the phase of the corresponding velocity with a slightly smaller peak-to-valley amplitude as compared with those in the front case. Meanwhile, the P z curves in the front for case 23 and in the rear for case 30 have larger amplitudes of the negative peaks than those of the positive peaks, while the P z curves of the rear for case 23 and the front for case 30 have larger amplitudes of the positive peaks than those of the negative peaks. The negative power at the front for case 23 enlarges the peakto-valley amplitude for the front power, which causes a large fluctuation value of 32.37%. Meanwhile, the negative power of the rear for case 30 enlarges the peak-to-valley amplitude of the rear power close to that of the front power, which causes a small fluctuation of 25.53%.
For an in-depth analysis, the pressure and vortex contours close to the front foil for case 23 and rear foil for case 30 are presented in Figure 11. On the lower surface of the rear foil of case 30, a leading edge vortex (LEV) shown in red is generated at t/T = 0.4 and then moves from the middle position at t/T = 0.5 to the trailing edge at t/T = 0.6; subsequently, the LEV is shed at the trailing edge between t/T = 0.7 and 0.8, and an LEV on the upper surface (shown in blue) is then generated at t/T = 0.9, as shown in the vortex contours on the right side of Figure 11. Among the pressure contour plots on the left side of Figure 11, the highest pressure difference on the lower and upper surfaces, which results in high force, is observed at t/T = 0.8, and the peak in P y curve is located in the vicinity of t/T = 0.75, at which the highest V y occurs in Figure 10. Meanwhile, the peak of P x is located at t/T = 0.67 due to the high-pressure difference and the highest V x at t/T = 0.65; the negative peak of P z is appeared at t/T = 0.35 due to the creation of the red LEV at the leading edge before t/T = 0.4, which impedes the pitch rotating movement, and the highest θ̇at t/T = 0.25. For the front foil in case 23, similar trends of the force curves and pressure distributions with higher strength levels are observed as compared with those of the rear foil in case 30. In terms of the characteristics of the vortex shedding, three wake types, in this case an undulation wake, LEV, and reverse von Karman vortex street, were reported based on the reduced frequency and angle of attack (AOA) from experiments with the particle image velocimetry technique. 25 The characteristics of the front foil in case 23 that match the undulation wake differ from those of the rear foil in case 30 that matches the LEV, which is reported to exist at a higher AOA. 25 Moreover, it has been reported that the lift force decreases as the AOA increases from above an AOA of 15°from a water channel experiment of which the condition is close to the range of Reynolds numbers used in our simulation. 26 The maximum AOAs for both cases are above 15°as well; thus, the smaller force at the rear foil in case 30 is speculated to be due to vortex shedding from the front foil, which causes the AOA of the rear foil in case 30 to increase. The trends of the forces, pressure distributions, and vortex characteristics of the front foil in case 30 and of the rear foil in case 23 can be explained in the manner similar to that of the front foil for case 23 and the rear foil for case 30. Accordingly, corresponding plots and detailed descriptions are omitted.
With regard to this primary optimal sets, to investigate the effect of the increased pitch angle, the condition sets at pitch angles of 65°, 70°, and 75°were selected as additional conditions. The condition sets at the medium separation distance of 4c, of which the global phase difference becomes close to the high-efficiency zone, were also added. These additional condition sets are listed in Table 3; thus, the effect of the global phase difference can be explored as well.
According to the results of the 6c cases in Figure 12, the efficiency level decreases but the fluctuation level increases sharply from the pitch angle of 65°; thus, the performance indicator becomes worse compared with that with a 65°pitch angle. For the separation distance of 4c, as depicted in Figure 12, the global phase difference is close to 25% of Φ 1-2 /360, and the efficiency level increases from 60°to the peak angle of 70°while the fluctuation level decreases as the pitch angle increases. Therefore, the 70°pitch angle in the 4c case is optimal due to the highest performance indicators.
The first important finding through these parametric studies is that the efficiency of the dual system is affected by the global phase difference, as reported in previous studies, 14,23 as is by the pitch angle when f* is fixed. In addition, the effect of increased pitch angle is closely related to other conditions, specifically the phase difference and the separation distance. Next, among all F I G U R E 10 Continued cases assessed here, the 70°pitch angle of (4c, front-lead), case 40 is optimal. In terms of efficiency, case 40 obtained a result of 0.595, showing 97% of 0.614 in the optimal case (70°pitch angle, 6.3c Lx, 180°phase difference, 92.2°g lobal phase difference, 0.12f*) of a previous CFD study, 14 showing 132% of 0.450 in the optimal case (70°p itch angle, 3c Lx, 130°global phase difference, 0.115f*) of another previous CFD study, 23 and showing approximately 93% of the Newman limit of 0.640 for tandem turbines. 27 Moreover, the power balance for case 40, in which front power to rear power is 29.90-29.58, is better than that of case 30, in which front power to rear power is 29.11-21.35, as listed in Tables 2 and 3.
Although case 40 may not be the global optimum, it can be one of the best local optima due to its moderate power fluctuation level, high power efficiency level, medium separation distance, and good power balance level, as found in these parametric studies involving CFD simulations of only 41 cases. The entire procedure and extracted optima in this study are presented in Figure 13.
If numerous cases with small increments in the reduced frequency, separation distance, and pitch angle are used, or if more sophisticated optimization methods are adopted, a better outcome than that of case 40 may be extracted. Meanwhile, it is shown that a simple parametric procedure based on the optimal range of the pitch angle, reduced frequency, and the global phase difference with a small number of cases could provide a good local optimum based on the suggested performance indicator, of which the efficiency is as high as 93% as compared with the Newman limit.

| CONCLUSION
In this study, to design a dual FHT system that mimics the leg structure and swimming scheme of ancient marine dinosaurs or turtles, a numerical parametric study was carried out to determine the separation distance and phase difference between the front and rear hydrofoils.  As the analysis tool, a CFD in-house code based on the Navier-Stokes equation was used, and 10 cases were initially chosen to investigate the effects of variations in the separation distance and phase difference. For evaluating the performance of each case, a mixed function of the power efficiency and fluctuation was suggested as an indicator. As a result of the parameter analysis of the 10 cases, the coefficients in the function were determined, and the rear-lead in the 2c Lx condition and the front-lead in the 6c Lx condition were chosen for investigating the effects of the pitch angle and the reduced frequency; subsequently, the optimal case of the first parametric analysis from 34 cases was obtained according to the value of the performance indicator. For the second parametric analysis, seven additional cases were selected based on the first optimum to examine the effect of increased pitch angle and the medium separation distance.
As a result of the sequential parameter analyses, the efficiency is found to be affected by the global phase difference as well as the pitch angle when the reduced frequency is fixed, while the effect of the increased pitch angle is closely related to other conditions, such as the F I G U R E 13 Flowchart of the parametric procedure. phase difference and separation distance. Consequently, among all 41 cases, the 70°pitch angle of (4c Lx, frontlead) was acquired as a final optimum due to the highest performance indicators as well as acceptable system length and power balance. It was also found that to derive an optimal design, it is necessary to analyze the effect of parameters, as in this study, in advance within the optimal ranges of design parameters, such as the reduced frequency, global phase difference, pitch angle, and others, as in previous studies.
In the future, we plan to conduct a study to verify the results of this parametric study through experiments and to derive a final design of the suggested dual FHT.