On the approximation of D.I.Y. water rocket dynamics including air drag
Laura Fischer, Thomas Günther, Linda Herzig, Tobias Jarzina, Frank Klinker, Sabine Knipper, Franz-Georg Schürmann, Marvin Wollek
PPreprint
On the approximation of D.I.Y. water rocket dynamicsincluding air drag
L. Fischer T. G¨unther ,a L. Herzig T. Jarzina F. Klinker , ,b S. Knipper F.-G. Sch¨urmann M. Wollek Marie-Curie-Gymnasium, Billy-Montigny-Platz 5, 59199 B¨onen, Germany TU Dortmund, Fakult¨at f¨ur Mathematik, 44221 Dortmund, Germany Eduard-Spranger-Berufskolleg, Vorheider Weg 8, 59067 Hamm, Germany
Abstract.
If you want to get accurate predictions for the motion of water andair propelled D.I.Y rockets, neglecting air resistance is not an option. Butthe theoretical analysis including air drag leads to a system of differentialequations which can only be solved numerically. We propose an approximationwhich simply works by the estimate of a definite integral and which is evenfeasible for undergraduate physics courses. The results only slightly deviatefrom the reference data (received by the Runge-Kutta method). The motion isdivided into several flight phases that are discussed separately and the resultingequations are solved by analytic and numeric methods. The different resultsfrom the flight phases are collected and are compared to data that has beenachieved by well explained and documented experiments. Furthermore, wetheoretically estimate the rocket’s drag coefficient. The result is confirmed bya wind tunnel experiment.
The classical rocket equation [15, eq. 2.12] is attributed to Konstantin EduardovichTsiolkovsky [23]. If the propellant is exhausted with constant speed the rocketequation can easily be integrated. In case of water and air propelled rockets theexhaust speed decreases together with the internal pressure and mass of the rocket. a [email protected] (corresponding author) b [email protected] Int. J. of Sci. Research in Mathematical and Statistical Sciences (2019) no. 6, 1-13We suggest the reader to choose this preprint version because the journal text contains severalformatting errors. a r X i v : . [ phy s i c s . pop - ph ] J a n reprint L. Fischer et. al. Physics of water rockets have been the subject of several investigations. In an earlypaper [17, p. 152], Nelson and Wilson claim that rocket thrust and mass as a functionof time, as well as the drag coefficient must be determined experimentally. Later,Finney assumes in [9] that “ the air in the rocket behaves as an ideal gas and [...] expands isothermally ”. He uses Bernoulli’s equation to determine an equation for thepropellant’s mass flow rate, see [9, eq. 3]. Furthermore, Finney deduces an equationfor the pressure as a function of time and calculated the rocket’s burn time, see [9,eq. 7, 8]. Therewith, he proposes the estimate h = 18 g t − D M g t (1)for the rocket’s height, where D = ρ air c D A R contains the parameters of the dragforce F D and M is the mass of the empty rocket, see [9, eq. 14]. As mentioned inthe appendix of his paper “ an adiabatic approximation would seem more natural .”Gommes slightly improves the thrust prediction for water rockets: “ the gas expansionhas to be modeled as an adiabatic process. ”, see [10]. He also includes the fact, that“ Air expansion [...] is accompanied by vapor condensation ”. A more meticulousinvestigation of the thermodynamics of the water rocket’s thrust phase was publishedby Romanelli, Bove and Madina in [19]. Indeed, the value of the polytropic exponent n in pV n = constant affects the rocket’s maximum altitude. There are several othereffects that raise the inaccuracy of the predictions more than that. Prusa used n = 1 . n . Since n is a constant, this doesn’t cause additional difficulties.Now let us get back to the rocket’s movement: Prusa derived an equation for therocket’s acceleration, cf. [18, eq. 2.2], and proposed a numerical algorithm to solve it,see [18, p. 723]. Prusa neglected additional thrust from compressed air which is leftat the end of the water propulsion phase. An even more accurate consideration of thewater rocket physics was given by Barrio-Perott et. al. in 2010, see [3]. They workedout a set of differential equations for the water thrust, cf. [3, eq. 18-21], as well asequations for the air thrust [3, eq. 28-31]. Both articles, [18] and [3] use numericalmethods to achieve the solution.In our paper the rocket’s ascent is divided into water thrust phase, air thrust phaseand upward coasting phase. We also set up the equations of motion for the rocketin Section 2. The equation that governs the gas expansion inside the rocket tankis solved analytically in Subsection 3.2, but it is not possible to analytically solvethe complete system of differential equations that describe the rocket’s launch. Firstand foremost, we study the water thrust phase in Section 4. On the one hand, weapply the Runge-Kutta method to the corresponding initial value problem to havesome proper reference data. On the other hand, we deduce a simple method toapproximate the results with high accuracy: Our approximation simply works by reprint On the approximation of D.I.Y. water rocket dynamics the calculation of a definite integral. One advantage is that this calculation canbe done with a simple graphing calculator without a computer algebra system andactually is feasible for undergraduate physics courses. Nonetheless, one receives veryaccurate results which only slightly deviate from the reference data (received bythe Runge-Kutta method). Section 5 is concerned with the air thrust phase. D.I.Y.rockets are chaotic systems and there is no way to receive analytic results. Therefore,we make reasonable assumptions to simplify the underlying equations and introducean efficiency factor for the air thrust. In Section 6 the upward coasting phase isdiscussed. A collection of the previous results is given in Section 7. Subsequently,our paper contains a detailed theoretical discussion of the drag coefficient and weverify the results by experimental data from a wind tunnel experiment, see Section8. Finally, we compare our formula of the rocket’s maximum altitude with theexperimental data from our D.I.Y. rocket launch experiments. Acknowledgments.
We are grateful to Prof. Dr. Andreas Br¨ummer of the Depart-ment of Fluidics at TU Dortmund University for kindly providing us access to thedepartment’s wind tunnel.
Consider a rocket moving with velocity v in vertical direction. Let m be the massof the rocket including its propellant at a given time t . During the interval of time dt , the rocket ejects the mass element dm with the exhaust speed v e . According tothe law of conservation of momentum, the velocity v of the rocket thereby increasesby dv = − v e m − dm , see e.g. [15, eq. 2.11]. Consequently, the thrust acceleration isgiven by a I = − dmdt · v e m . (2)Since the total mass is shrinking, i.e. dmdt <
0, it is a I >
0. A rocket within theterrestrial atmosphere additionally experiences deceleration from gravity and airdrag. Both of them decrease with the rocket’s height. However, a simple waterrocket does only reach a maximum altitude of a few meters . Consequently, thealtitude dependence of gravity and air drag can be neglected. For the magnitude It should be mentioned that a group of scientists at the University of Cape Town builta water and air propelled rocket that reached 830 m reprint L. Fischer et. al. g = GM earth r − of barycentric gravitational acceleration we use g = 9 . ms − duringnumerical calculations. Air resistance is modeled by the drag force F D = ρ air c D A R v where ρ air is the density of air, c D the drag coefficient, A R the reference area (later weuse the rocket’s cross sectional area A cs perpendicular to the direction of movement),and v the rocket’s velocity. Let us adopt the abbreviation D := 12 ρ air c D A R (3)from [9]. The drag force F D = Dv leads to a deceleration a D = F D m = Dm v . (4)An estimate of the drag coefficient of our model rocket is given in Section 8. Deceler-ation from gravity and drag point on the one hand and thrust acceleration on theother hand point in opposite directions. Therefore, the rocket’s total acceleration isgiven by a = a I − g − a D = − v e dmm dt − g − Dm v . (5)For a water rocket, the working mass (water) represents the major part of the rocket’stotal mass m . Anyway, the total mass strongly decreases with time. In order tointegrate (5) it is necessary to know the dependence of mass m and time t . Let (cid:101) m ( t ) be the time dependent working mass, namely the mass water in the rockettank. Further let M denote the constant mass of the empty rocket. The timedependent mass of the compressed air doesn’t contribute significantly to the totalmass of the rocket. Despite this fact, we will take it into account. Section (5) isconcerned with the additional air propulsion after the water thrust phase. Duringthe water thrust phase we treat the mass of compressed air as a constant, say m air . This approach is based on the idea that all water is expelled before the airescapes. Of course, this is somehow physically unrealistic. Indeed, an mixture ofwater and air is expelled, which yields a non-calculable chaos, in particular towardsthe end of the water thrust phase. Within our model the rocket’s mass is given by m ( t ) = M + m air + (cid:101) m ( t ) during the water thrust phase. Let ρ denote the densityof the working mass, in case of water this is ρ ≈ kg m − , and (cid:101) V ( t ) its volume. V b is the volume of the rocket’s tank and V ( t ) = V b − (cid:101) V ( t ) the part which is filledwith air (or another gas, maybe water vapor etc.). Therefore, the mass m and thedifferential dm read m = M + m air + ρ ( V b − V ( t )) ⇒ dm = − ρ dV. (6) reprint On the approximation of D.I.Y. water rocket dynamics Thereby, the dependence of mass m and time t is determined by the gas expansionduring the water thrust phase. Let p a denote the atmospheric pressure and p = p ( V )the pressure of propellant and gas inside the rocket. Initial values at rocket launch forgas pressure and volume are denoted by p and V , respectively. Bernoulli’s equation p = p a + v e ρ/
2, see [22], relates the exhaust speed v e to the pressure difference p − p a : v e = (cid:115) p ( V ) − p a ) ρ . (7)Indeed, the barometric pressure p a slightly depends on the altitude, see [12]. Ataltitudes that can be reached by a water rocket, p a decreases by about 12 P a/m =12 − bar/m . Since the pressure of the propellant will usually be about a few bar wecan neglect the altitude dependence of the barometric pressure. Another point is thatthe fluid’s velocity before leaving the nozzle is assumed to be zero in (7). Strictlyspeaking, we have to take into account the rejuvenation of the bottle. Consideran incompressible fluid (like water) which laminar flows from a point inside thebottle with cross-section A R through the nozzle with cross-section A at velocity w e .Let (cid:101) w denote the velocity of the fluid inside the bottle. The continuity equationleads to (cid:101) w = w e A/A R = w e r /R where R and r are the radii of bottle and nozzle,respectively. Usually, the radius of the nozzle will be small against the rocket’sradius, that is r (cid:28) R . From Bernoulli’s equation p + (cid:0) w e r /R (cid:1) ρ/ p a + w e ρ/ w e = (cid:115) p ( V ) − p a ) ρ · (cid:113) − (cid:0) rR (cid:1) = v e (cid:26) r R + O (cid:18)(cid:104) rR (cid:105) (cid:19)(cid:27) . As shown above , equation (7) represents the third order Taylor approximation w.r. t. r/R for the exhaust velocity w e . But what is the error while using (7) insteadof the latter equation? The radius of a commercially available PET bottle can beestimated to about 4 to 5 cm . Our nozzle has a diameter of about 9 mm . That leadsto R ≈ r and we get w e ≈ . v e . The error for the exhaust velocity is about0 . r = R/
3, which seems to be unrealistic large, theerror will be less than 1%. Obviously, approximation (7) is good enough for ourconcern. The working mass pressure and therewith the pressure difference depend onthe gas volume V . The propellant is expelled through a nozzle with a cross sectionalarea A at speed v e . Let d (cid:126)A be the vector surface element normal to A . In our case,the velocity vector (cid:126)v e is perpendicular to A as well. Therefore, the volumetric flow The corresponding Taylor series is 1 √ − x = 1 + 12 x + O (cid:0) x (cid:1) with x = r/R . reprint L. Fischer et. al. rate through the plane surface A reduces to dVdt = (cid:90) (cid:90) A (cid:68) (cid:126)v e , d (cid:126)A (cid:69) = Av e (8)and from (6) we deduce the mass flow rate dmdt = ddt [ M + m air + ρ ( V b − V ( t ))] = − ρ dVdt (8) = − ρAv e . (9)From (7) one gets ρv e = 2 ( p ( V ) − p a ). With the above considerations the rocket’sacceleration (5) takes the form a = ρAv e m − g − Dm v = 2 A ( p ( V ) − p a ) − Dv M + m air + ρ ( V b − V ) − g (10)where V is the gas volume as a function of time t . In order to solve the latter equationit is also necessary to know the relation p ( V ) of pressure and volume. “ Many real processes undergone by gases or vapours are approximately polytropic witha polytropic index typically between 1.0 and 1.7 [...]”, cf. [8]. It is argued in [10] thatthe gas expansion in a model rocket can also be described by a polytropic process.Pressure p and Volume V are related by pV n = constant . If p and V denote thecorresponding initial values this is p ( V ) = p (cid:18) V V (cid:19) n . (11)Volumetric flow rate (8) and exhaust speed (7) lead to dVdt = A v e = A (cid:115) p ( V ) − p a ) ρ (12)which together with (11) yields dVdt = A (cid:115) p V n V − n − p a ) ρ . (13)In case of V = V b we receive the rocket’s burn time by numerical integration of t b = 1 A (cid:114) ρ (cid:90) V b V dV (cid:112) p V n V − n − p a . (14) reprint On the approximation of D.I.Y. water rocket dynamics In fact, there is still left some compressed air after the whole water is expelled. Thisprovides an additional air thrust. Strictly speaking, (14) gives the duration of thewater-thrust phase. However, this time will be very short. We will discuss the airboost in further detail later in Section 4.
Specific investigation of the gas expansion requires the value of the polytropicexponent. Frequently, the polytropic index is chosen as n = 1 .
4, see [14, 18] forexample. This corresponds to the assumption that the gas in the rocket is dryair. As mentioned in [10] “ the air expansion in the rocket is accompanied by watervapor condensation, which provides an extra thrust ”. The water vapor pressure atwhich water vapor is in thermodynamic equilibrium with the water depends on thetemperature. For moist air the relation of saturation vapor pressure and temperaturecan be well approximated by an equation given by Arden Buck in [5]. If T denotesthe air temperature in ◦ C , the saturated vapor pressure p v is given by p v = K · exp (cid:20)(cid:18) . − T . (cid:19) (cid:18) T .
14 + T (cid:19)(cid:21) ,K = 611 . Nm = 6 . · − bar . (15)In [19, 10] an approximation of the polytropic index that depends on the water vaporpressure is given by n = 1 .
15 + (1 . − .
15) exp (cid:18) − p v p (cid:19) . (16)For example, p = 3 bar and T = 15 ◦ C yield n ≈ .
35. Figure 1 shows the polytropicindex in dependence of temperature and pressure, i.e. the function n ( T, p ) given by(16) together with (15). For the expansion of moist air in water rockets the polytropicindex is about 1 . ≤ n ≤ . First, the question arises if there is an analytical solution for (13). If that is not thecase, the corresponding equation can be solved numerically. The following approachis based on the Gaussian hypergeometric function. A little manipulation of theordinary differential (13) for the volume yields V n dVdt = A (cid:115) p V n ρ (cid:114) − p a p V n V n . reprint L. Fischer et. al. Figure 1: Polytropic index n ( T, p ) where temperature varies from 0 ◦ C to 40 ◦ C .The pressure range is 0 ≤ p ≤ bar and 0 ≤ p ≤ bar , respectively.Together with ddt (cid:104) V n +22 (cid:105) = n +22 V n dVdt the latter equation leads to ddt (cid:34)(cid:18) p a p V n (cid:19) n +22 n V n +22 (cid:35) = ( n + 2) A (cid:18) p a p V n (cid:19) n +22 n (cid:115) p V n ρ (cid:118)(cid:117)(cid:117)(cid:117)(cid:116) − (cid:34)(cid:18) p a p V n (cid:19) n +22 n V n +22 (cid:35) nn +2 . By using u := (cid:16) p a p V n (cid:17) n +22 n V n +22 , k := ( n + 2) A (cid:16) p a p V n (cid:17) n +22 n (cid:115) p V n ρ , α := 2 nn + 2we finally get du √ − u α = k dt . (17)This equation admits the solution u · H (cid:0) , α ; α + 1; u α (cid:1) = kt + ξ (18)where ξ = u · H (cid:0) , α ; α + 1; u α (cid:1) and u = u (0). Here H ( a, b ; c ; x ) = ∞ (cid:80) i =0 Γ( a + i )Γ( b + i )Γ( c )Γ( a )Γ( b )Γ( c + i ) i ! x i is a special generalized hypergeometric function aka Gausshypergeometric function. More precisely, this function is a special solution of Euler’shypergeometric differential equation, see [4, eq. 1.498], x ( x − w (cid:48)(cid:48) + (cid:0) ( a + b + 1) x − c ) w (cid:48) + ab w = 0 reprint On the approximation of D.I.Y. water rocket dynamics which can also be written as (cid:0) x ddx + a (cid:1)(cid:0) x ddx + b (cid:1) w = (cid:0) x ddx + c (cid:1) w (cid:48) . The function H ( a, b ; c ; x ) obeys ddx (cid:0) x c − H ( a, b ; c ; x ) (cid:1) = ( c − x c − H ( a, b ; c − x ) , (19) H ( a, b ; b ; x ) = (1 − x ) − a , (20)see [21, eq. 1.4.1.6] and [21, eq. 1.5.1]. Therefore, by writing g ( x ) = x α H (cid:0) , α ; α +1; x (cid:1) and ˜ g ( u ) = g ( x ( u )) = uH (cid:0) , α ; α + 1; u α (cid:1) with x ( u ) = u α we get d ˜ gdu = dgdx (cid:12)(cid:12)(cid:12) x = u α · dxdu (19) = 1 α x α − H (cid:0) , α ; α ; x (cid:1)(cid:12)(cid:12)(cid:12) x = u α · αu α − = H (cid:0) , α ; α ; u α (cid:1) (20) = 1 √ − u α as stated. Although a solution of (13) is represented by (18), its implicit form isinexpedient for our concerns. Of course, the burn time t b can be calculated from(18). But we also need a closed-form expression for the gas volume as a functionof time. Hence, the expansion equation (13) will be included into the system ofequations of motion for our rocket. However, let us take a look at the separatenumerical solution of (13) now. Within our calculations with the Runge Kuttamethod the temperature of water and vapor is chosen to be 15 ◦ C . The volume ofthe bottle and the initial water volume are 1 dm and 0 . dm , respectively. Thecorresponding Figure 2 was created with wxMaxima . The source code is available athttps://github.com/tguent/code.Figure 2: Expansion of water vapor from 0 . dm to 1 dm at 15 ◦ C . The left sideshows the expansion with an initial pressure of 3 bar . This leads to a polytropicexponent of n ≈ .
35. In the other case, the initial pressure is 10 bar , which leads to n ≈ .
39. The polytropic exponent was determined by (16). reprint L. Fischer et. al. The flight of the rocket can be divided into three parts: Thrust phase, coasting phasewith upwards motion, and free-fall phase. We model the thrust phase by dividing itagain into two parts: The main thrust is provided by the expelled water. Duringthe water thrust phase the acceleration of the water rocket is modeled by (5). Butusually there is still some compressed air left at the end of the water thrust phase.This gives rise to an additional nonzero air thrust. We refer to this as the air thrustphase in the following. After all propellant (water and air) is exhausted, the rocketcan be regarded as an object that is thrown vertically upwards. The correspondinginitial velocity is determined by the velocity at the end of the thrust phase. Afterreaching the maximum altitude, the rocket enters the free-fall phase. During theupward coasting and free-fall phase air drag has a crucial influence on the rocket’smovement. As we will see in the following, air drag is negligible during the thrustphase of a water rocket.
Remark 1 ( Example data).
For our numerical calculations we use the followingdata: The empty rocket has a mass of 1 / kg and a volume of 1 dm , its nozzle has adiameter of 9 mm . The rocket radius is 4 cm . Initial values for propellant volumeand pressure are V = 0 . dm and p = 3 bar at a temperature of 15 ◦ C . The massof the compressed air is approximately 1 . g . The air drag coefficient was set to c D = 1. As mentioned in Section 3, pressure p and volume V are related by polytropicexpansion (11) where V as a function of time t is determined by (13). The rocket’svelocity can be expressed as the rate of change of its position h by v = dhdt . Analogously,the acceleration (10) can be expressed as the rate of change of its velocity a = dvdt .Finally, this leads to the following first order system of differential equations for gasvolume V , altitude h and velocity v : dVdt = A (cid:115) p V n V − n − p a ) ρ , dhdt = v ,dvdt = 2 A ( p V n V − n − p a ) − Dv M + m air + ρV b − ρV − g . (21)The latter can be solved numerically. We use the fourth order Runge Kutta method,see [1], which is implemented in the computer algebra system wxMaxima , see [24].The step-size was set to 0 .
01. The corresponding numerical results for the referencedata given in Remark 1 are presented in Figure 3. In order to create comparability reprint On the approximation of D.I.Y. water rocket dynamics we use the following reference values v b ≈ . ms at an altitude of h b ≈ . m atthe end of the water thrust phase. For the corresponding wxMaxima source code seeagain https://github.com/tguent/code.Figure 3: Velocity and altitude of the rocket during the water thrust phase for thereference data, see Remark 1. At the end of the water thrust phase the rocket has avelocity of v b ≈ . ms at an altitude of h b ≈ . m . Solving the equations of motion (21) require advanced numerical calculations. In thefollowing we present a more simple estimation of altitude and velocity at the endof the water thrust phase. The first equation of (21) determines the gas expansion.Figure 2 shows the gas expansion for 0 ≤ t ≤ t b , where t b is determined by (14).This indicates that the gas volume as a function of time is not too far from beinglinear, at least as a rough estimate. If we use V ( t ) = V b − V t b t + V (22)for our calculations the corresponding equation of motion predicts a velocity ofabout 15 . ms at an altitude of 2 . m at the end of the thrust phase (again for thedata given in Remark 1). This is only a minor deviation downwards from the moreaccurate reference values of v b ≈ . ms and h b ≈ . m . The relative deviationof velocity and altitude at the end of the thrust phase is about 0.9% and 0.2%, With (22) the equations of motion (21) reduce to dhdt = v , dvdt = A p V n (cid:18) Vb − V tb t + V (cid:19) n − p a − Dv M + m air + ρV b − ρ (cid:16) Vb − V tb t + V (cid:17) − g . reprint L. Fischer et. al. respectively. Thus, the linear approach (22) can be used as a good approximation.Now consider the last equation of (21) which incorporates thrust and air resistance.Although air drag has crucial influence on the maximum altitude of the rocket ithas few influence during the thrust phase. During the subsequent free flight phasewhen the propellant is exhausted it is very important to incorporate air drag again.Next we calculate the data at the end of the water thrust phase in case of c D = 0.Together with (22), equation (21) leads to v n ( t ) = (cid:90) t A (cid:34) p V n (cid:16) Vb − V tb x + V (cid:17) n − p a (cid:35) M + m air + ρV b − ρ (cid:16) V b − V t b x + V (cid:17) dx − gt (23)at a given time t ≤ t b . Integration over the burn time t b , see (14), gives the velocity v nb at the end of the water thrust phase. The altitude at the end of a water rocket’sthrust phase is small. We use the mean acceleration ¯ a = v nb /t b and h nb := 12 ¯ at b = v nb t b v nb ≈ . ms for h nb ≈ . m . Theresults exceed the reference values v b ≈ . ms and h b ≈ . m by 1 .
8% and 4 , As mentioned before, the mass of air inside the rocket is small, i.e. m air (cid:28) M . Evenif we neglect the mass of the compressed air during the water thrust phase the errorwill be very small. Using (11) the mass of air inside the rocket can be calculatedfrom m air = ρ air n (cid:114) p p a V (25)where ρ air ≈ . gdm is the density of air at atmospheric pressure p a ≈ bar . Forexample: An initial air volume of V = 0 . dm at p = 3 bar contributes to thetotal mass of the rocket with about 2 g . In our case, the rocket has a curb weightof about 1 / kg plus 350 g propellant (water). The additional mass of 2 g leads to a reprint On the approximation of D.I.Y. water rocket dynamics deviation of about some centimeters of calculated altitude. I.e., the additional massof the compressed air doesn’t has a noticeable influence on the water thrust phasewhere the rocket is propelled by the water. But what effect has the expelled air afterthe water thrust phase? On the one hand, there is only a tiny working mass leftto induce a change in momentum. On the other hand, the air escapes very rapidly.Indeed, it already was a matter of discussion if the air boost can be neglected or not.Prusa notes: “ Due to its low density, the thrust provided by air alone is negligible,and in a launch without water, the rocket is barely able to lift off of the air pumpseal. ”, cf. [18, p. 719]. On the other hand, Barrio-Perotti et. al. state: “
This air isexpelled through the nozzle causing an additional increase in the rocket momentumthat sometimes cannot be neglected: for example, a rocket with air pressurized at 2bars can reach a distance of about 10 m [...]”, cf. [3, p. 1138]. In [3], the air propulsionis also described by a system of differential equations, see [3, eq. 29, 30, 31] togetherwith the algebraic equation [3, eq. 28] which have to be solved numerically. In orderto decide whether the air boost has to be taken into account or not we filled ourmodel rocket with compressed air (no water) at about 2 bar pressure. Actually, therocket lifted off and reached a significant altitude of at least two meters. Thus wedecided not to neglect the air propulsion. But since the effect is small compared tothe water thrust, we choose a straightforward approximation of the air thrust phasein the following.
The gas volume equals the bottle’s volume V b now. According to (11), initial pressure p and initial air density have reduced to p b = p (cid:18) V V b (cid:19) n and ρ b = m air V b (25) = ρ air V V b n (cid:114) p p a (26)The change in momentum can be calculated from the rocket equation (5). In orderto avoid complicated numerical methods we neglect the air drag during the air thrustphase. Barrio-Perotti et. al. evaluate water and air thrust in [3]. In comparison withthe duration of the water thrust phase, the air thrust phase is short. Beyond that,the rocket acceleration is strongly decreasing during the air propulsion, see [3, Fig.12]: If the air propulsion lasts about 0 .
025 seconds, the major contribution may beover within a hundredth of a second. Compared with this, the water thrust in [3, Fig.12] lasts six times longer (and leads to higher acceleration anyway). Based on theshort duration of the air propulsion, we assume that the exhaust velocity v e can betreated as a constant. With these simplifications the rocket equation dvdt = − v e dmm dt − g reprint L. Fischer et. al. can be integrated. (cid:90) v max v b dv = − v e (cid:90) MM + m air dmm − (cid:90) duration of air propulsion g dt. (27)Now let us rate the two terms on the right side of (27). Because the change ofmass can be huge the integral over the mass may have a significant contributionto the result, despite of the fact that the duration of the air propulsion is veryshort. Now consider the last integral. The product of gravitational acceleration andduration of air propulsion will be also small. So we neglect this part. The remainingterms result in v max ≈ v e ln (cid:16) M + m air M (cid:17) + v b where v b is the rocket velocity at the endof the water thrust phase. Indeed, the above considerations are further based onthe assumption that no air escapes before all water is expelled. This seems to beunrealistic. Particularly towards the end of the water thrust phase the propellantwill be a mixture of water and air. Therefore, the mass of the remaining air at theend of water thrust will be less than the initial amount of air inside the rocket givenby (25). In order to take into account the loss of air during the water thrust phasewe include an efficiency factor η and receive v max = v e ln (cid:18) M + η · m air M (cid:19) + v b . (28)Amongst other things, the efficiency factor η presumably depends on the initialproportion of water and air in the rocket, that is V /V b . If the initial amount of airequals the bottle volume, there is no water propulsion and the air provides the wholethrust. Obviously the air thrust efficiency should be η = 1 in this case. If the wholebottle is filled with water, η = 0 should indicate that there will be no air thrust.Furthermore, the efficiency strongly depends on the construction of the D.I.Y rocket.Since we consider a highly chaotic system, there should be a parameter to readjustour calculation to the experimental data of a special model rocket. Based on thepreceding considerations we choose η = (cid:18) V V b (cid:19) µ . (29)As a first estimate one can use µ = 1. In order to fine-tune this estimate for a specialD.I.Y. rocket, the exponent µ ≥ µ ≈ v e . We cannot simply calculate v e from (7),since it is based on Bernoulli’s equation for incompressible flow. The density of acompressible flow is not constant. But the derivation of Bernoulli’s equation from reprint On the approximation of D.I.Y. water rocket dynamics dpρ + v dv + g dz = 0 can be modified for compressible flow. With p ∝ ρ n , see [7,eq. 3.20], we receive the relationship nn − · pρ + v = constant which applies tocompressible adiabatic flows, see [7, eq. 3.21]. We use again the above notation ρ air for the air density at atmospheric pressure p a . Density and pressure at the beginningof the air thrust phase are denoted by ρ b and p b , respectively. Let us assume thatthe speed inside the bottle can be neglected. Therewith, we receive for the exhaustvelocity: v e = (cid:115) nn − (cid:18) p b ρ b − p a ρ air (cid:19) (26) = (cid:118)(cid:117)(cid:117)(cid:116) n p a ( n − ρ air (cid:34)(cid:18) p p a (cid:19) − n (cid:18) V V b (cid:19) n − − (cid:35) . (30)Combining (28) and (30) results in v max ≈ (cid:118)(cid:117)(cid:117)(cid:116) n p a ( n − ρ air (cid:34)(cid:18) p p a (cid:19) − n (cid:18) V V b (cid:19) n − − (cid:35) ln (cid:18) M + η · m air M (cid:19) + v b . (31)Since the duration of the air propulsion is very short, we neglect the additional heightwhich the rocket attains during this phase. But we take into account the enhancedvelocity from the air thrust. This significantly affects the upward coasting. By the time the working mass is expelled the rocket has already reached its topspeed. From this point we can regard the rocket as an object that is thrown verticallyupwards. Air resistance has a significant impact on the movement after the thrustphase. Luckily, the corresponding equations of motion can be solved analytically.The kinematics of an upward movement including air resistance are well known. Asthe rocket’s propellant is exhausted the remaining mass M of the rocket is constant.With the notation ψ = DM = ρ air c D A R M , (32)the air drag deceleration (4) takes the form a D = ψv during the free flight phase,see (3). For v ≤ v max the equation of motion dvdt = − (cid:0) g + ψv (cid:1) leads to the integralequation (cid:90) vv max dxg + ψx = − (cid:90) t dt where t = 0 corresponds to the end of the thrust phase. The solution is given by v ( t ) = (cid:114) gψ tan (cid:32) arctan (cid:32)(cid:115) ψg v max (cid:33) − (cid:112) ψg t (cid:33) . (33) reprint L. Fischer et. al. With v = 0 the duration of the upward coasting phase is given by t c = 1 √ ψg arctan (cid:32)(cid:115) ψg v max (cid:33) . (34)The upward coasting begins after the thrust phase and ends when the rocket hasreached the maximum altitude. The covered distance during the free flight phase isgiven by h c = (cid:90) t c (cid:114) gψ tan (cid:32) arctan (cid:32)(cid:115) ψg v max (cid:33) − (cid:112) ψg t (cid:33) dt = 1 ψ (cid:40) ln (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) cos (cid:32) arctan (cid:32)(cid:115) ψg v max (cid:33) − (cid:112) ψg t c (cid:33)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) − ln (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) cos (cid:32) arctan (cid:32)(cid:115) ψg v max (cid:33)(cid:33)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:41) = − ψ ln (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) cos (cid:32) arctan (cid:32)(cid:115) ψg v max (cid:33)(cid:33)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) . The first term in the second row cancels out by using (34) for t c . Usingcos (arctan x ) = √ x , i.e. ln | cos (arctan ( x )) | = ln (cid:18) √ x (cid:19) = −
12 ln (cid:0) x (cid:1) , one finally gets h c = 12 ψ ln (cid:18) ψg v (cid:19) (35)where ψ = ρ air c D A R M , see (32). Let us summarize the previous results and therewith derive a method to estimatethe maximum altitude of the rocket. First we have to calculate the burn time of thewater thrust phase by t b = 1 A (cid:114) ρ (cid:90) V b V dV (cid:112) p V n V − n − p a , see (14). Here A is the nozzle cross sectional area, ρ the water density, n thepolytropic exponent, p a the atmospheric pressure and p the initial pressure in the The absolute value function can be abandoned since x > reprint On the approximation of D.I.Y. water rocket dynamics rocket tank. The tank has the volume V b . At launch, V < V b is filled with air and V b − V is the initial water volume. The water thrust phase can be modeled by thesystem of differential equations (21). Its numerical solution gives velocity v b andaltitude h b at the end of the water thrust phase. Alternatively, for a rough estimatewe receive these values from v b ≈ v nb = (cid:90) t b A (cid:34) p V n (cid:16) Vb − V tb x + V (cid:17) n − p a (cid:35) M + m air + ρV b − ρ (cid:16) V b − V t b x + V (cid:17) dx − gt b ,h b ≈ h nb = v b t b , (36)see (23) and (24). M is the rocket’s curb mass, m air the initial mass of air in thetank and g the gravitational acceleration. After the water thrust phase, the rocketexperiences an additional acceleration due to the air propulsion. We only considerthe enhanced velocity (31) within our model v max ≈ (cid:118)(cid:117)(cid:117)(cid:116) n p a ( n − ρ air (cid:34)(cid:18) p p a (cid:19) − n (cid:18) V V b (cid:19) n − − (cid:35) ln (cid:18) M + η · m air M (cid:19) + v b where η is the efficiency factor of air propulsion, see (29). The rocket enters theupward coasting phase with the velocity v max at the altitude h b . During the upwardcoasting the rocket reaches an additional height of h c = Mρ air c D A R ln (cid:18) ρ air c D A R M g v (cid:19) , see (35). Finally, the maximum altitude of the rocket can be calculated from h max = h b + h c . The drag coefficient can be divided into its components which arise from pressuredrag and friction drag, see [25, eq. 7.63]. Friction drag is caused by the viscosity ofthe surrounding air in our case. Pressure drag is the “ difference between the highpressure in the front stagnation region and the low pressure in the rear separatedregion [...]”, cf. [25, p. 448]. Basically, the drag coefficient varies with the Reynoldsnumber Re = v L/ν where v is the rocket’s velocity, L is its characteristic length, reprint L. Fischer et. al. and ν is the kinematic viscosity of the surrounding atmosphere, see [25, eq. 7.61].A simple home-made rocket built from a PET bottle achieves a maximum speed ofabout 20 ms . From Re = v L/ν we see that the Reynolds number will not exceed 5 · in this case. Of course, a more professional constructed water rocket might receive amaximum Reynolds number which is about ten times higher. The relation of dragcoefficient and Reynolds number is usually obtained from laboratory experiments,see for example [25, fig. 7.16]. Since such precise aerodynamic considerations arebeyond the scope of this paper we try to get an estimate of a constant drag coefficientof our rocket in the following. “ The drag analysis of rockets [...] is usually simplifiedby considering the rocket to be made up of several simple basic components ” [11]. Wewill confine our considerations to the drag analysis of nose cone, body tube, base,fins, and a (small) constant value for the launch lugs. The latter segmentation of amodel rocket is shown in [11, fig. 17]. Due to interference drag the total drag of therocket amounts to more than the sum of the components: “[...] additional amountof drag is caused by the joining of the fins to the rocket body. [...]
Interference dragcan be as much as above the sum of the fin and body tube drag ”, cf. [11, p. 9].In this paper, the drag coefficients of nose cone, body tube, base, fins, interferenceand launch lugs are denoted by c nose , c tube , c base , c fin , c int , and c lau respectively. Thetotal drag coefficient c D = c nose + c tube + c base + c fin + c int + c lau represents the casethat the rocket is moving at zero angle to the wind direction. Any nonzero angleleads to an additional induced drag. The nose cone is exposed to pressure drag and skin friction drag. As mentioned in[11, p. 10] a flat nose cone (solely) would result in a drag coefficient of c N = 0 . A cs be the cross-sectional area of the body tube and A ws the wetted The kinematic viscosity of air at 15 ◦ C is ≈ , · − m s It is worth mentioning that the drag force F D = ρ air c D Av increases with increasing speed ofthe rocket. Although if the drag coefficient c D usually decreases. reprint On the approximation of D.I.Y. water rocket dynamics surface area of the rocket. Therewith, [11, eq. 8] takes the form C l := c nose + c tube = 1 . c f (cid:0) Ld (cid:1) A ws A cs (37)where c f is the skin friction coefficient and L/d the length to diameter ratio of therocket.
Flow separation causes low pressure at the rear of the rocket which results in basedrag. An equation that estimates the base drag is given by c base = 0 . √ c nose + c tube , (38)see [11, eq. 9]. Fins cause additional friction drag, pressure drag and induced drag.But they substantially enhance the flight stability of our rocket. Generally, the findrag depends on various parameters like the thickness to chord ratio, planform areaand cross-sections. We present a rough estimate in this paper. In order to get someupper limit we use fairly pessimistic assumptions for the zero lift fin drag coefficient.Following [11, p. 17], the fin thickness to chord ratio will rarely be greater than0 . cross-section for a pessimistic estimate.Let τ denote the thickness to chord ratio and c ∗ fin the zero lift fin drag coefficient inthe following. The plotted data in [11, fig. 40] suggests that c ∗ fin depends linearly on τ in for 0 . < τ < . c ∗ fin ≈ . τ (rectangular cross-section) , (39) c ∗ fin ≈ . τ (rounded cross-section) (40)The left side of figure 4 shows that the graphic illustration of (39) and (40) reproducesthe corresponding lines in [11, fig. 40]. Let us assume in the following that the finshave rectangular cross-section. The zero lift fin drag coefficient c ∗ fin is based onthe fin area A fin whereas the other drag coefficients are based on the body tubecross-sectional area A cs . Therefore, we have to adjust the coefficient c ∗ fin by c fin = c ∗ fin A fin A cs (39) = 0 . A fin A cs τ. (41) According to [11, fig. 41], the zero lift fin drag coefficient for steamlined cross-sections at ratio0 . .
019 (laminar) to 0 .
024 (turbulent) for 30 ms (100 ftsec ). reprint L. Fischer et. al. Figure 4: Left side: Zero lift fin drag coefficients for fins with rectangular and roundedcross-sections. Right side: Skin friction coefficient as a function of Reynolds numbergiven by Prandtl’s law, see [25, eq. 7.43].As above mentioned, rocket body and fins jointly cause additional intereferenz drag.Due to [11, eq. 21], the interference drag can be estimated by c int = c ∗ fin C R d A cs · N (39) = 0 . C R d A cs N τ , (42)where C R is the root chord of the fin, d the diameter of the body tube, A cs again itscross-sectional area, and N is the number of fins. The more accurate drag calculationgiven in [11] also includes drag from launch lugs. The examples given in [11, p. 44,49] lead to a small launch lugs drag coefficient of 0 .
02 to 0 .
03. For our rough estimatewe incorporate this kind of drag by adding c lau = 0 . Collecting the above discussion, we receive a base value for the drag coefficient ofthe model rocket by summation of c nose , c tube , c base , c fin , c int and c lau , see (37), (38),(39), (41), and (42). However, we have not considered the surface texture of therocket so far, “ roughness can cause drag to increase by about percent ”, c.f. [11,p. 43]. Since our rocket’s nose consists of half a tennis ball, fins and launch lugsare fixed with hot glue and tape, we presuppose that the surface might be ratherrough. This leads to a significant uncertainty which we incorporate as an error in thefollowing. Finally, we receive a more or less adequate formula for the drag coefficientof our D.I.Y. rocket by c D = (cid:20) C l + 0 . √C l + 0 . τA cs (cid:18) A fin + C R d (cid:19) N + 0 . (cid:21) · (1 . ± . reprint On the approximation of D.I.Y. water rocket dynamics with C l from (37). Beside the ingredients that enter into C l , τ is the thickness tochord ratio of the fins, C R the root chord of the fin, and N the number of the fins.The skin friction coefficient c f depends on the kind of air flow and varies greatlywith the Reynolds number, see right side of Figure 4. As c f depends on the rocket’sspeed it would be the best to include the skin friction into the system of differentialequations (21). However, our estimate of the total drag coefficient seems to be toovague to justify this approach. Therefore, we decided to use a constant value for theskin friction coefficient in the following calculations. We get an estimate of the drag coefficient by (43). The maximum speed of our rocketis about 17 ms . The characteristic length of our rocket is 0 . m . Hence, a kinematicair viscosity of 1 , · − m s (at 15 ◦ C ) leads to a maximum Reynolds number ofabout 4 · . Based on Prandtl’s law [25, eq. 7.43], a mean value for the skin frictioncoefficient on this scale is given by¯ c f = 14 · (cid:90) · . Re / x dRe x ≈ . . Indeed, this value is in accordance with the skin friction coefficient at a Reynoldsnumber of 4 · given in [11, fig. A-1]. Our rocket has a diameter of 8 cm which leadsto a cross-sectional area of A cs ≈ . m . The rocket’s surface area (taken as acylinder with patched-up half of a tennis ball ) amounts to A ws ≈ . m . The rockethas N = 3 fins. Each has a root chord of 7 . cm and a fin area of A fin ≈ . m .The fin-thickness-to-chord ratio is about 0 .
08. The drag coefficient (43) is c D = 0 . ± . , for a wxMaxima source code see again https://github.com/tguent/code. Prof. Dr. Andreas Br¨ummer of the chair of Fluidics at TU Dortmund Universitykindly provided us the opportunity to check our estimate of the drag coefficient in awind tunnel. The drag force was measured with a dynamometer with an error of0 . N . We always had to do an additional measurement for the calibration. Thus,in the worst case scenario our error adds up to ∆ F D = 0 . N . The wind speed wasgiven to the first decimal. Hence, we use an error of ∆ v = 0 . ms in the followingconsiderations. The drag coefficient c D is related to drag force F D and wind speed v Drag coefficients of a non spinning tennis balls have been measured to 0 . ± .
05, see [13]. reprint L. Fischer et. al. by c D = F D ρ air A R v . Therefore, the error of the drag coefficient can be calculated by∆ c D = c D (cid:16) ∆ F D F D + 2 ∆ vv (cid:17) .Wind speed v/ ms Drag Force F D /N Drag coefficient c D . ± .
05 0 . ± . . ± . . ± .
05 0 . ± . . ± . . ± .
05 0 . ± . . ± . . ± .
05 0 . ± . . ± . . ± .
05 0 . ± . . ± . . ± .
05 0 . ± . . ± . . ± .
05 0 . ± . . ± . . ± .
05 0 . ± . . ± . . ± .
05 0 . ± .
10 0 . ± . . ± .
05 0 . ± . . ± . . ± .
05 0 . ± . . ± . c D = 0 . ± . , see Table 1. In fact, this is a remarkable low value for our D.I.Y. rocket. However, fromFigure 5 we see that the stream line flow pattern suggests very good aerodynamicsproperties.Figure 5: Stream line flow pattern made visible by smoke in the wind tunnel at theDepartment of Fluidics at TU Dortmund University. reprint On the approximation of D.I.Y. water rocket dynamics The rocket launch was repeated with different initial values for pressure, amountof water, and temperature (of water). The model rocket which was used for theexperiment has weight 143 g . It was constructed from a bottle with a volume of1000 ml . At launch the rocket was filled with an amount of water V l and air V atpressure p and temperature T . The difference to atmospheric pressure p a ≈ bar is p − p a . The pressure p was measured with an error of ∆ p = 0 . bar . The erroris based on the manometer’s scale. The rocket has a radius of 40 mm . We drew ascale for the water level which is legible up to about ± mm . Within the relevantscope the rocket is cylindric. Accordingly, the initial amount of water V l – andtherewith the air volume V – is determined up to an error of ∆ V = 10 ml . Due tothe precision of temperature measurements the error of temperature doesn’t affect theresults. Determining the rocket’s altitude h mes makes up the major part in measuringinaccuracy within our experiment. We placed a measuring rod beside the launchingdevice. The whole flight phase was filmed and we provide links to the correspondingvideos in Table 4. By evaluating the videos we receive the rocket’s maximum altitudeby comparing to the length of the measuring rod. The extrapolation leads to anerror of ∆ h med = 0 . m . The measured altitudes can be confirmed by the given videolinks.It is the aim of this section to compare the experimental data with the theoreticalresults, see Table 3. We use the mean values of the experimental data for thetheoretical results. These are calculated as described in Section 4.1. Thereby, h calc is calculated by using the system of differential equations (21) for the thrust phase.We get h est by calculating the part of the altitude that is obtained from the waterthrust phase from (36). This can be done with a simple graphing calculator. Therocket’s drag coefficient was set to c D = 0 .
6. For the air thrust efficiency factor weuse µ = 3 in (29).A wxMaxima source code for the calculations from Table 3 is available athttps://github.com/tguent/code. T / ◦ C V l /ml V /ml p /bar h mes /m ±
10 650 ±
10 3 . ± .
05 17 . ± . m ±
10 675 ±
10 3 . ± .
05 17 . ± . m ±
10 745 ±
10 3 . ± .
05 14 . ± . m ±
10 745 ±
10 3 . ± .
05 17 . ± . m . ± . > m ±
10 750 ±
10 2 . ± .
05 7 . ± . m Table 2: The results from the rocket experiment reprint L. Fischer et. al. T / ◦ C V l /ml V /ml p /bar h calc /m h est /m . . m . m . . m . m . . m . m . . m . m . . m . m . . m . m Table 3: Theoretical predictionsvideo link1 https://youtu.be/YXGb75eOqbI2 https://youtu.be/9XY 6iSHXPE3 https://youtu.be/3l5JJHO0Tp84 https://youtu.be/RUkJwDCq8uU5 https://youtu.be/MifjHZH3Z7Y6 https://youtu.be/8MRufGWCUjMTable 4: Video links
10 Conclusion
The method introduced in this paper yields accurate theoretical predictions. Thepredictions in launches 1, 2 and 4 from Table 4 lie within the experimental error,see Table 2. Launch 5 was hard to evaluate because the rocket crashed into theceiling. Indeed, the main goal of launch 5 was to show that there is significant thrustif the rocket is filled with pressured air only. Launch 3 yielded an altitude of 13 . m at the lower bound which differs from the predicted value of 13 . m . The questioncomes up whether the error in altitude has to be enlarged. On the one hand, thecamera position has to be far enough from the experiment that the elevation angleremains small. On the other hand, huge distance leads to more blurred pictures.It stands to reason that one has to find a more suitable method to evaluate themaximum altitude. In case of launch 6, theoretical and experimental results don’tcoincide. Indeed, we launched the rocket several times without getting exploitabledata because the rocket didn’t lift off correctly. In some cases the rocket stuck toolong to the launching device due to friction. In other cases the rocket lurched throughthe air. Furthermore, our model D.I.Y. rocket began to leak after a huge amount ofexperiments.Putting it all together, our method is suitable to predict the altitude in case that therocket lifts off perfectly. But even in this case the water rocket physics represents a reprint On the approximation of D.I.Y. water rocket dynamics highly chaotic system. On this basis, the simple estimation proposed in this paperyields amazingly good results. References [1] M. Abramowitz, I. A. Stegun:
Handbook of mathematical functions with formulas,graphs, and mathematical tables . Dover Publications, 1965[2] R. Barrio-Perotti, E. Blanco-Marigorta, K. Arguelles-Diaz, J. Fernandez-Oro: Experi-mental Evaluation of the Drag Coefficient of Water Rockets by a Simple Free-Fall Test.
European Journal of Physics no. 5, 1039-1048 (2009)[3] R. Barrio-Perotti, E. Blanco-Marigorta, J. Fernndez-Francos, M. Galdo-Vega: Theo-retical and experimental analysis of the physics of water rockets. European Journal ofPhysics no. 5, 1131-1147 (2010)[4] I. N. Bronstein: Taschenbuch der Mathematik . B. G. Teubner Stuttgart-Leipzig, 1996[5] A. L. Buck: New equations for computing vapor pressure and enhancement factor.
J. Appl. Meteorol. no. 12, 1527-1532 (1981)[6] T. A. Campbell, M. Okutsu: Model Rocket Project for Aerospace EngineeringCourse: Trajectory Simulation and Propellant Analysis. Preprint 2017, arXiv:1708.01970[physics.ed-ph][7] L. J. Clancy: Aerodynamics . John Wiley & Sons, 1975[8] M. Clifford (ed.):
An Introduction to Mechanical Engineering . CRC Press Taylor &Francis Group, 2006 (Part 1), Hodder Education, An Hachette UK Company, 2010(Part 2)[9] G. A. Finney: Analysis of a water-propelled rocket: A problem in honors physics.
American Journal of Physics no. 3, 223-227 (2000)[10] C.J. Gommes: A more thorough analysis of water rockets: Moist adiabats, transientflows and inertial forces in a soda bottle. American Journal of Physics no. 3, 236-243(2010) DOI: 10.1119/1.3257702[11] G. M. Gregorek, Aerodynamic Drag of Model Rockets . Estes Industries, Penrose, CO,1970[12] International Civil Aviation Organization:
Manual of the ICAO Standard Atmosphere .Doc 7488/3, 3rd ed. 1993[13] R. Mehta, F. Alam, A. Subic: Review of tennis ball aerodynamics.
Sports Technology no. 1, 7-16 (2008) DOI: 10.1080/19346182.2008.9648446[14] D. Kagan, L. Buchholtz, L. Klein: Soda-bottle water rockets. Phys. Teach. Raumfahrtsysteme . Springer Vieweg, Berlin, Heidelberg, 2017, DOI: 10.1007/978-3-662-49638-1 2[16] J. Moran,
An Introduction to Theoretical and Computational Aerodynamics . John Wiley& Sons, New York, 1984[17] R. A. Nelson, M. E. Wilson: Mathematical analysis of a model rocket trajectory Part I:The powered phase.
Phys. Teach. no. 3, 150-161 (1976) reprint L. Fischer et. al. [18] J. M. Prusa: Hydrodynamics of a Water Rocket. Siam Rev. no. 4, 719-726 (2000)[19] A. Romanelli, I. Bove, F. G. Madina: Air expansion in the water rocket. AmericanJournal of Physics no.frm[o]–0, 762-766 (2013) DOI:10.1119/1.4811116[20] Simple Drag Tests for Water Rockets - seeds2lrn.com. http://fliphtml5.com/rftx/obtf(visited 05.03.2019)[21] L. J. Slater: Generalized Hyperbolic Functions . Cambridge University Press, 1966[22] H. St¨ocker:
Taschenbuch der Physik . Verlag Harri Deutsch, 1998[23] K. E. Tsiolkovsky: The Exploration of Cosmic Space by Means of Reaction Devices (inRussian).
The Science Review (1903)[24] A. Vodopivec: wxMaxima 18.02.0. http://andrejv.github.io/wxmaxima/[25] F. M. White: Fluid Mechanics . McGraw-Hill Education Ltd, 7th ed., 2011. McGraw-Hill Education Ltd, 7th ed., 2011