3.1.1. Calculation of membrane structure parameters

In the process of bat flapping, the membrane will produce a certain degree of spatial deformation due to its movement and elasticity, which can be divided into two parts: The first part is that because the wing skeleton flaps around the flapping axis, but the tail remains stationary, there will be an angle between the two, but the membrane is stuck to the wing skeleton and tail at the same time, so there will be a height difference between the leading edge and the trailing edge of the membrane relative to the initial plane. The second part is the elastic deformation of the membrane because of the aerodynamic force. To facilitate understanding, the deformation of the membrane is decomposed for illustration, as shown in Fig. 7.

Figure 7. The decomposition of membrane deformation.

Figure 7(a) represents the three-dimensional model diagram of the membrane in the downward process, including the deformation of the membrane and the projection of the membrane on the initial plane (projection is represented by dotted lines). Figure 7(b) shows the model in the downward process of the membrane, which is used to explain the deformation of the membrane, including the deformation diagram of the membrane and the projection diagram of the initial plane; Fig. 7(c) and (d) are schematic diagrams of the deformation of the membrane in Fig. 7(b) after decomposition (separated along the plane P ₁ P ₂ P ₃ P ₄ P ₅ P ₆ P ₁₀), where Fig. 7(c) represents the elastic deformation of the membrane because of the aerodynamic force; Fig. 7(d) represents the height difference between the torsional deformation of the leading and trailing edges of the membrane relative to the initial plane (the plane composed of dotted lines).

In Fig. 7, $\varphi$ (r, t) is the angle of the strip (The strip with a distance r from the flapping axis P ₁ P ₁₀ at time t) relative to the initial plane, which is defined as the passive torsion angle. Its physical significance is the angle between the strip string and the flapping axis after the membrane twisted around the leading edge axis. h _Pr1 is the height of the leading edge due to flapping; h _Pr2 is the height of the trailing edge due to flapping; Δh _head is the height of the leading edge due to the flexible deformation of the membrane; Δh _back is the height of the trailing edge due to the flexible deformation of the membrane. The above deformation parameters ( $\varphi$ (r, t), h _Pr1, h _Pr2, Δh _head, Δh _back, etc.) all represent the deformation of the membrane, so these variables are located on the membrane plane and in the membrane coordinate system(O _m -x _m y _m z _m ) .

Due to the particularity of the bat wing structure, the deformation of the side membrane fixed to the tail, the arm membrane of the middle section, and the hand membrane of the wing hand part are different, so the research needs to be divided into three parts [Reference Joshi, Jaiman, Li, Breuer and Swartz22], as shown in Fig. 8.

Figure 8. The membrane is divided into three parts.

The first part (P ₁ P ₁₀ P ₆ P ₁₁) is the part connected with the membrane and the body, in which P ₁ P ₁₀ is connected with the body and P ₁₀ P ₆ is connected with the tail, both of which do not produce deformation, so this part only needs to calculate the deformation of the leading edge. The second part(P ₁₁ P ₆ P ₅ P ₂) is the middle part of the wing; in this part, two equivalent deformations should be considered and calculated simultaneously. The third part(P ₂ P ₅ P ₄ P ₃) is the part where the membrane is connected with the wing hand; this part only needs to calculate the deformation of the trailing edge.

First, for the first part, analytical diagram of bat membrane deformation is shown in Fig. 9. In the figure, r represents the distance between the strip and the flapping axis P ₁ P ₁₀; $\varphi$ (r, t) is the passive torsion angle of the strip at this moment where the distance from the flapping axis is r; r _P6 represents the distance between the strip P ₆ P ₁₁ and the flapping axis; $\varphi$ _P6(t) represents the passive torsion angle of the strip where P ₆ P ₁₁ is located, which is also the largest passive torsion angle (absolute value) of the first part; θ ₀ is the flapping angle; Δh _head is the height formed by the flexible deformation of the membrane from the leading edge of the strip at a distance r from the flapping axis.

