Quasi-steady aerodynamic theory under-predicts glide performance in flying snakes Free

Animal flight requires producing and controlling aerodynamic lift and drag forces via neuromuscular control and/or specialized morphology. Flapping fliers use paired beating wings, while most vertebrate gliders deploy pairs of stretched skin between limbs (Byrnes et al., 2008; Bahlman et al., 2013; Khandelwal and Hedrick, 2022), extendable ribs (McGuire, 2003; McGuire and Dudley, 2005), webbed feet (McCay, 2001) or modified fins (Park and Choi, 2010; O'Dor et al., 2013). By contrast, flying snakes flatten their bodies and undulate through the air, turning the whole animal into an ‘S'-shaped morphing wing-body (Socha, 2002). The wing-body changes configuration continuously throughout the glide as the body forms straight segments connected by tight lateral bends that travel posteriorly (Fig. 1D) (Socha et al., 2005). The snake's aerial undulation is composed of horizontal and vertical waves which are 90 deg out of phase and whose frequencies differ by a factor of 2 (Yeaton et al., 2020). While airborne, the posture of the body is never symmetric about any axis, which has functional implications for force production and balance, and distinguishes flying snakes from most other flapping flyers or gliders (Khandelwal et al., 2023).

A fundamental component of snake flight is how effective the wing-body is at producing lift and drag forces, and how these forces change as the animal accelerates and the glide shallows (Socha et al., 2015). Understanding lift and drag production requires knowing the detailed kinematics of the morphing wing, as well as the lift and drag characteristics of the wing at each configuration. Given the orientation of the wing-body and the lift and drag coefficients of the cross-section, first-order approximations of the lift and drag forces can be calculated using the quasi-steady assumption (Ellington, 1984; Yeaton et al., 2020).

The quasi-steady assumption states that the forces acting on a wing are a function of its instantaneous speed and orientation to the flow; its acceleration and the time history of the flow are not considered. The validity of quasi-steady theory can be checked by examining whether the resultant aerodynamic force matches the center of mass acceleration of the animal. Although the quasi-steady assumption has proved insufficient to explain force production during flapping flight of insects and vertebrates (Hedenström et al., 2007; Warrick et al., 2005; Videler et al., 2004; Ellington et al., 1996; Muijres et al., 2008; Maxworthy, 1981; Dudley and Ellington, 1990; Dickinson et al., 1999; Chin and Lentink, 2016), it has not been tested explicitly for animal gliders (sensu stricto), with their nominally static body postures (Socha, 2011; Socha et al., 2015). Despite flying snakes displaying large postural changes involved in aerial undulation, a previous study found that undulation frequency was not significantly correlated to any glide performance variables, suggesting that undulation itself has a minor role in aerodynamic force production and that quasi-steady theory could be applicable (Socha and LaBarbera, 2005). Additionally, the fast forward speed compared with the relatively slow speed of undulation, quantified as the advance ratio, also suggests that quasi-steady theory might be applicable (Holden et al., 2014). However, the precise aerodynamic forces produced during snake flight have never been calculated because the detailed wing-body kinematics and the center of mass trajectory were unknown.

A quasi-steady analysis requires known lift and drag coefficients for the flyer of interest. The lift and drag characteristics of the 2D cross-sectional shape of Chrysopelea paradisi have been studied previously using load cell measurements, particle image velocimetry and computational fluid dynamics (Miklasz et al., 2010; Holden et al., 2014; Krishnan et al., 2014). The roughly triangular cross-section acts as a lifting bluff body (Figs 1D and 2A), with the concave ventral surface and protruding lips serving to improve force production (Miklasz et al., 2010). The shape produces appreciable lift over a large range of angle of attack, with the lift-to-drag ratio peaking at an angle of attack of 35 deg (Fig. 2A). The profile has gentle stall characteristics, with drag increasing slowly after 35 deg. The cross-section is also effective at producing lift at high angles of attack, as the lift coefficient at 60 deg is roughly the same as that at 15 deg. Additional physical modeling also determined that when two snake-like profiles are offset, there is a modification of the lift and drag through wake interaction effects (Miklasz et al., 2010). More recent physical modeling investigated such ‘tandem’ effects, finding appreciable lift-to-drag ratio changes for particular wing configurations of gap (horizontal spacing) and stagger (vertical spacing) that may be realizable during gliding (Jafari et al., 2021).

In addition to the lift and drag characteristics, quasi-steady theory requires knowledge of the body's kinematics: the instantaneous wing configuration, its orientation to the flow and each body segment's motion. Previous kinematic measurements of snake flight recorded either three or five landmarks on the body (Socha, 2002; Socha et al., 2005, 2010). Although these data enabled a detailed center of mass trajectory analysis, the time-varying body posture and its orientation to the airflow had to be estimated from the horizontal body posture and glide angle. However, in a real flying snake, different sections of the body might experience different angles of attack, and portions of the body are swept relative to the airflow, which will affect the forces and moments acting on the body. Additionally, the out-of-plane bending of the body causes the snake's wing-body to twist, thereby changing the local angles of attack (Yeaton et al., 2020). These complexities of the aerial snake made it impossible to estimate the locomotor forces from a few landmark points alone, and proved insufficient for guiding modeling studies. Additionally, the center of mass position, velocity and acceleration could only be estimated, as the true configuration of the body was not known.

Here, we tested the effectiveness of quasi-steady theory at approximating lift and drag forces along a flying snake's body during gliding. To do so, we did two things. (1) We conducted a new analysis of 3D kinematics data from our recent study (Yeaton et al., 2020). Specifically, we used 11–17 landmark points along the snake's body to recover the snake's time-varying body posture recorded in an experimental setting. This new analysis enabled us to estimate the center of mass position with high fidelity and subsequently calculate velocity and acceleration. Knowing the mass of the snake, we then used the center of mass acceleration to calculate the total time-varying aerodynamic force on the snake, which we took as the ‘true’ value. (2) Then, using as few assumptions as possible, we estimated the time-varying aerodynamic force produced across the entire body as predicted by quasi-steady theory. To do so, we used the 3D kinematics to define the local wing orientation, and calculated quasi-steady lift and drag forces using the previously measured 2D force coefficients, the blade element method and simple sweep theory. This calculation of forces is a sum of infinitesimal elements, which analogizes the snake as a collection of particles that connect with one another, but do not behave as a rigid body. We then tested the quasi-steady assumption by comparing the absolute errors between the center of mass force (‘true’) and the estimated theoretical aerodynamic force (‘model’). This study provides the first detailed analysis of the total aerodynamic load on flying snakes during gliding. Additionally, we used the kinematics data to estimate how Reynolds number, angle of attack and sweep angle vary along the body, which can be used to inform future physical and computational modeling of snake flight.

This study newly analyzes 3D kinematics data from an experiment originally presented in Yeaton et al. (2020). The main focus of that paper was to develop a quasi-steady analysis of gliding to examine the role of aerial undulation on gliding. Data related to undulation were presented therein, and the remainder of the un-analyzed kinematics data are addressed in the current study. A total of 43 glide trials from 7 snakes of the species Chrysopelea paradisi are presented here. We refer the reader to the previous paper for the full details of the experiment; for convenience, a brief summary is provided below.

We recorded the position of between 11 and 17 infrared (IR) tape markers placed along the dorsal surface of flying snakes using a 23-camera motion capture system with a sampling rate of 179 Hz. Experiments were conducted in ‘The Cube’, a four-story-tall black-box theatre located in the Moss Arts Center at Virginia Tech, modified for use as a large indoor glide arena (Fig. 1A–C) under IACUC protocol #15-034. Snakes were allowed to jump and glide under their own volition from a height of 8.3 m. From the landmark points, we exported the marker trajectories and filled gaps using Kalman filters. We then smoothed the marker time series using two passes of a second-order Butterworth filter, with cutoff frequencies selected separately for each coordinate time series for each marker (Winter, 2009); these values ranged from 7 to 17 Hz, which are greater than the nominal undulation frequencies observed in gliding snakes (∼1–2 Hz).

To reconstruct the time-varying 3D body, we rotated the previously determined snake airfoil segment (Socha, 2011) to lie within the local – plane at each position along the body. The segment is then scaled to a specified width [which is also the chord length, c(s)] based on its arc-length distance along the body, s (Fig. 3F). Qualitatively, the resulting wing-body model is visually very similar to the observed body while gliding.

For each trial, we compared the non-dimensional force errors as a function of the height of the animal's center of mass above the ground, as this proceeds with the progress of the glide in a manner which can be compared across all trials. Lastly, linear interpolation was used to interpolate the accelerations and forces such that each glide was sampled on a uniform height grid. This enabled us to calculate the averages and standard deviations of the forces, accelerations and errors for each individual.

In our aerodynamic analysis, we focused on quantifying the discrepancies between our model predictions and experimental observations. Specifically, we observed consistent deviations in the force estimates in the forward (⁠⁠) and vertical (⁠⁠) directions. Our model tended to overestimate drag force and underestimate lift force when compared with experimental data.

To quantify these discrepancies, we introduced time-varying multiplier factors b_L(t) and b_D(t) for the lift and drag forces, respectively. These factors represent how much we need to adjust the model's force estimates to align them with experimental observations, at each time t. Ideally, if the model perfectly matched the experiments, these multiplier factors would be 1, indicating no discrepancy.

Gap and stagger are defined as the relative horizontal and vertical position of the downstream body segment relative to the upstream segment. To quantify the relative spacing of the body for possible wake interaction effects, we calculated the time-varying gap and stagger values of the snakes throughout the trajectory. Gap and stagger are calculated in the trajectory reference frame, defined by rotating the body from the straightened frame such that the center of mass velocity is in the forward direction. In the trajectory reference frame, gap is defined as the horizontal distance between the upstream and downstream segments and stagger as the vertical distance (Fig. 2C).

Once in the trajectory reference frame, we found intersections of the body spline with a vertical plane at the center of mass (Fig. 2D,E). The locations of the intersections defined where the gap and stagger values were calculated. For each intersection, we took the relative displacement between the upstream and downstream airfoils as the gap and stagger, respectively. Sometimes, the body intersected the vertical plane multiple times. During these instances, we defined two separate gap and stagger values from the first and second intersections and the second and third intersections.

Analysis of the center of mass trajectory revealed steep glides compared with previous kinematic measurements of snake glides (Socha et al., 2010, 2005). This observation prompted us to revisit data from these other studies to extract gap and stagger information. For this analysis, we used data from Socha et al. (2010), which were the previous best kinematics data. This previous experiment recorded late-phase gliding in two individuals, with four trials per snake analyzed, of glides originating from a height of 15 m. The data comprise trajectories of five landmark points (head, 1/4 SVL, 1/2 SVL, 3/4 SVL and vent), rotated into the trajectory reference frame. Gap and stagger values were estimated from these data as described above, using marker pairs – head–1/2 SVL, 1/4–3/4 SVL and 1/2 SVL–vent – at time points when the lateral displacement of the marker pairs was the same. Time points with the same lateral displacement are an approximation to intersections of the body with the vertical plane passing through the center of mass; this approximation was made because the center of mass is not precisely known from the five-point data.

The analysis presented in this study is based on several layers of assumptions, beginning with the measured infrared marker trajectories. From the marker trajectories, the spline curve was used to represent the body as it moves through space. Next, the mass and width distributions were overlaid on the spline; the mass distribution was measured from anatomical snake sectioning, and the width distribution from photographs (Yeaton et al., 2020). The gap and stagger calculations were based on the spline, while the relative angle-of-attack and sweep angles of the wing segments both incorporate assumptions discussed below.

A necessary assumption for the quasi-steady force calculations is the orientation of the airfoil coordinate system. The airfoil coordinate system directly affects the aerodynamic force calculations, because it is one of two components that determine the angle of attack (the other being the average forward velocity vector of the snake along its trajectory). The need for an assumption arises from the fact that our kinematics data are limited: the IR system is assumed to have identified the centroid of the marker on the snake, which is along the dorsal backbone. As such, we have no quantitative information about the lateral positions (i.e. the width) of the snake's body. As an analogy, this situation is akin to identifying the root and tip of a bird's wing during flight, but not knowing the location of the leading and trailing edges. Without that information, it is impossible to specify the local angle of attack on any part of the snake's body. We selected a method to calculate the airfoil coordinate system that resulted in a reconstructed wing-body that was qualitatively similar to images of the snakes while airborne (Fig. 3G,H). The similarity included the twisting of the body at the ‘U’-bends, the overall orientation of the straight segments, and the calculated forces varying continuously along the body. Because we lack the full kinematic data, we have no way of calculating the error in angle of attack in our model, but the visual similarity of model to imagery of the real gliding snake provides some degree of confidence that the model can be informative and provide a basis to build upon.

The other aerodynamic assumptions are (1) extrapolation of lift and drag coefficients between angles of attack of 60–90 deg, (2) the use of simple sweep theory to account for the non-perpendicular body segments, and (3) quasi-steady aerodynamics. The lift extrapolation is based on the fact that the snake's cross-sectional shape is left–right symmetrical, and this shape at 90 deg will produce zero lift. The drag extrapolation was based on the similarity of the shape to other triangular shapes, with experimental coefficients reported to be around 2 (Hoerner, 1965). The simple sweep theory assumption enabled the use of our current best understanding of flying snake aerodynamics, as no data exist on the relative orientation of the airflow to the body. Likewise, the quasi-steady aerodynamics theory had not been tested on flying snake locomotion given our previously poorer understanding of the whole-body motion during a glide. Overall, the data presented in this study should help guide future work to address the assumptions used here.

The center of mass trajectory from the side and overhead views is shown in Fig. 4A,B for 43 trials from seven individuals. For the center of mass analysis, multiple glides from three individuals are shown separately, chosen because these individuals gave the greatest number of usable trials to analyze, and happened to span the observed mass range of animals tested. The side views indicate relatively steep glides compared with previous studies. The largest animal (snake 81, mass: 107.2 g) only covered about 4 m horizontally from an 8.3 m jump height. The lightest animal (snake 95, mass: 37.3 g) had more variation in the horizontal glide distance, but generally performed better and covered horizontal distances up to 5 m. Trial 618 (Fig. 4Aiii,Biii) was the best-performing glide, as it shallowed fastest and covered the most horizontal distance. However, the snake landed on the target tree placed in the glide arena (Fig. 1A), so the maximum potential horizontal distance was not recorded; it is estimated that this snake could have glided 7 m horizontally. The overhead view of the center of mass trajectories (Fig. 4B) shows that the glide paths were generally not straight. The center of mass motion has broad arcs, but no discernible turning events were observed. The overhead view shows no oscillations of the center of mass due to undulation.

The dimensional velocity polar diagrams (Fig. 4C) indicate initial forward velocities at or below 2 m s⁻¹ for all snakes. The velocity trajectories initially move vertically downwards on the diagram as the snakes accelerated downwards while keeping the horizontal velocity constant during the ballistic phase. The glide angle during this phase increases from roughly 0 deg to 60–75 deg before the velocity trajectory arcs upward. The shape of the arc upwards is similar across all glides and is consistent with motion onto the ‘terminal velocity manifold’ in theoretical models of gliding (Yeaton et al., 2017; Nave and Ross, 2019). For the smallest snake (Fig. 4Ciii), the trajectories curve upwards the most, indicating it progressed farthest through the glide. The curve upwards is more apparent when viewing the non-dimensional and rescaled velocity polar diagrams (Fig. 4D), as all effects of animal size have been removed. The glides that performed best, trial 618 from snake 95 and trial 505 from snake 88, progressed farthest on the velocity polar diagram and had higher initial horizontal velocities. The higher velocity likely caused the glide to transition more quickly, while slower horizontal velocities took longer to transition.

The quasi-steady forces and center of mass accelerations are shown in Fig. 5 for the three individuals with the greatest number of usable glides. If quasi-steady theory properly accounts for the aerodynamic forces, the left and right plots in Fig. 5 should be the same. We see the same trends across individuals in regard to the forces and accelerations, but quasi-steady theory is insufficient to explain the aerodynamic forces on the snake. In general, the force estimates are smoother than the accelerations, with smaller standard deviations during the trajectory. The lower standard deviation is due to only needing to perform one numerical derivative to calculate velocities, while the center of mass acceleration requires two derivatives, as well as integrating the forces having a filtering effect.

Both the average lateral force and lateral accelerations are near zero, indicating fairly good agreement between the accelerations and quasi-steady theory (Fig. 5A). The deviations are larger about the forward and vertical directions. The forward forces, which are the aerodynamic forces responsible for horizontal motion over the ground, are too low compared with the accelerations. The average accelerations peak at approximately 0.5 body weight (BW), while the forces peak at approximately 0.25 BW (Fig. 5B). By contrast, the vertical forces are higher than the accelerations. The vertical forces are responsible for offsetting the animal's weight against gravity. The time series of force indicate that more vertical force is produced than is needed to support the weight (forces are above the dashed line in Fig. 5C). The accelerations are lower and only cross the supporting body weight threshold for the lightest snake. Therefore, the quasi-steady force estimates show qualitatively the same trends as the accelerations, but the forces are underestimated in the forward direction and overestimated in the vertical direction.

Across individuals, the accelerations about the lateral direction are similar. However, there is a qualitative difference in the forward and vertical directions between the lightest and heaviest individuals (snake 95, mass: 37.3 g and snake 81, mass: 107.2 g) with a difference in mass of approximately 2.9×. The forward acceleration of the lighter snake peaks and then decreases near the end of the glide, whereas the heavier snake does not show this behavior. The vertical accelerations show the lighter snake accelerating upward near the end of the trajectory, whereas the heavier snake does not. The medium-heavy snake (snake 91, mass: 71 g, mass ratio of 1.9× compared with the lighter snake) appears to be accelerating slightly upwards at the end of the trajectory (Fig. 5C).

The absolute force errors, defined as the absolute difference between the calculated center of mass acceleration and the resultant aerodynamic force, are show in Fig. 6A for all of the seven snakes analyzed. Errors are near zero about the lateral direction. The forward force error is approximately 0.25 BW, while the vertical force error is approximately −0.5 BW. The physical interpretation of these errors is that quasi-steady theory overestimates drag and underestimates lift, as the drag force acts predominantly in the vertical direction, while the lift force acts predominantly in the forward direction. Decreasing drag and increasing lift will rotate the resultant aerodynamic force vector forward and reduce the error.

The required multipliers to the lift and drag forces to reduce the total error to zero are shown in Fig. 6B. The lift force needs to be multiplied by 1.35× and the drag force needs to be multiplied by 0.6× the quasi-steady values. Note that the change in drag is approximately constant throughout the glide, while the change in lift decreases from approximately 2× to 1× as the glide progresses. The force multipliers thus indicate that quasi-steady theory consistently overestimates the total drag force by approximately 40% and underestimates the total lift force by approximately 35%.

The time-varying lift and drag distributions for the best-performing glides are shown in Fig. 7 for different heights as the animals glide. The force distributions have the same orientation as predicted by quasi-steady theory, but have been scaled using the lift and drag multipliers such that the total force errors are near zero. Therefore, the true local lift and drag force may be quite different from those shown. By 10% through the glide (non-dimensional height z/h₀=0.9, the first row in Fig. 7), the snake bodies have already formed the ‘S’-shaped glide posture, and lift and drag forces are being produced. As the glides progress, the smaller animal (snake 95, mass: 37.3 g, SVL: 64.4 cm) forms a tight body posture, with the posterior body dropping below the head. The larger animal (snake 88, mass: 71.9 g, SVL: 88.8 cm) has a more open body posture, with fewer spatial periods of bending, and a more horizontal orientation.

There are commonalities in the force distributions for both animals. The drag force is continuous along the body, including along the straight segments and along the ‘U’-bends. By contrast, consistent with the simple sweep theory assumption, the lift force decreases to zero at the ‘U’-bends, and is largest along the straight segments. The transition region to zero lift is small and does not extend over the entire ‘U’-bend region. The small transition region is likely related to low sweep angles before and after the ‘U’-bends (Fig. 8E,F). The estimated force produced by the tail is small because of its small width. This small width also results in lower Reynolds numbers (Fig. 8A,B). For snake 88, the body is initially highly swept, with the straight segments angled roughly 45 deg relative to the forward direction. However, the model still indicates that these areas produce appreciable lift force. The side views of the glides indicate that the drag force acts in the vertical direction and against the direction of forward motion. The lift force is angled upwards relative to the horizontal and along the direction of forward motion.

For the glides shown in Fig. 7, the time histories of the Reynolds numbers, angles of attack, sweep angles and dynamic pressure fraction distributions are shown in Fig. 8. The Reynolds number (Fig. 8A,B) broadly tends to increase throughout the glide because the snake's airspeed increases as it accelerates. The larger animal has higher Reynolds numbers due to its greater width and speed. The Reynolds numbers peak midway along the body, where the animal is widest. The angles of attack are high, ranging from 60 to 90 deg at various points along the body. The ‘U’-bends have the highest angles of attack (shown as dashed lines in Fig. 8A,B), which results in the locations of zero lift production in Fig. 7. The ‘U’-bends have the highest sweep angles, and locations along the body near the ‘U’-bends have low sweep values. Even along the straight segments, the sweep angle is generally greater than 30 deg. Lastly, the dynamic pressure fraction is shown in Fig. 8G,H, with the 75% contour highlighted. Values of 100% indicate no decrease in lift or drag due to the swept wing, whereas values of 0% indicate no force being produced because of sweep. The straight segments of the snake maintain the most dynamic pressure, although there are regions where this is not the case (Figs 8G and 7 for snake 95 at 0.7 and 0.6 of the height fallen).

The relative spacing of the perpendicular body segments for all glides and at different progressions through the glide are shown in Fig. 9A and summarized in Table 1. Initially, the distributions of gap (horizontal spacing) and stagger (vertical spacing) show high spread, which decreases later in the glide. During the first quarter of the glides, the median gap is 3.0 c (chord) and the median stagger is 3.7 c. The second quarter of the glides exhibits the same median stagger of 3.7 c, but the median gap decreases to 1.0 c. During the third and fourth quarters of the glides, the median gap is effectively 0 c, while the median stagger is 4.2 c. As a measure of distribution spread, we use the interquartile range (IQR). Initially, the gap and stagger IQR is 3.8 c and 2.4 c, which decreases to 2.7 c and 1.8 c during the second quarter of the glide. The IQRs decrease further to 2.1 c and 1.7 c, 2.0 c and 1.3 c during the second half of the glides. The gap and stagger distributions indicate gaps that are near zero and even negative. A gap of zero indicates the forward airflow contacts the anterior and posterior body simultaneously, while a negative gap indicates the rear airfoil leads the front airfoil. Staggers are generally positive, indicating that the rear airfoil is below the front airfoil relative to the airflow.

Overlaid on the joint distributions are the gap and stagger measurement locations for wake interaction effects (Jafari et al., 2021). Gap and stagger combinations are observed in the aerodynamic interaction region, although these combinations are relatively rare. Wake interaction effects are most prominent along the top row, with the rear airfoil directly behind the front airfoil. There are also wake interaction effects along the second row, but the effect is smaller. The gap and stagger measurement locations initially overlap with the observed gap and stagger distributions, but the overlap is less later in the glide, as the observed distribution moves to a gap of zero.

The angles of attack of the front and rear airfoils are shown in Fig. 9B. The angle of attack of the rear airfoil is correlated with the angle of attack of the front airfoil. Initially, angles of attack are very high, upwards of 70–90 deg, with the median angle of attack being 74 deg for both the front and rear airfoils. As the glides progress, the spread of the angle of attack distribution decreases, and the median angle of attack ultimately decreases to 56 deg for the front airfoil and 65 deg for the rear airfoil. All angles of attack were high, indicating large drag coefficients and small lift coefficients (Fig. 2A). Unlike the angles of attack, the sweep angles of the front and rear airfoils (Fig. 9C) do not follow a clear trend. Initially, there is a large cluster of high sweep angles, which may be due to the jump and body formation phase of the glide. Later in the glide, the higher sweep angles are not seen, but there is no clear trend relating the sweep angle of the front and rear airfoils.

Lastly, gaps and staggers from this study and estimated gaps and staggers from a previous study with only five markers on the snake (Socha et al., 2010) are shown in Fig. 9D. The gaps from the five marker trials are much greater than observed in this study. The steep glides from this study have a gap of approximately 0 c, while the shallow glides from the previous study have a gap greater than 5 c, although the staggers are similar. The gap and stagger combinations from the five-point trials are generally in the wake interaction region, although few time points from the current study are located in the wake interaction region. The right panel of Fig. 9D indicates that although uncommon, gap and stagger configurations in the wake interaction region are observed in the present study.

This study of flying snake glide trajectories provides the first measurements of the snake's center of mass informed by the whole body. The overhead and side views of the glide trajectory provide our first measurements of the center of mass motion of flying snake glides. The recorded glides were steeper than would be anticipated from previous studies, with a horizontal glide distance of only 4 m from a jump height of 8.3 m, whereas maximum glides from a previous experiment averaged 10.1 m horizontally from a jump height of 9.6 m (Socha et al., 2005). There are several possible explanations for this difference, including the more massive snakes used in the present study, physical markers being placed on the body, lack of visual cues from the indoor glide arena, potential colony effects, or because wild-caught snakes are more accustomed to gliding. Alternatively or additionally, only the best-performing glides from previous studies were analyzed in detail, whereas all trials, regardless of performance, were analyzed in the present study. The larger mass directly affects the wing loading of the animal, with predictable results in the velocity polar diagram (see Fig. 4). To address visual cues, a target tree was placed in the arena (see Fig. 1), as has been done for previous studies. There was no discernible difference when handling the animals and encouraging them to jump. Additionally, the physical tape markers, as opposed to painted markers, were placed carefully on the animal's dorsal surface to minimize disturbance to their flattening. As a control, one snake was left unmarked and did not show noticeably different glide behavior or performance from that of the marked snakes. Although the glides were steep, the increased spatial and temporal resolution enabled us to fully quantify the position of the body in much greater detail than in previous studies. This improved understanding enabled us to calculate the center of mass, as well as estimate aerodynamic forces, and measure the gap and stagger of different body segments.

The overhead view of the glide path (Fig. 4B) did not show obvious center of mass deviations related to undulation. These data show the unfiltered center of mass, so the smoothing effects of digital filtering are not present. We did see broad arcing turns, although no distinct turning events were observed during the trials. This arcing motion could be due to the initial jump conditions and result from stabilization of the rotational motion (Yeaton et al., 2020). There is one trajectory that followed a fairly straight course, trial 413 from snake 91 (Fig. 4Bii), from the launch branch to the target tree. We recorded only a few glides that landed on the tree.

The velocity polar diagrams (Fig. 4C,D) have been discussed previously from a theoretical modeling perspective (see Yeaton et al., 2017; Nave and Ross, 2019). The empirical velocity polar diagrams provided here support several theoretical predictions from those previous studies. First, the wing loading rescaling enables comparisons between individuals that vary in mass. The theoretical studies predicted that smaller individuals would progress further through the diagram during a glide and show better glide performance, which was seen in this study, with individuals varying by a factor of 3 in mass. Additionally, a greater initial forward velocity was predicted to more quickly reach a steady-state glide and perform better. Here, we found the best-performing glides did indeed have greater initial forward velocities. All diagrams show a characteristic turn that signifies the transition from the ballistic to the shallowing phases of gliding. The theoretical model predicted that a trajectory in the velocity polar diagram, starting from a horizontal jump, will accelerate quickly downwards before curving upwards, transitioning onto the terminal velocity manifold, a higher-dimensional analog to the terminal velocity speed. This transition point occurs for different absolute airspeeds for different individuals, However, the rescaled diagrams show that the transition occurs for a non-dimensional speed of approximately 1, which is approximately the terminal velocity speed.

Quasi-steady theory was insufficient to explain the aerodynamic forces produced during short glides in flying snakes. When compared with the experimentally determined center of mass acceleration, the predicted quasi-steady lift forces were 35% lower and the quasi-steady drag forces were 40% higher. We found that the difference in drag was constant throughout the glides, whereas the difference in lift decreased throughout the glide (see Fig. 6).

One possibility for the lift multiplier changing with progress through the glide is that initially the angles of attack along the wing-body are high. As the glide progresses, the angles of attack decrease from 70–90 deg at the start to 50–70 deg at the end of the glide. The lift and drag curves are substantially different at high angles of attack; the lift coefficient is near zero, while the drag coefficient is at a maximum (see Fig. 3A). The relative change in the lift and drag coefficients also varies as the angle of attack is decreased. From 90 to 60 deg, the extrapolated lift coefficient increases from 0 to about 1, while the extrapolated drag coefficient only decreases from approximately 2 to 1.6. Therefore, the modification to the lift force is more sensitive during the initial portion of the glide as the angle of attack decreases from a high value. While the lift and drag extrapolated values at α=90 deg are based on theory for geometrically similar airfoils (Yeaton et al., 2020), the true lift and drag values in this regime of angle of attack are unknown, as previous experiments only measured up to 60 deg. The most likely difference in lift or drag coefficients in this regime is likely to be in drag, because the estimate at 90 deg is less certain than that in lift (which is zero at 90 deg as a result of the symmetry of the airfoil with respect to the airflow at this configuration). However, the results of this study show that the quasi-steady model over-estimates drag, suggesting that our extrapolation was not the major factor in the difference. More generally, it is likely that the discrepancy between model and reality stems from the nature of the kinematic data: in the trials analyzed here, the snake was not only undulating but continuously accelerating throughout the trajectory, yet the model takes a quasi-static approach, which intrinsically does not include acceleration. It would be interesting to apply the model in future studies to non-accelerating parts of a glide, which could in theory be obtained from snakes launched from a greater height. Such work would build on recent theoretical advances of our understanding of gliding across all animals including flapping flyers (Cheney et al., 2020; Usherwood et al., 2020; KleinHeerenbrink et al., 2022), beyond the gliding-restricted taxa such as flying snakes.

Tandem airfoil effects, due to wake interactions between a leading airfoil and a trailing airfoil, are an unsteady aerodynamic effect. They were considered here as a potential reason for the insufficiency of quasi-steady theory. The distributions of gap and stagger, and how they change depending on progress through the glide (Fig. 9A), revealed zero or even negative gap values at the end of the glide. The gap and stagger results indicate that flying snakes use a wide range of body configurations relative to the airflow. Notably, during all glides, the gap and stagger values spend little time in the region where wake interaction effects are greatest (Fig. 9D), from which it can be concluded that this effect is not the reason for the insufficiency of quasi-steady theory.

While the configurations may not be advantageous from an aerodynamics and glide performance perspective, they may serve ecologically relevant functions. For example, snakes executing short glides, or falling near vertically, may do so to escape predators. Falling straight downwards still requires controlling aerodynamic and inertial forces to ensure stability, but will not result in substantial horizontal travel. Steep glides may allow the snakes to fall away from the predator while staying on the same tree, or in the nearby area.

Flying snakes therefore have a large performance envelope within which to operate. Some glides can be viewed as predominantly falling, in which the animal moves vertically downward or only covers a few meters horizontally, whereas other glides cover significant horizontal distance, possibly to escape to a different tree. It is known that other gliding animals (Khandelwal et al., 2023), such as flying squirrels and Draco lizards, have higher shallowing rates and shallower glides than flying snakes (Socha et al., 2015), but they may not be able to execute very steep glides. The morphing wing-body of flying snakes may be particularly well suited to stabilize short glides, as it possibly allows the animal to correct for rotational torques. Additionally, the multi-wing configuration may enable wake interactions that increase glide performance (Fig. 9D), although this may not be volitional and may simply be an artifact of the shallowing glide.

Future studies are needed to address the assumptions used in this work. The biggest kinematic assumption is the airfoil coordinate system that is overlaid on the spline fit of the body. This assumption was used because the orientation of the body was not resolvable from the marker time series. The airfoil coordinate system directly affects the force estimates, as well as angle of attack and sweep angle estimates. Future experiments, with greater camera coverage, may be able to measure body orientation with higher resolution.

The next assumptions to be addressed are quasi-steady theory and simple sweep theory. Both assumptions were applied so that previously measured lift and drag coefficients could be used. An alternative to using quasi-steady force coefficients and simple sweep theory would be to perform anatomically accurate computational fluid dynamics (CFD) simulations, with a moving mesh derived from the kinematics analysis. This analysis would account for unsteady fluid phenomena, and provide details about wake interaction, vortex shedding and flow at the ‘U’-bends of the snake. This analysis would also provide insights into span-wise flow along the straight body segments. In fact, a full-body CFD analysis has been conducted recently, but it treated the snake in an idealized fashion, with the motion of the body sinusoidal, confined to a single, horizontal plane, in a non-accelerating regime (Gong et al., 2022). This CFD study provides new insight into potential aerodynamic mechanisms including enhancement by leading and trailing edge vortices, and it would be useful to conduct a similar force comparison with experimental data as done here. As an alternative to computational models, physical or robotic models could be tested in air or water tunnels, as in Holden et al. (2014), or in a tow tank, in which the model could be accelerated. Flow and force measurements can then be taken of the whole snake model, or of straight segments at different attack and sweep angle combinations (Fig. 8C,D). Finally, it would be informative to compare model results with trajectories that include shallower gliding seen in prior studies (Socha, 2002; Socha et al., 2005, 2010; Socha and LaBarbera, 2005).

In some ways, the reciprocal progression of experiment and model development suggested by this study recapitulates the early history of flight studies of birds, bats and insects (Alexander, 2003; Biewener and Patek, 2018; Dudley, 2000). Each of these groups shared common features of locomotion, but more detailed understanding of kinematics and aerodynamics was required to refine the theoretical understanding of their flight. In a parallel way, this study similarly represents a foundational step forward in understanding how flying snakes glide, providing new inspiration for future work.

Using a new kinematic analysis of the center of mass trajectory and body orientation, we tested whether quasi-steady theory can predict the time-varying lift and drag forces on flying snakes. Quasi-steady theory was insufficient to fully explain the forces acting on the body. Our results indicate that unsteady and 3D aerodynamic effects are likely important for snake flight, even during short glides. During some glides, we did find body configurations where the anterior body was located where unsteady wake interaction effects are possible, but this is unlikely to explain the shortcomings of the applied quasi-steady theory. However, wake interaction effects may be more prominent during late phase gliding. The time-varying body posture and orientation can be used in future computational dynamics and physical modeling studies to elucidate unsteady and 3D aerodynamic phenomena during snake flight.

We thank the Virginia Tech Institute for Creativity, Arts, and Technology (ICAT) for facilities access, crew support and funding. In particular, we thank Tanner Upthegrove, Nate McGowan and Benjamin Knapp. We thank Gary Nave, Talia Weiss, Michelle Graham, Jack Whitehead, Joel Garrett and Grant Baumgardner for assisting with the kinematics experiments. We thank Hodjat Pendar for help with developing the mathematical model of snake flight. We thank Barry Robert, Julie Settlage and Amy Rizzo for help with IACUC. Portions of this paper are reproduced from I.J.Y.’s PhD dissertation (Yeaton, 2018).