Numerical Analysis of the Effect of Flexibility on the Propulsive Performance of a Heaving Hydrofoil Undergoing Sinusoidal and Non-Sinusoidal Motions


 Numerical simulations of fluid-structure interaction (FSI) on an elastic foil heaving with constant amplitude in freestream flow are carried out at a low Reynolds number of 20,000. The commercial software STAR-CCM+ is employed to solve the flow field and the large-scale passive deformation of the structure. The results show that introducing a certain degree of flexibility significantly improves the thrust and efficiency of the foil. For each Strouhal number St considered, an optimal flexibility exists for thrust; however, the propulsive efficiency keeps increasing with the increase in flexibility. The visualisation of the vorticity fields elucidates the improvement of the propulsive characteristics by flexibility. Furthermore, the mechanism of thrust generation is discussed by comparing the time-varying thrust coefficient and vortex structure in the wake for both rigid and elastic foils. Finally, in addition to sinusoidal motions, we also consider the effect of non-sinusoidal trajectories defined by flattening parameter S on the propulsive characteristics for both rigid and elastic foils. The non-sinusoidal trajectories defined by S=2 are associated with the maximum thrust, and the highest values of propulsive efficiency are obtained with S=0.5 among the cases considered in this work.


INTRODUCTION
The fluid-structure interaction (FSI) phenomenon exists widely in nature and engineering systems [1][2][3] such as aerospace engineering, ocean engineering and biomedical engineering. These engineering fields always involve complex geometric shapes and large-scale active or passive flow-induced deformations. The nonlinear deformation of a structure is still challenging in FSI problems. Birds [4,5], insects [6,7] and aquatic animals [8] can gain efficient propulsion by flapping their flexible or elastic wings and fins. To simplify the model, in previous experiments or numerical simulations to analyse thrust generation and how to achieve efficient propulsion, the wings/fins have usually been considered as a rigid plate or foil [7,9,10]. However, in fact, birds and aquatic animals will interact with the surrounding viscous fluid to generate large-scale deformations when they use their wings/fins to fly/ swim. Inspired by this, the FSI problem in the propulsion of underwater bionic vehicles has been a hot topic in recent years [11]. Therefore, it is necessary to study the effect of flexibility on the thrust generation and effective propulsion of the foil.
In fact, many studies have reported the propulsion characteristics of rigid plates/foils immersed in free flow in the past few years. These plates/foils are set to undergo pitching, heaving and combined pitching and heaving motions in the fluid. Within a specific range of pitching and heaving frequency, the vortices on the surface of the plate/foil are shed and form a reverse von Kármán vortex street which is * Corresponding author: yupengyao@dlmu.edu.cn (P. Yu) a signature of thrust generation on the plate/foil [12][13][14][15]. In most of the literature, the dimensionless Strouhal number St is used to characterise the frequency of vortices shedding. Triantafyllou et al. [16] conducted experimental measurements on oscillating foils immersed in a uniform oncoming flow and found that the optimal propulsive efficiency was always attained within a narrow range of St (0.2≤St≤0.4). Eloy [17] reported similar conclusions to Triantafyllou et al. [16]. Pedro et al. [14] numerically investigated the influence of the Strouhal number and pitch angle on the efficiency of a foil undergoing the combined pitching and heaving motions. Similarly, Lewin and Haj-Hariri [18] numerically studied the power coefficient and fluid characteristics of a heaving foil at low Reynolds number. The effects of asymmetric non-sinusoidal trajectories were also discussed in this paper. The results showed that for a fixed St, the propulsive efficiency attained by the intermediate heaving frequency was always the largest, which was contrary to the conclusion of the ideal model of Wang et al. [19]. In addition to thrust force and efficiency, the wake patterns behind the rigid foil were also necessary for propulsion research [20]. Schnipper et al. [21] conducted an experimental study on a symmetrical pitching foil immersed in a vertically flowing soap film to reveal the transition mechanism of the wake by varying the oscillation frequency and amplitude of the flapping foil. Anderson et al. [22] compared the wake characteristics of pitching and heaving foils by a combination of experimental and numerical methods proposed in their study. The results showed that pitching and heaving foils had qualitatively similar wake maps spanned by the amplitude and flapping frequency. And the drag-thrust transition had a strong correlation with the changes in wake structure at low frequency and high amplitude.
Previously, we observed that the wings/fins of birds/ fish are simply modelled as a rigid foil to discuss bionic propulsion. However, these wings/fins are flexible and it is necessary to discuss the effect of that flexibility on the propulsive characteristics and vortex structure of the foil [23][24][25]. Heathcote and Gursul [26] experimentally measured the thrust and efficiency of a heaving flexible teardrop/flat plate at Reynolds numbers of 9,000 to 27,000. The study reported that the maximum thrust and efficiency were obtained in the case of intermediate stiffness for all Reynolds numbers considered. For a fixed Reynolds number and Strouhal number, the vortex structures corresponding to the maximum thrust coefficient and maximum efficiency were analysed and compared. The stronger vortex structures behind the foil were associated with higher thrust, and weaker leading edge vortices were associated with higher efficiency according to the visualisation of the flow fluid. Soon after, Alben [27] analysed the thrust generation and vortex patterns of a pitching plate with a new formulation of motion at a small amplitude. Michelin and Smith [28] used the potential flow theory to investigate the effect of flexibility on the propulsive performance of a twodimensional heaving wing in the inviscid limit. Recently, Zhang et al. [29] studied the effect of the mass ratio on the thrust generation of a two-dimensional elastic panel immersed in free flow. At low or intermediate mass ratios, the system resonance was beneficial for the improvement of the thrust force and propulsive efficiency. However, for the cases of high mass ratio, the inertia of the panel was dominant. Thus, the system resonance was not conducive to thrust generation and significantly reduced the efficiency. Regarding the effect of spanwise flexibility, Heathcote et al. [30] performed a water tunnel experiment on the effect of spanwise flexibility on the propulsion characteristics of a rectangular heaving wing. The results showed that a certain degree of spanwise flexibility can improve the thrust and efficiency. The Strouhal number for which flexibility was beneficial for propulsive performance was consistent with that observed in aquatic animals and birds in nature, namely 0.2≤St≤0.4.
Furthermore, most of the plates/foils mentioned in the previous literature experienced simple harmonic motions.
Recently, it has also been observed that non-sinusoidal trajectories have a significant effect on the propulsion characteristics. The experiment by Read et al. [31] made a sinusoidal variation of the angle of attack by introducing higher harmonics in the heaving motion. The thrust coefficient was greatly improved at high Strouhal numbers. Xiao and Liao [32] numerically investigated the effect of the harmonic cosine angle of attack profile on the propulsive performance. The results showed that the degradation of thrust and efficiency could be greatly alleviated or removed at a high Strouhal number when the angle of attack was imposed as a cosine function. The stronger reverse von Kármán vortex street in the wake was considered to induce the improvement of the propulsive performance. A numerical study by Lu et al. [33] reported that non-sinusoidal trajectories affected the aerodynamics of the two-dimensional foil by affecting the instantaneous force coefficient and flow structure. The author reported that non-sinusoidal motions could also enhance the reverse von Kármán vortex street in the wake, which induced an increase in thrust force. More recently, Boudis et al. [34] proposed a new type of non-sinusoidal function to calculate the aerodynamics of an NACA0012 foil undergoing combined pitching and heaving motions at a Reynolds number of 11,000. The authors reported that non-sinusoidal motions greatly improved the thrust force; however, the propulsive efficiency attained with sinusoidal trajectories was always higher than that with non-sinusoidal motions.
While several previous researchers have reported the thrust generation, flow structure around the plate/foil and the effect of kinematic parameters, the effects of large-scale deformation due to the interaction with the fluid have rarely been reported. Recently, Manjunathan et al. [35] reported the thrust generation mechanism of an elastic flat plate which was set to experience sinusoidal motions. However, the effect of non-sinusoidal trajectories on the thrust generation is not discussed. And the hydrodynamic performance of the NACA airfoil is better than that of flat plates. Therefore, it is necessary to investigate the effect of non-sinusoidal trajectories on the thrust characteristics of the flexible NACA airfoil. The main aim of this present work is to investigate the effect of flexibility on the propulsive characteristics, vortex shedding and wake pattern behind the foil at a low Reynolds number of 20,000. We also explore the effect of non-sinusoidal trajectories on the propulsive performance of an elastic foil. The rest of the paper is arranged as follows: First, we describe the problem of bionic propulsion and give the definition of related parameters in section 2. And we present the results and discuss the influence of flexibility on the propulsive performance in Section 3. Finally, we summarise the paper in Section 4.