Figure 9. Analytical diagram of bat membrane deformation.

P ₁ P ₂ is the leading edge of the membrane, which will be deformed by aerodynamic force. For calculation, P ₁ P ₂ is assumed to be a beam structure (a beam with a large amount of elastic deformation) with both ends fixed instantaneously. The beam is subjected to uniform force (assuming the aerodynamic force as uniform force), which will produce bending deformation, and its height due to elastic deformation is:

(2) \begin{equation}\Delta h_{\text{lead}}=\frac{\frac{q}{2}\left(\frac{1}{12}r^{4}-\frac{r_{P1P2}}{6}r^{3}\right)+\frac{q{r_{P1P2}}^{3}}{24}r}{\text{EI}_{\text{lead}}}\end{equation}

where EI _lead is the bending rigidity of the leading edge; q is the uniform aerodynamic load; r _P1P2 is the distance between P ₁ and P ₂.

Through the geometric relationship in Fig. 9, the expressions of chord length and passive torsion angle can be obtained:

(3) \begin{equation}c\left(r,t\right)=\sqrt{\left(2r\cdot \sin \frac{\theta _{o}}{2}-\text{sign}\left(\theta _{o}\frac{d\theta _{o}}{dt}\right)\Delta h_{\text{lead}}\right)^{2}+\left(P_{1}P_{10}+r\tan 30^{\circ}\right)^{2}}, 0\leq r\leq r_{P6}\left(t\right)\end{equation}

(4) \begin{equation}\varphi \left(r,t\right)=\pm \arctan \frac{2r\cdot \sin \frac{\theta _{o}}{2}-\text{sign}\left(\theta _{o}\frac{d\theta _{o}}{dt}\right)\Delta h_{\text{lead}}}{P_{1}P_{10}+r\tan 30^{\circ}}, 0\leq r\leq r_{P6}\left(t\right)\end{equation}

Next, calculate the second part, that is r _P6(t) ≤r ≤P ₀ P ₁. The schematic diagram of the structure of this part is shown in (c)and (d) in Fig. 7. Similar to the deformation of the leading edge, the trailing edge part P ₅ P ₆ is also assumed to be a beam structure with both ends fixed instantaneously, and then, the height of the leading and trailing edges due to elastic deformation are:

(5) \begin{equation}\Delta h_{\text{lead}}=\frac{\frac{q}{2}\left(\frac{1}{12}r^{4}-\frac{r_{P1P2}}{6}r^{3}\right)+\frac{q{r_{P1P2}}^{3}}{24}r}{\text{EI}_{\text{lead}}}\end{equation}

(6) \begin{equation}\Delta h_{\text{back}}=\frac{\frac{q}{2}\left(\frac{1}{12}\left(r-r_{P6}\right)^{4}-\frac{r_{P5P6}}{6}\left(r-r_{P6}\right)^{3}\right)+\frac{q{r_{P5P6}}^{3}}{24}\left(r-r_{P6}\right)}{\text{EI}_{\text{back}}}\end{equation}

where EI _back is the bending rigidity of the trailing edge; r _P5P6 is the distance between P ₅ and P ₆.; r _P6 is the distance between P ₆ and flapping axis.

The height of any point on the leading edge and the height of the corresponding trailing edge point:

(7) \begin{equation}h_{Pr1}=\frac{h_{P2}}{r_{P1P2}}r_{Pr1}=2r_{Pr1}\cdot \sin \frac{\theta _{o}}{2}\end{equation}

(8) \begin{equation}h_{Pr2}=\frac{r_{Pr2}-r_{P6}}{r_{P1P2}-r_{P6}}h_{P2}\end{equation}

The projection on horizontal plane of the chord length formed by the connecting point P _r1 and P _r2:

(9) \begin{equation}P'_{\!\!r1}P'_{\!\!r2}=P_{2}P_{5}+\left(\frac{r_{P1P2}-r_{Pr1}}{r_{P1P2}-r_{P6}}\right)\left[\left(r_{P6}\cdot \tan 30^{\circ}+P_{1}P_{10}\right)-P_{2}P_{5}\right]\end{equation}

The height difference between P _r1 and P _r2 is:

(10) \begin{equation}h_{Pr2Pr1}=\left(h_{Pr2}-\text{sign}\left(\theta _{o}\frac{d\theta _{o}}{dt}\right)\Delta h_{\text{back}}\right)-\left(h_{Pr1}-\text{sign}\left(\theta _{o}\frac{d\theta _{o}}{dt}\right)\Delta h_{\text{lead}}\right)\end{equation}

So, the expressions of chord length and passive torsion angle can be obtained:

(11) \begin{equation}\varphi \left(r,t\right)=\mathrm{atan} \left(\frac{h_{Pr2Pr1}}{P'_{\!\!r1}P'_{\!\!r2}}\right),r_{P6}\leq r\leq r_{P1P2}\end{equation}

(12) \begin{equation}c\left(r,t\right)=P_{r1}P_{r2}=\sqrt{{h_{Pr2Pr1}}^{2}+\left(P'_{\!\!r1}P'_{\!\!r2}\right)^{2}},r_{P6}\leq r\leq r_{P1P2}\end{equation}

Finally, for the third part, that is the torsion angle of the membrane between r≥ P ₁ P ₂, because this part is connected with the wing hand part, and the skeleton layout of the wing hand is fixed, the chord length of this part is related to the shape of the wingtip

For this part, the height of the trailing edge due to elastic deformation is:

(13) \begin{equation}\Delta h_{\text{hand}}=\frac{\frac{q}{2}\left(\frac{1}{12}\left(\frac{r-r_{P1P2}}{\cos \varphi _{\text{hand}}}\right)^{4}-\frac{1}{6}\frac{r_{P2P3}}{\cos \varphi _{\text{hand}}}\left(\frac{r-r_{P1P2}}{\cos \varphi _{\text{hand}}}\right)^{3}\right)+\frac{q}{24}\left(\frac{r_{P2P3}}{\cos \varphi _{\text{hand}}}\right)^{3}\left(\frac{r-r_{P1P2}}{\cos \varphi _{\text{hand}}}\right)}{\text{EI}_{\text{hand}}}\end{equation}

where

(14) \begin{equation}\cos \varphi _{\text{hand}}=\frac{r_{P2P3}}{\sqrt{\left(r_{P2P3}\right)^{2}+\left(r_{P2P5}\right)^{2}}}\end{equation}

So, the expressions of chord length and passive torsion angle can be obtained:

(15) \begin{equation}c\left(r,t\right)=\sqrt{\left(P_{2}P_{5}\cdot \left(1-\left(\frac{r-r_{P1P2}}{r_{\max }-r_{P1P2}}\right)^{I}\right)\right)^{2}+\left(\Delta h_{\text{hand}}\right)^{2}}, r\geq r_{P1P2}\end{equation}

where I is the wingtip shape index (I reflects the wingtip shape).

(16) \begin{equation}\varphi \left(r,t\right)=\mathrm{atan} \left(\frac{\Delta h_{\text{hand}}}{P_{2}P_{5}\cdot \left(1-\left(\frac{r-r_{P1P2}}{r_{P2P3}}\right)^{I}\right)}\right), r\geq r_{P1P2}\end{equation}

In summary, the passive torsion angle of the wing membrane is:

