The Effect of Torsional and Bending Stiffness on the Aerodynamic Performance of Flapping Wing

: For large bird-like flapping wing aircraft, the fluid–structure coupling problem is very important. Through the passive torsional deformation of the wing, sufficient thrust is generated and propulsion efficiency is ensured. Moreover, spanwise bending deformation will affect lift and thrust. The flow field on the surface of the wing and the geometric nonlinearity caused by the large deformation of the wing should be considered during the design process. Existing research methods do not solve this problem accurately and efficiently. This paper provides a method to analyze the fluid–structure coupling problem of the flapping wing which adopts the three-dimensional unsteady panel method to solve the aerodynamic force, and adopts the linear beam element model combined with the corotational formulation method to consider the geometric nonlinear deformation of the wing beam. This article compares the performance of the flapping wing with different torsional and bending stiffness, and analyzes the airfoil surface pressure coefficients at different portions of the wing during the period. The results show that torsional stiffness has a large influence on the lift coefficient, thrust coefficient and propulsion efficiency. Meanwhile, the torsional stiffness of the wing beam and the initial geometric twist angle of the wing need to be well coordinated to achieve high efficiency. Moreover, appropriate bending stiffness of the wing is conducive to improving propulsion efficiency.


Introduction
Compared with traditional fixed-wing aircraft and helicopters, bird-like flapping wing aircrafts have their own characteristics.The flexibility of the flapping wing significantly affects the aerodynamic force of the wing [1], and the passive flexible deformation can be exploited to improve efficiency.We hope that the fluid-structure coupling effect of the wing will cause the wing to periodically twist and bend, thus achieving an increase in thrust and propulsive efficiency.This may not result in optimal efficiency compared to active control, but it can eliminate the need for control devices, greatly simplifying the wing structure and reducing weight.
Previous research has been done on the fluid-structure coupling effect of flapping wings.DeLaurier [2] designed a large flapping wing that utilized the passive twisting of the wings to improve aerodynamic performance.The two-dimensional analytical method was used to solve the unsteady aerodynamic force, the wing mass was discretized into concentrated mass points, and the dynamic problems were jointly solved by combining the deformation of the wing.Smith [3] established a fluid-structure interaction model by the panel method with unsteady wake and the beam finite element model to simulate the moth's flight.The results showed that the model is consistent with the real flight state of the moth.Isogai [4] employed lifting surface theory to solve unsteady aerodynamic forces, and utilized the superposition of states of beams to solve dynamical problems to optimize the design of Delaurier's flapping wing.However, the influence of wing airfoil and camber was not considered.The results showed that to achieve optimal efficiency, the beam should have a large bending deformation.Larijani [5] considered airflow separation in the aerodynamic model based on Delaurier's model, adopted the finite element beam model to simulate the deformation of the beam, and achieved better accuracy than Delaurier.Reichert [6] used a simple surface element method to solve the aerodynamic forces, and a three-dimensional beam model was used to account for the geometric nonlinearity of the beam.Finally, an optimization study of the motion parameters was conducted.Wang [7] used the N-S equations method to solve the aerodynamic forces of a flapping rotor vehicle, which utilizes a preset wing deformation to simulate the fluid-structure coupling problem without taking into account the interaction effects.Wu [8] used the unsteady vortex lattice method and the modal superposition method to optimize the planar shape of a large aspect ratio flapping wing.Heathcotel [9] used experiments to study the impact of spanwise flexibility on flapping wings, and found that some spanwise flexibility can improve the propulsion efficiency and thrust of flapping wings.Aono [10] studied the effect of spanwise stiffness on flapping wings using the N-S equation solver and nonlinear beam model.He pointed out that changes in spanwise flexibility will affect the phase of motion and angle of attack, thereby affecting the thrust of the wing.Kodali [11] used the two-dimensional analytical method to solve the unsteady aerodynamic force and the Euler Bernoulli beam to consider the spanwise flexibility, and found that the flapping frequency is near the first spanwise natural frequency of the wing.Lankford [12] conducted a fluid-structure coupling analysis on a flapping wing that can fly vertically and compared it with experiments.The aerodynamic calculation utilized a solver based on the N-S equations, and the structural calculation utilized a universal solver that can calculate large deformations.Yang [13] studied a highly flexible flapping wing in hover.The aerodynamic force was calculated by an improved vortex lattice method and the structure was modeled using a fully nonlinear beam and shell finite elements.Yang [14] studied the force generation mechanism of the dovelike flapping micro aircraft through numerical simulation and experiments.The digital image correlation technology was utilized to measure the dynamic deformation of the wing, and aerodynamic forces were solved by a CFD solver based on N-S equations.
There are some shortcomings in the existing studies.Most studies on the fluidstructure interaction problem of the flapping wing adopted analytical models or singlesurface panel models to calculate aerodynamic forces.Parameters such as the pressure coefficient on the wing surface cannot be calculated, and the working status of the airfoil surface is unclear.In fact, in the Reynolds number range of 10 5 , the airfoil and surface pressure coefficient will significantly affect the aerodynamic results [15].Moreover, the high precision method of solving N-S equations which consumes too much time is not suitable for the initial design calculations and analysis of fluid-structure interaction problems.Few studies have considered the geometric nonlinearity caused by the large bending deformation of the flapping wing.The default assumption is that the wings deform linearly, which is obviously inconsistent with reality.There are few studies on flapping wing under passive torsion and bending coupling at the same time, which is inconsistent with the fact that actual wings are in a bending-torsion coupling state.Most studies only consider the influence of motion parameters and wing stiffness, but do not consider the influence of the geometric parameters, such as the initial geometric twist angle.
To address the problems of the current research, this paper develops a method for solving the fluid-structure coupling problem of the flapping wings.The unsteady threedimensional panel method is used to solve the pressure coefficient of the wing surface, and the linear beam element combined with the corotational formulation method is used to consider the geometric nonlinear problem of large deformation of the wing.A dynamic analysis is conducted on the passive deformed wings with different torsional and bending stiffness under different initial geometric twist angles, and aerodynamic coefficients such as lift coefficient, thrust coefficient and propulsion efficiency are compared.Various aerodynamic coefficients and geometric deformation parameters throughout the period are compared under the bending-torsion coupling state, and the wing surface pressure coefficients at different sections are analyzed.The influence of the torsional and bending stiffness of the wing on the aerodynamic performance is discussed, the mechanism is analyzed and design guidance is given.

