The water rocket: Thrust from air

After the water is expended, the tank (bottle) still contains substantial pressure that releases quickly through the nozzle. This final burst of air can impart a significant boost to velocity (at least 30% depending on the mass of the rocket), so we shouldn’t ignore this contribution to thrust. For convenience, however, we ignore any further affects on air temperature due to water vapor although we still calculate pressure changes adiabatically using our adjusted value of \(\lambda\). Any water vapor would likely condense by the time airflow begins.

Thrust from choked airflow

When the ratio of ambient pressure to total absolute tank pressure is less than the "choke ratio" $$\alpha_c = \frac{P_a}{P_e} = \left(\frac{2}{\lambda+1}\right)^\frac{\lambda}{\lambda-1} \tag{A1}$$ then the outflow is choked, or limited, to the speed of sound: $$c = \sqrt{\lambda RT} \tag{A2}$$ where

  • \(c\) = speed of sound in air, approximately 343 m/s at 20°C, 331 m/s at 0°C
  • \(R\) = ideal gas constant divided by gas molecular weight (8.3145 / 0.02896 for air)
  • \(T\) = constant gas temperature in Kelvin (e.g. 0°C = 273 K)

The choke ratio is \(\alpha_c = 0.455\) for moist air, using \(\lambda=1.34\); the traditional value is 0.528 using \(\lambda=1.4\) for dry air.

We can use a value of \(T\) from equation (W6), resulting from the expansion of the air that ejected all of the water. The presence of water vapor would have a moderating effect on the final temperature and pressure, so we use the modified value of \(\lambda\) (1.34 instead of 1.4) to calculate the speed of sound. The range of possible values of \(c\) turns out to make just a small difference to apogee altitude, however, typically less than a meter.

To figure out the mass flow rate of air during choked flow, we first need to find the density of compressed humid air in the tank. To do this, we first get the saturation vapor pressure \(P_\text{sat}\) from the Arden Buck equation already given previously in equation B0.

We assume that the air in the bottle is saturated to 100% relative humidity. The density of humid air is given by $$\rho_\text{ humid~air} = \frac{P_\text{d}}{R_\text{d} T} + \frac{P_\text{v}}{R_\text{v} T} \tag{A3}$$ where \(\rho_\text{ humid~air}\) has units of kg/m3, and all of the terms on the right are known:

  • \(P_\text{d} = P_0 - P_\text{v}\) = partial pressure of dry air, or the difference between normal dry air pressure (atmospheric pressure) and water vapor pressure (Pa)
  • \(R_\text{d} = 287.058\) = specific gas constant for dry air in J/(kg·K)
  • \(T\) = air temperature (K)
  • \(P_\text{v} = P_\text{sat}\) = pressure of water vapor (Pa)
  • \(R_{v} = 461.495\) = specific gas constant for water vapor in J/(kg·K)

From this, we can calculate the total initial mass of air and water vapor available for thrust: $$m_{a0} = \frac{P_e}{P_a}\rho_\text{ humid~air}V\tag{A4}$$ where \(V\) is the total volume of the bottle, and \(P_e\) from equation (W5), the absolute pressure of the air in the tank emptied of water, is already adjusted for adiabatic cooling during the expansion that provided the water thrust.

The density of the air in the bottle is simply the mass/volume ratio: $$\rho_0 = \frac{m_a}{V}\tag{A5}$$ where \(m_a = m_{a0}\) initially. This is the mass of air remaining in the bottle at each time step.

The mass flow rate for choked flow is: $$m'=C_dA\sqrt{s \rho_0 P_e} \tag{A6}$$ where in addition to the previous definitions, \(s = \text{constant} = \lambda \left( \frac{2}{\lambda+1} \right)^\frac{\lambda+1}{\lambda-1} = 0.4548\) for \(\lambda = 1.34.\)

\(C_d\) is a discharge coefficient, a correction factor to account for the ratio of the actual discharge from the nozzle to the ideal discharge. We could use our value of \(C_v\) here, but it may be simpler to approach the problem of mass flow rate during choked flow a bit differently, using water-thrust equations we already have.

