Aerodynamics of the flying snake Chrysopelea paradisi: how a bluff body cross-sectional shape contributes to gliding performance

‘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 large-amplitude 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 cross-sectional 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 cross-sectional 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 cross-sectional 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 mid-body, where the width in the air increases the most relative to its normal condition, the cross-sectional shape consists of a semi-triangular top (the dorsal surface) and a concave bottom (the ventral surface), with two downward-projecting ‘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

vortex shedding and significant profile drag due to separation. Generally, bluff body flows exhibit interaction of the separating boundary layer with the free shear layer and the wake (Williamson, 1996). The vortex-shedding process begins as the separated boundary layer supplies vorticity to the shear layer, which in turn rolls into a vortex. The vortex grows in strength and size until it reaches a point when the supply of vorticity from the shear layer is interrupted, and, through interaction with the opposite side vortex, the vortex sheds downstream and a vortex streak ensues (Williamson, 1996). Studies on cylinders with rectangular cross-section (Yoon et al., 2010), asymmetric bluff bodies (Hu and Zhou, 2009) and thick airfoils (Sarraf, 2010) have shown that shape asymmetry, thickness and angle of attack affect the flow structure, forces and vortex-shedding characteristics.

The previous study of flying snake aerodynamics used simple semi-circular tubing to model the snake's cross-section, 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 cross-sectional shape, but the results suggested that the snake possesses certain favorable aerodynamic characteristics. The most snake-like geometries sustained increases in lift up to an angle of attack of 30 deg, with a maximum lift coefficient of 1.5 and maximum lift-to-drag 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 cross-sectional 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 take-off 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 cross-sectional shape of C. paradisi Boie 1827. We used a high-fidelity model of the snake's shape, based on stereo imaging of the snakes in flight (Socha, 2011). We used load cell measurements and time-resolved (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 vortex-shedding characteristics of the snake's cross-sectional geometry across a behaviorally relevant range of angles of attack and Reynolds numbers. The experiments were performed in a low-speed, low-turbulence 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 steady-state 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.

The snake's cross-sectional 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 lift-to-drag 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.

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), half-cylinders (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 counter-rotating 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 vortex-induced 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.

The spectra of the force measurements collapsed into two peaks (Fig. 7), which indicates the presence of two different vortex-shedding 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 v-component 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 x-axis and spanwise (Y/C) location on the y-axis. 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 vortex-shedding frequency of St≈0.24. This angle also showed the greatest symmetry in the wake about Y/C=0. Moreover, by observing the high-intensity 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.

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 vortex-formation and -shedding processes and the downstream evolution of the wake. The initial stage of the trajectories for state I, which corresponds to the vortex-formation 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.

This study investigated the aerodynamic characteristics of the cross-sectional 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 near-maximum 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 low-pressure 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 cross-sectional 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 semi-circular 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 vortex-induced 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 vortex-induced lift-generation 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 high-pressure 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 vortex-induced lift processes.

The congruence of results from different methods of analysis clearly indicates that the snake's cross-sectional shape displayed different vortex-shedding 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 vortex-shedding 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 vortex-shedding 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 cross-sectional 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.

We compared the aerodynamic performance of the snake's cross-sectional shape with that of some well-known bluff bodies and airfoils (Figs 11, 12). Surprisingly, the snake's geometry generally out-performed the other airfoils and bluff bodies within the snake's low Reynolds number range. Compared with the previous semi-circular models tested as a first-order approximation (Miklasz et al., 2010), the anatomically similar snake shape exhibited a larger maximum lift coefficient and higher drag coefficients. The semi-circular models showed a sharp increase in L/D at α=30 deg, producing a maximum L/D of 2.8, and also maintained near-maximum L/D values over a wide range of angles of attack. The semi-circular 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 cross-sectional shape appears to play an important role in aerodynamic force production. Specifically, the concave ventral surface with ventrally protruding lips is an important lift-producing feature that may allow the snake's cross-sectional 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), half-full and empty semi-circular 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 0015-0035 (Sarraf, 2010), which produced maximum C_L values of 1 to 1.2, or an insect-inspired 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 semi-circular 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 cross-sectional 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 semi-circular model, which had a thickness of 29% (Miklasz et al., 2010).

The lift curve of the snake's cross-sectional 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 force-producing 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.

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 cross-sectional 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 flow-relevant 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 near-maximum lift production range. Considering the poor aerodynamic performance of the snake's cross-sectional 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 cross-sectional 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 cross-sectional shape maintained high lift values and near-maximum lift-to-drag 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 intermediate-sized 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 intermediate-sized snakes glide more effectively than the largest snakes, the smallest snakes also exhibit lower performance. It is possible that the intermediate-sized 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 three-dimensional and unsteady effects, may augment the snake's aerodynamic force production beyond what is delivered by its cross-sectional 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 whole-body 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.

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 cross-sectional shape, including a spectral analysis, a PIV-based estimation of the aerodynamic forces, POD and coherent structure identification and tracking.

In this study, we modeled the snake's body using a large aspect ratio (AR) model in order to isolate the effects of cross-sectional 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 end-effects, three-dimensionality and lift-induced 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 S-like 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 quasi-steady approach is a reasonable first approximation for an initial investigation of the effects of the aerial snake's cross-sectional 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.

The experimental model is an anatomically similar replica of the snake's cross-sectional 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 closed-loop water tunnel at Virginia Tech. The tunnel has a test section of 0.6×0.6×1.5 m and can produce a low-turbulence freestream of less than 3% at speeds up to 1 m s⁻¹ and flow uniformity at the test location better than 99% (Weiland and Vlachos, 2009).

The lift and drag produced by the snake's cross-sectional 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 (LCFD-5, accuracy: ±0.15% FSO, DMD-465WB, strain-gage 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 super-elliptical 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.

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 high-speed video camera (XS-3, Integrated Design Tools, Tallahassee, FL, USA) was used to record double-pulsed images of particle motion at a pair-sampling rate of 500 Hz, which was sufficient to resolve a wide range of frequencies in the wake, with a maximum expected vortex-shedding 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⁻¹, 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.

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 omni-directional 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 (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 vortex-shedding behaviors. We chose POD because it provides the ability to quantitatively compare the spatio-temporal 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 spatio-temporal dynamics of the flow, and therefore a spatial- or temporal-only type of analysis would not suffice.

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.

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 point-based 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 direction-spanning 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 iso-contours of the swirling strength. The vortex-tracking algorithm described previously (Gifford, 2011) was used to compile the trajectory of each vortex.