Aerodynamic Model
According to Prandtl's boundary layer theory, the flow around the wing can be divided into viscous flow around the wing surface and potential flow.The panel method can be used to calculate the potential flow.Many programs adopted the panel method to calculate the surface pressure of the three-dimensional wing.Magnus [16] adopted the high-order panel method to solve subsonic and supersonic potential flow.Ashby [17] adopted low-order surface source and surface dipole to solve three-dimensional wing aerodynamic forces, and used the integral boundary layer method to estimate the frictional drag.Maskew [18] and Johnson [19] also used the three-dimensional surface element method to calculate wing aerodynamic forces.Vest [20] adopted the low-order three-dimensional panel method combined with the free wake model to solve the aerodynamic force of pigeons.This paper adopts the method based on Hess [21] and improves it for unsteady problems.
As shown in Figure 1, the wing surface adopts two types of panels: surface source and surface dipole, and the panels are all in the form of planes.In its two-dimensional form [22], the surface adopts a constant linear vortex.In three dimensions, it is difficult to deal with the coordination problem of the surface vortex.Considering the equivalence of the surface dipole and surface vortex [23], a linearly changing surface dipole is used instead.The wing is divided into M parts in the chord direction, and the spanwise direction is divided into N parts by N-lines.At each time step, the distribution of each surface source is uniform, and the surface dipole distribution between every two N-lines changes linearly in both spanwise and chordwise directions.After the wing is deformed, the panels will warp, making it no longer a plane.However, by projecting the space quadrilateral onto the plane unit for calculation, there will be no noticeable calculation error [21].The wake is composed of surface dipoles shed from the trailing edge of the wing.According to Helmholtz's law, the total circulation is conserved in the flow field.The wake circulation in the first row of the trailing edge is that of the surface circulation of the wing at k−1 time steps minus the surface circulation of the wing segment at the same spanwise position at k time steps.Over time, the dipole panels continue to fall off from the trailing edge of the wing, producing surface dipoles of varying strengths in the wake.Since the reduced frequency of a large flapping wing aircraft is below 0.2 during cruising flight [24], the unsteadiness is not strong.Similar to Reichert's [6] "frozen-wake" method, the wake in the model used in this paper smoothly moves backward at the speed of the free flow after the trailing edge, without free deformation with the induced velocity.This will not cause perceptible errors and can greatly improve the calculation speed.
Since the dipole changes linearly along the N-lines, the M surface dipoles on the wing between the two N-lines and the  dipoles in the wake are considered together in the calculation, which is equivalent to N unknown numbers for dipoles on the entire wing.Adding the unknown M × N sources, there are a total of M * N + N unknowns.According to the Neumann boundary condition, the velocity in the normal direction at the control point of the wing surface panel is 0, so M * N equations can be established.The velocity of the surface is affected by the induced velocity of the surface panels, the free flow, the flapping motion of the wing, and the elastic deformation motion.The remaining N equations are established through Kutta conditions on the trailing edge.According to the Kutta condition, the pressure on the upper and lower surfaces of the trailing edge of the wing must be equal.For the unsteady form, it is expressed as where  and  are the tangential velocity and velocity potential of the upper and lower surfaces of the trailing edge at k time steps.

Structural Model
For the problem in the article, the wing undergoes large bending deformation, which exceeds the small deformation assumption.Patil [25] found that the geometric nonlinear effects caused by large deformations are very important for the accurate calculation of aeroelastic problems.There have been many studies on large deformation beams.Wang The Neumann boundary condition is used, so it is necessary to compute the induced velocity of the jth panel at the control points of the ith panel.The velocity due to the constant strength source density is denoted → V ij (source).For velocities due to dipoles, the velocity potential is first required.The potential due to the dipole distribution on the panel at points of space is obtained by integrating over the panel: where ξ and η are panel local coordinates, µ(ξ, η) is the dipole distribution on the panel and ϕ unit is the potential of a unit point dipole with axis normal to the panel.The velocity due to the dipole distribution is Since the dipole changes linearly along the N-lines, the M surface dipoles on the wing between the two N-lines and the M wake dipoles in the wake are considered together in the calculation, which is equivalent to N unknown numbers for dipoles on the entire wing.Adding the unknown M × N sources, there are a total of M * N + N unknowns.According to the Neumann boundary condition, the velocity in the normal direction at the control point of the wing surface panel is 0, so M * N equations can be established.The velocity of the surface is affected by the induced velocity of the surface panels, the free flow, the flapping motion of the wing, and the elastic deformation motion.The remaining N equations are established through Kutta conditions on the trailing edge.According to the Kutta condition, the pressure on the upper and lower surfaces of the trailing edge of the wing must be equal.For the unsteady form, it is expressed as where V t and ϕ are the tangential velocity and velocity potential of the upper and lower surfaces of the trailing edge at k time steps.

Structural Model
For the problem in the article, the wing undergoes large bending deformation, which exceeds the small deformation assumption.Patil [25] found that the geometric nonlinear effects caused by large deformations are very important for the accurate calculation of aeroelastic problems.There have been many studies on large deformation beams.Wang [26] adopted nonlinear geometric equations to establish a nonlinear geometric model of the beam and solved the coupling problem of aeroelasticity and flight mechanics of high-aspectratio wing aircraft.Crisfield [27] and Garcia [28] adopted the corotational formulation method to solve the large deformation problem of the three-dimensional beam.The method assumes that the structural deformation is in a linear state in the local reference system, and the nonlinearity is completely introduced by the overall large deformation.
This paper adopts the finite element method based on Kuo-Mo's [29] method to solve the large deformation problem of the beam.Since the beam in this article is in the form of a slender beam, the deformation is small relative to the length, and the strain is still within the linear elastic range, small strain and small angle assumptions are adopted in each element, and a linear beam element is used.Large displacement deformation is caused by the rigid body movement and rotation of the element, and the overall deformation shows large deformation geometric nonlinearity [30].This method is relatively simple to accurately solve the exact beam equations [31].As shown in Figure 2, each element in the method is associated with a local Cartesian coordinate system X E i (i = 1, 2, 3).When deformation occurs, the coordinate system translates and rotates with the element.The solution is performed in the global coordinate system X B i (i = 1, 2, 3), which remains fixed to the wing.
Aerospace 2023, 10, x FOR PEER REVIEW 5 of 25 method is associated with a local Cartesian coordinate system  ( = 1,2,3).When deformation occurs, the coordinate system translates and rotates with the element.The solution is performed in the global coordinate system  ( = 1,2,3), which remains fixed to the wing.According to the above discussion, a three-dimensional linear beam element is used, each unit has 12 degrees of freedom, and each point has 3 translational and 3 rotational degrees of freedom.The local element stiffness matrix is  : The method requires an iterative solution.After p iterations, the point coordinates  ⃗ and node displacements  ⃗ in the global coordinate system are known.Through a series of transformations, the element transformation matrix between the local coordinate system and the global coordinate system can be obtained as  , as well as the strain of the local element  ⃗ = 0,0,0,  ,  ⃗ ,  ⃗ ,  ⃗ , 0,0,  ,  ⃗ ,  ⃗ .
Once the strain of the local element is known, the internal stress of the local element can be calculated: According to the above discussion, a three-dimensional linear beam element is used, each unit has 12 degrees of freedom, and each point has 3 translational and 3 rotational degrees of freedom.The local element stiffness matrix is [k E ]: The method requires an iterative solution.After p iterations, the point coordinates and node displacements → Q B j in the global coordinate system are known.Through a series of transformations, the element transformation matrix between the local coordinate system and the global coordinate system can be obtained as [T] B E , as well as the strain of the local element Once the strain of the local element is known, the internal stress of the local element can be calculated: For the current beam configuration, the element stiffness matrix and internal stress in the local coordinate system needs to be transformed into the global coordinate system: Finally, the element stiffness matrices are aggregated into the overall stiffness matrix: The internal stresses of the element are integrated into the overall internal force vector: For existing external load, it is assumed that the converged solution corresponding to the load has been obtained in p iterations.Now, an external load increment is added and the external load changes to → F B ext .We solve the displacement increment generated by the external load increment in step p + 1 and update the calculation results: Then, the new node displacement is of the nodes can be obtained.New node coordinates are available.
This starts the next iteration until

Aeroelastic Coupling
During the wing flapping process, the unsteady aerodynamic force, the inertial force of the wing mass, the structural damping force and the deformation internal force of the beam reach a dynamic balance.The influence of aerodynamic force causes the deformation of the beam, and the deformation of the beam will affect the aerodynamic force.The iterative calculation is repeated until convergence.Starting at time step k, the equilibrium equations are solved for time step k + 1: where → q k+1 is the deformation displacement, → F k+1 is the force and moment applied to the beam, [K] is the stiffness matrix, [M] is the mass matrix and [D] is the structural damping matrix which can be ignored here.In order to discretize the second order dif-ferential equation, the Newmark-β method is adopted [5].It is an implicit method that is unconditionally stable.

