Influence of posture during gliding flight in the flying lizard Draco volans

The agamid lizards of the genus Draco are undoubtedly the most renown reptilian gliders, using their rib-supported patagial wings as lifting surfaces while airborne. Recent investigations into these reptiles highlighted the role of body posture during gliding, however, the aerodynamics of postural changes in Draco remain unclear. Here, we examine the aerodynamics and gliding performances of Draco volans using a numerical approach focusing on three postural changes: wing expansion, body camber, and limb positioning. To this aim, we conducted 70 three-dimensional steady-state computational fluid dynamics simulations of gliding flight and 240 two-dimensional glide trajectory calculations. Our results demonstrate that while airborne, D. volans generates a separated turbulent boundary layer over its wings characterized by a large recirculation cell that is kept attached to the wing surface by interaction with wing–tip vortices, increasing lift generation. This lift generating mechanism may be controlled by changing wing expansion and shape to modulate the generation of aerodynamic force. Furthermore, our trajectory simulations highlight the influence of body camber and orientation on glide range. This sheds light on how D. volans controls its gliding performance, and conforms to the observation that these animals plan their glide paths prior to take off. Lastly, D. volans is mostly neutral in pitch and highly maneuverable, similar to other vertebrate gliders. The numerical study presented here thus provides a better understanding of the lift generating mechanism and the influence of postural changes in flight in this emblematic animal and will facilitate the study of gliding flight in analogous gliding reptiles for which direct observations are unavailable.


Introduction
Life in the trees comes with a series of challenges for animals, including moving in a three-dimensional (3D) environment with inclined surfaces and substrates of variable width and compliance [1,2].Moreover, arboreal habitats are discontinuous, making it a necessity for animals to employ gap-crossing behaviors in order to move around between branches and trees [3,4].In this context, several vertebrate and invertebrate taxa are capable of aerial gliding locomotion, whereby they perform a controlled descent by converting gravitational potential energy to aerodynamic work, using their body as a wing [5][6][7].This results in net horizontal movement and allows them to move efficiently in their 3D habitat.
Among arboreal reptiles, the agamid flying lizards of the genus Draco, representing around 45 species endemic to South-East Asia [8], are undoubtedly the most renown gliders.These gracile squamates show a suite of extreme anatomical specializations to gliding locomotion, including a broad patagium spanning the entire length of the body as well as a pair of hyoid folds that serve as canards to increase stability in flight [9][10][11].The patagium of Draco is supported by five to seven ribs depending on the species [12,13], which can be extended laterally through muscular contraction and active manual handling to increase the area of the lift-generating surface [9,10,14,15].
The gliding prowess of Draco was recognized early on, but for a long time, our understanding of the aerodynamics of gliding flight in these reptiles remained obscured by their cryptic lifestyle and limitations of experimental setups [9,11,16,17].New information has come to light recently from direct observations [15,18,19], wind-tunnel experiments [20,21], and computational fluid dynamics (CFD) [21].These highlight that Draco: (i) usually keeps its wings open using its forelimbs, forming a 'composite wing' [15] that increases lift generation [21]; (ii) negotiates obstacles and plans glide paths using visual cues [18]; (iii) actively changes its body posture, and in turn wing shape, as well as its tail orientation while airborne, allowing for efficient maneuvering [18,20] and control over aerodynamic forces and moments [19]; (iv) takes advantage of wing-tip vortices to keep the flow attached over its wings and delay stall [21].
However, despite these recent investigations, much remains unclear about the effects of postural changes on the aerodynamics of gliding flight in Draco.In this study, we thus address the aerodynamics of Draco volans through CFD numerical simulations, with a focus on three postural changes: wing expansion, body camber, and limb positioning.Such postural changes are indeed commonly observed in Draco during gliding, allowing this animal to control its aerodynamic forces and moments, and negotiate its glide path [15,18,19].
The calculated aerodynamic coefficients are then fed into a two-degrees-of-freedom flight dynamics model and used to simulate glide trajectories under various conditions.This work thus aims to: (1) describe the lift generation mechanism employed by D. volans for gliding; (2) determine the aerodynamics of the highly specialized morphology of D. volans; (3) identify the impact of postural changes on gliding performance.The new data presented here aim to complement published experimental results and further our understanding of the gliding behavior of this emblematic animal.

