The dynamics of passive feathering rotation in hovering flight of bumblebees

The fluid-structure interaction problem of the flapping wings of bumblebees is considered, with focus on the action of elastic joints between wings and body. Morphological measurements and kinematic reconstruction of the wing motion using synchronized high-speed video recordings are described. They provide the necessary input data for numerical modelling. In particular, for the first time, the moments of inertia of bumblebee’s wing are determined using realistic mass distribution. A computational fluid dynamics solver is combined with a dynamical model that describes the wing motion. The model consists of the wings approximated as flat plates, connected with the body by elastic hinges. The results of high-resolution numerical simulations are presented. The hinged plate model produces realistic feathering motion and accurate time-average estimates of the aerodynamic performance in hover, despite some discrepancy in the instantaneous values of aerodynamic forces compared with the fully prescribed model. A parameter sweep reveals that the hinge is not exactly tuned to maximum efficiency during hovering flight, but slightly offset away from the maximum. ∗Corresponding author: dkolom@gmail.com; dkolomenskiy@jamstec.go.jp Preprint submitted to Journal of Fluids and Structures March 22, 2019

: Bumblebee wing morphology. The wing outline shape, including the mean contours, ±1 and ±2 standard deviation intervals calculated using 20 forewing and 18 hindwing samples, the outline shapes of selected intact samples, and the closed-contour approximation used in the CFD simulations. Red, green and blue lines show the position of the veins on the wing. Three different colors are used for visual distinction between different lines that have different numbers, each having its distinct constant diameter in the model. The corresponding biological classification (Michener, 2007, page 50) is shown in the right, for reference. The black and white marker shows the center of mass situated at x c /R = 0.379 and y c /R = −0.019. Some of the wing samples were relatively intact and others were worn.
where R is in meters and J xx , J yy , J xy are in kg m 2 . Afterwards, the values obtained from (3) are used as input data for CFD simulations.

122
Since veins account for more than 60% of the wing mass (Kolomenskiy 123 et al., 2019), and our calculation assumes circular piece-wise constant cross-124 section, the approximation error may be significant.   : Sample frames from two synchronized video recordings: camera 1 (top row) and camera 3 (bottom row). Frames 517 and 520 correspond to downstroke, frames 525 correspond to upstroke. Theoretical rigid wing contour lines, shoulder points (plus signs) and body markers (dots) are superposed on the images. proximated as a single solid flat plate that can rotate about the hinge point 157 at the shoulder, therefore its orientation with respect to the body is fully 158 described with three angles. The body is also assumed rigid, therefore, it is 159 straightforward to relate the position of the shoulder points in the labora-    In total, 7 video sequences have been analyzed, that correspond to hov-185 ering flight of different individuals. These videos have been published in an 186 online repository . The measured parameters in-187 clude the wing length R, the wing beat frequency f , the body inclination 188 angle β and the anatomical stroke plane angle η, see Table 1, as well as the where AR = R 2 /S = 3.66 is the aspect ratio evaluated using the area S = 193 wing dxdy = 0.273R 2 of the intact wing in Fig. 1, the radius of the second 194 moment of area calculated using the same intact wing contour, the reference average wing speed and the Reynolds number The average wing-tip speed is equal to U t = 1.75U 2 and the wing-tip Reynolds 198 number is equal to Re t = 1.75Re 2 . The air density and kinematic viscosity at   (FZ-300i, A&D, Japan) and measured the wing length using a digital caliper.  In addition, we weighed individual #4 immediately after recording its 212 flight. The result is included in Table 1. It differs by 18.6% from the isometric 213 fit (8), i.e., falls within one standard deviation interval above the isometric 214 fit. This is shown using a large cross marker in Figure 2  The computational approach followed in the present study is in continuity 220 with our previous work (Engels et al., 2016b,a;Ravi et al., 2016). We employ 221 FluSI 1 , a Fourier pseudo-spectral solver with volume penalization (Engels 222 et al., 2016b). In the simulation, the bumblebee is approximated by three isometrically scaled such as to match the wing length R as given in Table 1.

233
In all numerical simulations in the present study, the body is fixed in the laboratory reference frame. Prescribed time-periodic functions φ(t) and θ(t), as shown in Fig. 6, determine the position of the wing tip. Rotation of the wing about its longitudinal axis is described by α(t) which is determined from an elastic hinge model similar to that proposed by Whitney and Wood (2010). The model employs an equivalent linear torsional spring-damper element as an abstraction for the combined effect of the compliance and structural damping of the muscles, shoulder joints and proximal parts of the wings. We infer the model coefficients from minimizing the discrepancy between simulated and measured time series of α(t), as we explain later in Section 3.3. In the numerical simulation, the feathering angle α(t) is determined from the equation of passive feathering motion, where M aero is the aerodynamic pitching moment, K is the elastic hinge the aerodynamic ground effect is negligible .

255
Strong fluid-structure coupling is used. Equation (9) is transformed 256 in a system of two first-order differential equations and discretized using 257 the second-order Adams-Bashforth scheme. Since the same time marching 258 scheme is also used in the Navier-Stokes solver, it is straightforward to eval-  Quasi-periodic regime is in our case is reached after three wing beat cycles. as shown in Fig. 6 (b). In the second, α(t) is passive, i.e, it is modelled using 282 equation (9)  for the individual #2, but almost zero for #1 and small negative for #3.

