The aerodynamics of flying snake airfoils in tandem configuration

Snakes in the Southeast Asian genus Chrysopelea have evolved the ability to glide (Socha, 2002, 2011; Socha et al., 2015). As they lack specialized surfaces to produce aerodynamic forces, they instead use their entire body as a wing. Once airborne, the snake splays its ribs to transform a rounded cylindrical body to a roughly triangular cross-sectional shape (Fig. 1A), which has the effect of increasing the projected area and lowering the profile in the dorsoventral axis. In addition to this body flattening, flying snakes undulate laterally by sending traveling waves posteriorly down the body. At any instant, the snake appears to assume an S-like posture (Fig. 1B), but the undulating motion continuously repositions the body, with substantial displacement in the vertical axis (Fig. 1C). As the snake moves through its glide trajectory, its glide path shallows, and the relative orientation of the oncoming airflow changes accordingly. Although the body can be considered as a single, reconfigurable wing, it can also be viewed as a composite of many segments, whose specific orientation and positioning relative to the airflow varies throughout the glide. In certain configurations, downstream segments may interact with the wake created by upstream segments. How does this continuous repositioning of the body during aerial undulation influence the snake's aerodynamics?

The cross-sectional shape of the body plays a major role in the snake's glide performance. In Chrysopelea paradisi, this shape is symmetrical in the fore–aft axis, with a semi-triangular dorsal surface, a flat ventral surface, and a pair of ventrally oriented ‘lips’ on each lateral edge (Socha, 2011) that appear to run along most of the length of the body between the snout and vent, producing a concave shape. Although this cross-section can be modeled as a bluff body, its aerodynamic performance is superior to some airfoils in steady flows at similar Reynolds numbers (Holden et al., 2014). Two studies provide our current understanding of the aerodynamics of this two-dimensional shape (Holden et al., 2014; Krishnan et al., 2014). Positive lift is created at angles of attack (α) greater than 10 deg, with lift coefficients ranging from approximately 1.0 to 1.5 at Reynold numbers between 3000 and 15,000. The lift coefficient tends to be slightly higher at lower Reynolds numbers within this range. Lift coefficients peak sharply at α=35 deg and gradually decrease up to the maximum studied angle of 60 deg. Drag coefficients also increase dramatically after α=35 deg. Vortex trajectory analysis of the experimental velocity field and subsequent computational investigations (Krishnan et al., 2014) suggest that vortex-induced suction on the dorsal surface of the snake may underlie the surprisingly high lift coefficients at 35 deg. Overall, near-maximum lift-to-drag ratios (L/D) ranging from 1.3 to 2.7 occur between α=15 deg and α=40 deg, and positive L/D occurs beyond α=10 deg, providing a physical rationale for the snake's glide performance.

However, the aerodynamic performance of single two-dimensional airfoils alone is not sufficient to explain the gliding abilities of real snakes. For instance, the shallowest glide angle reported for C. paradisi, 13 deg (Socha et al., 2005), indicates that flying snakes are capable of achieving an average L/D of 4.5, which is almost twice as large as the maximum value of 2.7 found in the modeling study by Holden et al. (2014). As first-order approximations, initial studies neglected a range of complexities of the system, including undulation, 3D shape, possible variation of the snake profile along the body, and unsteady effects. Clearly, some aerodynamic mechanisms that contribute to glide performance remain unexplained. Unsteady and 3D effects have previously been shown to augment lift in some biological systems (e.g. Bahlman et al., 2013; Wang et al., 2014; Shyy et al., 2009; Song et al., 2014; Muijres et al., 2008; Lentink and Dickinson, 2009; Bomphrey et al., 2006), but the snake's large advance ratio (the ratio of the forward motion to the reciprocating motion) compared with that of insects and birds (Vogel, 2003; Dickson and Dickinson, 2004; Holden et al., 2014) suggests that unsteady mechanisms are less likely to produce a significant aerodynamic contribution. Recent studies (Jafari et al., 2017; Yeaton et al., 2020) found that undulation does have a functional influence on glide performance, acting to increase the stability of the snake. However, the theoretical modeling of the snake only incorporated single-foil aerodynamic coefficients, and did not examine interactive aerodynamics.

