Abstract
A prominent feature of gliding flight in snakes of the genus Chrysopelea is the unique crosssectional shape of the body, which acts as the lifting surface in the absence of wings. When gliding, the flying snake Chrysopelea paradisi morphs its circular crosssection into a triangular shape by splaying its ribs and flattening its body in the dorsoventral axis, forming a geometry with fore–aft symmetry and a thick profile. Here, we aimed to understand the aerodynamic properties of the snake's crosssectional shape to determine its contribution to gliding at low Reynolds numbers. We used a straight physical model in a water tunnel to isolate the effects of 2D shape, analogously to studying the profile of an airfoil of a more typical flyer. Force measurements and timeresolved (TR) digital particle image velocimetry (DPIV) were used to determine lift and drag coefficients, wake dynamics and vortexshedding characteristics of the shape across a behaviorally relevant range of Reynolds numbers and angles of attack. The snake's crosssectional shape produced a maximum lift coefficient of 1.9 and maximum lifttodrag ratio of 2.7, maintained increases in lift up to 35 deg, and exhibited two distinctly different vortexshedding modes. Within the measured Reynolds number regime (Re=3000–15,000), this geometry generated significantly larger maximum lift coefficients than many other shapes including bluff bodies, thick airfoils, symmetric airfoils and circular arc airfoils. In addition, the snake's shape exhibited a gentle stall region that maintained relatively high lift production even up to the highest angle of attack tested (60 deg). Overall, the crosssectional geometry of the flying snake demonstrated robust aerodynamic behavior by maintaining significant lift production and nearmaximum lifttodrag ratios over a wide range of parameters. These aerodynamic characteristics help to explain how the snake can glide at steep angles and over a wide range of angles of attack, but more complex models that account for 3D effects and the dynamic movements of aerial undulation are required to fully understand the gliding performance of flying snakes.
INTRODUCTION
‘Flying’ snakes are a group of five species of arboreal snakes in Southeast Asia known for their ability to glide. When gliding, a flying snake undulates from side to side, sending largeamplitude traveling waves posteriorly down the body (Socha, 2002; Socha et al., 2005). Upon becoming airborne, the snake first splays its ribs to the side, flattening out in the dorsoventral axis and creating an unconventional crosssectional shape (Socha, 2011) (Fig. 1). This combination of undulation and body morphing produces a glider with symmetries drastically different from all other biological or engineered flyers: the undulation results in a body that does not exhibit bilateral symmetry, and yet its rib splaying creates an airfoil with fore–aft symmetry. These features make the gliding system of flying snakes strikingly unique compared with any other animal that aerially locomotes, from gliding ants (Yanoviak et al., 2005; Yanoviak et al., 2011) to flapping bats (Tian et al., 2006; Hedenström et al., 2007; Riskin et al., 2008; Hubel et al., 2010).
Although many open questions remain in the study of animal flight, the aerodynamic basis of flight in birds, bats and insects has received focused attention, leading to extensive understanding of kinematics, steady and unsteady effects, 2D and 3D dynamics, wake patterns and fluid–structure interactions (Lehmann, 2004; Wang, 2005; Tobalske, 2007; Hedenström and Spedding, 2008; Song et al., 2008; Usherwood and Lehmann, 2008; Lehmann, 2009; Spedding, 2009; Johansson et al., 2010). In comparison, there have been far fewer studies on the aerodynamics of animals that can only glide (Emerson and Koehl, 1990; McCay, 2001; Bishop, 2007; Alexander et al., 2010; Miklasz et al., 2010; Park and Choi, 2010; Bahlman et al., 2013), despite their large morphological and taxonomic diversity and the far greater number of independent evolutionary origins of gliding flight (Dudley et al., 2007; Dudley and Yanoviak, 2011). In particular, little is known about the physical mechanisms of gliding in snakes, with previous studies focusing on kinematics (Socha, 2002; Socha and LaBarbera, 2005; Socha et al., 2005; Socha, 2006; Miklasz et al., 2010) and only one study to date that considered aerodynamics (Miklasz et al., 2010). Here, we focus on the aerodynamics of crosssectional shape of the paradise tree snake, Chrysopelea paradisi. This species is the most proficient glider of the flying snakes, capable of glide ratios as high as 4.5 and maneuverability by active turning (Socha, 2002; Socha et al., 2005; Socha, 2011).
The crosssectional shape of an airborne C. paradisi varies along its length. Although the total 3D shape is not precisely known, for much of the body between the head and the tail the shape is roughly similar to a rounded triangle (Socha, 2011). At midbody, where the width in the air increases the most relative to its normal condition, the crosssectional shape consists of a semitriangular top (the dorsal surface) and a concave bottom (the ventral surface), with two downwardprojecting ‘lips’ at each edge (Fig. 1). This profile more closely resembles a bluff body than a streamlined airfoil, and the aerodynamic implications of this shape are unknown.
Compared with conventional airfoils, bluff bodies exhibit considerably different aerodynamic performance, characterized by
List of symbols and abbreviations
 AR
 aspect ratio
 C_{D}
 coefficient of drag
 C_{L}
 coefficient of lift
 D
 drag force
 DPIV
 digital particle image velocimetry
 EF
 energy fraction
 FSE
 fluctuating structures energy
 J
 advance ratio
 L
 lift force
 L′
 characteristic length used to calculate St number
 POD
 proper orthogonal decomposition
 PSD
 power spectral density
 Re
 Reynolds number
 RPC
 robust phase correlation
 St
 Strouhal number
 TR
 time resolved
 U
 freestream veolcity
 α
 angle of attack
 λ
 modulation amplitude
 ω
 modulation frequency