Specimen data
The 3D models (hereafter 'geometries') used in this study derive from the Draco volans specimen Louisiana State Museum of Natural History (LSUMZ) herp 81 750 housed in the LSUMZ, Baton-Rouge, Louisiana, USA.This specimen was CTscanned in the University of Florida Nanoscale Research facility (70 kV, 200 µA, voxel size 0.063 mm) and is accessible on the MorphoSource 2.0 repository (http://n2t.net/ark:/87602/m4/M77362).All vertical slices were converted from 16 bits to 8 bits using ImageJ v.2.1.0to save storage space.The dataset thus comprises 1862 vertical slices with a resolution of 351 × 800 pixels and a voxel size of 0.063 mm.The slices were then imported into the 3D segmentation software Mimics v.21.0 (Materialise, Leuven, Belgium), and 3D surface models of the skin and skeleton of the specimen were generated in individual masks using the threshold tool for given gray levels and exported as STL files.

3D geometries
The segmented geometry obtained is highly irregular and comprises numerous intersecting faces.We thus modified it through various processing steps in Blender v.2.90.1 in order to obtain a corrected CTscan-derived geometry.In addition, we created an armature that mimics the skeleton of Draco so that individual body segments of the geometry could be moved in the regions of the animal's joints.This allowed us to generate different geometries of Draco under different postures.See supplementary information for more details on geometry reconstruction (figure S1).
We thus constructed geometries of D. volans that mimic gliding postures known for this animal (inspired by the sequential short-exposure photographs of [15]).First, these postures represent increasing degrees of body camber, defined as the angulation between the anterior and posterior trunk portions of the geometry.We selected three degrees of camber, P0, P7, and P15, where the anterior and posterior trunk portions are set at ca. 0 • , 7 • , and 15 • to each other (figure 1).For P15, the tail is set at ca. 60 • to the horizontal as observed in living Draco individuals during gliding flight [15].Second, these postures differ in the orientation of the forelimbs, which can be oriented anteriorly (hereafter referred to as 'armsfirst' (AF) postures), or oriented transversely over the leading edge of the wings (referred-to as 'armstransverse' (AT) postures), mimicking the composite wing employed by Draco [15].Lastly, one of the postures has partially folded wings ('wings-folded' (WF) posture) as preserved on specimen LSUMZ herp 81750.
In summary, we constructed seven geometries of D. volans: the 'resting posture' P0 WF which is based on specimen LSUMZ herp 81750 (figure S1(e)), postures P0 AF (or 'standard gliding' posture), P7 AF , and P15 AF representing AF postures with increasing degrees of body camber, and their counterpart AT postures P0 AT , P7 AT , and P15 AT .All of these rigid-body posture geometries were exported in STL format.
We measured the planform area A (the area of the wings and trunk, 'total wing area' of [16]), wing area A wing , and total area A tot by isolating the upper surface cells of each geometry (i.e.removing the lower surface cells to exclude volume information) and projecting them on the horizontal plane in Paraview v.5.9.0.P0 WF and P0 AF have an A of 18.41 cm 2 and 23.87 cm 2 respectively, both of which fall within the range known for D. volans [10,22].We use A as the reference area for all our geometries following common practice for airplane-like biological models [23].

Measurements
The mass m of a geometry was estimated under P0 WF and P0 AF postures by measuring the volume of the geometry prior to smoothing (figure S1(d)) and multiplying it by a density value of 1000 kg.m −3 , similar to the volumetric mass estimation methodology described by Allen et al [24].We did not attempt to reconstruct the air-filled volume of the mouth, trachea, and lungs as they could not be reliable identified on the CT slices.
The estimated mass for P0 WF and P0 AF was 7.25 g and 8.29 g respectively.Both estimations fall within the range of mass known for D. volans (between 2 g and 13 g [10]; mean value of 8.475 g in [25]; see also [22]).The difference in mass between both postures results from the fact that our geometries do not maintain a constant volume when deformed (see supplementary information).Thus, the volume and mass of P0 AF are likely overestimated due to the extension of the wing.As P0 is a direct reproduction of the original specimen LSUMZ herp 81750, its mass of 7.25 g will be used for all subsequent analyses.
The wing loading (WL) of P0 WF and P0 AF is 38.63 N m −2 and 29.80 N m −2 respectively.Both values are within the range of WL known for D. volans [22] and are comparable to those reported for the similar-sized Draco maculatus as well (Draco whiteheadi in [9]).Conversely, those WL values are much higher than those reported by McGuire and Dudley [17] for Draco species of similar mass (Draco quinquefasciatus: m = 6.5 g and WL = 10.5 N m −2 ; Draco formosus: m = 8.8 g and WL = 14.1 N m −2 ), and even surpass that of the largest measured species (Draco fimbriatus: m = 18.7 g and WL = 23.5 N m −2 ).As both our mass and planform area measurements fall within the ranges for D. volans published in other studies, we are unsure why our WL differ from those reported for McGuire and Dudley [17] (they might result from different reference area calculations), and thus refrain from comparing absolute values of WL here.
Table 1.Key measurements of the resting (P0WF) and standard gliding (P0AF) postures.mc is the mean value of rc, tc, and ten regularly spaced chords along the width of the wing; R is the aerodynamic center, taken at 25% of mc; AR equals s 2 A −1 ; WL equals mg A −1 using m of P0WF for both geometries, with g = 9.81 m s −2

Numerical setup
Animal gliding has traditionally been described through steady-state aerodynamics under the assumption that animals glide at equilibrium, whereby the aerodynamic forces counteract weight so that the animal glides in a straight path at a fixed glide angle and velocity [6,23,26].However, there is now a consensus that animals actually rarely reach equilibrium gliding, if at all, and are instead actively maneuvering during gliding [18,19,26,27].Nevertheless, equilibrium gliding is theoretically possible at least in some taxa [27], especially for medium-to long-range glides (height loss > ca. 25 m [28]), and appears relatively more common in Draco lizards than in any other gliders [6,17], provided there are no obstacles in their path.Such glide paths may not accurately represent the complex real glide paths of Draco lizards in forest environments [18], but indicate that steady-state aerodynamics are a relevant first step to the study of gliding in Draco.Similarly, although Draco changes is posture and wing shape while gliding, the study of the aerodynamics of static postures has proven to be a relevant first step, as shown by Lau et al [21] regarding forelimb position.
An unobstructed glide trial in which the animal does not change its posture is essentially identical to a wind tunnel or CFD setup where a Draco geometry is immersed in an incoming uniform airflow, which is what is reconstructed here.Moreover, as equilibrium gliding implies a constant glide angle, it can be described using static simulations rather than dynamic ones which incorporate motion but are very computationally demanding.
Given that the total length of Draco ranges from 15 to 40 cm (60-150 mm in snout-vent length (SVL) [16], and considering an SVL:total length ratio of ca.2.5 as in table 1), and that these animals glide at mean velocities ranging between 3 and 10 m s −1 [17], we estimate that these animals glide in a range of Reynolds numbers (Re) of 10 4 -10 5 , similar to previous estimations [19,21].This range is typical of vertebrate gliders and flyers and is indicative of transitional to turbulent flows [29][30][31].Given the nature of this flow, we thus conduct three-dimensional steadystate Reynolds averaged Navier-Stokes (RANS) CFD simulations of gliding flight.RANS simulations were selected over Detached Eddy Simulations (DES), as done by Lau et al [21] in order to reduce computational cost and allow for the conduction of the large number of simulation cases required for this study.
All CFD cases were run using the open-source software OpenFOAM v.2012 on a Gnu/Linux Debian 9 workstation (Intel(R) Xeon(R) CPU X5650, 2.67 GHz, 12 cores, 24 Go Ram).All cases were run in parallel on six processors.
We set up a rectangular (2.6 × 2 × 2 m) computational domain extending three times the length of our D. volans geometry upstream, ten times downstream, and five times in all other directions (figure 2(a)).We then discretized the computational domain into a series of hexahedral and split-hexahedral cells using the snappyHexMesh mesh generation utility of OpenFOAM.We specified a refinement box (where the mesh is iteratively refined) extending four times the length of the D. volans geometry downstream and half this length in all other directions for all cases (figure 2(b)) except for those derived from postures P15 AF and P15 AT where the refinement box extends the length of the geometry downward to encompass the strongly downturned tail.The mesh was clustered near the body to resolve the turbulent boundary layers (figure 2(c)).Meshing parameters, and sensitivity analyses regarding mesh resolution are presented in the supplementary information.
The flow field around D. volans is simulated by solving the RANS equations (which express conservation of time-averaged mass and momentum) in conjunction with the k − ω shear-stress transport turbulence model, a widely used and satisfactorily accurate model in order to predict fluid flow [32,33].The solution of the discretized governing equations is computed by using the Semi-Implicit Method for Pressure Linked Equations (SIMPLE) algorithm, commonly used for steady-state RANS simulations [34].The convective terms in the equations are approximated using a bounded second-order accurate upwind finite-volume scheme for the convective terms and a second-order Gauss scheme for the viscous terms.The iterative Gauss-Sidel model is used to solve the pressure and velocity equations, as well as the transport equations for the turbulence model.
The upstream and top faces of the rectangular domain relative to the D. volans geometry were designated as the inlet (i.e. the inflow) with a turbulence intensity of 5% and a turbulence length scale of 7% of the characteristic length of the geometry (here taken at the total length l, table 1).At those boundaries, the inlet velocity (equal to the flight velocity) is assigned along with the inflow angle.The latter can be modified to simulate different angles of attack (AoA).The downstream and bottom faces were designated as the outlet (i.e.outflow), where a zero-pressure outlet was imposed.Lastly, a no-slip boundary condition was applied to the D. volans geometry.
Each case was assigned a specific velocity and AoA, defined as the angle between the main axis of Draco (i.e. the line between the tip of the snout and the midpoint of the body), and its velocity direction (identical to that of the simulated airflow, figure 3(a)).Draco has been observed to glide at relatively high AoAs [19], and previous computational studies have simulated relatively lower AoAs to determine the stall angle in Draco, which corresponds to an AoA of ca.20 • [21].Nevertheless, previous studies have mostly focused on tracking AoA in the mid-glide and landing phases in Draco, but not take off [19,21].Yet, it has been shown that Draco can launch at a wide range of negative take off angles [18], suggesting a similarly wide range of AoAs.In addition, it has been suggested that Draco may use its gliding ability to escape predators using much steeper take offs to increase its initial velocity (e.g. to escape the gliding colubrid snake Chrysopelea [35]).Therefore, we selected a set of 6 angles of AoAs:   .This range does not extend to the high AoAs sampled by Lau et al [21] as these cases give rise to highly unsteady separated flow conditions, and did not converge using our steady RANS numerical setup.However, it does cover the stall angle of 20 • recovered by this study using DES simulations.Lastly, this range covers some low negative AoAs to account for steeper glides, similar to Zhao et al's [36] simulated range for the flying squirrel Glaucomys.
In total, 70 3D steady-state CFD simulation cases were run at a uniform freestream velocity U∞ of 5 m s −1 , corresponding to seven different postures (P0 WF , P0 AF , P0 AT , P7 AF , P7 AT , P15 AF , P15 AT ) and to ten AoAs, taken in the set (−10 and 25 • ).Each case is run for a number of iterations (ranging from 140 to 900) tailored to allow convergence of the results to a steady-state solution.Sixty-two additional cases were set up for sensitivity analyses (see Supplementary Information).
For a freestream velocity of 5 m s −1 , and using the mean patagial chord mc as a reference as done by Lau et al [21], we estimate a Reynolds number (Re) of 1.07 × 10 4 for P0 AF .This is the same order of magnitude as the values calculated by Lau et al [21] for Draco dussumieri, a slightly larger species (based on data from Khandelwal and Hedrick et al [19], mean chord 56 mm), and slightly lower than that obtained from their CFD analyses, but this is likely due to their much larger geometry (mean chord 79 mm) and higher freestream velocity of 7 m s −1 .

Model validation
This work relies on a numerical approach to assess the gliding performances of D. volans, namely the lift, drag and pitching moment coefficients, defined as: We confronted our calculated aerodynamic coefficients to those of Lau et al [21], showing a very good correspondence in our C D results, and with some degree of caution for higher AoA, with our C L results as well (detailed in supplementary information; figure S2).Thus, we consider that our model can be reliably applied to the study of gliding locomotion in D. volans.

Trajectory simulations
The calculated values of the aerodynamic coefficients for each geometry and AoA were used to simulate 240 two-dimensional trajectories using different jump heights, postures and angles of attack by resolving equations of motion through time, following similar studies [37].The following equations are solved for each time step: where V is the glide velocity magnitude (m s −1 ), V x and Vy are the horizontal and vertical glide velocities (m s −1 ) respectively, γ is the glide angle, F x and F y are the horizontal and vertical aerodynamic forces (N), a x and a y are the horizontal and vertical accelerations (m s −2 ), and g is the gravitational acceleration (taken at 9.81 m s −2 ).
For each time step ∆t (set equal to 1 ms), the velocity and location for the following time step are determined by: where t is the time (s) and X and Y are the horizontal and vertical locations (in meters) of the center of mass, respectively.
We compute four sets of simulations, starting from different heights: 10 m, 15 m, 20 m, and 30 m.Although perch height does not directly affect the trajectories, higher perches offer more vertical space to develop glide phases.In turn, simulating different perch height allows us to assess the gliding behaviors and performances of Draco, and compare the efficiency of different postures in various settings (e.g.short range versus long range jumps).
First, we selected 10 m to encompass the recorded glides of Khandelwal and Hendrick [19] (max observed take-off height 4-5 m) as well as that of Khandelwal and Hendrick [18] (mean take-off height 7.4 m), and the majority of the experimental glide trials conducted for the gliding snake Chrysopelea which inhabits the same habitats as Draco [35].
We then conducted simulations from 15 m perches (as also done for Chrysopelea [38]), as well as 20 m and 30 m perches to encompass a broad array of perch heights.30 m represents an estimate of tropical rainforest canopy height worldwide (although the Indo-Malayan rainforest inhabited by Draco are typically tens of meters higher [39]), and remains within the range of perch heights observed in in Draco [12].This height is also characteristic of medium-to longrange glides in vertebrate gliders, for which optimal glide efficiency can be reliably approached by maximizing glide range using steady-state aerodynamics [28].
All glides started with an initial horizontal velocity V x (0) of 2.5 m s −1 and no vertical initial velocity (i.e. a strictly horizontal take off), as used by Clark et al [20] for similar trajectory simulations.This is close to the velocity magnitude reported by Khandelwal et al [40] by the time Draco fully extends its wings.The posture and AoA are also fixed throughout the glide to enable the use of the aerodynamic coefficients obtained from the CFD analyses.While this is a simplification of the actual gliding behavior of Draco, which actively controls its glide trajectory [15,18,19], it allows for a systematic assessment of posture on gliding performance.
The trajectory simulations were run for one of two AoAs: that which maximizes the lift-to-drag ratio C L /C D (case-dependent), or that which maximized C L .We first selected the AoA of (C L /C D ) max , which corresponds to a glide angle γ which maximizes glide range and is a standard measure of gliding performance [23,41].However, animals rarely maximize glide range during gliding, and adopt different AoAs depending on the situation [18,28].For shorter glide paths (e.g. between neighboring trees), the animal might maximize lift generation instead to ensure a shorter ballistic dive.Thus, we selected a second AoA of 25 • , corresponding to the maximum computed C L for all postures (see below), and conforms to the AoAs maximizing C L in other studies on Draco [20,21].

Influence of posture on the aerodynamic coefficients
For all postures, the lift coefficient C L increases with AoA (figure 3(b)).There is no trace of a plateau or decrease in C L , indicating that D. volans does not experience stalling for AoAs below 25 • , the highest AoA considered here.Conversely, Lau et al [21] reported that stalling occurred respectively at ca. 20 • or 25 • for CFD and wind tunnel experiments in Draco, which conforms to the value obtained by Clark et al [20].The higher stall angle suggested here likely accounts for our more detailed reconstruction of the wings, or for intraspecific variation, however comparisons at high AoA must be interpreted with caution due to the uncertainties associated with the turbulence model.P0 WF shows very similar C L values compared to P0 AF and P0 AT .The P0 AF , P7 AF , and P15 AF cases show that C L markedly increases with camber (figure 3(b)).This is particularly striking for negative AoAs (P0 AF : −0.43 and P15 AF : −0.12 at −10 • ) and for high positive AoAs (P0 AF : 0.81 and P15 AF : 1.05 at 25 • ).However, there is hardly any difference between counterpart AF and AT geometries (figure 3(b)).In contrast, Lau et al [21] reported a slightly higher C L for AT postures than AF postures in Draco.However, this may result from our detailed reconstruction of the hand, especially the palm and digits that may function as a secondary airfoil in an AF posture.
The drag coefficient C D increases when the AoA departs from 0 • for positive AoA, and to a lesser degree for negative AoAs (figure 3(b)).P0 WF shows very similar C D values compared to P0 AF and P0 AT .The P0 AF , P7 AF , and P15 AF cases show that C D strongly increases with camber for high positive AoAs (P0 AF : 0.42 and P15 AF : 0.63 at 25 • ).This difference decreases with shallower AoAs, and there is only a small difference in C D for high negative AoAs (P0 AF : 0.15, P7 AF : 0.15, and P15 AF : 0.21 at −10 • ; figure 3(b)).For all geometries, the counterpart AF and AT geometries have similar C D , but with slightly higher values for the latter for low positive AoAs (5 • to 15 • ).
The pitching moment coefficient C M decreases overall with AoA for all geometries despite a slight

Influence of angle of attack on the flow fields
Considering an AT uncambered gliding posture P0 AT , the behavior of the velocity (and pressure) field varies markedly with AoA (figure 4).For an AoA of 0 • , as seen on the parasagittal (streamwise) plane at mid-span, the airflow separates in the anterior half of the wing (as shown by the isoline of null horizontal velocity), forming a separation bubble immediately downstream of the leading edge, before reattaching in the posterior half of the wing.This separation bubble comprises a region of lower pressure and velocity, which is transversely restricted over the posterior wing and in the wake of the animal, covering only the middle portion of each wing's span (figure 4).
The frontal surface of the animal (i.e. that facing the airflow, figure 5(a)) is gradually displaced dorsally or ventrally for steeper AoAs (figure 4).As a result, the velocity and pressure difference around the wing becomes progressively stronger, such that starting from an AoA of 5 • , the airflow above the wing does not reattach, forming a continuous shear layer (figure 4).This indicates that the airflow is fully separated along the dorsal surface of the wings of the animal.This shear layer also becomes thicker as AoA increases (figure 4).In addition, as AoA increases, the legs become more exposed to the incoming airflow and generate a secondary region of low velocity in their wake (figure 4).
As shown by the velocity streamlines and vortical structures (visualized through an iso-surface of the ) with ∥ Ω ∥ the vorticity magnitude and ∥ S ∥ the strain rate magnitude), both reported in figure 5, the above-wing shear layer occurring for positive AoAs is characterized by the development of a large cell of turbulent flow recirculation.This cell starts at the leading edge and spans the entire length of the wings, and recirculates the air in a mostly anteroposterior (chordwise) loop.It generates dorsally a leading-edge suction that increases in strength for positive AoAs (figures 4 and 5).A smaller, slightly posteromedially oriented wing tip vortex also occurs.This vortex swirls medially (from the wing tip to the dorsal wing surface) just lateral to the recirculation cell.As described by Lau et al [21] using time-resolved simulations, the wing-tip vortices likely help in keeping the flow energized and attached over the wing (figure 5).