Here, we focused on the aerodynamic effects of the relative positioning and orientation of the body. Specifically, we hypothesized that the staggered S-shaped configuration of the gliding snake produces a complex interaction between the fore and aft body segments, whereby wake formed upstream is intercepted by downstream segments. Such tandem interactions have been observed in flapping wings (particularly in dragonflies) with an often favorable contribution to the aerodynamics (Akhtar et al., 2007; Lehmann, 2008, 2009; Maybury and Lehmann, 2004; Usherwood and Lehmann, 2008; Wang and Russell, 2007; Warkentin and DeLaurier, 2007; Weimerskirch et al., 2001). Static tandem airfoils (with Re≈10⁴–10⁵) have been shown to experience mixed aerodynamic effects (Scharpf and Mueller, 1992; Michelsen and Mueller, 1987; Husain et al., 2005), demonstrating that tandem arrangements can also produce poorer performance. A preliminary study by Miklasz et al. (2010) on flying snakes suggested that tandem effects may be substantial in certain configurations. They found that certain configurations of gap (defined as the horizontal spacing difference) and stagger (defined as the vertical spacing difference) resulted in changes in L/D relative to a single airfoil, increasing as much as ∼50% and decreasing by ∼20%. However, that study used simple models that lacked an anatomically accurate shape, kept the angles of attack fixed at a single value (of 25 deg), measured forces only on one airfoil, and only explored a small portion of the gap–stagger parameter space.

In this study, we used anatomically accurate 3D-printed models to examine the interactions of two flying snake airfoils, representing upstream and downstream body segments. The angles of attack and the gap and stagger between the two airfoils were varied to reproduce a wide range of conditions experienced by the snake during aerial undulation. Load cells were used to measure the individual lift and drag forces acting on each airfoil, and in follow-up trials, digital particle image velocimetry (DPIV) provided flow fields to examine configurations of particular interest. By determining the physical bounds on the aerodynamic capabilities of the tandem airfoils, this study furthers our knowledge about how flying snakes produce forces for gliding and control, and allows us to predict the limits on the snake's aerial performance.

First, load cells were used to measure the lift and drag forces on each of the airfoils in the tandem arrangement. Next, we used DPIV to obtain the flow field measurements for select cases that exhibited significant lift enhancement compared with single airfoil models. The velocity field results were used to estimate the aerodynamic forces and verify the load cell measurements, to look at the dynamical interactions between the fore and aft sections, and to investigate the specific sources of the high lift coefficients.

In designing our experiments, we aimed to minimize the number of independent variables so that the number of test cases would not become impractically large. Previous measurements (Holden et al., 2014) showed that the aerodynamic performance of a single airfoil is only mildly dependent on the Reynolds number. So, we kept the Reynolds number constant at Re=13,000, which corresponds to an adult snake in a fully developed glide (Socha et al., 2005). We addressed the effects of variation in spacing (gap and stagger) between the airfoil sections of the snake, and the airfoil's angle of attack. Thus, the parameter space consisted of the distance between the airfoils parallel to the flow (gap, Δx), the distance perpendicular to the flow (stagger, Δy), and the upstream and downstream angles of attack, α_u and α_d. Sampling of this large parameter space entailed a trade-off between being able to capture the important physical phenomena and generating a feasible test matrix. To reduce the volume of the parameter space, we focused on the set of parameters that are most typical of snakes during flight based on the previously determined kinematic data (Socha et al., 2010). Guided by these data, the gap and stagger distances were varied over the ranges 2c<Δx<8c (with increments of 2c) and 0c<Δy<5c (with increments of 1c), where c is the chord length of the airfoil. As no data exist on the snake's exact angles of attack along the body, we used the range of glide angles across the whole body reported in Socha et al. (2005) as the basis for the range of angles of attack in the experiments: 0≤α_u, α_d≤60 deg. This approximation is justified by the kinematic data indicating that the snakes glide with the ventral portion of the body facing the ground (Socha et al., 2010). This means that typically the body cross-sectional shape is approximately aligned with the plane of mean body posture. We also considered the possibility that the fore and the aft body may not have the same angle of attack, but we assumed that the snake cannot exhibit extreme values of twist (Moon, 1999); therefore, we rejected combinations with difference of greater than 30 deg between α_u and α_d. The test matrix with a total of 192 (=4×6×8) trials is summarized in Table 1. The combinations of the angles of attack were sorted in two groups: (i) the upstream angle of attack was kept fixed at α_u=30 deg and the downstream angle of attack was varied from α_d=0 deg to 60 deg, and (ii) the downstream angle of attack was kept fixed at α_d=30 deg and the upstream angle of attack was varied from α_u=0 deg to 60 deg.