(17) \begin{equation}\varphi \left(r,t\right)=\begin{cases} \pm \arctan \frac{2r\cdot \sin \frac{\theta _{o}}{2}-\text{sign}\left(\theta _{o}\frac{d\theta _{o}}{dt}\right)\Delta h_{\text{lead}}}{P_{1}P_{10}+r\tan 30^{\circ}}, 0\leq r\leq r_{P6}\left(t\right)\\[9pt] \text{atan}\left(\frac{h_{Pr2Pr1}}{P'_{\!\!r1}P'_{\!\!r2}}\right),r_{P6}\leq r\leq r_{P1P2}\\[9pt] \text{atan}\left(\frac{\Delta h_{\text{hand}}}{r_{P2P5}\cdot \left(1-\left(\frac{r-r_{P1P2}}{r_{P2P3}}\right)^{I}\right)}\right), r\geq r_{P1P2} \end{cases}\end{equation}

The chord length distribution of the wing surface c(r, t) can be summarized as:

(18) \begin{equation}c\left(r,t\right)=\begin{cases} \sqrt{\left(2r\cdot \sin \frac{\theta _{o}}{2}-\text{sign}\left(\theta _{o}\frac{d\theta _{o}}{dt}\right)\Delta h_{\text{lead}}\right)^{2}+\left(P_{1}P_{10}+r\tan 30^{\circ}\right)^{2}}, 0\leq r\leq r_{P6}\left(t\right)\\[5pt] \sqrt{{h_{Pr2Pr1}}^{2}+\left(P'_{\!\!r1}P'_{\!\!r2}\right)^{2}},r_{P6}\leq r\leq r_{P1P2}\\[5pt] \sqrt{\left(r_{P2P5}\cdot \left(1-\left(\frac{r-r_{P1P2}}{r_{\max }-r_{P1P2}}\right)^{I}\right)\right)^{2}+\left(\Delta h_{\text{hand}}\right)^{2}}, r\geq r_{P1P2} \end{cases}\end{equation}

The initial parameters are used to calculate the chord length and passive torsion angle, and the curve of chord length and torsion angle in one cycle is obtained, as shown in Figs. 10 and 11. The brown area in Figs. 10 and 11 is the real-time change curve of the chord length along the span in one cycle. For the convenience of explanation, take the chord length and the passive torsion angle curve at any time for explanation, as shown in the red curve in the figure. The four points marked in the red curve correspond to the chords P ₁ P ₁₀, P ₆ P ₁₁, P ₂ P ₅ and the wing tip in Fig. 8. The three curves enclosed by these four points represent the chords in the three regions of the wing membrane in Fig. 8. It is easy to understand by comparing Figs. 10, 11 and Fig. 8.

Figure 10. Chord length change in one cycle. (The abscissa should be the span length of the membrane, but the active deformation of the wings will cause the span length to change in real time. In this paper, the BEM is used to divide the span into 1000 strips, and strip algebras are used to replace the span length).

Figure 11. Passive torsion angle change in one cycle.

3.1.2. Calculation of relative velocity

The relative velocity on the strip is composed of three parts: the passive torsion velocity caused by the flexible deformation of the membrane, the flapping velocity caused by the reciprocating flapping of the wing, and the freestream velocity.

The passive torsional velocity v _nt(r, t) is:

(19) \begin{equation} v_{\textrm{nt}}\left(r,t\right)=-\frac{1}{4}c\left(r,t\right)\frac{\partial \varphi (r,t)}{\partial t}\end{equation}

The flapping velocity of the wing v _pf(r, t) at the strip with the distance from flapping axis r is:

(20) \begin{equation}v_{\text{pf}}\left(r,t\right)=-r\frac{d}{dt}\theta _{o}\left(t\right)=-\pi fA_{f}r\cdot \cos (2\pi {ft})\end{equation}

The projection velocity v _lf (t) of the freestream velocity on the surface x _m O _m z _m can be shown in Fig. 12, which can be obtained:

Figure 12. Projection velocity of freestream velocity.

(21) \begin{equation}v_{\text{lf}}\left(t\right)=v_{f\infty }\cdot \cos \beta _{2}\end{equation}

(22) \begin{equation}\cos \beta _{2}=\sqrt{1-\left[\sin \left(\frac{\pi }{2}-\beta \right)\sin {\theta _{o}}\right]^{2}}=\sqrt{1-\sin ^{2} \theta _{o}\cos ^{2} \beta }\end{equation}

