Leaky Flow through Simplified Physical Models of Bristled Wings of Tiny Insects during Clap and Fling

In contrast to larger flight-capable insects such as hawk moths and fruit flies, miniature flying insects such as thrips show the obligatory use of wing–wing interaction via “clap and fling” during the end of upstroke and start of downstroke. Although fling can augment lift generated during flapping flight at chord-based Reynolds number (Re) of 10 or lower, large drag forces are necessary to clap and fling the wings. In this context, bristles observed in the wings of most tiny insects have been shown to lower drag force generated in clap and fling. However, the fluid dynamic mechanism underlying drag reduction by bristled wings and the impact of bristles on lift generated via clap and fling remain unclear. We used a dynamically scaled robotic model to examine the forces and flow structures generated during clap and fling of: three bristled wing pairs with varying inter-bristle spacing, and a geometrically equivalent solid wing pair. In contrast to the solid wing pair, reverse flow through the gaps between the bristles was observed throughout clap and fling, resulting in: (a) drag reduction; and (b) weaker and diffuse leading edge vortices that lowered lift. Shear layers were formed around the bristles when interacting bristled wing pairs underwent clap and fling motion. These shear layers lowered leakiness of flow through the bristles and minimized loss of lift in bristled wings. Compared to the solid wing, peak drag coefficients were reduced by 50–90% in bristled wings. In contrast, peak lift coefficients of bristled wings were only reduced by 35–60% from those of the solid wing. Our results suggest that the bristled wings can provide unique aerodynamic benefits via increasing lift to drag ratio during clap and fling for Re between 5 and 15.


Introduction
Miniature flying insects with body lengths less than 1 mm, such as thrips, fairyflies, and some parasitoid wasps, have to contend with generating lift in the face of large viscous resistance at chord-based Reynolds number (Re) on the order of 10 or lower [1]. In contrast to our understanding of the aerodynamics of insect flight at the scale of fruit flies and above [2][3][4][5][6][7], the fluid dynamic mechanisms that enable the smallest flying insects to generate lift remain unclear. At low Re∼O(10), drag forces substantially peak and hinder aerodynamic performance [8,9]. Tiny insects show two distinct physical adaptations when compared to larger insects such as fruit flies and hawk moths, including: (1) obligatory use of "clap and fling" wing-wing interaction as a part of free-flight wingbeat kinematics [10][11][12]; and (2) wings consisting of a thin solid membrane with long bristles on the fringes. Recent studies by Santhanakrishnan et al. [13] and Jones et al. [14] have proposed that bristled wings can reduce drag forces needed to fling wings apart at low Re. However, the physical mechanisms underlying drag reduction by bristled wings and the impact of bristles on lift generation via clap and rotation at Re∼10. No clear benefit in aerodynamic performance was observed with bristled wings when compared to geometrically equivalent solid (non-bristlesd) wings. Another recent study by Lee et al. [33] experimentally investigated the aerodynamics and flow structure of comb-like wings with varying gap sizes and angle of attack at Re∼10. Their results showed that comb-like wings generate larger aerodynamic force per unit area in comparison to solid wings. However, both of the above studies [32,33] did not address wing-wing interaction commonly seen in free-flight kinematics of these insects. Santhanakrishnan et al. [13] approximated the wing bristles as a homogeneous porous layer and conducted 2D numerical simulations of clap and fling. Drag and lift forces were lowered by porous wing pairs in both clap and fling, when compared to solid wing pairs. Jones et al. [14] modeled wing bristles in two dimensions as a cylinder array and examined forces generated during single wing in translation and in two-wing fling. Bristles were found to reduce drag forces required to fling the wings apart. In both these studies [13,14], flow through the bristles was not examined.
In this study, we used a dynamically scaled robotic model to comparatively examine aerodynamic force coefficients and flow structures generated by idealized solid and bristled wing physical models during clap and fling. Scaled-up physical models of solid and bristled wing pairs were programmed to move following the kinematics described in a previous study of low Re clap and fling by Miller and Peskin [28]. Chordwise and inter-bristle flows were visualized using 2D particle image velocimetry (PIV). Spanwise PIV measurements were used to examine the proportion of reverse flow through the bristles via quantifying leakiness [34]. Leakiness was used to provide a non-dimensional estimate of volumetric flow that was leaked through the bristle array along the wing span, in the direction opposite to the wing motion (due to flow leaking through the gaps between bristles), as compared to the theoretical maximum volumetric flow that can pass through the inter-bristle gaps with no viscous resistance. Leakiness of zero represents a solid/non-bristled wing, as there is no effective wing surface porosity for fluid to pass through the wing span. Leakiness of one represents the inviscid case wherein reverse flow through the bristles located along the wing span can be obtained as the product of wing tip velocity and the total area occupied by inter-bristle gaps. Time-varying lift and drag forces were measured using strain gauges. Aerodynamic forces were non-dimensionalized to obtain coefficients of lift (C L ) and drag (C D ). Force coefficients for the solid and bristled wing models were compared to the vorticity at the leading edge (LE) and trailing edge (TE) to examine how lift generation was impacted by the presence of bristles. Leakiness obtained from spanwise PIV was used to develop a physical explanation for larger drag reduction by bristled wings (relative to lift reduction) during clap and fling, as compared to solid wings.

Physical Models
In this first experimental study of flow through bristled wing pairs performing clap and fling, we did not examine the diversity in bristled wing structure that is seen among tiny insects. We fabricated and comparatively tested three bristled wing models of varying inter-bristle spacing (herein referred to as gap width between bristles, W) and a solid wing that had no bristles ( Figure 1). Our choice of gap widths were based on a recently published study by Jones et al. [14], in which the authors performed a computational study of 2D bristled wings in fling (bristles modeled as a cylinder array with spacing W). Jones et al. [14] performed morphological analysis of a range of tiny insect wings and determined ratios of gap width between bristles (W) to the bristle diameter. We used the same range of ratios as the models that were examined in Jones et al. [14], but report the values of gap width between bristles divided by wing span (W/S) instead, as our wings were three-dimensional. The gap width (W) between two bristles was given by the relation: where S represents the wing span (90 mm in all models), n represents number of bristles on both sides of the membrane (see Figure 1 for values), and diameter of bristle (D) ≈ 1.3 mm. Four wing models were fabricated for this study (see Figure 1) with: W = 0 (solid wing), W = 4.7 mm (W/S∼0.05), 13.7 mm (W/S∼0. 15), and 21.2 mm (W/S∼0. 23). For the solid wing model, gap width between pair of bristles was zero, such that W/S = 0. Idealized physical models of solid and bristled wings were fabricated using a 3.17 mm thick polycarbonate sheet, such that the model wings had a rectangular cross-section with 90 mm span and 45 mm chord length ( Figure 1). Our model wings were thus nearly two orders of magnitude larger than the typical length scales of tiny insect wings (wing length∼0.5 mm). Bristled wings were fabricated using a 3.17 mm thick polycarbonate membrane (width, c s = 8.5 mm) that was surrounded on either side by cylindrical glass rods of 1.30 mm diameter to represent bristles ( Figure 1). Glass rods were chosen for fabricating bristles so as to provide optical access needed for PIV measurements of inter-bristle flow. The chord length for bristled wings were measured from one bristle tip to the opposing bristle tip on the other side of the solid acrylic membrane (Figure 1).
The wing models ( Figure 1) were tested using the dynamically scaled clap and fling robotic platform ( Figure 2, described below) to obtain measurements of aerodynamic forces and visualization of flow structures. The overall chord and span lengths of the solid and bristled wings were unchanged, and were selected to minimize confinement effects from the walls of the tank. The bristled wings tested in this study were simplified as compared to those seen in tiny insects, and we did not match the following to biologically relevant values: number of bristles (n), percentage of solid membrane area, and relative bristle length (c h /S). (E) The bottom view of a single bristled wing model with directions of measured tangential (F T ) and normal forces (F N ) experienced during rotation by angle α. Lift and drag forces were measured using lift and drag bracket, respectively, by taking components of F T and F N in the vertical (F L ) and horizontal (F D ) directions. F L = lift force; F D = drag force; LE = leading edge; TE = trailing edge; x,y,z are fixed coordinate definitions. The views in (A,E) are referred to in this paper as "spanwise" and "chordwise" views, respectively.

Dynamically Scaled Robotic Model
The model used for replicating wing-wing interaction ( Figure 2) consisted of two wings that were programmed to move symmetrically in opposite directions along a straight line (no change in elevation and stroke plane angles). Translation and rotation of each wing was achieved using two programmable 2-phase hybrid stepper motors with integrated encoders (model ST234E, National Instruments Corporation, Austin, TX, USA). Each pair of rotational and translational stepper motors were rigidly mounted onto an aluminum base plate, which in turn was allowed to slide with minimal friction within two T-slotted extrusions (80/20 Inc., Columbia City, IN, USA) using bearing pads (80/20 Inc., Columbia City, IN, USA). Each wing was mounted onto a stainless steel D-profile shaft of diameter 6.35 mm via a custom-made aluminum L-bracket. Each D-profile shaft was coupled to one rotational stepper motor using a pair of nylon miter gears for transmitting motion from the rotational stepper motor shaft at a 90-degree angle. Each translational stepper motor was coupled to a nylon pinion gear that was in turn coupled to a 0.30 m long nylon rack. The rack for each translational stepper motor was rigidly mounted onto an aluminum bar that was coupled to the T-slotted extrusion. Each of the two D-profile shafts passed through the aluminum base plate and was supported by a steel ball bearing that was press-fitted onto the aluminum base plate. All these motors were controlled using a multi-axis motion controller (PCI-7350, National Instruments Corporation, Austin, TX, USA) via a custom program in LabVIEW software (National Instruments Corporation, Austin, TX, USA). One revolution of each stepper motor was divided into 20,000 steps using a stepper motor drive (SMD-7611, National Instruments Corporation, Austin, TX, USA). The entire assembly was mounted on the top of a 10.6 × 10 −2 m 3 acrylic tank measuring 0.51 m × 0.51 m in cross-section and 0.41 m in height. The tank was filled with 99% glycerin such that the wings were completely immersed in the fluid medium.