Influence of posture on the flow fields
We now compare the aerodynamic performance for various body postures.The size of the low-velocity regions is much smaller in P0 WF than in the other geometries due to this posture's partially folded wings (figure 6).Similarly, these regions strongly increase in size from P0 AF to P15 AF and P0 AT to P15 AT (figure 7).At the same time, the leading-edge suction increases in strength for positive AoA.As a result, the size of the recirculation cell increases markedly with camber, forming more regular vortices above the wings (figures 8 and S3).
The flow fields behave similarly between counterpart AF and AT geometries with respect to AoA.However, the diameter of the wing recirculation cell for positive AoAs is always greater for AT geometries than for their AF counterparts (figure 8).At the same time, the above-wing vortices are much larger and more regular in AT geometries (figure S3).In agreement with Lau et al [21], we suggest that the vortices and recirculation cells may be smaller and more irregular in AF postures because the forelimbs perturb the fluid flow upstream of the leading-edge (figure S3), promoting transition to a turbulent state shortly upstream of the leading edge of the wing (figure 7).However, determining the exact location of the transition point requires computationally intensive high fidelity simulations, which are beyond the scope of this work.In addition, for AT geometries, the leadingedge of the wing is more blunt, increasing the leadingedge suction and the size of the separation layer at this level, and thus the diameter of the recirculation cell.
Excluding the wings, the body of the animal also has influence on the scalar fields and vortices.For all geometries, a narrow vortex forms at the back of the head and runs along the back of the animal, increasing in diameter with camber (figure S3).A separation layer is also visible in the sagittal plane of the animal, occurring as a separation bubble at the base of the tail for P0 AT , and more anteriorly as a fully separated layer for more cambered geometries (figure 9).In addition, a boundary layer develops over the tail, especially for highly cambered geometries (P15 AF , P15 AT ) where it is angled downwards.For high positive AoAs (e.g.20 • ), a separation bubble occurs at the midlength of the tail in P15 AF and P15 AT , which does not occur for less cambered geometries (figure 9), indicating that the tail functions as a secondary airfoil independent  from the wings and the orientation of the body.Lastly, as mentioned above, the legs of Draco also function as secondary airfoils, generating secondary separation and vortices (figures 5-8 and S3).This is further exacerbated for more cambered geometries (P15 AT ), where the legs contribute more to the frontal surface and displace a further amount of air while gliding (figure 8 and S3).