(23) \begin{equation}\cos \beta _{3}=\frac{\cos \left(\frac{\pi }{2}-\beta \right)}{\cos \beta _{2}}=\frac{\sin \beta }{\cos \beta _{2}}\end{equation}

where v _f∞ is freestream velocity, β ₂ is the angle between freestream velocity and surface x _m O _m z _m , θ ₀ is the flapping angle, β is the angle between the flapping plane and the horizontal plane, β ₃ is the angle between v _lf(t) and axis O _m x _m .

In the plane x _m O _m z _m , the combination of the three velocities is shown in Fig. 13:

Figure 13. Composite of relative velocity.

The relative velocity and the components of velocity on the x_m and z_m axes are shown below:

(24) \begin{align} \begin{cases} v_{Rx}\left(r,t\right)=-v_{\text{lf}}\cos \beta _{3}+v_{\text{nt}}\sin \varphi \\[5pt] v_{Rz}\left(r,t\right)=v_{\text{pf}}+v_{\text{lf}}\sin \beta _{3}-v_{\text{nt}}\cos \varphi \\[5pt] v_{R}\left(r,t\right)=\sqrt{v_{Rx}^{2}+v_{Rz}^{2}} \end{cases} \end{align}

3.2.1. Translational circulation aerodynamic force

The translational circulation generated by the membrane consists of two parts, one is the circulation generated by the freestream flow acting on the membrane. The other part is the change of the velocity caused by the flapping of the wing, resulting in the circulation. The quasi-steady translational aerodynamic force can be obtained through the vector synthesis of lift and drag in the membrane coordinate system. The tangential component of translational aerodynamic force is very small and can be ignored. The two-dimensional quasi-steady constant lift and drag coefficient of the membrane can be given according to the empirical equation by Dickinson [Reference Dickinson24]. However, the equation summarized by Dickinson is an empirical equation based on insect experiments, and the model in this paper is a bat. Due to the large difference in size, this equation cannot be directly used to calculate the bat model, and the coefficient needs to be corrected according to the experimental data of bats. Using the data of the bat lift coefficient and drag coefficient obtained by Colorado [Reference Colorado, Barrientos, Rossi and Breuer25], the empirical equation summarized by Dickinson is modified. The revised Eqs. (25) and (26) are as follows:

(25) \begin{equation}C_{L}=-8.487-9.588\,\sin\, (0.0233\alpha +4.304)\end{equation}

(26) \begin{equation}C_{D}=0.2844+0.2262\cos\, (0.09426\alpha +3.652)\end{equation}

The lift and drag on each strip are:

(27) \begin{equation}dL_{\text{trans}}=\frac{1}{2}\rho v_{\text{trans}}^{2}C_{L}{cdr}\end{equation}

(28) \begin{equation}dD_{\text{trans}}=\frac{1}{2}\rho v_{\text{trans}}^{2}C_{D}{cdr}\end{equation}

(29) \begin{equation}v_{\text{trans}}=\sqrt{v_{f}^{2}+v_{\infty }^{2}}\end{equation}

(30) \begin{equation}\alpha _{\text{trans}}=\mathrm{a}\tan \left(\frac{v_{\infty }\cdot \sin \alpha +v_{f}}{v_{\infty }\cdot \cos \alpha }\right)\end{equation}

Aerodynamic force on each strip is:

(31) \begin{equation}dN_{\text{trans}}=dL_{\text{trans}}\cdot \cos \alpha _{\text{trans}}+dD_{\text{trans}}\cdot \sin \alpha _{\text{trans}}\end{equation}

3.2.2. Rotational circulation aerodynamic force

