Aerodynamic performance of a bristled wing of a very small insect

Aerodynamic force generation capacity of the wing of a miniature beetle Paratuposa placentis is evaluated using a combined experimental and numerical approach. The wing has a peculiar shape reminiscent of a bird feather, often found in the smallest insects. Aerodynamic force coefficients are determined from a dynamically scaled force measurement experiment with rotating bristled and membrane wing models in a glycerin tank. Subsequently, they are used as numerical validation data for computational fluid dynamics simulations using an adaptive Navier–Stokes solver. The latter provides access to important flow properties such as leakiness and permeability. It is found that, in the considered biologically relevant regimes, the bristled wing functions as a less than 50%\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$50\%$$\end{document} leaky paddle, and it produces between 66 and 96%\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$96\%$$\end{document} of the aerodynamic drag force of an equivalent membrane wing. The discrepancy increases with increasing Reynolds number. It is shown that about half of the aerodynamic normal force exerted on a bristled wing is due to viscous shear stress. The paddling effectiveness factor is proposed as a measure of aerodynamic efficiency.


3
194 Page 2 of 13 thrips (Thysanoptera). They are active fliers, which implies that bristled wings produce enough force to support the animal body weight and propel it through the air (Yavorskaya et al. 2019;Cheng and Sun 2018;Zhao et al. 2019). Considering the bristled morphology as a biological adaptation, it is important to look into the potential benefits and penalties from the mechanical standpoint. The present study investigates into the aerodynamic aspect of this problem. Horridge (1956) conjectured that the smallest insects may have forsaken their larger relatives' airfoil action exploiting the lift force perpendicular to the direction of motion. Instead, they use a mechanism by which the drag on the upstroke is made less than that on the downstroke. Horridge's analysis built upon earlier experiments with a ten percent thick airfoil (Thom and Swart 1940) having the drag at low Reynolds number independent of the angle of attack. Then logically, Horridge hypothesized that bending of the bristles could be critical for producing the necessary aerodynamic asymmetry. For instance, single cells with cilia such as Paramecium use asymmetric power-and recovery strokes based on bending. Kuethe (1975) extended upon the idea about the key role of elastic deformation. However, recent high-speed videography in free flight, e.g., of a featherwing beetle Nephanes titan (Yavorskaya et al. 2019), a chalcid wasp Encarsia formosa (Cheng and Sun 2018) and a thrips (Zhao et al. 2019;Lyu et al. 2019) shows large variation in the wing orientation during one stroke cycle, and relatively little bending of the setae.
It is, therefore, plausible that the aerodynamic asymmetry is primarily achieved by changing the angle of attack, instead of bristle bending deformation. An idealized twodimensional computational study (Jones et al. 2015) demonstrated such possibility of drag-based mean force generation by a cyclic motion of a flat plate. Three-dimensional numerical simulations of Encarsia formosa (Cheng and Sun 2018) and thrips Frankliniella occidentalis (Lyu et al. 2019) with realistic wing kinematics also show cycle-average forces sufficient for body weight support, under the assumption that the wings are impermeable plates. However, until now, the force-generation capacity and aerodynamic function of the bristled wing morphology remains a largely open question.
Outwardly similar bristled appendages may function as virtually impermeable paddles or as leaky rakes, depending on the Reynolds number and on the geometrical parameters. This was pointed out in Cheer and Koehl's theoretical study (Cheer and Koehl 1987) based on Umemura's matchedasymptotic solution for a pair of cylinders (Umemura 1982), and confirmed in a recent numerical investigation of twodimensional linear arrangements of multiple cylinders by Lee et al. (2020). Three-dimensional numerical computations of the forces of a bristled wing have been carried out by Barta and Weihs (2006) using the Stokes flow approximation for a linear array of slender ellipsoids. However, later Navier-Stokes simulations by Davidi and Weihs (2012) showed that, in the relevant for insects range of the Reynolds number between 10 and 100, the Stokes approximation is inaccurate by a factor greater than two. Despite that complication, both studies agreed that suitably spaced sparse arrays of rods can approach an impermeable wing in terms of aerodynamic force generation.
A linear array of cylindrical rods has been studied experimentally by Sunada et al. (2002). Their experimental apparatus implemented protocols of rectilinear translation and rotation about the vertical axis in a water-glycerin tank. The array of rods in all cases produced less force than a solid rectangular plate with the same external dimensions. However, the force per wing surface area was larger for the array of rods. Sato et al. (2013) downscaled this configuration to the size of one millimeter and tested it in a bench-top wind tunnel. Zhao et al. (2019) performed wind-tunnel force measurements on a thrips wing, paralleled by computational fluid dynamics (CFD) analysis.
A separate line of research focused on the bristled wing leakiness as a mechanism to facilitate the clap-andfling interaction. Thus, Santhanakrishnan et al. (2014) and Jones et al. (2016b) concluded from two-dimensional Navier-Stokes simulations that bristles reduce the force to move the wings apart. Experiments with mechanical models in a glycerin solution by Kasoju et al. (2018) and Ford et al. (2019) showed that bristled wings can have higher lift to drag ratio during clap and fling than solid plate wings.
Virtually all previous computational and experimental estimates of the aerodynamic forces of bristled wings were based on highly simplified morphological representation. They covered a range of operating conditions and showed a variety of dynamical effects. However, it is not self-evident which effects are actually realized in the biological wings.
The objectives of our study are to implement and crossvalidate an experimental facility and a numerical simulation software for studying the aerodynamics of bristles wings of bio-realistic shape. We constructed a dynamically scaled model that accurately matches the shape of a featherwing beetle Paratuposa placentis in terms of the bristle size and orientation, and the shape of the central membrane (blade). The diameter of the bristles was selected taking into consideration the secondary outgrowths on the setae measured in that species, see Appendix 1. Likewise, we implemented Navier-Stokes simulations of this wing.
We consider a revolving motion, which is a useful simplified kinematic protocol for studying the aerodynamics of biological flapping wings (Usherwood and Ellington 2002;Wolfinger and Rockwell 2015;Jones et al. 2016a). A revolving setup cannot replace a flapping setup as an aerodynamic test bed, but it offers a convenience of having a smaller kinematic parameter space while preserving the important spanwise gradient of the velocity that is also present in flapping (wing tip moves faster than the wing root). The two essential parameters are the angle of attack and the Reynolds number. Thus, a two-dimensional parameter sweep is performed with a bristled and with an equivalent membrane wing, the aerodynamic forces are measured and the flow properties are analyzed.