COMPUTATIONAL APPROACH PROBLEM DESCRIPTION
The fluid-structure interaction (FSI) of a two-dimensional elastic NACA0012 foil undergoing pure heaving motion is considered at low Reynolds numbers. The chord length c and dimensionless relative thickness are 0.1 m and 0.12 c, respectively. The heaving motion of the foil is governed by the following equation: where h 0 denotes the dimensionless amplitude of the heaving motion, ω = 2πf denotes the angular frequency. f is the oscillation frequency of the flapping foil. It is observed that the leading edge of bird wings and fish fins has greater rigidity and produces hardly any flexible deformation [36]. Therefore, a new bionic deformation shape of the flexible foil is introduced in this paper. Only the trailing edge produces flow-induced deformation, as shown in Fig. 1. The distance between the leading edge and the deformation position of the trailing edge X S is 1/3c. The pitching motion is introduced by the flexibility of the foil, unlike most rigid foils undergoing combined heaving and pitching motion in the previous literature. The dimensionless bending stiffness parameter is expressed as the following equation: In the present work, the relative dimensionless bending parameter λ* is defined as the ratio of each foil to the most flexible foil. The details are shown in Table 1. The Reynolds numbers Re is expressed as follows: where ρ 1 denotes the fluid density. μ denotes the dynamic viscosity of the fluid. U∞ denotes the freestream velocity. The Reynolds number is fixed at 2×10 4 in this paper. In the works of Boudis et al. [34], the flattening parameter S was used to define the non-sinusoidal trajectories of the foil as follows: where S denotes the flattening parameter. When S = 1, the trajectory is sinusoidal. Eq. (4) is used to achieve the nonsinusoidal trajectories of the flexible foil. Fig. 2 shows the non-sinusoidal trajectories; when S = 1, the trajectory of the oscillating foil is harmonic. The Strouhal number is an important dimensionless parameter to characterise the kinematics in the study of the hydrofoil propulsion, defined as follows: To evaluate the propulsive performance of the elastic foil, we define the following quantities. The time-varying thrust coefficient and the time-varying power-input requirement coefficient are calculated by the following equations: where F x is the instantaneous force of the foil in the thrust direction. The thrust coefficient C T consists of two components, namely, the pressure component and shear component, expressed as C T,P and C T,S . F y is the instantaneous force of the foil in the lift direction. v(t) is the velocity of the leading edge of the foil in the lift direction. The time-averaged coefficients of thrust and power are calculated from the average force per unit period and they can be expressed as follows: Moreover, the propulsive efficiency can be defined as follows:

GOVERNING EQUATIONS
The flow around the two-dimensional elastic foil is assumed to be unsteady, viscous, and incompressible. It can be solved by calculating the Navier-Stokes equations as follows: where p is the pressure, F is the volume forces, u 1 is the fluid velocity, and ν is the viscosity.
The commercial computational fluid dynamic solver STAR-CCM+ is employed to solve the flow field. The Finite Volume Method is used to discretize the N-S equations. The SIMPLE (Semi-Implicit Method for Pressure-Linked Equations) algorithm is used to solve pressure-velocity coupling of the flow field. The pressure and momentum are discretised by a second-order implicit scheme. The turbulent model SST k-ω is used in all numerical simulations. The deformation of the flexible foil is governed by the structural momentum conservation equation, which can be expressed follows: where σ is the stress tensor, ρ 2 is the foil density, b is the body force, d is the displacement of the structure. The structure is considered as a linear constitutive material and can be expressed as follows: where σ is the stress tensor, ε is the strain tensor, I is the identity tensor, λ is Lamé's first parameter, μ is the shear modulus. The Finite Element Method is used to solve the large-scale deformation of the foil. Implicit coupling is employed. The fluid domain and structure domain exchange the geometric and force information at each time step. The pressure and shear force calculated are enforced as boundary conditions of the structure solver. The displacements of the structure are enforced as boundary conditions of the fluid solver. The interface of the fluid and structure needs to ensure that the velocity and displacement are continuous functions. The equations governing these boundary conditions are expressed as follows: where subscript f denotes the fluid, subscript s denotes the structure. Fig. 3 shows the computational domain and boundary conditions. The computational domain consists of two parts, a background zone and an overset zone. The background zone is a rectangle with 60c × 40c, and the overset zone is a circle with a radius of 1c. A no-slip condition is imposed on the elastic foil surface, which means that the velocity of the fluid over the foil is zero. At the left, top and bottom boundaries, the boundary condition is set to the velocity inlet, which means that the pressure is zero gradient and the velocity is equal to freestream velocity U∞ = 0.2 m/s. The right boundary is set to the pressure outlet. The mesh around the foil is refined to accurately simulate the fluid structure. To make the y + ≤1, the height of the first layer grid is 10 -5 c. The overall situation of the grid is shown in Fig. 4.

