Launcher with empty 1.3 liter bottle.

In my Soda-bottle rocket simulation post, I posted a buggy version of a water rocket simulation that my students had written, and that I had helped them extend. I asked for help on both this blog and the AP-physics mailing list, since it was clear to me that the simulation was not doing a good job of capturing the phenomena even qualitatively.

Several people pointed me to an old article:

Soda-Bottle Water Rockets

David Kagan, Louis Buchholtz, and Lynda Klein

*The Physics Teacher* 33:150-157, March 1995

I read that article using the university library subscription (there may be copies on the web that aren’t behind the paywall, but I didn’t look for them). I also did a lot of poking around on the web (mainly in Wikipedia). I found a number of bugs in the simulation, the most important of which was that the students had computed the momentum change of the rocket based on the velocity of the exhaust relative to the ground, rather than relative to the rocket.

I reworked the water jet to include a correction for the cross-sectional area of the bottle, though this makes less than a half-a-percent difference for the bottles we were using. I also followed the lead of Kagan et al. and used a fixed exhaust volume for the water, rather than a fixed time step in the simulation. They had used far too large a step (10–30 cc) in their simulation, resulting in over-estimates of the speed and “burnout” height.

The biggest problems with their simulation (which I’ll call the KBK simulation), though, were that they did not simulate air drag or the air jet that occurs once the water has been emptied. These two corrections change the max height from about 136m to about 19m, a much more reasonable number.

The air jet is essential in modeling what happens when an empty bottle is launched. Their simulation suggests that it doesn’t move, but experimentally it goes quite high—maybe half as high as an optimally filled bottle. When I duplicated their simulation, I saw that they still had an extra atmosphere or more of pressure in the bottle at their reported “burnout”. I had the simulation record data at this point (when the bottle emptied of water) and at the real burnout (when the pressure inside dropped to the ambient). If the bottle is overfilled, there may not be enough air pressure to empty it. The point where bottle is full enough that the air pressure is just enough to empty the bottle is an interesting one, as it seems to maximize the burnout height (though not the max height, as the speed at burnout is fairly low then). There also seems to be a bit of glitch in the simulation just above that fill level—I’ve not tracked down the bug yet, but I think it may be from the last packet of water dropping the internal pressure below ambient.

The air-jet simulation that we had written was based on crude modifications to the water-jet code, and it did not seem to work at all. Because the *Matter and Interactions* textbook we’ve been using does not cover anything on fluid flow (not even Bernoulli’s principle for incompressible fluids), I read a number of Wikipedia articles. This is one area of Wikipedia where the articles are far from being standalone—they look like they were extracted from the middle of physics book without providing the links to the required prerequisite material.

I did finally find some good approximations I could use in Orifice plate and Choked flow, which turned out to be better explanations than in Rocket nozzle and more detailed than the article on Bernoulli’s Principle. I ended up modeling both choked flow and subsonic flow, as both seem to be relevant for the soda-bottle rocket.

The bottle doesn’t have an expansion nozzle, but the simulation suggests that adding one might get better performance since the air-jet contributes a lot and we should be able to get at least the first part of the air jet to be supersonic with an appropriate nozzle. I’ve not yet investigated what the specs should be for the nozzle. Rigging up a nozzle that doesn’t interfere with the launch release on the Aquaport launcher might be tough, and I don’t want to build another launcher.

We now get a reasonable result for the optimal fill, at around 350 mL for a 1.3 L bottle (27% full) or 600mL for a 2L bottle (30% full). The empty and full bottles behave reasonably, and the 14g paper and masking tape “rocket” fired with compressed air off a 30cm long launching tube is simulated to go about 4.5 times as high as a 2L soda bottle.

Plot of maximum height vs. the amount of water in the bottle, for both the new simulation and the simulation in the paper by K, Buchholtz, and Klein. Note that adding the drag force makes the maximum height almost independent of fill for a wide range, and adding the air jet gives a reasonable value for the launch of an empty bottle.

I was going to post the data files from the simulations, but WordPress.com does not support any reasonable plain-text formats—a major irritation that I’ve had with them repeatedly. If people want to look at the simulation outputs and have too much trouble re-running them, I can try putting the outputs elsewhere and linking to them.