Materials and methods
This section describes, both, the mechanical apparatus for the dynamically scaled experiment and the numerical simulation setup.

Geometrical model and kinematics of the wing
The mechanical wing model and the CFD model share the same morphological features that mimic the real insect wing. Modelling artifacts such as the external boundaries of the fluid domain and attachment at the root of the wings are different in the experiment and in the simulation.

Geometrical model
The wing is modelled after one of the tiniest beetles, Paratuposa placentis Deane, 1931. Adult individuals were collected in Vietnam in the Cát Tiên National Park. Morphological measurements were acquired using light microscopes (BX43, Olympus Corporation, Tokyo, Japan and SMZ168, Motic China group, Ltd., China) and a scanning electron microscope (SEM) JSM-6380 (JEOL, Tokyo, Japan), following the same protocol as described in an earlier study . The external morphology of one of the samples is shown in Fig. 1a.
The samples were measured in AutoCAD (Autodesk, Inc., USA) using images taken with the light microscopes and SEM. Wing lengths and distances between tips and bases of setae (bristles) were measured using the light microscopic photographs. The wing length was measured as the distance between the base of the wing and the apex (the most distant point on the setae fringe). These measurements were made on ten wings. Then, diameters of setae were measured using SEM images. We measured both the diameters of the stems of the setae and the external diameters of the setae including the lengths of the outgrowths in middle regions of the setae. The diameter measurements were made on 20 setae. The scaled mechanical model (Fig. 1b) and the CFD model ( Fig. 1c) have similar major morphological features such as the number of long setae, their orientation, position on the wing blade, and the shape of the wing blade contour. The model wing length, defined as the distance from the axis of rotation to the apex, is scaled up to R = 93 mm.
The blade is approximated as a flat plate with a uniform thickness h = 1 mm , since the SEM measurements show that its out-of-plane deviation is less than 3% of the wing length. It is cut from a steel sheet using a printed photographic image of the insect wing as a template. Long setae are modelled as long and straight circular cylinders, all having the same diameter b = 0.36 mm , fabricated from hardened steel wire. Thus, the ratio b∕R = 0.00388 is fixed equally in the experiment and in the simulation. Small setae near the wing base (root) are neglected in the experiment as they are unlikely to contribute to the fluid-dynamic forces. Note that the surface of each seta has a complex micro-relief composed of more-or-less regularly spaced outgrowths . Therefore, the bristle diameter b in this study corresponds to an effective diameter, which is about two times as large as the diameter a seta having the outgrowths removed. The drag of a cylinder with the effective diameter b is the same as for a thinner cylinder covered with outgrowths. The value of the effective diameter used in this study is based on the results of a dedicated towing tank experiment described in Appendix 1.
In the mechanical model experiments and numerical simulations alike, the wing was flat. We did not have an objective to account for the wing deformation. No visually noticeable deformation has been observed during the experiment (also see static bending tests in the Supplementary material 1, section S1). In the simulation, the wing was perfectly rigid. However, there were some minor differences in the position of the bristles. They were soldered on one side of the blade in the experiment, but protruded along the mid-plane of the blade in the simulation. Besides that, in the experiment, the shape of the blade near the root was modified for the anchorage mechanism, and the root sections of the model were relatively wide to provide sufficient structural stiffness.
An equivalent membrane wing was fabricated from the bristled model by gluing adhesive tape sheets on either side and cutting the membrane out along the line connecting the bristle tips. Thus, the membrane only covers gaps between long bristles. The shape of the proximal part of the blade remains unchanged. The membrane wing outline is shown in Fig. 1b. This equivalent membrane wing shape definition is consistent with previous experiments (Sunada et al. 2002), but it is not the only one possible. The reader should keep in mind that the equivalent membrane wing is merely an analytical construct providing an intuitive reference for physical interpretation of the bristled wing data. It shows what happens when the gaps between bristles are perfectly sealed.
Since membrane wings are more common than bristled wings, they have received more attention in the previous research. In particular, it became customary to define the aerodynamic force coefficients using as reference quantities the projected area of the wing and the velocity at the radius of gyration determined from the second moment of area (Ellington 1984a, b), which is straightforward for a membrane wing. For a bristled wing, however, the projected area is extremely small, and using it for the force normalization can perplex the comparison with membrane wings (Sunada et al. 2002). In this study, therefore, we always use the geometrical parameters of the equivalent membrane wing non-dimensionalization. Thus, the reference area is equal to the projected area of the membrane wing, S ref = S membrane = 0.52R 2 . The mean chord length is then equal to c mean = S ref ∕R = 0.52R , hence the aspect ratio is as small as R 2 ∕S ref = 1.9 (cf. Chen et al. (2018): for a fruit fly, 3.64 for a bumblebee, 2.78 for a hawkmoth). The radius of the second moment of area (geometric radius of gyration) is calculated using the same outline as shown in Fig. 1b, yielding r g = 0.63R . For comparison, the projected area of the bristled wing is equal to S bristled = 0.08R 2 (i.e., as small as 0.15S membrane ), and 27% of it belong to the blade.