All experiments were conducted with models that were geometrically accurate replicas of the cross-sectional shape of the airborne snake (Fig. 1A). Each model was fabricated in a stiff acrylate photopolymer (VeroWhite, Sculpteo, San Leandro, CA, USA) using a material-jetting additive manufacturing system (Connex 350, Stratasys, Eden Prairie, MN, USA). The models had a chord width of 25.4 mm and an aspect ratio of 20 to simulate the approximately straight sections of the snake's body normal to the flow (Fig. 1B). The models spanned the entire width of the tunnel to eliminate end effects, as in Holden et al. (2014). The experiments were conducted in a closed-loop water tunnel at Virginia Tech with a 0.6×0.6×1.5 m test section. The freestream turbulence intensity at the test location was less than 3% at speeds up to 0.5 m s⁻¹ (Gifford et al., 2011). The speeds were chosen to maintain Reynolds number similarity to the range exhibited by the snakes in the air (Socha, 2011).

The experimental setup was a modified version of the one used in Holden et al. (2014). Two grids of holes were drilled in the sidewalls to accommodate the second model, either upstream or downstream from the first (Fig. 2A). With this configuration, the forces on one of the models were measured while the other model was fixed to the sidewalls. Four 5 lb (22 N) load cells (LCFD-5, accuracy: ±0.15% FSO, DMD-465WB, strain-gage amplifiers; Omega Engineering, Stamford, CT, USA) were used to simultaneously measure the vertical and horizontal forces on the model. LabVIEW software (National Instruments, Austin, TX, USA) and a DAQ board (NI PCI-6251) were used to record the data for each trial at a sampling rate of 1000 Hz for 30 s. The weight of the experimental setup and model was accounted for by zeroing the force measurement before each trial.

The flow field was spatiotemporally resolved using planar DPIV measurements. To minimize the undesirable optical effects in the captured images, the acrylic sidewalls of the force-measuring setup were modified to accommodate detachable walls of acrylic specifically cut for each combination of gap and stagger (Fig. 2B). The experimental models were rigidly mounted to the detachable walls for the DPIV measurements. To seed the water tunnel, small, neutrally buoyant, hollow glass spheres (mean diameter, 126.4 μm; Potters Industries, Valley Forge, PA, USA) were used to act as flow tracers. The particles were illuminated using a dual-head laser system (Pegasus PIV, New Wave Research, Portland, OR, USA). Optical elements were used to spread the laser beam into a thin sheet with an approximately 1 mm thickness. To eliminate the shadow cast by the models, two mirrors were placed in a windowed box submerged into the water above the experimental setup to reflect the laser plane back through the test section. The desired region of interest included the two models, the distances between them (maximum of 6 chords of gap and 3 chords of stagger; Table 1), and the downstream wake. We used two cameras to capture this large region with sufficient resolution, including adequate overlap between the field of view of the two cameras. To capture the entire height of the desired field of view, the cameras were placed at one chord length vertical offset for the cases with stagger of Δy=2c or 3c, and at the same level for the cases with stagger of Δy=0c or 1c (Fig. 2C). High-speed video cameras (XS-5, Integrated Design Tools, Tallahassee, FL, USA) were used to record 1024×1280 pixel images at a sampling rate of 1000 Hz. Using 105 mm lenses with an aperture of f/2.0, an average image magnification of 108.1 μm/pixel was achieved. Before taking each set of data, a calibration grid (an aluminium plate with a grid of 3.2 mm dots at 12.7 mm spacing) was placed at the laser sheet plane and recorded. The calibration data were used to correctly overlap and align the images of the two cameras.