Influence of posture on glide phases
All gliding trajectories from a 30 m perch can be divided in three or four informal phases (figure 10): (1) a ballistic dive where the animal falls rapidly while the glide angle steeply increases, and where the acceleration is essentially downwards.This phase ends once the glide angle reaches its maximum, in agreement with similar definitions for real glide trials [6].(2) A shallowing glide where the animal follows a parabolic trajectory while the glide angle decreases.The total velocity first increases, reaching its global maximum, while the acceleration vector rapidly orients upwards, then both the total velocity and acceleration decrease in intensity, with the latter reorienting backward.This phase ends when the glide angle and the vertical velocity both reach their local minimum, which coincides with the point where the acceleration vector reorients downwards (figure 10).(3) A descent phase where the trajectory is markedly straighter.
Both the glide angle and velocity increase to a local maximum, which coincides with the upward reorientation of the acceleration vector.The acceleration then gradually decreases in intensity while the glide angle and velocity tend to stabilize.(4) An equilibrium glide, where the acceleration is essentially null (i.e.<5% of maximum acceleration for each posture) such that the animal glides in a straight path at a constant glide angle and velocity, as defined by aerodynamic theory [6,26].This latter phase occurs only for postures P7 AF , P15 AF , and P15 AT for a 30 m perch (figure 11(a)).For all other postures, the acceleration decreases sharply but retains a significant intensity until the animal reaches the ground.Note that the glide phases described here characterize simulated glides of a rigid-body geometry falling to the ground, and thus do not accurately reflect the glide phases previously described for real glide trials in vertebrates [6] and Draco in particular [18].
Given the higher (C L /C D ) max of P7 AF , this geometry is expected to cover the greatest glide range.P7 AF indeed produces the longest glide from a 30 m perch (86.07 m), with P15 AF falling slightly shorter (74.99 m), and P0 AF even more so (68.61 m) (table 2; figure 10(a)).P7 AF also produces the longest glide range for 15 m or 20 m perches, but P15 AF and P15 AT cover a greater glide range for a 10 m perch (table 2, figure 12).Moreover, AT geometries cover shorter or similar glide distances than their AF counterparts for all jump-heights and postures, which conforms to their lower (C L /C D ) max (table 2; figures 10-12).This is particularly striking for a 30 m perch, where P7 AT covers a much shorter distance (66.46 m) than P7 AF (86.07 m).
Considering a 30 m perch, P7 AF , P15 AF and P15 AT have the shortest ballistic dive phases of all geometries, in line with their lower maximum glide angles (table 2; figures 10 and 11(a)).As a result, the animal accumulates less downwards velocity before the start of the shallowing glide phase, and the maximum velocity achieved by these postures is thus much lower than that of other geometries (P0 AF : 13.37 m s −1 and P15 AF : 9.43 m s −1 ; table 2; figure 10(c)).In turn, the duration of the shallowing glide phase is much shorter such that the descent phase begins earlier (P0 AF : 3.92 s and P15 AF : 2.94 s) and after a much shorter fall (P0 AF : 19.61 m and P15 AF : 10.89 m).Lastly, the gliding trajectories of these geometries are the only ones to reach an equilibrium glide phase (figure 11(a)).The animals thus fall at a constant velocity that is much lower than that of other geometries (figure 10(c)) and are able to stay airborne for a longer time and to cover a much longer glide range (table 2; figure 10(a)).; smaller markers indicate all other α.Note that the differences between the trajectories reflect only the vertical space available, as the computed trajectories are all identical relative to perch height but stop once the animal reaches the ground.
The duration and range of all glides increases markedly with perch height (figure 12).However, each glide trajectory being identical in our rigid-body simulations, these differences only reflect the capacity of each phase to develop in the time and space available before the animal reaches the ground.For 20 m perches, only postures P15 AF and P15 AT reach equilibrium gliding, whereas postures P0 AT and P7 AT do not reach the descent phase (table 2; figure 12).Only postures P7 AF , P15 AF , and P15 AT reach the descent phase for 15 m perches, none of them reaching equilibrium.Lastly, none of the geometries reach the  descent or equilibrium phases for 10 m jumps, which are solely described by the ballistic dive and shallowing glide (table 2).