CODE VALIDATION
To confirm the accuracy of the FSI solver in STAR-CCM+, we carried out two validations and the results were compared with the experimental data and numerical simulation data reported in the literature. The first validation is a heaving NACA0012 foil and the time-averaged thrust coefficient is calculated using the kinematic parameters h 0 = 0.175, Re = 2×10 4 , and the dimensionless reduced frequency k varying in the range of (2-8), (k = 2πfc/ U∞). In Fig. 6, the results attained in the present work are compared with the experimental data of Heathcote et al. [30] and the numerical data of Young and Lai [37] dB. It can be seen that all the simulation results show similar trends and are in good agreement.
The second validation is the numerical simulation of a flowinduced passive deformation of a flexible plate connected behind the cylinder at low Reynolds numbers. The fluid-structure interaction of a flexible plate added behind the cylinder proposed by Turek and Hron [38] has been a typical problem in the study of FSI in recent decades. The cylinder and the elastic plate are placed in a rectangular laminar field. The fluid flow through the cylinder and the unsteady vortex is shed, which stimulates the flexible plate to periodically oscillate. The computational configuration is shown in Fig. 7. The dimensionless length and thickness of the plate are 3.5D and 0.2D, respectively. The structure to fluid density ratio, dimensionless Young's Modulus (E/ρf U∞ 2 ), Reynolds number and Poisson ratio are set to 10, 1400, 100, and 0.4, respectively. The no-slip boundary condition is imposed at the top and bottom of the computational domain. The left is set to the velocity inlet, and a parabolic velocity profile is set as follows: v 1 (0, y) = 1.5U (16) where U is the average inflow velocity and the maximum inflow velocity is 1.5U. The left of the plate is fixed on the cylinder surface. The no-slip boundary condition is applied on the surface of the plate. To realise accurate simulations of the flow structure, the grid around the cylinder and plate is refined. The minimum grid size is set as 0.03D to ensure the grid converges. Fig. 8 shows the periodic vibration of the splitter plate with multiple bending modes caused by alternate shedding of vortices over the cylinder. The time-varying displacement of the reference point A (Fig. 7) in the y-direction when the plate achieves periodic self-sustaining oscillation is shown in Fig. 9. The result shows great agreement with the data reported by Turek and Hron [38]. The validations performed in this section confirmed that the FSI solver in STAR-CCM+ is able to simulate the deformation of the elastic hydrofoil.

