The Effect of Spanwise Folding on the Aerodynamic Performance of a Passively Deformed Flapping Wing

The wings of birds exhibit multi-degree-of-freedom motions during flight. Among them, the flapping folding motion and chordwise passive deformation of the wings are prominent features of large birds in flight, contributing to their exceptional flight capabilities. This article presents a method for the fast and accurate calculation of folding passive torsional flapping wings in the early design stage. The method utilizes the unsteady three-dimensional panel method to solve the aerodynamic force and the linear beam element model to analyze the fluid–structure coupling problem. Performance comparisons of folding flapping wings with different kinematics are conducted, and the effects of various kinematic parameters on folding flapping wings are analyzed. The results indicate that kinematic parameters significantly influence the lift coefficient, thrust coefficient, and propulsion efficiency. Selecting the appropriate kinematic and geometric parameters is crucial for enhancing the efficiency of the folding flapping wing.


Introduction
Unlike fixed-wing aircraft, flapping wings have numerous degrees of freedom, encompassing various motions and deformations [1], such as active flapping and twisting, passive torsional and bending deformations, and spanwise folding.These dynamic features enable flapping wings to efficiently generate lift and thrust, thereby enhancing flight performance [2].Consequently, comprehending the impact of wing motions on flight performance is pivotal in flapping-wing design.Given that the spanwise folding motion is a primary characteristic of large birds in flight [3], this paper predominantly investigates the influence of the folding motion combined with the wing's passive torsional deformation on performance.
Flapping wings with folding motions have been studied.Send [4] designed the famous smart bird and utilized thin-plate theory and CFD software to study the unsteady two-dimensional airfoil.Finally, experiments were applied to the analysis.Kim [5] built a folding bird-like flapping wing, analyzed its performance using simple methods, and conducted flight tests.Huang [6] employed the two-dimensional airfoil method to solve the folding flapping wing's aerodynamic forces and optimized the structural parameters in the linkage mechanism.Han [7] applied the panel method to study the unsteady aerodynamic force of the flapping wing and found that, compared with simple up-and-down flapping, the wing combined with the folding motion can obtain greater lift and thrust, but the wing twist was not considered.Yang [8] employed the unsteady vortex lattice method to calculate the aerodynamic force of the folding ornithopter and employed an optimization algorithm to optimize the shape and motion parameters to maximize the average lift.Karimian [9] applied a simple two-dimensional lift line model and a linear beam element model to study the aeroelastic model of a flapping wing with folding.He found that a wing with folding motion can improve flight performance and propulsion efficiency, and that the phase of the movement of the outer wing has a great influence on aerodynamic performance.Verstraete [10] utilized the unsteady vortex lattice method and a linear finite element model to simulate a seagull-like flapping wing with folding and studied the effect of folding on the flutter speed.Chang [11] applied the method of solving the Navier-Stokes equations to calculate the foldingflapping-wing model.The torsional deformation in the model was preset without the fluid-structure interaction calculation.He found that lift and thrust were mainly generated in the downward-flapping stage.Lang [12] applied commercial software to solve N-S equations to study the effects of the folding parameters of a flat-panel wing on its performance.Bie [13] also applied a commercial N-S equation solver to study a bat-like folding flapping wing and analyzed the effects of the inner-to-outer wing ratio.Ryu [14] experimentally studied the folding flapping wing with a four-bar linkage mechanism and concluded that the folding motion of the wing is conducive to the generation of lift.Chen [15] employed an experimental method to study the influence of wing-folding amplitudes on the vertical aerodynamic force of a small flapping wing.Qin [16] employed a PIV experiment to study the evolution and physical characteristics of the tip vortex of a folding flapping wing.
The existing research exhibits numerous shortcomings.The design and manufacture of folding flapping wings have basically relied on experiments and simple analysis methods, which makes it difficult to guide the design.Regarding analytical approaches, the single-layer vortex lattice method is unable to compute the surface pressure coefficient; it can only assess changes in lift and thrust, leaving the wing surface state unclear.However, the surface pressure coefficient significantly influences the airfoil lift-drag ratio [17].Although the method of solving the N-S equations can accurately determine the aerodynamic force, it is time-consuming and unsuitable for the initial design stage.Considering the fluid-structure coupling problem further exacerbates the computational time, rendering it impractical.Regarding research problems, the impact of folding on propulsion efficiency has not been considered.Propulsion efficiency is, in fact, crucial for achieving efficient flight.Studies seldom consider the effects of wing torsional deformation and the geometric parameters of the wing itself, concentrating solely on the impact of the folding motion.This oversight results in widespread separation on the wing surface, deviating from real-world scenarios [18].
To address the shortcomings of current research, this article proposes a method for solving the fluid-structure coupling problem of folding flapping wings.The unsteady threedimensional panel method is employed to calculate the wing surface pressure coefficient and aerodynamic parameters, while the linear beam element model is utilized to account for the outer wing's passive torsional deformation, as shown in Figure 1.A kinematic analysis of folding flapping wings with varying kinematic parameters is conducted.A dynamic analysis of folding flapping wings with different parameters undergoing passive twisting due to aerodynamic forces is performed.Performance parameters such as the lift coefficient, thrust coefficient, and propulsion efficiency are compared and analyzed.Various aerodynamics throughout the period are compared, and the wing surface pressure coefficients are analyzed.The influence of different parameters on the aerodynamic performance of folding flapping wings is discussed, their mechanism is analyzed, and design guidance is provided in the conclusion.

Aerodynamic Model
This article employs the unsteady panel method to calculate the 3D aerodynamic force.The panel method is an inviscid method that has been studied by many researchers.Magnus [19] and Ashby [20] utilized high-order and low-order surface element methods to calculate 3D aerodynamic forces, respectively.Similarly, Maskew [21] and Johnson [22] utilized the panel method to establish a general program for solving wing aerodynamic forces.Vest [18] utilized the panel method to study the flight performance of birds.The method in this article originates from Hess [23], and Hess's method is extended to unsteady situations.
As shown in Figure 2, the panel model consists of a wing panel and a wake panel.The wing panel consists of a source panel and a dipole panel, and the wake panel consists of only the dipole panel.The wing is divided into M portions in the chord direction and N portions by N-lines in the spanwise direction.To save calculation time, the free-wake model is not used, but the fixed-wake model described by Reichert [24] is employed.The velocity (Neumann) boundary condition is employed as the boundary condition, which means that the velocity normal to the wing surface is zero.For unsteady conditions, the trailing edge of the wing must satisfy the unsteady Kutta condition [25], which is

Aerodynamic Model
This article employs the unsteady panel method to calculate the 3D aerodynamic force.The panel method is an inviscid method that has been studied by many researchers.Magnus [19] and Ashby [20] utilized high-order and low-order surface element methods to calculate 3D aerodynamic forces, respectively.Similarly, Maskew [21] and Johnson [22] utilized the panel method to establish a general program for solving wing aerodynamic forces.Vest [18] utilized the panel method to study the flight performance of birds.The method in this article originates from Hess [23], and Hess's method is extended to unsteady situations.
As shown in Figure 2, the panel model consists of a wing panel and a wake panel.The wing panel consists of a source panel and a dipole panel, and the wake panel consists of only the dipole panel.The wing is divided into M portions in the chord direction and N portions by N-lines in the spanwise direction.To save calculation time, the free-wake model is not used, but the fixed-wake model described by Reichert [24] is employed.

