Abstract In this paper, we analyze the trajectory and body attitude data of Felix Baumgartner’s supersonic free fall through the atmosphere on October 14, 2012. As one of us (UW) was scientific advisor to the Red Bull Stratos team, the analysis is based on true body data (body mass, wetted pressure suit surface area) and actual atmospheric data from weather balloon measurements. We also present a fully developed theoretical analysis and solution of atmospheric free fall. By matching the flight data against this solution, we are able to derive and track the drag coefficient C D from the subsonic to the transonic and supersonic regime, and back again. Although the subsonic drag coefficient is the expected C D = 0.60 ± 0.05, surprisingly the transonic compressibility drag coefficient is only 19% of the expected value. We provide a plausible explanation for this unexpected result.

Citation: Guerster M, Walter U (2017) Aerodynamics of a highly irregular body at transonic speeds—Analysis of STRATOS flight data. PLoS ONE 12(12): e0187798. https://doi.org/10.1371/journal.pone.0187798 Editor: Christof Markus Aegerter, Universitat Zurich, SWITZERLAND Received: April 21, 2017; Accepted: October 26, 2017; Published: December 7, 2017 Copyright: © 2017 Guerster, Walter. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited. Data Availability: All relevant data are within the paper. Funding: This work is the result of a bachelor thesis. Therefore, the authors received no specific funding for this work. Competing interests: The authors have declared that no competing interests exist. Abbreviations: A ⊥ , wetted pressure suit surface area; C D , drag coefficient; g, Earth’s gravity; m, mass; Ma, Mach number; Re, Reynolds number; R S , specific gas constant; t, time; T, temperature; v, velocity; h, altitude (height); α, angle of attack (AOA); φ, degree of latitude

I. Introduction On October 14, 2012, and as part of the Red Bull Stratos Mission [1], Felix Baumgartner performed a record-breaking supersonic free fall from the stratosphere. Other than from earlier jumps, namely by Joe Kittinger in 1960 [2] and a supersonic jump by Alan Eustace on October 24, 2014 [3], only for Felix’ extensive flight data are available. It is the objective of this paper to understand the physics of the free fall, in particular the aerodynamic behavior in the transonic regime. The first analysis of a high-altitude free fall dates back to 1996, when Mohazzabi and Shea [4] were able to analytically solve the equation of motion for v(h) by power series including the atmospheric friction term and a standard barometric law. From a linear approximation of an implicit solution they were able to roughly analyze Kittinger’s jump confirming his key achievements. However, owing to the linear approximations made, the result becomes less reliable for jumps with extensive free fall segments as with Felix Baumgartner’s and Alan Eustace’s jumps. In addition, the power series approximation works only with a strict barometric law and a constant aerodynamic friction term. In 2010 Benacka [5] revised Mohazzabi’s high-altitude free fall analysis by taking into account a linear temperature gradient as about actual in the troposphere. For a drag coefficient that does not depend on the Mach number he was able to solve the equation of motion analytically by power series, again for v(h). Later in 2011 Vadim et al. [6] and 2014 Rygalov [7] published two papers. The earlier deals with the theoretical analysis of the maximal drag deceleration in free fall and the corresponding altitude, whereas the latter shows that an emergency egress at an altitude around 100 km can be survived and it provides a relation between the initial altitude and the altitude of the transonic transition. His analysis is based on a drag coefficient independent on the Mach number and on a temperature-independent barometric law. J. M Colino and A. J Barbero in 2013 [8] published a course quantitative data analysis based on a spreadsheet and intended as an introductory physics course of Felix’ free fall. Although they use a similar approach as described in our paper, they do not have precise data and use only a standard atmosphere model to derive v(h) and v(t). Therefore they were not able to extract the key aerodynamic parameter, the drag coefficient C D , from the computed product C D A ⊥ , with A ⊥ being the wetted surface area. Since the insight into the physics of the supersonic jump is gained only through a well-defined drag coefficient C D , a full-fledged numerical investigation with exact data is needed, which is the objective of this paper.