Burnout speed as a function of water fill. Note that the air jet increases the burnout speed substantially for nearly empty bottles, and that there is a small boost throughout from treating the liftoff from the launcher like a very short-barrel gun. The air drag has very little effect for the short duration of the powered phase of the flight—it mainly affect the ballistic phase.

The KBK simulation is done with finer-grain simulation steps than in the paper, avoiding their 10% overshoot from too large a step.

Height at which the bottle stops expelling water, either because the water is all gone, or because the air pressure has dropped to the ambient pressure. Note that the KBK claim that the this is well fit by a straight line corresponding to a cylinder of water from the nozzle does not match even their simulation, except around the 35% fill region.

I find the cusp near 1.35 L interesting—this is where the rocket just finishes emptying the water when it runs out of pressure. That’s the setting to maximize the height from which water rains down on people, I think. I suspect that with exactly the right settings, this height is the same as the max height, with the water slowly dribbling out as the rocket it traveling almost ballisticly.

We still need to 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, and I was pointed to a number of methods. The most common ones measure height with trigonometry, measuring angles with a couple of theodolites. In a comment on the previous post, TW pointed me to two good tutorials: “a single station method, http://www.apogeerockets.com/education/downloads/newsletter92.pdf Or you can use a two station method.http://www.apogeerockets.com/education/downloads/newsletter93.pdf” I don’t want to build theodolites for just this experiment, though.

The KBK article mentioned using a thread wrapped around a Pasco super pulley ($24) which is similar to the Vernier “Ultra Pulley” ($24). The thread-and-pulley approach seems good for getting finely spaced measurements for the first couple of meters (validating the simulation up to burnout), but not so good for measuring max height.

The other approach is to use video, but that again is best at measuring the first few meters of the launch.

I considered measuring time of flight (which is simpler to measure than height), but the variation in time of flight is fairly small.

The drag model is still a very crude one from first principles, with a guess at the drag coefficient. We should probably do some experiments towing a soda bottle on a string with a force gauge. There is a long enough hill on campus that we should be able to make measurements up to 15 m/s, if we don’t crash trying to stay on the bike path, read the bike speed, and read the force gauge all at once. Still, that 15 m/s is nowhere near the 55 m/s maximum predicted by the simulation. I’ll need to think a bit more about how to measure drag safely.

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 would run much faster without the Unum units in every computation.