295
Therefore, this small residual phase lag seems to be a measurement error 296 rather than a modelling artifact.

302
The vertical aerodynamic force, Fig. 7 (b), shows two distinct peaks around 303 the middle of each translation phase. They are larger but narrower in the case 304 of passive rotation, which is consistent with the differences in α(t) discussed 305 above. The pointwise difference in F az (t) is particularly large on upstroke, 306 when α differs by as much as 20 • . However, this overestimate cancels out 307 with the underestimates before and after the peak, such that the wingbeat 308 time-average force, as shown with the dashed lines, is almost identical in the 309 two cases and it is 8% less than the body weight (the last column in Table 1).

310
The horizontal force, see Fig. 7(c), shows similar trends with the difference 311 that the peak associated with the upstroke is negative, such that the time- kinematics.

318
In summary, the passive wing rotation model (9)

328
Let us now discuss sensitivity of the results to the elastic hinge parameters 329 K, C and α 0 . Since these parameters cannot be measured directly, they will 330 be evaluated by fitting the model to the experiment data. It will be insightful 331 to see the influence of each parameter separately before solving the full opti-332 mization problem. Contributions of each term in the passive feathering equation (9) are illustrated in Fig. 10. The aerodynamic pitching moment M aero is small during the reversals near t/T = 0 and 0.5, when the translation velocity of the wing is small. M aero is large positive in the middle of downstroke, it is large negative in the middle of upstroke. The inertial pitching moment In this section, we look for the optimal values of K, C and α 0 that 355 minimize the cost function where α(t) is the time evolution of the feathering angle obtained from the 357 numerical simulation using (9), and α exp (t) is the feathering angle measured 358 in the experiment.

359
It may be expected that the hinge parameters vary among different in-  Supplementary material S3 also supports this assumption.

374
In view of the isometric scaling (1) of the wing mass versus length, it is 375 reasonable to introduce a similar isometric scaling for the hinge stiffness, that holds for flexible-plate hinges with thickness, width and length scaled 377 linearly with R (Whitney and Wood, 2010). We refer to K * as the hinge 378 stiffness factor. It can be regarded as a composite material property. The

379
cost function e for all individuals of the test set is plotted in Fig. 11 with 380 respect to K * . The optima K * opt,i , where i = 1, ..., 4 is the individual index, 381 are indicated with circles and included in Table 2. The average stiffness 382 factor plus/minus standard deviation is equal to

418
The energetic efficiency of hovering flight can be measured using the figure   419 of merit F M = P ideal /P , where P ideal is the ideal power determined by the 420 Rankine-Froude momentum theory as P ideal = 2ρw 3 0 A 0 , where A 0 = 2ΦR 2 421 is actuator disc area and w 0 = L/2ρA 0 is the induced velocity, and P = 422 1 T T 0 P a (t)dt is the mean aerodynamic power from the numerical simulation.

423
After arithmetic simplification, we obtain a short formula The numerical results are shown in Fig. 13. For all individuals, the plots of

438
In the present study, it is shown that a single hinged plate model, orig-439 inally designed for diptera, also provides a reasonably accurate approxima- shows no clear trend, but the mean value is close to zero, α 0 ≈ 0.

453
From the aerodynamic perspective, it is found that passive feathering 454 provides the required lift-generation capacity for a realistic energetic cost. It 455 is interesting that the hinge parameters are not exactly tuned to maximum 456 efficiency during hovering, but the stiffness is slightly larger than the optimal 457 value. We conjecture that the difference may be interpreted as a safety factor 458 that helps to avoid abrupt decrease in the efficiency when flight conditions  Table 1.

527
As a next step, we track the wing tips, reconstruct the wing tip trajec-528 tories, convert them to the body reference frame and best-fit a plane, in the 529 least-mean-square sense. The morphological stroke plane angle is determined 530 as the angle between the normal to that plane and the body longitudinal axis.

531
The stroke plane, in our definition, is inclined by the same angle to the body, periodic, with small deviations that may be due to actuation, wing-wake 558 interaction and measurement errors. The next processing step consists in 559 low-pass filtering the data at 450 Hz using the 4th order Butterworth filter 560 and upsampling the result on a 100-times finer grid using spline interpolation.

561
Thus we discard those points that produce unrealistically large accelerations.

562
The resulting time profiles are shown with dotted and dashed lines that 563 correspond to the left and the right wing, respectively.  The time evolution of φ, α and θ is described with less than 1 degree error 577 using, respectively, 4, 5 and 4 harmonics. These coefficients are used as input 578 data for the CFD simulation. tions of the passive rotation model (9). The flapping motion in the direction 583 of φ was driven by a piezoelectric actuator, while passive rotation in α was 584 allowed by an elastic hinge. These angles, as well as the small out-of-plane 585 deviation θ due to compliance, were measured simultaneously with the ver-  The 'split-cycle' case is presented in Fig. B.19. Similarly, the peaks of α are slightly overestimated, which leads to overall lower lift of the passive 642 feathering model. The time evolution of F shows oscillatory behavior due to