Performance Parameters
In unsteady flows, according to the unsteady Bernoulli theorem [32], the pressure coefficient is expressed as where V move is the motion velocity on the wing relative to the fixed coordinate system, V is the airflow velocity on the wing surface, and ϕ is the velocity potential.The wing performs a pitching motion around the elastic axis, accompanied with a flapping motion around the wing root flapping axis.The moving wing generates horizontal thrust F x (t), vertical lift F z (t), torsional moment M t (t) and bending moment M b (t).M t (t) and M b (t) are the combination of aerodynamic moment and inertial force moment.The power output from the wing is where β(t) is the flapping angle, and M b_root (t) is the bending moment at the wing root.F x , F z , M t (t), M b (t) and P are the time-average amount within a period: For the convenience of discussion and research, dimensioned quantities can be transformed into dimensionless coefficients: Propulsion efficiency is defined as the ratio of the thrust coefficient to power coefficient: 2.5.Solver Validation 2.5.1.Aerodynamic Force Static 3D aerodynamic force is verified by McAlister's experiment [33].The experiment is conducted in a wind tunnel to study a straight rectangular wing of airfoil NACA0015.The wing has an aspect ratio of 6.576 and a square wingtip.The experimental Reynolds number of the wing is 2,500,000.In the article, the data at an angle of attack of 12deg are used for comparative studies.A fixed wake model is used here, and the wake flows along the free flow.A total of 30 panels are used in the span direction and 90 panels are used in the chord direction.As shown in Figure 3, the pressure coefficient in most portions of the wing is in good agreement with the experiment.The boundary layer at a higher Reynolds number is relatively thin, so the pressure coefficient is closer to the inviscid case.There is an obvious error at the trailing edge near the wing tip.This may be due to the influence of the wing tip vortex on the wing tip boundary layer, resulting in an obvious three-dimensional boundary layer effect, and the boundary layer thickens at the wing tip trailing edge.The area of this effect is very small and does not affect the overall validation of the method.portions of the wing is in good agreement with the experiment.The boundary layer at a higher Reynolds number is relatively thin, so the pressure coefficient is closer to the inviscid case.There is an obvious error at the trailing edge near the wing tip.This may be due to the influence of the wing tip vortex on the wing tip boundary layer, resulting in an obvious three-dimensional boundary layer effect, and the boundary layer thickens at the wing tip trailing edge.The area of this effect is very small and does not affect the overall validation of the method.The verification of unsteady aerodynamic force is compared with Lin's calculation.Lin [34] studied the flapping motion of a rectangular wing with the airfoil of NACA0014 by solving the 3-dimensional Euler equations.The aspect ratio of the wing is 8, the flapping amplitude angle is 15 degrees, and the reduced frequency is 0.1.The wing has active The verification of unsteady aerodynamic force is compared with Lin's calculation.Lin [34] studied the flapping motion of a rectangular wing with the airfoil of NACA0014 by solving the 3-dimensional Euler equations.The aspect ratio of the wing is 8, the flapping amplitude angle is 15 degrees, and the reduced frequency is 0.1.The wing has active twisting that changes linearly along the span and the maximal twisting angle at the wing tip is 4 degrees.The phase shift between flapping and twisting is set as 90 degrees, which means that the pitching angle is maximal when the flapping angle is zero.
As shown in Figure 4, it can be seen that the results of the method in this article are in good agreement with the results of the unsteady 3-dimensional global method based on the Euler equations.
Aerospace 2023, 10, x FOR PEER REVIEW 9 of 25 twisting that changes linearly along the span and the maximal twisting angle at the wing tip is 4 degrees.The phase shift between flapping and twisting is set as 90 degrees, which means that the pitching angle is maximal when the flapping angle is zero.As shown in Figure 4, it can be seen that the results of the method in this article are in good agreement with the results of the unsteady 3-dimensional global method based on the Euler equations.

Geometric Nonlinearity
Dowell's experiment [35] is used to verify the geometric nonlinearity of the beam.In the experiment, a beam with a rectangular cross-section was used.The root was fixed on the device and could be rotated to different pitch angles.A concentrated downward load of P = 1 − 4lb is applied to the tip.Under different root pitch angles, the beam will produce displacements in the flatwise direction, edgewise direction and twist angle.The flatwise direction W and edgewise direction V are perpendicular to the long and short sides of the rectangular section, respectively.The twist angle is a characteristic of the nonlinear beam model.In the linear model, the twist angle would be zero due to the absence of the coupling between the twist and bending displacements.As shown in Figure 5, it can be seen that the calculated results are in good agreement with the experimental results.

Geometric Nonlinearity
Dowell's experiment [35] is used to verify the geometric nonlinearity of the beam.In the experiment, a beam with a rectangular cross-section was used.The root was fixed on the device and could be rotated to different pitch angles.A concentrated downward load of P = 1 − 4lb is applied to the tip.Under different root pitch angles, the beam will produce displacements in the flatwise direction, edgewise direction and twist angle.The flatwise direction W and edgewise direction V are perpendicular to the long and short sides of the rectangular section, respectively.The twist angle is a characteristic of the nonlinear beam model.In the linear model, the twist angle would be zero due to the absence of the coupling between the twist and bending displacements.As shown in Figure 5, it can be seen that the calculated results are in good agreement with the experimental results.

Fluid-Structure Coupling
In order to evaluate the accuracy and reliability of present method for calculating the aeroelastic response of flapping wings, the present method has been applied to the flapping wing studied by DeLaurier [2].DeLaurier designed and built a large unmanned flapping wing, and conducted calculations and experiments on the aircraft at an angle of attack of 6 degrees and a speed of 45 ft/s.The wings of DeLaurier's flapping wing are passively twisted and bended, which is similar to the flapping wing in this article.All aerodynamic and structural parameter information is provided in [2], where the model has been made according to DeLaurier's model and was calculated using the method of this article, and the results have been compared with DeLaurier's calculation and experiment.It takes about 15 minutes to calculate this example.
As shown in Figure 6a, the lift calculated by the present method is higher than the experimental result.The reason for this difference is that if the experiment was conducted at a Reynolds number of approximately 280,000, the lift will be reduced compared to the inviscid calculation results due to the influence of the boundary layer.As shown in Figure 6b, compared with the experimental results, the thrust calculated by the present method is bigger at high frequencies.The reason is that in high-frequency flapping motion, there will be a moment in the cycle when the angle of attack of the outer wing is relatively large, so airflow separation will occur and the drag of the airfoil will increase.The present method in this article is an inviscid method and cannot simulate this phenomenon.Since the calculation modeling is more accurate and complex, it has some advantages compared to the calculation results of DeLaurier.In order to evaluate the accuracy and reliability of present method for calculating the aeroelastic response of flapping wings, the present method has been applied to the flapping wing studied by DeLaurier [2].DeLaurier designed and built a large unmanned flapping wing, and conducted calculations and experiments on the aircraft at an angle of attack of 6 degrees and a speed of 45 ft/s.The wings of DeLaurier's flapping wing are passively twisted and bended, which is similar to the flapping wing in this article.All aerodynamic and structural parameter information is provided in [2], where the model has been made according to DeLaurier's model and was calculated using the method of this article, and the results have been compared with DeLaurier's calculation and experiment.It takes about 15 minutes to calculate this example.
As shown in Figure 6a, the lift calculated by the present method is higher than the experimental result.The reason for this difference is that if the experiment was conducted at a Reynolds number of approximately 280,000, the lift will be reduced compared to the inviscid calculation results due to the influence of the boundary layer.As shown in Figure 6b, compared with the experimental results, the thrust calculated by the present method is bigger at high frequencies.The reason is that in high-frequency flapping motion, there will be a moment in the cycle when the angle of attack of the outer wing is relatively large, so airflow separation will occur and the drag of the airfoil will increase.The present method in this article is an inviscid method and cannot simulate this phenomenon.Since the calculation modeling is more accurate and complex, it has some advantages compared to the calculation results of DeLaurier.
For the spanwise dynamic response distributions, as shown in Figure 6c, the twist angle amplitudes Θ of the inner wings are almost the same, but there are some differences in the outer wing sections.Their phase angles are slightly larger than the calculation results of DeLaurier.As shown in Figure 6d, the flapping elastic amplitudes H along the spanwise are flatter than the calculation results of DeLaurier, and their phase angles are slightly smaller.These deviations may be due to a combination of differences in aerodynamic and structural models.