The time-resolved DPIV images were processed with the publicly available code PRANA (Eckstein and Vlachos, 2009a,b; Eckstein et al., 2008). We used the Robust Phase Correlation (RPC) algorithm with three iterations of multigrid iterative deformable windows (Scarano, 2001) to correlate the images and determine the instantaneous velocity fields. A window deformation algorithm was used as it is known to produce higher accuracy velocity field measurements by preventing the loss of particle-image pairs as the result of in-plane shear (Scarano and Riethmuller, 2000). Gaussian image masking was applied to the image windows to decrease induced high-frequency noise in the cross-correlations. Between passes, spatial filtering of the velocity fields was used to stabilize the iterative process (Schrijer and Scarano, 2008). The DPIV processing was carried out with three passes, one pass with 64×64 pixel interrogation windows, and two passes with 32×32 pixel interrogation windows. After each pass, the results were validated using velocity thresholding and universal outlier detection to replace vectors calculated from incorrectly matched correlation peaks (Westerweel and Scarano, 2005). The final processed velocity fields had a uniform vector grid spacing of 8 pixels.

Next, we used proper orthogonal decomposition (POD) with the method of snapshots (Smith et al., 2005; Berkooz et al., 1993) to decompose the dataset into orthogonal basis functions that capture the energetically important modes of the data. As the lower energy modes were more likely to be overwhelmed by experimental noise, the influences from the noise could be minimized by removing those modes (Charonko et al., 2010). Here, we reconstructed the velocity field by keeping the modes that accounted for 80% of the total energy.

Contour plots of the overall L/D of the tandem airfoils in the Δx−Δy plane show that when (α_u=20 deg, α_d=30 deg) or vice versa, the tandem L/D was generally larger than other combinations of the angles of attack, and it was fairly insensitive to the relative position of the airfoils (Fig. 4C,D). At these combinations of angles of attack, the maximum value of the tandem L/D exceeded 2.2 and was almost 10% larger than the maximum lift-to-drag ratio of a single airfoil (L/D=2 at α=0 deg) and 16% larger than L/D_single (the average of the lift-to-drag ratios of single models with the same angles of attack, found using Eqns 5a and 5b). With the other combinations of the angles of attack, the tandem configuration generally had a positive effect; however, the maximal L/D, if present, was confined to a few specific postures (e.g. the case with α_u=30 deg, α_d=60 deg; Fig. 4G). We also observed some negative effects; for example, a 12% decrease compared to L/D_single with (α_u=60 deg, α_d=30 deg) (Fig. 4H).

To investigate the physics underlying the aerodynamic properties of the tandem model, DPIV flow field data were collected for a select few configurations. As indicated in Fig. S1 by circled data points, the DPIV configurations included those showing significant changes in the tandem L/D with respect to single airfoil data along with some of their ‘adjacent’ positions without a large tandem effect. For instance, by just changing the downstream angle of attack, the significant tandem effect in one configuration (Δx=6c, Δy=2c, α_u=30 deg, α_d=0 deg) nearly vanished with α_d=20 deg. In addition to visualizing the flow dynamics, DPIV data were used for calculating the lift and drag coefficients of each airfoil in order to verify those obtained from the load cell data. Fig. S2 shows that, in general, the results obtained from the two methods agreed well, with the DPIV results mostly lying within the 95% confidence interval of the load cell data, and with the maximum relative error being about 18%.

The two important subsets in the DPIV configurations were: (i) the biologically relevant cases – with (α_u=30 deg, α_d=0 deg), the largest increase in L/D was observed at 6 chords (6c) of gap and 1 to 2 chords (1–2c) of stagger, intriguingly coincident with the most probable relative position of the fore and the aft body in airborne snakes (Socha et al., 2010); and (ii) the extreme tandem effect cases – with (α_u=30 deg, α_d=0 deg), and two chords (2c) of gap, only one chord of change in stagger from Δy=0c to 1c caused by far the biggest change in the lift-to-drag ratio. Table 2 summarizes the individual lift and drag coefficients of the two airfoils for these two groups of configurations. Within each group, the upstream and downstream angles of attack and the gap were the same, but the stagger varied. For the biologically relevant group, the upstream airfoil (30 deg angle of attack) showed similar aerodynamic performance to a single model at 30 deg angle of attack. By contrast, the downstream airfoil produced substantially larger lift and smaller drag with Δy≤2c; however, the enhanced performance was lost at Δy=3c. For the extreme tandem effect group, at Δy=0c, the lift and drag of both models changed with respect to their single model counterparts, and the overall lift nearly summed to zero, resulting in a drastic negative change in the combined lift and combined L/D (Fig. 3; Fig. S1). At Δy=1c, the upstream airfoil with 0 deg angle of attack produced positive lift, as opposed to the negative lift of a single model, making the most important contribution to the drastic enhancement of the aerodynamic performance. At Δy=2c, the upstream lift became negative again but remained larger than the single model lift; therefore, the tandem system experienced a smaller positive change in aerodynamic performance.