II. Theoretical analysis Let h be the altitude of the skydiver above ground. We define the downwards velocity v to be positive, i.e. dh = −v⋅dt. Then the equation of motion of a body with mass m subject to gravitational and aerodynamic forces is known to be (1) where C D is the aerodynamic drag coefficient, including the pressure drag, skin friction drag and interference drag. The coefficient quite generally depends on the type of flow (laminar or turbulent) and hence on the Mach number Ma and Reynolds number Re. We can safely neglect a Knudsen number dependency because at stratospheric altitudes we are in the continuum flow regime. The quantity ρ is the atmospheric density depending on altitude and temperature T. A ⊥ is the aerodynamically wetted pressure suit area depending on the angle of attack α as defined later. We assume g to be the gravity of Earth and not Earth’s gravitational acceleration. This is because Earth’s gravity results jointly from Earth’s gravitational force plus the centrifugal force due to the rotation of the Earth and therefore is dependent on the altitude h and geographical latitude φ. Because of Earth’s atmosphere up to roughly 100 km altitude on average co-rotates with Earth’s surface, we have to take centrifugal forces into account. We first define the instantaneous terminal velocity (2) With this the acceleration from Eq (1) can be written as (3) If v t would be independent of altitude and velocity, the acceleration would cease at v = v t . So, v t is the terminal velocity at instantaneous flight conditions. Analytical solutions for v t = const Let v 0 = v(t 0 = 0) and h 0 = h(t 0 = 0). If v t is considered constant, then Eq (3) can be integrated directly We finally get or equivalently (4) We rewrite this equation as By taking v = −dh/dt into account this equation can directly be integrated to deliver under the initial precondition (5) For a small initial time interval, gt ≪ v t , we derive after some lengthy approximation This simplifies for an initial zero speed v 0 = 0 to the convenient expression and By the same token we now derive v(h). We first substitute In the equation of motion (3). We hence get and and finally (6) Eq (6) is a result of Mohazzabi and Shea [4] with generalized initial conditions. Analytical solution for v t ≠ const Yet, in reality v t is not constant. In particular, and according to the barometric formula the atmospheric density decreases drastically with altitude and also somewhat with air temperature via the scale height where R s = 286.91 J kg−1K−1 is the specific gas constant of standard atmosphere and g 0 = 9.798 m s−2. With this we rewrite the equation of motion Eq (3) as (7) where (8) We note for later purposes that in free fall air drag always counteracts gravitational force and hence (9) Eq (7) can no longer be integrated fully analytically. However, for flight data analysis we only need to consider the change in velocity Δv for small time intervals Δt. In this case the analytical solution is given by the Taylor series By Eq (7) we have not only , but with the approximation h = h 0 − v 0 t we derive additionally So, with definitions t 0 = 0, v(t 0 ) = v 0 we obtain and , which inserted into the Taylor series delivers (10) Extraction of aerodynamic parameters from flight data To derive aerodynamic parameters from flight data, we have to define a small-time interval Δt, then read from the fight data the measured change in velocity Δv in that interval, and with this finally solve Eq (10) for to derive C D A ⊥ from via Eq (8). Because the Δt2-term is much smaller than the Δt-term, the Banach fixed-point theorem applies and we can solve Eq (10) for iteratively by contraction mapping. So, in a first step we set the Δt2-term to zero and find as a first approximation The latter follows from Eq (9). This is the trivial finding that attributes any deviation of constant acceleration from g to atmospheric drag. For a second-order approximation, we insert this result into the Δt2-term of Eq (10), which yields Again, solving for delivers (11) This is the second order solution of Eq (9) for a given Δt and a measured Δv. Analysis of the above derivation reveals that the Δv/v 0 -term in square brackets accounts any variation of acceleration to an adjustment of the aerodynamic parameter v t0 from its trivial determination as given in the round brackets, and the v 0 Δt/(2H)-term, which accounts for air density increase in the interval, makes an exception from this. We insert this advanced solution into Eq (8), and with the notation that the index 0 indicates values at the beginning of a time interval, we finally derive C D A ⊥ for any small time interval (12) with Transonic regime The objective of our work is to understand the aerodynamic behavior of a human body with pressure suit in the transonic regime by deriving the drag coefficient C D with Stratos’ flight data from Eq (12). The transonic regime is the velocity range around the speed of sound a or Mach number Ma = 1 with (13) Here κ air = 1.403 is the adiabatic index and R s = 286.91 J K−1kg−1 is the specific gas constant of the standard atmosphere. Since Ma is temperature dependent (but not density-dependent as one might assume intuitively) an atmospheric temperature profile T(h) is essential to determine the correct Mach number.

