1 Introduction

As is well-known,a flying animal must generate lift to support its weight and thrust to overcome the drag in forward flight. A fixed wing at an angle of attack (AOA) in the range from 10^◦ to 15^◦ can produce sufficient high lift but cannot provide thrust. This is why an animal cannot move for long time or long distance by gliding. For flying species,the aerodynamic forces are produced mainly by the flapping wings in the flight. In the past decades,the flow physics of flapping motion has been investigated extensively,especially the mechanism of high lift and thrust generation in insect flight^{[1, 2]}. However,in bat flight,the deformation of a flapping wing may have new mechanism due to its special structure.

A bat wing consists of flexible muscularized membrane and upper limbs. First,the membrane attaches all the lateral border along its body,from the neck to the ankle^[3]. This arrangement limits the supination/pronation motion of the bat wings,which has been commonly observed and investigated in insect flapping flight. In fact,the AOA at the root,which is also the body tilt angle,is very small and changes little during the wing stroke,while the AOA at the wing tip is large and varies greatly during the wing beat^{[4, 5]}. The variation of the AOA distribution from the wing root to the wing tip is called twist-morphing or twisting. Second,the twist-morphing is complex and controlled actively by the upper limbs,similar to human’s. The variation of wing twist or camber can also be caused by the joint motion controlled by the wing muscle^[3]. The unique bat wing morphology suggests a high potential to adjust the wing morphing according to the aerodynamic demand^[5].

There have been many experimental studies on bat flight. The flapping parameters,such as the amplitude,the frequency,and the stroke plane angle versus the flight velocity have been investigated broadly^{[6, 7, 8, 9]}. Many steady aerodynamic and momentum theories have been used to predict the lift and drag in bat flight^{[6, 7, 8]}. With the development of the particle image velocimetry (PIV) technology,the wake of a bat can be measured to estimate the forces^{[10, 11, 12]}. Wolf et al.^[4] and Busse et al.^[5] investigated the deformation of bat wings. However,up to now,few studies have described the wing morphing and the mechanism of aerodynamic force generation^[13].

The purpose of this paper is to find out the mechanism of aerodynamic force generated by the twist-morphing wing. As a preliminary study,an unsteady panel method is developed to predict the aerodynamic forces,and a three-dimensional (3D) plate with the aspect ratio AR of 3 is used to model the bat wing with the linearly distributed AOA along the span-wise direction. In the discussion,the twisting motion is compared with the pitching motion of a flapping rigid-plate. The aerodynamic mechanism of the twisting wing is concluded finally. 2 Numerical methods

Yu et al.^{[2, 14]} proposed that for highly unsteady flow,the Strouhal number (St) is the dominant parameter governing the character of the flow and the viscous effect becomes in- significant (St>=1/Re,where Re is the Reynolds number). In bat flight,St ∈ [0.25,1.25] and Re ∈ [10³,10⁵]. Furthermore,it is reasonable to consider the flow excited by the bat wing to be incompressible. Therefore,the flow field is potential and incompressible,and the panel method can be used to analyze the problem,by imposing the Kutta conditions on the leading,lateral, and trailing edges. In this section,the panel method is introduced briefly,and the verification and validation are shown subsequently. 2.1 Panel method

The panel method is a highly developed approach,which has been used successfully to inves- tigate the plunging and/or pitching airfoil,bird/insect wing flapping,and fish swimming^{[15, 16, 17, 18, 19, 20]}. Here,we just give an overview of the method. The reader is refereed to Ref.^[21] for more details on the theory and implementation.

One of the challenges in the panel method is how to treat the separation. Fortunately,the bat wing is so thin that the separation position is confined near to the leading edge. A zero thickness plate,divided into the finite number of linear panels N_b,is used to model the bat wing. The constant strength distribution of doublet μ is situated at the midpoint of each panel. At each time step,a new row of wake panels is added to the wake sheet from the leading, lateral,and trailing edges. The doublet strength of the new wake panels is determined by the Kutta condition,i.e.,the instantaneous strength equals that of the nearby doublet at the leading,lateral,or trailing edge. At subsequent time steps,all of the wake panels are convected downstream with the local flow velocity. The “zero normal flow” boundary condition is used to determine the doublet strength on the plate.

Another challenge is to avoid the numerical singularity when the wake panels are too close or they are near the plate,as the induced velocity becomes too large to lead the computation to divergence. Especially,it happens frequently near the leading edge. The desingularization factor δ is introduced in our in-house panel method program. The desingularized velocity induced by the vortex segment l is

where d is the distance from the point to the segment,and u is the induced velocity determined by the Biot-Savart law. This method can also avoid the wake panels from the leading edge going through the plate.

When the doublet strength has been determined,the pressure can be resolved by the un- steady Bernuolli’s equation. Then,the aerodynamic force can be derived by