The previous study of flying snake aerodynamics used simple semicircular tubing to model the snake's crosssection, and aimed to address the effects of degree of body concavity and presence or absence of a backbone protrusion (Miklasz et al., 2010). These models were coarse and only roughly approximated the snake's true crosssectional shape, but the results suggested that the snake possesses certain favorable aerodynamic characteristics. The most snakelike geometries sustained increases in lift up to an angle of attack of 30 deg, with a maximum lift coefficient of 1.5 and maximum lifttodrag ratios of 2.4–2.8 across a wide range of angles of attack. Stall was gentle, with lift production occurring even at angles of attack of 60 deg. When the ventral concavity was fully filled, the model displayed the smallest lift coefficients, suggesting that the ventrally projecting lips are important features for the snake's force production.
Although the study by Miklasz and colleagues (Miklasz et al., 2010) provided some insight into flying snake aerodynamics, it represents only a preliminary analysis. (1) The simple geometry did not reflect the snake's true crosssectional shape, and small differences in shape may incur large differences in flow patterns. (2) Only one Reynolds number was examined. Chrysopelea paradisi exhibits significant variation in body size and glide speed, resulting in a range of Reynolds numbers that characterize its gliding flight (Re ~5000–15,000) (Socha et al., 2005). Furthermore, in a glide trajectory the forward speed can increase by a factor of six relative to the initial takeoff velocity, meaning that the snake must experience increasing Reynolds number within the course of a single glide. (3) The fluid dynamics were not studied, and therefore the physics that underlie the snake's aerodynamic behavior are not known.
Here, we comprehensively characterized the aerodynamic performance of the 2D crosssectional shape of C. paradisi Boie 1827. We used a highfidelity model of the snake's shape, based on stereo imaging of the snakes in flight (Socha, 2011). We used load cell measurements and timeresolved (TR) digital particle image velocimetry (DPIV), referred to hereafter as TRDPIV, to determine lift and drag coefficients (C_{L} and C_{D}, respectively), wake dynamics and vortexshedding characteristics of the snake's crosssectional geometry across a behaviorally relevant range of angles of attack and Reynolds numbers. The experiments were performed in a lowspeed, lowturbulence water tunnel (flume) at velocities that satisfy Reynolds number dynamic similarity. Our study is the first to quantify the aerodynamic performance of the anatomically similar shape of a flying snake, and specifically evaluates whether the aerodynamic performance of the bluff body alone is sufficient to enable gliding. This is the first study focused on the steadystate aerodynamic characteristics of the flying snakes and will serve as a foundation for future work that will incorporate body movement and unsteady effects. More broadly, this study contributes to our growing understanding of biological lifting surfaces and flight in low Reynolds regimes.
RESULTS
Force measurements: lift and drag coefficients
The snake's crosssectional geometry produced negative lift coefficients for angles of attack of −10 to 5 deg, and positive coefficients thereafter (Fig. 2). The largest increase in lift coefficient occurred between angle of attack α=5 deg and α=35 deg, during which the drag coefficient gradually increased. At α=35 deg, the lift rapidly increased, producing a maximum lift coefficient of 1.9. For higher angles of attack, the lift gradually decreased while the drag rapidly increased. A strong Reynolds number dependence was observed at α=35 deg. At this angle of attack, a peak in lift occurred at higher Reynolds numbers, but did not occur at Reynolds numbers of 3000–7000. The lifttodrag ratio exhibited similar behavior, producing a maximum lift/drag ratio (L/D) of 2.7 at α=35 deg due to the peak in lift for the higher Reynolds numbers.
The results from the TRDPIV force estimates are shown in Fig. 3 and are compared with the force measurements for Re=13,000. These two different methods of determining the aerodynamic forces show good agreement and are within the uncertainty ranges. The TRDPIV lift and drag coefficients were decomposed into the three components: the pressure coefficients and the parallel and perpendicular momentum flux coefficients. Fig. 3 shows that the model produced a large pressure force associated with lift for angles of attack between 20 and 60 deg. At α=35 deg, we observed a reduction in the pressure component accompanied by a negative increase in the parallel momentum, which was offset by preserved normal momentum up to 45 deg followed by a subsequent reduction. For angles of attack greater than 35 deg, the large increase in drag corresponded to a large rise in the pressure coefficient as well as drastic changes in both the parallel and normal momentum components representing the momentum deficit of the wake and signifying the onset of stall. These plots clearly show that the pressure predominantly governs both lift and drag, with counteracting effects between the contributions of the parallel and normal momentum components, and also that increased convective acceleration losses in the streamwise direction result in increased drag and reduced lift.
Fig. 4 compares the force coefficients computed from the different techniques described above for an angle of attack of 35 deg as a function of Reynolds number. There was a steep increase in the pressure force associated with lift between Reynolds numbers of 7000 and 9000.
Proper orthogonal decomposition modes and flow field measurements
The energy fractions and entropy calculations from the proper orthogonal decomposition (POD) analysis are shown in Fig. 5 for Re=13,000. The energy fractions shown in Fig. 5A clearly indicate that the first two modes, which correspond to the Kármán vortex streets produced by the model, are the most energetic modes and contain nearly 70% of the fluctuating structures energy (FSE). POD studies on the wake produced by circular cylinders (Ma, 2000), halfcylinders (Cruz, 2005) and a live swimming fish (Epps, 2010) have also shown that the dynamics of Kármán vortex streets can be accurately reconstructed using the two most energetic modes. The entropy plots in Fig. 5B show that an angle of attack of 35 deg produced the lowest entropy, meaning that this case is more ordered in space and time, and that the dynamics of the flow field can be reconstructed with fewer modes than in the other cases. Although a small number of modes is sufficient to capture the majority of the energy for all angles of attack, Fig. 5A,B reveals that the organization of these vortices is significantly increased for α=35 deg. This observation provides insight that helps to explain the lift enhancement observed in Figs 2 and 3. The entropy gradually increased between α=0 and 20 deg and significantly decreased at 35 deg. The entropy then increased again at α=45 deg and decreased at α=60 deg.
Spatial eigenmodes 1 and 2 are plotted in Fig. 5C for angles of attack of 0, 35 and 60 deg and Re=13,000. These plots show dimensionless vorticity overlaid with the velocity vectors of the mode shapes. Eigenmodes 1 and 2 for α=0 deg exhibit counterrotating vortices that decrease in strength with downstream position. The velocity field shapes also suggest that there may be a vortex structure that formed in the concave ventral surface of the model at this angle of attack, but this region in the flow field was blocked from the camera and was not resolved. The eigenmodes for α=35 deg show similar behavior, but with stronger vortices that were deflected downwards. The vortex patterns of the highest angle of attack are significantly larger than those for the lower angles of attack. Eigenmode 2 contains a pair of smaller vortical structures located directly behind the model, which is then followed by a large vortical structure with a concentrated core and long branches of vorticity that extend from the perimeter of the vortex towards the beginning of the separated boundary layer. The vortices located further downstream show similar structures with decreased strength.
The mean recirculation region of a wake is defined as the region in the flow within which the mean local streamwise velocity is opposite in direction with respect to the freestream velocity, U; namely, u×U<0. Hence, the boundary enclosing the recirculation region produced by the snake model is defined by the locations in the wake along the streamwise direction where u×U=0. These are calculated across the wake by first finding the points where u×U changes sign and then performing a weighted linear interpolation between these two points to find the exact location. This results in the recirculation regions shown in Fig. 6 for Re=13,000.
For a cylinder, the length of the mean recirculation is known to decrease with increasing Re, producing a larger base suction and larger drag coefficients (Williamson, 1996). The snake model showed similar behavior (Fig. 6B), as the length of the recirculation region decreased with Reynolds number and the drag coefficient gradually increased. As the angle of attack increased from 0 to 35 deg, the recirculation region length decreased and was deflected downwards. At the highest angles of attack, the recirculation region enveloped a larger region on the dorsal surface. The length and size of the recirculation region increased between 45 and 60 deg. The decrease of the mean recirculation for 0 to 35 deg implies that the wake vortices are forming closer to the dorsal surface, augmenting the vortexinduced suction and thereby increasing lift. Once the wake elongates for the higher angles of attack, the added suction is lost, reducing lift and increasing drag.
Spectral analysis
The spectra of the force measurements collapsed into two peaks (Fig. 7), which indicates the presence of two different vortexshedding behaviors. Lower angles of attack (up to α=30 deg) produced the higher Strouhal number peak of 0.24, whereas the higher angles of attack (α=35 deg and above) produced the lower Strouhal number peak of 0.16, which also had larger power densities. At low Reynolds numbers, α=35 deg exhibited a Strouhal number slightly less than 0.24, and at higher Reynolds numbers it showed Strouhal numbers larger than 0.24.
Plots of the spatially resolved power spectra of the vcomponent velocity fields are shown in Fig. 9 for Re=13,000. The spectra are calculated along all the TRDPIV velocity measurement points at the X/C=2 location, as shown in Fig. 8B. For every point, the signal spectral energy as a function of frequency is calculated and we construct contours of spectral energy as a function of dimensionless frequency (St number) plotted along the xaxis and spanwise (Y/C) location on the yaxis. Y/C=0 corresponds to the centerline of the wake.
The dominant spectral components observed in these spectral maps (Fig. 9) show good agreement with the power spectral densities (PSDs) of the force measurements (Fig. 7). The angle of attack of 0 deg showed a narrow band in frequency and vertical position, with most of the power concentrated at the vortexshedding frequency of St≈0.24. This angle also showed the greatest symmetry in the wake about Y/C=0. Moreover, by observing the highintensity spectral band, the width of the wake can be directly estimated to be between Y/C=−0.6 and 0.6. As the angle of attack increases to 20 and 35 deg, the energy spectrum spreads, and the wake shifts downwards, becoming asymmetric around the centerline and expanding in width. The angle of attack of 35 deg showed the greatest degree of asymmetry about Y/C=0 as the wake was deflected further downward. Angles of attack of 45 and 60 deg exhibited the highest power magnitudes but here the dominant shedding frequency is reduced to approximately St≈0.16. These observations are consistent with the shedding frequencies estimated from the force spectra (Fig. 7). However, the TRDPIV spectra further reveal that for the higher angles of attack, additional dominant frequencies emerge, shown as high power peaks at higher frequencies. This yields a broad spectrum and suggests a more complex and turbulent wake.
Trajectories and strength of coherent structures
The mean vortex trajectories for the upper and lower vortices for Re=13,000 are plotted in Fig. 10A. The analysis of the vortex trajectories provides insight into the vortexformation and shedding processes and the downstream evolution of the wake. The initial stage of the trajectories for state I, which corresponds to the vortexformation process, is much more compact with a difference in initial vertical position of ~0.6 chords. The initial points of the state II trajectories covered a much wider range with a vertical difference of over 1 chord. The trajectories for the vortex pairs of state I converged on nearly the same trajectory at ~1.5 chord lengths downstream of the model's centroid. The trajectories at higher angles of attack did not converge. The point of maximum circulation for angles of 0 to 35 deg occurred at ~1 chord downstream of the model. The maximum circulation of the upper vortex for α=45 deg occurred closer to the model than for the other angles of attack, and at α=60 deg it was located the farthest downstream.
DISCUSSION
This study investigated the aerodynamic characteristics of the crosssectional shape of one species of flying snake (C. paradisi) across a behaviorally relevant range of Reynolds number and angle of attack. Because this shape occurs along a large percentage of the length of the body, it must contribute greatly to the snake's force production while gliding. Overall, this shape exhibited a robust aerodynamic performance, maintaining high lift coefficients between α=20 and 60 deg and nearmaximum L/D values over a range of angles of attack between 15 and 40 deg. As shown in Fig. 3, the pressure term from the TRDPIV force decomposition had the greatest influence on the total aerodynamic forces. Low angles of attack (α=−10 to 5 deg) produced negative lift coefficients because the concave ventral surface produced a lowpressure zone, perhaps due to the formation of a vortex ‘trapped’ in the cavity. These angles of attack also exhibited the smallest drag coefficients, likely because this is the range where the crosssectional shape presents the minimum perpendicular projection.
The angle of attack of 35 deg produced the largest lift coefficient of 1.9, which occurred due to a spike in the pressure term for the higher Reynolds numbers (greater than 7000). Initially, we were concerned that this peak was an experimental artifact, but its occurrence in two different load cell setups, in the independently derived TRDPIV force measurements and in the previous study of flying snake shape (i.e. the coarse semicircular model) (Miklasz et al., 2010) all confirmed that the peak is real. The vortex dynamics and vortex trajectory analysis help to explain the mechanism responsible for this behavior (Figs 5, 6, 10). The POD entropy shown in Fig. 5B reveals that for α=35 deg the vortices are more organized. The locations at which the vortex trajectories merge reveal the approximate wake closure point, which in turn is related to the base pressure and form drag of the body. The combination of more coherent vortices with the closure of the wake in closer proximity to the suction side of the airfoil helps explain the vortexinduced lift augmentation observed for α=35 deg. Moreover, the degree of asymmetry of the wake and the strength of the formed vortices in the proximity of the suction side of body (the top surface) suggest that a vortexinduced liftgeneration process is at play. For the case of α=35 deg, we notice that peak circulation and trajectory convergence points are closer to the body suction side and that the wake exhibits strong downwash (higher deflection of the trajectories downwards). As a result, the higher Reynolds numbers produced a larger highpressure zone in the ventral region and a more significant negative pressure peak at the dorsal surface. This also suggests that the spike in lift and L/D observed at α=35 deg is the result of vortexinduced lift processes.
Vortexshedding behavior
The congruence of results from different methods of analysis clearly indicates that the snake's crosssectional shape displayed different vortexshedding behaviors, which depend on the angle of attack. Two different Strouhal numbers were observed in the spectral analysis, indicating that the lower angles of attack (α=0–30 deg) and higher angles of attack (α=40–60 deg) demonstrated different shedding characteristics as the separation point on the dorsal surface shifted from the rounded apex to the leading edge. Lower angles of attack, corresponding to a Strouhal number of 0.24, produced smaller vortices (as shown by the POD analysis). Lower angles of attack showed greater spectral organization, with the majority of the power concentrated at the dominant Strouhal number (supplementary material Movies 1, 2). Higher angles of attack (St=0.16) produced much larger vortical structures and circulation magnitudes, corresponding to increased drag and fluctuations in the fluid dynamic forces (supplementary material Movie 3).
The vortexshedding behavior of the snake model is similar to previously studied National Advisory Committee for Aeronautics (NACA) airfoils. In particular, thick, symmetric airfoils of 15–35% thickness also produced two different vortexshedding mechanisms that depend on the angle of attack and thickness (Sarraf, 2010). A sharp decrease in Strouhal number (from 0.3 to 0.2) was similarly observed when these airfoils reached a critical angle of attack, whereby the separation point shifted on the suction side towards the leading edge. Previous studies on rounding corners of square cylinders showed that the Strouhal number increased with increasing rounding radius (Hu and Zhou, 2009). Similarly, the snake's crosssectional shape has a large radius of curvature near the apex when oriented at lower angles of attack and produced a higher Strouhal number. At high angles of attack, the shape presents a sharp radius of curvature at the leading edge and produced correspondingly lower Strouhal numbers.
Comparison with bluff bodies and lowspeed airfoils
We compared the aerodynamic performance of the snake's crosssectional shape with that of some wellknown bluff bodies and airfoils (Figs 11, 12). Surprisingly, the snake's geometry generally outperformed the other airfoils and bluff bodies within the snake's low Reynolds number range. Compared with the previous semicircular models tested as a firstorder approximation (Miklasz et al., 2010), the anatomically similar snake shape exhibited a larger maximum lift coefficient and higher drag coefficients. The semicircular models showed a sharp increase in L/D at α=30 deg, producing a maximum L/D of 2.8, and also maintained nearmaximum L/D values over a wide range of angles of attack. The semicircular model produced a more rapid increase in L/D between α=0 and 10 deg, but our anatomically similar model exhibited a higher stall angle and larger lift coefficients after stall.
The dorsoventral (top–bottom) asymmetry of the snake's crosssectional shape appears to play an important role in aerodynamic force production. Specifically, the concave ventral surface with ventrally protruding lips is an important liftproducing feature that may allow the snake's crosssectional shape to generate large lift coefficients over a wide range of angles of attack by maintaining a high pressure zone within the cavity. Other airfoils and bluff bodies with camber, such as the circular arc (Sunada, 1997), halffull and empty semicircular flying snake models (Miklasz et al., 2010), and the S1223 airfoil (Selig, 2005) generated significantly larger maximum lift coefficients than symmetric airfoils, such as the NACA 0012 (Alam et al., 2010) and NACA 00150035 (Sarraf, 2010), which produced maximum C_{L} values of 1 to 1.2, or an insectinspired corrugated airfoil, which produced maximum C_{L} values of 0.8 to 1.0 (Murphy and Hu, 2009). Miklasz and colleagues (Miklasz et al., 2010) found that filling in the concave ventral surface of the semicircular model decreased the maximum lift coefficient by more than 20%. Thickness also appears to be an important characteristic, as the thicker airfoils and bluff bodies produced higher stall angles. Studies on thick, symmetric airfoils (Sarraf, 2010) found that increasing airfoil thickness from 15 to 35% decreased the slope of the lift curve at low angles of attack, but also increased the stall angle from 21 to 40 deg. Our model of the snake's crosssectional shape possessed a thickness of 39% (relative to the chord) and stall angle of 35 deg. This greater thickness may have led to an increased drag coefficient because the perpendicular projection of the airfoil is larger. If true, this may also explain why the anatomical snake model produced larger drag coefficients than the semicircular model, which had a thickness of 29% (Miklasz et al., 2010).
The lift curve of the snake's crosssectional geometry displays significantly different behavior compared with more conventional streamlined airfoils (e.g. the S1223 airfoil) (Selig and Guglielmo, 1997; Selig, 2005) (Fig. 12). Streamlined airfoils with positive camber, which are typically operated at higher Reynolds numbers, produce positive lift coefficients at α=0 deg and rapid increases in lift up to the stall angle, which often occurs at α=20 deg or less. After this stall angle, the lift coefficient decreases rapidly, possibly resulting in a catastrophic stall for an aerial vehicle. The snake model did not produce positive lift coefficients until α=10 deg, but maintained increases in lift up to α=35 deg and exhibited a gradual stall region. These characteristics may explain the snake's favorable forceproducing mechanisms despite its continuously changing posture during a trajectory, in which the body's orientation in the pitch axis may vary by over 60 deg (Socha et al., 2005). In particular, the gradual stall behavior at high angles of attack may allow the snake to develop lift during the ballistic dive phase of its trajectory, during which the body likely encounters very high angles of attack (i.e. >45 deg). This may be critical for the process of transition between falling and gliding in any glide event.
Implications for gliding in flying snakes
Although the glide angle and relative body orientation of C. paradisi has been identified in previous kinematic studies (Socha et al., 2005; Socha et al., 2010), the actual local Reynolds number and angle of attack along the snake's body are unknown, and these may vary significantly given the snake's complex kinematics (Socha et al., 2005; Miklasz et al., 2010). The maximum lift coefficient and L/D of the snake's crosssectional shape occurred at an angle of attack of 35 deg. Our results suggest that, for maximum horizontal gliding, the snake should optimize its glide angle and body posture to produce a 35 deg angle of attack over the majority of its body. While the actual variation in angle of attack along the snake's body is unknown, the snake has been shown to orient its whole body roughly at an average angle of 25 deg (in a flowrelevant reference frame) while gliding (Socha et al., 2010). If this angle is reflective of the actual angle of attack along the body, the snake may be operating in a nearmaximum lift production range. Considering the poor aerodynamic performance of the snake's crosssectional shape in the negative lift range, it seems unlikely the snake uses low angles of attack (<10 deg) during gliding. However, it could possibly use this negative lift range for purposes of control, such as to change its pitch angle or to induce a turn. The fore–aft symmetry of the snake's crosssectional shape appears to be an important characteristic, because it means that when the leading and trailing edges reverse positions during the snake's undulatory motion, similar aerodynamic performance is maintained despite the radically different orientations.
The snake's crosssectional shape maintained high lift values and nearmaximum lifttodrag ratios over a wide range of angles of attack, which may help explain the snake's performance as the glide angle changes within a glide trajectory (Socha et al., 2005). After the snake jumps into the air and begins the ballistic dive, it falls at a steep angle (as great as 62 deg) as a result of low lift production. As the snake accelerates, it produces more lift and its glide angle shallows to as little as 13 deg, and it can even produce a net positive upwards acceleration (Socha et al., 2010). Our results suggest that snakes using α=35 deg would experience a boost in lift as the Reynolds number increased across the threshold of Re=7000, enabling the snake to produce more shallow glides and cover greater horizontal distances.
Previous studies on flying snake kinematics used assumptions of equilibrium gliding to estimate the maximum L/D values during the shallowing phase of the trajectory; these values are 14–60% larger than those reported here (3.08–4.33 versus 2.7) (Socha et al., 2005; Socha et al., 2010). Based on the results from this study (L/D_{max}=2.7) and assuming equilibrium gliding, C. paradisi should not be able to produce glide angles of less than 20 deg, but in fact glide angles as shallow as 13 deg have been observed (Socha et al., 2005). This demonstrates that our 2D models are insufficient to fully capture the performance of real snakes, or that the snakes are not using equilibrium gliding. Additionally, it has been observed that intermediatesized C. paradisi are better gliders than the largest snakes, capable of covering greater horizontal distances at lower minimum glide angles (Socha and LaBarbera, 2005). These smaller snakes operate at lower Reynolds numbers (resulting from their smaller body width and lower glide speeds). Their superior performance relative to larger snakes stands in contrast to our study of 2D snake geometry, which indicates that the best aerodynamic performance should occur at higher Reynolds numbers. Considering that smaller snakes appear to be morphologically and physiologically similar to the larger ones (Socha and LaBarbera, 2005; Socha, 2011), this contrasting evidence suggests that smaller snakes may compensate for the aerodynamic deficits of the lower Reynolds number by exploiting other mechanisms for lift augmentation and drag reduction. Alternatively, smaller snakes may exhibit subtle differences in shape that have not yet been recorded. Lastly, the Re=7000 transition may help to explain a further feature of body size on glide performance. Although intermediatesized snakes glide more effectively than the largest snakes, the smallest snakes also exhibit lower performance. It is possible that the intermediatesized snakes cross the Reynolds number threshold, and thus experience enhanced lift when at angles of attack of 35 deg.
In conclusion, our analysis of a 2D geometry only partially explains the snake's gliding performance. Other mechanisms, including threedimensional and unsteady effects, may augment the snake's aerodynamic force production beyond what is delivered by its crosssectional shape. In addition, as the snake undulates, its body posture is constantly reconfigured such that the spacing between upstream and downstream body segments is continuously changing. Thus, the wake from upstream may interact with downstream body segments, which may have a significant effect on the snake's wholebody aerodynamic performance. Sweep angle effects, which were also neglected in this study, may produce complex 3D flow structures that could also change the aerodynamic performance. Overall, the aerodynamic forces produced by the snake's unique mode of aerial locomotion likely result from a complex interaction between its morphological characteristics and its dynamic gliding behavior, utilizing multiple physical mechanisms. Future studies on more complex models, including those that incorporate dynamic movements, are required to fully understand the aerodynamic performance of flying snakes.
MATERIALS AND METHODS
Two types of experimental measurement techniques were used for this study. Load cell measurements were performed to determine the fluid dynamic forces, and TRDPIV was used to obtain flow field measurements. From these measurements, several different types of analyses were employed to examine the aerodynamic performance of C. paradisi's crosssectional shape, including a spectral analysis, a PIVbased estimation of the aerodynamic forces, POD and coherent structure identification and tracking.
Justification for 2D and quasisteady modeling
Unsteady 3D effects are well known to affect the aerodynamics of flapping flyers (Wang et al., 2004; Wang, 2005; Warrick et al., 2005; Johansson et al., 2010), and conventional, steadystate aerodynamics are insufficient for fully explaining how insects, birds and bats fly. A useful index for determining whether a quasisteady approach can be used for studying animal flight is the advance ratio, J, which compares the animal's forward motion with the reciprocating motion of its wings. For flapping flyers, J is of the order of 1, with values for insects extending lower and large birds reaching as high as J=4 (Vogel, 2003; Dickson and Dickinson, 2004). Using typical values of gliding (Socha et al., 2005), J can be estimated in flying snakes as: (1) where U is glide speed, ω is undulation frequency and λ is undulation amplitude. Here, we assume that the cyclic sidetoside and vertical movement of the body involved in aerial undulation is the reciprocating motion. The high value of J for the snake is more than an order of magnitude greater than that of flapping flyers, suggesting that the snake's undulation is less important for force generation than is flapping. A second approach is to compare the glide speed with the speed of the posteriorly directed traveling wave that moves down the body: (2) where U_{w} is the traveling wave speed. This value of J is even higher. Together, these indices strongly suggest that the snake's aerodynamics are dominated by the overall forward motion of the snake and not its aerial undulation. Therefore, the forces on the snake can be reasonably approximated by the steadystate conditions of any particular body posture. This dimensional analysis suggests that a quasistatic approach is reasonable for a firstorder study of the snake's aerodynamics.
In this study, we modeled the snake's body using a large aspect ratio (AR) model in order to isolate the effects of crosssectional shape and angle of attack (α) on aerodynamic performance. Generally for wings with low AR, the lift and drag coefficients can be strongly influenced by the AR due to endeffects, threedimensionality and liftinduced drag (Abbott and Von Doenhoff, 1959). However, for higher AR wings, the lift and drag coefficients are affected less by AR and instead depend more directly on α. During its glide, the snake assumes an Slike shape composed of two to three straight sections connected with curved segments. Each straight section can have an AR of the order of 8–10 (Socha and LaBarbera, 2005; Socha, 2011), which can be considered in the medium to high AR range. This AR, in combination with the high advance ratio, suggests that a 2D and quasisteady approach is a reasonable first approximation for an initial investigation of the effects of the aerial snake's crosssectional shape and angle of attack. This argument only applies for portions of the snake's body that are orthogonal to the freestream flow. Furthermore, this does not suggest that 3D effects are not important, but simply that the general approach used in this study may yield insight into one component of the snake's gliding apparatus.
Experimental model and test facility
The experimental model is an anatomically similar replica of the snake's crosssectional geometry, designed from an analysis of aerial images of the snake's flattened body (Socha, 2011). The geometry was scaled to give the experimental model a chord (c) of 25.4 mm, which represents the maximum body width of the snake while gliding. As in other studies, the chord was used as the characteristic length to determine the Reynolds number (Socha et al., 2005; Miklasz et al., 2010; Socha et al., 2010). The model was rapid prototyped from ABS plastic with an aspect ratio of 20 to simulate a single 2D segment of the snake's body when gliding and to reduce end effects. Aluminium rods were embedded in the model to increase its strength and stiffness. The experiments were conducted in a closedloop water tunnel at Virginia Tech. The tunnel has a test section of 0.6×0.6×1.5 m and can produce a lowturbulence freestream of less than 3% at speeds up to 1 m s^{−1} and flow uniformity at the test location better than 99% (Weiland and Vlachos, 2009).
Aerodynamic force measurements
The lift and drag produced by the snake's crosssectional geometry were measured as a function of Reynolds number and angle of attack to determine the snake's overall aerodynamic performance. To encompass a large range of possible snake flight kinematics, we explored a parameter space of seven different Reynolds numbers ranging from 3000 to 15,000, and 15 angles of attack (α) ranging from −10 to 60 deg (Fig. 13). The experimental setup was contained in a rigid truss system mounted in a water tunnel (Fig. 13). Lift and drag forces were measured simultaneously using four 5 lb (22 N) load cells (LCFD5, accuracy: ±0.15% FSO, DMD465WB, straingage amplifiers; Omega Engineering, Stamford, CT, USA) mounted to the model support arms. These support arms were completely encased by acrylic sidewalls to prevent interaction with the flow in the water tunnel. The leading edge of the acrylic sidewalls was designed with superelliptical geometry (Narasimha, 1994), which promoted the formation of a stable boundary layer and prevented separation. To eliminate 3D effects, the model spanned the entire width between the sidewalls with a clearance of ~1 mm.
The data were recorded using LabVIEW software (National Instruments, Redmond, WA, USA) and a DAQ board (NI PCI6251) at a sampling rate of 1000 Hz. Each trial consisted of 30 s of recording. The weight of the experimental setup and model was accounted for by zeroing the force measurement before each trial. The force coefficients were defined as: (3) (4) where L and D are the measured lift and drag forces, ρ is the density of water, U is the freestream velocity, l is the span of the model and c is the chord length of the model.
To investigate how vortex shedding varied with Reynolds number and angle of attack, we performed spectral analysis on the force time series to determine the dominant spectral components of the force measurements. The frequencies from the spectral analysis were nondimensionalized by calculating the corresponding Strouhal number: (5) where f is the vortexshedding frequency and L′ is a characteristic length, defined as the projection of the model perpendicular to the freestream (Miranda, 2005) (Fig. 8A). This dimension (L′) essentially scales the characteristic length scale of the wake. Spectral analysis was also preformed on the timeresolved velocity field measurements (described in ‘Flow field measurements’, below) and the results were nondimensionalized the same way and were compared against the force spectral analysis. The velocity field spectra were calculated along the red dashed line shown in Fig. 8B. Several different streamwise locations were analyzed, and a position of 2 chords downstream of the model's centroid was selected because it showed the highest energy and significant variation between angles of attack. It should be noted that the mean closure point of the wake is also expected to be at approximately X/C=2.
Flow field measurements
We used TRDPIV (Willert and Gharib, 1991; Adrian, 2005) to measure instantaneous, spatially resolved velocity vectors to quantify flow field characteristics including vorticity, strain rate and field pressure. This analysis was conducted on a subset of the total experimental parameter space, with case selection informed by the force measurement results. Five angles of attack were tested at a Reynolds number of 13,000 and four Reynolds numbers were tested for an angle of attack of 35 deg. This selection of the TRDPIV test conditions was based on observations from the load cell measurements. For all angles of attack except 35 deg, we observed negligible Reynolds dependence. The angle of attack variation at Re=13,000 was selected because it was more representative of gliding by an adult snake.
In these tests, the experimental model was rigidly mounted to acrylic sidewalls, as shown in Fig. 13B. A concave mirror was mounted above the snake model at the top of the experimental setup to reflect the laser plane back through the test section, eliminating the shadow cast by the model. The mirror was mounted on an acrylic ‘boat’, which slightly depressed the free surface at the top of the water tunnel to create a stable interface between the surface and acrylic for the laser plane to pass.
The water tunnel was seeded with small, neutrally buoyant, hollow glass spheres (mean diameter, 126.4 μm) that acted as flow tracers and were illuminated by a planar laser sheet (~1 mm thickness) using a 527 nm wavelength, 10 mJ dual head laser system (Pegasus PIV, New Wave Research, Portland, OR, USA). A highspeed video camera (XS3, Integrated Design Tools, Tallahassee, FL, USA) was used to record doublepulsed images of particle motion at a pairsampling rate of 500 Hz, which was sufficient to resolve a wide range of frequencies in the wake, with a maximum expected vortexshedding frequency of 15 Hz. The interrogation region for angles of attack up to 35 deg was 1280×512 pixels and the interrogation region of the higher angles of attack was 1280×1024 pixels, which was required to capture the full width of the wake. The magnification was 109.34 μm pixel^{−1}, producing a particle image size of 2–4 pixels.
Robust phase correlation (RPC) (Eckstein, 2008; Eckstein, 2009a; Eckstein, 2009b) was used with multigrid iterative deformable windows (Scarano, 2002) to correlate the image pairs and to measure the instantaneous velocity fields. Gaussian image masking and spatial filtering of the velocity fields between iterations were employed to stabilize the iterative process (Schrijer, 2008). All TRDPIV measurements were processed using two passes with 64×64 pixel interrogation windows and three passes with 32×32 pixel interrogation windows. After each pass, the vectors were validated using velocity thresholding and universal outlier detection (Westerweel, 2005). The final velocity fields had a uniform vector grid spacing of 8 pixels.
Force calculations from TRDPIV measurements
A velocity field momentumbased approach (e.g. Unal, 1997; Oudheusden, 2006; Ragni, 2009) was used to estimate the lift and drag coefficients from TRDPIV measurements for comparison with the direct load cell force measurements. The application of the integral form of the momentum equation to a control volume that is statistically steady yields the following integral: (6) where is the aerodynamic force acting on the model, ρ is the density of water, is the velocity vector, is the outward pointing unit normal vector, p is the pressure and τ is the viscous stress tensor (Ragni, 2009). The boundary of the control volume, shown in Fig. 8B, was taken sufficiently far from the experimental model so that viscous stresses and turbulent stresses are negligible. Eqn 6 can be simplified and rewritten in Cartesian coordinates for the control volume to yield the aerodynamic forces (per unit length): (7) (8) where u and v are the x and y components of velocity (Ragni, 2009). The forces resulting from the shear stresses along the vertical and horizontal surfaces of the control volume were computed and found to be negligible, and were therefore neglected in Eqn. 6.
The first terms, which will be referred to as the parallel momentum fluxes, and the second terms, referred to as the perpendicular momentum fluxes, can be calculated directly from the flow field measurements. The last term is the pressure force acting on the control surface boundaries, which was computed using the omnidirectional virtual boundary integration technique (Liu, 2006; Charonko et al., 2010). The use of this approach enables a direct calculation of the flow components that contribute to the generation of forces exerted on the body. In particular, we can determine the relative significance of changes in the parallel and normal convective acceleration around the body with respect to each other and with respect to the pressure components. Such physical insight could not be attained by force load cell measurements alone.
POD
POD (also known as Karhunen–Loéve decomposition) is a mathematical technique used to decompose an ensemble of signals into orthogonal modes that optimally capture the average energy content (Smith, 2005). POD has been used in a wide variety of problems in fluid mechanics to decompose complicated flow fields into their energetically dominant eigenmodes (Berkooz, 1993; Ma, 2000; Cruz, 2005; Smith, 2005; Epps, 2010). Here, we used POD to compare the dominant, orthogonal modes of the different vortexshedding behaviors. We chose POD because it provides the ability to quantitatively compare the spatiotemporal characteristics of the flow across angle of attack and Reynolds number. In contrast, direct observation of the flow field can provide only qualitative comparisons, and other approaches such as spectral analysis are limited to the temporal domain only. For the purposes of this work, it is important to analyze the combined spatiotemporal dynamics of the flow, and therefore a spatial or temporalonly type of analysis would not suffice.
The POD method was implemented for the mean subtracted flow fields using the method of snapshots (see Sirovich and Kirby, 1987; Smith, 2005). The fluctuating structures energy (FSE) and energy fraction (EF) of eigenmode i are defined as: (9) (10) where λ_{i} is the eigenvalue of mode i, and N is the total number of modes. FSE corresponds to the cumulative energy content of N modes and EF provides the ratio of the energy of an individual mode with respect to this total. Hence, EF will take values between 0 and 100%, with the most dominant modes having high fractions and the lowest ones asymptotically approaching zero. A very organized wake would have a small number of very dominant modes while the vast majority of the modes will have very low energy content. In contrast, a very turbulent flow will have a more uniform mode energy distribution.
Entropy is a measure of the complexity and order of a spatiotemporal structure, and was computed as follows: (11) The entropy can range from 0 to 1, indicating that the energy is contained in the first mode or the energy is evenly distributed between all modes, respectively (Cruz, 2005). The use of entropy allows us to consolidate the information contained across the modes to one single number for each tested condition. Similar to the interpretation of EF, an entropy value of near zero would indicate a very organized flow, whereas an entropy value near one would imply that the energy in the flow is more uniformly distributed across all modes, indicating a more turbulent flow.
However, it should be noted that a mode is not necessarily a vortical structure, because a vortex can be composed of the superposition of multiple modes whose strengths vary as a function of time.
Identification and tracking of coherent structures
Bluff body flows are dominated by the formation and shedding of vortices. Analysis of these coherent structures can give insight into complex wake dynamics. Two different methods were used to identify vortex cores: a pointbased velocity ‘labeling’ algorithm and the swirling strength (λ_{ci}) criterion (Zhou, 1999; Jiang et al., 2002). Two passes were performed on the velocity fields to identify possible vortex cores. The first pass labeled each grid point in the flow field depending on which quadrant the velocity vector pointed to. The second pass checked each point's surrounding neighbors for the directionspanning property to see whether these neighbors contained each of the four possible direction labels and was therefore at the center of swirling motion. Next, the swirling strength of these possible vortex core locations, computed using the λ_{ci} criterion, was compared against a threshold value [2.6% of the maximum (Chakraborty et al., 2005)] to eliminate incorrect identifications. The swirling strength of each grid point was calculated by computing the imaginary part of the complex eigenvalue of the velocity gradient tensor. The circulation (Γ) was calculated around isocontours of the swirling strength. The vortextracking algorithm described previously (Gifford, 2011) was used to compile the trajectory of each vortex.
FOOTNOTES

Author contributions
D.H., J.J.S. and P.P.V. designed the experiments; D.H. conducted the experiments; D.H., J.J.S., N.D.C. and P.P.V. analysed the data and wrote the manuscript.

Competing interests
The authors declare no competing financial interests.

Funding
This research was partially supported by Defense Advanced Research Projects Agency (DARPA) grant W911NF1010040 to J.J.S. and P.P.V.

Supplementary material
Supplementary material available online at http://jeb.biologists.org/lookup/suppl/doi:10.1242/jeb.090902//DC1
 © 2014. Published by The Company of Biologists Ltd