Computational Model
In this article, the passive torsion of the wing is driven by the aerodynamic force and inertia force.Therefore, in order to utilize the aerodynamic force to produce torque on the beam as much as possible, the beam is placed at 10% of the chord line.In order to allow the wings to be passively twisted by aerodynamic force, the flapping wing in this article For the spanwise dynamic response distributions, as shown in Figure 6c, the twist angle amplitudes Θ max of the inner wings are almost the same, but there are some differences in the outer wing sections.Their phase angles are slightly larger than the calculation results of DeLaurier.As shown in Figure 6d, the flapping elastic amplitudes H max along the spanwise are flatter than the calculation results of DeLaurier, and their phase angles are slightly smaller.These deviations may be due to a combination of differences in aerodynamic and structural models.

Computational Model
In this article, the passive torsion of the wing is driven by the aerodynamic force and inertia force.Therefore, in order to utilize the aerodynamic force to produce torque on the beam as much as possible, the beam is placed at 10% of the chord line.In order to allow the wings to be passively twisted by aerodynamic force, the flapping wing in this article does not apply a structure similar to a D-box to prevent excessive torsional stiffness.It is assumed that the torsional stiffness is only provided by the beam.The wing profile is assumed to be rigid and is perpendicular to the profile of the beam element node.As shown in Figure 7a, the distribution of the beam elements is consistent with the aerodynamic panel, and is divided into N parts along the span direction by N-lines.The mass of the wing is assumed to be distributed by volume, discretely divided into mass points on the chord line, located at 10% of the chord behind the beam.Figure 7b shows the forces in the 2D wing profile.The airfoil performs plunge motion and pitching motion around the elastic axis, with aerodynamic and inertial forces as driving forces.The angle of attack is the angle of incidence plus the effects of wake-induced velocity and self-motion.

Structural and Kinematic Parameters
The flapping wing is modeled after the Canadian goose, weighing 3 kg and having a wing mass of 0.15 kg.The half wing span is 0.8m, the wing aspect ratio is 12, the wing root chord length is 0.3 m and the root-to-tip ratio is 2. The forward flight speed is 15 m/s, the range of wing Reynolds number is 367,500 to 183,750.The flapping period is 0.33 s, so the corresponding wingtip reduction frequency k is 0.0952 and the corresponding wingtip Strouhal number St is 0.167.The flapping motion is a simple harmonic motion, and the flapping angle is () = 30().As shown in Figure 6a, the wing root airfoil type is NACA8412, the wing tip airfoil type is NACA0012, and the rest of the wing has an interpolated transition from the wing root to the wing tip.The angle of incidence of the root airfoil is 7 degrees assuming the vehicle is performing steady level flight.If the flapping wing is designed to generate sufficient lift and thrust and have high propulsion efficiency, the airfoil must operate within a suitable angle of attack [20,36], so the wing needs to have a suitable deformation law during flapping.In order to balance manufacturing difficulty and flight efficiency, after extensive calculations, it is assumed that the cross-sectional diameter of the beam changes linearly with the wing span, which means the torsional stiffness changes with the wing span to the fourth power.As shown in Figure 8a, torsional stiffness GJ1 to GJ3 represent different torsional stiffnesses.As shown in Figure 8b, bend-

Structural and Kinematic Parameters
The flapping wing is modeled after the Canadian goose, weighing 3 kg and having a wing mass of 0.15 kg.The half wing span is 0.8m, the wing aspect ratio is 12, the wing root chord length is 0.3 m and the root-to-tip ratio is 2. The forward flight speed is 15 m/s, the range of wing Reynolds number is 367,500 to 183,750.The flapping period is 0.33 s, so the corresponding wingtip reduction frequency k is 0.0952 and the corresponding wingtip Strouhal number St is 0.167.The flapping motion is a simple harmonic motion, and the flapping angle is β(t) = 30cos(ωt).As shown in Figure 6a, the wing root airfoil type is NACA8412, the wing tip airfoil type is NACA0012, and the rest of the wing has an interpolated transition from the wing root to the wing tip.The angle of incidence of the root airfoil is 7 degrees assuming the vehicle is performing steady level flight.If the flapping wing is designed to generate sufficient lift and thrust and have high propulsion efficiency, the airfoil must operate within a suitable angle of attack [20,36], so the wing needs to have a suitable deformation law during flapping.In order to balance manufacturing difficulty and flight efficiency, after extensive calculations, it is assumed that the cross-sectional diameter of the beam changes linearly with the wing span, which means the torsional stiffness changes with the wing span to the fourth power.As shown in Figure 8a, torsional stiffness GJ1 to GJ3 represent different torsional stiffnesses.As shown in Figure 8b, bending stiffness EI1 to EI2 and the case of spanwise rigid represent different bending stiffnesses.Torsional and bending stiffnesses that are too small are of little practical significance and are not discussed.

Variation of Parameters with Initial Geometric Twist Angle
The twisting and bending moments of the wing are smallest at the wing tip and gradually increase toward the wing root.For passively deformed wings, the initial geometric twist angle of the wing will significantly affect wing deformation rules and aerodynamic performance parameters.The initial geometric twist angle is the twist change in the airfoil angle of incidence from the wing tip to the wing root during initial manufacturing, and the initial geometric twist angle varies linearly along the wing span.As shown in Figure 8, the aerodynamic coefficients of wings with different torsional and bending stiffnesses as the initial geometric twist angle changes are given.The biggest problem with the inviscid model is that it cannot simulate the trailing edge separation that will lead to a rapid increase in drag at large angles of attack.But this problem can be discussed through the pressure coefficient on the wing surface.Let the torsional stiffnesses of the wing be GJ1, GJ2, and GJ3, respectively, and the bending stiffnesses be rigid, EI1, and EI2, respectively.Combinations of different torsional and bending stiffnesses are presented in Case 1 to Case 9, as shown in Table 1.As shown in Figure 9a,b, the lift coefficient and thrust coefficient increase as the torsional stiffness of the beam increases.This is because the passive twist angle of the wing will decrease when the torsional stiffness increases, the range of angle of attack of the wing will increase, and therefore the thrust will increase.And due to the influence of aeroelasticity and the camber of the airfoil, the twist angle decreases while the angle of attack increases during downward flapping more than during upward flapping, so the lift increases.Under the same torsional stiffness, the lift will increase as the initial geometric twist angle increases, because it is equivalent to increasing the average angle of attack of the wing.Under the same torsional stiffness, as the initial geometric twist angle increases, the thrust coefficient first increases and then decreases.This is because thrust is mainly

