Gas station without pumps

2012 June 4

Soda-bottle rocket simulation: take 2

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]))

10 Comments »

  1. […] would like us to continue on the rocket experiment, now that the simulation seems to be working.  Launching next to a tall building and making a video recording may be our best bet for measuring […]

    Pingback by Physics homework chapter 12 « Gas station without pumps — 2012 June 5 @ 09:33 | Reply

  2. […] Soda bottle rocket simulation take 2, I mentioned that the authors of  one paper used “a thread wrapped around a Pasco super […]

    Pingback by Homemade super pulley « Gas station without pumps — 2012 June 10 @ 23:34 | Reply

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

    Pingback by Soda-bottle rocket simulation « Gas station without pumps — 2012 June 14 @ 14:06 | Reply

  4. […] while the rocket is on the launcher—it is modeled more like a gun than like a rocket. (But my soda-bottle rocket simulator can model these paper bullets also.)  It would probably best to have a shorter barrel for doing […]

    Pingback by Nerf gun prototype 1 « Gas station without pumps — 2012 July 3 @ 21:12 | Reply

  5. I was going to download the latest draft of the python code. The article says “Here is the latest draft of the simulation code…” but I didn’t see a link to the code.

    Comment by Aaron Titus — 2012 September 12 @ 20:26 | Reply

    • You should see a “show source” button on the line just below the “Here is the latest draft …” paragraph. If not, try a different browser, as I see it with Firefox, Chrome, and Safari. I don’t have Internet Exploder, so I can’t test with that, nor with any of the rarely-used Linux browsers.

      Comment by gasstationwithoutpumps — 2012 September 12 @ 21:58 | Reply

      • Thank you. I am using Chrome. The “show source” link was light gray and I just didn’t see it. Thanks for the tip. I see it now

        Comment by Aaron Titus — 2012 September 13 @ 04:14 | Reply

  6. […] Soda-bottle rocket simulation: take 2 […]

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

  7. […] Soda-bottle rocket simulation: take 2 […]

    Pingback by Blogoversary 3 | Gas station without pumps — 2013 June 1 @ 20:00 | Reply

  8. […] 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, with […]

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


RSS feed for comments on this post. TrackBack URI

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out /  Change )

Google photo

You are commenting using your Google account. Log Out /  Change )

Twitter picture

You are commenting using your Twitter account. Log Out /  Change )

Facebook photo

You are commenting using your Facebook account. Log Out /  Change )

Connecting to %s

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

%d bloggers like this: