# Gas station without pumps

## 2012 May 28

### Soda-bottle rocket simulation

In my Soda-bottle rockets post, I showed a picture of the very simple soda-bottle rocket launcher I’ve used for the past 10 years (with pointers to the plans for making it), and in Soda-bottle rockets used, I mentioned the possibility of doing a physics lab using a \$26 triggered launcher to get more repeatable launches. I bought the launcher, along with other toys (like the Astro-blaster and a syringe for pressure-volume experiments), so I had the students design a lab for it.

The trigger is a simple design: a spring-loaded lever that hooks over the rim of the bottle and is pulled off by a tug on a long cord.  The force needed to release the trigger is pretty high, so the launcher has to be staked to the ground (or a few heavy bricks placed in front of the legs, to keep the launcher from sliding). Close up of the launcher release. The stake can be seen just to the right of the bottle. Note that we had to add a paperclip to the launch lever as a jump ring, since the hook provided for the cord would not reach the hole in the lever intended for it. The bottle here has no pressure in it, so did not leave the launcher, even though the lever is pulled back to release position.

First, I insisted on them coming up with equations to model the rocket and writing a simulation, so that they could optimize the amount of water to put in the rocket. Then they did some preliminary tests to see whether the simulation was anywhere near correct (with 0g of water and with 400g of water).

The first simulation, which was done almost entirely by the students, predicted much too low a max height, because it only included two sources for thrust: the initial slide off the launcher (using a simple energy model and constant pressure), and conservation of momentum as the water was expelled.

The model was then improved by adding an air jet after the water was expelled, using essentially the same equations as for the water jet, but with lower density and a slightly different update for the pressure on each time step.  By this point I was heavily involved in the modeling, because I wanted a model that we could use. The model with both the water jet and the air jet predicted too high a max height, so one more term was added: a drag force proportional to the velocity squared.

All the parameters except the drag coefficient (and frontal area, if the bottle tumbles) are directly measured, except for the local gravitational field for which we used a Wolfram Alpha widget.  We didn’t measure barometric pressure, temperature, and humidity ourselves, but relied on reports on the web from a local weather station a couple blocks away.

We now have a simulation that has roughly the right values at the extremes (empty bottle and full bottle), and that peaks at an intermediate value, but I’m still not convinced that the model works properly.  We consistently get that the amount of water for max height is around 1/9th the volume of the bottle, but I seem to remember getting the best results in the past with amounts around 1/3rd.  One thing that the model does get right is that it predicts a much higher flight for a 14g paper and masking tape “rocket” fired with compressed air off a 30cm long launching tube than for a 42g bottle off the triggered launcher.

So now we need to do two things:

1. Do a lot of launches with various amounts of water and measure the flights.  My main concern here is coming up with an accurate measurement.  Lots of people do water and compressed-air rocket labs of various sorts—what do they measure and how?  I considered measuring time of flight (which is simpler to measure than height), but the variation in time of flight is fairly small.  For height, we are probably limited to a trigonometric measurement, but I have my doubts about the accuracy of measurement with a protractor and plumb-bob.
2. Look for phenomena we haven’t properly included in the simulation.  For example, we have not simulated any friction for the water or air leaving the bottle, nor checked to see whether the air flow is subsonic (I think it is).  If the simulation predicts super-sonic air flow, then we need to change the model, since the lack of any expansion nozzle means that there is definitely a limitation to Mach 1 for the air jet.

I’d like help from the physics (or rocketry) experts, both on the simulation the students wrote (with help from me) and on how to do the empirical measurements to test the fit of the simulation.

Here is the latest draft of the simulation code, which uses the VPython library for graphing and the Unum library for handling units.  The code written by the students did not include explicit units, but did everything in standard SI units.  I added the explicit units to check to see if we had any conversion errors (there were none).  The code would run much faster without the Unum units in every computation.

[UPDATE 2012 June 14: please see Soda-bottle rocket simulation: take 2 for a better simulation.]

```#################################
# Bottle Rocket Simulation      #
# 24-28 May 2012                #
#################################

### Module Imports

from __future__ import division
from argparse import ArgumentParser
from math import sqrt, pi
import vis
from vis import crayola as color
from vis.graph import gcurve, gdisplay

from unum.units import *
from unum import Unum

# derived units

cc = cm*cm*cm
liter = 1000.*cc
kPa = 1000.*Pa

### Constant Initialization

g_grav = 9.80665* m/s/s	# standard gravity
gravity_acceleration = 9.7995 * m/s/s # from Wolfram Alpha's gravitational field widget
# for the gravitational field in Santa Cruz, CA
# http://www.wolframalpha.com/widgets/view.jsp?id=d34e8683df527e3555153d979bcda9cf
water_density = 1 * g/cc

specific_gas_constant_for_dry_air = 287.058 * J/(kg * K)
# http://en.wikipedia.org/wiki/Density_of_air
specific_gas_constant_for_water_vapor = 461.5 * J/(kg * K)
# http://www.engineeringtoolbox.com/density-air-d_680.html

inHg = 3386.389 * Pa	# from http://en.wikipedia.org/wiki/Inch_of_mercury

def air_density_from_pressure_temp_humidity(p,t,h):
return p /(specific_gas_constant_for_dry_air* t)  * (1+h)/(1+h*specific_gas_constant_for_water_vapor/specific_gas_constant_for_dry_air)

# Debugging check that air density ok
# standard_air_density = air_density_from_pressure_temp_humidity(100*kPa, 273.15*K, 0.0)
# print "Standard air density should be 1.2754 kg/m3, is {}".format(standard_air_density)

### Command-line Argument Parsing

parser = ArgumentParser(description='Bottle rocket simulator.')
parser.add_argument('--temperature', default=18, type=float, help="temperature in Celsius")
parser.add_argument('--humidity', default=0.6, type=float, help="relative humidity (between 0 and 1.0)")
parser.add_argument('-w', '--water_step', default=50, type=int, help="increment for water amount in grams")
parser.add_argument('-t', '--timestep', default=0.001, type=float, help="time step in seconds")
parser.add_argument('-p', '--pressure', default=3, type=float, help="gauge pressure in bar")
parser.add_argument('-b', '--barometric', default=30.23, type=float, help="barometric pressure in inches of Mercury")
args = parser.parse_args()

temperature = (args.temperature+ 273.15)* K

gauge_pressure=args.pressure * bar
barometric_pressure = args.barometric*inHg
barometric_density = air_density_from_pressure_temp_humidity(barometric_pressure,temperature,args.humidity)

print "Gauge pressure=", gauge_pressure.asUnit(Pa)
print "barometric pressure=", barometric_pressure.asUnit(Pa)

def air_density_from_pressure(p):
# Air density for 100% humidity at default temperature
return air_density_from_pressure_temp_humidity(p,temperature, 1.00)

# Rocket and launcher parameters
rocket=args.rocket

if rocket=="1.3liter":
launcher_length = 2.2*cm #, 0.022, 0.024, 0.036 # measured (from O-ring, from O-ring, from base of launcher)
nozzle_diameter = 2.13*cm # measured
bottle_diameter = 9.19*cm # measured
bottle_volume = 1.302*liter # measured

bottle_mass = 48*g  # measured
drag_coefficient = 1.0	# guess at drag coefficient for bottle 	(0.82 for long cylinder, 0.47 for sphere, 0.75 for model rocket)
# http://en.wikipedia.org/wiki/Drag_coefficient
elif rocket=="paper":
bottle_mass = 14*g	#measured
launcher_length = 30*cm #measured
nozzle_diameter = 2.125*cm # measured
bottle_diameter = nozzle_diameter+2*mm #measured
bottle_volume = pi/4 *nozzle_diameter*nozzle_diameter*launcher_length # cylinder
drag_coefficient = 0.8	# guess at drag coefficient (0.82 for long cylinder, 0.47 for sphere, 0.75 for model rocket)

else:
raise ValueError("Error: unknown rocket type: {}".format(rocket))

frontal_area = pi* bottle_diameter*bottle_diameter/4

def sign(x):
return 1 if x>0 else 0 if x==0 else -1

def air_drag(v):
return -sign(v.asNumber(m/s))*0.5*barometric_density*v*v*drag_coefficient*frontal_area

def velocity_sqrt(v_squared):
return sqrt(v_squared.asNumber(m*m/s/s))*m/s

if args.graph:
position_window =gdisplay(title='Position (m)', y=0, foreground=color.black, background=color.white)
velocity_window =gdisplay(title='Velocity (m/s)', y=100, foreground=color.black, background=color.white)
pressure_window =gdisplay(title='Air pressure (kPa)', y=200, foreground=color.black, background=color.white)
#   water_window =gdisplay(title='Water mass (g)', y=300, foreground=color.black, background=color.white)
#   energy_window =gdisplay(title='Kinetic energy (J)', y=400, foreground=color.black, background=color.white)
#   momentum_window =gdisplay(title='|Momentum| (kg m/s)', y=500, foreground=color.black, background=color.white)

### Simulation Function

def single_run(gauge_pressure, water_volume, use_rate=True, use_graph=True):
### Graph Setup

if args.graph:
graphs = dict()
graphs["position"] = gcurve(gdisplay=position_window, color=color.red)
graphs["velocity"] = gcurve(gdisplay=velocity_window, color=color.green)
graphs["pressure"] = gcurve(gdisplay=pressure_window, color=color.green)
#       graphs["water"] = gcurve(gdisplay=water_window, color=color.blue)
#       graphs["energy"] = gcurve(gdisplay=energy_window, color=color.black)
#       graphs["momentum"] = gcurve(gdisplay=momentum_window, color=color.black)

def update_graph(graphs):
graphs["position"].plot(pos=(time/s, rocket_position/m))
graphs["velocity"].plot(pos=(time/s, rocket_velocity/(m/s)))
graphs["pressure"].plot(pos=(time/s, air_pressure/kPa))
#   graphs["water"].plot(pos=(time/s, water_mass/g))
#   graphs["energy"].plot(pos=(time/s, kinetic_energy/J))
#   graphs["momentum"].plot(pos=(time/s, abs(momentum.asNumber(kg*m/s))))

# initialization
time = 0*s
delta_time = 0.1* args.timestep*s # use smaller time-step while powered
rocket_position = 0*m
air_pressure = barometric_pressure + gauge_pressure
air_volume = bottle_volume - water_volume
air_gasconst = air_pressure * air_volume
air_density = air_density_from_pressure(air_pressure)
air_mass = air_volume *air_density

water_mass = water_volume * water_density
nozzle_area = pi * nozzle_diameter * nozzle_diameter / 4
max_height = 0*m

# while on launcher, gains initial speed
starting_energy = (air_pressure - barometric_pressure) * nozzle_area * launcher_length
starting_velocity = velocity_sqrt(2*starting_energy/(bottle_mass + water_mass + air_mass))

rocket_velocity = starting_velocity

# water jet phase
while water_volume > 0*liter and rocket_position >= 0*m and air_pressure > barometric_pressure:
if use_rate:
vis.rate(1*s/delta_time)
time += delta_time
air_pressure = air_gasconst / air_volume
exhaust_pressure = max(air_pressure - barometric_pressure, 0*Pa)
exhaust_force_mag = exhaust_pressure * nozzle_area
exhaust_impulse_rocket_mag = exhaust_force_mag * delta_time
exhaust_velocity_rocket_mag = velocity_sqrt(exhaust_impulse_rocket_mag / (nozzle_area * water_density * delta_time))
exhaust_volume = nozzle_area * exhaust_velocity_rocket_mag * delta_time
if exhaust_volume > water_volume:
exhaust_volume = water_volume
water_volume -= exhaust_volume
water_mass = water_volume * water_density
exhaust_mass = exhaust_volume * water_density
air_volume = bottle_volume - water_volume
rocket_mass = bottle_mass + water_mass+ air_mass

exhaust_impulse_earth_mag = exhaust_mass * abs(rocket_velocity - exhaust_velocity_rocket_mag)
rocket_velocity += delta_time * air_drag(rocket_velocity) / rocket_mass
rocket_velocity -= gravity_acceleration * delta_time
rocket_velocity += exhaust_impulse_earth_mag / rocket_mass
kinetic_energy = 0.5 * rocket_mass * rocket_velocity * rocket_velocity
momentum =  rocket_mass * rocket_velocity

rocket_position += rocket_velocity * delta_time
if rocket_position > max_height:
max_height = rocket_position
if use_graph: update_graph(graphs)

empty_after = time
empty_velocity = rocket_velocity

air_volume=bottle_volume
# air jet phase
while air_pressure > barometric_pressure and rocket_position >= 0*m:
if use_rate:
vis.rate(1*s/delta_time)
time += delta_time
exhaust_pressure = max(air_pressure - barometric_pressure, 0*Pa)
exhaust_force_mag = exhaust_pressure * nozzle_area
exhaust_impulse_rocket_mag = exhaust_force_mag * delta_time
air_density = air_density_from_pressure(air_pressure)
exhaust_velocity_rocket_mag = velocity_sqrt(exhaust_impulse_rocket_mag / (nozzle_area * air_density * delta_time))
exhaust_volume = nozzle_area * exhaust_velocity_rocket_mag * delta_time
exhaust_mass = exhaust_volume * air_density
air_mass = air_volume *air_density
air_pressure *= (air_mass-exhaust_mass)/air_mass
rocket_mass = bottle_mass + air_mass

#       print "DEBUG:", time, "air_mass=",air_mass, "exhaust_mass=", exhaust_mass, "exhaust_v=", exhaust_velocity_rocket_mag

exhaust_impulse_earth_mag = exhaust_mass * abs(rocket_velocity - exhaust_velocity_rocket_mag)
rocket_velocity += delta_time * air_drag(rocket_velocity) / rocket_mass
rocket_velocity -= gravity_acceleration * delta_time
rocket_velocity += exhaust_impulse_earth_mag / rocket_mass
kinetic_energy = 0.5 * rocket_mass * rocket_velocity * rocket_velocity
momentum =  rocket_mass * rocket_velocity

rocket_position += rocket_velocity * delta_time
if rocket_position > max_height:
max_height = rocket_position
if use_graph: update_graph(graphs)

peak_velocity=rocket_velocity

# ballistic phase
delta_time = args.timestep*s # use normal time step
while rocket_position >= 0*m:
if use_rate:
vis.rate(1*s/delta_time)
time += delta_time

rocket_velocity += delta_time * air_drag(rocket_velocity) /rocket_mass
rocket_velocity -= gravity_acceleration * delta_time
rocket_position += rocket_velocity * delta_time
kinetic_energy = 0.5 * rocket_mass * rocket_velocity * rocket_velocity
momentum =  rocket_mass * rocket_velocity

if rocket_position > max_height:
max_height = rocket_position
if use_graph: update_graph(graphs)

return starting_velocity, empty_after, empty_velocity, time, peak_velocity, max_height

### Simulation Caller

print("water, init_velocity, empty_after, empty_velocity, peak_velocity, flight_time, max_height")
max_volume = 1 if rocket=="paper" else int(bottle_volume.asNumber(cc))

for water in range(0,max_volume,args.water_step):
starting_velocity, empty_after, empty_velocity, flight_time, peak_velocity, max_height = \
single_run(gauge_pressure, water*cc, use_rate=args.realtime, use_graph=args.graph)
print("{}\t{}\t{}\t{}\t{}\t{}\t{}".format(water*cc, starting_velocity, empty_after, empty_velocity, peak_velocity, flight_time, max_height))
```

1. How high are your rockets going? You could launch them next to a building that you can premeasure every 1/2 meter or so and video tape the launch. This would also give you speed parameters like using the picket fence with photogates.

Comment by Joanna Evans MacDonald — 2012 May 29 @ 14:12

• I had thought of trying to use the video camera to measure speed. There aren’t many tall buildings near here, so I doubt that we could measure height that way (the bottles go well higher than our house), but we could probably measure speed for the first 5m, if we could mark up the side of the house. The simulation indicates that the rocket is only powered for the first 1–2m, and the peak speed is no more than about 25 m/s, so we should be able to get 3 or 4 frames of the bottle as a ballistic projectile at 29.97 fps. I’m a bit worried about parallax effects though, as the camera is going to be near the ground, and the bottle will be some distance from the wall.

Comment by gasstationwithoutpumps — 2012 May 29 @ 14:33

2. I can’t know for sure whether it will help, but if it were mine, the first place I would look to upgrade this model is line 102: the implementation of air_drag(v). I would expect air drag to be a major factor, and the ideal drag equation is such a rough estimator that, in my limited exposure to projects involving aerodynamics, I’ve only seen it used as the basis for an interpolation function between as many measured drag values as the project could collect. I don’t think any professional uses it to calculate air drag from first principles, even with a measured drag coefficient. The drag coeff ends up being a function of v, the relevant “transverse area” isn’t the cross-section if there’s any sort of concavity in the direction of motion (such as you have in the bottom of a soda bottle), the exponent on v ends up not being quite 2 due to non-laminar air flow, etc. Everyone just uses wind tunnels, models hanging under aircraft, etc., measures, and interpolates.

If you haven’t already done so, you might try to think of clever ways to measure the pressure on the bottle when pushed bottom-first into wind of various speeds. If you could get someone to take you out for a few drives (ideally in different temp/humidity conditions) while you sit in the back seat and make lots of measurements at lots of different speeds, you could curve fit it, empirically validate (or invalidate) the drag equation, and get a reliable implementation of air_drag(v) (for a specific rocket, that is.) Make sure to hold it far enough away from the car to avoid the car’s own laminar flow unless you have an accurate anemometer that can measure the air speed directly.

Comment by Glen — 2012 May 29 @ 19:30

• Yes, I’m pretty sure that the drag force is awful. I’ll discuss possible measurements with the students.

I was also worried about friction for the exhaust. Is that likely to be a significant problem?

Comment by gasstationwithoutpumps — 2012 May 29 @ 19:51

• Good point about the exhaust. I don’t know enough about fluid dynamics to answer, but the flow rate is such a significant factor that I would want to measure it, too. So, how to measure the rate of flow as a function of pressure?

I don’t see an easy way to measure this using air pressure (you may have one), but the pressure ultimately becomes water pressure at the nozzle anyway, so maybe you could use water pressure as a substitute for air pressure. Cut the nozzle end off a soda bottle and attach it to some sort of tall tube that’s open at the top. The pressure at the nozzle will depend solely on the depth of the water column above it, not on friction of any sort at the nozzle or on the sides of the tube (since water is essentially incompressible) or on the shape or capacity of the container above it. That would give you a very accurate measure of nozzle pressure.

Open the nozzle, let it run for a few (carefully measured) seconds into a bucket, quickly close it, and see how much water went through the nozzle in that time. Since pressure is linear with depth, the average pressure could be determined by the midpoint between the depth of the column at start and the depth at end. You could then gather a lot of good data points showing rate of flow past the nozzle as a function of pressure.

You’d need a water column 33 ft above the nozzle to get a full atmosphere of overpressure, but you might be able to do that with a 50 ft section of PVC pipe running down a 45 degree incline down at the beach. Just a thought.

I assume these rockets are pumped up to more than 1 atm of overpressure, but I’m guessing that this test would produce clean enough data that you could confidently extrapolate to the full pressure of your bottle at start.

I think it’s great that you’re teaching physics this way. I MUST do something like this when my kids are old enough.

Comment by Glen — 2012 May 30 @ 14:48

• We’re launching at 3bar, which would be 30.6m of water. I don’t think I want to try building a tower and plumbing that high! It is close to the municipal water pressure, though. Unfortunately, the plumbing in my house is narrower than the soda-bottle neck, so the flow rate would be limited by some other part of the plumbing.

We could measure the water flow rate roughly by filming a captive soda bottle, and seeing how quickly the water emptied.

Measuring the flow rate of the exhaust air is tougher, but the simulation suggests that the air is more important than the water.

Comment by gasstationwithoutpumps — 2012 May 30 @ 17:02

3. There is an article that will help you in The Physics Teacher magazine. It is called “Soda-Water Bottle Rockets” by David Kagan, Louis Buckholtz, and Lynda Klein. It is here, http://tpt.aapt.org/resource/1/phteah/v33/i3/p150_s1?isAuthorized=no

Comment by TW — 2012 June 2 @ 08:22

• Thanks, an AP Physics teacher on the AP-physics mailing list pointed me to the Kagan, Buckholz, and Klein article yesterday, and I downloaded it. There are two things there of note: 1) they use Bernoulli’s equation to compute the exhaust velocity of the water, which is probably better than what we did in our simulation, 2) they did not model the air jet, so their model would imply no height for the no-water case, which is clearly very wrong.

They used the superpulley and string method, which requires having a superpulley (about \$24 from Pasco or Vernier) plus a photogate big enough to work with pulley (another \$45, unless I jury-rig something). The Vernier “Ultra Pulley” has 15cm circumference of the inner groove and 1.16E-6 kg m^2 moment of inertia, according to their documentation. I don’t know what sort of bearings it has. Pasco tells us about their bearings and the mass of the pulley, but not the moment of inertia or the circumference. They do say that it mounts on a 1/4″x20 thread (the same as a camera tripod), so it should be easy to mount the pulley (Vernier fails to specify the threads on theirs).

I’ve been thinking about the possibility of making an pulley using Lego parts that would work with the tiny photogate I have, but I’m not sure it is feasible. The largest Lego pulley we have has only about an 11cm circumference, it has 2 spokes instead of 10, and Lego shafts in holes are lousy bearings with a very high drag. Keeping the electronics dry could be a problem also. I understand that winding several meters of thread around the pulley can be a problem also. Thy only used about 2m of thread, so that only the power phase of the launch was measured.

I’ve also been trying to figure out whether we could jury rig something with a couple of old CDs and the magnetic bearing from an old mousetrap car—although magnetic bearings are very low friction, they can’t take much shear force, and I suspect that the jerk on the string would just pull the whole thing apart.

I’ve now downloaded and looked at the apogeerockets.com newsletters also—they provide a nice tutorial on measuring height. Unfortunately, I don’t have 2 theodolites, and I’m not sure I want to make them. (The less accurate single-station method in

Comment by gasstationwithoutpumps — 2012 June 2 @ 09:28

4. After reading “Soda-Water Bottle Rockets” by David Kagan, Louis Buckholtz, and Lynda Klein, I see a couple of problems with our simulation.

1. We tried accelerating just the water leaving the nozzle, but it can’t move without the rest of the water moving with respect to the bottle. This is an error which will overestimate the speed of the exhaust. The article uses Bernoulli’s principle, and ends up with
exhaust_velocity=sqrt(2 *exhaust_pressure / (water_density *(1 – (nozzle_diameter/bottle_diameter)**4)))
while the formula we used was
exhaust_velocity=sqrt(exhaust_pressure / water_density)

The diameter ratio to the 4th power is about 0.003 and so makes almost no difference, but there is the factor of 2 that we need to explain.

2. The other error is that we assumed that air_pressure*air_volume is constant, which would be the case if the air was at constant temperature, but the expansion is so rapid that there isn’t time for much heat flow, and the air cools rapidly. One can see that the temperature drops below the dew point, because the bottle has a little cloud inside when it lands.

The article uses the assumption that the expansion is adiabatic (no heat flow) and gets
air_pressure = initial_air_pressure * (initial_volume/volume)**1.4
where the 1.4 comes from assuming that air is a diatomic gas. This is also a big difference.

Comment by gasstationwithoutpumps — 2012 June 2 @ 15:54

• I found where our factor of 2 error came from. We had assumed that the impulse was equal to the ejected mass times its velocity, but forgot that the mass had to be accelerated up to the exhaust velocity, so the impulse was mass * velocity/2. (more correctly (internal velocity+external velocity)/2, but we’ve been assuming that the velocity of the water inside the bottle is small enough to be negligible, which is pretty much true for these diameters.

I’ve not yet looked at Bernoulli’s principle for compressible flow, which would apply to the air jet.

Comment by gasstationwithoutpumps — 2012 June 2 @ 17:16

5. […] the pressure down to around 20–40 psi and the volume is small.  I have no idea what plastic the commercial soda-bottle rocket launcher uses (maybe ABS, which is not as brittle as PVC), and they have a pressure release around 60psi (4 […]

Pingback by NASA releases reasons for removing paper rocket activity « Gas station without pumps — 2012 October 14 @ 10:38

6. […] Soda-bottle rocket simulation […]

Pingback by 2012 in review « Gas station without pumps — 2012 December 31 @ 11:18

7. […] that I believe I bought from Arbor Scientific (which can be seen in the superpulley post and the water-rocket simulation posts 1 and 2).  Nonetheless, several families did take instructions for making the simple PVC launchers, […]

Pingback by Soda bottle rockets used again | Gas station without pumps — 2016 May 12 @ 20:48

This site uses Akismet to reduce spam. Learn how your comment data is processed.