The vortex shedding behavior of the tandem airfoils was observed from the temporal spectra of the load cell measured forces. Fig. 5A,B shows the data for the biologically relevant configurations and a few of their adjacent configurations with no significant tandem effects. The upstream airfoils, all of which had an angle of attack of 30 deg, showed power spectral densities quite similar to that of a single airfoil at the same angle of attack (Holden et al., 2014) with a distinct dominant frequency at f*≈0.23. However, the behavior of the downstream airfoils varied with the tandem configuration. For the only two configurations with significant tandem effects (Δx=6c, Δy=1c,2c, α_u=30 deg, α_d=0 deg), the downstream airfoil deviated from its single model counterpart to be ‘locked’ to the upstream dominant frequency. The proximity of the upstream and downstream peak frequencies for (Δx=6c, Δy=2c, α_u=30 deg, α_d=20 deg) could not be interpreted in the same way because the peak frequencies of single airfoils at 20 deg and 30 deg angles of attack were already close to each other. In the cases (Δx=4c,6c, Δy=3c, α_u=30 deg, α_d=0 deg), the tandem arrangement had negligible influence on the downstream dominant frequency (f*≈0.34), which remained close to that of a single airfoil. For the last configuration (Δx=4c, Δy=2c, α_u=30 deg, α_d=0 deg), the downstream power spectral density exhibited both previously observed peak frequencies, which indicated some influence from the upstream wake without the vortex shedding frequencies being completely locked. The temporal spectra of the extreme tandem effect cases and their adjacent configurations (Fig. 5C,D) were less orderly, with the peak frequency varying from f*≈0.23 to f*≈0.42, and neither of the upstream and downstream groups showed a preferred dominant frequency. In particular, no clear peak could be identified for the case (Δx=2c, Δy=0c, α_u=0 deg, α_d=30 deg); however, increasing the stagger to Δy=1c resulted in two locked dominant frequencies. When the stagger was increased to Δy=2c, the secondary peak vanished and the primary dominant frequencies no longer matched. For all other configurations, no secondary peak was observed, but the upstream and downstream primary peaks matched.

Similar spectral behavior was observed from the spatially resolved power spectra of the vertical velocity component one chord downstream of each airfoil (Figs S3 and S4). In addition to the dominant frequencies, the velocity spectra revealed the behavior of the wake structure behind each model, indicating that the position and distribution of the wake varied with the configuration. For the biologically relevant cases (Fig. S3), the upstream high-power band was relatively narrow and symmetric relative to the airfoil position, with the maximum energy contained between y=−c and y=c. In contrast, the downstream wake was asymmetric and shifted upward, except for the (Δx=6c, Δy=2c, α_u=30 deg, α_d=20 deg) case, which had a symmetric spectrum. The extreme tandem effect cases (Fig. S4) showed a strong influence of the tandem configuration on the vertical position, distribution width and frequency of the high-power band. For the case (Δx=2c, Δy=0c, α_u=0 deg, α_d=30 deg), the power was spread over a large range of frequencies and no dominant frequency could be identified. By increasing the stagger to Δy=1c, a clear upstream band emerged at f*≈0.42, which was also observed in the downstream spectra at y=2c. In addition, the downstream spectra showed a distinct band at f*≈0.23 that was shifted downward. Increasing the stagger to Δy=2c caused the downstream band to become stronger and symmetric, with no trace of the upstream dominant frequency.