Kinematic protocol
The spatial orientation of a rigid wing model is commonly described using three Euler angles. In the present setup, one of the angles-the one that describes the deviation of the spanwise axis from the horizontal plane-is identically equal to zero. Then, let denote the angle of attack, which in our setup is the angle between the wing and the horizontal plane. Note that, in the present rotating wing setup, the geometrical angle of attack (i.e., relative to the horizontal plane) and the kinematic angle of attack (i.e., relative to the wing inflow direction) are equal. To vary , the wing is rotated about its longitudinal axis conventionally defined as a line connecting the root and one of the most distal bristle tips, see Fig. 1b. In every test, is set before the beginning of motion, and remains constant during the test. Our tests are at = 0 • , 15 • , 30 • , 45 • , 60 • , 75 • and 90 • .
The only angle that varies in time is the positional angle , which determines the rotation with respect to the vertical axis. The motion starts from rest at t = 0 , = 0 and the wing accelerates until it reaches = ∕8 , then continues rotation with its ultimate constant angular speed = 2 f . In the experiment, we tested at five different values of the frequency of rotation: f = 0.04 rps , 0.08 rps, 0.12 rps, 0.16 rps, 0.2 rps, 0.4 rps, 0.7 rps. T h i s corresponds to the range of the angular speed = 0.25, … , 4.40 s −1 . The Reynolds number based on the wing length and the wing tip speed takes the values Re R = R 2 ∕ = 6.0 , 12.1, 18.1, 24.2, 30.2, 60.4 and 105.7. Perhaps a more commonly used in studies on animal flight definition of the Reynolds number is based on the mean chord length c mean and the circumferential velocity of the wing at the radius of gyration r g , yielding Re = r g c mean ∕ = 2.0 , 4.0, 5.9, 7.9, 9.9, 19.8 and 34.6, respectively. We will mainly use this definition through the paper, unless otherwise is stated. In the numerical simulations, we considered Re = 2.0 , 9.9 and 39.6. See Section 4.1.2 in Shyy et al. (2007) for a discussion of the Reynolds number definitions in flapping-wing aerodynamics.