Influence of angle of attack on glide trajectories
The trajectory simulations differ markedly according to the AoA at which the geometry is set.The ballistic phase is extremely steep for negative AoAs such that the animal falls roughly vertically (figure 13).For higher AoAs, the slope of the ballistic phase shallows as AoA increases such that the animal falls less and transitions more quickly to the shallowing glide phase.This is particularly striking for an AoA of 25 • , which has a much shallower slope than the AoA of (C L /C D ) max for all postures (figures 12 and 13).As a result, the shallowing glide phase is very short such that the descent or equilibrium glide phases begin after only a 10-12 m fall for this AoA (figure 11(b)).The descent phase is also extremely short compared to (C L /C D ) max cases, and this phase is even absent for postures P15 AF and P15 AT , where the equilibrium phase immediately follows the shallowing glide (figure 11(b)).Overall, the acceleration is much smaller than for (C L /C D ) max cases, which allows all postures to reach equilibrium gliding.As a result, the cases set at 25 • provide a much shallower glide, allowing the animal to cover some the longest glide ranges of all AoA considered here when jumping from a 10 m perch (figure 13).However, for higher perch heights, an AoA of 25 • provides shorter glide ranges relative to other AoAs, including that of (C L /C D ) max (figure 12).For a 30 m perch, the (C L /C D ) max cases recover the longest glide ranges of all AoAs for all postures.Lastly, the minimum perch height for which the (C L /C D ) max cases recover the longest glide ranges appears linked to camber.(C L /C D ) max cases indeed recover the longest glide ranges starting from 15 m perches for P7 AF , P15 AF and P15 AT whereas this occurs only for 30 m perches for P0 AF and P0 AT (figure 12).