The membrane is made of flexible material. In the process of reciprocating flapping, it will be affected by aerodynamic force, which will cause passive torsion, resulting in an angular velocity along the leading edge axis. However, the translational aerodynamic model cannot capture the amount of circulation caused by the continuous change of torsional angular velocity. Especially when the area of the membrane is large and the material is soft, this rotational circulation cannot be ignored [Reference Sane and Dickinson26]. Therefore, aerodynamic force caused by the passive rotation needs to be considered. The aerodynamic force of passive rotation on any strip is:

(32) \begin{equation}dF_{\mathrm{rot}}=\frac{1}{2}\rho \omega _{\theta }\omega _{t}C_{r}{rc}^{2}{dr}\end{equation}

(33) \begin{equation}C_{r}=\pi \!\left(0.75-\hat{x}_{0}\right)\end{equation}

(34) \begin{equation}\omega _{\theta }=\frac{d}{dt}\theta =\pi fA_{f}\cos\!(2\pi ft)\end{equation}

(35) \begin{equation}\omega _{t}=\frac{d}{dt}\varphi\end{equation}

where $\omega$ _θ is the flapping angular velocity, $\omega$ _t is the passive rotation angular velocity, C _r is the rotational circulation aerodynamic force coefficient, x ₀ is the dimensionless length between the pitch axis and the leading edge axis.

3.2.3. Added-mass force

The above is to use the attachment vortex mechanism to calculate the circulation aerodynamic force (translation and rotation). But in many cases, especially when the wing is large, the non-circulation mechanism also has a great influence on the aerodynamic force. Therefore, it is necessary to consider the factors of non-circulation mechanism in aerodynamic calculation.

In aerodynamic research, non-circulation mechanism generally refers to added-mass force. The method mentioned by Wood [Reference Whitney and Wood23] is used to derive and calculate the normal added-mass force on the plane of the membrane:

(36) \begin{equation}dF_{\mathrm{add},N}=-m_{\mathrm{add},\mathrm{mid},N}\alpha _{N}-m_{\mathrm{add},\varphi,N}\ddot{\varphi }\end{equation}

And:

(37) \begin{equation}\alpha _{N}=-r\left(\ddot{\theta }+\dot{\varphi }\dot{\phi }\right)=-r\ddot{\theta }\end{equation}

(38) \begin{equation}m_{\mathrm{add},\mathrm{mid},N}=\frac{1}{4}\pi \rho c^{2}\end{equation}

(39) \begin{equation}m_{\mathrm{add},\varphi,N}=\frac{1}{4}\pi \rho c^{2}c_{\mathrm{dis}}\end{equation}

where c _dis is the distance between the pitch axis and the midpoint of the strip.

So, the normal added-mass force on the plane of the membrane is:

(40) \begin{equation}dF_{\mathrm{add},N}=\frac{1}{4}\pi \rho c^{2}\ddot{\theta }r^{2}dr-\frac{1}{4}\pi \rho c^{2}c_{\mathrm{dis}}\ddot{\varphi }rdr\end{equation}

3.2.4. Inertia force

Some scholars have conducted experimental studies on aerodynamic forces and found that inertia force plays an important role in aerodynamic test experiments, especially in the vicinity of wing pitch inversion, where the inertial force is almost equivalent to the additional mass force [Reference Ke, Zhang, Shi and Chen27]. Therefore, it is necessary to consider the contribution of inertia force.

Because the membrane is connected to the wing skeleton and the shape is irregular, the BEM still needs to be used to calculate the inertial force on each strip. The inertial force of each chord is:

(41) \begin{equation}dF_{\text{inert},x}=-dm_{\text{strip}}\cdot \dot{v_{Rx}}\!\left(r,t\right)\end{equation}

(42) \begin{equation}dF_{\text{inert},z}=-dm_{\text{strip}}\cdot \dot{v_{Rz}}\!\left(r,t\right)\end{equation}

And:

(43) \begin{equation}dm_{\text{strip}}=\rho _{\mathrm{mem}}\cdot h_{\mathrm{mem}}c\!\left(r,t\right)dr\end{equation}