Page 5 of 13 194
The time profiles of the instantaneous angular speed (t) are prescribed as where t is the time from start in seconds. The positional angle (t) varies in time accordingly, in radians, Supplementary material 2 contains a video that illustrates the wing motion.

Fluid-dynamic force normalization
The reference wing speed is obtained as U g = r g = 0.63 R . The force coefficients are determined as where L (lift) is the force in the vertical direction and D (drag) is the force in the direction opposite to the instantaneous velocity of the wing tip (apex).

Rotation mechanism
The dynamically scaled mechanical model of the wing (Fig. 1b) is mounted on a support holder as shown in Fig. 2a. An adjustable anchorage allows to set the geometrical angle of attack to any fixed value from 0 • to 90 • with ± 0.5 • precision, using a digital angle finder. The support holder can rotate about the vertical axis, which passes through the wing root, and the rotation is driven by a NEMA 17 stepper motor 42HS60-1704A (Changzhou Jinsanshi Mechatronics Co., Ltd., Jiangsu, China), transmission belt and a gearbox with the transmission ratio that can be set as 2.5:1 or 12.5:1. The stepper motor is controlled using a TB6560 V2 driver connected to an Arduino Uno R3 controller, which enables gradual constant acceleration in the beginning of rotation. The control program is composed using the AccelStepper library and it allows prescribing the desired angular speed and acceleration of the stepper motor.

Fluid conditions
The wing model is fully immersed in the water-glycerin solution (which is a Newtonian fluid) filling a 50 cm × 80 cm × 25 cm rectangular container (aquarium), see Fig. 2b. The volume fraction of water in the solution is small, and constant temperature is maintained to within ± 0.1 • C from the target value corresponding to the kinematic viscosity of 360 mm 2 s −1 , calibrated using a capillary viscometer (VPJ-2, Ecroskhim, Saint Petersburg, Russia).

Force measurement
The fluid-dynamic forces exerted on the model are measured using a strain gauge load cell (Beijing XNQ Electric Co., Ltd., Beijing, China; for specifications see Supplementary material 1, section S2) with 0.0098 N accuracy. The sampling rate of the force measurement is 1000 Hz . The signal is processed using a custom-built amplifier, a resistor-capacitor circuit low-pass filer, and an L-Card E20-10 analog-to-digital converter (L-Card Ltd., Moscow, Russia). In addition, low-pass biquadratic digital filtering at 5 Hz is applied for denoising. The drag and the lift are measured, respectively, in two different sets of experiments, using the load cell oriented as shown in Fig. 2a to measure the drag, or rotated by 90 • about the longitudinal axis to measure the lift. Each time series obtained for the model at any given and Re is the average of three repetitions. In addition to that, calibration runs have been performed at equal Re, but with no wing model attached. To find the time evolution of the fluid dynamic force, calibration load cell voltage signal is then subtracted from the force measurement signals.