Fig. 6 shows the probability density of the upstream and downstream dominant frequencies for different combinations of the angles of attack. Each probability density curve was obtained from the compilation of all peaks in the spectra of any combination of gap and stagger. The curves show that the upstream model experienced a single peak, which was always close to f*≈0.23 when the upstream airfoil had an angle of attack of 30 deg; otherwise, it was a function of the configuration and varied with the angle of attack. The upstream dominant frequency was generally shared by the downstream model, which exhibited one or more additional peaks that could be smaller or larger than the upstream peak frequency. All in all, the tandem configuration influenced the spectra of the models by changing the dominant frequency, which was usually shared by both models, and by creating multiple peak frequencies for the downstream model.

The velocity fields obtained from the DPIV were used to visualize the flow topology and better understand the mechanisms through which the tandem configuration modified the aerodynamic performance of the system. Fig. S5 shows the mean velocity fields for the extreme tandem effect cases of Table 2. The background color indicates the velocity magnitude normalized with respect to the freestream velocity and the gray shaded regions indicate where the flow field data could not be captured because of the shadows in the laser sheet. With zero stagger, the downstream airfoil interfered with formation of the upstream wake and caused it to become excessively elongate, a feature that was not observed with Δy=1c or 2c. Furthermore, the mean flow speeds in the vicinity of both airfoils were altered by the change in the stagger. For instance, the flow speed close to the dorsal surface of the downstream model at Δy=0c was notably smaller than that at Δy=1c or 2c, whereas the flow speed close to the ventral surface of the upstream body decreased when the stagger was increased from Δy=0c to 1c.

For the configurations included in Table 2, the lift of the 0 deg angle-of-attack airfoils showed the strongest dependence on the tandem arrangement. The details of the flow structure around those airfoils were considered using the mean streamlines in their vicinity (Fig. 7). It was revealed that flow generally separated from the concave ventral surface, creating a ‘trapped’ vortex in the cavity. A similar behavior has been previously observed with a single snake-like airfoil at 0 deg angle of attack, for which the low pressure region created by the flow separation from the lower surface of the model dominated the dorsal wake and a negative lift was produced (Holden et al., 2014). However, with the tandem arrangement, the size of the dorsal wake and the strength of the trapped vortex varied with the configuration. In particular, for (Δx=2c, Δy=1c, α_u=0 deg, α_d=30 deg), the gap flow did not allow separation from the lower surface, and no ventral vortex could be formed (Fig. 7E).

Finally, the POD modes of the velocity field in each configuration were used to study the vortex structure in the flow. Santa Cruz et al. (2005) and Epps and Techet (2010) have shown that the first two modes can be used to accurately reconstruct the dynamics of von Kármán vortex streets. Furthermore, another index of the complexity and order of the spatio-temporal velocity field could be obtained by calculating ‘entropy’ from the energies of the POD modes (Santa Cruz et al., 2005). Entropy ranges from 0 to 1, with a small value indicating that most of the energy is contained in just a few modes (an ordered flow), and a large value indicating that the energy is distributed among several modes (a complex flow). Fig. S6 shows the entropies for the tandem configurations of Table 2. The majority of the tandem model entropies were larger than the maximum entropy of H=0.425 for a single model (Holden et al., 2014), meaning that the tandem model flow was less ordered and entailed greater turbulence. Also, the largest entropies belonged to the cases with the largest interaction between the two airfoils. In particular, a very large entropy for the case (Δx=2c, Δy=0c, α_u=0 deg, α_d=30 deg) was observed because the downstream airfoil interfered with the formation of upstream wake, which in return altered the flow around the downstream model. This is evident from Fig. 8A,B in which a strong and organized von Kármán vortex street could not be formed. With Δy=1c, the entropy decreased but the interaction between the models still existed as the vortices shed from the two models were coupled (Fig. 8C,D). With Δy=2c, the entropy decreased even more as the interaction between the two models considerably diminished, and the more energetic vortex shedding mode produced by the downstream model was only slightly influenced by the presence of the upstream model (Fig. 8E,F).