EFFECT OF FLEXIBILITY ON THRUST AND EFFICIENCY
We study the effect of flexural rigidity on the propulsive performance of a flexible foil undergoing sinusoidal motion in this section. The Reynolds number Re, amplitude of the heaving motion h 0 , the structure to fluid density ratio and Poisson ratio are set to 2 × 10 4 , 0.3, 7.8 and 0.3, respectively. The Strouhal number St is varied in a narrow range from 0.2 to 0.4, consistent with the St of animals in nature (such as insects, birds, and fish) observed in Taylor et al. [39]. Fig. 10 shows the overall propulsive performance as a function of the dimensionless relative bending stiffness, including the tip amplitude of the trailing edge (h Tip ), the time-averaged thrust coefficient (C T,mean ), the time-averaged power coefficient (C P,mean ) and the propulsive efficiency (η). As expected, the trends of hTip with λ* are very similar within the range of St considered in this work and it decreases rapidly with the increase of bending stiffness of the foil at λ*≤10. Then, it decreases steadily at λ*≥10. h Tip for higher St is larger than that for lower St over the whole range of bending stiffness. This is attributed to the increase in the reaction force exerted by the fluid on the foil, as the oscillation frequency of the foil increases. The variation trends of C T,mean and C P,mean with λ* are plotted in Fig. 10(b) and Fig. 10(c). The trends are similar and show a peak in C T,mean and C P,mean for all St considered, which means that an optimal bending stiffness exists for thrust. The peak moves to a higher value of bending stiffness with the increase of the Strouhal number St. For St = 0.2, 0.25 and 0.3, C T,mean and C P,mean gradually increase until λ* = 2 and start decreasing beyond it. C T,mean and C P,mean increase sharply until λ* = 6 and start falling beyond it for St = 0.35 and 0.4. Fig. 11 shows the comparison of the averaged velocity profile in the wake of the rigid and elastic foils for St = 0.4. The three foils correspond to the rigid foil, the one that produces the largest thrust, and the most flexible foil. For the elastic foil of λ* = 6, a higher positive additional x component of velocity is observed than in the most flexible foil (λ* = 1). The magnitude of jet-like velocity is considered to be related to the generation of thrust, so a suitable flexibility can produce larger thrust than the most flexible foil (λ* = 1). The far greater flexibility is found to be disadvantageous in terms of thrust generation. The Strouhal number St has a relatively large effect on the thrust improvement especially for high St. For St=0.4, the maximum C T,mean obtained by the flexible foil is around 1.42, almost 1.6 times larger than that by the rigid foil. However, it is 1.28 times that for St = 0.2. C T,mean obtained with high St is always higher than that obtained with low St and shows a monotonical increase with the increase of St. The maximum C T,mean obtained with St = 0.4 (at λ* = 6) is around 6.6 times larger than that with St = 0.2 (at λ* = 2). It can be predicted that a higher St may yield even higher thrust. Fig. 10(d) shows the propulsive efficiency η as a function of bending stiffness λ*. For St = 0.3, 0.35 and 0.4, η reaches the peak at λ* = 2, then steadily decreases at λ*≥2. For St = 0.2 and 0.25, there is no visible optimal value for η, and it decreases in the overall range of bending stiffness. It is assumed that introducing greater flexibility of the foil may produce higher efficiency. Unlike the thrust coefficient, a lower St seems to be more conducive to efficiency improvement.
The wake signature behind the foil has a significant effect on the propulsive performance. Fig. 12 shows the vorticity fields for three foils with different bending stiffness. The leading edge of the foil is moving downwards through the mean position at t = 1/4T. The rotational direction and position of the vortex in the wake have a decisive effect on the thrust generation of the foil. As shown in Fig. 12, the vortex pairs are arranged in two rows, with the counterclockwise vorticity on the top and the clockwise vorticity on the bottom. The arrow indicates the direction of the flow velocity. Induced by the vortex in the wake, the flow direction is the same as the direction of freestream velocity. Therefore, the direction of the reaction force exerted by the fluid on the foil is opposite to that of freestream. The formation of a reverse von Kármán vortex street in the wake generates thrust on the foil. It is noted that the vortex pairs from the most flexible foil (λ* = 1) are weaker and more compactly distributed in the x-direction than the other two foils. This vortex structure leads to a severe reduction in the thrust force of the most flexible foil. It is also seen that the deformation of the trailing edge of the elastic foil reduces the flow separation of the leading-edge vortex and makes it weaker. Minimising the strength of the leading-edge vortex is beneficial for achieving high efficiency [40,41]. The flexibility of the foil can improve both the thrust and the efficiency, whereas the bending stiffness for maximising thrust and efficiency is not same. The foil maximising the efficiency seems to be more flexible than that maximising the thrust [ Fig. 10(b) and Fig. 10(d)]. The vortex structure for the two foils oscillating in one typical period for St = 0.4 is shown in Fig. 13. Fig. 13(a)-(d) correspond to the stiffness (λ* = 6) that produces the maximum thrust and Fig. 13(e)-(h) correspond to the stiffness (λ* = 2) that produces the maximum efficiency. The vortices form on the leading edge of the foil and are swept into the trailing edge vortex system in both cases. The vortices in the wake from the stiffer foil (λ* = 6) are seen to be stronger, which induces high thrust. And the more flexible foil is better at making the leading edge vortices weaker, which is consistent with the higher efficiency.

