Aerodynamics of a highly irregular body at transonic speeds—Analysis of STRATOS flight data

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 CD from the subsonic to the transonic and supersonic regime, and back again. Although the subsonic drag coefficient is the expected CD = 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.


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]  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 welldefined 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 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 With this the acceleration from Eq (1) can be written as 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 rewrite this equation as By taking v = −dh/dt into account this equation can directly be integrated to deliver under the initial precondition 1=v 2 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 By the same token we now derive v(h). We first substitute In the equation of motion (3). We hence get Eq (6) is a result of Mohazzabi and Shea [4] with generalized initial conditions.
Analytical solution for v t 6 ¼ 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 −1 K −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 We note for later purposes that in free fall air drag always counteracts gravitational force and hence 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 _ v, but with the approximation h = h 0 − v 0 t we derive additionally and v 0 , which inserted into the Taylor series delivers 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 v 2 0 =v 2 t0 to derive C D A ? from v 2 t0 via Eq (8). Because the Δt 2 -term is much smaller than the Δt-term, the Banach fixed-point theorem applies and we can solve Eq (10) for v 2 0 =v 2 t0 iteratively by contraction mapping. So, in a first step we set the Δt 2 -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 Δt 2 -term of Eq (10), which yields 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 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 Here κ air = 1.403 is the adiabatic index and R s = 286.91 J K −1 kg −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.

III. Data analysis procedure
For deriving the drag coefficient C D from Eq (12) we need to have the following data

Flight trajectory data v(t) and h(t)
2. Atmospheric data, i.e. air density ρ(h) and air temperature data T(h).  1). Velocity and time were measured by an on-body GPS system. Atmospheric data. On October 12, at 4h 51min GMT, a radiosonde was launched from the airfield measuring the profile data p(h), T(h), wind direction, and wind speed. Based on this measurement a FIM forecast for flight day October 14 was generated. From this and with the ideal gas law, the atmospheric density profile was derived as with R s = 286.91 J K −1 kg −1 . Both profiles, ρ(h) and T(h), are depictured in Fig 1 with errors of about 0.1%.

Flight attitude data and wetted area
The angle of attack (AOA), α(t), is defined (see Fig 2) such that α = 0 is a regular "belly down" attitude for skydivers. We derived α(t) (depicted in Fig 1) from the full story video published by GoPro [1]. The error of the AOA such determined is estimated to be 15-20%.
The wetted pressure suit surface area is determined from the AOA and roll angle ϕ (angle around the body z-axis) quite generally as Because Felix' attitude after 17 seconds of free fall did not show any significant roll, i.e. ϕ = 0, we have This confirms the expected result that for a belly-down attitude A ? (α = 0) = A x .

Pressure suit data
In order to determine the pressure suit data, on September 17, 2012, during a test run, pictures of the pressurized suit were taken (see Fig 3) along the three body frame axes as defined in Fig  2, always together with a sheet of paper of size 17@ × 11@ = 0.1206 m 2 for area reference. From these the following effective suit cross sections were derived The total mass of Felix and the fully dress-up suit was determined to be

Initial interval considerations
Eq (12) is applicable to flight data only if in the equation the aerodynamic term is significantly larger than the velocity data error δv, i.e.
For our digitalization technique δv % 0.05 Á v and hence This is particularly important for the initial free fall segment, where ρv Á Δt is extremely small. In this regime where Δt = t and v % g Á t we get as the condition of the beginning t 0 of the first relevant time interval With the given pressure suit data and assuming an initial C D = 0.65 and belly-down flight attitude, i.e. A ? = A x , we get rðt 0 Þ t 2 0 ¼ 1:6 kg s 2 m À 3 Assuming h 0 % 1 2 gt 2 0 and from the atmospheric density profile ρ(h), we finally get for the first data analysis interval t 0 ¼ 17 s

Sampling interval estimator
Because the contribution of the drag term to the velocity is strongly varying over altitude and hence time, we need to adapt the sampling interval width. It needs to be chosen such that the drag deceleration δv shall be minimum 3 times as much as the data imprecision equaling 1 m/s dv À 3 m=s ð19Þ We chose Eq (10) for interval estimation. For that we assume v 2 0 =v 2 t0 as given and specify Δv = δv. With this, and Δh % v 0 Á Δt, and discarding at high drag the gravitational term we obtain In first-order approximation this yields In second-order approximation we have

And hence for time interval estimation
with v 2 t0 given by Eq (8). If the time interval becomes so short that Δh < 200 m we set Δh = 200 m.

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 trustregion 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.
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 for Ma 0:6 C D ¼ 0:50 þ 1:23ðMa À 0:25Þ 2 for 0:6 Ma 1:1 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 / ρv 2 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: In all three cases, the dashed line is the data function as given by Eqs (22) with (23), i.e. as derived from method C. The data shows that although all three methods provide qualitatively the same results, method C obviously delivers the best result.

Initial interval width
Next, we study the effect of the width of the omitted initial interval (see section Data Analysis Procedure -Initial Interval Considerations) on the results. As shown in Fig 6 we have varied the initial interval in a reasonable range. Besides minor differences, the data points are qualitatively the same (particularly the downstream data point, to which Eq (22) was fitted).

Sampling interval width
The sampling interval is expected to have more effect on the results than all the other parameters. Fig 7 shows the data points for three different velocity increase intervals δv. According to Eq (20) a changing δv results in different time intervals. Obviously, the choice δv = −10m/s from theoretical considerations seems optimal.

AOA uncertainty
To study the impact of the relatively high AOA uncertainty on the resulting drag coefficient as given by Eq (22) with Eq (23), we performed a sensitivity analysis by randomly varying the AOA for each time step by up to 20% around its value shown by in Fig 1 (AOA limited to 0å nd 90˚). We ran 1000 samples and calculated for each value the fitting parameters A, B, and C. A histogram of them is shown in Fig 8 with Note, that the mean of the sensitivity analysis matches our calculated fitting parameters for A and B. For the slope C they do not match exactly. This is due to the limitation of the AOA to 0˚and 90˚and the high sensitivity of this fitting parameter. Compared to theoretical expectations for subsonic speeds as given by Eq (21) we have with C D = 0.60 ± 0.05 a perfect match between our data and theoretical expectations. In addition, our overall pattern, is compliant with an increase of C D below the sonic speed and a decrease above it. However, for transonic speeds we see an increase of the drag coefficient of just ΔC D = 0.14, which is only 19% of the theoretically expected ΔC D = 0.74.

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.
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 Á 10 5 − 12.2 Á 10 5 for the considered time zone 29 s < t < 68 s.

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 % 10 5 . 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 crosssection) 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 > 10 4 . 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.

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 (boundarylayer) separation. This separation empirically causes a roughly quadratic increase of the wavedrag 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.