Navier-Stokes solver
The task of setting up a numerical simulation presents difficulties that are unique to bristled wings. The bristles are extremely thin in comparison with the wing length. A direct numerical simulation requires to resolve the fluid motion on these two different scales. This motivates us to employ an adaptive method. We use a wavelet-based incompressible Navier-Stokes solver for fully adaptive computations in time-varying geometries (Engels et al. 2020). It is an opensource software (https ://githu b.com/adapt ive-cfd/WABBI T) optimized for distributed-memory computer systems. The solver is based on Cartesian multi-block decomposition of the computational domain. Spatial derivatives in the governing equations are evaluated on locally uniform Cartesian grid blocks using second-order finite-difference schemes. The grid is adapted by adding or removing blocks of grid points, depending on the magnitude of wavelet coefficients used as refinement indicators. A Runge-Kutta-Chebyshev time marching scheme (Verwer et al. 2004) is employed, which allows explicit integration of all terms in the momentum equation with a reasonably large time step. The velocity-pressure coupling is enforced through the artificial compressibility method. No-slip boundaries are modelled using the volume penalization method.

Flow configuration in numerical simulations
The computational domain is an 8R × 8R × 8R periodic cube, it contains the fluid and the wing model excluding the attachment base and driving mechanism. The wing rotates about the vertical central axis of the domain. The computational domain is block-wise Cartesian. Each block contains 23 × 23 × 23 grid points. The grid is dynamically adapted to the solution, with the limitation of maximum nine levels of refinement, which determines the minimum grid spacing as Δx min = 0.00071R . The threshold value for thresholding wavelet coefficients is fixed as = 10 −3 . The time step Δt is adapted according to the CFL condition with the Courant number equal to 1. The artificial speed of sound is set as c 0 = 30.38 R and the volume penalization parameter is C = 7.82 × 10 −6 −1 .

Results and discussion
The experiment and the CFD results are presented together in this section in a topic-oriented manner focusing on different physical effects. Figure 3 shows an example time profile of the two measured components of the force. This particular case corresponds to the angle of attack of = 60 • , and the Reynolds number is Re = 9.9 , but the time evolution is similar in all cases (see the full data set in Supplementary material 3 and 4). Our choice of = 60 • for illustrative purposes instead of a more conventional 45 • angle is inspired by the large kinematic angle of attack during downstroke found in the smallest bristled-wing insects [e.g., thrips Frankliniella occidentalis (Lyu et al. 2019)].

Time evolution of the forces
The process begins with the acceleration reaction producing a peak of the force. As the acceleration ends, the force relaxes to what we call a quasi-steady value. The transient is slower at the higher Re than at the lower Re but, in all cases considered in this study, the force is nominally constant within the interval of t ∈ [ ∕3, ∕2] . The corresponding time interval is gray-shaded in Fig. 3. We exclude the later part of the sequence when the wing first . The gray shaded rectangle shows the time averaging interval for calculation of the quasi-steady forces encounters the wake of the rotating support holder and then its own wake. We are not interested neither in the transient effects in this study, which may be complicated by the inertia of the mechanical device and wake interactions. Even though the numerical simulation shows similar peaks as in the experiment, these peaks ultimately have little relevance to the forces on flapping wings of insects. Note that the large lift in the experiment during the first 0.2 s of the process is an artifact of low-pass filtering at 5 Hz. Also note that the forces in the numerical simulation show no oscillation after the angular acceleration discontinuity. Hereafter, let us focus on the quasi-steady values of the forces that can offer a crude approximation for a flapping wing in the mid downstroke and upstroke.

Time-average quasi-steady forces
The average force coefficients over the interval t ∈ [ ∕3, ∕2] are displayed in Fig. 4 for a selected combination of Re and . Similar plots for all regimes realized in the experiment are provided in sections S4 and S5 of Supplementary material 1. Fig. 4a portrays a typical variation with respect to the angle of attack , while Fig. 4b elucidates the Reynolds number dependence. In the experiment, the bristled and the membrane models, both, show similar variation of C D and C L with and Re. The values for the bristled wing are confirmed, in addition, by the results of the numerical simulation. The discrepancy between the experiment and the simulation is less than 5% of the maximum C D and 9% of the maximum C L at Re = 9.9 . It increases to 14% for C D at Re = 2.0 , probably because of the wall effect.
In all cases, C L is substantially less than C D , which means that this wing, in the range of Re considered here, is better suited for drag-based flight (force in the opposite direction to the wing motion) than the lift-based (force perpendicular to the direction of wing motion). This is usual for the smallest fliers (Lyu et al. 2019;Jones et al. 2015). Low lift-to-drag ratio ( < 1 ) can be explained by a combination of factors including the low Reynolds number < 100 and the small aspect ratio of the wing < 2 . The difference between the magnitudes of C D and C L increases as Re decreases.
The values of Re in the lower end are sufficiently small to observe transition to a Stokesian regime in which C D and C L are both decreasing functions of Re, which has been reported previously for nominally two-dimensional foils (Jones et al. 2015;Thom and Swart 1940). The values of Re in the higher end are sufficiently large to see the membrane wing C L increase with Re due to the leading-edge vortex, also in agreement with translating wing experiments and numerical simulations (Jones et al. 2015;Miller and Peskin 2004). In contrast, C L of the bristled wing continues to decrease with increasing Re, which means that the bristled wings does not benefit from the leading-edge pressure peak as much as the membrane wing does.
It is instructive to see how large force the bristled wing can generate, in per cent of the membrane wing force. A color map of such quantity, C Dbristled ∕C Dmembrane , is plotted in Fig. 5 versus and Re. We focus on C D here, since we know from the previous discussion that C L is much smaller.
The bristled wing produces less drag than the membrane wing under equal conditions, but no less than 60% . This ratio consistently increases as or Re decrease. At the lowest Reynolds number considered, Re = 2.0 , the ratio is above 90% regardless of . It should be noted that the typical range of Re of P. placentis is supposedly below 20. For example, for a ptiliid beetle Nephanes titan that has the wing length R = 0.66% , which is only slightly larger than P. placentis, kinematic measurements (Yavorskaya et al. 2019) suggest that the flapping frequency is equal to f = 207 Hz and the amplitude is Φ = 195 • . Taking the kinematic viscosity of air at the temperature 25 • C as = 1.56 × 10 −5 m 2 s −1 , we find that Re R varies around 2ΦfR 2 ∕ = 39.3 during one flapping cycle. A wing of N. titan has a similar outline shape to P. placentis, therefore, Re based on r g = 0.63R and c mean = 0.52R varies around 12.9. This means that the bristled wings produce no less than 75% of the force that the equivalent impermeable membrane wings could produce.

Paddling effectiveness factor
The drag-based flight mode in insects (Jones et al. 2015) requires reciprocating motion of the wings: the drag should be larger when the wings move down than when they move up. It can be referred to as the "rowing mechanism" (Cheng and Sun 2018;Lyu et al. 2019), for its similarity with swimming or paddling action that was already noticed in an early work by Horridge (1956).
The present experiment focuses on steady rotation, but these force measurements can be used for quick estimation of the forces acting on flapping wings as well, under the quasi-steady assumption. As the Reynolds number decreases, viscous diffusion of the vorticity becomes more efficient. This means that the flow field adapts faster to the instantaneous position change of the boundary, and the aerodynamic forces depend less on the time history of the flow. Further, if the wing motion is periodic in time, the acceleration reaction makes zero net contribution to the time-average forces. Hence, at low Reynolds number, cycle-average forces can be expressed in terms of the force coefficients from the rotating wing experiment.
To make the quasi-steady calculation simple, let us consider an idealized "vertical stroke" (cf. Jones et al. (2015)). To maximize the net vertical force at a given Re, as the wings move down, the drag coefficient is held at its maximum C Dmax , but when the wings move up, it is kept at its minimum C Dmin . C D is controlled by changing the angle of attack using feathering rotation. The lift coefficient is zero in this scenario, C L = 0 . The real motions of insect wings are, in general, more complex. In the case of P. placentis, precise kinematic reconstruction is yet to be done. Nevertheless, in view of the known trend in the smallest insects to make deep U-shaped strokes to produce large upward pointing drag (Lyu et al. 2019), the simple vertical stroke model, which relies on the drag in a similar way, can provide useful insights. Jones et al. (2015) introduced an aerodynamic performance metric suitable for wing kinematics that predominantly use drag for hovering: the ratio of the net vertical force to the net total force that the muscles have to produce: C V ∕C T , when written in terms of cycle-average force coefficients. Overlines denote time averaging over one wing beat cycle in this discussion. If the left and the right wings move symmetrically, the forces can be taken per wing.
Let us derive an approximate formula quantifying the effectiveness of the vertical stroke gait and analyze its Reynolds number dependence. For simplicity of calculation, we assume a triangle wave time profile of the up an down elevation angle, The feathering angle is piecewise constant in time such that the kinematic angle of attack is equal to 90 • during the downstroke (flat-on) and 0 during the upstroke (edge-on). As Fig. 4a suggests, the maximum drag coefficient C Dmax is attained at = 90 • and the minimum C Dmin is at = 0.
The total force on the wing is entirely due to drag, it acts in the direction perpendicular to the wing, its magnitude can be calculated as F T (t) = 1 2 C D (t)̇2(t)r 2 g S ref , and its vertical component is F V (t) = −F T (t)sigṅ(t) cos (t) . The time averaging yields for the total force perpendicular to the wing, and for the vertical force. A representative approximate value of the flapping amplitude is Φ = 2 (cf. Yavorskaya et al. (2019) for N. titan, 195 ± 4 • in the frontal and 187 ± 3 • in the dorsal projection). We obtain which is a dimensionless quantity that we term the "paddling effectiveness factor" (PEF). It can be taken as a proxy for the efficiency of a wing as a drag-based propulsor, when precise kinematic information is not available. The numerator is approximately the time average useful vertical force carrying the insect aloft, and the denominator is approximately the time average parasite resistance force counteracted by muscles. In the high Reynolds number limit, C Dmin is much smaller than C Dmax , therefore, PEF approaches 2∕ . In the range of Re considered in this study, PEF varies with Re as shown in Fig. 6. It agrees with the general trend for propulsors to lose efficiency at low Re. At Re = 2.0 , the useful vertical force of the flapping propulsor is only about 15% of the total resistance force.
The bristled wing never outperforms the membrane one. At Re = 34.6 , the difference is about 13% . It reduces to 7% as Re decreases to 2.0. It can be speculated from the Stokes flow considerations that this trend can be extrapolated, i.e., PEFs of the bristled wing and of the membrane wing should converge in the limit of Re ≪ 1.

Pressure force and shear force
For membrane wings on the fruit fly scale and above, it is known that pressure forces dominate the shear viscous forces (Roccia et al. 2013). An extreme illustrative example is the fluid-dynamic force acting on a thin flat plate in the normal direction to it. Even at a low Reynolds number, this force is practically due to the surface pressure only, for geometrical reasons.
Another useful idealized example is the Stokes drag on a sphere. In the low Reynolds number limit, the pressure forces account for 1/3 of the total drag on the sphere, and the shear resultant contribution is 2/3 of the total, respectively. Likewise, for a circular cylinder at any Re < 1 , the pressure and the shear stress components are practically equal in magnitude (Dennis and Shimshoni 1965). Even though the bristled wing in our study produces in total almost the same amount of aerodynamic force as the membrane wing, the force-generating structural elements are different. The bristles are circular cylinders, some parts of their surface is perpendicular to the flow direction and some part is parallel. It is, therefore, reasonable to expect the shear force be of the same order of magnitude as the pressure force.
Using the pressure distribution over the surface of the bristled wing model in the numerical simulations at = 60 • , by numerical integration of the pressure gradient multiplied with the volume penalization mask function over the entire computational domain, we calculated the pressure component of the total force in the direction normal to the wing plane, F Np , at time t = 2∕ . The shear component was evaluated as F N = F N − F Np , where F N is the total normal force. Table 1 contains the values of relative contributions F Np ∕F N and F N ∕F N , where the total normal force is related with the lift L and the drag D as F N = L cos + D sin . In the range of Re between 2.0 and 39.6, the pressure force accounts for 44% of the total normal force, and the shear accounts for 56% . This force breakdown reflects the fact that the bristles are bluff bodies. However, it should be interpreted with caution because, while the resolution of 15.14 grid points per blade thickness and 5.46 points per bristle diameter in our numerical simulations is sufficient for accurate evaluation of the total force, the relative contribution of the pressure force may be underestimated by as much as 30%, see Supplementary material 1, section S3. Cheer and Koehl (1987) pointed out that bristled appendages can operate as highly permeable rakes or as virtually impermeable paddles. Which operational mode is realized depends on the combination or morphological parameters such as the bristle diameter and spacing, and on the Reynolds number of the flow. Our force measurements suggest that, in the range of Re considered, the bristled wing of P. placentis can produce aerodynamic forces of the same order of magnitude as membrane wings of the same size. To work as a paddle, the bristled rim must effectively block the air flow through the wing. Let us describe this flow.

Flow field and leakiness
In the following analysis, we consider the wing at a typical rowing angle of attack = 60 • . Let us begin by Paddling effectiveness factor (7) as a function of the Reynolds number depicting the flow velocity component in the direction normal to the wing plane, relative to the wing, at Re = 9.9 . Figure 7a shows its distribution over the wing plane, in the dimensionless form u n ∕ R . The region occupied by the solid structure of the wing is masked with the white color. A black line connects the tips of the bristles, and the region exterior to it is also masked. The color map visualizes the velocity distribution. The velocity is very small, u n ∕ R < 0.1 , close to the proximal side of the wing and near the blade where the bristles are densely packed. As the bristles extend radially, the gap increases and the flow velocity through the gap also increases to u n ∕ R ≈ 0.5 . Further increase of the gap would entail larger throughflow velocity. The flux calculated as an integral of u n over the inter-bristle space Σ evaluates as Dividing it by the "ideal" flux based on the volume across which the wing sweeps we obtain the "leakiness" of the wing (cf. Cheer and Koehl (1987); Kasoju et al. (2018)) which equals 0.24 at Re = 9.9 . As a matter of comparison with earlier related results (Cheer and Koehl 1987), the leakiness is also equal to 0.24 for a pair of 1 μm cylinders with 15 μm spacing at the diameter-based Reynolds number 0.01. These values are close to the typical parameters of setae in P. placentis. shows how effectively the viscous shear stresses block the flow through gaps between the bristles. It depends on Re. The values of obtained using the same procedure for different Re are shown in Table 1. increases with the Reynolds number, as expected (Cheer and Koehl 1987).
An alternative measure to the leakiness is the permeability, as used in porous media flows. The Darcy law relates the fluid velocity u n and the pressure gradient ∇p driving the fluid flow through a permeable continuous medium with permeability k and dynamic viscosity , where p∕ n = ∇p ⋅ , is the flow direction unit vector. For a non-homogeneous anisotropic structure such as a bristled wing, the local permeability can be defined as k = − u n ∕( p∕ n) . The pressure gradient in the entire computational domain was calculated using the same finitedifference scheme and the same grid as in the CFD simulation. Then, it was sampled on the wing plane using linear interpolation in ParaView 5.8.0, and its scalar product with the normal vector of the wing plane was calculated. Thus  obtained spatial distribution of k, normalized by ∕ , is displayed in Fig. 7b. Similar to the throughflow velocity, it increases with the gap between the bristles, but it is less sensitive to the distance from the axis of rotation, because the increase of the velocity with x is compensated by a similar increase in the magnitude of ∇p.
A convenient single-valued criterion of the permeability of the entire wing is the non-dimensionalized average permeability where Σ is the sub-region of the wing plane that includes the inter-bristle space and the interior of the wing, and it defines the virtual membrane area occupied by an equivalent porous medium. This expression evaluates to = 4 × 10 −3 for the data in Fig. 7b. The values for other Re are provided in Table 1. Note that ∫ Σ dxdy∕ ∫ Σ dxdy = 1 − S bristled ∕S membrane = 0.85.

Conclusions and perspectives
The bristled wing model always produces less aerodynamic force than the equivalent membrane model although the difference varies with the Reynolds number. Thus, the drag differs between the two models by less than 10% in the lower end of Re tested, but the discrepancy increases to 40% for the highest Re. The lift shows a similar trend, except that the variation of C L with Re is smaller. The normal force of the bristled wing is evaluated as 44% due to pressure and 56% due to shear. However, our calculations of the those contributions separately are less accurate than the total force calculations.
Like many miniature insects, P. placentis may use dragbased "rowing" to fly. The effectiveness of this aerodynamic mechanism can be estimated knowing the minimum and the maximum drag coefficient. The paddling effectiveness factor obtained from such calculations is lower for the bristled wing than for the membrane one, but the difference is small when Re is in the typical range of the tiniest insects.
Although our study shows no net aerodynamic benefit for an isolated bristled wing model of P. placentis as compared with an equivalent membrane wing, the handicap is small. Positive effect of permeability during clap and fling of a pair of wings (Kasoju et al. 2018) may help equalize the aerodynamic performance. Thus, from the aerodynamic point of view and under the assumptions used in this study, the membrane appears merely unnecessary for flight at this low Re. Advantages of not having a membrane may be sought beyond the aerodynamics. For example, bristled wings are (12) = ∫ Σ kdxdy ∫ Σ dxdy , likely to be much lighter than their membrane counterparts and have lower cost of formation. Our model does not account for the wing deformation, unsteady aerodynamic effects and aerodynamic interaction between two wings, though it is known that these factors can be important [e.g, flexible clap and fling (Miller and Peskin 2009)]. To include them in the model is a promising direction for the future work.
cylinder bar was b base = 5 mm in diameter and 112.5 mm long.
An assembly of six parallel 3D printed models was mounted on a light alluminium frame, attached on a force sensor. The spacing between the bar axes was 82 mm, which corresponds to the typical spacing between setae in the distal part of P. placentis wing (Fig. 8c). During the experiment, the models were dipped in glycerin and the frame slided above the surface over the distance of 200 mm with constant velocity U = 3.37 mm s −1 in the direction perpendicular to the bars. With the kinematic viscosity of glycerin = 360 mm 2 s −1 , the Reynolds number based on the diameter was Re b = b base V∕ = 0.047 . The force sensor signal was recorded similarly as in the rotating wing experiment. The drag was time-averaged over the middle third of the time series.
Similarly, assemblies of cylindrical bars were tested under the same conditions. Circular cylinders were manufactured in an assortment of sizes from the same neutrallybuoyant polymer. The diameter varied from b base to the outer diameter of outgrowth: 5, 7, 9, 10, 11, 12, and 16.5 mm. It was found that the drag of the model with outgrowths is situated between the values corresponding to the 10 and 11 mm cylinder models, and it is significantly different from both. Figure 8d illustrates this point graphically. Thus, the effective diameter is evaluated as b eff = 2.1b base . Supplementary material 7 contains the underlying numerical data and statistical tests.