Variation of Parameters with Initial Geometric Twist Angle
The twisting and bending moments of the wing are smallest at the wing tip and gradually increase toward the wing root.For passively deformed wings, the initial geometric twist angle of the wing will significantly affect wing deformation rules and aerodynamic performance parameters.The initial geometric twist angle is the twist change in the airfoil angle of incidence from the wing tip to the wing root during initial manufacturing, and the initial geometric twist angle varies linearly along the wing span.As shown in Figure 8, the aerodynamic coefficients of wings with different torsional and bending stiffnesses as the initial geometric twist angle changes are given.The biggest problem with the inviscid model is that it cannot simulate the trailing edge separation that will lead to a rapid increase in drag at large angles of attack.But this problem can be discussed through the pressure coefficient on the wing surface.Let the torsional stiffnesses of the wing be GJ1, GJ2, and GJ3, respectively, and the bending stiffnesses be rigid, EI1, and EI2, respectively.Combinations of different torsional and bending stiffnesses are presented in Case 1 to Case 9, as shown in Table 1.As shown in Figure 9a,b, the lift coefficient and thrust coefficient increase as the torsional stiffness of the beam increases.This is because the passive twist angle of the wing will decrease when the torsional stiffness increases, the range of angle of attack of the wing will increase, and therefore the thrust will increase.And due to the influence of aeroelasticity and the camber of the airfoil, the twist angle decreases while the angle of attack increases during downward flapping more than during upward flapping, so the lift increases.Under the same torsional stiffness, the lift will increase as the initial geometric twist angle increases, because it is equivalent to increasing the average angle of attack of the wing.Under the same torsional stiffness, as the initial geometric twist angle increases, the thrust coefficient first increases and then decreases.This is because thrust is mainly generated by the outer section of the wing, and the ability of the airfoil to generate thrust is strongest when the average angle of attack is near 0 degrees.gles of attack during the flapping process, and the lift-to-drag ratio of the airfoil at large angles of attack decreases significantly.Under the same torsional stiffness, as the initial geometric twist angle increases, the propulsion efficiency first increases and then decreases.This occurs because an excessively positive initial geometric twist angle leads to an overly high positive angle of attack in the outer wing section during downward flapping, while an excessively negative initial geometric twist angle results in an overly large negative angle of attack in the outer wing section during upward flapping.Excessively large positive and negative angles of attack will significantly increase drag and reduce thrust.Especially, the large negative angle of attack state should be avoided, because the airfoil is designed to work at a positive angle of attack, and the lift-to-drag ratio is very bad at the negative angle of attack state.In actual conditions, there will be a large range of airflow separation at the trailing edge at a large angle of attack, and the performance will be worse than that calculated using the inviscid method.It can be seen that under the same torsional stiffness, the lift coefficient, thrust coefficient and propulsion efficiency of a wing with smaller bending stiffness are stronger than those of a wing with rigid bending stiffness.The variation rules of different bending stiffnesses with initial geometric twist angle are basically the same.Bending stiffness has little impact on lift coefficient, but has a large impact on thrust coefficient and propulsion efficiency.The speed generated by the bending deformation and flapping are superimposed, As shown in Figure 9c, as the torsional stiffness increases, the propulsion efficiency decreases.This is because the increase in stiffness causes larger positive and negative angles of attack during the flapping process, and the lift-to-drag ratio of the airfoil at large angles of attack decreases significantly.Under the same torsional stiffness, as the initial geometric twist angle increases, the propulsion efficiency first increases and then decreases.This occurs because an excessively positive initial geometric twist angle leads to an overly high positive angle of attack in the outer wing section during downward flapping, while an excessively negative initial geometric twist angle results in an overly large negative angle of attack in the outer wing section during upward flapping.Excessively large positive and negative angles of attack will significantly increase drag and reduce thrust.Especially, the large negative angle of attack state should be avoided, because the airfoil is designed to work at a positive angle of attack, and the lift-to-drag ratio is very bad at the negative angle of attack state.In actual conditions, there will be a large range of airflow separation at the trailing edge at a large angle of attack, and the performance will be worse than that calculated using the inviscid method.
It can be seen that under the same torsional stiffness, the lift coefficient, thrust coefficient and propulsion efficiency of a wing with smaller bending stiffness are stronger than those of a wing with rigid bending stiffness.The variation rules of different bending stiffnesses with initial geometric twist angle are basically the same.Bending stiffness has little impact on lift coefficient, but has a large impact on thrust coefficient and propulsion efficiency.The speed generated by the bending deformation and flapping are superimposed, so that the angle of attack range of the wing becomes larger, and the result is an increase in thrust.The speed generated by the bending deformation of the wing is greater at the outer wing section, resulting in an increase in the angle of attack range of the outer wing, and the thrust weight moves to the outer wing.Since the outer wing's small camber airfoil generates thrust more efficiently, the total propulsion efficiency increases.
For a flapping wing aircraft that adopts passive twisting wings, it can achieve an effect similar to active twisting.Of course, the wing torsional stiffness and initial geometric twist angle must match to work efficiently.The initial geometric twist angle of the wing will also be different depending on the different torsional stiffness.It is best to make the outer wing section work in a larger angle of attack range where separation is not serious, so that sufficient thrust and propulsion efficiency can be obtained.When designing the torsional stiffness of the wing, we must consider the design strategy and purpose, whether it emphasizes large thrust or more efficiency.