Our results showed that the tandem configuration of snake-like airfoils could produce both positive and negative aerodynamic interactions, depending on the arrangement. At best, an almost 10% enhancement of the lift-to-drag ratio can be achieved: a maximum of 2.2 for the tandem airfoils (Fig. 4) compared with a maximum of 2.0 for the single airfoil, excluding the special case of α=35 deg (Holden et al., 2014). In the other direction, the tandem lift-to-drag ratio could be 12% smaller than L/D_single. The results also indicated strong dependence of aerodynamic performance of the tandem snake-like airfoils on their relative position as well as their angles of attack (Fig. 3). In general, the results were more sensitive to stagger, and the tandem effects almost universally diminished when the stagger was larger than three cord lengths. This could be explained, for the most part, by the vortex–blade interaction (interaction of vortical unsteady flows with solid bodies; e.g. Rockwell, 1998), which was the main mechanism of aerodynamic force production in this system. With smaller stagger, the vortices shed from the upstream airfoil passed closer to the downstream airfoil, and therefore could produce a larger influence on it. This mechanism, however, could not explain the few cases for which the tandem effects increased beyond three chord lengths of stagger.

Flow field data helped to elucidate the details of how the lift and drag coefficients varied with the tandem configuration. Here, we focused on the cases with the extreme tandem effects (Fig. S5). The mean velocity fields showed that for Δy=0c, the interference of the downstream airfoil with the upstream wake resulted in a distorted wake. Moreover, the mean velocities close to the dorsal surface of the downstream airfoil did not reach the large magnitudes observed for Δy=1c and 2c. Therefore, the circulation-based downstream lift was considerably smaller. However, adding one chord of stagger modified the gap flow such that: (1) the upstream wake could form but was pushed by the gap flow closer to the model, (2) the mean velocity close to the ventral surface of the upstream airfoil decreased, and (3) the mean velocity close to the dorsal surface of the downstream airfoil increased. These mechanisms augmented both the upstream and downstream lift through vortex–blade interaction and circulation modification. Increasing the stagger to Δy=2c diminished the effects of the abovementioned mechanisms, and the combined lift and combined L/D decreased. Nonetheless, the downstream lift was a maximum at Δy=2c because the mean wake vortices were formed closest to the dorsal surface of the downstream airfoil (Fig. 7F; and also deduced from the smaller downstream wake in Fig. S5C).

The 0 deg angle-of-attack airfoils were particularly examined because their lift production was dominated by the separation of flow from the model's lower surface originating from the leading edge. Holden et al. (2014) observed negative lift for a single airfoil at 0 deg angle of attack and argued that it was the result of a low-pressure zone on the concave ventral surface created by a vortex trapped in the cavity. In this study, the negative lift was observed for almost all of the 0 deg angle-of-attack airfoils, but the magnitude of the negative lift was dependent on the tandem configuration. The most important cases were considered in Table 2 and Fig. 7. As shown by Fig. 7A, separation occurred at the leading edge and the resulting vortex created a strong downforce that resulted in negative lift. By increasing the stagger, the trapped vortex gained more strength (Fig. 7B,C) and created larger downforces, which explains why the magnitude of the negative downstream lift increased with stagger for these cases. Fig. 7D suggests the same negative lift mechanism; however, when the stagger was increased to Δy=1c, the gap flow changed in a way that did not allow the trapped vortex to form below the airfoil (Fig. 7E). As a result, no downforce was created and the upstream lift became positive. Further increasing the stagger to Δy=2c allowed a weak vortex to reappear (Fig. 7F), and the upstream lift became negative again, although with a smaller magnitude than the case with Δy=0c.

The vortex–blade interaction was also confirmed as the primary factor in the lift production of the tandem system by the spectral analysis results, which showed that the upstream and downstream dominant frequencies were similar for the configurations with significant effect on the lift coefficients. For instance, in Fig. 5A, the dominant frequency of the upstream airfoils (30 deg angle of attack) matched that of a single model with the same angle of attack, likely because the upstream flow was not influenced by the downstream model. The dominant frequency of the downstream airfoil (0 deg angle of attack) with Δy=1c and 2c matched that of the upstream airfoil (f*≈0.23), and a substantial lift enhancement was observed for these cases. However, at Δy=3c, the downstream dominant frequency became ‘unlocked’ from the upstream wake and its value suddenly changed to the single model value (f*≈0.34). Thereby, the lift enhancement was entirely lost (Figs 4 and 8). Another key feature of the tandem wake was that it had a larger spatial distribution than a single model. For instance, the spatial distribution of the wake of the model at 0 deg angle of attack in tandem configuration could be comparable to that of a single model at 60 deg angle of attack (compare Figs S3 and S4 with the results of Holden et al., 2014).