ANALYSIS OF MECHANISM OF THRUST GENERATION
The mechanism of thrust generation is analysed for a rigid and an elastic foil by comparing the temporal variation of C T , C P and the wake structure in this section. The amplitude of heaving motion is h 0 = 0.3 for both cases. The Strouhal number and Reynolds number based on the chord length are 0.3 and 2×10 4 , respectively. The dimensionless relative bending stiffness is fixed at λ* = 6 corresponding to the maximum thrust coefficient and efficiency for St = 0.3 (Fig. 10). Fig. 14(a) shows the time-varying displacement of the tip of the trailing edge h T and instantaneous thrust coefficient C T with two components of pressure and shear force for a rigid foil. The time-varying power coefficient for a rigid foil is plotted in Fig. 14(b). It is noted that the variation of h T is sinusoidal and three peaks are observed in three typical cycles as expected, whereas six peaks are shown in the temporal variation of C T . In other words, C T attains its peaks twice in one oscillation cycle. C T is minimum and negative at t/T = 1.0 when the foil is at the top position. Then, C T increases sharply to the maximum value when the foil moves downwards to reach the position between the top and mean position at t/T = 1.2. At t/T = 1.5, the foil moves to the bottom position and C T decreases from the maximum to the minimum value. When the foil returns Another pair of vortices also develop on the lower surface, which will shed from the trailing edge in the next cycle of downward movement. It is also observed that although the vortices form and are shed in pairs, the strength of the pair of vortices is not equal. The clockwise vortex will eventually vanish in the wake when the foil moves downwards, while the counterclockwise vortex will vanish when the foil moves upwards.

Elastic foil
According to the analysis of the rigid foil, the time-varying displacement of the tip of the trailing edge h T and time-varying thrust coefficient C T with two components of pressure and shear force for an elastic foil in three typical cycles are shown in Fig. 16(a). The temporal variation of displacement of trailing edge h T is sinusoidal, similar to the rigid foil. Due to the passive deformation of the elastic foil, the time at which h T reaches its maximum value is slightly later than that of the rigid foil. The peak value of h T is 0.38, which is 26.7% larger than that of the rigid foil. The temporal variation of C T is sinusoidal. The variation frequency of C T is twice as large as the oscillation frequency of the foil. It is noted that the deformation of the elastic foil has little effect on the shear force component of thrust; however, it greatly improves the positive pressure component. The time-averaged thrust coefficient C T,mean is significantly improved and is 0.56, 55.6% larger than the rigid foil. The time-averaged power coefficient C P,mean is 2.53, 7.2% larger than the rigid foil. The flexibility of the foil greatly increases C T,mean while having little effect on C P,mean . Based on this point, the propulsive efficiency η is notably improved. η for the elastic foil at St = 0.3 is around 0.22, 47.6% larger than the rigid foil.
In Fig. 17, the vortex shedding pattern for the elastic foil is similar to the rigid foil described previously in detail. The vortices form and develop on the leading edge of the foil and are then shed from the trailing edge. The reverse von Kármán vortex street is observed in the trailing edge vortex system behind the foil. The shedding of the vortices minimises the thrust when the foil moves to the top or bottom position.