Variation of Parameters during the Period
In order to further study the state of the wing and the effects of torsional and bending stiffness during the entire period, Case 1, Case 3, Case 6 and Case 9 were selected for comparison.Among them, Case 1 and Case 3 have the same torsional stiffness but different bending stiffnesses.Case 3, Case 6 and Case 9 have the same bending stiffness but different torsional stiffnesses.The initial geometric twist angle is taken as 0 degrees.
As shown in Figure 10, the changes in lift coefficient and thrust coefficient within a cycle are compared.It can be seen that the lift coefficient changes with time basically as a sinusoidal curve, while the flapping angle β(t) is a cosine curve.Due to the influence of unsteadiness and inertial force of the wing mass, the lift coefficient curve has a phase lag ϕ relative to the sinusoidal curve.The lift coefficient is positive in the downward flapping process and negative in a small section in the upward flapping process.The thrust coefficient is positive in the downward phase and is basically positive in the upward process, and the thrust coefficient in the downward process is greater than that in the upward phase.Comparing Case 1 with Case 3, as the bending stiffness decreases, the peak and average values of the lift coefficient and thrust coefficient will increase during the entire cycle, and the phase lag ϕ 3 is greater than ϕ 1 .Comparing Case 3, Case 6 and Case 9, it can be seen that as the torsional stiffness increases, the lift and thrust coefficient will increase, and the difference in phase lag is very small.
As shown in Figure 11a,c, the wingtip twist angles and wing root torsional moment coefficients during the period are basically sinusoidal curves.Since the lift center and the mass center are behind the torsion center, the twist angle is negative during the downward flapping process and positive during the upward phase.Because of the angle of incidence and camber of the wing, the positive lift is much greater than the negative lift during the entire cycle, so the negative twist angles are much larger than the positive twist angles and the torsional moment coefficients of wing root are negative during the period.Figure 11b,d show the maximum of the twist angles and torsional moment coefficient along the wing span.It can be seen that the twist angles are basically linear twist, with the largest twist angles at the wing tip and largest torsional moment coefficients at the wing root.By comparing Case 1 and Case 3 in Figure 10, as the bending stiffness decreases, the range of twist angle and torsional moment coefficient increases.This occurs because the combination of bending deformation and flapping motion produces greater lift and torsional moment, and phase lag occurs due to the combined motion.Comparing Case 3, Case 6 and Case 9, it shows that, as the torsional stiffness decreases, the twist angles increase and the torsional moment coefficients decrease.
process and negative in a small section in the upward flapping process.The thrust coefficient is positive in the downward phase and is basically positive in the upward process, and the thrust coefficient in the downward process is greater than that in the upward phase.Comparing Case 1 with Case 3, as the bending stiffness decreases, the peak and average values of the lift coefficient and thrust coefficient will increase during the entire cycle, and the phase lag  is greater than  .Comparing Case 3, Case 6 and Case 9, it can be seen that as the torsional stiffness increases, the lift and thrust coefficient will increase, and the difference in phase lag is very small.As shown in Figure 11a,c, the wingtip twist angles and wing root torsional moment coefficients during the period are basically sinusoidal curves.Since the lift center and the mass center are behind the torsion center, the twist angle is negative during the downward flapping process and positive during the upward phase.Because of the angle of incidence and camber of the wing, the positive lift is much greater than the negative lift during the entire cycle, so the negative twist angles are much larger than the positive twist angles and the torsional moment coefficients of wing root are negative during the period.Figure 11b,d show the maximum of the twist angles and torsional moment coefficient along the wing span.It can be seen that the twist angles are basically linear twist, with the largest twist angles at the wing tip and largest torsional moment coefficients at the wing root.By comparing Case 1 and Case 3 in Figure 10, as the bending stiffness decreases, the range of twist angle and torsional moment coefficient increases.This occurs because the combination of bending deformation and flapping motion produces greater lift and torsional moment, and phase lag occurs due to the combined motion.Comparing Case 3, Case 6 and Case 9, it shows that, as the torsional stiffness decreases, the twist angles increase and the torsional moment coefficients decrease.As shown in Figure 12, since the wing in Case 1 is spanwise rigid, there is no bending deformation.Comparing Case 3, Case 6 and Case 9 in Figure 12a,c, the wing with higher torsional stiffness has a larger bending deformation and bending moment coefficient at the wing root.This occurs because the large torsional stiffness makes the twist angle range smaller and the angle of attack range and lift coefficient range larger, resulting in a larger bending moment coefficient and bending deformation.As shown in Figure 12b, the maximum deformation of the wing tip in the vertical direction is more than 0.2 times the wing span, and the displacement occurs in the span direction, reflecting the characteristic of geometric nonlinearity.As shown in Figure 12d, the bending moment coefficient gradually increases from the wing tip to the wing root.As shown in Figure 12, since the wing in Case 1 is spanwise rigid, there is no bending deformation.Comparing Case 3, Case 6 and Case 9 in Figure 12a,c, the wing with higher torsional stiffness has a larger bending deformation and bending moment coefficient at the wing root.This occurs because the large torsional stiffness makes the twist angle range smaller and the angle of attack range and lift coefficient range larger, resulting in a larger bending moment coefficient and bending deformation.As shown in Figure 12b, the maximum deformation of the wing tip in the vertical direction is more than 0.2 times the wing span, and the displacement occurs in the span direction, reflecting the characteristic of geometric nonlinearity.As shown in Figure 12d, the bending moment coefficient gradually increases from the wing tip to the wing root.Figure 13 shows the average lift coefficient and average thrust coefficient along the span direction during the period.The lift coefficient gradually decreases from the wing root to the wing tip.The thrust coefficient first gradually increases from the wing root, and due to the three-dimensional effect, the thrust coefficient of the wing tip will decrease.At the wing root, the thrust is even negative.Generally speaking, the outer wing mainly produces thrust, and the inner wing mainly produces lift.Comparing Case 1 and Case 3, it can be seen that as the bending stiffness decreases, the effect on the lift coefficient is not as obvious as on the thrust coefficient.The maximum value of the thrust coefficient will become larger and the center of thrust generation will be closer to the wing tip.Since the outer wing's small camber airfoil is more suitable for generating thrust efficiently, the propulsion efficiency will also increase.Comparing Case 3, Case 6 and Case 9, changes in torsional stiffness have no obvious impact on the thrust coefficient, but the lift coefficient of a wing with large torsional stiffness will increase significantly throughout the span.Figure 13 shows the average lift coefficient and average thrust coefficient along the span direction during the period.The lift coefficient gradually decreases from the wing root to the wing tip.The thrust coefficient first gradually increases from the wing root, and due to the three-dimensional effect, the thrust coefficient of the wing tip will decrease.At the wing root, the thrust is even negative.Generally speaking, the outer wing mainly produces thrust, and the inner wing mainly produces lift.Comparing Case 1 and Case 3, it can be seen that as the bending stiffness decreases, the effect on the lift coefficient is not as obvious as on the thrust coefficient.The maximum value of the thrust coefficient will become larger and the center of thrust generation will be closer to the wing tip.Since the outer wing's small camber airfoil is more suitable for generating thrust efficiently, the propulsion efficiency will also increase.Comparing Case 3, Case 6 and Case 9, changes in torsional stiffness have no obvious impact on the thrust coefficient, but the lift coefficient of a wing with large torsional stiffness will increase significantly throughout the span.

Pressure Coefficient at Different Portions during the Period
In order to further study the state of the airfoil surface during the period, the airfoil surface pressure coefficients of different span positions at a specific time in Case 1, Case 3, Case 6, and Case 9 were compared.The time selected should include the moment of maximum lift coefficient so that more features will be revealed.As shown in Figure 14, in Case 1, the pressure coefficients at 25% and 75% the wing span for t/T = 0.25 +  /2, t/T = 0.5 +  /2, t/T = 0.75 +  /2, and t/T = 1 +  /2 are selected.It can be seen that at 25% span of the wing, the airfoil camber and average angle of attack are large, and the unsteadiness is weak.During the entire period, the airfoil basically operates within the positive angle of attack range and mainly generates lift.At 75% span, the airfoil is at a large positive angle of attack when t/T = 0.25 +  /2, and when t/T = 0.75 +  /2, the airfoil is at a large negative angle of attack.At these moments, a sharp suction peak is generated at the leading edge of the airfoil, and the airfoil has the strongest ability to generate thrust.Due to the camber of the airfoil, the pressure coefficient curve is narrower at negative angles of attack, indicating that the lift-

Pressure Coefficient at Different Portions during the Period
In order to further study the state of the airfoil surface during the period, the airfoil surface pressure coefficients of different span positions at a specific time in Case 1, Case 3, Case 6, and Case 9 were compared.The time selected should include the moment of maximum lift coefficient so that more features will be revealed.As shown in Figure 14, in Case 1, the pressure coefficients at 25% and 75% the wing span for t/T = 0.25 + ϕ 1 /2π, t/T = 0.5 + ϕ 1 /2π, t/T = 0.75 + ϕ 1 /2π, and t/T = 1 + ϕ 1 /2π are selected.It can be seen that at 25% span of the wing, the airfoil camber and average angle of attack are large, and the unsteadiness is weak.During the entire period, the airfoil basically operates within the positive angle of attack range and mainly generates lift.