IV. Results Analysis methods To extract the aerodynamic parameter C D A ⊥ from the data, we applied three different methods. Method A is straightforward in that starting with a measured initial vertical speed it solves the equation of motion (1) and C D A ⊥ is adjusted such that the vertical speed at the end of each time interval fits the measured value. Method B makes use of the equation preceding Eq (6) For a given set of measured interval data h 0 ,v 0 ;h,v this equation is solved for v t by the trust-region method with dogleg step strategy. We recall that this method B assumes that v t and hence the atmospheric density ρ is constant over each interval. Method C makes use of Eq (12), which is explicit in C D A ⊥ and in addition takes into account also linear variations of air density and scale height variations over each interval. Overall, the results of all three methods did not show any grave qualitative differences. However, because method A and method B featured some numerical instabilities in particular at high altitudes while method C was not only stable throughout the data range, but is also physically more elaborated, we discuss here only the results derived by method C. For a detailed comparison of the three methods see section Sensitivity Analysis–Different Methods below. Method C is based on the approximation v⋅Δt = Δh ≪ 2H ≈ 14 km. A lower Δh-limit arises from increasing relative velocity gain errors with decreasing Δh. We hence settled at an interval width of Δh = 200 m. According to the above section Initial Interval Considerations, the first interval started at t 0 = 17 s into flight corresponding to an altitude of h 0 = 37 640 m. Thereafter we used for the time interval estimator, Eq (20), with δv = −10 m/s down to h end = 13040 m. Fig 4 displays the resulting aerodynamic drag coefficient as a function of Mach number. We note that the absolute inaccuracy of the data is mainly driven by the inaccuracy of the area A(AOA). Therefore, we performed a sensitivity analysis described in Section V. PPT PowerPoint slide

PowerPoint slide PNG larger image

larger image TIFF original image Download: Fig 4. Drag coefficient as a function of Mach number. Red circles for increasing and blue crosses for decreasing velocities. The dashed line is the final result as given by Eq (26) and the gray region the total uncertainty as given by Eq (25). The green line is the theoretical transonic drag coefficient for a rotating cube according to Hoerner [11] and as given by Eq (21). https://doi.org/10.1371/journal.pone.0187798.g004 There are two distinct phases (cf. full line in Fig 1): A velocity increase phase for 17s < t < 50s (red circles) and a velocity decrease phase 50 s < t < 120 s (blue crosses). In addition, a green line is shown, which according to Hoerner [11] is the transonic drag coefficient for a rotating cube, namely (21) This empirical dependency was used by the Red Bull team to predict the maximum velocity at a given exit altitude, or the required altitude to ensure that his maximum speed would achieve Mach 1.15 vice versa and hence achieve supersonic speed of a free falling human body for the first time. The scatter of data point and equivalently the data error obviously is much bigger for the increasing velocity data points than for the decreasing velocity data points. This is because the drag D ∝ ρv2 affects the measured speed with decreasing altitude and increasing speed. Nevertheless, the Mach dependency of the up-velocity and down-velocity data points coincide. We therefore now focus on the more accurate down-velocity data (blue crosses). To derive an empirical drag coefficient dependence on the Mach number, we use the same intervals and dependencies as Eq (21), but with the three variable fit parameters A, B, and C: (22) From fitting the downstream data, we derive the fit parameter as (23)

VI. Discussion In order to grasp this unexpected result, we apply fluid dynamics theory (Standard literature for further information on aerodynamics can be found in reference [10] and [11].). The key elements are the overall shape of the body and the relevant Reynold numbers. Relative to the wind flow, Baumgartner’s body is quite blunt, where the pressure suit with its folds adds roughness, while the technical equipment and cameras cause an unevenness of the surface (see Fig 9). Moreover, the backpack and Baumgartner’s limbs form even a quite complex body shape. We therefore identify Baumgartner’s irregular blunt body very approximately by a cylinder with surface roughness on all length scales. PPT PowerPoint slide

PowerPoint slide PNG larger image