Aerodynamic Model
This article employs the unsteady panel method to calculate the 3D aerodynamic force.The panel method is an inviscid method that has been studied by many researchers Magnus [19] and Ashby [20] utilized high-order and low-order surface element methods to calculate 3D aerodynamic forces, respectively.Similarly, Maskew [21] and Johnson [22 utilized the panel method to establish a general program for solving wing aerodynamic forces.Vest [18] utilized the panel method to study the flight performance of birds.The method in this article originates from Hess [23], and Hess's method is extended to un steady situations.
As shown in Figure 2, the panel model consists of a wing panel and a wake panel The wing panel consists of a source panel and a dipole panel, and the wake panel consists of only the dipole panel.The wing is divided into M portions in the chord direction and N portions by N-lines in the spanwise direction.To save calculation time, the free-wake model is not used, but the fixed-wake model described by Reichert [24] is employed.The velocity (Neumann) boundary condition is employed as the boundary condition which means that the velocity normal to the wing surface is zero.For unsteady conditions the trailing edge of the wing must satisfy the unsteady Kutta condition [25], which is The velocity (Neumann) boundary condition is employed as the boundary condition, which means that the velocity normal to the wing surface is zero.For unsteady conditions, the trailing edge of the wing must satisfy the unsteady Kutta condition [25], which is where V t is the tangential velocity, and Φ is the velocity potential of the trailing edge at time step k.More details can be found in Reference [26].

Structural Model
To simulate the torsional deformation of the wing beam, a linear beam element model is utilized.The local element stiffness matrix is [K]: where EI is the bending stiffness and GJ is the torsional stiffness of an element, and l is the length of an element.

Aeroelastic Coupling
The unsteady aerodynamic force, the wing mass's inertial force, and the beam's deformation internal force all achieve a dynamic balance during the flapping of the wings.There is an interaction between the deformation of the wing and aerodynamic forces.Given the data at time step k, the equilibrium equation at time step k + 1 can be solved: where [K] is the stiffness matrix, [M] is the mass matrix, → q k+1 is the deformation displacement, and → F k+1 is the force and moment on the beam.The Newmark-β method [27] is utilized to solve this second-order equation.This method is unconditionally stable.

Performance Parameters
The unsteady Bernoulli theorem [28] states that the pressure coefficient in unsteady flows can be written as follows: where V is the velocity of the wing, V ∞ is the far-field velocity, Φ is the velocity potential, and V move is the motion velocity of the wing with respect to the global coordinate system.The outer wing performs a pitching motion around the outer wing beam, accompanied by a flapping motion around the axis of the tip of the inner wing, while the inner wing flaps around the axis of the inner wing root.The torsional moment M t (t), bending moment M b (t), vertical lift F z (t), and horizontal thrust F x (t) are all produced by the moving wing.The wing's output of power is as follows: As shown in Figure 1, the wing consists of two parts: the inner wing and the outer wing.β in (t) is the flapping angle of the inner wing, and β out (t) is the flapping angle of the outer wing.The bending moment at the wing root is M b_in (t), while the bending moment at the outer wing root is M b_out (t).
F x , F z , and P are the following time-averaged quantities: The dimensionless forms of these quantities are Propulsion efficiency is defined as The unsteady aerodynamic force obtained with the present method was validated with the numeral calculation described by Lin [29] on a flapping rectangular wing with the airfoil NACA0014.It is solved by using the 3-dimensional Euler equations.The wing's maximum twisting angle is 4 degrees at the wing tip, and active twisting varies linearly along the span with a phase lag of 90 degrees.As shown in Figure 3, the method employed in this article yields results that closely match Lin's calculations.
The dimensionless forms of these quantities are Propulsion efficiency is defined as 2.5.Solver Validation

Aerodynamic Force
The unsteady aerodynamic force obtained with the present method was validated with the numeral calculation described by Lin [29] on a flapping rectangular wing with the airfoil NACA0014.It is solved by using the 3-dimensional Euler equations.The wing's maximum twisting angle is 4 degrees at the wing tip, and active twisting varies linearly along the span with a phase lag of 90 degrees.As shown in Figure 3, the method employed in this article yields results that closely match Lin's calculations.

Fluid-Structure Coupling
The aeroelastic response obtained with the present method was validated with the flapping wing studied by DeLaurier [30].DeLaurier performed computations and experiments with the flapping wing he designed.Similar to the flapping wing in this article, DeLaurier's flapping wings are also passively twisted.All parameters are provided in Reference [30] to facilitate the comparison of methods.
As shown in Figure 4a,b, compared with the calculations used by DeLaurier, the lift coefficient and thrust coefficient calculated using this method are in good agreement with the experiment.Since the non-stick model is used, there are still some differences, especially when trailing-edge separation occurs on the outer wing after the period becomes small.
The aeroelastic response obtained with the present method was validated with the flapping wing studied by DeLaurier [30].DeLaurier performed computations and experiments with the flapping wing he designed.Similar to the flapping wing in this article, DeLaurier's flapping wings are also passively twisted.All parameters are provided in Reference [30] to facilitate the comparison of methods.
As shown in Figure 4a,b, compared with the calculations used by DeLaurier, the lift coefficient and thrust coefficient calculated using this method are in good agreement with the experiment.Since the non-stick model is used, there are still some differences, especially when trailing-edge separation occurs on the outer wing after the period becomes small.
As shown in Figure 4c, the maximum twist angle amplitudes Θ  along the inner wing are almost the same, but there are some differences in the outer wing sections.The phase angles Φ  between the twisting angle and flapping motion are slightly larger than the results of DeLaurier.As shown in Figure 4d, compared to DeLaurier's results, the flapping elastic amplitudes H  are larger in the inner wing but smaller in the outer wing, and their phase angles Φ  are marginally smaller.These deviations could result from a mix of variations in the structural and aerodynamic models.As shown in Figure 4c, the maximum twist angle amplitudes Θ max along the inner wing are almost the same, but there are some differences in the outer wing sections.The phase angles Φ θ between the twisting angle and flapping motion are slightly larger than the results of DeLaurier.As shown in Figure 4d, compared to DeLaurier's results, the flapping elastic amplitudes H max are larger in the inner wing but smaller in the outer wing, and their phase angles Φ H are marginally smaller.These deviations could result from a mix of variations in the structural and aerodynamic models.

Wing Calculation Model
The wings of the folding flapping wing are divided into two sections, namely, the inner wing and the outer wing.The inner section is torsionally rigid and does not twist with flapping, while the outer wing passively twists with flapping.At the same time, the outer wing section performs a flapping motion relative to the inner section.
In this article, inertia force and aerodynamic force are responsible for the wing's passive torsion.The lift center on the airfoil is around 25% of the chord length to utilize the aerodynamic force to produce torque on the beam as much as possible, and considering the issues of the wing structure, the beam is placed at 10% of the chord line.It is assumed that the beam is the only source of torsional stiffness.As shown in Figure 5a, the beam elements are distributed into N portions along the span direction, which is consistent with the aerodynamic panel.It is assumed that the wing's mass is distributed by volume and is discretely separated into mass locations at 10% of the chord behind the beam on the chord line.As shown in Figure 5b, the forces and moments in the 2D wing profile, which conducts pitching and plunging motions along the elastic axis, are depicted.The effective angle of attack equals the angle of incidence plus the self-motion.

Wing Calculation Model
The wings of the folding flapping wing are divided into two sections, namely, the inner wing and the outer wing.The inner section is torsionally rigid and does not twist with flapping, while the outer wing passively twists with flapping.At the same time, the outer wing section performs a flapping motion relative to the inner section.
In this article, inertia force and aerodynamic force are responsible for the wing's passive torsion.The lift center on the airfoil is around 25% of the chord length to utilize the aerodynamic force to produce torque on the beam as much as possible, and considering the issues of the wing structure, the beam is placed at 10% of the chord line.It is assumed that the beam is the only source of torsional stiffness.As shown in Figure 5a, the beam elements are distributed into N portions along the span direction, which is consistent with the aerodynamic panel.It is assumed that the wing's mass is distributed by volume and is discretely separated into mass locations at 10% of the chord behind the beam on the chord line.As shown in Figure 5b, the forces and moments in the 2D wing profile, which conducts pitching and plunging motions along the elastic axis, are depicted.The effective angle of attack equals the angle of incidence plus the self-motion.

Structural and Kinematic Parameters
The shape of the wing is inspired by the Canadian goose, which has a half wingspan of 0.8 m, an aspect ratio of 12, a root chord length of 0.3 m, and a root-to-tip ratio of 2. The flight velocity is 15 m/s, and the Reynolds number of the wing is from 367,500 to 183,750.The flapping period is 0.33 s, the wing-tip Strouhal number St = 0.167, and the wing-tip reduced frequency k = 0.0952 can be calculated.As shown in Figure 5a, the wing-root and wing-tip airfoils are NACA8412 and NACA0012, respectively, and the remaining portion of the wing has a linearly interpolated transition from the wing root to the wing tip according to the parameters of the NACA airfoil family.The root airfoil's angle of incidence is 4 degrees.The airfoil must function within an appropriate angle of attack if the flapping wing is intended to generate enough lift and thrust and have good propulsion efficiency [18,31].To balance the manufacturing difficulty and ensure that the wing has a suitable torsion law, it is assumed that there is a linear relationship between the beam's cross-sectional diameter and wingspan; thus, the torsional stiffness is supposed to vary with the wingspan to the fourth power.As shown in Figure 6b, the torsional stiffness GJ represents torsional stiffnesses along the outer wing.Since this article discusses the effect of spanwise folding on a flapping wing, the bending stiffness is assumed to be rigid to exclude interference.

Structural and Kinematic Parameters
The shape of the wing is inspired by the Canadian goose, which has a half wingspan of 0.8 m, an aspect ratio of 12, a root chord length of 0.3 m, and a root-to-tip ratio of 2. The flight velocity is 15 m/s, and the Reynolds number of the wing is from 367,500 to 183,750.The flapping period is 0.33 s, the wing-tip Strouhal number St = 0.167, and the wing-tip reduced frequency k = 0.0952 can be calculated.As shown in Figure 5a, the wing-root and wing-tip airfoils are NACA8412 and NACA0012, respectively, and the remaining portion of the wing has a linearly interpolated transition from the wing root to the wing tip according to the parameters of the NACA airfoil family.The root airfoil's angle of incidence is 4 degrees.The airfoil must function within an appropriate angle of attack if the flapping wing is intended to generate enough lift and thrust and have good propulsion efficiency [18,31].To balance the manufacturing difficulty and ensure that the wing has a suitable torsion law, it is assumed that there is a linear relationship between the beam's cross-sectional diameter and wingspan; thus, the torsional stiffness is supposed to vary with the wingspan to the fourth power.As shown in Figure 6b, the torsional stiffness GJ represents torsional stiffnesses along the outer wing.Since this article discusses the effect of spanwise folding on a flapping wing, the bending stiffness is assumed to be rigid to exclude interference.
where  _ is the flapping amplitude of the inner wing,   is the mean folding angle,  _ is the folding angle amplitude of the outer wing, and   is the phase difference between the flapping angles of the inner and outer wings.In order to show the characteristics of folding flapping wings,  _ = 30 ° is used, and the other parameters are studied with different values.

Results and Discussion
To study the impact of different parameters on flight performance, a kinematic analysis of folding flapping wings with different kinematic parameters was conducted based on real bird flight postures.The flapping phase angles   are 0, 45, and 90 degrees, the inner wing flapping angles  _ are 20 and 30 degrees, and the average folding angles   are 0 and 30 degrees.The different combinations are Case 1 to Case 9, as shown in Table 1.The trajectories of the folding flapping wing during flapping are different due to different kinematic parameters.Figure 7a shows the trajectory of the downward-flapping The lengths of the inner and outer sections are, respectively, 0.45 and 0.55 of the total span.For the convenience of discussion, the flapping motion of the inner/outer wing is a simple harmonic motion.As shown in Figure 6a, the flapping motion law of the inner wing is and the flapping motion law of the outer wing relative to the inner wing is where β in_max is the flapping amplitude of the inner wing, β mean is the mean folding angle, β out_max is the folding angle amplitude of the outer wing, and ϕ β is the phase difference between the flapping angles of the inner and outer wings.In order to show the characteristics of folding flapping wings, β out_max = 30 • is used, and the other parameters are studied with different values.

Results and Discussion
To study the impact of different parameters on flight performance, a kinematic analysis of folding flapping wings with different kinematic parameters was conducted based on real bird flight postures.The flapping phase angles ϕ β are 0, 45, and 90 degrees, the inner wing flapping angles β in_max are 20 and 30 degrees, and the average folding angles β mean are 0 and 30 degrees.The different combinations are Case 1 to Case 9, as shown in Table 1.It can be seen from the density of the trajectory that the Case 1 and 2 velocities will be greater in the rear part of the downward and upward flapping than that in the front part.
When the maximum flapping angle  _ is large, the overall flapping angle is large, and the unsteadiness of the entire wing is enhanced.

Velocity of Wing Tip
The wing is composed of two sections, and the outer section of the wing rotates around the end of the inner section of the wing.As shown in Figure 8a, the flapping velocity at the wing tip is the sum of the migration velocity  ⃗⃗  and the relative velocity  ⃗⃗  .The coordinate of the wing tip at time t is  ⃗⃗  = (  ,   ), and the migration velocity of the wing tip linked to the local coordinate system of the inner wing in the current wing configuration is as follows: It can be seen that the average folding angle β mean affects the configuration of the folded wing during flapping.When β mean = 0 • , the configuration of the wing at the highest and lowest points is symmetrical with respect to half of the inner wing flapping angle, as in Cases 1, 3, and 5.It is asymmetric when β mean is 30 • , as in Cases 2, 4, and 6.We expect that the trajectory will be different during upward and downward flapping.If the outer wing folds downward during upward flapping, the upward speed of the wing can be reduced, and the loss of lift during upward flapping can be reduced.The configurations of real large birds are asymmetric during flight.
As shown in Cases 1, 5, and 9, when the flapping angular phase ϕ β is small, the range of movement and the unsteadiness of the outer wing become large.And the smaller ϕ β is, the more similar the trajectories of the wing are during upward and downward flapping.For example, the upward and downward trajectories in Cases 9 and 10 are exactly the same, while the upward and downward trajectories in Cases 1 and 2 are different.It can be seen from the density of the trajectory that the Case 1 and 2 velocities will be greater in the rear part of the downward and upward flapping than that in the front part.
When the maximum flapping angle β in_max is large, the overall flapping angle is large, and the unsteadiness of the entire wing is enhanced.

Velocity of Wing Tip
The wing is composed of two sections, and the outer section of the wing rotates around the end of the inner section of the wing.As shown in Figure 8a, the flapping velocity at the wing tip is the sum of the migration velocity → V m and the relative velocity → V r .The coordinate of the wing tip at time t is → L tip = (y tip , z tip , and the migration velocity of the wing tip linked to the local coordinate system of the inner wing in the current wing configuration is as follows: The relative velocity of the outer wing relative to the inner wing is where The sum of

Aerodynamic Performance
For Case 1 to Case 10, the changes in the lift coefficient, thrust coefficient, and propulsion efficiency with the initial geometric twist angle of the wing were studied.The initial geometric twist angle is the twist change in the airfoil's angle of incidence from the wing tip to the wing root during initial manufacturing, and the initial geometric twist angle varies linearly along the wingspan.Figure 8b shows the flapping velocity of the wing tip with different kinematic parameters, where the positive and negative are based on the z-direction velocity.In general, the curve of the velocity change with time is similar to a sinusoidal curve.It can be seen that the larger the phase angle ϕ β , the smaller the velocity peak, which is similar to the trajectory analysis.And the larger the ϕ β , the more lagging the phase.The smaller ϕ β is, the closer it is to a simple synthesis of the inner and outer wing velocities.According to the comparison between Cases 1 and 2, when β mean = 30 • , the velocity peak value is slightly less than or equal to that at β mean = 0 • , but the difference is tiny.Obviously, when β n_max is larger, the wing-tip velocity is larger.In general, the phase angle ϕ β has the greatest impact on the wing-tip velocity.

Aerodynamic Performance
For Case 1 to Case 10, the changes in the lift coefficient, thrust coefficient, and propulsion efficiency with the initial geometric twist angle of the wing were studied.The initial geometric twist angle is the twist change in the airfoil's angle of incidence from the wing tip to the wing root during initial manufacturing, and the initial geometric twist angle varies linearly along the wingspan.
As shown in Figure 9a, the lift coefficient increases as the initial geometric twist angle increases; this occurs because when the geometric twist angle increases, the average angle of attack of the airfoil increases.In Cases 1, 3, 5, 7, and 9, it can be seen that when β mean = 0 • , except for the initial geometric twist angle, the parameters have little impact on the lift coefficient.This occurs because when β mean = 0 • , the configuration of the two-section wing in the upward and downward phases is symmetrical.On the contrary, it can be seen in Cases 2, 4, 6, and 8 that the asymmetry of the wing configuration during the flapping period has a greater impact on the lift coefficient.The lift coefficient of the asymmetric configuration is larger, which is consistent with the analysis in Figure 7. Due to the difference in the inner wing flapping angle β in_max , the flapping trajectories during the cycle are different due to the asymmetry caused by β mean .When β in_max is larger, the lift coefficient will be larger, as shown in Cases 2 and 4. Phase ϕ β has no obvious effect on the lift coefficient.
As shown in Figure 9b, the thrust coefficient first increases and then decreases as the initial geometric twist angle increases.This occurs because the airfoil has the strongest thrust generation capability near an average angle of attack of 0 degrees.The thrust coefficients are obviously divided into three levels, with the highest in Cases 7 and 8.The thrust coefficients in Cases 5 and 6 are slightly smaller than those in Cases 3, 4, 9, and 10, and those in Cases 1 and 2 are the smallest.It can be seen that β mean has a small impact on the thrust coefficient; the flapping angle β in_max and phase ϕ β have a greater impact on the thrust coefficient.This is consistent with the discussion of the wing-tip flapping velocity.The thrust coefficients in Cases 3 and 4 are larger than those in Cases 1 and 2 with the same ϕ β , which means that the thrust coefficient is proportional to the β in_max angle.This occurs because a large flapping angle indicates a high flapping velocity and strong unsteadiness.With the same β in_max , the thrust coefficients in Cases 5 and 6 are larger than those in Cases 1 and 2. Similar to the wing-tip flapping velocity analysis, a smaller phase angle ϕ β of the outer wing means a higher flapping velocity.The thrust coefficient is generally proportional to the flapping velocity of the outer wing.
As shown in Figure 9c, the propulsion efficiency first increases and then decreases as the initial geometric twist angle increases.The propulsion efficiency distribution is similar to the distribution of the thrust coefficient.The propulsion efficiency is high when the thrust coefficient is big.A high thrust coefficient means a large maximum angle of attack of the airfoil during the period.At high angles of attack, the lift-to-drag ratio of the airfoil will drop sharply because of the airflow separation.The aerodynamic calculation method in this article is an inviscid method that cannot simulate separation, but with the help of the airfoil surface pressure, this can be discussed.
to the distribution of the thrust coefficient.The propulsion efficiency is high when the thrust coefficient is big.A high thrust coefficient means a large maximum angle of attack of the airfoil during the period.At high angles of attack, the lift-to-drag ratio of the airfoil will drop sharply because of the airflow separation.The aerodynamic calculation method in this article is an inviscid method that cannot simulate separation, but with the help of the airfoil surface pressure, this can be discussed.

Parameters within a Cycle
In order to further study the influence of kinematic parameters on folding flapping wings, the changes in aerodynamic parameters during the period were examined.The initial geometric twist angle of the wing was set to −4 degrees.
As shown in Figure 10a, the lift coefficient curve during the entire period is similar to a sinusoidal curve.The lift coefficient is generally positive, but it is negative for a short period during the upward phase.According to Cases 1, 5, and 9 or Cases 2, 6, and 10, the larger the phase angle   , the more lagging the phase of the lift coefficient, and the smaller the peak value of the lift coefficient.This is consistent with the wing-tip velocity shown in Figure 8b.Comparing Cases 1 and 2, Cases 5 and 6, or Cases 9 and 10, it can be seen that when   = 30 °, the phase of the maximum lift coefficient is advanced over that when   = 0 °, but the phase of the minimum lift coefficient lags behind that when   = 0 °.And it can be seen that the lift coefficient is larger when   = 30 °, which is consistent with the conclusion in Figure 10a.The   significantly affects the lift coefficient during the period.During the upward-flapping phase, the lift loss is reduced due

Parameters within a Cycle
In order to further study the influence of kinematic parameters on folding flapping wings, the changes in aerodynamic parameters during the period were examined.The initial geometric twist angle of the wing was set to −4 degrees.
As shown in Figure 10a, the lift coefficient curve during the entire period is similar to a sinusoidal curve.The lift coefficient is generally positive, but it is negative for a short period during the upward phase.According to Cases 1, 5, and 9 or Cases 2, 6, and 10, the larger the phase angle ϕ β , the more lagging the phase of the lift coefficient, and the smaller the peak value of the lift coefficient.This is consistent with the wing-tip velocity shown in Figure 8b.Comparing Cases 1 and 2, Cases 5 and 6, or Cases 9 and 10, it can be seen that when β mean = 30 • , the phase of the maximum lift coefficient is advanced over that when β mean = 0 • , but the phase of the minimum lift coefficient lags behind that when β mean = 0 • .And it can be seen that the lift coefficient is larger when β mean = 30 • , which is consistent with the conclusion in Figure 10a.The β mean significantly affects the lift coefficient during the period.During the upward-flapping phase, the lift loss is reduced due to the wing's folding motion.It can be seen from Cases 1 and 3 or Cases 5 and 7 that a larger β in_max means a larger range of the lift coefficient curve, and β in_max has little impact on the change in phase.
when   = 0 ° in the downward phase but is smaller than that when   = 0 ° in the upward phase.This shows that the asymmetric motion trajectory reduces the thrust during the upward-flapping phase and concentrates the generation of thrust in the downward-flapping phase.It can be seen from Cases 1 and 3 or Cases 5 and 7 that a larger  _ means a larger range of the thrust coefficient curve due to the increased overall unsteadiness.And  _ has little impact on the thrust coefficient curve phase.As shown in Figure 11a, there is some relationship between the wing-tip twist angle and lift coefficient curves.For different cases, when the range of the lift coefficient is large, the range of the negative wing-tip twist angle will also be large.And the phase of the wing-tip twist angle is consistent with the lift coefficient.The twist angle is positive during the upward phase of flapping and negative during the downward phase for the wing's lift center is located behind the torsion center.The negative twist angles are significantly bigger than the positive twist angles throughout the period because of the angle of incidence and camber of the wing, which causes the positive lift to be substantially greater than the negative lift.The outer wing's maximum twist angles are displayed along its wingspan in Figure 11b, where it can be seen that the twist angles are almost linear.Figure 11c shows the maximum effective angle of attack along the wingspan.It can be seen that the effective angle of attack of the inner wing increases linearly along the span, and due to the torsional deformation of the outer wing, there is a sudden change in the effective angle of attack of the outer wing and the inner wing.As shown in Figure 8b, there exists a correlation between the maximum equivalent angle of attack and the wing-tip flapping velocity, where a higher wing-tip flapping velocity corresponds to a higher effective angle of attack.The maximum effective angle of attack does not exceed 15 degrees, which shows that the method in this article is effective.As shown in Figure 10b, the thrust coefficient curves have two peaks during the period.The peak value in the downward phase is greater than the peak value in the upward phase.The thrust coefficient is generally positive.Comparing Cases 1, 5, and 9 or Cases 2, 6, and 10, it can be seen that, similar to the case of the lift coefficient, the smaller the flapping phase angle ϕ β , the larger the thrust coefficient, and the more lagging the phase of the thrust coefficient.Comparing Cases 1 and 2, Cases 5 and 6, or Cases 9 and 10, it can be seen that when β mean = 30 • and ϕ β ̸ = 0 • , the thrust coefficient is slightly larger than that when β mean = 0 • in the downward phase but is smaller than that when β mean = 0 • in the upward phase.This shows that the asymmetric motion trajectory reduces the thrust during the upward-flapping phase and concentrates the generation of thrust in the downward-flapping phase.It can be seen from Cases 1 and 3 or Cases 5 and 7 that a larger β in_max means a larger range of the thrust coefficient curve due to the increased overall unsteadiness.And β in_max has little impact on the thrust coefficient curve phase.
As shown in Figure 11a, there is some relationship between the wing-tip twist angle and lift coefficient curves.For different cases, when the range of the lift coefficient is large, the range of the negative wing-tip twist angle will also be large.And the phase of the wing-tip twist angle is consistent with the lift coefficient.The twist angle is positive during the upward phase of flapping and negative during the downward phase for the wing's lift center is located behind the torsion center.The negative twist angles are significantly bigger than the positive twist angles throughout the period because of the angle of incidence and camber of the wing, which causes the positive lift to be substantially greater than the negative lift.The outer wing's maximum twist angles are displayed along its wingspan in Figure 11b, where it can be seen that the twist angles are almost linear.Figure 11c shows the maximum effective angle of attack along the wingspan.It can be seen that the effective angle of attack of the inner wing increases linearly along the span, and due to the torsional deformation of the outer wing, there is a sudden change in the effective angle of attack of the outer wing and the inner wing.As shown in Figure 8b, there exists a correlation between the maximum equivalent angle of attack and the wing-tip flapping velocity, where a higher wing-tip flapping velocity corresponds to a higher effective angle of attack.The maximum effective angle of attack does not exceed 15 degrees, which shows that the method in this article is effective.

Pressure Coefficient
In order to further study the influence of kinematic parameters on flight performance, the surface pressure coefficient of the airfoil at 75% of the span during the period was studied.Since the inner section of the wing does not twist, the pressure coefficient of the inner wing is not discussed.In order to further show the characteristics of the airfoil surface pressure coefficient, the pressure coefficient is studied at the maximum and minimum points of the lift coefficient curves, that is,   = ⁄ 0.25 +   /2 and   = ⁄ 0.75 +   /2.Due to the unsteadiness and the asymmetry of the folded-wing motion, the maximum and minimum lift coefficient points have phase angle deviations relative to the inner wing flapping angle, which are   and   , respectively, as shown in Figure 10a.The other two moments in the cycle are set to   = ⁄ 0.5 + (  +   ) 2 ⁄  and   = ⁄ 1 + (  +   ) 2 ⁄ .As shown in Figure 12, based on kinematic laws and the pressure coefficient distribution, the suction peak of the pressure coefficient distribution means that the positive effective angle of attack of the airfoil is large when   = ⁄ 0.25 +   /2 in Cases 1 and 2, and the negative effective angle of attack is large when   = ⁄ 0.75 +   /2.At the other two moments, the effective angle of attack of the airfoil is very small, generating almost no lift and thrust.Comparing Figure 12a,b, it can be seen that due to the different   , the effective angle of attack of the airfoil in the cycle is greatly affected.Combined with Figures 8b and 9a, it shows that although   has a small impact on the velocity, it can affect the effective angle of attack of the outer wing.Moreover, the average angle of

Pressure Coefficient
In order to further study the influence of kinematic parameters on flight performance, the surface pressure coefficient of the airfoil at 75% of the span during the period was studied.Since the inner section of the wing does not twist, the pressure coefficient of the inner wing is not discussed.In order to further show the characteristics of the airfoil surface pressure coefficient, the pressure coefficient is studied at the maximum and minimum points of the lift coefficient curves, that is, t/T = 0.25 + ϕ max /2π and t/T = 0.75 + ϕ min /2π.Due to the unsteadiness and the asymmetry of the folded-wing motion, the maximum and minimum lift coefficient points have phase angle deviations relative to the inner wing flapping angle, which are ϕ max and ϕ min , respectively, as shown in Figure 10a.The other two moments in the cycle are set to t/T = 0.5 + (ϕ max + ϕ min )/2π and t/T = 1 + (ϕ max + ϕ min )/2π.
As shown in Figure 12, based on kinematic laws and the pressure coefficient distribution, the suction peak of the pressure coefficient distribution means that the positive effective angle of attack of the airfoil is large when t/T = 0.25 + ϕ max /2π in Cases 1 and 2, and the negative effective angle of attack is large when t/T = 0.75 + ϕ min /2π.At the other two moments, the effective angle of attack of the airfoil is very small, generating almost no lift and thrust.Comparing Figure 12a,b, it can be seen that due to the different β mean , the effective angle of attack of the airfoil in the cycle is greatly affected.Combined with Figures 8b and 9a, it shows that although β mean has a small impact on the velocity, it can affect the effective angle of attack of the outer wing.Moreover, the average angle of attack during the period when β mean = 30 • in Case 2 is biased toward the positive angle of attack compared to that when β mean = 0 • in Case 1, which makes the lift coefficient larger.
attack during the period when   = 30 ° in Case 2 is biased toward the positive angle of attack compared to that when   = 0 ° in Case 1, which makes the lift coefficient larger.
(a) (b) As shown in Figure 13, the effective positive angle of attack at   = ⁄ 0.25 +   /2 and the negative angle of attack at   = ⁄ 0.75 +   /2 in Cases 3 and 4 both increase significantly compared with those in Cases 1 and 2. The problem in Cases 3 and 4 is that there is always a moment when the effective angle of attack is excessively large during the upward or downward phase.This will aggravate the separation of the trailing edge of the airfoil, and the lift-to-drag ratio will drop sharply, thus reducing the propulsion efficiency.This situation should be avoided by adjusting the initial geometric twist angle.As shown in Figure 14, the airfoil surface pressure coefficient in Cases 5 and 6 represents that the effective-angle-of-attack range of the airfoil is larger than that in Cases 1 and 2. Due to the smaller flapping phase angle   , the flapping velocity of the outer wing is larger than that in Cases 1 and 2, which is consistent with the analysis in Figure 8b.And it can also be seen in Figure 9 that the thrust coefficient and propulsion efficiency are much larger than those in Cases 1 and 2. Comparing Cases 5 and 6 with Cases 1 and 2, the impact of   is the same.As shown in Figure 13, the effective positive angle of attack at t/T = 0.25 + ϕ max /2π and the negative angle of attack at t/T = 0.75 + ϕ min /2π in Cases 3 and 4 both increase significantly compared with those in Cases 1 and 2. The problem in Cases 3 and 4 is that there is always a moment when the effective angle of attack is excessively large during the upward or downward phase.This will aggravate the separation of the trailing edge of the airfoil, and the lift-to-drag ratio will drop sharply, thus reducing the propulsion efficiency.This situation should be avoided by adjusting the initial geometric twist angle.As shown in Figure 13, the effective positive angle of attack at   = ⁄ 0.25 +   /2 and the negative angle of attack at   = ⁄ 0.75 +   /2 in Cases 3 and 4 both increase significantly compared with those in Cases 1 and 2. The problem in Cases 3 and 4 is that there is always a moment when the effective angle of attack is excessively large during the upward or downward phase.This will aggravate the separation of the trailing edge of the airfoil, and the lift-to-drag ratio will drop sharply, thus reducing the propulsion efficiency.This situation should be avoided by adjusting the initial geometric twist angle.As shown in Figure 14, the airfoil surface pressure coefficient in Cases 5 and 6 represents that the effective-angle-of-attack range of the airfoil is larger than that in Cases 1 and 2. Due to the smaller flapping phase angle   , the flapping velocity of the outer wing is larger than that in Cases 1 and 2, which is consistent with the analysis in Figure 8b.And it can also be seen in Figure 9 that the thrust coefficient and propulsion efficiency are much larger than those in Cases 1 and 2. Comparing Cases 5 and 6 with Cases 1 and 2, the impact of   is the same.As shown in Figure 14, the airfoil surface pressure coefficient in Cases 5 and 6 represents that the effective-angle-of-attack range of the airfoil is larger than that in Cases 1 and 2. Due to the smaller flapping phase angle ϕ β , the flapping velocity of the outer wing is larger than that in Cases 1 and 2, which is consistent with the analysis in Figure 8b.And it can also be seen in Figure 9 that the thrust coefficient and propulsion efficiency are much larger than those in Cases 1 and 2. Comparing Cases 5 and 6 with Cases 1 and 2, the impact of β mean is the same.As shown in Figure 15, as  _ becomes larger, the pressure coefficient shows that the maximum effective angle of attack in Cases 7 and 8 becomes excessively large compared to that in Cases 5 and 6.From a practical perspective, such a large effective angle of attack of the outer wing is not allowed during the flapping process [18].As shown in Figure 16, in Cases 9 and 10, the range of the effective angle of attack is large when the phase angle   = 0 ° because the flapping velocity is large, according to Figure 8b.But according to Figure 9, its thrust and propulsion efficiency are smaller than those in Cases 3 and 4. Therefore, after paying the price of a large airfoil angle of attack and severe airflow separation, Cases 9 and 10 have no advantages over Cases 3 and 4 in the lift coefficient and thrust coefficient.When the flapping phase angle   = 0 °, although it is beneficial to improve the flapping velocity of the outer wing, the outer wing is prone to have an excessive angle of attack, and the thrust coefficient and propulsion efficiency are not guaranteed.As shown in Figure 15, as β in_max becomes larger, the pressure coefficient shows that the maximum effective angle of attack in Cases 7 and 8 becomes excessively large compared to that in Cases 5 and 6.From a practical perspective, such a large effective angle of attack of the outer wing is not allowed during the flapping process [18].As shown in Figure 15, as  _ becomes larger, the pressure coefficient shows that the maximum effective angle of attack in Cases 7 and 8 becomes excessively large compared to that in Cases 5 and 6.From a practical perspective, such a large effective angle of attack of the outer wing is not allowed during the flapping process [18].As shown in Figure 16, in Cases 9 and 10, the range of the effective angle of attack is large when the phase angle   = 0 ° because the flapping velocity is large, according to Figure 8b.But according to Figure 9, its thrust and propulsion efficiency are smaller than those in Cases 3 and 4. Therefore, after paying the price of a large airfoil angle of attack and severe airflow separation, Cases 9 and 10 have no advantages over Cases 3 and 4 in the lift coefficient and thrust coefficient.When the flapping phase angle   = 0 °, although it is beneficial to improve the flapping velocity of the outer wing, the outer wing is prone to have an excessive angle of attack, and the thrust coefficient and propulsion efficiency are not guaranteed.As shown in Figure 16, in Cases 9 and 10, the range of the effective angle of attack is large when the phase angle ϕ β = 0 • because the flapping velocity is large, according to Figure 8b.But according to Figure 9, its thrust and propulsion efficiency are smaller than those in Cases 3 and 4. Therefore, after paying the price of a large airfoil angle of attack and severe airflow separation, Cases 9 and 10 have no advantages over Cases 3 and 4 in the lift coefficient and thrust coefficient.When the flapping phase angle ϕ β = 0 • , although it is beneficial to improve the flapping velocity of the outer wing, the outer wing is prone to have an excessive angle of attack, and the thrust coefficient and propulsion efficiency are not guaranteed.According to the above discussion, Cases 3 and 4 have advantages over Cases 9 and 10.Now, compare Cases 3 and 4 with Cases 5 and 6.When the initial geometric twist angle is −4 degrees, there are some differences between the pressure coefficients in Figures 13  and 14.In order to make the comparison more convincing, the geometric twist angle in Case 4 is set to −7 degrees.As shown in Figures 14b and 17, the range of the airfoil pressure coefficient is almost the same, which means that the maximum effective angle of attack in the cycle is almost the same.According to Figure 9, it can be seen that the lift coefficient, thrust coefficient, and propulsion efficiency in Case 4 when the initial geometric twist angle is −7 degrees are greater than those in Case 6 when the initial geometric twist angle is −4 degrees.It can be concluded that a larger flapping phase angle   has benefits, which can reduce the lift loss in the upward-flapping process of the folding flapping wing.The above-discussed results can be compared in Table 2 It can be seen that, excluding Cases 7 and 8, whose equivalent angle of attack is excessively large, Case 4, with a geometric twist angle of −7 degrees, has advantages over the other cases in terms of the lift coefficient, thrust coefficient, and propulsion efficiency.The same conclusion can also be learned from Figure 9.According to the above discussion, Cases 3 and 4 have advantages over Cases 9 and 10.Now, compare Cases 3 and 4 with Cases 5 and 6.When the initial geometric twist angle is −4 degrees, there are some differences between the pressure coefficients in Figures 13  and 14.In order to make the comparison more convincing, the geometric twist angle in Case 4 is set to −7 degrees.As shown in Figures 14b and 17, the range of the airfoil pressure coefficient is almost the same, which means that the maximum effective angle of attack in the cycle is almost the same.According to Figure 9, it can be seen that the lift coefficient, thrust coefficient, and propulsion efficiency in Case 4 when the initial geometric twist angle is −7 degrees are greater than those in Case 6 when the initial geometric twist angle is −4 degrees.It can be concluded that a larger flapping phase angle ϕ β has benefits, which can reduce the lift loss in the upward-flapping process of the folding flapping wing.According to the above discussion, Cases 3 and 4 have advantages over Cases 9 and 10.Now, compare Cases 3 and 4 with Cases 5 and 6.When the initial geometric twist angle is −4 degrees, there are some differences between the pressure coefficients in Figures 13  and 14.In order to make the comparison more convincing, the geometric twist angle in Case 4 is set to −7 degrees.As shown in Figures 14b and 17, the range of the airfoil pressure coefficient is almost the same, which means that the maximum effective angle of attack in the cycle is almost the same.According to Figure 9, it can be seen that the lift coefficient, thrust coefficient, and propulsion efficiency in Case 4 when the initial geometric twist angle is −7 degrees are greater than those in Case 6 when the initial geometric twist angle is −4 degrees.It can be concluded that a larger flapping phase angle   has benefits, which can reduce the lift loss in the upward-flapping process of the folding flapping wing.The above-discussed results can be compared in Table 2 It can be seen that, excluding Cases 7 and 8, whose equivalent angle of attack is excessively large, Case 4, with a geometric twist angle of −7 degrees, has advantages over the other cases in terms of the lift coefficient, thrust coefficient, and propulsion efficiency.The same conclusion can also be learned from Figure 9.The above-discussed results can be compared in Table 2 It can be seen that, excluding Cases 7 and 8, whose equivalent angle of attack is excessively large, Case 4, with a geometric twist angle of −7 degrees, has advantages over the other cases in terms of the lift coefficient, thrust coefficient, and propulsion efficiency.The same conclusion can also be learned from Figure 9.

Conclusions
This article proposes a method for the rapid and accurate calculation of folding passive torsional flapping wings during the early design stage.The unsteady three-dimensional surface element method is employed to compute aerodynamic forces, the linear beam element method is utilized to calculate torsional deformation, and the fluid-structure coupling calculation is conducted using the Newmark-β method.Using various kinematic parameters, this article investigates the changing rules of the lift coefficient, thrust coefficient, and propulsion efficiency as the initial geometric twist angle changes.It analyzes the lift coefficient, thrust coefficient, and torsional deformation of the wing throughout the entire period, along with the wing surface pressure coefficient at different moments.This article draws the following conclusions: 1.
Kinematic parameters significantly impact performance.The flapping phase angle ϕ β between the inner and outer wings affects the movement velocity of the outer wing.A larger ϕ β implies a smaller flapping velocity for the outer wing, resulting in a smaller thrust coefficient and a more pronounced lag in the lift coefficient phase.A Larger ϕ β also indicates asymmetry in the flapping trajectory.

2.
The average folding angle β mean significantly influences the lift coefficient as it impacts the trajectory of the flapping wing.When β mean is positive, the folding motion is noticeable during the upward-flapping process, resulting in small lift loss and, consequently, a larger lift coefficient during this phase.It has a minor impact on the thrust coefficient and propulsion efficiency but greatly affects the effective angle of attack of the outer wing.β mean should be designed to align with the initial geometric twist angle.Otherwise, the airfoil's angle of attack throughout the period might become excessive, leading to severe airflow separation and reduced propulsion efficiency.

3.
The flapping angle β in_max of the inner wing primarily influences the overall wing's unsteadiness.An increase in β in_max results in an elevation in the thrust coefficient, but it may lead to an excessive effective angle of attack for the outer wing, reducing the propulsion efficiency.4.
For folding flapping wings, there are principles for selecting kinematic parameters.
From the perspective of the lift coefficient, the folding motion should be applied to reduce lift loss during the upward-flapping phase.The configuration of the wing during the period should be asymmetrical, and a large flapping phase angle ϕ β and a positive average folding angle β mean should be selected.Regarding the thrust coefficient and propulsion efficiency, the generation of thrust should be concentrated in the downward-flapping phase when the lift-to-drag ratio is high.So, a positive β mean should be selected.The inner wing flapping angle β in_max can adjust the overall unsteadiness of the wing and the thrust.The initial geometric twist angle can be matched with the average folding angle β mean to maintain the airfoil's effective angle of attack within a reasonable range during upward and downward flapping.

Figure 4 .
Figure 4. Comparison of calculation results with Delaurier's calculation and experimental results.(a) Average lift; (b) average thrust; (c) torsion amplitude and phase angle; (d) flapping elastic amplitude and phase angle.

Figure 4 .
Figure 4. Comparison of calculation results with Delaurier's calculation and experimental results.(a) Average lift; (b) average thrust; (c) torsion amplitude and phase angle; (d) flapping elastic amplitude and phase angle.

Figure 5 .
Figure 5. Wing model.(a) Structure of the wing; (b) forces and moments on the airfoil.

Figure 5 .
Figure 5. Wing model.(a) Structure of the wing; (b) forces and moments on the airfoil.

Figure 6 .
Figure 6.Structural and kinematic parameters of folding flapping wing.(a) Flapping angle of folding wing; (b) torsional stiffness of outer wing.

Figure 6 .
Figure 6.Structural and kinematic parameters of folding flapping wing.(a) Flapping angle of folding wing; (b) torsional stiffness of outer wing.

4. 1 .
Kinematic Analysis 4.1.1.Flapping Trajectory of the Folding Flapping Wing The trajectories of the folding flapping wing during flapping are different due to different kinematic parameters.Figure 7a shows the trajectory of the downward-flapping stage in Cases 1 to 10, and Figure 7b shows the trajectory of the upward-flapping stage in Cases 1 to 10.
Biomimetics 2024, 9, 42 9 of 20  flapping.For example, the upward and downward trajectories in Cases 9 and 10 are exactly the same, while the upward and downward trajectories in Cases 1 and 2 are different.

Figure 7 .
Figure 7.The trajectory of folding wing during the period.(a) The trajectory of downward flapping from Case 1 to Case 10; (b) the trajectory of upward flapping from Case 1 to Case 10.

Figure 7 .
Figure 7.The trajectory of folding wing during the period.(a) The trajectory of downward flapping from Case 1 to Case 10; (b) the trajectory of upward flapping from Case 1 to Case 10.

Figure 11 .
Figure 11.Twist angle and effective angle of attack along the span.(a) Wing-tip twist angle during the period; (b) maximum twist angles along the wingspan; (c) maximum effective angle of attack along the wingspan.

Figure 11 .
Figure 11.Twist angle and effective angle of attack along the span.(a) Wing-tip twist angle during the period; (b) maximum twist angles along the wingspan; (c) maximum effective angle of attack along the wingspan.

Figure 12 .
Figure 12.Pressure coefficient at 75% of the wingspan in Case 1 and Case 2. (a) Pressure coefficient in Case 1; (b) pressure coefficient in Case 2.

Figure 13 .
Figure 13.Pressure coefficient at 75% of the wingspan in Case 3 and Case 4. (a) Pressure coefficient in Case 3; (b) pressure coefficient in Case 4.

Figure 12 .
Figure 12.Pressure coefficient at 75% of the wingspan in Case 1 and Case 2. (a) Pressure coefficient in Case 1; (b) pressure coefficient in Case 2.

Figure 13 .
Figure 13.Pressure coefficient at 75% of the wingspan in Case 3 and Case 4. (a) Pressure coefficient in Case 3; (b) pressure coefficient in Case 4.

Figure 13 .
Figure 13.Pressure coefficient at 75% of the wingspan in Case 3 and Case 4. (a) Pressure coefficient in Case 3; (b) pressure coefficient in Case 4.

Figure 14 .
Figure 14.Pressure coefficient at 75% of the wingspan in Case 5 and Case 6.(a) Pressure coefficient in Case 5; (b) pressure coefficient in Case 6.

Figure 15 .
Figure 15.Pressure coefficient at 75% of the wingspan in Case 7 and Case 8. (a) Pressure coefficient in Case 7; (b) pressure coefficient in Case 8.

Figure 14 .
Figure 14.Pressure coefficient at 75% of the wingspan in Case 5 and Case 6.(a) Pressure coefficient in Case 5; (b) pressure coefficient in Case 6.

Figure 14 .
Figure 14.Pressure coefficient at 75% of the wingspan in Case 5 and Case 6.(a) Pressure coefficient in Case 5; (b) pressure coefficient in Case 6.

Figure 15 .
Figure 15.Pressure coefficient at 75% of the wingspan in Case 7 and Case 8. (a) Pressure coefficient in Case 7; (b) pressure coefficient in Case 8.

Figure 15 .
Figure 15.Pressure coefficient at 75% of the wingspan in Case 7 and Case 8. (a) Pressure coefficient in Case 7; (b) pressure coefficient in Case 8.

Figure 16 .
Figure 16.Pressure coefficient at 75% of the wingspan in Case 9 and Case 10.(a) Pressure coefficient in Case 9; (b) pressure coefficient in Case 10.

Figure 17 .
Figure 17.Pressure coefficient at 75% of the wingspan in Case 4 with initial geometric twist angle of −7 degrees.

Figure 16 .
Figure 16.Pressure coefficient at 75% of the wingspan in Case 9 and Case 10.(a) Pressure coefficient in Case 9; (b) pressure coefficient in Case 10.

Figure 16 .
Figure 16.Pressure coefficient at 75% of the wingspan in Case 9 and Case 10.(a) Pressure coefficient in Case 9; (b) pressure coefficient in Case 10.

Figure 17 .
Figure 17.Pressure coefficient at 75% of the wingspan in Case 4 with initial geometric twist angle of −7 degrees.

Figure 17 .
Figure 17.Pressure coefficient at 75% of the wingspan in Case 4 with initial geometric twist angle of −7 degrees.
the wing β in flapping angle of inner wing β out flapping angle of outer wing β in_max flapping amplitude of the inner wing β out_max flapping amplitude of the outer wing β mean mean folding angle U tip flapping velocity of the wing tip GJ torsional stiffnesses of the beam EI bending stiffnesses of the beam ϕ β phase difference between the flapping angles of the inner and outer wings ϕ max phase angle deviations relative to the inner wing flapping angle at maximum lift coefficient point ϕ min phase angle deviations relative to the inner wing flapping angle at maximum lift coefficient point θ geo initial geometric twist angle of the wing θ e f f effective angle of attack of wing θ max maximum passive twisting angle of the wing θ tip passive twisting angle of the wing tip H max flapping elastic amplitude of the wing Φ θ phase angles between twisting angle and flapping motion Φ H phase angles between bending response and flapping motion Φ up velocity potential of the upper surface of the trailing edge Φ down velocity potential of the upper surface of the trailing edge k reduced frequency St Strouhal number

Table 1 .
Combinations of different kinematic parameters.

Table 2 .
Comparison of performance parameters.