EFFECT OF NON-SINUSOIDAL TRAJECTORIES ON PROPULSIVE PERFORMANCE
In this section, we consider the effect of flexibility on the propulsive performance of an elastic foil undergoing two different non-sinusoidal trajectories defined by the flattening parameters S. The mathematical equation of the non-sinusoidal trajectories is based on the work of Boudis et al. [34]. The settings of the Reynolds number Re and amplitude of the heaving motion h 0 are the same as the previous sinusoidal motion, which are 2×10 4 and 0.3c, respectively. The Strouhal number St is varied in a narrow range from 0.2 to 0.4. Fig. 18(a)-(b) shows the time-averaged thrust coefficient C T,mean in one typical cycle as a function of bending stiffness λ* for S = 2 and S = 0.25. The trends of C T,mean for the two nonsinusoidal trajectories are similar to sinusoidal motion. A peak value in the C T,mean curves for all St considered in this section is also observed. The peak values move to a larger λ* with the increase of St in two cases. For St =0.2 and 0.25, CT,mean increases until λ* = 6 and steadily decreases beyond λ* = 6. For St = 0.3 and 0.35, C T,mean increases rapidly at λ*≤10 and decreases beyond it. C T,mean increases sharply until λ* = 15 and starts decreasing at λ*≥15. To clearly evaluate the effect of nonsinusoidal trajectories on C T,mean , the maximum time-averaged thrust coefficient MAX(C T,mean ) for each St in three motions is plotted in Fig. 18(c). It is observed that MAX(C T,mean ) increases monotonically with the increase of St. The MAX(C T,mean ) obtained with non-sinusoidal trajectories is much larger than that with sinusoidal trajectories, especially for high St. At St = 0.4, the values of MAX(C T,mean ) for S = 2 and 0.5 are 2.27 and 2.06, respectively, which are around 1.6 and 1.45 times larger than that for S = 1.
The time-averaged power coefficient C P,mean for S = 2 and S = 0.5 as a function of bending stiffness is plotted in Fig. 19(a)-(b). The variation trends of C P,mean are similar to C T,mean and show a peak value in the curves for two non-sinusoidal motions considered in this section. The peak value moves to high λ* with the increase of St. For S = 2, for each St, the variation trends of C P,mean with λ* are extremely steady at high bending stiffness and C P,mean increases sharply with the increase of λ* at low stiffness. By contrast, for S = 0.5, C P,mean reduced rapidly at high λ*. Fig. 19(c) shows the comparison of the maximum time-averaged power coefficient MAX(C P,mean ) attained by nonsinusoidal trajectories and sinusoidal trajectories. For each St, especially for high St, MAX(C P,mean ) for S = 2 is far greater than that for S = 1 and S = 0.5. For St = 0.4, the values for MAX(C P,mean ) for S = 1 and S = 0.5 are 6.35 and 8.33, respectively. For S = 2, the value is 14.8, which is 2.33 and 1.78 times larger than that for S = 1 and S = 0.5. An excessive power-input requirement seems to indicate low propulsive efficiency.  The maximum efficiency obtained with S = 0.5 is much higher than that with S = 1 and S = 2. The maximum efficiency obtained with S = 0.5 at St=0.2 is 0.4, 1.53 times larger than that with S = 1. This is due to the substantial increase in thrust and the limited power input. For S = 2, the efficiency is lower than that with S = 1, hence the great power input.
To analyse the improvement of non-sinusoidal motions on the thrust and efficiency, we compare the vortex structure around the elastic foil in one typical cycle at St = 0.3 as shown in Fig. 21 and Fig. 22. The bending stiffness in the two cases is λ* = 15, corresponding to the maximum thrust coefficient for non-sinusoidal trajectories. For S = 2, it is observed that two pairs of vortices in opposite rotating directions are shed from the trailing edge of the elastic foil in one cycle. While the foil is moving downwards from t/T = 0 to t/T = 0.5, two counterclockwise vortices are sequentially shed, while during the reverse movement from t/T = 0.5 to t/T = 1.0, two clockwise vortices are shed. This is different from the way of shedding vortices of sinusoidal trajectories, in which only one vortex is shed per half-cycle. It is also noted that the first shedding vortex forms and develops on the leading edge, whereas the second shedding vortex is on the trailing edge. This will further increase the pressure differential between the upper and lower surfaces of the foil, thereby increasing the lift force compared to that of sinusoidal trajectories. The mode of two pairs of vortices shedding is associated with the case of the maximum thrust force. As shown in Fig. 22, for S = 0.5, the mode of shedding vortices is similar to S = 2. However, the second vortex is much weaker and vanishes rapidly without completely shedding. It is observed that the vortices have not fully developed on the surface of the foil. In the previous section, we pointed out that the leading edge vortex has a significant effect on the propulsive efficiency. This undeveloped leading edge vortex seems to greatly improve the efficiency of the foil for S = 0.5, which is consistent with the conclusion in Fig. 20(c).