Although our results show that large lift enhancement could be obtained with some tandem configurations, it is not possible to conclude that the snakes take the staggered posture to exploit such aerodynamic enhancements. During a fully developed glide (i.e. at glide angles less than ∼40 deg), the body of flying snakes does pass through these configurations, indicating that they experience the tandem aerodynamic effects. But the snakes also take postures corresponding to other configurations with Δx>6c and Δy>3c (Socha et al., 2010), for which little or no tandem effects were found. Moreover, the snakes rarely exhibit the small gaps needed for the largest tandem effects we observed. However, the configurations with both models at angles of attack between 20 deg and 30 deg produced a near-maximum lift-to-drag ratio for almost any combination of gap and stagger (Fig. 4). So, it is possible in theory for the snakes to modulate their body posture to exploit such optimal angles of attack, where they can assume vastly different body postures without losing the maximal L/D. In fact, during the shallower part of the glide, the snakes maintain a body angle of 25 deg from the glide path (Socha et al., 2010), but no further inferences could be made owing to the lack of data for the exact variation of the airfoil angle along the body.

Although we observed a 10% enhancement of the overall lift-to-drag ratio in certain tandem configurations, this effect was not sufficient to explain the smallest glide angles observed in snakes. None of the configurations, even those with the largest aerodynamic enhancements, could increase the overall L/D to the observed values of 3.1–4.3 exhibited in the extreme glide performances (Socha et al., 2005, 2010). However, the results of this study are not sufficient to reject the hypothesis, as a large portion of the parameter space remains unexplored. In particular, we did not examine the special angle of attack of 35 deg, for which a single airfoil produced narrow dominant peaks in the lift coefficient (C_L=1.9) and the lift-to-drag ratio (L/D=2.7) with more than 30% increase with respect to the neighboring values (Holden et al., 2014). It might be possible that the tandem interactions augment the peak even more and produce sufficiently large lift-to-drag ratios. Additionally, this study did not address the effects of the sweep angle or the 3D shape of the snake body, both of which may have significant effects on the aerodynamics of the snake.

Generalizing our findings to flying snakes suggests that they can achieve relatively large changes in the overall L/D with slight adjustments in their posture. Therefore, we postulate that it is possible for the snakes to use the staggered configuration to control their glide trajectory. A simple mechanism for this control strategy was described in a previous theoretical study (Jafari et al., 2014), in which the trajectory of the flying snake was restricted to remain in the longitudinal plane. In this control hypothesis, the snakes would have to actively keep their body in a favorable posture in which the sensitivity of the aerodynamic forces to tandem configuration is high. The biologically relevant cases being within this region of high sensitivity conforms to this mechanism of control, but follow-up experimental work is needed to test the hypothesis. It is also possible that sensitivities have negative implications, and that small control errors lead to large rotational failures. However, the fact that in hundreds of observed glides no snake has lost control (J.J.S., personal observation) provides some anecdotal evidence that the system is not unstable. Furthermore, a recent study (Yeaton et al., 2020) that incorporated a more accurate representation of the snake's 3D motion, body shape, mass distribution and aerodynamics demonstrated that the aerial undulatory motion functions as a control mechanism for flight stability.

This study determined the aerodynamic performance of a tandem system with airfoils having the cross-sectional shape of the flying snake C. paradisi. As previously observed in studies with a single airfoil model, vortex shedding and flow separation dominate the aerodynamics of the tandem system. The main mechanism through which the tandem configuration modifies the overall aerodynamic performance is vortex–blade interaction. Overall, our findings provide new insight into the underlying physics of how the snake could maintain stability or maneuver in mid-air, despite the lack of specialized morphology serving as control surfaces. To test the existing hypotheses, these results should be incorporated into theoretical studies of the dynamics of flying snakes. Previous efforts (Jafari et al., 2014, 2017) have only considered the aerodynamic properties of a single, isolated airfoil, but clearly tandem interactions must be considered for a more accurate understanding of real snake gliding.

We thank Chris Williams and the DREAMS lab for fabrication of the snake model.