where p_u and p_l are the pressures of the upper and lower surfaces of the plate,respectively. S_k and n_k refer to the area and the unit normal vector of the panel k,respectively. Then,the lift coefficient C_L and the drag coefficient C_D can be obtained as follows:

where F_z is the vertical upward force (lift),Fx is the horizontal backward force (drag),pis the air density,U is the inflow velocity,and S is the plate area. 2.2 Numerical convergence

The numerical parameters such as the desingularization factor,the grid resolution,and the time step may affect the accuracy and the convergence of the results. It is necessary to assess their effects before the validation.

A 3D pitching and plunging flat plate with A_R = 2 is a typical sample to examine the numerical convergence. The kinematic equations are given by

where h(t) is the plunging motion,while α(t) is the pitching motion with pitch pivot on the 0.25 chord from the leading edge^[22].

Figure 1 shows the time-dependent lift coefficient C_L versus the desingularization factor δ, where T is the period and t is the time. It is illustrated in Fig. 1 that the numerical oscillation can be eliminated when δ≥0.35 with the surface grid 20 × 20 and the time step dt = 0.05.

Figure 2 shows the effects of the grid resolution and the time step. It indicates that the results are convergent. For the subsequent computation,the desingularization factor δ is set to be 1.0,the time step dt is set to be 0.05,and the grid is set to be 20×30 (for A_R greater than 2). 2.3 Validation

To validate the method,we perform quantitative comparisons with previous experimental data^[23] and numerical results^[24].

Firstly,we consider a flat plate with A_R = 2,pitching from 0^◦ to 90^◦ with a constant angular velocity and the pitch pivot point at the leading edge. The motion is depicted as α= 2κt,where k is the reduced frequency of the pitch,α is the pitching angle,and t is the time. Figure 3 shows the profiles of C_L and C_D versus the angle and the time. It can be seen that the pitching motion agrees well with the experimental data (the Reynolds number is 20 000 in the experiment) obtained by Kenneth et al.^[23].

Secondly,a comparison with a computational fluid dynamics (CFD) sample^[24],where a wing with A_R = 1 revolves in the quiescent flow,is performed. The wing root is extended out 0.52 chord from the axis of rotation with the chord oriented at an angle of 60^◦ relative to the sweeping plane. The wing is initially at rest,and then accelerates to a constant rotational rate as it revolves around the axis. The motion of the revolving angle Φ follows^[24]

where Τ is the non-dimensional time,a = 2.6,Ω₀^* = 0.98,Τ1 = 2,Τ2 = 8.4,and Φ0 = 2π.

Figure 4 shows the results of the time-dependent lift and drag coefficients obtained by the present method and Ref. ^[24]. It shows that the present panel method can give a good prediction (Re = 30 000)^[24].

3 Kinematic models

As mentioned in Section 1,the morphing of the bat wing is very complex. To investigate the aerodynamic forces generated by the wing,a simplified flapping model and a twisting model are given in this section in mathematical description. As a preliminary study,a rectangular plate is used to model the bat wing whose membrane is thin enough. 3.1 Flapping model

Flapping is the basic motion in bat flight,which can be described in a body-fixed coordinate system (x,y,z). In Fig. 5,the plate with the chord length c and the span length b flaps around the wing root (the flapping pivot),the AOA is defined as the angle between the chord and the incoming flow (along the x-axis),and the flapping angle θ is defined as that between the plate and the plane determined by the flapping pivot and the y-axis.

The flapping motion with respect to the angular movement is almost sinusoidal in bat flight^{[6, 25, 26]},based on which we can get

where θ0 is the average flapping angle,θA is the flapping amplitude,! is the angular frequency, and t is the time. For symmetric flapping,θ0 ω is set to be zero. 3.2 Twisting model

Because of the limitation of the special arrangement of the bat wing,the AOA at the wing root is fixed as the tilted angle of the body,and it changes slightly in the forward flight. The changing AOA is performed with the twisting motion by the flexible wing membrane and the hand of the bat. To simplify the problem,the twisting is defined as the AOA varying from the root to the tip with the linear distribution,which is illustrated in Fig. 6. A local coordinate system (ξ,η,ζ) is set on the shoulder of the twist-morphing plate,and the twisting pivot is fixed at the leading edge,the η-axis. Then,the AOA distribution can be expressed by

where α_tip is the tip AOA,and r is the distance from the root.

    Since αroot changes slightly during flapping and its value is very small,we set it to be a constant obtained by experiments^{[4, 5]}. We assume that αtip follows a simple harmonic motion based on the periodic deformation of the flapping wing,i.e., where αtip0 is the average value,αA is the morphing amplitude,and Φ is the phase difference relative to the flapping motion. The greater αA is,the greater the twist-morphing is. 3.3 Kinematical parameter selection