Kinematics
The focus of this study was to quantify forces and flow structures generated during clap and fling of bristled wings. At present, there are no published data on clap and fling wingbeat kinematics in tiny insects that describes chordwise rotation and translation in time. This is primarily due to the incredibly small sizes of these insects (body lengths < 1 mm, chord∼0.5 mm). Even with modern high-speed cameras, resolving chordwise rotation and translation at the end of upstroke and downstroke is challenged by lack of control of animal position and orientation during free-flight. As a result, we relied on using prescribed kinematics based on a computational study by Miller and Peskin [28]. These kinematics have been used in other published studies of clap and fling aerodynamics [13,14,30]. The original kinematics developed by Miller and Peskin [28] was intended to provide smooth wing motion throughout the cycle.
The motion profiles (velocity in terms of steps/second) for the stepper motors used in the robotic platform ( Figure 2) were developed based on the wing kinematics used in computational study of clap and fling by Miller and Peskin [28]. The motion profile for the stepper motors (velocity in terms of steps/second) was obtained from the variation of non-dimensional velocity with non-dimensional time (τ) during a clap and fling cycle. Figure 3 shows variation of non-dimensional velocity with non-dimensional time during rotation and translational motion of a wing for two continuous cycles. Non-dimensional time was defined as: where t represents instantaneous time and T represents time taken to complete one cycle (see Table 1 for values of T) that includes one fling half-stroke, stroke reversal, and one clap half-stroke ( Figure 4). Non-dimensional velocity was represented as the ratio of instantaneous wingtip velocity (U) to steady translational velocity (U ST ) of the wing ( Figure 3). U ST was used as the velocity scale for non-dimensionalization of the instantaneous wingtip velocity identical to previous studies of low Re clap and fling aerodynamics [13,14,28,30]. Each wing underwent both rotation and translation, as shown by the two curves in Figure 3. U during wing rotation (U rot ) was calculated as the product of chord (c) and angular velocity of the wing (ω rot ), represented by the thick-lined curve in Figure 3. U during wing translation (U trans ) is shown in the trapezoidal (thin-lined) curve in Figure 3. The motion profiles for the second wing of the wing pair were made identical to the first wing, but driven in opposite directional sense with respect to the first wing. Both wings moved along a straight line without changes to the elevation and stroke angles. Figure 3 shows overlap in translational and rotational motion during both clap and fling half-stroke. While the computational study of clap and fling by Miller and Peskin [28] used symmetric clap and fling strokes, Arora et al. [30] varied the amount of overlap between rotation and translation. In this study, the wings were made to translate along with rotational motion (100% overlap) during clap half-stroke, while the wings started to translate after 50% of wing rotation was completed (50% overlap) during fling half-stroke. We were limited to these asymmetric kinematics between clap and fling due to the design constraints in our robotic platform. Both wings were tethered at their trailing edges to vertical shafts coupled to driving motors ( Figure 2). This did not permit change in physical point of rotation during clap from trailing edge to leading edge. Note that this did not preclude us from realizing effective leading edge rotation of the wing pair during clap-this was achieved via modification of the control program driving the wings.  [28]. Thin line represents non-dimensional linear translational velocity and the thick line indicates non-dimensional rotational velocity. Translational and rotational velocities were both non-dimensionalized by the steady translational velocity of the wing (U ST ). Shaded gray region for τ = 0.8-1.2 indicates the motion profile of a continuous clap and fling cycle where experimental data were acquired. During fling, wing translation started at 50% of wing rotation (50% overlap). During clap, wing translation occurred throughout the wing rotation (100% overlap). Motion profile is shown for only one wing. The motion of the other wing of the pair was identical to this profile but in the opposite direction. Both wings moved along a straight line during both rotation and translation, with no changes to the stroke plane and elevation angles. U = instantaneous wingtip velocity; U during rotation motion is given by U rot and was calculated as the product of chord (c) and angular velocity of the wing (ω rot ); U during translation is given by U trans ; τ = dimensionless time defined by Equation (2).  The motion was started with both wings being parallel to each other and separated by a spacing of roughly 10% chord length, prior to performing fling half-stroke (τ ≈ 0-0.2 in Figures 3 and 4). During fling half-stroke, wings were purely rotated with respect to their trailing edges for 50% of this half-stroke and then translational motion was started with the remaining 50% of rotation (in fling) occurring simultaneously. Once the wings individually rotated to 45 degrees, they were purely translated at 45 degrees angle for about 1 chord length (τ ≈ 0.2-0.4). The wings subsequently underwent stroke reversal (90 degrees rotation in the direction opposite to individual rotation in fling) to position the wings for start of clap motion (τ ≈ 0.4-0.7 in Figures 3 and 4). During clap half-stroke, the wings were made to rotate and translate simultaneously (τ ≈ 0.8-1 in Figures 3 and 4). The maximum angle of rotation for each wing in both clap and fling was 45 degrees, such that the inter-wing angle in between a wing pair in clap and fling never exceeded 90 degrees ( Figure 5). Force and flow measurements reported in this paper were acquired only during rotational portions of clap and fling continuously, starting with clap at the end of 1st cycle to the end of fling of 2nd cycle (τ ≈ 0.8-1.2 in Figures 3 and 4). Table 1. Experimental conditions used for PIV and force measurements comparing solid and bristled wing pairs during clap and fling. Camera recording rates used for time-resolved PIV (TR-PIV) measurements in chordwise plane (FR TR-PIV ) and phase-locked PIV (PL-PIV) measurements in spanwise plane (FR PL-PIV ) were varied to obtain data at identical non-dimensional time points (τ) across every Re tested. Cut-off frequencies used for filtering strain gauge data ( f cut-off ) were varied for every Re. U ST = steady translational velocity of the wing, Re = chord-based Reynolds number given by Equation (3), T = time taken to complete one cycle, which includes one fling half-stroke, stroke reversal, and one clap half-stroke ( Figure 4). As noted above, the degree of separation between the wings used in our study was ∼10% of chord length. There is variability in the spacing between the wings during clap and fling of tiny insects. For example, free-flight recordings of thrips published by Santhanakrishnan et al. [13] show that the wings do not come in full contact but rather have roughly the same separation as considered in this manuscript. Other studies such as that of Kolomenskiy et al. [29] have computationally modeled full contact between the wing pairs in clap and fling. The latter was not attempted in our study. Inter-wing angle for a particular half-stroke is specified both in terms of degrees (left ordinate) and percentage of wing rotation (right ordinate).

A B
In addition to the above simplifications with respect to prescribed wing kinematics, we did not consider complete 3D wing beat kinematics seen in flapping flight of most insects, including: flapping translation, wing supination, and wing pronation [2]. Finally, our prescribed motion does not account for changes to the stroke plane and elevation angle. A horizontal stroke plane is commonly used in hovering studies, and the first paper documenting clap and fling in the tiny wasp E. formosa [10] was based on hovering observations as well.

Test Conditions
The solid and bristled wing models were comparatively tested under varying chord-based Reynolds number (Re) from 5 to 15, defined as: where ν is the kinematic viscosity of the fluid medium, c is the chord length of the wing, and U ST is the steady translational velocity of the wing. U ST was used as the velocity scale in calculating Re, identical to previous studies of low Re clap and fling aerodynamics [13,14,28,30]. U ST was calculated using Equation (3) for a desired Re and motion profiles were created for clap and fling, as shown in Figure 3. A 99% glycerin solution was used as the fluid medium in all the experiments. The kinematic viscosity of 99% glycerin solution used in this study was measured using Cannon-Fenske routine viscometer (size 400, Cannon Instrument Company, State College, PA, USA) to be 790 × 10 −6 m 2 ·s −1 .
The density of 99% glycerin solution (ρ) used in this study was determined to be 1261 kg·m −3 , via the measurement of mass of known volume of the fluid using a standard mass balance (Scout Pro SP401, Ohaus Corporation, Parsippany, NJ, USA). For constant ν and a fixed value of U ST , Re did not change between bristled to solid wing models as the chord length c was maintained constant between them (c = 45 mm, see Figure 1). Re was varied between 5 and 15 in this study by only altering the steady translational velocity U ST (values of U ST are in Table 1). Starting from the non-dimensional velocity versus time profile (Figure 3), different rescaling values were used to vary dimensional U ST to obtain Re in the range of 5-15 (Table 1).

Particle Image Velocimetry (PIV)
PIV was used to visualize the flow structures formed along the chordwise (horizontal plane) and spanwise (vertical plane) directions during the wing motion. Inter-bristle flow visualization data were used to quantify the volume of fluid leaked in between the bristles of the bristled wing models.

PIV Along Chordwise Direction
Two-dimensional time-resolved PIV (2D TR-PIV) was used to visualize chordwise flow field generated during clap and fling of a wing pair at a particular Re ( Figure 6). Hollow glass spheres of 10-micron diameter (110P8, LaVision GmbH, Göttingen, Germany) were used as seeding particles in the fluid medium. Seeding particles were mixed in the fluid medium contained in the aquarium tank at least a day before PIV data acquisition to provide adequate time for settling and homogenous initial distribution. One horizontal PIV plane (HP located at mid-span of the wing, Figure 6) was illuminated using a single cavity Nd:YLF laser (Photonics Industries International, Inc., Bohemia, NY, USA) that provides a 0.5 mm diameter beam of 527 nm in wavelength. A cylindrical lens of 10 mm focal length was used to make a planar laser sheet from the laser beam. A high-speed complementary metal-oxide-semiconductor (CMOS) camera with a spatial resolution of 1280 × 800 pixels, maximum frame rate of 1630 frames/s, and pixel size of 20 microns × 20 microns (Phantom Miro 110, Vision Research Inc., Wayne, NJ, USA) was positioned at the bottom of the tank and focused onto the seeding particles in the plane HP using a 60 mm constant focal length lens (Nikon Micro Nikkor, Nikon Corporation, Tokyo, Japan). The aperture of the camera lens was set to 2.8 for all TR-PIV measurements. Two trigger signals were generated using a custom LabVIEW program (National Instruments Corporation, Austin, TX, USA) each at the beginning of clap and fling half strokes. These trigger signals were provided as an input to the high speed controller unit of the TR-PIV hardware. These triggers enabled the acquisition of TR-PIV data from the starting point for 10 clap and 10 fling cycles. The PIV particle size was in the range of 1.5-3 pixels. Average particle displacements in the measurement plane ranged 4-8 pixels, depending on Re. PIV data were acquired for clap and fling half-stroke across all the conditions tested (varying Re and wing models). For each test condition, 100 raw TR-PIV images were recorder during both clap and fling by varying CMOS camera frame rate based on Re (see Table 1 for values). 2D TR-PIV data from 10 continuous clap and fling cycles were recorded for each experimental condition (Re and wing design).

