Figures
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 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.
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; CD, drag coefficient; g, Earth’s gravity; m, mass; Ma, Mach number; Re, Reynolds number; RS, 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 CD, from the computed product CDA⊥, 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 CD, 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 CD 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 vt would be independent of altitude and velocity, the acceleration would cease at v = vt. So, vt is the terminal velocity at instantaneous flight conditions.
Analytical solutions for vt = const
Let v0 = v(t0 = 0) and h0 = h(t0 = 0). If vt is considered constant, then Eq (3) can be integrated directly
We finally get
or equivalently
(4)
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 ≪ vt, we derive after some lengthy approximation
This simplifies for an initial zero speed v0 = 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 vt ≠ const
Yet, in reality vt 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 Rs = 286.91 J kg−1K−1 is the specific gas constant of standard atmosphere and g0 = 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 = h0 − v0t we derive additionally
So, with definitions t0 = 0, v(t0) = v0 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 CDA⊥ 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/v0-term in square brackets accounts any variation of acceleration to an adjustment of the aerodynamic parameter vt0 from its trivial determination as given in the round brackets, and the v0Δ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 CDA⊥ 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 CD 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 Rs = 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.
III. Data analysis procedure
For deriving the drag coefficient CD from Eq (12) we need to have the following data
- Flight trajectory data v(t) and h(t)
- Atmospheric data, i.e. air density ρ(h) and air temperature data T(h).
- Pressure suit data, i.e. the total mass m, incl. body weight, and the effective surfaces along his bodies symmetry axes Ax,Ay,Az
- Flight attitude data, i.e. angle of attack α(t) of his body against flight direction, to derive the wetted area A⊥(t)
Flight trajectory data
Felix jumped from the capsule at an altitude of 38969.4 m. At 50.0 s into flight and at an altitude of 28833 m he obtained his maximum vertical velocity of 377.1 ms−1 equaling Mach 1.25. His free fall lasted down to an altitude of 2566.8 m where he pulled the drag chute. The detailed flight trajectory data v(t) and h(t) were extracted from the graphical plots as given in the Summary Report of the Red Bull Stratos Team [9] by digitalization (see Fig 1). Velocity and time were measured by an on-body GPS system.
The gray area on the left indicates that for t < 17 s, owing to insignificant drag, no drag coefficient can be derived. The gray area on the right indicates times when Felix experienced a dangerous flat spin at a reduced AOA.
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
(14)
with Rs = 286.91 J K−1kg−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
(15)
This confirms the expected result that for a belly-down attitude A⊥(α = 0) = Ax.
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 m2 for area reference. From these the following effective suit cross sections were derived
(16)
The ledge is of size 17”x11” to determine its cross-section. (Credit: Art Thompson / Stratos team).
The total mass of Felix and the fully dress-up suit was determined to be
(17)
Gravity of earth
For the gravity of Earth, we apply the WGS 84 ellipsoidal gravity formula in combination with a FAC altitude correction factor
with altitude h and geographical latitude φ. For φ = 33.39° at the jump site this reduces to
(18)
Hence, the altitude-dependent gravity introduces a relative g error of only 0.01%.
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 t0 of the first relevant time interval
With the given pressure suit data and assuming an initial CD = 0.65 and belly-down flight attitude, i.e. A⊥ = Ax, we get
Assuming and from the atmospheric density profile ρ(h), we finally get for the first data analysis interval
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
(19)
We chose Eq (10) for interval estimation. For that we assume as given and specify Δv = δv. With this, and Δh ≈ v0 ⋅ Δ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
(20)
with
given by Eq (8). If the time interval becomes so short that Δh < 200 m we set Δh = 200 m.
IV. Results
Analysis methods
To extract the aerodynamic parameter CDA⊥ 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 CDA⊥ 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 h0,v0;h,v this equation is solved for vt by the trust-region method with dogleg step strategy. We recall that this method B assumes that vt and hence the atmospheric density ρ is constant over each interval.
Method C makes use of Eq (12), which is explicit in CDA⊥ 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 t0 = 17 s into flight corresponding to an altitude of h0 = 37 640 m. Thereafter we used for the time interval estimator, Eq (20), with δv = −10 m/s down to hend = 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.
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).
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)
V. Sensitivity analysis
Different methods
In Fig 5 the drag coefficients are provided as derived from the three different methods A, B, C (see Section Results –Analysis Methods above). 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.
The drag coefficient as derived from the three different methods A, B, and C.
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).
t0 = 17.1 s was finally chosen from theoretical considerations.
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.
Intervals are given by the different speed changes δv within the interval (see Eq (20)).
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° and 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 the standard variation σ. The value of the fitting parameter is depicted on the abscissae and their quantity on the ordinate. The gray area in Fig 4 shows the ±σ uncertainty of the distribution function. It can be seen from the error in Eqs (24) that for increasing Mach numbers the error due to the AOA uncertainty increases, in particular for the slope C of the line for Ma > 1.1, because it shows with σ = 0.26 the highest variation. For the derived ±σ the corresponding error of the fitting parameters are
(24)
AOA sensitivity analysis for fit parameters A, B, and C. For each trajectory time step the AOA was randomly sampled 1000 times with variations up to 20%. The abscissae show the value of a fitting parameter binned in seven value ranges. The vertical bar to each range measures the number of values obtained from the variations that fall into this range. From the resulting histogram, the standard variation σ for each of the fitting parameters A, B, and C is derived.
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.
The combined error of the data statistics Eq (23) and the AOA sensitivity analysis Eq (24) is
(25)
This continuous function together with the total error as given by Eq (25) is displayed in Fig 4 as a dashed line and the gray region.
Compared to theoretical expectations for subsonic speeds as given by Eq (21) we have with CD = 0.60 ± 0.05 a perfect match between our data and theoretical expectations. In addition, our overall pattern, is compliant with an increase of CD below the sonic speed and a decrease above it. However, for transonic speeds we see an increase of the drag coefficient of just ΔCD = 0.14, which is only 19% of the theoretically expected ΔCD = 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.
(Credit: Art Thompson / Stratos team).
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.
The representative physical length scale of 1 m and a constant dynamic viscosity of 15×10−6 kg s−1 m−1.
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 CD from the Reynolds number below the critical Mach number. In our analysis, the constant A in Eq (22) models this constant pressure drag.
(from S.F. Hoerner [11]).
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.
References
- 1.
GoPro: Red Bull Stratos—The Full Story. Available: https://www.youtube.com/watch?v=dYw4meRWGd4. Cited 1 September 2016
- 2.
Sky Dive From The Edge Of Space (1960). Available: https://www.youtube.com/watch?v=Z8RRkMsHOMU. Cited 1 September 2016
- 3.
Markoff J, Parachutist’s Record Fall: Over 25 Miles in 15 Minutes, in The New York Times, ed, 2014. Available from https://mobile.nytimes.com/2014/10/25/science/alan-eustace-jumps-from-stratosphere-breaking-felix-baumgartners-world-record.html, Cited 18 April 2017
- 4. Mohazzabi P, Shea JH. High-altitude free fall. American Journal of Physics 1996; 64:1242–1245.
- 5. Benacka J. High-altitude free fall revised. American Journal of Physics 2010; 78: 616–619.
- 6.
Vadim R, John J, Sarah F. High Altitude Free Fall: Theoretical Analysis of G-Forces Impacting Human Subjects. in 41st International Conference on Environmental Systems, ed: American Institute of Aeronautics and Astronautics, 2011
- 7.
Rygalov VY, Greene JK, Jurist JM. Supersonic Survival Envelope Analysis for High Altitude Free Fall. 44th International Conference on Environmental Systems, 13–17 July 2014, Tucson, Arizona, ICES-2014-108.
- 8. Colino JM, Barbero AJ. Quantitative model of record stratospheric freefall. European Journal of Physics 2013; 34: 841–848.
- 9.
Stratos RB. Summary Report: Findings of the Red Bull Stratos Scientific Summit, California Science Center, Los Angeles, California, USA, 2013. Available:https://issuu.com/redbullstratos/docs/red_bull_stratos_summit_report_final_050213
- 10.
Anderson JD. Fundamentals of Aerodynamics, ed: McGraw-Hill; 2011
- 11.
Hoerner SF. Fluid-dynamic drag practical information on aerodynamic drag and hydrodynamic resistance, 2nd ed. Midland Park, NJ; 1965.