larger image TIFF original image Download: Fig 9. The complex suit shape evoking an equally complex aerodynamics. (Credit: Art Thompson / Stratos team). https://doi.org/10.1371/journal.pone.0187798.g009 The Reynolds number Re (see Fig 10) describes the proportion between the inertial and viscous forces and hence the type of flow, laminar or turbulent. It is between 2.3 ⋅ 105 − 12.2 ⋅ 105 for the considered time zone 29 s < t < 68 s. PPT PowerPoint slide

PowerPoint slide PNG larger image

larger image TIFF original image Download: Fig 10. Reynolds number as a function of time into flight. The representative physical length scale of 1 m and a constant dynamic viscosity of 15×10−6 kg s−1 m−1. https://doi.org/10.1371/journal.pone.0187798.g010 Drag crisis Fluid dynamic tells us that blunt bodies, i.e. bodies which are not streamlined, are subject dominantly to pressure drag and only little to skin friction drag. Because they are blunt they undergo a so-called drag crisis (see Fig 11) for Re ≈ 105. Drag crisis is a drop of drag by up to a factor of 10 for smooth surfaces like a table tennis ball. The drop is due to the fact that the laminar flow from the stagnation point to the separation point, located at maximum body cross-section) becomes turbulent, which moves the separation point of the flow backward (delayed separation). This reduces recirculation behind the blunt body and hence lowers drag. However, if the surface is rough (sand-grain size k increases), turbulence in the upstream flow is induced already at lower Reynolds number, causing a less pronounced drag drop (see Fig 11). If the surface is very rough, uneven, and the body’s shape is even not well defined, as in our case, we do not expect any drag crisis because we have turbulence all over the surface at any Reynolds numbers Re > 104. We therefore do not expect, and in fact do not see, any dependency of C D from the Reynolds number below the critical Mach number. In our analysis, the constant A in Eq (22) models this constant pressure drag. PPT PowerPoint slide

PowerPoint slide PNG larger image

larger image TIFF original image Download: Fig 11. Drag crisis of cylinders as a function of surface sand-grain roughness k against cylinder diameter d. (from S.F. Hoerner [11]). https://doi.org/10.1371/journal.pone.0187798.g011 Wave drag At the critical Mach number local supersonic flow and hence shock waves start to emerge from major surface discontinuities (edges) where the flow is forced to accelerate. This erratic shock formation and general flow instabilities correspond to pressure discontinuities, which transform a considerable part of the impinging flow energy into heat. This loss in energy corresponds to a counteracting force, the so-called wave drag (a.k.a. transonic compressibility drag). However, for Mach number up to 0.7 or 0.8 compressibility effects have only minor effects on the flow pattern and drag. For increasing Mach numbers the surface extent that generates shock wave increases and the shock waves intensify. They now interact with the boundary layer so that a separation of the boundary layer occurs immediately behind the shock. This condition accounts for a large increase in drag, which is known as shock-induced (boundary-layer) separation. This separation empirically causes a roughly quadratic increase of the wave-drag coefficient with Mach number. Generally, the term B(Ma − 0.6)2 in Eq (22) models the wave drag. Once a supersonic flow has been established, however, the flow stabilizes and the drag coefficient is reduced. We then essentially have a shock wave that builds up at the upstream front of the body and another one at the unsteady geometry at the downstream tail of the body. It is usually modeled by the term – C(Ma − 1.1) in Eq (22). The shock waves generate the characteristic double boom of a supersonic body flying past, which was indeed captured for Baumgartner’s supersonic free fall in this video recording provided here http://www.youtube.com/watch?v=yZFz6y4UCuo. Given this fluid dynamics effects it becomes clear why drag crisis and the wave-drag coefficient is so insignificant: Because the body is quite blunt, the pressure drag seems to be dominant all the way up to sonic speed. The surface roughness and the complex body shape at all speeds seem to induce a highly turbulent flow even close to the surface. This destroys any boundary layer, thus strongly suppressing a delayed flow separation effect and shock-induced boundary-layer separation.

Acknowledgments The authors thank the whole STRATOS team, in particular Art Thompson, for support and making the flight data available to us, and Albert Pernpeintner for valuable comments on the data interpretation.