Numerical Investigation of the Semi-Active Flapping Foil of the Wave Glider

A numerical investigation is conducted to study the propulsive performance of the semi-active flapping foil of the wave glider, where the heaving smotion is fully prescribed, and the pitching motion is determined by the hydrodynamic force and torsion spring. A mesh for two-dimensional NACA0012 foil with the Reynolds number Re = 42000 is produced, and a dynamic mesh and sliding interface are used in the computation. The influences of reduced frequency, spring stiffness, and critical pitching amplitude on the hydrodynamic characteristics of semi-active flapping foil are systematically investigated. We find that there is a critical reduced frequency: When the reduced frequency is lower than the critical value, the propulsive performance of flapping foil can be improved exponentially, and when the reduced frequency is higher than the critical value, the semi-active flapping foil cannot provide an effective thrust. For a greater reduced frequency, there is an optimal spring stiffness value, which corresponds to the maximum value of the output power coefficient. For a lower reduced frequency, the mean value of the output power coefficient monotonically decreases as the spring stiffness increases. We also notice that the propulsive efficiency of flapping foil monotonically decreases as the spring stiffness increases. Finally, we find that the appropriate critical pitching amplitude can improve the propulsive performance of semi-active flapping foil, especially for greater heaving amplitudes.


Introduction
The wave glider is a new class of wave-propelled, persistent ocean vehicle [1,2]. It consists of a surface boat and a submerged glider, which are connected by an umbilical cable. Due to the unique architecture, the wave glider can capture the wave energy from the ocean and convert it into forward thrust. The principle of wave propulsion is shown in Figure 1 [3][4][5]. The surface boat moves with the wave, and the submerged glider, with flapping hydrofoil, follows the movement. The hydrofoil flaps under hydrodynamic forces and continues to generate thrust [5,6]. Therefore, it is essential to study the hydrodynamic characteristics of the wave propulsion system of the wave glider, especially the propulsion performance of the flapping foil, in order to optimize the wave glider. The factors affecting the wave propulsion performance of the wave glider include the coupling effect between the surface boat and the submerged glider and the propulsion performance of the flapping hydrofoil [6][7][8][9][10]. Several experiments were performed to examine the ability of the wave glider, such as its velocity and longest voyage [11,12]. Many researchers use theoretical analysis to study the performance of the wave glider, establish a six-degrees-of-freedom dynamic equation, identify the key hydrodynamic parameters, predict the wave glider velocity using linear regression analysis, and analyze the motion characteristics of the wave glider, based on the D-H (Denavit and Hartenberg) approach [7,8,13,14]. Numerical simulations were also applied to investigate the hydrodynamic characteristics of the wave glider, analyze the lift coefficient, and drag coefficients of flat and NACA foils using steady-state simulation, and study the propulsive performance of fullactive flapping foil as well as the coupling of the foil and the hull using the unsteady Unsteady-Reynolds-Averaged-Navier-Stokes (URANS) method [5,10,15]. A FINE/Marine and STAR-CCM+ cosimulation has also been used to analyze the coupling effect of the wave glider and to point out that the passive flapping hydrofoil plays a decisive role in the wave glider's propulsion performance [6]. It should be noted that there is still a lack of numerical simulations of flapping foil, especially in the determination of the propulsive performance of the semi-active flapping foil of the wave glider.
Flapping hydrofoil has complex unsteady fluid-structure coupling effects, and research on its hydrodynamic characteristics is still lacking [16,17]. Flapping hydrofoil has a better maneuverability, environmental adaptability, and efficiency than the conventional propeller propulsion mode. In particular, the semi-active flapping propulsion mode has a more straightforward mechanical structure than the fully active flapping foil [5,[18][19][20][21]. Many authors carried out theoretical, simulation, and experimental studies on the hydrodynamic performance of active and semi-active flapping hydrofoil [19,20,[22][23][24]. The form of motion for fully active flapping hydrofoil includes the heaving motion, the pitching motion, and the coupled heaving-pitching motion, which is fully prescribed [24,25]. The factors affecting the propulsion performance of fully active flapping hydrofoil mainly include a reduced frequency, the pitching profile, the Strouhal number, the pivot location, the phase between heaving and pitching motions, and the foil shape [5,23,[26][27][28][29].
The disadvantage of fully active flapping foil is that the heaving and pitching motions require complex mechanical structures, which reduce the usability and reliability of the system [21]. This shortcoming can be alleviated by using semi-active flapping foil, which can reduce one degree of freedom to simplify the mechanical structure [30]. Semi-active flapping hydrofoil can only choose the heaving motion as the active motion and the pitching motion as the passive motion in order to obtain the thrust [20,[31][32][33]. As for the difference between the principle of semi-active and fully active flapping foils, the factors influencing the propulsive performance of semi-active flapping hydrofoil are not the same as those of the fully active flapping one. Additional influencing factors include the spring stiffness, added mass, and critical pitching amplitude of the pitching motion [20,32,34]. Conventional semi-active flapping foil has a rigid body, and some scholars also carried out research work on the hydrodynamic performance of flexible flapping hydrofoil [35,36]. Limiting the magnitude of the pitching motion may reduce the propulsion efficiency of flapping foil [26,33]. The factors affecting the wave propulsion performance of the wave glider include the coupling effect between the surface boat and the submerged glider and the propulsion performance of the flapping hydrofoil [6][7][8][9][10]. Several experiments were performed to examine the ability of the wave glider, such as its velocity and longest voyage [11,12]. Many researchers use theoretical analysis to study the performance of the wave glider, establish a six-degrees-of-freedom dynamic equation, identify the key hydrodynamic parameters, predict the wave glider velocity using linear regression analysis, and analyze the motion characteristics of the wave glider, based on the D-H (Denavit and Hartenberg) approach [7,8,13,14]. Numerical simulations were also applied to investigate the hydrodynamic characteristics of the wave glider, analyze the lift coefficient, and drag coefficients of flat and NACA foils using steady-state simulation, and study the propulsive performance of full-active flapping foil as well as the coupling of the foil and the hull using the unsteady Unsteady-Reynolds-Averaged-Navier-Stokes (URANS) method [5,10,15]. A FINE/Marine and STAR-CCM+ co-simulation has also been used to analyze the coupling effect of the wave glider and to point out that the passive flapping hydrofoil plays a decisive role in the wave glider's propulsion performance [6]. It should be noted that there is still a lack of numerical simulations of flapping foil, especially in the determination of the propulsive performance of the semi-active flapping foil of the wave glider.
Flapping hydrofoil has complex unsteady fluid-structure coupling effects, and research on its hydrodynamic characteristics is still lacking [16,17]. Flapping hydrofoil has a better maneuverability, environmental adaptability, and efficiency than the conventional propeller propulsion mode. In particular, the semi-active flapping propulsion mode has a more straightforward mechanical structure than the fully active flapping foil [5,[18][19][20][21]. Many authors carried out theoretical, simulation, and experimental studies on the hydrodynamic performance of active and semi-active flapping hydrofoil [19,20,[22][23][24]. The form of motion for fully active flapping hydrofoil includes the heaving motion, the pitching motion, and the coupled heaving-pitching motion, which is fully prescribed [24,25]. The factors affecting the propulsion performance of fully active flapping hydrofoil mainly include a reduced frequency, the pitching profile, the Strouhal number, the pivot location, the phase between heaving and pitching motions, and the foil shape [5,23,[26][27][28][29].
The disadvantage of fully active flapping foil is that the heaving and pitching motions require complex mechanical structures, which reduce the usability and reliability of the system [21]. This shortcoming can be alleviated by using semi-active flapping foil, which can reduce one degree of freedom to simplify the mechanical structure [30]. Semi-active flapping hydrofoil can only choose the heaving motion as the active motion and the pitching motion as the passive motion in order to obtain the thrust [20,[31][32][33]. As for the difference between the principle of semi-active and fully active flapping foils, the factors influencing the propulsive performance of semi-active flapping hydrofoil are not the same as those of the fully active flapping one. Additional influencing factors include the spring stiffness, added mass, and critical pitching amplitude of the pitching motion [20,32,34]. Conventional semi-active flapping foil has a rigid body, and some scholars also carried out research work on the hydrodynamic performance of flexible flapping hydrofoil [35,36]. Limiting the magnitude of the pitching motion may reduce the propulsion efficiency of flapping foil [26,33]. Therefore, this factor is rarely considered when semi-active flapping hydrofoil is used as a propeller. The heaving motion of the wave glider does not require an additional energy input. It is more concerned with the mean value of the propulsion power and is not sensitive to propulsion efficiency. Based on this feature, this paper studies the propulsive performance of wave glider flapping foil and further explores how to improve it.
The purpose of this paper is to study the wave driving mechanism of the wave glider. Furthermore, we analyze the propulsive characteristics of semi-active flapping foil, which is driven by sinusoidal waves. The heaving motion of semi-active flapping foil is prescribed, and the pitching motion is determined by the hydrodynamic force and movement acting on the foil. In the present work, the influence of a reduced frequency, spring stiffness, and critical pitching amplitude on the propulsive performance of the wave glider are systematically investigated.

Problem Description
To study the wave propulsion mechanism of the wave glider, the coupling between the surface boat and the submerged glider is simplified into a semi-active flapping foil, as shown in Figure 2. c is the chord length of the foil, d is the thickness of the foil, and the pivot point is located at c/3 of chord. The mass ratio r is the density ratio between the foil and the fluid, and its value is kept constant in this study. The heaving motion of the semi-active flapping foil of the wave glider is driven by the waves from the ocean [3,37,38]. For the sake of the present research, the heaving motion is simplified into a harmonic motion. Ignoring the gravity and buoyancy of the hydrofoil, the pitching motion is completely determined by the hydrodynamic force and torsion spring. Therefore, this factor is rarely considered when semi-active flapping hydrofoil is used as a propeller. The heaving motion of the wave glider does not require an additional energy input. It is more concerned with the mean value of the propulsion power and is not sensitive to propulsion efficiency. Based on this feature, this paper studies the propulsive performance of wave glider flapping foil and further explores how to improve it. The purpose of this paper is to study the wave driving mechanism of the wave glider. Furthermore, we analyze the propulsive characteristics of semi-active flapping foil, which is driven by sinusoidal waves. The heaving motion of semi-active flapping foil is prescribed, and the pitching motion is determined by the hydrodynamic force and movement acting on the foil. In the present work, the influence of a reduced frequency, spring stiffness, and critical pitching amplitude on the propulsive performance of the wave glider are systematically investigated.

Problem description
To study the wave propulsion mechanism of the wave glider, the coupling between the surface boat and the submerged glider is simplified into a semi-active flapping foil, as shown in Figure 2. c is the chord length of the foil, d is the thickness of the foil, and the pivot point is located at c/3 of chord. The mass ratio r is the density ratio between the foil and the fluid, and its value is kept constant in this study. The heaving motion of the semi-active flapping foil of the wave glider is driven by the waves from the ocean [3,37,38]. For the sake of the present research, the heaving motion is simplified into a harmonic motion. Ignoring the gravity and buoyancy of the hydrofoil, the pitching motion is completely determined by the hydrodynamic force and torsion spring. The prescribed heaving motion ( ) and heaving velocity ( ) can be written as: where is the heaving amplitude, f is the heaving frequency, is the heaving phase( a constant value of zero is used in this study), is the free stream velocity, * is the reduced frequency, which is the non-dimensional heaving frequency, defined by * = / . The differential equation of the pitching motion of flapping foil can be modelled as follows: ( ) + ( ) + ( ) = ( ) (2) where is the moments of inertia of the foil, is the pitching damping factor, is the spring stiffness factor of the foil and * is the normalized spring constant of the spring stiffness, defined * = /( ). ( ) represents the total hydrodynamic torque acting on the hydrofoil. To consider the influence of the critical pitching amplitude, the pitching damping factor can be written as The prescribed heaving motion y(t) and heaving velocity . y(t) can be written as: where y 0 is the heaving amplitude, f is the heaving frequency, Φ is the heaving phase( a constant value of zero is used in this study), U ∞ is the free stream velocity, f * is the reduced frequency, which is the non-dimensional heaving frequency, defined by f * = f c/U ∞ . The differential equation of the pitching motion of flapping foil can be modelled as follows: 4 of 20 where I is the moments of inertia of the foil, ζ is the pitching damping factor, k is the spring stiffness factor of the foil and k * is the normalized spring constant of the spring stiffness, defined k * = k/ 1 2 ρU 2 ∞ . M(t) represents the total hydrodynamic torque acting on the hydrofoil. To consider the influence of the critical pitching amplitude, the pitching damping factor ζ can be written as where θ c is defined as the critical pitching amplitude, which can be used to limit the pitching amplitude [36]. The structure of the critical pitching amplitude can be found in a previous work [5]. The effective angle of attack (AOA) α e (t) can be calculated by where α h (t) represents the heaving induced angle, when considering the influence of the upstream flow. The instantaneous input power P I (t) and the thrust power P O (t) can be calculated by where . θ(t) is defined as the pitching velocity, while F x (t) and F y (t) represent the horizontal and vertical hydrodynamic forces of the flapping foil, respectively. The dimensionless coefficients of the hydrodynamic forces and the moment are defined as where C X (t), C Y (t), and C M (t) are the dimensionless coefficients of the thrust, lift, and moment. ρ is the fluid density. The mean value of the input power coefficient C PI and the output power coefficient C PO are given by where the over-line represents the average value over one period. The propulsive efficiency η is written as

Numerical Method and Validation
We use the CFD software, FLUENT version 19.1, to simulate the flow around the semi-active hydrofoil, and the URANS equations are solved by the Finite Volume Method (FVM) [39,40]. As shown in Figure 3, a mesh for two-dimensional NACA0012 foil is produced, and the concepts of dynamic mesh and sliding interface are used in the computation. The vertical motion of the inner zone is used to simulate the prescribed heaving motion of the hydrofoil, and the rotation of the inner zone is used to simulate the pitching motion of the hydrofoil. The grid of the inner zone is fixed, and the grid of the outer zone is changed dynamically. The inner zone and the outer zone exchange data through the sliding interface. A Reynolds number of Re = U ∞ c/γ = 42000 is selected for all cases. To address the choice of turbulence model for the simulation of the semi-active flapping foil, the following options are used : Spalart-Allmaras, k − ω SST and k − ω SST with low-Re corrections. The comparison of their performance in the case under consideration is shown in Appendix A. The one-equation Spalart-Allmaras turbulence model is chosen for the simulation of semi-action flapping foil, considering the benefits of a lower requirement in terms of the quality of the mesh and a reliable result [39]. The SIMPLE algorithm is selected for the pressure-velocity coupling. The second-order scheme is used for the pressure, momentum, and turbulent viscosity resolution. The second-order implicit scheme is used for the unsteady formulation. An absolute convergence criterion of 10 −5 is set for the x-velocity, y-velocity, continuity, and nut. Further simulation parameter settings are described in detail in a previous study [5]. Based on the dynamic-hydrodynamic coupling method, the pitching motion can be calculated and controlled by User Defined Functions (UDFs) based on the FLUENT. The heaving motion is prescribed and controlled by UDFs. The flow chart of the dynamic-hydrodynamic coupling method for semi-active flapping foil is shown in Figure 4. First, the initial conditions are defined (the heaving and pitching motion and physical parameters of the foil). Secondly, the differential Equation (2) is estimated using the Newton-Euler method (fourth-order Runge-Kutta method) to obtain the pitching angle and velocity of oscillating foil. Then, depending on Equation (1), the position and velocity of the heaving motion are updated. We complete the numerical simulation process by repeating the above steps [20]. To validate the simulation strategy, grid independence studies are carried out using four different grid levels (Coarse grid 24000 cells, Medium-1 grid 48000 cells, Medium-2 grid 64000 cells, and Fine grid 96000 cells), as well as three different corresponding time-steps (1000 ts/2000 ts/5000 ts/10000 ts). A periodic solution for simulation requires 10-20 flapping periods, depending on the initial conditions. The stopping criterion is that the variation of the mean value of the thrust coefficient is less than 1% between successive periods. The evolution of the thrust coefficient for the independence tests is shown in Figure 5. The mean value and mean error of the thrust coefficient for the different grid levels and time-steps are provided in Table 1. Considering the result quality and computation efficiency, the medium-1 grid and medium-2 time-step are selected for the simulations. To further verify the credibility of the simulation strategy, the present research is validated against the previous study of Kinsey and Dumas under the same conditions (fully active flapping foil) [39]. The time variations over one cycle of the thrust and lift coefficients are illustrated in Figure 6. The result indicates that the present study agrees well with Kinsey and Dumas's research, which confirms the credibility of the simulation strategy. To validate the simulation strategy, grid independence studies are carried out using four different grid levels (Coarse grid 24000 cells, Medium-1 grid 48000 cells, Medium-2 grid 64000 cells, and Fine grid 96000 cells), as well as three different corresponding time-steps (1000 ts/2000 ts/5000 ts/10000 ts). A periodic solution for simulation requires 10-20 flapping periods, depending on the initial conditions. The stopping criterion is that the variation of the mean value of the thrust coefficient is less than 1% between successive periods. The evolution of the thrust coefficient for the independence tests is shown in Figure 5. The mean value and mean error of the thrust coefficient for the different grid levels and time-steps are provided in Table 1. Considering the result quality and computation efficiency, the medium-1 grid and medium-2 time-step are selected for the simulations.    To further verify the credibility of the simulation strategy, the present research is validated against the previous study of Kinsey and Dumas under the same conditions (fully active flapping foil) [39]. The time variations over one cycle of the thrust and lift coefficients are illustrated in Figure 6. The result indicates that the present study agrees well with Kinsey and Dumas's research, which confirms the credibility of the simulation strategy.

Results and discussion
In this section, the propulsive performance of semi-active flapping foil is discussed in detail. The effects of reduced frequency, spring stiffness and critical pitching amplitude on the hydrodynamic characteristics of the foil are further systematically investigated. The simulation parameters are selected, depending on the structure of the wave glider and the wave characteristics in typical sea conditions [5,7], which are illustrated in Table 2.

Results and Discussion
In this section, the propulsive performance of semi-active flapping foil is discussed in detail. The effects of reduced frequency, spring stiffness and critical pitching amplitude on the hydrodynamic characteristics of the foil are further systematically investigated. The simulation parameters are selected, depending on the structure of the wave glider and the wave characteristics in typical sea conditions [5,7], which are illustrated in Table 2.

The Influence of the Reduced Frequency
First, the influence of the reduced frequency on the propulsive characteristic of semi-active flapping foil is studied. The simulation parameters of y 0 = 0.2, 0.4, k = 0.05-0.2, and f * = 0.1-0.68 are selected. The mean value of the output power coefficients, against the variation of the reduced frequency, is shown in Figure 7. It is found that the mean value of the output power coefficients increases monotonically and significantly for y 0 = 0.2 by increasing the reduced frequency, as shown in Figure 7a. The same growth law for y 0 = 0.4 can be found, as shown in Figure 7b. We also find that there is a critical reduced frequency: When the reduced frequency is lower than the critical reduced frequency, the mean value of the output power coefficient increases exponentially, and when the reduced frequency is higher than the critical reduced frequency, the semi-active flapping foil cannot acquire a steady thrust. As the reduced frequency increases, the pitching amplitude of the flapping foil increase under the action of hydrodynamic forces. Excessive pitching amplitude prevents the hydrofoil producing effective propulsion. This is one reason to choose a reduced frequency range of 0.1-0.68. Meanwhile, we investigate the effect of the critical pitching amplitude on the propulsive performance to improve the adaptation of the flapping foil at a higher reduced frequency.

The influence of the reduced frequency
First, the influence of the reduced frequency on the propulsive characteristic of semi-active flapping foil is studied. The simulation parameters of = 0.2, 0.4, = 0.05 − 0.2, and * = 0.1 − 0.68 are selected. The mean value of the output power coefficients, against the variation of the reduced frequency, is shown in Figure 7. It is found that the mean value of the output power coefficients increases monotonically and significantly for = 0.2 by increasing the reduced frequency, as shown in Figure 7a. The same growth law for = 0.4 can be found, as shown in Figure 7b. We also find that there is a critical reduced frequency: When the reduced frequency is lower than the critical reduced frequency, the mean value of the output power coefficient increases exponentially, and when the reduced frequency is higher than the critical reduced frequency, the semi-active flapping foil cannot acquire a steady thrust. As the reduced frequency increases, the pitching amplitude of the flapping foil increase under the action of hydrodynamic forces. Excessive pitching amplitude prevents the hydrofoil producing effective propulsion. This is one reason to choose a reduced frequency range of 0.1-0.68. Meanwhile, we investigate the effect of the critical pitching amplitude on the propulsive performance to improve the adaptation of the flapping foil at a higher reduced frequency. The propulsive efficiency against the reduced frequency is illustrated in Figure 8. The propulsive efficiency increases as the frequency is reduced, except in the cases = 0.2, * = 0.27, and = 0.4, * = 0.09.  The propulsive efficiency against the reduced frequency is illustrated in Figure 8. The propulsive efficiency increases as the frequency is reduced, except in the cases y 0 = 0.2, k * = 0.27, and y 0 = 0.4, k * = 0.09. To illustrate the effect of the reduced frequency on the propulsive performance, the time evolution of the thrust coefficient and pitching motion over one cycle are shown in Figure 9. The variation of the thrust coefficient in the case of * = 0.17 is much lower than that of * = 0.68. The maximum values of the instantaneous thrust coefficient are 15.87 and 0.24, respectively. To make it visible, the time variation of the thrust coefficient of * = 0.17 was locally amplified, as shown in Figure 9a. Therefore, it is evident that the reduced frequency has a significant effect on the propulsion performance of the flapping hydrofoil, and the same conclusion can be found in [41,42]. We also note To illustrate the effect of the reduced frequency on the propulsive performance, the time evolution of the thrust coefficient and pitching motion over one cycle are shown in Figure 9. The variation of the thrust coefficient in the case of f * = 0.17 is much lower than that of f * = 0.68. The maximum values of the instantaneous thrust coefficient are 15.87 and 0.24, respectively. To make it visible, the time variation of the thrust coefficient of f * = 0.17 was locally amplified, as shown in Figure 9a. Therefore, it is evident that the reduced frequency has a significant effect on the propulsion performance of the flapping hydrofoil, and the same conclusion can be found in [41,42]. We also note that there is a phase difference between the evolution of the thrust coefficient for f * = 0.17 and f * = 0.68. The maximum thrust coefficients of f * = 0.17 and f * = 0.68 occur at t/T = 0 and t/T = 0.1, respectively. It should be noted that due to the huge difference in the hydrodynamic force and torque action on flapping foil, there is a significant difference between the pitching motion profile and pitching amplitude between f * = 0.17 and f * = 0.68. The pitching amplitudes of f * = 0.17 and f * = 0.68 are ±14.5 • and ±56 • , respectively, as shown in Figure 9b. To illustrate the effect of the reduced frequency on the propulsive performance, the time evolution of the thrust coefficient and pitching motion over one cycle are shown in Figure 9. The variation of the thrust coefficient in the case of * = 0.17 is much lower than that of * = 0.68. The maximum values of the instantaneous thrust coefficient are 15.87 and 0.24, respectively. To make it visible, the time variation of the thrust coefficient of * = 0.17 was locally amplified, as shown in Figure 9a. Therefore, it is evident that the reduced frequency has a significant effect on the propulsion performance of the flapping hydrofoil, and the same conclusion can be found in [41,42]. We also note that there is a phase difference between the evolution of the thrust coefficient for * = 0.17 and * = 0.68. The maximum thrust coefficients of * = 0.17 and * = 0.68 occur at / = 0 and / = 0.1, respectively. It should be noted that due to the huge difference in the hydrodynamic force and torque action on flapping foil, there is a significant difference between the pitching motion profile and pitching amplitude between * = 0.17 and * = 0.68. The pitching amplitudes of * = 0.17 and * = 0.68 are ±14.5° and ±56°, respectively, as shown in Figure 9b.  Figure 10a. Thus, there is a negative pressure region near the lower surface. Since the foil moves upward, and the heaving velocity is relatively large, a positive pressure region appears on the upper surface. The distribution of the negative and positive pressure regions is shown in Figure 10b, and there exists a large beneficial effective angle of attack (AOA) at that moment, as shown in Figure 9b. Therefore, the maximum thrust coefficient is generated at / = 0.1,  Figure 10a. Thus, there is a negative pressure region near the lower surface. Since the foil moves upward, and the heaving velocity is relatively large, a positive pressure region appears on the upper surface. The distribution of the negative and positive pressure regions is shown in Figure 10b, and there exists a large beneficial effective angle of attack (AOA) at that moment, as shown in Figure 9b. Therefore, the maximum thrust coefficient is generated at t/T = 0.1, as illustrated in Figure 9a. As shown in Figure 10a, the trailing vortex of the foil sheds, and the wake of the foil appears as a reversed Karman Vortex Street. The leading-edge vortex (LEV) and reversed Karman vortex street are beneficial for improving the thrust performance of flapping foil [18,24]. Comparing Figure 10c with Figure 10a, it is found that the vorticity intensity of f * = 0.17 at t/T = 0 is weaker than that of f * = 0.68 at t/T = 0.1, and the wake region of the hydrofoil is rapidly dissipated. The corresponding pressure contours are presented in Figure 10d. Since the leading-edge vortex (LEV) and trailing-edge vortex (TEV) are formed on the lower surface of the foil, a negative pressure region is generated. Due to the upward motion of the foil, a positive pressure region is formed in the upper surface. Comparing Figure 10d with Figure 10b, we find that the intensity of pressure of the case f * = 0.17 is significantly weaker than that of f * = 0.68, and the effective angle of attack (AOA) is lower. These factors further weaken the propulsion of the flapping hydrofoil.
dissipated. The corresponding pressure contours are presented in Figure 10d. Since the leading-edge vortex (LEV) and trailing-edge vortex (TEV) are formed on the lower surface of the foil, a negative pressure region is generated. Due to the upward motion of the foil, a positive pressure region is formed in the upper surface. Comparing Figure 10d with Figure 10b, we find that the intensity of pressure of the case * = 0.17 is significantly weaker than that of * = 0.68, and the effective angle of attack (AOA) is lower. These factors further weaken the propulsion of the flapping hydrofoil.

The influence of the spring stiffness
According to the previous section, the reduced frequency has a significant impact on the hydrodynamic characteristics of semi-active flapping foil. As shown in Figures 7 and 8, it is obvious that the stiffness of the torsion spring affects the propulsion performance of the hydrofoil. In this section, the influence of the spring is investigated in more detail. The comparison of the mean value of the output power coefficient and propulsive efficiency for different torsion spring stiffness is illustrated in Figure 11

The Influence of the Spring Stiffness
According to the previous section, the reduced frequency has a significant impact on the hydrodynamic characteristics of semi-active flapping foil. As shown in Figures 7 and 8, it is obvious that the stiffness of the torsion spring affects the propulsion performance of the hydrofoil. In this section, the influence of the spring is investigated in more detail. The comparison of the mean value of the output power coefficient and propulsive efficiency for different torsion spring stiffness is illustrated in Figure 11. When the heaving amplitude and reduced frequency are relatively large (y 0 = 0.2, f * = 0.68 and y 0 = 0.4, f * = 0.34), the mean value of the output power coefficient firstly increases and then decreases, when the spring stiffness further increases, as shown in Figure 11a. When the heaving amplitude and reduced frequency decrease (y 0 = 0.2, f * < 0.68 and y 0 = 0.4, f * < 0.34), the mean value of the output power coefficient decreases monotonically as the spring stiffness increases. Figure 11b shows that except in the case of k * = 0.09, y 0 = 0.4, f * = 0.34, the propulsive efficiency decreases monotonically as the spring stiffness increases. To further investigate the above phenomenon, the influence of the torsion spring on two cases of the performance of the flapping foil is studied, namely for a larger reduced frequency (y 0 = 0.2, f * = 0.68) and smaller one (y 0 = 0.2, f * = 0.17). The time evolution over one cycle of the thrust, lift and moment coefficients, and pitching profile, with different torsion spring stiffness (k * = 0.09, 0.27, 0.37), are shown in Figure 12a-d. There exists a phase difference in the evolution curve of the thrust coefficient of the 1/10 cycle between k * = 0.09 and k * = 0.27, 0.37. It is noted that the maximum value of the thrust coefficient appears at t/T = 0.1 and t/T = 0.6, t/T = 0 and t/T = 0.5, respectively, as shown in Figure 12a. The maximum instantaneous thrust coefficient of k * = 0.27 is greater than that of k * = 0.09 and k * = 0.37. Thus, the mean value of the output power coefficient of k * = 0.27 is larger than that of k * = 0.09 and k * = 0.37, as shown in Figure 11a. We also find that the flapping foil of k * = 0.09 experiences a large drag force at t/T = 0.3-0.4. The reason for this phenomenon is that the heaving velocity decreases to zero at this time interval, and an overshoot of pitching motion of the flapping foil occurs due to the relatively soft spring stiffness, as shown in Figure 12d. When the torsion spring stiffness increases, the flapping foil can quickly be brought back to a horizontal plane, thereby suppressing the overshoot of the pitching motion. Furthermore, we note that as the torsion spring stiffness coefficient increases, the lift and moment coefficients increase significantly. The increase in the amplitude of the moment coefficient is higher than that of the thrust coefficient, which leads to a decrease in the propulsive efficiency, as shown in Figure 12b,c.
(a) (b) Figure 11. Comparison of (a) the mean value of the output power coefficient and (b) the propulsive efficiency for varying spring stiffness at = 0.2 − 0.4, * = 0.09 − 0.68.
The time evolution over one cycle of the thrust, lift and moment coefficients, and pitching profile, with different torsion spring stiffness ( * = 0.09, 0.27, 0.37), are shown in Figure 12a-d. There exists a phase difference in the evolution curve of the thrust coefficient of the 1/10 cycle between * = 0.09 and * = 0.27, 0.37. It is noted that the maximum value of the thrust coefficient appears at / = 0.1 and / = 0.6 , / = 0 and / = 0.5 , respectively, as shown in Figure 12a. The maximum instantaneous thrust coefficient of * = 0.27 is greater than that of * = 0.09 and * = 0.37. Thus, the mean value of the output power coefficient of * = 0.27 is larger than that of * = 0.09 and * = 0.37, as shown in Figure 11a. We also find that the flapping foil of * = 0.09 experiences a large drag force at / = 0.3 − 0.4. The reason for this phenomenon is that the heaving velocity decreases to zero at this time interval, and an overshoot of pitching motion of the flapping foil occurs due to the relatively soft spring stiffness, as shown in Figure 12d. When the torsion spring stiffness increases, the flapping foil can quickly be brought back to a horizontal plane, thereby suppressing the overshoot of the pitching motion. Furthermore, we note that as the torsion spring stiffness coefficient increases, the lift and moment coefficients increase significantly. The increase in the amplitude of the moment coefficient is higher than that of the thrust coefficient, which leads to a decrease in the propulsive efficiency, as shown in Figure 12b and 12c. To further understand the influence of the torsion spring stiffness on the flow around the flapping foil, the five typical instant vorticity contours and corresponding pressure contours for different spring stiffness ( * = 0.09, * = 0.27) are shown in Figure 13 and Figure 14. It is noted that a leading-edge vortex (LEV) is formed on the lower surface of the foil at different instants, i.e., t/T = 0, 0.1 for * = 0.09, * = 0.27, respectively. A negative region is found on the lower surface of the foil, as shown in Figure 14a and 14b. At this time, the foil moves upward, and a positive pressure region is found on the upper surface of the foil. Thus, the maximum value of the thrust coefficient is produced at this moment, as shown in Figure 12a. Comparing Figure 14a with Figure 14b, we find that due to the different intensities of the leading-edge vortex (LEV) and the trailing-edge vortex (TEV), the intensity of the negative pressure zone, formed on the lower surface of the foil for * = 0.27, is higher than that formed on the lower surface of the foil for * = 0.09; therefore, the maximum value of the thrust coefficient for * = 0.27 is larger than that for * = 0.09. Owing to the effect of the tough torsion spring ( * = 0.27), a new trailing-edge vortex (TEV) is generated and separated at To further understand the influence of the torsion spring stiffness on the flow around the flapping foil, the five typical instant vorticity contours and corresponding pressure contours for different spring stiffness (k * = 0.09, k * = 0.27) are shown in Figures 13 and 14. It is noted that a leading-edge vortex (LEV) is formed on the lower surface of the foil at different instants, i.e., t/T = 0, 0.1 for k * = 0.09, k * = 0.27, respectively. A negative region is found on the lower surface of the foil, as shown in Figure 14a,b. At this time, the foil moves upward, and a positive pressure region is found on the upper surface of the foil. Thus, the maximum value of the thrust coefficient is produced at this moment, as shown in Figure 12a. Comparing Figure 14a with Figure 14b, we find that due to the different intensities of the leading-edge vortex (LEV) and the trailing-edge vortex (TEV), the intensity of the negative pressure zone, formed on the lower surface of the foil for k * = 0.27, is higher than that formed on the lower surface of the foil for k * = 0.09; therefore, the maximum value of the thrust coefficient for k * = 0.27 is larger than that for k * = 0.09. Owing to the effect of the tough torsion spring (k * = 0.27), a new trailing-edge vortex (TEV) is generated and separated at t/T = 0.35, and a corresponding negative pressure zone is created at the tail of the upper surface of the foil, as shown in Figures 13e and  14e. For a low spring stiffness (k * = 0.09), no negative pressure zone is formed at the tail of the upper surface of the foil. Thus, the tough torsion spring improves the propulsive performance of the flapping foil at t/T = 0.35. Above all, we conclude that under greater reduced frequency conditions, proper spring stiffness can increase the maximum value of the instantaneous thrust coefficient and suppress the drag force caused by the pitching overshoot.
flapping foil, the five typical instant vorticity contours and corresponding pressure contours for different spring stiffness ( * = 0.09, * = 0.27) are shown in Figure 13 and Figure 14. It is noted that a leading-edge vortex (LEV) is formed on the lower surface of the foil at different instants, i.e., t/T = 0, 0.1 for * = 0.09, * = 0.27, respectively. A negative region is found on the lower surface of the foil, as shown in Figure 14a and 14b. At this time, the foil moves upward, and a positive pressure region is found on the upper surface of the foil. Thus, the maximum value of the thrust coefficient is produced at this moment, as shown in Figure 12a. Comparing Figure 14a with Figure 14b, we find that due to the different intensities of the leading-edge vortex (LEV) and the trailing-edge vortex (TEV), the intensity of the negative pressure zone, formed on the lower surface of the foil for * = 0.27, is higher than that formed on the lower surface of the foil for * = 0.09; therefore, the maximum value of the thrust coefficient for * = 0.27 is larger than that for * = 0.09. Owing to the effect of the tough torsion spring ( * = 0.27), a new trailing-edge vortex (TEV) is generated and separated at t/T = 0.35, and a corresponding negative pressure zone is created at the tail of the upper surface of the foil, as shown in Figure 13e and Figure 14e. For a low spring stiffness ( * = 0.09), no negative pressure zone is formed at the tail of the upper surface of the foil. Thus, the tough torsion spring improves the propulsive performance of the flapping foil at t/T = 0.35. Above all, we conclude that under greater reduced frequency conditions, proper spring stiffness can increase the maximum value of the instantaneous thrust coefficient and suppress the drag force caused by the pitching overshoot.    Figure 15. It is pointed out that the maximum thrust coefficient of * = 0.09 is higher than that of * = 0.27 and * = 0.37, while the instantaneous thrust coefficient of * = 0.09 is lower than the other two cases at almost all times. As shown in Figure 15b, due to a weaker torque spring ( * = 0.09), the effective angle of attack (AOA) is more sensitive to the hydrodynamic forces and movement of the flapping foil, and the pitching amplitude is higher than that of * = 0.27 and * = 0.37. According to Figure 15a, five typical snapshots in the range at t/T = 0 -0.4 are selected. The instantaneous pressure contours are illustrated in Figure 16. It could be observed that the intensity of the positive and negative pressure regions, formed on the upper and lower surfaces of the foil for * = 0.37, is slightly higher than that of * = 0.09. However, the pitching angle seems to play a more influential role in the thrust generation of the flapping foil for a lower reduced frequency. Figure 11a also illustrates that for a lower reduced frequency (y 0 = 0.2, f * < 0.68 and y 0 = 0.4, f * < 0.34), the propulsive performance of the flapping foil decreases as the spring stiffness increases, and both the mean values of the output power coefficient and the propulsive efficiency decrease. The time history of the thrust coefficient and pitching profile, during one period, is illustrated in Figure 15. It is pointed out that the maximum thrust coefficient of k * = 0.09 is higher than that of k * = 0.27 and k * = 0.37, while the instantaneous thrust coefficient of k * = 0.09 is lower than the other two cases at almost all times. As shown in Figure 15b, due to a weaker torque spring (k * = 0.09), the effective angle of attack (AOA) is more sensitive to the hydrodynamic forces and movement of the flapping foil, and the pitching amplitude is higher than that of k * = 0.27 and k * = 0.37. According to Figure 15a, five typical snapshots in the range at t/T = 0-0.4 are selected. The instantaneous pressure contours are illustrated in Figure 16. It could be observed that the intensity of the positive and negative pressure regions, formed on the upper and lower surfaces of the foil for k * = 0.37, is slightly higher than that of k * = 0.09. However, the pitching angle seems to play a more influential role in the thrust generation of the flapping foil for a lower reduced frequency. illustrated in Figure 15. It is pointed out that the maximum thrust coefficient of = 0.09 is higher than that of * = 0.27 and * = 0.37, while the instantaneous thrust coefficient of * = 0.09 is lower than the other two cases at almost all times. As shown in Figure 15b, due to a weaker torque spring ( * = 0.09), the effective angle of attack (AOA) is more sensitive to the hydrodynamic forces and movement of the flapping foil, and the pitching amplitude is higher than that of * = 0.27 and * = 0.37. According to Figure 15a, five typical snapshots in the range at t/T = 0 -0.4 are selected. The instantaneous pressure contours are illustrated in Figure 16. It could be observed that the intensity of the positive and negative pressure regions, formed on the upper and lower surfaces of the foil for * = 0.37, is slightly higher than that of * = 0.09. However, the pitching angle seems to play a more influential role in the thrust generation of the flapping foil for a lower reduced frequency.

The influence of the critical pitching amplitude
According to the results of the previous section, the stiffness of the torsion spring has different effects at different reduced frequencies. However, for the wave propeller of the wave glider, the torsion spring set on the hydrofoil cannot be changed according to different sea conditions, once it is selected. Therefore, we prefer to choose a unique spring stiffness to ensure a higher propulsion performance of the flapping foil at different reduced frequencies. According to Figure 12 and Figure  13, it is found that the decrease of the hydrofoil propulsion performance is related to the overshoot of the pitching motion, when * = 0.09. Therefore, we wish to improve the overshoot of the hydrofoil by setting the critical pitching amplitude and thereby improving the propulsive performance of the hydrofoil under the soft spring stiffness coefficient.
The comparison of the mean value of output power coefficients and propulsive efficiency for varying critical pitching amplitudes is presented in Figure 17. The selection of the critical pitching amplitudes is determined according to the magnitude of the pitching motion of the flapping foil under a free condition. It is found that the mean value of the output power coefficient first increases as the critical pitching amplitude increases ( ≤ 50) and then basically does not change ( > 50). When the critical pitching amplitude is greater than 20, the mean value of the output power coefficient is already higher than that under a free condition, as shown in Figure 17a. Moreover, compared with Figure 11a, the mean value of the output power coefficient for = 0.2, = 50, * = 0.09, * = 0.68 is higher than the maximum value in the free state. However, as shown in Figure 17b, the propulsive efficiency increases with the increase of the critical pitching amplitude, except in the case of = 40. The maximum value of the propulsive efficiency, considering the critical pitching amplitude, is less than that under the free condition for = 0.2, * = 0.09. On the other hand, the maximum value of the propulsive efficiency is slightly higher than that under the

The Influence of the Critical Pitching Amplitude
According to the results of the previous section, the stiffness of the torsion spring has different effects at different reduced frequencies. However, for the wave propeller of the wave glider, the torsion spring set on the hydrofoil cannot be changed according to different sea conditions, once it is selected. Therefore, we prefer to choose a unique spring stiffness to ensure a higher propulsion performance of the flapping foil at different reduced frequencies. According to Figures 12 and 13, it is found that the decrease of the hydrofoil propulsion performance is related to the overshoot of the pitching motion, when k * = 0.09. Therefore, we wish to improve the overshoot of the hydrofoil by setting the critical pitching amplitude and thereby improving the propulsive performance of the hydrofoil under the soft spring stiffness coefficient.
The comparison of the mean value of output power coefficients and propulsive efficiency for varying critical pitching amplitudes is presented in Figure 17. The selection of the critical pitching amplitudes is determined according to the magnitude of the pitching motion of the flapping foil under a free condition. It is found that the mean value of the output power coefficient first increases as the critical pitching amplitude increases (θ c ≤ 50) and then basically does not change (θ c > 50). When the critical pitching amplitude is greater than 20, the mean value of the output power coefficient is already higher than that under a free condition, as shown in Figure 17a. Moreover, compared with Figure 11a, the mean value of the output power coefficient for y 0 = 0.2, θ c = 50, k * = 0.09, f * = 0.68 is higher than the maximum value in the free state. However, as shown in Figure 17b, the propulsive efficiency increases with the increase of the critical pitching amplitude, except in the case of θ c = 40. The maximum value of the propulsive efficiency, considering the critical pitching amplitude, is less than that under the free condition for y = 0.2, k * = 0.09. On the other hand, the maximum value of the propulsive efficiency is slightly higher than that under the free condition at y = 0.4, k * = 0.09. It is found that the critical pitching amplitude can improve the thrust power of flapping foil. The time history over one cycle of the thrust coefficient and the pitching profile for = NULL (the pitching amplitude is unlimited) and = 50 are shown in Figure 18. For = 50 , the instantaneous thrust coefficient has two peaks in half a cycle. Conversely, for = NULL, there is only one peak, as shown in Figure 18a. The first peak for = 50 appears at t/T = 0.32. Meanwhile, the pitching motion reaches the critical pitching amplitude, as shown in Figure 18b. Moreover, it is also noted that the second peak for = 50 is slightly larger than the peak for = NULL. Figure   19b also shows that due to the influence of the critical pitching amplitude, the pitching profile is modified, and the magnitude of the pitching motion decreases to 50 .
As shown in the partially enlarged view in Figure 18a, there is a drastic change in the instantaneous thrust coefficient at t/T = 0.32 − 0.36 . Therefore, the five corresponding instantaneous vorticity and pressure contours are illustrated in Figure 19 and Figure 20. Comparing the vorticity contours for = 50 with that for = NULL, a new trailing-edge vortex (TEV) is generated at t/T = 0.34, and a corresponding negative pressure zone is created at the tail of the upper surface of the foil, as shown in Figure 19c and Figure 20c. Due to the influence of the trailing-edge vortex (TEV), the magnitude of the second peak of the thrust coefficient for = 50 is slightly enhanced. Since the pitching motion of the flapping foil is limited by the critical pitching amplitude, a strong positive pressure region forms on the lower surface of the foil, and a negative pressure region forms on the upper surface of the foil, as shown in Figure 20b. Thus, a large magnitude of the thrust The time history over one cycle of the thrust coefficient and the pitching profile for θ c = NULL (the pitching amplitude is unlimited) and θ c = 50 • are shown in Figure 18. For θ c = 50 • , the instantaneous thrust coefficient has two peaks in half a cycle. Conversely, for θ c = NULL, there is only one peak, as shown in Figure 18a. The first peak for θ c = 50 • appears at t/T = 0.32. Meanwhile, the pitching motion reaches the critical pitching amplitude, as shown in Figure 18b. Moreover, it is also noted that the second peak for θ c = 50 • is slightly larger than the peak for θ c = NULL. Figure 19b also shows that due to the influence of the critical pitching amplitude, the pitching profile is modified, and the magnitude of the pitching motion decreases to 50 • . The time history over one cycle of the thrust coefficient and the pitching profile for = NULL (the pitching amplitude is unlimited) and = 50 are shown in Figure 18. For = 50 , the instantaneous thrust coefficient has two peaks in half a cycle. Conversely, for = NULL, there is only one peak, as shown in Figure 18a. The first peak for = 50 appears at t/T = 0.32. Meanwhile, the pitching motion reaches the critical pitching amplitude, as shown in Figure 18b. Moreover, it is also noted that the second peak for = 50 is slightly larger than the peak for = NULL. Figure   19b also shows that due to the influence of the critical pitching amplitude, the pitching profile is modified, and the magnitude of the pitching motion decreases to 50 . As shown in the partially enlarged view in Figure 18a, there is a drastic change in the instantaneous thrust coefficient at t/T = 0.32 − 0.36 . Therefore, the five corresponding instantaneous vorticity and pressure contours are illustrated in Figure 19 and  To further explore the influence of the critical pitching amplitude on the propulsive performance of the flapping foil, the evolution of the instantaneous pressure contours from 0.324T to 0.34T is presented in Figure 21a-i. We find that when the pitching motion reaches the critical pitching amplitude, the positive pressure region and the negative pressure region are formed on the upper and lower surface of the foil, as shown in Figure 21b. Then, the intensity of the pressure region is enhanced with the downward movement of the foil. The effect of the critical pitching amplitude is gradually reduced after t/T=0.332. The change process of the pressure contours corresponds to the evolution of the thrust coefficient, as shown in Figure 18a and Figure 21, and the influence of the critical pitching amplitude on the propulsion performance of flapping foil is clearly revealed. To further explore the influence of the critical pitching amplitude on the propulsive performance of the flapping foil, the evolution of the instantaneous pressure contours from 0.324T to 0.34T is presented in Figure 21a-i. We find that when the pitching motion reaches the critical pitching amplitude, the positive pressure region and the negative pressure region are formed on the upper and lower surface of the foil, as shown in Figure 21b. Then, the intensity of the pressure region is enhanced with the downward movement of the foil. The effect of the critical pitching amplitude is gradually reduced after t/T=0.332. The change process of the pressure contours corresponds to the evolution of the thrust coefficient, as shown in Figure 18a and Figure 21, and the influence of the critical pitching amplitude on the propulsion performance of flapping foil is clearly revealed. To further explore the influence of the critical pitching amplitude on the propulsive performance of the flapping foil, the evolution of the instantaneous pressure contours from 0.324T to 0.34T is presented in Figure 21a-i. We find that when the pitching motion reaches the critical pitching amplitude, the positive pressure region and the negative pressure region are formed on the upper and lower surface of the foil, as shown in Figure 21b. Then, the intensity of the pressure region is enhanced with the downward movement of the foil. The effect of the critical pitching amplitude is gradually reduced after t/T=0.332. The change process of the pressure contours corresponds to the evolution of the thrust coefficient, as shown in Figures 18a and 21, and the influence of the critical pitching amplitude on the propulsion performance of flapping foil is clearly revealed.

Conclusions
To study the propulsion performance of the wave glider's wave propeller, the propeller is simplified into a semi-active flapping hydrofoil, where the heaving motion is prescribed and the pitching motion is determined by the hydrodynamic force and movement. In the present study, the influence of the reduced frequency, spring stiffness, and critical pitching amplitude on the hydrodynamic characteristics of semi-active flapping foil are systematically investigated. It is noted that there is a critical value of reduced frequency: When the reduced frequency is lower than the critical value, the mean value of the output power coefficient increases exponentially, and when the reduced frequency is higher than the critical value, the semi-active flapping foil cannot acquire a steady thrust. Under different reduced frequencies, the spring stiffness has a different effect on the performance of the foil. For a higher reduced frequency, there is an optimal value of spring stiffness, which corresponds to the maximum value of the output power coefficient. As the spring stiffness increases, the mean value of the output power coefficient increases first and then decreases. For a lower reduced frequency, the mean value of the output power coefficient decreases monotonically as the spring stiffness increases. Except in special cases, concerning either higher or lower reduced frequencies, the propulsive efficiency of flapping foil monotonically decreases as the spring stiffness increases. Finally, we investigate the influence of the critical pitching amplitude on the propulsive performance of flapping foil and find that by using the critical pitching amplitude, the mean value of output power coefficient and the propulsive efficiency of semi-active flapping foil are improved, especially for higher heaving amplitudes. This is an important factor for improving the propulsive performance of the wave glider under high sea conditions.
The current study reveals the influence of the reduced frequency, spring stiffness and critical pitching amplitude on the propulsion performance of flapping hydrofoil and provides theoretical guidance for the optimal design of the wave glider. In the present research, the influence of damping and foil shape are neglected, and this problem requires further study in connection with the simulation of the propulsive characteristics of flapping foil.

Conclusions
To study the propulsion performance of the wave glider's wave propeller, the propeller is simplified into a semi-active flapping hydrofoil, where the heaving motion is prescribed and the pitching motion is determined by the hydrodynamic force and movement. In the present study, the influence of the reduced frequency, spring stiffness, and critical pitching amplitude on the hydrodynamic characteristics of semi-active flapping foil are systematically investigated. It is noted that there is a critical value of reduced frequency: When the reduced frequency is lower than the critical value, the mean value of the output power coefficient increases exponentially, and when the reduced frequency is higher than the critical value, the semi-active flapping foil cannot acquire a steady thrust. Under different reduced frequencies, the spring stiffness has a different effect on the performance of the foil. For a higher reduced frequency, there is an optimal value of spring stiffness, which corresponds to the maximum value of the output power coefficient. As the spring stiffness increases, the mean value of the output power coefficient increases first and then decreases. For a lower reduced frequency, the mean value of the output power coefficient decreases monotonically as the spring stiffness increases. Except in special cases, concerning either higher or lower reduced frequencies, the propulsive efficiency of flapping foil monotonically decreases as the spring stiffness increases. Finally, we investigate the influence of the critical pitching amplitude on the propulsive performance of flapping foil and find that by using the critical pitching amplitude, the mean value of output power coefficient and the propulsive efficiency of semi-active flapping foil are improved, especially for higher heaving amplitudes. This is an important factor for improving the propulsive performance of the wave glider under high sea conditions.
The current study reveals the influence of the reduced frequency, spring stiffness and critical pitching amplitude on the propulsion performance of flapping hydrofoil and provides theoretical guidance for the optimal design of the wave glider. In the present research, the influence of damping and foil shape are neglected, and this problem requires further study in connection with the simulation of the propulsive characteristics of flapping foil. Acknowledgments: The authors would like to thank the anonymous reviewers for their constructive comments and suggestions.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A
The turbulence model can impact the instantaneous forces and moment when flow separation and vortex shedding occur. To address the choice of turbulence model of simulation for semi-active flapping foil, we test the following turbulence models: Spalart-Allmaras, k − ω SST and k − ω SST with low-Re corrections. The simulation parameters of y 0 = 0.4, θ 0 = 50 • , and f * = 0.21 are selected. All numerical studies are repeated using the Fluent version 19.1 software.
The time variation over one cycle of the thrust coefficient and the lift coefficients is shown in Figure A1, and the corresponding vorticity contours are shown in Figure A2. The comparison of results from the three turbulence models for semi-active flapping foil is shown in Table A1. The instantaneous force coefficients are superposed nicely, and the simulations of different turbulence models reached satisfactory periodicity. We find that the Spalart-Allmaras simulation compares well with the SST simulation and the SST simulation with low-Re corrections, and the relative difference of the thrust coefficient between the Spalart-Allmaras and SST simulation with low-Re corrections is 1.8 %.
From the evolution of force coefficients of Figure A1 and vorticity contours of Figure A2, we note that the different characters of propulsive performance for different turbulence models are observed, those differences are concentrated near the top and bottom positions of heaving motion, where the thrust force is nearly minimal. These differences have little effect on the propulsive performance of the semi-active flapping foil, as shown in Table A1. The above results illustrate that the choice of turbulence model for this study slightly affects the propulsive performance of the semi-active flapping foil. Therefore, considering the benefits of a lower requirement in terms of the quality of the mesh and a reliable result, the turbulence model of Spalart-Allmaras is chosen for the simulation of semi-action flapping foil.

Acknowledgments:
The authors would like to thank the anonymous reviewers for their constructive comments and suggestions.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A
The turbulence model can impact the instantaneous forces and moment when flow separation and vortex shedding occur. To address the choice of turbulence model of simulation for semi-active flapping foil, we test the following turbulence models: Spalart-Allmaras, − SST and − SST with low-Re corrections. The simulation parameters of = 0.4 , = 50°, and * = 0.21 are selected. All numerical studies are repeated using the Fluent version 19.1 software.
The time variation over one cycle of the thrust coefficient and the lift coefficients is shown in Figure A1, and the corresponding vorticity contours are shown in Figure A2. The comparison of results from the three turbulence models for semi-active flapping foil is shown in Table A1. The instantaneous force coefficients are superposed nicely, and the simulations of different turbulence models reached satisfactory periodicity. We find that the Spalart-Allmaras simulation compares well with the SST simulation and the SST simulation with low-Re corrections, and the relative difference of the thrust coefficient between the Spalart-Allmaras and SST simulation with low-Re corrections is 1.8 %.
From the evolution of force coefficients of Figure A1 and vorticity contours of Figure A2, we note that the different characters of propulsive performance for different turbulence models are observed, those differences are concentrated near the top and bottom positions of heaving motion, where the thrust force is nearly minimal. These differences have little effect on the propulsive performance of the semi-active flapping foil, as shown in Table A1. The above results illustrate that the choice of turbulence model for this study slightly affects the propulsive performance of the semiactive flapping foil. Therefore, considering the benefits of a lower requirement in terms of the quality of the mesh and a reliable result, the turbulence model of Spalart-Allmaras is chosen for the simulation of semi-action flapping foil.