We first calculate a snapshot of the speed of sound based on the current air temperature inside the bottle using (A2), with \(\lambda=1.34\). Then, using (W2) given previously for water, we substitute this speed of sound for \(v\) and the current air density \(\rho\) inside the bottle to calculate the mass flow \(m'\).

The choked airflow thrust is then the same formula as (W3) given previously, the mass flow rate multiplied by the velocity: $$F_a = m'c \tag{A7}$$ where the speed of sound \(c\) from (A2) is a function of temperature that changes with each time increment.

Thrust from non-choked airflow

As choked airflow thrust is calculated for new time increments, the internal pressure decreases until the choke condition is no longer met. Then we can use similar formulas as (W1), (W2), and (W3) for compressible fluid flow: $$v=C_v \sqrt{2\frac{\lambda}{\lambda-1}\frac{p_e}{\rho_0}} \tag{A8}$$ $$m'=C_c \rho_0 A v \tag{A9}$$ $$F_a=m'v \tag{A10}$$

Remember, here \(p_e\) is the excess (measured) pressure, not absolute pressure \(P_e\). Also we'll use the conservative value \(C_c = 0.9\) as we did in (W2).

The value \(C_v\) for air can be found by setting \(v=c\) and using (A1) to set \(p_e=P_a/\alpha_c − P_a\), then solving for \(C_v\) in (A8). This value yields the speed of sound when the internal pressure is right at the choke threshold. \(C_v\) falls in a range between 0.45 and 0.6 depending on the initial launch conditions of temperature, pressure, and fill volume. \(C_v\) needs to be recalculated at each time step during choked flow to determine when non-choked flow begins.

Final pressure at the end of the time interval

We need to calculate the pressure at the end of the time interval, which would be the pressure to use for calculations at the beginning of the next time interval. To do this, we need to get the new air density for the next interval \(i\): $$\rho_{a~i} = \frac{m_a - m'\Delta t}{V} \tag{A11}$$ then from this relation $$\frac{P_1}{\rho_1^\gamma}=\frac{P_2}{\rho_2^\gamma}=\text{constant}$$ we get the new pressure for the next time step: $$P_{e~i} = P_e \left(\frac{\rho_{a~i}}{\rho_a}\right)^\gamma \tag{A12}$$ where the subscript \(i\) denotes the start of the next time interval, and values without the subscript \(i\) are for the current time step.

To finish up initializing the next time step:

  • Use (A13) to update \(C_v\) in (A8),
  • Use (W6) to calculate the temperature in the bottle from the new pressure,
  • Get a new speed of sound value in (A2) using this new temperature.

Comments

  1. AnonymousSaturday, 01 April, 2023

    Hello, what is C_d in equation A6?

    ReplyDelete
    Replies
    1. AnachronistSunday, 02 April, 2023

      Good catch. \(C_d\) is the discharge coefficient. It looks like I skipped a step there. I ended up not actually using \(C_d\) in my calculations, but determined the mass flow rate using a different approach. I have made a correction to the document. Thanks for pointing this out.

      Delete
    2. AnonymousMonday, 03 April, 2023

      Thank you! Another question, in equation A6 we use $P_d = P_0 - P_v$, here $P_0$ is the atmospheric pressure correct?

      Delete
    3. AnachronistMonday, 03 April, 2023

      You mean A3. Yes \(P_0\) is the atmospheric pressure in Pascals.

      Delete
    4. AnonymousTuesday, 04 April, 2023

      I've been using your equations to write a MATLAB script to simulate a water rocket for a class project. Your website is amazing, it helped me tremendously !! Infinite thanks from Switzerland !!

      Delete
    5. AnachronistTuesday, 04 April, 2023

      I am glad it has been helpful. The spreadsheet I provide in the last chapter of this series may also be helpful, although the formulas may be easier to trace if you download it and load it into Excel instead of using the Google Sheets version.

      Delete

Post a Comment

Popular posts from this blog

Syncing Office 365 Outlook to Google calendar using Power Automate

Image
When my employer switched over from the Google office suite to Microsoft's Office 365, it caused unexpected disruption in my family. I had been sharing my work calendar with my wife (sharing no details, just busy times). It became an integral part of her calendar, to help her coordinate (between her job and my job) who picks up and drops off my son from school, and to let her know when I was free to take a phone call during my work day. Since the switchover to O365, my ability to share a calendar externally was eliminated. Based on documentation I could find, it should be possible, but the option doesn't appear in my employer's implementation of Outlook. I can view my personal Google calendar in my work Outlook calendar, but I can't share it the other way. Microsoft's Power Automate provides a way to solve this problem. It lets you create flows of data, triggered manually or automatically by events, to move data around, not only among the Microsoft products but
Read more

Whose hands are biggest? You may be surprised.

Image
A recent project to create an ergonomic handle for 3D printing led me down a path that introduced me to anthropometric measurements of the human hand, which in turn revealed some interesting facts about hand sizes. I had no idea that hand size is dependent on nationality, but I found a number of research articles on this topic. For those who don't want to read further, the answer to the title question is: Based on the data found, Filipino males have the biggest hands, by a significant margin. In fact, the margin is large enough that this population may be considered an outlier from the general human population with respect to hand size. Even Filipina females have larger hands than most males of other nationalities. On the other end of the spectrum, Vietnamese females and Indian females have the smallest hands. Now you can read on for the whole story, or you can skip down to the chart at the end. Idea for a toy My unintentional investigation into hand sizes started while de
Read more

New approach to screw threads in OpenSCAD

Image
This is an intellectual exercise with a practical end-result. There are better approaches to this problem, but the approach I describe here is interesting. A screw thread from a stack of discs Imagine a stack of 50 coins, stacked so that each coin is slightly offset from the one below it, with the offset moving around in a circle with each new coin. Here is what it looks like for fifty coins 20 mm in diameter and 0.5 mm thick, offset by 2 mm from a center of rotation, and stacked with a rotation offset of 30 degrees per coin ( click to enlarge any image in this article ): Obviously, this is a spiral stack of coins. But it also approximates screw threads. The thread profile looks like a sinewave rather than the trapezoidal profile of machine screw threads. Indeed, others have observed this, and used the freeware parametric CAD program OpenSCAD to make screw threads by extruding an offset circle while twisting it, as in this model on Thingiverse , for example. From a distan
Read more