PIV along Spanwise Direction
Two-dimensional phase-locked particle image velocimetry (2D PL-PIV) was used to visualize flow through the bristles of the physical models during clap and fling. Hollow glass spheres of 10 microns in diameter (110P8, LaVision GmbH, Göttingen, Germany) were used as seeding particles in the fluid medium. A vertical PL-PIV plane located at approximately 10% of chord length away from the leading edge (VP1 in Figure 7) was illuminated using a double-pulsed, single-cavity Nd:YAG laser (Gemini 200-15, New Wave Research, Fremont, CA) with wavelength of 532 nm, maximum repetition rate of 15 Hz and pulse width in the range of 3-5 ns. The laser beam was converted to a 2D planar sheet of thickness 5-6 mm using a 10 mm focal length cylindrical lens. A scientific CMOS (sCMOS) camera with spatial resolution of 2600 × 2200 pixels at a frame rate of 50 frames/s, and a maximum pixel size of 6.5 microns × 6.5 microns (LaVision Inc., Ypsilanti, MI, USA) was used for recording raw PL-PIV image pairs in frame-straddling mode [35] during fling half-stroke. Similarly, another vertical PL-PIV plane located at approximately 10% of chord length away from the trailing edge (VP2 in Figure 7) was used to record raw PL-PIV image pairs during clap half-stroke. Seeding particles in the laser sheet plane were focused using a 60 mm constant focal length lens (Nikon Micro Nikkor, Nikon Corporation, Tokyo, Japan). The aperture of camera lens was set to 2.8 for all PL-PIV measurements. A trigger signal was generated for PL-PIV using a custom LabVIEW program (National Instruments Corporation, Austin, TX, USA) at the beginning of both clap and fling and was provided as an input to the programmable timing unit (PTU) of the PL-PIV hardware. The trigger signal was used as a reference to offset PL-PIV image acquisition to occur at selected phase-locked time points for 10 fling and 10 clap cycles. The PIV particle size was in the range of 1.5-3 pixels. Laser pulse separation intervals between the two images in an image pair (dt) was in the range of 2-9 ms, depending on Re. Average particle displacements ranged 4-8 pixels, depending on Re.
2D PL-PIV data for every Re and wing design test condition were collected at eight different time points of fling half-stroke defined in terms of percentage of fling duration (τ ≈ 1.06, 1.07, 1.09, 1.1, 1.11, 1.13, 1.14, and 1.2 in Figure 3). These time points were selected based on distributing the maximum fling angle into eight equally spaced angular points (12.5%, 25%, 37.5%, 50%, 62.5%, 75%, 87.5%, and 100%). Similarly, 2D PL-PIV data were collected for eight different time points of clap half-stroke defined in terms of percentage of clap duration (τ ≈ 0.8, 0.86, 0.87, 0.89, 0.9, 0.91, 0.93, and 0.94 in Figure 3). These time points were selected based on distributing the maximum clap angle into eight equally spaced angular points (0%, 12.5%, 25%, 37.5%, 50%, 62.5%, 75%, and 87.5%) for every Re and wing design test condition. The angular points were obtained using the clap and fling kinematics during rotational motion ( Figure 5A,B). The tank containing the fluid was placed on top of a rotating turn table with eight equally spaced angular points, in order to collect PL-PIV data along spanwise direction of the wing at the above clap and fling angles without moving camera and laser. For every test condition (Re and wing design), 10 image pairs were acquired (one image pair per clap and one image pair per fling cycle) from 10 continuous clap and fling cycles. Each image pair was phase-locked at one time point in clap or fling cycle. The recording rate of the camera was varied based on Re (Table 1). Inter-bristle flow during fling was visualized for 8 equally spaced angular points ranging 12.5-100% fling time at vertical plane VP1 located at 10% chord length from the LE. Inter-bristle flow during clap was visualized for eight equally spaced angular points ranging 0-87.5% clap time at vertical plane VP2 located at 10% chord length from the TE. The setup was rotated using a turn table (located at the base of the acrylic tank) so that the wing span being illuminated was always normal to the fixed laser sheet and fixed camera focal plane (across every angular position where PL-PIV data were acquired). LE = leading edge; TE = trailing edge; x, y, z are fixed coordinate definitions.

Post Processing
Raw data obtained from chordwise and spanwise PIV recordings were processed using DaVis 8.3.0 software (LaVision GmbH, Göttingen, Germany). The glass rods used in making bristled wings were masked in raw PL-PIV images prior to processing. Multi-pass cross-correlation was performed on the PIV image pairs with an initial window size of 64 × 64 pixels (two passes) and a final window size of 32 × 32 pixels (two passes), each with 50% overlap. Post-processing was performed by rejecting velocity vectors with peak ratio Q less than 1.2, and interpolation was used to replace empty vectors. The processed velocity vector fields were cycle-averaged across 10 fling cycles and 10 clap cycles. 2D velocity components, including (u, v) from chordwise TR-PIV measurements along the x-y plane and (u, w) from spanwise PL-PIV measurements along x-z plane, were obtained following cycle-averaging. Figure 2A shows the coordinate axes used for spanwise PL-PIV and Figure 2E shows the coordinate axes used for chordwise TR-PIV.

Force Measurements
Time-varying forces on solid and bristled wings were measured using uni-axial linear strain gauges of grid size 3.2 mm × 2.5 mm and nominal resistance of 350 ohms (SGD-3/350-LY13, Omega Engineering, Inc., Stamford, CT, USA). Lift (F L ) was defined as the force in the vertical direction while the drag (F D ) was defined as the force along the horizontal direction ( Figure 2E). Two different custom aluminum L-brackets were fabricated for individual/non-simultaneous lift and drag measurements ( Figures 2C,D). Wing models were attached to these L-brackets, which were in turn connected to D-shafts of the robotic model ( Figure 2). The thickness of the L-brackets used for drag and lift measurements were 2 mm and 1 mm, respectively. The variation in thicknesses was necessary to resolve the lower values of lift (in comparison to drag) expected for low Re in the range of 5-15. The thicker drag bracket that was used showed similar response (signal to noise ratio) across all bristled wing models. The response of the drag bracket was evaluated based on maintaining a minimum signal to noise ratio of 5. In addition, we wanted to be consistent by using the same drag bracket for all the wing models. Similar threshold was maintained for lift bracket also across all bristled wing models tested. Strain gauges were mounted on both sides of the L-brackets ( Figure 2). The surfaces of the L-brackets were smoothened and cleaned thoroughly so that the strain gauges could be bonded to surfaces with minimal air gap. Although a wing pair model was driven to mimic wing-wing interaction (Figure 2), lift and drag forces were measured only on one wing of the pair. This was based on the assumption that lift and drag forces would be identical in magnitude in each wing of the wing pair, as the motion profile used for prescribing motion of both wings was symmetric.

Data Acquisition
Raw strain gauge data in the form of analog voltage signals were acquired using a data acquisition board (NI USB-6210, National Instruments Corporation, Austin, TX, USA) via a custom program written in LabVIEW (National Instruments Corporation, Austin, TX, USA). The strain gauge data were acquired at a sampling frequency ( f s ) of 100 kHz at all Re (5-15) for solid and bristled wing models. The number of samples collected during clap and fling for a particular Re can be determined as the product of f s and corresponding clap and fling duration (Table 1). Strain gauge data were acquired for 30 continuous cycles. Prior to acquisition, the setup was driven for 10 complete cycles to establish a periodic steady state condition in the tank. Starting from 11th cycle, data were acquired from the start of clap of that cycle (τ ≈ 0.8-1 in Figures 3 and 4) to end of fling of the next cycle (τ ≈ 1-1.2 in Figures 3 and 4). This enabled us to capture the raw force data only during wing-wing interaction (clap and fling). For every set of strain gauge data, we also recorded the: (a) angular position of the wing using the integrated encoders housed within the stepper motors; and (b) voltage signal before the start of wing motion (when both wings are at rest) for baseline shift/offset correction. The latter raw voltage signal is referred to herein as "zero offset".

Processing
Raw force signals in terms of voltage were processed using a custom script written in MATLAB (The Mathworks, Inc., Natick, MA, USA). The first step in processing was filtering the raw voltage data using a zero-phase delay, third order, low-pass digital Butterworth filter with cutoff frequency varying from 12 to 35 Hz. The cutoff frequency for a particular Re was approximately 10 times of the sampling frequency ( f s ) for that specific Re, the rationale for which was based on a previous study [3]. Table 1 presents the cutoff frequency for filtering at every Re tested. The second step in processing was removing the zero offset from filtered voltage data obtained at the end of first step. The zero offset data was filtered similar to the first step, and then subtracted from the filtered voltage data obtained at end of first step (across 30 cycles of clap and fling). To isolate the aerodynamic forces, we subtracted the inertial forces generated by the wings using the motion profile from recorded data. The inertial forces were obtained from the prescribed motion profile (Figure 3) as the product of wing mass times wing acceleration, and amounted to less than 5% of the measured forces (O(1 mN)). The third step in processing was to use the slope obtained from strain gauge calibration to convert the filtered force signal from Volts to Newtons (N). The fourth step was to calculate non-dimensional coefficients of lift (C L ) and drag (C D ) using equations presented in the next section. The last step in processing was cycle-averaging the calculated non-dimensional lift (C L ) and drag (C D ) coefficients across 30 cycles of clap and fling. We used the Shapiro-Wilk normality test to verify that both force coefficients were normally distributed about the mean (across 30 cycles of clap and fling). The standard deviations of C L and C D were evaluated for each test condition (Re and W/S) about cycle-averaged coefficient values from 30 cycles of clap and fling. For the entire range of Re and W/S tested, the standard deviations (σ) of cycle-averaged C L and C D were less than or equal to 0.16 and 0.9, respectively. This corresponded to the following range of standard errors of the cycle-averaged force coefficients (σ/ √ 30) across all conditions: ±0.03 for C L and ±0.16 for C D .

Vorticity
Out-of-plane vorticity components were calculated from the 2D velocity vector fields in chordwise plane HP (ω z ; Figure 6) and spanwise planes VP1 and VP2 (ω y ; Figure 7) as: where u, v, and w refer to the velocity components along the x, y, and z directions, respectively ( Figure 2A,E).

Lift and Drag Coefficients
From the third step of force data processing, the filtered force signal (after subtraction of filtered zero offset data and inertial force) was converted from Volts to Newtons (N). As shown in Figure 2E, the lift L-bracket ( Figure 2D) measured tangential force (F T ) and drag L-bracket ( Figure 2C) measured normal force (F N ) acting on the wing. Note that each specific L-bracket had different orientation of strain gauge elements so as to effectively sense the normal force acting on the gauges.
Lift (F L ) was defined as the vertical component of tangential force (F T in Figure 2E) and determined using the following relation: Drag (F D ) was defined as the horizontal component of normal force (F N in Figure 2E) and determined using the following relation: where α represents the instantaneous angle of a single wing with the vertical ( Figure 2E). Non-dimensional coefficients of lift (C L ) and drag (C D ) below were used to quantify the forces generated during clap and fling: where U ST represents the steady translational velocity of the wing (values in Table 1). U ST was used as the velocity scale in calculating both lift and drag coefficients, identical to previous studies of low Re clap and fling aerodynamics [28,30]. Using instantaneous wingtip velocity for calculating lift and drag coefficients would render force coefficients to be unrealistically large during periods of wing deceleration (or unrealistically small during periods of wing acceleration). To take into account the effects of the flow through the bristles, an identical "characteristic surface area" (A) was used in calculating C L and C D for both solid and bristled wings. The characteristic surface area (A) was defined as chord c times span S (see Figure 1 for values of c and S).
The following peak ratio was calculated based on maximum values of force coefficients for the various test conditions examined in this study (Re, W/S): where C F refers to a particular coefficient (i.e., C F = C D for drag coefficient, C F = C L for lift coefficient), and the denominator refers to the solid wing (W/S = 0). While peak lift over peak drag ratios were used to examine aerodynamic benefits of bristled wing pairs as compared to solid wings, ∆C F,max was used to characterize the extent of lift and drag reduction by bristled wings as compared to the solid wing.
In addition to the characteristic surface area described above, we also examined calculation of force coefficients using a scaled surface area for the bristled wing (A ) following [32]. This definition removed the areas of gaps between bristles from solid wing characteristic surface area A, and was given by: where c s is the solid membrane width of the bristled wing, S is the wing span, n denotes the number of bristles (n = 30, 12, 8) on both sides of the membrane, D is the diameter of an individual bristle (D = 1.3 mm), and c h is the length of bristle on either side of the solid membrane (see Figure 1). Force coefficients for the bristled wing model that were calculated using A (in place of A) in Equations (8) and (9) are denoted as C L (lift) and C D (drag), respectively.

Leakiness
Depending on Re, inter-bristle spacing, and angle of attack, an array of bristles immersed in an incompressible, highly viscous flow can act as a leaky rake or a solid paddle. Cheer and Koehl [34] defined a non-dimensional term known as leakiness (Le) to characterize the amount of fluid leakage occurring in viscous flow through a cylinder array. Leakiness is the ratio of the volumetric flow rate of fluid that is displaced in the gaps between bristles in the direction opposite to the large-scale viscous flow, to the maximum volumetric flow rate that can be leaked through the inter-bristle spacing under no viscous resistance. Le was calculated for the bristled wing models from spanwise PL-PIV data in vertical planes VP1 and VP2 ( Figure 7B) during fling and clap, respectively, as: where Q viscous is the volumetric flow rate per unit width calculated from the experiments (Figure 8). Q ideal denotes the volumetric flow rate per unit width of leaky, inviscid flow through the same geometry, calculated as: where S denotes the wing span, n denotes the total number of bristles on both sides of the membrane of a bristled wing model, and U tip denotes the tip velocity of the wing in the direction normal to the instantaneous wing position (obtained from rescaling the motion profile in Figure 3). U tip was used because we rotated our setup during 2D PL-PIV measurements using a turn table, so as to quantify the flow through the bristles at different angular positions during clap and fling ( Figure 7) without moving the laser sheet and camera. Therefore, the horizontal component of the 2D PL-PIV velocity vector field (used for volumetric flow rate calculations) is in the direction normal to the wing, i.e., velocity component in the direction of wing motion. U tip was defined as the sum of the rotational (U rot ) velocity and component of translational (U trans ) velocity in the direction normal to the instantaneous wing position. U = U rot for rotation given by product of chord length (c) and angular velocity of the wing (ω rot ) and U = U trans for translation. U tip was calculated from instantaneous rotational (U rot ) and translational (U trans ) velocities as: where α is the instantaneous angle of a single wing with respect to vertical ( Figure 6C).  (14), where U rot represents instantaneous rotational velocity, U trans represents instantaneous translational velocity, and α represents the instantaneous angle made by a single wing with respect to vertical. U tip denotes the wingtip velocity in the direction normal to the instantaneous wing position. The setup was rotated using a turn table located at the base of the tank to acquire PL-PIV data of inter-bristle flow at different positions (% of θ c , θ f ). LE = leading edge; TE = trailing edge.
The viscous volumetric flow rate per unit width (Q viscous ) for a bristled wing model was defined based on the amount of fluid passing through the bristles along the wing span, in the direction opposite to the wing motion. Q viscous was calculated using the equation: (15) where Q solid is the volumetric flow rate per unit width of fluid displaced by a solid wing in the direction of wing motion, and Q bristled is the volumetric flow rate per unit width of fluid displaced by a bristled wing model (with the same span as solid wing) in the direction of wing motion. Both Q solid and Q bristled were calculated along spanwise line L located at a distance of 10% chord length from the right side face of the solid membrane ( Figure 8A), when viewing the wing along the x-z plane ( Figure 2). The length of line L used is equal to the vertical extent of field of view used in PL-PIV which is about 1.5 times wing span. The volumetric flow rate per unit width displaced by solid wing was calculated from experimental data, as opposed to theoretical estimation, to: (a) account for change in volumetric flow rate per unit width along span due to reverse flow structures at the LE and TE (LEV and TEV); and (b) use the same field of view as bristled wing models.
Volumetric flow rate per unit width (Q wing ) was calculated using the equation below in a custom script in MATLAB (The Mathworks, Inc., Natick, MA, USA).
where u denotes the horizontal component of velocity at line "L" (Figure 8A), in the direction normal to the instantaneous wing position.

Average Shear Stress
Average shear stress τ xz was used to quantify the strength of the shear layers formed around the bristles along the span of the wing. τ xz was calculated along spanwise line "L" located at a distance of 10% chord length (identical location to where leakiness Le was calculated; Figure 7A) as follows: where τ xz represents the x-z component (see Figure 7A for coordinate system) of shear stress at point (x,z) along line "L", and µ represents the dynamic viscosity of the fluid medium. u and w are velocities along x and z directions, respectively. x, y, z are fixed coordinate definitions shown in Figure 2A,B,E.  Figure 9). This was in direct contrast to the solid wing, where the LEV and TEV were observed to be concentrated near the leading and trailing edges. With increasing clap angle from 12.5% to 62.5%, both the LEV and TEV were observed to grow in both the solid and bristled wing models. With further increase in clap angle (from 62.5% to 87.5%), the strengths of both the LEV and TEV were observed to decrease. A downward jet was observed from the cavity between the trailing edges.

Chordwise Flow in Clap and Fling: LEV and TEV
The flow generated by wing rotation during fling using the robotic model developed in this study was qualitatively similar to the observations reported by Lighthill [27] and Maxworthy [24]. During fling, as the wings rotated about their trailing edges, two large LEVs formed over each wing at 25% fling angle ( Figure 10). Weak circulating flow with diffuse vorticity was seen at the trailing edges. With increasing fling angle from 50% to 75%, the LEV grew stronger as compared to the TEV in both the solid and bristled wing models. Suction of fluid into the cavity formed between the wings in fling was observed. The onset of wing translation at 50% fling time (τ = 1.1 in Figure 3), equivalent to roughly ∼45% fling angle ( Figure 5B), did not noticeably enhance vorticity near the TE of solid wing until at 100% angle ( Figure 10A). However, the translation motion was found to make the flow through the chord more diffuse for bristled wing models ( Figure 10B-D).
The velocities generated by the bristled wing pair were markedly smaller when compared to the solid wing pair, showing generation of weaker chordwise flow by bristled wings. Throughout clap and fling, both the LEV and TEV were weaker and highly diffuse for the bristled wing when compared to those of the solid wing (row V in Figures 9 and 10). Within the bristled wing models, both the LEV and TEV strengths (proportional to ω z ) were found to decrease with increasing spacing (W/S) between the bristles. Unlike solid wings, vorticity near the TE of the bristled wing model with largest inter-bristle spacing (W/S = 0.23) was close in magnitude to LE vorticity (see row V of Figures 9  and 10). Symmetric LEV-TEV has been previously noted to not be conducive for lift generation at low Re [9].

Force Generation during Clap and Fling
For a solid wing model at Re = 10, the drag coefficient reached maximum value during fling half-stroke ( Figure 11A) at around τ ≈ 1.08 (40% angle in fling half-stroke, see Figure 5B). This is prior to the start of overlapping translational motion of the wings at τ = 1.1 (50% angle in fling half-stroke; Figure 5B), as seen in the motion profile prescribed to drive wing motion (Figure 3). This increase in the drag force could potentially be due to added mass effects at the start of fling half-stroke [36]. In contrast to the solid wing, lower fluctuations were observed in the time-variation of drag coefficients of all bristled wing models ( Figure 11A). This could be due to flow passing through the bristles instead of over the wing. However, for the bristled wing with least spacing (W/S = 0.05), we found the drag force to peak at similar angular position (40% angle in fling half-stroke) to solid wing model, but was smaller in magnitude. Comparing the drag coefficients between solid and bristled wing models ( Figure 11A) in fling at Re = 10, roughly 2-3 times reduction in C D was observed across all the bristled wing models. This suggests that bristles can reduce the forces needed to fling the wings, which is supported by findings of previous computational studies [13,14]. Another point of interest in Figure 11A is a small amount of negative drag for solid wing at τ ≈ 1.14. At this time point, rotational deceleration of the wings occurred simultaneous with the onset of translational acceleration (Figure 3). Rotational deceleration started prior to this time point (τ = 1.1) and decreased drag force as well as TEV strength (proportional to TEV vorticity) across all wing models ( Figure 10). However, the sudden onset of translational acceleration at τ = 1.1 was nearly coincident with negative drag force production in the solid wing. This suggests that added mass effects expected with this sudden onset of translational motion could be responsible for the negative drag. Interestingly, no negative drag force was observed at this time point for any of the bristled wing models. With further progress in fling half-stroke (τ = 1.14-1.2), translational acceleration increased the drag coefficients of all models. This is supported by TR-PIV visualization of chordwise flow (Figure 10), where TEV increased in vorticity at the end of fling due to translational acceleration.
Negative drag force during clap indicates that the force was acting opposite in direction to that of fling. During clap, drag force peaked at about τ = 0.94 ( Figure 11A). This was prior to the wings coming close together, where the fluid between the wings is expected to resist the wings from coming in close proximity. During both clap and fling, peak drag coefficients were found to decrease with increasing inter-bristle spacing (W/S).
In contrast to the drag reduction seen in the bristled wings, lift generated by solid and bristled wings during clap and fling at Re = 10 was roughly in the same order of magnitude ( Figure 11B). The lift force generated was found to reach peak value during clap at around τ = 0.88 (40% angle in clap half-stroke; Figure 5A) for all the wing models. This could be potentially due to wing-wake interaction during stroke reversal (90 degrees rotation before start of clap half-stroke ; τ ≈ 0.4-0.7 in Figures 3 and 4). This is also prior to start of decelerating portion of translational motion at τ = 0.9 (50% angle in clap half-stroke in Figure 5A). During clap, the lift coefficient was observed to decrease from solid to bristled wing model of least spacing (W/S = 0.05) by a small amount. A sudden drop in lift was observed at the end of the clap ( Figure 11B) due to decreased vorticity at the leading edge ( Figure 9). In addition, a small amount of negative lift was observed at the end of clap. This could be due to deceleration of wings and flow induced by the other wing, as both wings come close together. During most of the fling half-stroke, lift coefficient was found to only change slightly from solid to any bristled wing model. A small amount of negative lift was observed during fling for the solid wing at Re = 10 ( Figure 11B). Little to no negative lift was observed in all bristled wing models during fling at Re = 10, showing a modest improvement from the solid wing purely via the inclusion of bristles. Peak value of lift coefficient at Re = 10 was slightly lower for the bristled wing as compared to the solid wing ( Figure 11B). Increasing spacing between the bristles (W/S) did not noticeably affect the lift coefficients, as compared to the impact of changing W/S on the drag coefficients ( Figure 11A).
Peak drag (C D,max ) and lift coefficients (C L,max ) were calculated from the absolute values of drag (C D ) and lift coefficients (C L ) throughout the clap and fling cycle. Peak drag and lift coefficients during clap and fling of the solid wing model were nearly doubled when decreasing Re from 10 to 5 ( Figure 12A,B). This shows the inherent challenge with generating sufficient lift relative to drag at Re∼O(1), corresponding to the lowermost limit of flapping flight by insects [13,14,28]. Within Re range of 5-10, all bristled wing models (W/S = 0.05, 0.15, 0.23) showed an order of magnitude reduction in peak drag coefficient when compared to the solid wing. Looking at only bristled wing models in the Re range of 5-10, large reduction in drag was found with increasing W/S from 0.05 to 0.15. With further increase in spacing between the bristles (W/S = 0.15 to W/S = 0.23), the peak drag coefficients did not noticeably change. A similar trend was observed for peak lift coefficient of bristled wings in Re range of 5-10. Increasing Re from 10 to 15 did not show a marked impact in peak force coefficients of solid and bristled wings. However, the reduction in peak lift was smaller compared to reduction in peak drag ( Figure 12A,B) when comparing solid and bristled wings. Both peak force coefficients decreased with increasing Re from 5 to 15 irrespective of W/S ( Figure 12A,B). Decrease in peak force coefficients with increasing Re was most pronounced for the solid wing ( Figure 12A,B). Overall, the inclusion of bristles rendered the peak force coefficients to only slightly vary with Re in the range of 5-15. The most pronounced impact of the bristles in terms of large drag reduction (over solid wings) was observed in Re range of 5-10, corresponding to the biological Re range of flapping flight in tiny insects that posess bristled wings [14].
To quantify the reduction in peak drag and lift coefficients between a bristled wing pair (non-zero W/S) to the solid wing pair (W/S = 0), we calculated ∆C D,max and ∆C L,max using Equation (10). Across the range of Re tested, bristled wings reduced peak drag coefficient by 50%-90% and peak lift coefficient by 35%-60% when compared to the corresponding peak coefficients of the solid wing model ( Figures 13A,B). Thus, the inclusion of bristles can disproportionally affect lift and drag generated during clap and fling-even for the smallest W/S model examined here (W/S = 0.05). Across Re ranging from 5 to 15, peak lift to peak drag ratios ( Figure 13C) were larger in all the bristled wing models when compared to the solid wing model. While peak lift to peak drag ratio increased with increasing non-dimensional inter-bristle spacing (W/S) from 0.05 to 0.15, the ratio was nearly unaffected when W/S was further increased to 0.23. These results demonstrate that bristled wings can provide aerodynamic benefits over solid wings for Re ranging between 5 and 15 via: (a) augmentation in peak lift to peak drag ratio; and (b) disproportionally larger drag reduction than lift reduction.
From the above results, we observe that ∼55% increase in gap width (W) relative to span (S) from 0.15 to 0.23 produces small changes in the peak force coefficients. To examine the relative importance of even a sparse number of bristles (n = 8 in W/S = 0.23 model) on aerodynamic force generation, we acquired force measurements on a membrane-only model with no bristles at Re = 10.
The membrane-only model can be viewed as the limiting case of W/S > 1. The time-variation of lift and drag forces of the membrane-only model (c s = 8.5 mm and span S = 90 mm as in Figure 1) was observed to follow a similar trend as the other wing models in this study ( Figure A3). For the membrane-only model, peak drag was observed at τ ≈ 1.08 (40% angle in fling half-stroke; Figure 5B) and peak lift occurred at τ ≈ 0.88 (40% angle in clap half-stroke; Figure 5A)-both similar to the other wing models ( Figure A3). However, the magnitudes of both lift and drag forces were lower for the membrane-only model when compared to the sparsest bristled wing model (W/S = 0.23). In addition, the peak lift to peak drag ratio was roughly 53% lower in the membrane-only model (C L,max /C D,max = 0.16) when compared to the sparsest bristled wing of W/S = 0.23 (C L,max /C D,max = 0.3). This suggests that having even a small number of bristles for a given span (as in W/S = 0.23), such that W/S < 1, does indeed impact force generation during clap and fling at Re∼O (10). More experiments are needed on bristled wing models with 0.23 < W/S < 1 to identify the precise ratio of gap width to span (W/S) and number of bristles (n) where peak force coefficients are essentially unchanged. When force coefficients for the bristled wing were calculated using the scaled (or reduced) surface area given by Equation (11), a reduction in bristled wing drag coefficient (C D ) was observed as compared to the solid wing ( Figure A4A) during pure wing rotation in fling (τ ≈ 1-1.1, Figure 5B). This was similar to what was observed when using identical (or characteristic) surface areas between both wing models ( Figure 11A). The scaled drag coefficient of the bristled wing (C D in Figure A4A) increased when translation and rotation started to overlap (τ ≈ 1.1) until τ ≈ 1.12, after which no difference was observed between the bristled and solid wings. However, the scaled lift coefficients of the bristled wing (C L in Figure A4B) were larger than the solid model for most of the fling half-stroke (τ 1.05 in Figure A4). This is in exact contrast to the minor reduction in lift coefficient of the bristled wing when using identical surface areas ( Figure 11B). Examining dimensional lift forces for this particular test condition ( Figure A3B) shows that calculation of lift coefficient using reduced surface area can incorrectly suggest that lift augmentation occurred in bristled wings during fling as compared to the solid wing.

Inter-Bristle Flow during Clap and Fling: Leakiness
To examine the physical mechanism underlying drag reduction by inclusion of bristles, we examined flow fields in spanwise planes VP1 and VP2 ( Figure 7B) via PL-PIV. Horizontal velocity profiles were extracted at line "L", located at a distance of about 10% chord in the same direction as wing motion ( Figure 8A), for θ f = 37.5% and Re = 10 ( Figure 14A,B). Velocities decreased for the bristled wing models with smallest and largest gap widths (W/S = 0.05 and 0.23) when compared to the solid wing, suggesting that flow leaked through the gaps in these bristled wing models. The spanwise distribution of horizontal velocity in the bristled wing model with least spacing (W/S = 0.05) followed the same trend as that of the solid wing ( Figure 14A). In contrast, recirculating flow with both negative and positive velocities were observed in the bristled wing model with the largest spacing (W/S = 0.23; Figure 14B). Fluid leakage through the gaps between bristles was quantified using leakiness (Le) calculated from Equations (12)- (16). Irrespective of W/S, Le showed similar variation across the clap and fling cycle for Re = 10 ( Figure 14C). Le varied between 0.1 and 0.5 in clap and between 0.1 and 0.9 in fling. The time point where peak Le was observed for Re = 10 (θ f = 37.5% in Figure 5B; τ ≈ 1.08) remained invariant with changing spacing between the bristles (W/S).
Despite the ∼55% increase in W/S, Le of W/S = 0.15 and W/S = 0.23 models were quite similar at Re = 10 ( Figure 14C). It is intuitive to expect that the W/S = 0.23 model should permit a larger volume of fluid to leak through the bristles when compared to the W/S = 0.15 model with a smaller gap width (W). However, this expectation does not necessitate that Le of W/S = 0.23 model also be greater than Le of W/S = 0.15 model. As per Equation (12), Le is the ratio of flow rate leaked by a bristled wing model under viscous conditions (Q viscous ) to the theoretical/ideal flow rate that would be leaked by the same wing model in an inviscid flow (Q ideal ). A small change in Le, such as between W/S = 0.15 and 0.23 models at Re = 10 ( Figure 14C), can therefore be interpreted as both models leaking fluid by roughly the same fraction relative to their individual ideal volumetric flow rate.
While Q viscous was calculated from PL-PIV data, Q ideal depends on the bristled wing design and prescribed wing kinematics as per Equation (13). To maintain the same wing span and bristle diameter across all bristled wing models, W was increased in our models by lowering the number of bristles (n) based on Equation (1). The ideal volumetric flow rate, given by Equation (13), would increase with decreasing n. As U tip is identical across all wings for a specific value of Re, it will cancel out when considering ratio of Q ideal between two bristled wing models of different W/S. Thus, Equation (13) can be used to give the following ratios for a specific value of Re: (Q ideal | W/S=0.23 )/(Q ideal | W/S=0.15 ) = 84.8/82.2 = 1.03 (3% increase in Q ideal ). In contrast, (Q ideal | W/S=0.15 )/(Q ideal | W/S=0.05 ) = 82.2/70.5 = 1.17 (17% increase in Q ideal ); and (Q ideal | W/S=0.23 )/(Q ideal | W/S=0.05 ) = 84.8/70.5 = 1.20 (20% increase in Q ideal ). As Le is inversely related to Q ideal as per Equation (12), Le of W/S = 0.15 and 0.23 bristled wing models need not show as much of a difference as compared to the differences in Le of: W/S = 0.05 and 0.15 models, or W/S = 0.05 and 0.23 models. When W/S was increased from 0.15 to 0.23, the expected increase in Q viscous (due to larger inter-bristle spacing) in W/S = 0.23 was most likely countered by a small increase in Q ideal in W/S = 0.23-resulting in only a minimal increase in Le as per Equation (12). In this context, it is important to note that the definition of Le used in this study was originally proposed by Cheer and Koehl [34] to quantify leakage flow through an array of two cylinders. Despite use of the same definition to quantify Le through arrays of multiple cylinders (as in this study), the close similarity in Le between W/S = 0.15 and W/S = 0.23 wing models was also mirrored by the similarity of peak force coefficients of the same two wing models (Figure 12). An examination of the inter-bristle vorticity and velocity fields during clap and fling (Figures 15-17) shows that high-shear regions are formed around the bristles at Re = 10. For all bristled wing models, shear layers around bristles was observed and was of significant value throughout the clap half-stroke. During fling, the shear layers around bristles were found to be minimum for θ f < 37.5%. The vorticity of these shear layers increased until θ f = 62.5% and decreased for θ f > 62.5% of fling half-stroke. This was coincident with the time points where peak leakiness (θ f = 37.5% fling angle, τ ≈ 1.08) and minimum leakiness (θ f = 62.5% fling angle; τ ≈ 1.12) were observed. This suggests that rotation and translation can impact leakiness in opposite manner, such that pure rotation in fling increases leakiness while onset of translation decreases leakiness and aids in increasing vorticity of shear flow around the bristles. With increase in spacing between the bristles (W/S) from 0.05 to 0.23, the vorticity of shear layers formed around the bristles was found to decrease during both clap and fling for all θ c and θ f (Figures 15-17). For bristled wing model of W/S = 0.23, we observed the flow to recirculate around the bristles for θ f < 62.5% fling angle (Figures 17 and A5C), while this phenomenon was absent during clap. In addition, the magnified inter-bristle flow field for W/S = 0.15 at θ f = 37.5% ( Figure A5B) showed the onset of small amount of recirculation around each bristle. However, the regions of recirculating flow were found to be quite small for W/S = 0.15 due to large shear layers ( Figure A5B) formed around neighboring bristles. These shear layers could further reduce the effective gaps between bristles and thereby limit the extent of recirculation and reverse flow through the bristles. For W/S = 0.05, larger shear layers without any recirculation were observed around the bristles when compared to W/S = 0.15 and 0.23 ( Figure A5A). These factors could potentially contribute to lower leakiness in W/S = 0.05 as compared to other two bristled wing models of larger W/S.A larger gap between adjacent bristles (as in W/S = 0.23) would favor the development of recirculation regions around bristles. However, increasing the gap to such an extent that there is no interaction between shear layers of adjacent bristles may not favor the development of recirculation. Comparing the W/S = 0.23 model with the limiting case of membrane-only wing model with no bristles (Figure A3), we can observe that interactions between the shear layers formed around adjacent bristles can have a non-negligible impact on forces generated during clap and fling. In the case of W/S = 0.23 model, it appears that the gap width was large enough for more volume of flow to leak through than W/S = 0.05, but not so large so that shear layers of adjacent bristles in W/S = 0.23 did not interact. If the latter were true, we would not expect: (a) recirculation of flow around bristles in W/S = 0.23 model; and (b) there to be any appreciable variation in forces generated by W/S = 0.23 and membrane-only models (contradicted by Figure A3). Although W/S = 0.05 favors the formation of high-shear regions around the bristles on account of the small gap width, there is not adequate gap in between interacting shear layers to permit flow to recirculate. To quantify the strength of the shear layers formed around the bristles, we calculated average shear stress τ xz using Equation (17) at the same location "L" where Le was determined ( Figure 8A). Irrespective of W/S, average shear stress was larger in clap half-stroke as compared to fling at Re = 10 ( Figure 14D). The average shear stress decreased with increasing W/S. The average shear stress for W/S = 0.15 and 0.23 models showed only a small amount of change, similar to Le for these two wing models ( Figure 14C).
Comparing Figure 14C,D for a particular bristled wing model, it is clear that leakiness is inversely related to the average shear stress. This reinforces our previous suggestion that interacting high-shear regions formed around the bristles (e.g., W/S = 0.05 in Figures 15 and A5A) can generate "inter-bristle blockage", lowering the effective gap through which fluid can leak through, which consequently lowers Le. The recirculating flow around the bristles for W/S = 0.15 and 0.23 wing models (Figures 16  and 17) started to dissipate after θ f > 37.5% fling angle ( τ ≈ 1.08). This time point was just prior to start of translation motion during fling ( τ ≈ 1.1; Figure 5B). The shear layers around the bristles were found to increase in strength with the onset of translation motion, as evidenced by an increase in τ xz for τ ≥ 1.1 in Figure 14D. This was particularly evident during clap (100% overlap of translation and rotation), where the shear layers around the bristles were stronger than during fling ( Figure 14D). The high-shear flow around the bristles thus directly contribute d to decrease in Le. While pure rotation in fling increased Le due to lower shear stress near the wing bristles (τ= 1-1.1 in Figure 14C,D), combining translational motion with rotation in clap increased shear stress and lowered Le (τ = 0.8-1.0 in Figure 14C,D).