# Soda Bottle Rocket Simulation
# Mon Jun 4 08:30:18 PDT 2012
# Kevin Karplus, Abe Karplus, Milo B.
# This 1D simulation of soda-bottle water rockets was started as
# a homework exercise in a home-school high school physics class.
# After the initial programming was done, the program was compared to the
# simulation described in
# Soda-Bottle Water Rockets
# David Kagan, Louis Buchholtz, and Lynda Klein
# The Physics Teacher 33:150-157, March 1995
#
# A few bugs in the program were fixed based on that comparison,
# but this simulation simulates several phenomena that the KBK
# simulation did not---including drag and the jet of air that occurs
# after the water has run out. The air jet is particularly important
# for simulating what happens when an empty bottle is launched.
#
# The drag model is extremely crude and the drag coefficient is a guess.
# Experiments towing empty soda bottles at different speeds would be needed to
# construct and validate a better drag model.
# Tumbling of the bottle (easily observed) has not been incorporated into the
# the drag model.
#
# The model has not yet been compared to experimental data.
### Module Imports
from __future__ import division
from argparse import ArgumentParser
from math import sqrt, pi
# The VPython package is used only for graphs, and can
# be replaced by other graphing packages (or the graphs can be
# removed from the simulator, and external plotting programs used).
import vis
from vis import crayola as color
from vis.graph import gcurve, gdisplay
# The Unum package is used to be sure that units are consistently
# handled throughout the simulation.
# The program could be sped up substantially by using some convention
# about units for each variable, at some risk of introducing
# hard-to-detect errors.
# An ideal compromise would be to have compile-time unit checking and conversion,
# but Python does not offer an easy way to do that.
from unum.units import *
from unum import Unum
Unum.VALUE_FORMAT = "%9.4f" # format for printing units
# derived units
cc = cm*cm*cm
kPa = Unum.unit("kPa", 1000.*Pa, "kilopascal")
# non-standard pressure units (definitions form www.wolframalpha.com)
psi = Unum.unit("psi", 8896442230521./1290320000. *Pa, "pounds per square inch")
inHg = Unum.unit("inHg", 514731./152. * Pa, "inches of mercury")
## Debugging check that pressure units ok
# print "DEBUG: 1 psi=", psi.asUnit(kPa), "1 inch of mercury=", inHg.asUnit(kPa)
def unit_sqrt(squared_unum):
"""sqrt function for the Unum package, retaining appropriate units"""
square_exps = squared_unum._unit
exps = {unit: square_exps[unit]/2 for unit in square_exps}
new_unit = Unum(exps,1)
return sqrt(squared_unum.asNumber(new_unit*new_unit))*new_unit
### Physical constants
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
# We'll ignore any slight changes in water density due to temperature or impurities.
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
air_adiabatic_index = 1.4
# http://en.wikipedia.org/wiki/Heat_capacity_ratio
# based on the approximation that air is a diatomic ideal gas
air_index_ratio = air_adiabatic_index/(air_adiabatic_index-1)
# ratio needed for compressible flow equation
# http://en.wikipedia.org/wiki/Bernoulli%27s_equation
choked_flow_limit = (0.5*(air_adiabatic_index+1))**air_index_ratio
choke_multiplier = (0.5*(air_adiabatic_index+1))**(-(air_adiabatic_index+1)/(air_adiabatic_index-1))
# flow is choked if the ratio of pressures is greater than the choked_flow_limit
# http://en.wikipedia.org/wiki/Choked_flow
def air_density_from_pressure_temp_humidity(p,t,h):
# Compute the density of the air from absolute pressure, absolute temperature, and relative humidity
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 "DEBUG: Standard air density should be 1.2754 kg/m3, is {}".format(standard_air_density)
def sign(x):
"""return the sign of a number"""
return 1 if x>0 else 0 if x==0 else -1
### Class for defining rocket parameters
class rocket_type:
def __init__(self,name,
nozzle_diameter=2.13*cm, # measured for soda bottles
launcher_length=2.2*cm, # measured for AquaPod launcher (2.2cm or 2.4cm to O-ring, 3.6cm full lenth)
drag_coefficient=1.2, # guess at drag coefficient
bottle_diameter=None, # compute if not provided
bottle_volume=None, # compute if not provided
bottle_mass=48*g, # measured for Crystal Geyser bottle
do_air_jet=True):
self.name=name
self.nozzle_diameter=nozzle_diameter
self.nozzle_area= pi/4 * nozzle_diameter*nozzle_diameter
self.launcher_length = launcher_length
self.drag_coefficient=drag_coefficient
# The default value is a guess at drag coefficient for a possibly tumbling bottle
# Use 0.82 for long cylinder, 0.47 for sphere, 0.75 for model rocket
# http://en.wikipedia.org/wiki/Drag_coefficient
if bottle_diameter is None:
bottle_diameter= nozzle_diameter+2*mm # assume cylinder, with 1mm walls (measured for paper rocket)
self.bottle_diameter=bottle_diameter
self.bottle_area = pi/4 * bottle_diameter *bottle_diameter
# Note: the bottle diameters are measured on the outside, to get frontal area for the
# drag calculations, but the same values are used for the inside diameter to get the
# cross-sectional area of the bottle in computing the speed of the water jet.
if bottle_volume is None:
bottle_volume = self.nozzle_area * launcher_length # assume cylinder
self.bottle_volume =bottle_volume
self.bottle_mass =bottle_mass
self.do_air_jet=do_air_jet
self.speed2_scale = 1- (self.nozzle_area/self.bottle_area)**2
# scaling factor for square of speed of exhaust
# see Bernoulli's principle for incompressible flow
# sqrt(speed2_scale) is called the "velocity of approach factor" in
# http://en.wikipedia.org/wiki/Orifice_plate
def __str__(self):
print_items = ["name", "launcher_length", "nozzle_diameter", "nozzle_area",
"bottle_diameter", "bottle_area", "bottle_volume", "bottle_mass",
"drag_coefficient", "do_air_jet"]
return "\n".join(["# {}={}".format(x,self.__dict__[x]) for x in print_items])
def air_drag(self,v,density):
"""a simple drag model, returning 1D force as a function of velocity
and the density of the surrounding air.
This should probably be replaced by an empirically validated drag model.
"""
return -sign(v.asNumber(m/s))*0.5*density*v*v*self.drag_coefficient*self.bottle_area
### Define some standard rockets
# a 1.3 liter bottle for Crystal Geyser sparkling mineral water
rocket_cg = rocket_type("1.3liter",
bottle_volume = 1.302*L, # measured
bottle_diameter = 9.19*cm, # measured (outside diameter)
bottle_mass = 48*g) # measured
# a paper "rocket" made at Maker Faire 2012, launched with a blast of 50psi air
rocket_paper = rocket_type("paper",
launcher_length = 30*cm, #measured
nozzle_diameter = 2.125*cm, # measured
bottle_mass = 14*g, #measured
drag_coefficient = 0.8) # guess at drag coefficient (0.82 for long cylinder, 0.47 for sphere, 0.75 for model rocket)
# a 2 liter soda bottle, measurements from Kagan, Buchholtz, and Klein
rocket_2L = rocket_type("2liter",
bottle_diameter = 11.0*cm,
bottle_volume = 2.08*L,
bottle_mass = 48.4*g)
# a 2 liter soda bottle, measurements from Kagan, Buchholtz, and Klein
# Turn off air-jet simulation, drag, and initial boost from launcher
# in an attempt to match the simulation in Kagan, Buchholtz, and Klein paper
rocket_KBK = rocket_type("KBK",
bottle_diameter = 11.0*cm,
bottle_volume = 2.08*L,
bottle_mass = 48.4*g,
launcher_length=0*cm,
do_air_jet=False,
drag_coefficient=0)
# What rockets are allowed to be specified on the command line?
rockets= [rocket_cg, rocket_paper, rocket_2L, rocket_KBK]
### Command-line Argument Parsing
# The argparse module has been used to provide command-line options for specifying
# various parameters and control options. It also provides a simple help system
# invoked by running
# rocket_3.py --help
parser = ArgumentParser(description='Bottle rocket simulator.')
parser.add_argument('-r', '--realtime', default=False, action="store_true",
help="""if set, use "rate" function to try to make simulation run at real speed
""")
parser.add_argument('-g', '--graph', default=False, action="store_true",
help="""if set, graph various rocket measurements as functions of time.
""")
parser.add_argument('-t', '--temperature', default=18, type=float,
help="""ambient temperature in Celsius
default=%default
""")
parser.add_argument('--humidity', default=0.6, type=float,
help="""relative humidity (between 0 and 1.0)
default=%default
""")
parser.add_argument('--rocket', default="1.3liter", choices=[r.name for r in rockets],
help="""rocket types set several parameters:
1.3liter a Crystal Geyser water bottle
paper a paper "rocket" launched on compressed air
2liter a 2-liter soda bottle
KBK a 2-liter soda bottle, with no drag and no air jet,
to match parameters in paper by Kagan, Buchholtz, and Klein
default=%default
""")
parser.add_argument('-w', '--water', default=50, type=int,
help="""increment for initial water amount in grams.
Used for plotting rocket behavior with different amounts of water initially in bottle.
default=%default
""")
parser.add_argument('--step_time', default=0.001, type=float,
help="""time step in seconds for ballistic phase
default=%default
""")
parser.add_argument('--step_water', default=0.5, type=float,
help="""volume step in cc for water_jet phase
default=%default
Note: Kagan, Buchholtz, and Klein used much too large a step (10cc to 30cc)
""")
parser.add_argument('-p', '--pressure', default=(50*psi).asNumber(bar), type=float, help="gauge pressure in bar")
parser.add_argument('-b', '--barometric', default=29.936, type=float, help="barometric pressure in inches of Mercury")
args = parser.parse_args()
# We need temperature on an absolute scale, not Celsius
init_temperature = (args.temperature+ 273.15)* K
gauge_pressure= args.pressure * bar # pressure read from the bicycle pump
barometric_pressure = (args.barometric*inHg).asUnit(bar) # ambient air pressure
barometric_density = air_density_from_pressure_temp_humidity(barometric_pressure, \
init_temperature,args.humidity)
# density of the ambient air, used for the drag model
# Look for rocket in rockets (it must be there)
for rocket in rockets:
if rocket.name==args.rocket:
break
# Echo the parameters, so that the output is properly documented.
# It might be a good idea to output a program version number and a date as well,
# but that is left as a future enhancement.
print "# Gauge pressure=", gauge_pressure.asUnit(bar), gauge_pressure.asUnit(psi)
print "# barometric pressure=", \
barometric_pressure.asUnit(inHg), \
barometric_pressure.asUnit(bar), \
barometric_pressure.asUnit(psi)
print rocket
# The setup for the graph windows is rather crude,
# with which properties are plotted determined by commenting out code in 3 places.
# The code has not been improved, since we don't use the graphs often.
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, init_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_height/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
rocket_height = 0*m
init_air_pressure = barometric_pressure + gauge_pressure
init_air_density = air_density_from_pressure_temp_humidity(init_air_pressure,init_temperature,1.0)
water_volume=init_water_volume
air_volume = rocket.bottle_volume - water_volume
air_gasconst = init_air_pressure * (air_volume**air_adiabatic_index)
# http://en.wikipedia.org/wiki/Adiabatic_expansion
air_pressure= init_air_pressure
air_mass = air_volume *init_air_density
water_mass = water_volume * water_density
max_height = 0*m
# while on launcher, gains initial speed
starting_energy = (air_pressure - barometric_pressure) * rocket.nozzle_area * rocket.launcher_length
starting_speed = unit_sqrt(2*starting_energy/(rocket.bottle_mass + water_mass + air_mass))
rocket_velocity = starting_speed
rocket_mass = rocket.bottle_mass + water_mass+ air_mass
# water jet phase
while water_volume > 0*L and rocket_height >= 0*m and air_pressure > barometric_pressure:
air_pressure = air_gasconst / (air_volume**air_adiabatic_index)
exhaust_pressure = max(air_pressure - barometric_pressure, 0*Pa)
# from Bernoulli's principle for incompressible flow:
# u^2/2 + barometric_pressure/water_density = constant = v^2/2 + air_pressure/water_density
# Including in the model the internal water velocity v: v = u*nozzle_area/bottle_area
# u^2-v^2 = (1- (nozzle_area/bottle_area)^2)*u^2 = barometric_pressure/water_density - air_pressure/water_density
exhaust_speed = unit_sqrt(2 * exhaust_pressure / (water_density*rocket.speed2_scale))
if exhaust_speed < 1.e-5*m/s:
# for very low exhaust speeds, use a tiny time step, rather than a fixed volume
delta_time = 1.e-5*s
exhaust_volume = rocket.nozzle_area * exhaust_speed * delta_time
if exhaust_volume>water_volume: exhaust_volume=water_volume
else:
# Use steps of constant water leaving rather than fixed time step
# Note: Kagan, Buchholtz, and Klein used much too large a step (10cc to 30cc)
# resulting in overestimates of the speed of the rocket.
exhaust_volume = args.step_water*cc
if exhaust_volume>water_volume: exhaust_volume=water_volume
delta_time = exhaust_volume / (rocket.nozzle_area * exhaust_speed)
water_volume -= exhaust_volume
water_mass = water_volume * water_density
exhaust_mass = exhaust_volume * water_density
air_volume = rocket.bottle_volume - water_volume
rocket_mass = rocket.bottle_mass + water_mass+ air_mass
rocket_velocity += delta_time * rocket.air_drag(rocket_velocity,barometric_density) / rocket_mass
rocket_velocity -= gravity_acceleration * delta_time
rocket_velocity += exhaust_speed*exhaust_mass / rocket_mass
kinetic_energy = 0.5 * rocket_mass * rocket_velocity * rocket_velocity
momentum = rocket_mass * rocket_velocity
rocket_height += rocket_velocity * delta_time
time += delta_time
if use_rate:
vis.rate(1*s/delta_time)
if rocket_height > max_height:
max_height = rocket_height
if use_graph: update_graph(graphs)
empty_time = time
empty_height= rocket_height
empty_speed = rocket_velocity
empty_pressure = air_pressure
# air jet phase
if rocket.do_air_jet:
delta_time = 0.1* args.step_time*s # use smaller time-step while powered
while air_pressure > barometric_pressure and rocket_height >= 0*m:
if use_rate:
vis.rate(1*s/delta_time)
time += delta_time
# Air density is not a constant, because the temperature and pressure both change.
# Perhaps it is best to compute it from the mass and volume directly.
air_density = air_mass/air_volume
pressure_ratio= barometric_pressure/air_pressure
if air_pressure >= choked_flow_limit * barometric_pressure:
# choked flow, barometric pressure doesn't matter
# http://en.wikipedia.org/wiki/Choked_flow
# Because of the taper of the bottle, we'll pretend that the discharge coefficient is close to 1.
exhaust_mass_rate = rocket.nozzle_area* \
unit_sqrt(air_adiabatic_index*air_density*air_pressure*choke_multiplier)
else:
# subsonic flow
# http://en.wikipedia.org/wiki/Orifice_plate (Equation 4)
expansion_factor = sqrt(pressure_ratio**(2/air_adiabatic_index) \
*air_index_ratio \
*(1-pressure_ratio**(1/air_index_ratio)) \
/ (1-pressure_ratio) )
# http://en.wikipedia.org/wiki/Orifice_plate (Equation 3)
exhaust_mass_rate= expansion_factor*rocket.nozzle_area \
* unit_sqrt(2*air_density*(air_pressure-barometric_pressure))
## BUG: I don't believe that the exhaust_speed is right, since
## the density of the air in the jet is not being computed.
## Because we don't have an expansion nozzle,
## let's underestimate speed by assuming density is the same as in rocket.
exhaust_speed = exhaust_mass_rate/(air_density*rocket.nozzle_area)
exhaust_mass= exhaust_mass_rate*delta_time
exhaust_momentum = exhaust_mass*exhaust_speed
rocket_mass = rocket.bottle_mass + air_mass + water_mass
rocket_velocity += delta_time * rocket.air_drag(rocket_velocity,barometric_density) / rocket_mass
rocket_velocity -= gravity_acceleration * delta_time
rocket_velocity += exhaust_momentum / rocket_mass
kinetic_energy = 0.5 * rocket_mass * rocket_velocity * rocket_velocity
momentum = rocket_mass * rocket_velocity
# I don't think that this air_pressure update is quite right,
# as it assumes no change in temperature of the air.
air_pressure *= (air_mass-exhaust_mass)/air_mass
air_mass -= exhaust_mass
rocket_height += rocket_velocity * delta_time
if rocket_height > max_height:
max_height = rocket_height
if use_graph: update_graph(graphs)
burnout_speed=rocket_velocity
burnout_time = time
burnout_height = rocket_height
# ballistic phase
rocket_mass = rocket.bottle_mass + water_mass+ air_mass
delta_time = args.step_time*s # use normal time step
while rocket_height >= 0*m:
if use_rate:
vis.rate(1*s/delta_time)
time += delta_time
rocket_velocity += delta_time * rocket.air_drag(rocket_velocity,barometric_density) / rocket_mass
rocket_velocity -= gravity_acceleration * delta_time
rocket_height += rocket_velocity * delta_time
kinetic_energy = 0.5 * rocket_mass * rocket_velocity * rocket_velocity
momentum = rocket_mass * rocket_velocity
if rocket_height > max_height:
max_height = rocket_height
if use_graph: update_graph(graphs)
return dict(
starting_time=0*s,
water=init_water_volume,
starting_height=rocket.launcher_length,
starting_speed=starting_speed,
empty_time=empty_time,
empty_height=empty_height,
empty_speed=empty_speed,
empty_pressure=empty_pressure,
burnout_time=burnout_time,
burnout_height=burnout_height,
burnout_speed=burnout_speed,
max_height=max_height,
flight_time=time)
### Simulation Caller
# which values from the simulation do we want to print, and with what units
fields= [("water",L),
("empty_time",s),
("empty_speed",m/s),
("empty_height", m),
("empty_pressure", bar),
("burnout_speed", m/s),
("burnout_height", m),
("flight_time", s),
("max_height", m)]
print( "\t".join(["{:15}".format(f[0]) for f in fields]))
print( "\t".join([" {:11}".format(f[1].strUnit()) for f in fields]))
# We don't use water in the paper rocket, so only run the loop with the 0 value
volume_limit = 1 if rocket.name=="paper" else int(rocket.bottle_volume.asNumber(cc)+1)
for water in range(0,volume_limit,args.water):
values = single_run(gauge_pressure, water*cc, use_rate=args.realtime, use_graph=args.graph)
print("\t".join([ "{:9.4f}".format(values[f[0]].asNumber(f[1])) for f in fields]))

### Like this:

Like Loading...