CONCLUSIONS
The effect of flexibility on the propulsive performance of a heaving foil undergoing sinusoidal and non-sinusoidal motions at a low Reynolds number has been numerically investigated in this study. We employ the commercial software STAR-CCM+ to solve the flow field and large-scale passive deformation of the foil. Two validations have been performed and the accuracy of the solver in STAR-CCM+ has been verified by comparing results with the published data in the literature. We analyse the overall propulsive performance as a function of the dimensionless relative bending stiffness. Introducing a certain degree of flexibility can significantly improve the thrust and efficiency of the foil. A peak value is observed in the thrust coefficient, which means that an optimal bending stiffness exists for thrust. The velocity profile in the wake is used to illustrate the improvement. A higher positive additional x component of velocity is observed for the elastic foil. The higher magnitude of jet-like velocity will result in a higher thrust. For optimal stiffness, the maximum time-averaged thrust coefficient obtained with the elastic foil is almost 1.6 times larger than that with a rigid foil. For efficiency, there is no visible optimal value within the range of stiffness studied in this work. Introducing a greater degree of flexibility seems to be more beneficial for efficiency. From the analysis of the flow field and vortex structure in the wake, it can be concluded that the higher vortices intensity of the trailing edge is related to the case of high thrust, and the weaker leading edge vortices are associated with the case of high efficiency. A higher Strouhal number St is more beneficial for thrust generation and a lower St seems to be more conducive to efficiency improvement.
The thrust generation mechanism of rigid and elastic foils is analysed. The variations of h T and C T are similar for the rigid and elastic foils. C T attains its peaks twice in one oscillation cycle in both cases. The maximum C T of the elastic foil is much higher than that of the rigid foil. By analysing the vortex structure for both cases, it is noted that the vortices form and develop on the leading edge and then move backwards along the surface of the foil. Finally, a pair of counter-rotation vortices are shed from the trailing edge of the foil to form the reverse von Kármán vortex street, which induces the generation of thrust.
Furthermore, the effect of non-sinusoidal trajectories on the propulsive performance of elastic foils is studied. The nonsinusoidal trajectories are realised by the flattening parameters S. The results show that the non-sinusoidal trajectories significantly affect the propulsive performance and vortex structure of the foil. For non-sinusoidal motions of S = 2, the time-averaged thrust coefficient is much higher than that obtained with sinusoidal motions. However, since the power input has also been greatly increased, the efficiency is lower than that obtained with sinusoidal motions. For non-sinusoidal motions of S = 0.5, the time-averaged thrust coefficient is also greatly improved and the power input is only slightly increased compared to the sinusoidal motions. The optimal propulsive efficiency of the elastic foil is obtained at S = 0.5. The impact of non-sinusoidal trajectories on thrust generation and efficiency is discussed in more detail by visualisation of the flow field. For S = 2, the shedding vortices show the new mode. Two pairs of counter-rotating vortices are shed in one cycle, which induces a stronger jet to improve the thrust. For S = 0.5, the surface of the foil shows the undeveloped vortices, which are associated with the high efficiency.