Discussion
In this study of clap and fling wing-wing interaction at low Re ranging from 5 to 15, we observed that: (a) lift and drag coefficients were both lowered in the case of bristled wings when compared to those of solid wings (for most of the cycle); (b) with increase in W/S of bristled wings beyond 0.12, very little change was observed in lift and drag coefficients; and (c) peak drag reduction in a bristled wing as compared to a solid wing was disproportionally larger than peak lift reduction between bristled and solid wings. LEV and TEV were both smaller and diffuse in all the bristled wings tested as compared to the solid wing, during both clap and fling half-strokes. This lowered the lift generated by bristled wing models. Leaky inter-bristle flow in the direction opposite to wing motion was observed in all bristled wing models throughout clap and fling. This inter-bristle flow through the bristles primarily reduced drag forces needed to perform clap and fling. From the above findings, we propose that bristled wings can provide unique aerodynamic benefits via improved lift to drag ratios during clap and fling at Re∼O(10) as compared to solid wings.
Drag reduction is more significant when going from a solid wing to a bristled wing due to leaky, reverse flow through the bristles throughout clap and fling (see . This is expected to lower the resistance offered by the fluid to the wing motion. Leaky flow through the bristled wing models in fling was accompanied by a reduction in the average shear stress near the wing span ( Figure 14D), which in turn lowered the viscous drag experienced by bristled wing pairs during fling. Lift generated during rotation in fling is associated with the bound circulation of the leading edge vortices and asymmetric reduction of trailing edge vorticity (see Ellington [15] and Sane [2] for detailed explanation). From Figures 9 and 10, we see a reduction in leading edge vorticity in all the bristled wing models as compared to the solid wing model. This lowered the lift generated by bristled wings in clap and fling relative to the solid wing ( Figure 13B). From inter-bristle flow visualization (Figures 15-17), we observed the formation of shear layers around the bristles throughout clap and fling. The shear layers around adjacent bristles reduce the effective gap through which fluid leakage occurs, which aids in reducing the loss of lift in bristled wings. We also expect that the shear layers around the bristles would interact with leading edge vortices during fling, causing the vorticity to diffuse across a larger portion of the wing (as seen in Figures 9 and 10) and potentially also help in minimizing loss of lift.

Implications for Aerodynamic Performance
Jones et al. [14] conducted image analysis of a large variety of tiny insect wings and found that the non-dimensional inter-bristle spacing (based on bristle diameter) was in the range of 4-12. Our bristled wing physical models with W/S = 0.05 and 0.15 fall in the above biological range (W/D ≈ 4 and 11, respectively). We also included W/S = 0.23 (W/D ≈ 16) to examine how increasing inter-bristle spacing beyond the biologically observed range would impact lift and drag. Across the Re range 5-15 that was examined, all the bristled wing models showed larger ratio of peak lift to peak drag relative to the solid wing model ( Figure 13C). Interestingly, the ratio of peak lift to peak drag increased noticeably when increasing W/S from 0.05 to 0.15, but this ratio was not nearly not impacted with further increase in W/S from 0.15 to 0.23 ( Figure 13C). Moreover, the lowest lift coefficients were also generated by the same bristled wing (W/S = 0.23; Figure 12B). This shows that the increasing inter-bristle spacing beyond W/S = 0.15 model only minimally improves aerodynamic benefit in terms of peak lift relative to peak drag, but at the cost of further reduction in lift. The latter should be particularly viewed as a compromise to aerodynamic performance, as a minimum amount of lift is necessary for takeoff and supporting body weight. This could be one of the reasons why spacing between the bristles in tiny insect wings (W/D) is in the range of 4-11, as bristled wings with larger inter-bristle spacing (such as the W/S = 0.23 model in this study) may not generate sufficient lift force to support their body weight. The extent of variation in the number of bristles per unit span of wings in tiny flying insects warrants further investigation, especially to understand if the inter-bristle spacing is closely tailored with a particular number of bristles (or a narrow range of numbers) to obtain the most drag reduction for the least amount of loss in lift.
An accurate assessment of the lift requirements of tiny insects is difficult due to a scarcity of body mass and free-flight data. Similar to inter-bristle spacing, it is possible that there is a large variation in the weight of tiny insects and thus in the minimum lift force required by these insects. To understand if the loss of lift observed in our bristled wing models is indeed a compromise to achieve larger drag reduction and a subsequent increase in lift to drag ratio, it is instructive to look at an estimate of minimum lift coefficient needed for supporting the body weight of a tiny insect during hovering. We considered the following variables from published measurements on thrips for calculating C L : body mass ≈ 10 µg [37]; chord (c) = 0.3 mm [32]; span (S) = 0.81 mm [32]; wing tip velocity, U tip ≈ 0.5 m·s −1 [13]. Using Equation (8) with weight in the numerator, U tip in place of U ST , ρ air = 1.225 kg·m −3 , and a characteristic surface area (A = c·S) of 0.24 mm 2 , we obtain C L ∼O(1). The peak lift coefficient of our solid wing model was in this order of magnitude for Re between 5 and 15 ( Figure 12B). However, peak C L of the bristled wing model with W/S = 0.05 was reduced by 45% from the peak C L of the solid wing (∆C L,max at Re = 10, Figure 13B). This suggests there may indeed be a compromise to lift generated by bristled wings versus the lift required to support body weight. Further increase in inter-bristle spacing to W/S = 0.15 and 0.23 resulted in about 60% reduction of peak C L from the peak C L of the solid wing (∆C L,max , Figure 13B). However, this compromise in lift generated by bristled wings should not be viewed in isolation, as the drag reduction realized by bristled wings during fling is quite large (50-90% reduction in C D compared to solid wing; Figure 13A). This is particularly noteworthy when considering that the peak drag coefficient of the solid wing model at Re = 10 was more than twice of the bristled wing with W/S = 0.05 ( Figure 12A). Given the obligatory use of clap and fling among miniature flying insects, drag reduction benefits of bristled wings can thus be quite substantial.

Comparison with Previous Studies
An experimental study by Sunada et al. [32] found that forces acting on single bristled wings undergoing translation and rotation were lower when compared to single solid wing models. This was also supported by our finding in wing pairs. We validated the force measurement methodology used in our study through comparison of published data in Sunada et al. [32] for the case of a single wing in constant velocity translation (see Appendix A for details). A recent 2D computational study by Jones et al. [14] on bristled wings noted that there was substantial decrease in force required to fling the wings apart at low Re. From our results, we observed that drag decreased significantly for bristled wings for Re from 5 to 15. When compared to single wing translation, Santhanakrishnan et al. [13] showed via 2D numerical simulations that larger drag reduction was seen during fling of a bristled wing pair as opposed to lift reduction during fling of a wing pair at low Re. This is also supported by our finding that bristled wings can substantially reduce drag in fling without a proportional decrease in lift at Re∼O(10). Arora et al. [30] found via numerical simulations that for symmetrical clap and fling motion of a solid wing pair, increasing the overlap between the translational and rotational phases resulted in significant increase in lift and drag forces. For any particular overlap, peak lift was observed during clap half-stroke and peak drag was observed during fling. In our study, we had asymmetry in overlap between clap and fling motion profile (100% overlap of translation with rotation during clap; 50% overlap during fling). From our results, we also observed peak lift during clap half-stroke and peak drag during fling half-stroke.

Inter-Bristle Flow and Leakiness
To examine the fluid dynamic mechanism underlying the large drag reduction by bristled wings, we determined leakiness (Le) of flow through the bristles from phase-locked PIV data along the wing span. Le was found to vary throughout the entire clap and fling cycle and reach the peak value for Re > 8 at τ ≈ 1.08 and found to increase with increasing spacing between the bristles (W/S).
An observation of the inter-bristle vorticity fields (Figures 15-17) showed the formation of small-scale regions of high-shear flow around the bristles at all time points during clap and fling. A recent study by Lee et al. [33] comparing single solid and bristled wings in translation, observed these shear layers around the bristles during pure translation at different angles of attack. They suggested that these shear layers would in fact block the gaps between the bristles, making bristled wing models to behave similar to a solid wing model and producing nearly same amount of lift force as that generated by a solid wing model. Irrespective of bristled wing W/S ratio, we observed that the reduction of peak drag coefficient (C D,max ) by bristled wings (when compared to the solid wing) was larger than reduction of peak lift coefficient (C L,max ) relative to the solid wing ( Figure 12). Across the range of Re examined here, C D,max of bristled wings were reduced by 60-90% of the solid wing ( Figure 13A). In contrast, C L,max of bristled wings were reduced by only 35-60% ( Figure 13B). Thus, bristled wing pairs interacting via clap and fling can provide smaller reduction in lift forces as compared to a larger reduction in drag forces.The disproportionally larger drag reduction relative to lift reduction resulted in increasing peak lift to peak drag ratios for bristled wings (all W/S) beyond those of the solid wing for Re ranging from 5 to 15 ( Figure 13C).
In general, leakiness was observed to peak at τ ≈ 1.08 across Re and W/S conditions examined in this study (θ f = 37.5% fling angle; see Figure 5D). Peak leakiness (Le max ) for a particular W/S increased with increase in Re from 5 to 10, but did not appreciably change with further increase in Re from 10 to Re = 15 ( Figure 18). The bristled wing model with least inter-bristle spacing (W/S = 0.05) had the lowest value of Le max across the entire range of Re examined. While increasing W/S from 0.05 to 0.15 resulted in increasing Le max across all Re, further increase of W/S from 0.15 to 0.23 did not appreciably alter Le max (Figure 18). The latter trend mirrored the relative insensitivity of peak force coefficients when increasing W/S from 0.15 to 0.23 ( Figure 12).
To understand how aerodynamic force coefficients were related to leakiness, we examined C L and C D at the time point where Le max was observed ( Figure 18). Compared to Le max , drag coefficient at peak leakiness (C D | Le max ) showed an opposite trend with respect to its variation with Re. Irrespective of Re, C D | Le max decreased with increase in inter-bristle spacing (W/S) from 0.05 to 0.15, and to a lesser extent with further increase in W/S from 0.15 to 0.23. Thus, drag coefficient appears to be inversely related to leakiness (similar to how average shear stress was related to Le in Figure 14D). This observation supports our finding in this study, as well as those of previous studies [13,14], that bristled wings can reduce the drag force needed to fling wings apart. In direct contrast to C D | Le max , lift coefficient at peak leakiness (C L | Le max ) showed relatively less variation for the range of W/S and Re examined here. Further, C L | Le max was roughly an order of magnitude smaller than C D | Le max , suggesting that Le had a smaller impact on lift forces as compared to drag forces. From these results, we suggest that Le plays a larger role in the reduction of drag forces by bristled wings as compared to lift reduction.
An observation of the inter-bristle flow field at peak Le (corresponding to τ ≈ 1.08, θ f = 37.5%; Figure A5) showed the formation of recirculating flow from velocity vectors around the bristles. The recirculation of flow around the bristles reduces the overall volumetric flow rate of fluid swept by the wing, by maximizing the leakiness between the bristles. The additive volume of fluid that is fluxed into the cavity between two bristled wings in clap and fling (whenever leakiness is non-zero) is expected to change the pressure inside the cavity, reducing the drag force required to clap and fling the wings ( Figure 19A). While previous computational studies [13,14] have shown lower drag to be a primary result of clap and fling motion using a bristled wing pair, they did not observe the high-shear flow around the bristles and were unable to clarify the physical mechanism underlying the observed drag reduction (other than conjecturing that flow must be leaked through the bristles).
One interesting observation is the time point at which we see larger reduction in drag force (which is τ ≈ 1.08; see Figure 11A). This is the same time point ( τ ≈ 1.08) where we observed peak Le value for all the bristled wing models in Re range of 5-15. Wing rotation during fling reaches peak velocity value at 37.5% fling angle ( Figure 5D), and thus it is likely that pure rotational acceleration contributes to realizing maximum possible leakiness. For bristled wing models of W/S = 0.15 and 0.23, we observed the flow to recirculate around the bristles for θ f < 62.5% fling angle, while this phenomenon was absent during clap. These recirculating flow started to dissipate after θ f > 37.5% fling angle ( τ ≈ 1.08) for both wing models (W/S = 0.15, 0.23). This time point is just prior to start of translation motion ( τ ≈ 1.1; Figure 5B). The shear layers around the bristles were found to increase in strength with the onset of translation motion. This is particularly evident during clap (100% overlap of translation and rotation), where the shear layers around the bristles were stronger than during fling. The high-shear flow around the bristles with the onset of translational motion thus directly contributes to decrease in Le. The importance of the small-scale, high-shear regions in decreasing Le is also supported by the weaker vorticity seen in the case of wing models with larger inter-bristle spacing (W/S = 0.15, 0.23; Figures 16 and 17). Leaky flow through bristles was lowest at Re = 5 for any bristled wing model ( Figure 18). This suggests that bristles may not be efficacious in drag reduction at Re∼O(1). Bristled appendages have been previously observed to typically act as solid paddles at low velocities, but transition to acting as leaky rakes during moments of higher velocity or interactions with surfaces [34]. The lower tip velocity of the wing can thus limit the peak leakiness achieved for Re = 5.

Considerations for Physical Model Studies of Bristled Wings
The central difficulty with conducting in vivo experimental measurements on tiny insects is needing to resolve wing length scales of O(0.1 mm) when flapping at frequencies of O(100 Hz). Computational and experimental modeling studies can thus be highly useful in examining the underlying flow physics. However, computational fluid dynamics (CFD) studies examining clap and fling aerodynamics with bristled wings have been challenged by the large computational costs associated with modeling bristled wings while simultaneously resolving the large-scale flow field generated around the surface areas of an interacting wing pair. This is because tiny insects can have fifty to hundred or more bristles per wing. To date, bristled wings have been modeled as either: (a) a homogenous porous layer [13]; or (b) 2D array of cylinders [14]. While assumption (a) does not capture the spatial complexity and resolve inter-bristle flows,assumption (b) can only resolve inter-bristle flow per unit span (bristles at wing tip that are laterally distributed about the chord).
To model a conservative number of bristles (say 30-50), 300-500 spatial grid points would be required along one direction and along the length of any one wing. This is based on an assumption that at least 10 spatial grid points would be needed between each bristle to resolve inter-bristle flows. For the case of two wings in a domain large enough to avoid boundary effects, accurate simulations would require spatial grids that are at least on the order of 3000 3 -5000 3 . The time step requirements for most immersed boundary type methods (typically used to simulate these types of problems) are on the order of 10 −4 -10 −6 of a wingbeat. The situation is improved with adaptive mesh refinement and adaptive time stepping. Even with the latter capabilities, only simulations over short times are possible at the time of this study. Physical model based approaches, such as in this paper, can help to address the challenges with in vivo measurements and computational modeling. Further, experimental studies can provide data on more realistic bristled wing models for verification and validation of simplified porous or poroelastic models of bristled wings such as in Santhanakrishnan et al. [13].
We observed added mass effects during fling via abrupt decrease in the drag coefficient of the solid wing to negative values (see Figure 11A). Added mass effects such as these depend on fluid properties, including the density and dynamic viscosity [38]. We matched the Re to biologically relevant range of O(1)-O(10) for scaled-up physical models of bristled and solid wings by using glycerin as the fluid medium. As glycerin has substantially higher density and viscosity than air, added mass effects are expected to be far more prominent in this study as compared to what one would expect in air. This was one of the reasons underlying our choice of presenting non-dimensional force coefficients in the paper (as opposed to dimensional forces). However, we expect that the observed trends in peak force coefficients (Figures 12 and 13) and leakiness ( Figures 14C and 18) with respect to changing Re and W/S would not be impacted by change in fluid medium to air. Thus, the findings of this study are expected to be applicable to clap and fling flight at low Re∼O(10) observed in miniature flying insects such as thrips.
Although this study shows that leakiness can be impacted by changing spacing between the bristles as well as chord-based Reynolds number, it is important to note that alternative definitions of Re will need to be considered for a more complete understanding of flapping flight using bristled wings. Specifically, chord-based Re could be less relevant to use in bristled wing designs with large inter-bristle spacing (large W/S) as opposed to Re based on bristle diameter (Re b ). Jones et al. [14] matched Re b to those observed in tiny insects and reported lower leakiness values than what we observed. Despite not matching Re b , both our results and those of Jones et al. [14] showed decrease in C D with increasing inter-bristle spacing. However, their computational model examined flow through 2D arrays of cylinders in fling and not finite span bristled wings that were considered in this study. Further studies are needed to examine the differences in forces and flow patterns from matching bristle-based and chord-based Re to biologically-relevant values seen in flight of tiny insects.

Use of Reduced Surface Area for Evaluation of Bristled Wing Force Coefficients
When identical wing area was used in calculation of lift and drag coefficients of solid and bristled wings ("characteristic area" A defined as product of chord c and span S), peak coefficients of all the bristled wing models were lower than the corresponding values in solid wing model across the entire range of Re that was examined (Figures 12 and 13). However, this observation was contradicted when scaled (reduced) wing area A (defined as per Equation (11)) was used in the calculation of lift and drag coefficients of bristled wings ( Figure A4). Specifically, Figure A4 shows that peak C L calculated using A (denoted in Figure A4 as C L ) actually exceeded peak C L of the solid wing, while peak C D was slightly lower than peak C D of the solid wing. Comparing the scaled lift coefficients of the bristled wing model (C L ) with the corresponding dimensional lift forces ( Figure A3), it appears that exclusion of the gaps in between the bristles in the scaled surface area Equation (11) can potentially lead to an incorrect observation that lift augmentation in fling was observed in bristled wings as compared to solid wings.
Indeed, the simplest reason for the above discrepancy is purely mathematical, i.e., lowering the surface area in computing force coefficients of bristled wings would numerically increase their values. A more fundamental argument that could explain the above discrepancy is that use of the reduced surface area Equation (11) implicitly assumes that flow through the bristles is mostly identical to the ideal/inviscid case illustrated in Figure 8B, without allowing the possibility that flow can incompletely leak through the gaps between the bristles on account of viscous interactions. Although previous studies of single bristled wings in translation and/or rotation [32,33] have used reduced areas in calculation of force coefficients, we suggest that comparative studies of bristled and solid wings during clap and fling at low Re opt for using a characteristic surface area that is identical across bristled and solid wings. This is especially advisable when the bristled wing designs being compared have lower inter-bristle spacing, where viscous interactions in between adjacent bristles would be non-negligible (as seen in Figures 15-17 and A5).

Interaction of Inter-Bristle Flow with LEV and TEV
We sought to address a limitation of previous computational studies of wing-wing interaction with bristled wings [13,14], where the chordwise flow structures (LEV and TEV) were not examined in conjunction with inter-bristle flow during clap and fling. At Re = 10, we observed the following in fling: (1) LE and TE vorticity were larger for solid wings as compared to bristled wings (see Figure 10A,D); and (2) with increasing W/S, both LE and TE vorticity decreased ( Figure 10A-D) while Le increased ( Figure 14C). We observed similar results during clap (see Figure 9), with bristled wings showing lower vorticity at the LE and TE as compared to those of the solid wing.
The magnitude of vorticity around the individual bristles (ω y ), obtained from PL-PIV measurements (Figures 15-17), was larger than the magnitude of vorticity around the LE and TE of the bristled wing (ω z ; Figures 9 and 10). From chordwise PIV measurements (Figures 9 and 10), we observed both the LEV and TEV to be diffuse in all the bristled wing models as compared to more coherent LEV and TEV in the solid wing model. Interaction of the inter-bristle flow with the LEV and TEV could thus: (a) remove some of the momentum from the LEV and TEV ( Figure 19); and (b) promote diffusion of vorticity at LE and TE of bristled wings during both clap and fling. The lower vorticity of the LEV is expected to reduce the lift force generated in bristled wings, as lift generated in fling is associated with the bound circulation of the leading edge vortices during wing rotation [15]. However, shear layers were observed to form around the bristles for all of the bristled wing models tested in this study (Figures 15-17). These shear layers lower the effective gap between the bristles through which reverse flow can occur, resulting in the bristled wing behaving like a solid wing, as proposed by Lee et al. [33]. Therefore, the generation of shear layers around the bristles (particularly augmented in translation) can aid in minimizing the loss of lift generated by bristled wing models during clap and fling.

Limitations of the Study
We acknowledge the following simplifications used in our study with respect to the wing design: (1) use of rectangular cross section for wing design with bristles attached laterally; and (2) approximation of wings as rigid structures. In terms of wing planform, tiny insect wings are teardrop-shaped with bristles attached laterally along the span and also longitudinally at the tip. Free flight videos of thrips [13] indicate the wings are flexible, leading to "peel" of the wings as opposed to the fling considered in this study. Forewing images of thrips [31,39] show variation in membrane area and a large number of bristles (≈50-100), neither of which was considered in this study. The structural deformation of a wing, membrane area and number of bristles can be expected to alter the leakiness, strength of LEV and TEV differently than those observed in this paper. These aspects of bristled wing design need to be investigated in future studies.
In addition to aerodynamic performance, bristles may have other benefits that were not considered in this study. Reduction of forces needed to fling a pair of wings at low Re can reduce the muscular energetic workload necessary for a wingbeat, and this factor needs to be investigated further. Wing bristles have been shown to assist in folding and unfolding the wing [22]. Finally, previous studies have suggested that bristles may enhance electrostatic charges to help in dispersal [14,40], or act as airflow and vibration sensors [41,42]. These alternative applications will need to be examined in future studies.

Conclusions
This study shows that bristled wings, such as those seen in the smallest flying insects, can provide a disproportionally large drag reduction with a comparatively smaller loss of lift during clap and fling at chord-based Reynolds numbers on the order of 1-10. Leaky flow in the direction opposite to wing motion was observed in the gaps between bristles, which played a central role in drag reduction by lowering the shear stress along the wing span. Shear layers formed around closely-spaced bristles limited the maximum leaky flow rate by lowering the effective gap through which the fluid can leak through, thereby minimizing the loss of lift.
Supplementary Materials: The following are available online at www.mdpi.com/2311-5521/3/2/44/s1, Video S1: Prescribed wing pair motion based on kinematics shown in Figure 3, viewed from the front of the tank (as in Figure 2A). Three continuous cycles are shown for a representative bristled wing pair (W/S = 0.05). Video S2: Example of raw 2D TR-PIV data acquired on a representative bristled wing pair (W/S = 0.05) during clap (τ= 0.8-1.0 in Figure 3) at Re = 10, as viewed from the bottom of the tank (chordwise view as in Figure 9). Video slowed by 23 times of recorded frame rate for display purpose. Video S3: Example of raw 2D TR-PIV data acquired on a representative bristled wing pair (W/S = 0.05) during fling (τ= 1.0-1.2 in Figure 3) at Re = 10, as viewed from the bottom of the tank (chordwise view as in Figure 10). Video slowed by 23 times of recorded frame rate for display purpose.

Acknowledgments:
The authors would like to thank: Jamey D. Jacob for loaning equipment used for aerodynamic force measurements, Adrian Macias for his assistance in processing PIV data, and three anonymous reviewers for their constructive comments.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A. Validation of Force Measurement Methodology
To validate the method used to determine force coefficients presented in this study, we replicated one of the experiments described by Sunada et al. [32] that examined constant velocity linear translation of a single, rectangular solid wing at a fixed angle of attack (α). We compared the drag coefficients from our measurements, using the methodology described in this paper, to those published in Figure 5B of Sunada et al. [32].

Appendix A.1. Experimental Setup
We modified our clap and fling robotic model (Figure 2) to drive a single rectangular solid wing in constant velocity linear translation, as shown in Figure A1. Translation and rotation of the wing were achieved using two hybrid stepper motors with integrated encoders (same as those shown in Figure 2). In general, the validation setup was similar to the robotic model used for clap and fling motion ( Figure 2). The main modification in Figure 2 was that a single wing was attached to a D-profile shaft using an L-bracket that was identical in design as shown in Figure 2D. 5 chord lengths, and finally linearly decelerated for about 1.5 chord lengths before coming to rest. Another stepper motor was used to vary the angle attack (α) of the wing.

Appendix A.4. Test Conditions
The solid wing model was tested at Reynolds number (Re) of 12, defined using Equation (3). The tank was filled with water-glycerin mixture with kinematic viscosity (ν) of 6.33 × 10 −4 m 2 ·s −1 , measured using the same viscometer as described in Section 2.4 of the paper. The angle of attack α ( Figure A2A) was varied from −20 to 50 degrees, identical to the constant velocity translation experiments in Figure 5B of Sunada et al. [32]. The drag forces on the solid wing were measured during constant velocity stage of linear translation, using uni-axial linear strain gauges mounted on a custom aluminum L-bracket ( Figure A1B). The thickness of the L-bracket used was 0.65 mm. This was thinner than, but similar in geometric design to, the lift L-bracket that was used for clap and fling study ( Figure 2D). The force data that were collected represent the tangential force (F T ) acting on the wing ( Figure A2A). Figure A1. Experimental setup used for validation of aerodynamic force measurement methodology. A single rectangular solid wing (aspect ratio = 3) was prescribed to move in constant translation at a fixed angle of attack (Re = 12), to replicate one of the experiments reported in Sunada et al. [32]: (A) the top view; and (B) the front view of the setup, with coordinate definitions. Drag force acting on the wing was measured for angle of attack ranging from −20 degrees to 50 degrees, using strain gauges bonded to an L-bracket attached to the wing. The L-bracket design and assembly with the wing was identical to Figure 2D. x, y, z are fixed coordinates with respect to the tank; X', Y', Z' are relative coordinates that move with respect to the wing. α = angle of attack; LE = leading edge; TE = trailing edge. Figure A2. Drag force calculation and coefficient of drag (C D ) versus angle of attack (α): (A) Drag force (F D ) was defined as force acting in the direction opposite to wing motion, and calculated using measured forces in tangential direction (F T ). (B) Coefficient of drag (C D ) versus angle of attack (α), where filled squares represent C D obtained from this study, empty squares represent data points extracted from Figure 5B of Sunada et al. [32] using ImageJ software [43]. Axes limits in B are retained the same as in Figure 5B of Sunada et al. [32] to not alter uncertainty associated with data point extraction in ImageJ. LE = leading edge; TE = trailing edge.   (11). Scaled lift and drag coefficients for the solid wing model (W/S = 0) are the same as C L and C D given by Equations (8) and (9), respectively. W/S = 0 (-); W/S = 0.23(---).

Appendix A.6. Results
For a solid wing in constant translation at Re = 12, our study shows that drag coefficient increases with increasing α ( Figure A2B). Our results are in fairly good quantitative agreement as the previous results published in Figure 5B of Sunada et al. [32]. Note that the axes limits in Figure A2B are identical to those in Figure 5B of Sunada et al. [32]. This was because the data points in Figure 5B of Sunada et al. [32] were extracted using ImageJ software [43], and there is visual uncertainty associated with this process. Plotting our data with identical axes limits as in Figure 5B of Sunada et al. [32] would not artificially alter the uncertainty inherent to the image analysis. A slight increase in drag coefficients is seen in our data ( Figure A2B) as compared to those of Sunada et al. [32]. This could potentially be due to the larger thickness of the wings used in this study (3.17 mm) as compared to the 0.6 mm wing thickness in Sunada et al. [32], as an increase in wing thickness can increase the drag force.