Pressure Coefficient at Different Portions during the Period
In order to further study the state of the airfoil surface during the period, the airfoil surface pressure coefficients of different span positions at a specific time in Case 1, Case 3, Case 6, and Case 9 were compared.The time selected should include the moment of maximum lift coefficient so that more features will be revealed.As shown in Figure 14, in Case 1, the pressure coefficients at 25% and 75% the wing span for t/T = 0.25 +  /2, t/T = 0.5 +  /2, t/T = 0.75 +  /2, and t/T = 1 +  /2 are selected.It can be seen that at 25% span of the wing, the airfoil camber and average angle of attack are large, and the unsteadiness is weak.During the entire period, the airfoil basically operates within the positive angle of attack range and mainly generates lift.At 75% span, the airfoil is at a large positive angle of attack when t/T = 0.25 +  /2, and when t/T = 0.75 +  /2, the airfoil is at a large negative angle of attack.At these moments, a sharp suction peak is generated at the leading edge of the airfoil, and the airfoil has the strongest ability to generate thrust.Due to the camber of the airfoil, the pressure coefficient curve is narrower at negative angles of attack, indicating that the lift- At 75% span, the airfoil is at a large positive angle of attack when t/T = 0.25 + ϕ 1 /2π, and when t/T = 0.75 + ϕ 1 /2π, the airfoil is at a large negative angle of attack.At these moments, a sharp suction peak is generated at the leading edge of the airfoil, and the airfoil has the strongest ability to generate thrust.Due to the camber of the airfoil, the pressure coefficient curve is narrower at negative angles of attack, indicating that the lift-to-drag ratio of the airfoil is not high at that moment.When t/T = 0.25 + ϕ 1 /2π and t/T = 0.75 + ϕ 1 /2π, the maximum absolute value of the airfoil surface pressure coefficient is similar, so that the airfoil can thrust efficiently; otherwise, the excessive angle of attack during the upward or downward process may cause the lift-to-drag ratio of the airfoil to drop sharply.When t/T = 0.5 + ϕ 1 /2π and t/T = 1 + ϕ 1 /2π, the airfoil is in a small angle of attack range and basically only generates drag.Since it is an inviscid model, the phenomenon of laminar flow transition and air flow separation at the trailing edge cannot be seen.
To study the effect of wing bending stiffness, the pressure coefficients in Case 3 are compared with those in Case 1.In Case 3, the pressure coefficients at 25% and 75% of the wing span for t/T = 0.25 + ϕ 3 /2π, t/T = 0.5 + ϕ 3 /2π, t/T = 0.75 + ϕ 3 /2π, and t/T = 1 + ϕ 3 /2π are selected.As shown in Figure 15, at 25% wing span, the pressure coefficient in Case 3 is almost the same as that in Case 1, which can explain why the lift coefficient does not change much as the bending stiffness changes.At 75% span, the peak value of the pressure coefficient in Case 3 is greater than that in Case 1 as the wing is spanwise rigid.This shows that the wing with appropriate spanwise flexibility can provide a larger thrust coefficient at the outer wing.
Aerospace 2023, 10, x FOR PEER REVIEW 20 of 25 to-drag ratio of the airfoil is not high at that moment.When t/T = 0.25 +  /2 and t/T = 0.75 +  /2, the maximum absolute value of the airfoil surface pressure coefficient is similar, so that the airfoil can thrust efficiently; otherwise, the excessive angle of attack during the upward or downward process may cause the lift-to-drag ratio of the airfoil to drop sharply.When t/T = 0.5 +  /2 and t/T = 1 +  /2, the airfoil is in a small angle of attack range and basically only generates drag.Since it is an inviscid model, the phenomenon of laminar flow transition and air flow separation at the trailing edge cannot be seen.
To study the effect of wing bending stiffness, the pressure coefficients in Case 3 are compared with those in Case 1.In Case 3, the pressure coefficients at 25% and 75% of the wing span for t/T = 0.25 +  /2, t/T = 0.5 +  /2, t/T = 0.75 +  /2, and t/T = 1 +  /2 are selected.As shown in Figure 15, at 25% wing span, the pressure coefficient in Case 3 is almost the same as that in Case 1, which can explain why the lift coefficient does not change much as the bending stiffness changes.At 75% span, the peak value of the pressure coefficient in Case 3 is greater than that in Case 1 as the wing is spanwise rigid.This shows that the wing with appropriate spanwise flexibility can provide a larger thrust coefficient at the outer wing.By comparing Case 1 and Case 3, it can be seen that reducing the bending stiffness has a greater impact on the outer wing section, which increases its angle of attack range and increases its thrust coefficient.It should be noted that this conclusion is drawn under inviscid conditions.In fact, when the airfoil is at a larger angle of attack, the airflow will separate more strongly at the trailing edge, which will deteriorate the performance.The performance may not be so optimistic.
To study the effect of wing torsional stiffness, the pressure coefficients in Case 6 and Case 9 are used to compare with that in Case 3. In Case 6, the pressure coefficients at 25% and 75% wing span for t/T = 0.25 +  /2, t/T = 0.5 +  /2, t/T = 0.75 +  /2 and t/T = 1 +  /2 are selected.As shown in Figure 16, at 25% of the wing span, the pressure coefficient shows that Case 6 has a larger positive angle of attack than Case 3, indicating that as the torsional stiffness increases, the lift will increase.At 75% wing span, the pressure coefficient leading edge suction peak at t/T = 0.25 +  /2 and t/T = 0.75 +  /2 in Case 6 is larger than that in Case 3, indicating that increasing torsional stiffness can provide a larger thrust coefficient at the outer wing.Moreover, the maximum absolute value of the pressure coefficient in these two moments is quite different, which means that the angle of attack is too large during the downward flapping process.This shows that the By comparing Case 1 and Case 3, it can be seen that reducing the bending stiffness has a greater impact on the outer wing section, which increases its angle of attack range and increases its thrust coefficient.It should be noted that this conclusion is drawn under inviscid conditions.In fact, when the airfoil is at a larger angle of attack, the airflow will separate more strongly at the trailing edge, which will deteriorate the performance.The performance may not be so optimistic.
To study the effect of wing torsional stiffness, the pressure coefficients in Case 6 and Case 9 are used to compare with that in Case 3. In Case 6, the pressure coefficients at 25% and 75% wing span for t/T = 0.25 + ϕ 6 /2π, t/T = 0.5 + ϕ 6 /2π, t/T = 0.75 + ϕ 6 /2π and t/T = 1 + ϕ 6 /2π are selected.As shown in Figure 16, at 25% of the wing span, the pressure coefficient shows that Case 6 has a larger positive angle of attack than Case 3, indicating that as the torsional stiffness increases, the lift will increase.At 75% wing span, the pressure coefficient leading edge suction peak at t/T = 0.25 + ϕ 6 /2π and t/T = 0.75 + ϕ 6 /2π in Case 6 is larger than that in Case 3, indicating that increasing torsional stiffness can provide a larger thrust coefficient at the outer wing.Moreover, the maximum absolute value of the pressure coefficient in these two moments is quite different, which means that the angle of attack is too large during the downward flapping process.This shows that the initial geometric twist angle of the wing is not appropriate under this torsional stiffness, and the initial geometric twist angle needs to be reduced to prevent an excessively large angle of attack to improve propulsion efficiency.
initial geometric twist angle of the wing is not appropriate under this torsional stiffness, and the initial geometric twist angle needs to be reduced to prevent an excessively large angle of attack to improve propulsion efficiency.In Case 9, the pressure coefficients at 25% and 75% of the wing span for t/T = 0.25 +  /2, t/T = 0.5 +  /2, t/T = 0.75 +  /2, and t/T = 1 +  /2 are selected.As shown in Figure 17, at 25% of the wing span, the pressure coefficients show that Case 9 has a smaller positive angle of attack than Case 3, indicating that as the torsional stiffness decreases, the lift will decrease.At 75% wing span, the pressure coefficient leading edge suction peak at t/T = 0.25 +  /2 and t/T = 0.75 +  /2 in Case 9 is smaller than that in Case 3, indicating that decreasing torsional stiffness provides a smaller thrust coefficient at the outer wing.It can be seen from the pressure coefficient that, compared with the downward flapping process, the angle of attack is excessively large during upward flapping.This shows that the initial geometric twist angle needs to be increased to prevent an excessively large angle of attack during the upward flapping process to improve propulsion efficiency.In Case 9, the pressure coefficients at 25% and 75% of the wing span for t/T = 0.25 + ϕ 9 /2π, t/T = 0.5 + ϕ 9 /2π, t/T = 0.75 + ϕ 9 /2π, and t/T = 1 + ϕ 9 /2π are selected.As shown in Figure 17, at 25% of the wing span, the pressure coefficients show that Case 9 has a smaller positive angle of attack than Case 3, indicating that as the torsional stiffness decreases, the lift will decrease.At 75% wing span, the pressure coefficient leading edge suction peak at t/T = 0.25 + ϕ 9 /2π and t/T = 0.75 + ϕ 9 /2π in Case 9 is smaller than that in Case 3, indicating that decreasing torsional stiffness provides a smaller thrust coefficient at the outer wing.It can be seen from the pressure coefficient that, compared with the downward flapping process, the angle of attack is excessively large during upward flapping.This shows that the initial geometric twist angle needs to be increased to prevent an excessively large angle of attack during the upward flapping process to improve propulsion efficiency.
initial geometric twist angle of the wing is not appropriate under this torsional stiffness, and the initial geometric twist angle needs to be reduced to prevent an excessively large angle of attack to improve propulsion efficiency.In Case 9, the pressure coefficients at 25% and 75% of the wing span for t/T = 0.25 +  /2, t/T = 0.5 +  /2, t/T = 0.75 +  /2, and t/T = 1 +  /2 are selected.As shown in Figure 17, at 25% of the wing span, the pressure coefficients show that Case 9 has a smaller positive angle of attack than Case 3, indicating that as the torsional stiffness decreases, the lift will decrease.At 75% wing span, the pressure coefficient leading edge suction peak at t/T = 0.25 +  /2 and t/T = 0.75 +  /2 in Case 9 is smaller than that in Case 3, indicating that decreasing torsional stiffness provides a smaller thrust coefficient at the outer wing.It can be seen from the pressure coefficient that, compared with the downward flapping process, the angle of attack is excessively large during upward flapping.This shows that the initial geometric twist angle needs to be increased to prevent an excessively large angle of attack during the upward flapping process to improve propulsion efficiency.By comparing the pressure coefficients of Cases 3, 6, and 9, it is shown that increasing torsional stiffness can increase thrust and lift.At the same time, the torsional stiffness of the wing needs to match the initial geometric twist angle to prevent an excessively large angle of attack during the cycle, thereby reducing the airfoil lift-to-drag ratio.