where dm _strip is the mass of each strip, $\rho$ _mem is the density of the membrane, h _mem is the thickness of the membrane.

The inertial force of each chord is:

(44) \begin{equation}dF_{\text{inert}}=\left[dF_{\text{inert},x},0,dF_{\text{inert},z}\right]^{T}\end{equation}

(45) \begin{equation}dF_{\text{inert},N}=-dm_{\text{strip}}\cdot \dot{v_{R}}\left(r,t\right)\end{equation}

3.2.5. Resultant aerodynamic force

Based on the above aerodynamic forces, the resultant aerodynamic forces can be obtained as:

(46) \begin{equation}{dF_{\text{res}}}={{dN_{\text{trans}}}+dF_{\text{rot}}}+d{F_{\text{add},N}}+dF_{\text{inert},N}\end{equation}

The normal aerodynamic force obtained is only in the wing coordinate system, which needs to be transferred to the body coordinate system. The rotation matrix is used to transform the aerodynamic coordinate on the strip to the body coordinate system O _b -x _b y _b z _b . The rotation transformation matrix R is:

(47) \begin{equation}R=R\left(x,\theta _{o}\right)R\left(y,\varphi \right)=\left(\begin{array}{c@{\quad}c@{\quad}c} \cos \varphi & 0 & \sin \varphi \\[5pt] \sin \theta _{o}\sin \varphi & \cos \theta _{o} & -\sin \theta _{o}\cos \varphi \\[5pt] -\cos \theta _{o}\sin \varphi & \sin \theta _{o} & \cos \theta _{o}\cos \varphi \end{array}\right)\end{equation}

So, the aerodynamic force in the body coordinate system is:

(48) \begin{equation}dF_{\text{body}}=\left[\begin{array}{c} F_{\text{Drag}}\\[5pt] F_{\text{Lateral}}\\[5pt] F_{\text{Lift}} \end{array}\right]=R\cdot dF_{\mathrm{res}}\vec {k}=\left[\begin{array}{c} -\cos \theta _{o}\sin \varphi \cdot dF_{\mathrm{res}}\\[5pt] \sin \theta _{o}\cdot dF_{\mathrm{res}}\\[5pt] \cos \theta _{o}\cos \varphi \cdot dF_{\mathrm{res}} \end{array}\right]\end{equation}

The three forces in Eq. (48) are the components of the aerodynamic force of the single wing on the x, y, and z axes, which represent drag, lateral force, and lift. The lateral force is not of great significance due to the equivalent reverse effects of the two wings. This study only considers the lift (dF _Lift) and drag(dF _Drag). The lift and its components obtained according to basic parameters are shown in Fig. 14. The drag and its components obtained according to basic parameters are shown in Fig. 15.

Figure 14. Lift and its components.

Figure 15. Drag and its components.

As can be seen from the Fig. 14, aerodynamic force is mostly derived from translational circulation caused by freestream flow and wing flapping. The average lift in one period accounts for 75.35%. The aerodynamic force caused by passive rotation only accounts for 14.01% of the whole, and all of them are negative lift. The added-mass force has great influence on aerodynamic force, accounting for 10.06%, and produces positive lift. However, due to the reciprocating flapping of wings, the average inert lift in one period only accounts for 0.58% of the whole, which is almost zero. It has almost no effect on the average lift, but only affects the transient lift.

As can be seen from the Fig. 15, aerodynamic force is mostly derived from translational circulation, the added-mass force and rotational circulation, while the inertia force has little influence. Among them, the translational circulation mainly affects the average drag, and at the same time, due to the existence of the translational circulation, the positive average drag will be generated. The added-mass force mainly affects the instantaneous drag, but only 7% of the average drag, and produces negative drag. The rotational circulation mainly affects the average drag, accounting for 18% of the whole, but it produces a negative average drag. However, due to the reciprocating flapping of the wing, the average drag in one cycle only accounts for 3% of the whole and has almost no influence on the average drag and instantaneous drag.