The kinematical parameters used here are obtained from Refs. ^[5] and ^[9]. The detailed information is shown in Table 1 and Table 2,where S_P is the full wing span (tip-to-tip),and c is the mean chord. Here,A_R = S_P/c. It contains two wings and the body attributions,and is more than twice as large as that with one wing. ω is the non-dimensional angular frequency defined by ω = 2πfc/u. θ_A is the flapping amplitude.

    Based on the experimental data shown in Tables 1 and 2,we choose A_R to be 3,ω to be around 1,and θ_A to be around 40^◦. According to Refs. ^[4] and ^[5],the body tilt angle α_rootis estimated around 7^◦,the tip angle α_tip0 is estimated around 15^◦,and the amplitude of the twisting α_A is estimated around 28^◦. All the parameters are listed in Table 3. 4 Results and discussion

In this section,we will first study the lift/drag affected by the flapping amplitude and the twisting amplitude,and then discuss the aerodynamic mechanism of the twisting model.

According to Section 2,the resolution of the surface grid is 20×30 with 20 in the chord-wise direction and 30 in the span-wise direction (see Fig. 7).

4.1 Effects of flapping and twisting

There are six parameters in the twisting model. However,in this section,only two param- eters are considered for their main contribution in controlling the flapping and twisting. One is the flapping amplitude θ_A,and the other is the morphing amplitude α_A. Based on Table 3, the effects of these two parameters on the force generation are studied and the time-averaged lift/drag coefficient is defined by

where F represents C_Lor C_D,and T is the flapping period. 4.1.1 Forces generated by flapping wing

Flapping is a common motion in animal flight. Can a bat get sufficient lift and thrust by this motion in its forward flight? The rigid-wing with a basic twist deformation is used to learn how the flapping motion affects the forces. The flapping amplitude θA is set to vary from 0^◦ to 40^◦,and the twisting amplitude is set to zero,i.e.,α_A = 0. The other parameters are the same as those in Table 3.

As shown in Fig. 8,a fixed wing,whose flapping amplitude is zero,can produce lift and drag with the values about 0.62 and 0.12,respectively. Obviously,the lift cannot balance the weight of the bat and the drag cannot support the forward flight.

When θA increases,the lift and drag generated by the flapping wing are different from those of a fixed wing. During the downstroke (0≤t/T ≤ 0.5),the larger the flapping amplitude θA is,the greater the forces are generated. During the upstroke (0.5≤t/T ≤ 1),the lift decreases and the drag turns to the thrust. But the averaged lift and drag increase with the increase in the flapping amplitude (see Fig. 8(c)).

When the amplitude is close to 40^◦,the averaged lift coefficient is about 0.8 and the averaged drag coefficient is about 0.19. It means that,for a non-twisting wing,the flapping motion can increase the lift and drag simultaneously,but cannot provide the thrust in flight,although the thrust has appeared during the upstroke.

4.1.2 Forces generated by twisting flapping wing

In the last subsection,we learned that the thrust cannot be obtained by a flapping rigid-wing in the forward flight. How to generate sufficient thrust? For a bat,to search the effects of the wing morphing may be an effective way to answer this question.

In this subsection,the twisting motion is considered with the amplitude αA varying from 0^◦ to 53^◦. With the kinematic parameters listed in Table 3,the numerical results are shown in Fig. 9.

Obviously,the peak (during the downstroke) and the trough (during the upstroke) of the lift coefficient curve are decreasing with the increase in the twisting amplitude (see Fig. 9(a)), but the average value increases with its maximum ω close to 1.0 when the twisting amplitude is around 40^◦. Therefore,the twist-morphing can enhance the lift by as much as 25%,relative to the rigid flapping wing. When α_A is in the range from 20^◦ to 43^◦,the averaged lift increases slightly (see Fig. 9(c)).

The drag curve is much different from that of the flapping rigid-wing. When the morphing amplitude increases,the twisting wing produces higher negative drag. Then,thrust will occur during not only the upstroke but also the downstroke (see Fig. 9(b)). The time-averaged drag is negative when α_A is greater than 10^◦ (see Fig. 9(c)). Therefore,a flapping twisting-wing can generate the thrust,which has been surmised^[13]. Similar to the lift curve,the average thrust peak appears when the twisting amplitude is close to 40^◦.

Now,we can see that the twisting motion of the bat wing can not only enhance the timeaveraged lift but also decrease its fluctuation during the whole stroke. Especially,the motion can generate not the drag but the thrust,which is required in the forward flight.

4.2 Force generation mechanism of twisting wing

How does a twisting wing generate the thrust during the downstroke and the upstroke? And how does it decrease the fluctuation of the lift? As mentioned in the introduction,the mechanism of high lift and thrust generation in insect flapping flight have been investigated widely. The pronation and supination in flapping motion are important to generate the forces,especially the thrust. To reveal the aerodynamic mechanism of the twisting wing,the pitching motion of a rigid plate is introduced to compare the aerodynamic forces. 4.2.1 Pitching model (AOA at radius of gyration along span-wise direction)