Conclusions
This article establishes a method to quickly and effectively solve the fluid-structure interaction problem of flapping wings in the early design stage.In computing such problems, a lot of time can be saved compared to a method that requires days to solve the N-S equations in the global domain.In the method, a three-dimensional inviscid unsteady aerodynamic model is used to calculate the pressure coefficient of the wing surface and aerodynamic performance, a linear finite element model combined with the corotational formulation method is used to solve the geometric nonlinear problem of the beam, and the Newmark method is used to solve the dynamic coupling problem.The article analyzes the variation of lift coefficient, thrust coefficient and propulsive efficiency with the initial geometric twist angle of the wing under different torsion and bending stiffness.And the lift coefficient, thrust coefficient, torsional and bending deformation of the wing during the entire cycle, as well as the pressure coefficient of the wing surface at different moment were analyzed under different torsional bending stiffness.The article reaches the following conclusions: 1.
By comparing with experiments and calculations, the method proposed in this paper has relatively high accuracy in calculating three-dimensional inviscid unsteady aerodynamic forces and three-dimensional slender geometric nonlinear beams, and is suitable for solving the problem in this article.

2.
For a passively twisted wing, an increase in torsional stiffness leads to an increase in the lift and thrust coefficients and a concomitant decrease in propulsive efficiency.This illustrates the need to make trade-offs during design.

3.
Under the same torsional stiffness, as the initial geometric twist angle of the wing increases, the lift coefficient increases, and the thrust coefficient and propulsion efficiency first increase and then decrease.For a certain wing torsional stiffness, there is a suitable initial geometric twist angle that maximizes propulsion efficiency.An inappropriate initial geometric twist angle will lead to an excessive angle of attack at a certain time during the period, thereby reducing propulsion efficiency.4.
For a passively twisted wing, the wing is in a bending-torsion coupling state.Compared with a spanwise rigid wing, a wing with some spanwise flexibility can increase the thrust coefficient of the outer wing, making the thrust distribution along the wing span more reasonable, thereby increasing the overall thrust coefficient and propulsion efficiency.

Figure 2 .
Figure 2. Beam elements in local and global coordinate systems.

Figure 2 .
Figure 2. Beam elements in local and global coordinate systems.

Figure 3 .
Figure 3.Comparison of calculation results with McAlister's experiment.

Figure 3 .
Figure 3.Comparison of calculation results with McAlister's experiment.

Figure 4 .
Figure 4. Comparison of calculation results with Lin's calculations.

Figure 4 .
Figure 4. Comparison of calculation results with Lin's calculations.

Figure 5 .
Figure 5.Comparison of calculation results with Dowell's experiments.
(a) Average lift (b) Average thrust (c) Torsion amplitude and phase angle (d) Flapping elastic amplitude and phase angle

Figure 6 .
Figure 6.Comparison of calculation results with Delaurier's calculation and experiment.

Figure 6 .
Figure 6.Comparison of calculation results with Delaurier's calculation and experiment.

Figure 9 .
Figure 9. Variation of parameters with initial geometric twist angle.

Figure 9 .
Figure 9. Variation of parameters with initial geometric twist angle.

Figure 10 .
Figure 10.Lift and thrust coefficients during the period.

Figure 10 .Figure 11 .
Figure 10.Lift and thrust coefficients during the period.

Figure 11 .
Figure 11.Wing twist angle and torsional moment coefficient.

Figure 12 .
Figure 12.Wing bending displacement and bending moment coefficient.

Figure 13 .
Figure 13.Average lift coefficient and thrust coefficient along the wing span during the period.
(a) Pressure coefficient at 25% the wing span (b) Pressure coefficient at 75% the wing span

Figure 14 .
Figure 14.Pressure coefficient at different portions at various times in Case 1.

Figure 13 .
Figure 13.Average lift coefficient and thrust coefficient along the wing span during the period.

Figure 13 .
Figure 13.Average lift coefficient and thrust coefficient along the wing span during the period.
(a) Pressure coefficient at 25% the wing span (b) Pressure coefficient at 75% the wing span

Figure 14 .
Figure 14.Pressure coefficient at different portions at various times in Case 1.

Figure 14 .
Figure 14.Pressure coefficient at different portions at various times in Case 1.
(a) Pressure coefficient at 25% the wing span (b) Pressure coefficient at 75% the wing span

Figure 15 .
Figure 15.Pressure coefficient at different portions at various times in Case 3.

Figure 15 .
Figure 15.Pressure coefficient at different portions at various times in Case 3.
(a) Pressure coefficient at 25% the wing span (b) Pressure coefficient at 75% wing span

Figure 16 .
Figure 16.Pressure coefficient at different portions at various times in Case 6.
(a) Pressure coefficient at 25% the wing span (b) Pressure coefficient at 25% the wing span

Figure 17 .
Figure 17.Pressure coefficient at different portions at various times in Case 9.

Figure 16 .
Figure 16.Pressure coefficient at different portions at various times in Case 6.
(a) Pressure coefficient at 25% the wing span (b) Pressure coefficient at 75% wing span

Figure 16 .
Figure 16.Pressure coefficient at different portions at various times in Case 6.
(a) Pressure coefficient at 25% the wing span (b) Pressure coefficient at 25% the wing span

Figure 17 .
Figure 17.Pressure coefficient at different portions at various times in Case 9.Figure 17.Pressure coefficient at different portions at various times in Case 9.

Figure 17 .
Figure 17.Pressure coefficient at different portions at various times in Case 9.Figure 17.Pressure coefficient at different portions at various times in Case 9.

Table 1 .
Combinations of different torsional and bending stiffnesses.

Table 1 .
Combinations of different torsional and bending stiffnesses.