Discussion
Our CFD results demonstrate that while gliding, D. volans generates large recirculation cells over its wings that induce regions of low pressure and velocity over most of the wing and body of the animal and in its wake (figures [4][5][6][7][8][9].This generates a strong difference in the pressure and velocity distribution above and below the wing (figures 4, 7), which enhances lift generation (figure 3), as predicted by aerodynamic theory [29,31].Our RANS simulations indicate that the airflow separates above the wing in close proximity to the wing tip vortices (figures 5, S3).This likely indicates that the separated region interacts with the wing tip vortices which help in keeping the separated region attached to the wing surface, as described by Lau et al [21] using time-resolved DES simulations in D. volans.In turn, this maintains the lift mechanism at high angles of attack.This interaction is typical in low aspect-ratio wings [42,43], such as those of vertebrate gliders [6,21,44].In addition, postural changes result in changes in the pressure and velocity distribution, with the difference between the upper and lower wing surfaces increasing with camber (figures 4-7), as well as in changes in the diameter of the recirculation cells (figures 8, S3).In turn, these differences in velocity and pressure distributions result in an increase in the generation of lift, and to a lower degree of drag (figure 3(b)) as is typical of low aspectratio wings [45].Thus, active control of body posture in flight is expected to provide D. volans with a measure of control over its lift enhancing mechanism.This conforms with field observations that demonstrate that Draco actively controls its body posture during gliding, namely its body camber and orientation [18], forelimb position [15], tail orientation [20], and wing camber [19].
The leading edge of the wing always faces the incoming airflow in D. volans and is thus subject to a pressure that would push it backward and upward for positive angles of attack (figure 5).While it remains unknown whether the strength of thoracic musculature of D. volans is sufficient to compensate this increase in air pressure while gliding, it is evident that holding the leading edge of the wings with the forelimbs would provide further support to keep the wing expanded.This may explain why D. volans typically adopts an AT posture when gliding [15].From an aerodynamics standpoint, the expansion of the wing increases the planform area, which in turn increases the aspect ratio and decreases wing loading (table 1).As glide range decreases when wing loading increases in D. volans [11,16,17], keeping the wing fully expanded appears paramount to increase gliding efficiency.Furthermore, D. volans may benefit from modulating its wing area during gliding, for instance to maneuver around obstacles [7,19].The differences observed here between P0 WF and P0 AF indeed indicate that such modulation would reduce the size of the vortices and low pressure and velocity regions around the animal (figure 6), providing a means of reducing the instantaneous aerodynamic force.
Our trajectory simulations show that when the lift-to-drag ratio is maximized, that is, if D. volans is capable of maintaining its body orientation close to the corresponding angle of attack for a given posture, then the animal is theoretically capable of reaching equilibrium gliding (figure 11(a)), corresponding to its maximum glide range for a 30 m jump (figure 12).Glide range can also be further extended by adopting a more cambered posture, as the increase in pressure and velocity difference and, to a lesser degree, in the lift-to-drag ratio with camber allows these postures to stay airborne for a longer time than uncambered ones regardless of perch height (table 2; figure 12).
However, when considering a fixed angle of attack, maximizing glide range relies heavily on perch height and path obstruction, as all glide phases need sufficient time and vertical space to develop (figures [11][12][13].While this can be attained by maximizing the lift-to-drag ratio when jumping from high perches (e.g. 30 m), only highly cambered postures maintain this trend for medium-high perches and none of our simulations recover this strategy as maximizing glide range for short perches (e.g. 10 m) at which points other angles of attack provide longer glide ranges for all postures considered (table 2; figure 13).
Nevertheless, our CFD analyses also demonstrate that active body positioning would allow the animal to control its aerodynamic performance and the distribution of its glide phases.This is further supported by our trajectory simulations as the duration and falling height of the ballistic dive and shallowing glide phases decrease with both camber and angle of attack (figures 11 and 13).Based on the lift generating mechanism described above, we suggest that by increasing its camber and angle of attack, D. volans could increase the size of the recirculation cells and attached shear layers above its wings, as well as secondary layers developed around its legs and tail (figures 4, 7 and 8), which all contribute to lift and drag generation (figure 3(b)).In turn, this would decrease both the glide angle and the total velocity of the glide, which underpin the transition from the initial fall (i.e. the ballistic dive and shallowing glide phases) to the descent phase (figure 11).
In line with this, our trajectory simulations at various angles of attack show that for lower perch heights (e.g. 10 m), increasing the angle of attack significantly increases glide range while reducing height loss (figures 12 and 13).In addition, by increasing body camber and/or the angle of attack, and in turn lift and drag generation (figure 3), the glide angle, velocity and acceleration all decrease significantly overall, including at the end of the trajectory (figures 10 and 11).In turn, this reduces the duration of the descent phase in favor of the equilibrium glide phase, which occurs when the acceleration is very low (figure 11(b)).
These results conform to the observation that vertebrate gliders choose their glide paths prior to takeoff [18], and that their angle of attack during gliding does not necessarily maximize their lift-to drag ratio for a given situation [28].For example, when considering only the initial portion of their glide (i.e. the ballistic dive and beginning of the shallowing glide; figure 13), increasing the angle of attack (but keeping it below the stall angle) enables the animal to fall more slowly and for a shorter height, following a shallower trajectory than when the lift-to-drag ratio is maximized (figure 13).This would appear beneficial for an animal aiming for a neighboring branch rather than one further away [6,28], and conforms to the observation that Draco takes off at a shallower angle when planning a shorter glide [18].Conversely, an animal may choose to maximize its gliding speed at the cost of glide range by adopting a lower angle of attack (figures 12 and 13).Such behavior is indeed observed in vertebrate gliders when escaping from a predator [6].As suggested by Socha et al [35] this may be especially true for Draco as one of its predators, the flying colubrid Chrysopelea, is capable of gliding at higher speed than Draco.
However, contrary to the rigid-body gliding simulations presented here, Draco is known to change its posture and angle of attack while gliding [19].As a result, glide range appears independent from perch height in Draco (at least for relatively short perches) [18].Nevertheless, our systematic approach to the study of the aerodynamic and gliding performances of each posture in D. volans provides some insight into the influence of postural changes in flying lizards.An individual may for instance begin gliding by adopting an uncambered posture at a low or negative angle of attack to maximize its ballistic dive phase and vertical acceleration before increasing its camber and angle of attack to decrease its shallowing rate.Such an active gliding behavior has already been demonstrated by direct observations [18,19], and wind tunnel experiments [20], but is also supported by the very low pitching moment coefficients recovered here (figure 3(c)).The range of C M obtained here indeed suggests that D. volans is mostly neutral in pitch for the range of angles of attack considered, although cambered postures are slightly more stable.However, additional pitch stability may be provided by the hyoid folds, which are not considered here.In other words, D. volans, while not unstable, is highly maneuverable.It lacks a strong passive righting mechanism, which favors an active gliding behavior.Similarly low pitching moments have been reported in other vertebrate gliders [36,46], and an increase in maneuverability has long been considered to improve gliding or flying performance through active behavior in aerial vertebrates [46,47].In line with this, Khandelwal and Hedrick [19] reported that Draco continuously changes its pitch angle in order to maintain a relatively constant angle of attack.
In summary, our results provide new insight on the presumed gliding behavior of D. volans, such as: (1) The lift generating mechanism of D. volans is strengthened by body camber, and to a lesser degree by adopting an AT posture (figures 3, 7, 8 and S3).An AT posture also provides a means to hold the wings expanded and glide efficiently.This could explain why Draco typically adopts this posture when gliding [15].(2) D. volans could modulate its wing area by controlling their expansion or shape to reduce its instantaneous aerodynamic force (figure 6).This could provide a means to maneuver around obstacles [7,19].(3) Maximizing glide range depends on the length of each glide.For long glides (>30 vertical meters), glide range is maximal if the animal maintains the angle of attack that maximizes the liftto-drag ratio (table 2; figure 12).For shorter glides (e.g.ten vertical meters), higher angles of attack result in more efficient gliding (table 2; figure 13).A higher angle of attack would thus be beneficial for D. volans when aiming for a neighboring branch rather than one far away [6,28].This conforms with the observation that D. volans plans its glide path prior to take off [18] and that the angle of attacks observed in vertebrate gliders do not necessarily maximize their lift-to drag ratio for a given situation [28].(4) Conversely, D. volans may benefit from increasing its gliding speed and initial fall height rather than glide range by adopting a lower (possibly negative) angle of attack, a low degree of camber, and less expanded wings (figures 3, 4, 7, 8, 13 and S3).D. volans indeed uses gliding to escape predators, such as the flying colubrid snake Chrysopelea which is also a capable glider [35].(5) D. volans appears mostly neutral in pitch (figure 3(c)), and is thus highly maneuverable.In line with this, D. volans is known to continuously change its pitch angle, wing shape and posture [18,19].Vertebrate gliders are typically mostly neutral in pitch [36,46], and this improves their gliding performance and maneuverability [46,47].
The work presented here thus provides a better picture of the aerodynamics, gliding performances, and presumed gliding behavior of D. volans.Yet, many factors influencing gliding performance ought to be considered.Further studies may for instance incorporate unsteady CFD simulations to better reconstruct the influence of sequential postural changes in a single glide, or consider the compliance of the wing membranes, as highlighted by direct observations in Draco [19].Yet, the data and approach presented here provide a means of describing the aerodynamics of this enigmatic lizard, for which gliding behavior and performance have only recently started being quantified by experimental and numerical studies [18][19][20][21].Similar methods could also be employed to study gliding flight in other poorly known extant gliders, or even in extinct taxa where direct measurements and observations are impossible.There are indeed a variety of extinct gliding reptiles analogous to Draco [11], such as the Late Triassic kuehneosaurids [48] or the Late Permian weigeltisaurids [49][50][51], for which gliding performance and behavior remains poorly understood.

Conclusion
The flying agamid lizard Draco volans is one of the most renown vertebrate gliders, capable of using its patagial wings to generate lift and glide in its forest canopy habitat.The 3D steady-state CFD simulations presented here demonstrate that while airborne and holding its wings with its forelimbs, D. volans generates a turbulent boundary layer over its wings characterized by a large recirculation cell that is likely maintained attached above the body by the interaction with wing tip vortices.This mechanism induces a strong suction over the wing that enhances lift generation, especially for high camber and angles of attack, and may be modulated by changing wing expansion and shape to control the generation of aerodynamic force.Such a lift enhancing mechanism is typical of other vertebrate gliders, and of aircrafts with low aspect-ratio wings in general.Based on these CFD results, rigid-body gliding trajectory simulations highlight the importance of body camber and orientation in D. volans, as both influence glide range.For instance, increased camber and angle of attack provide a means to increase glide range for short jumps, which may be beneficial in certain situations.Thus, D. volans may actively control its posture to control its aerodynamic and gliding performance depending on the chosen glide path, as previously reported by direct observations.Lastly, D. volans is mostly neutral in pitch and highly maneuverable, similar to other vertebrate gliders.The work presented here thus provides an important first step to the understanding of the lift generating mechanism and the influence of body posture on aerodynamic performance in D. volans, and provides a means to address gliding flight in other extinct or extant gliders for which direct observations are impossible.Further studies may incorporate additional aspect of the morphology and behavior of Draco to provide a clearer picture of the gliding prowess of this emblematic animal.

Figure 1 .
Figure 1.Selected Draco volans geometries used for CFD cases in dorsal, anterior, and left lateral views.(a) P0AF, no camber arms-first geometry; (b) P7AF, low camber arms-first geometry; (c) P15AF, high camber arms-first geometry; (d) close-up view of P0AF geometry with arms flexed at the elbow and hands pointing anteriorly ('arms-first'), illustrating the anatomical axes and orientations relative to the sagittal axis used in the text; (e) close-up view of P0 AT geometry with arms oriented transversely and resting above the leading edge of the wing ('arms-transverse').Scale bars equal 2 cm.

Figure 2 .
Figure 2. Overview of computational domain and mesh: (a) computational domain in oblique view illustrating the position of the geometry relative to the incoming airflow; (b) detail of the mesh on the sagittal plane (Z = 0), illustrating the refinement box around the geometry; (c) close-up oblique view of (b) illustrating the refinement layers close to the body and the discretization of the geometry.

Figure 3 .
Figure 3. Aerodynamic performances for different postural geometries of Draco volans immersed in a unidirectional laminar airflow at U∞ = 5 m s −1 .(a) free-body diagram of gliding Draco volans illustrating glide angle γ, angle of attack α, forces and pitching moment; (b) lift and drag coefficients against α; (c) pitching moment coefficient CM against α; (d) lift-do-drag ratio CL/CD against α.

Figure 4 .
Figure 4. Velocity field around a Draco volans arms-transverse standard gliding posture P0 AT at various angles of attack α, with isoline of null horizontal velocity (Ux = 0 m s −1 ) illustrating airflow separation.Left column, velocity field in anterior view, on transverse (spanwise) plane set in the posterior half of the wing; right column, velocity field in left lateral view, on parasagittal (streamwise) plane set mid-span of left wing.Arrows indicate direction of freestream velocity U∞; gray color indicates freestream velocity U∞ = 5 m s −1 .All cases taken at last time step.

Figure 5 .
Figure 5. Oblique view of scalar fields and vortices around Draco volans arms-transverse standard gliding geometry P0 AT set at an angle of attack of 10 • .(a) Pressure distribution on the surface of the geometry; (b) streamlines representing the velocity field around the left half of the geometry.(c), Q-criterion isosurfaces (Q = 1000) colored by the velocity field.Transparent gray color indicates freestream pressure P∞ = 0 Pa and velocity U∞ = 5 m s −1 .All cases taken at last time step. classical

Figure 6 .
Figure 6.Velocity field around a Draco volans resting posture (P0WF) at different angles of attack, with isoline of null horizontal velocity (Ux = 0 m s −1 ) illustrating airflow separation.Scalar field in anterior view, on transverse (spanwise) plane set in the posterior half of the wing.Arrows indicate incoming airflow direction.Gray color indicates freestream velocity U∞ = 5 m s −1 .Both cases taken at last time step.

Figure 7 .
Figure 7. Velocity field around Draco volans geometries at an angle of attack of 0 • (a) and 20 • (b), with isoline of null horizontal velocity (Ux = 0 m s −1 ) illustrating airflow separation.Velocity fields in left lateral view on parasagittal (streamwise) plane set mid-span of left wing.Arrows indicate direction of freestream velocity U∞; gray color indicates freestream velocity U∞ = 5 m s −1 .All cases taken at last time step.

Figure 8 .
Figure 8. Left lateral view of streamlines representing the velocity field around selected Draco volans posture geometries set at an angle of attack of 20 • .Transparent gray color indicates freestream velocity U∞ = 5 m s −1 .All cases taken at last time step.

Figure 9 .
Figure 9. Velocity field around Draco volans geometries for selected cases, with isoline of null horizontal velocity (Ux = 0 m s −1 ) illustrating airflow separation.(a) Selected arms-transverse geometries at an angle of attack of 20 • ; (b) highly cambered arms-transverse geometry P15 AT set at various angles of attack.Velocity fields in left lateral view, on sagittal plane of animal.Arrows indicate direction of freestream velocity U∞; gray color indicates freestream velocity U∞ = 5 m s −1 .All cases taken at last time step.

Figure 10 .
Figure 10.Gliding simulations from a 30 m perch for different postural geometries of Draco volans starting at an initial horizontal velocity Vx(0) of 2.5 m s −1 .Each geometry set at a fixed angle of attack which maximizes CL/CD.(a) Simulated gliding trajectories; (b) glide angle; (c), velocity.Left column, comparison of all postural geometries; right column, details of the slightly cambered arms-first geometry P7AF illustrating the four gliding phases.Markers indicate time in seconds.Red arrows represent acceleration vectors every 0.2 s.V, total velocity; Vx, horizontal velocity; Vy, vertical velocity.

Figure 11 .
Figure 11.Gliding simulations from a 30 m perch for different postural geometries of Draco volans starting at an initial horizontal velocity Vx(0) of 2.5 m s −1 , illustrating the variation of acceleration through the gliding phases.Each geometry gliding at a fixed angle of attack α which maximizes CL/CD (a), or α = 25 • which maximizes CL for the studied range of α (b).Markers indicate time in seconds.Red arrows represent acceleration vectors every 0.2 s.

Figure 12 .
Figure 12.Variation in gliding simulations for postural geometries P0AF and P0 AT (a), P7AF and P7 AT (b), and P15AF and P15 AT (c) of Draco volans starting at an initial velocity Vx(0) of 2.5 m s −1 .Top, variation due to perch height (10 m, 15 m, 20 m, 30 m); Bottom, variation due to angle of attack α (−10 • to 25 • ) for all perch heights.Large filled markers and solid lines indicate α of (CL/CD)max; large empty markers and dotted lines indicate α of 25• ; smaller markers indicate all other α.Note that the differences between the trajectories reflect only the vertical space available, as the computed trajectories are all identical relative to perch height but stop once the animal reaches the ground.

Figure 13 .
Figure 13.Variation due to angle of attack over the first 10 m of the glide range for a 10 m perch for postural geometries P0AF and P0 AT (a), P7AF and P7 AT (b), and P15AF and P15 AT (c) of Draco volans starting at an initial velocity Vx(0) of 2.5 m s −1 .Large filled markers and solid lines indicate angle of attack α of (CL/CD)max; Large empty markers and dotted lines indicate α of 25 • .

Table 2 .
Aerodynamic and gliding performances for all postures for two-dimensional trajectory simulations starting from different perch heights.The glide phases are identical for all geometries provided sufficient vertical space and are thus only indicated for a 30 m perch.