3D flapping motion of the rigid wings has been investigated numerically^{[1, 22, 27, 28, 29]} and exper- imentally broadly^{[30, 31, 32]}. The pitching motion with the positive (or reverse) direction is called supination (or pronation),which is a mechanism to alter the AOA of the wing.

It is found that the kinematic parameters at the radius of gyration are very fantastic. The physical parameters such as the speeds and the chord length at this radius are usually used as reference quantities for normalization^{[1, 33, 34]}. In the present study,the AOA of the pitching motion is the same as that of the twisting motion at the radius of gyration.

The radius of gyration,of the second moment of the area,is determined by r₂ = (∫ ^r2ds/S)^0.5,where r is the radial distance from the fulcrum and S is the wing area. For a rectangle plate, r₂ = 0.58b. Therefore,the AOA of a rigid pitching wing,α_r2,is set as

and the pitching pivot is at the leading edge. 4.2.2 Comparisons of forces produced by twisting wing and pitching plate

The curves of the forces generated by the twisting and the pitching plates are demonstrated in Fig. 10. For these two models,the lift (or drag) curves coincide well with each other. The results indicate that the twisting in bat flight is of the same function with the pitching,which is necessary motion in insect flight.

Why are the aerodynamic forces generated by twisting wing and the pitching wing similar? For the linear distribution of the aerodynamic force along the span-wise direction,the radius of gyration,r₂,is at the theoretical position of the pressure center. In our study,it is selected as the kinematic reference of the pitching,although the pressure is not distributed linearly for its 3D effects. The position of the pressure center of the twisting model is determined by

where r is the vector of position,n is the outward unit normal vector of the upper surface, R_F is the position of the pressure center,Δ_p is the pressure difference of the lower and upper side,and F is the concentrated force. As a special case,the roll moment is considered,which approximates to M_x when the tilt angle of the body is very small. Therefore, where R₂ is the arm of the force,regarded as the center of the pressure approximately.

In the two cases where the twisting amplitudes are 15^◦ and 28^◦,we find that the value of the force arm is 1.85 (see Fig. 11),that means R₂ = 0.62b,which is very close to the radius of gyration for a rectangular flat plate,r₂ = 0.58b. Therefore,it is easy to understand why the lift/drag curves of the two types of motion have similar patterns.

Furthermore,it is indicated that the function of the twisting motion of the bat wing is analogous to that of the pitching motion of a rigid wing. In other words,the twisting motion changes the effective AOA of the bat wing like the pronation and supination of the rigid wing. During the downstroke,the “pronation” of the twisting wing makes the AOA become negative, but the effective AOA remains positive while it decreases. Thus,the thrust is generated as a component of the net force,and the lift is high enough although it is less than that without “pronation”. During the upstroke,the “supination” of the bat wing makes the AOA increase but the effective AOA is also smaller than that of the wing without “supination”. Therefore, the negative lift decreases and the thrust is generated.

In short,a bat can generate the sufficient thrust with the flexible twisting wing. Simultaneously,the lift is enhanced and its fluctuation is decreased,relative to the rigid flapping wing.

5 Conclusions

In this paper,the aerodynamic forces generated by the twisting model-bat-wing with the hypothesis that the AOA of the twisting is linearly distributed are investigated by the panel method. The results indicate that a flapping rigid-wing can produce higher lift than a fixed wing,but there is no time-averaged thrust. To overcome the limitation of the arrangement that the bat wing attaches its body side,the flexible wings must be twisted to generate the thrust during the downstroke and upstroke. Simultaneously,the mean lift is enhanced by as much as 25% while the fluctuation of the lift is decreased.

The amplitude of the twisting is a key parameter to control the lift and thrust generation. When the amplitude is from 20^◦ to 43^◦,the averaged lift increases slightly and the maximum appears around 40^◦. The averaged thrust is positive when the amplitude is greater than 10^◦ and the maximum appears when it is close to 40^◦.

To reveal the aerodynamic mechanism of the twisting wing,a pitching motion (with prona- tion and supination) of a rigid-plate is brought in as a comparison. It is found that the aero- dynamic forces (lift/drag) generated by a twisting plate-wing are in similar patterns to those by a pitching wing,i.e.,the twisting in bat flight is of the same function with the prona- tion/supination motion in insect flight.

In conclusion,the twisting motion is very important in bat flight,which can make the bat wing not only generate the sufficient thrust,but also enhance the lift,just like the prona- tion/supination in insect flapping flight.

Acknowledgements The authors thank Professor Bing-gang TONG and Associate Professor Lin BAO for their helpful comments.