Gas station without pumps

2016 May 12

Soda bottle rockets used again

Filed under: Uncategorized — gasstationwithoutpumps @ 20:48
Tags: , , ,

I had another chance on Tuesday this week to play with soda-bottle water rockets, which I have not done since my son was taking home-school physics and we wrote the timing program for measuring the ascent of the rockets that later turned into PteroDAQ to go along with the homemade Lego “superpulley”.

My wife volunteered me to help out at the Spring Hill School’s Family Science Night—a tradition I started 9 years ago (here are my notes from the 2 years I ran it).  She had thought I could set up the “Dr. K’s Applied Electronics” display I’d used at the Mini Maker Faire, but I didn’t have the time it would take to set up the display (over an hour, to do it right).  Instead, we agreed to revive the soda-bottle water rockets, which are always a hit with elementary-school kids.

I raced home from running the pulse-monitor lab for the Applied Electronics course, loaded up my bike trailer with empty mineral water and soda bottles, a couple of floor pumps, and some rocket launchers.  My wife also brought a few copies of the instructions for making the PVC launchers, though not the Spanish-language version, which we first created in 2001, when my son was in kindergarten (a post I wrote in 2011 talked a bit about the activities I did with the Kindergarten kids to explain how rockets worked).

The rockets were a big hit at Family Science Night last night, though the kids preferred the commercial launcher 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, with the intent of doing it as a fun project this summer.

It was fun to do the science outreach, though I only managed to talk with one family about how rockets worked, and it gave me a good excuse for not doing any grading that night (after 3 days of doing not much besides grading, teaching, and supervising the students in lab, a night supervising kids sending up soda-bottle water rockets was a welcome relief).

2012 October 14

NASA releases reasons for removing paper rocket activity

Filed under: Robotics — gasstationwithoutpumps @ 10:37
Tags: , , , , , ,

In a previous post, I had groused about NASA removing plans for a common compressed-air launcher from their educational web site for “safety reasons” without explaining what the hazards were.  (There are several possible hazards, including shooting the rockets at people, firing more massive projectiles, and exploding PVC pressure containers.)  I asked for a copy of the engineering report, and (10 weeks later) I’ve finally gotten a reply.

They have a FAQ page which they have sent out as PDF files to people who inquired, but they don’t seem to have put it on their website anywhere (a strange oversight—I would have put it up on the web site first and sent people links to it, rather than sending out PDFs).  They did not send me the engineering report I requested, insisting that it was only available through a Freedom of Information Act request.  Again, a very strange, anti-education, anti-safety approach—I would have put the engineering report on the web, since there was clearly public interest and a need for the information to be disseminated.  I get the impression that NASA is being run by lawyers and politicians, whose first instinct is to make everything require expensive intervention by lawyers and whose second instinct is to prevent the spread of information if at all possible.  There may have been a time when NASA was run by scientists and engineers, but they certainly aren’t now.

In any event the FAQ is quite clear on what hazard they were talking about:

NASA does not recommend use of PVC pipe with compressed gasses, including air. Under pressure, PVC can shatter or explode under pressure or from an external force.

At this time NASA has no plans to redesign an activity using PVC pipe to construct a launch system that utilizes compressed air. NASA will assess other materials and designs and may release a new high-powered paper rocket launcher at a future date.

NASA completed an inquiry into this activity and determined that the launcher, or design equivalents, should not be used. NASA does not recommend use of PVC pipe with compressed gasses, including air. PVC can shatter or explode under pressure or from an external force. NASA recommends that individuals and organizations should immediately discontinue use of these launchers.

Q8: Can I obtain a copy of the engineering report referenced in the discontinuation notice?
You may request a copy of the report under the Freedom of Information Act. (http://www.hq.nasa.gov/office/pao/FOIA/agency)

Exploding PVC is one of the possible hazards I had considered, and perhaps the most dangerous one, since not all people who build with PVC are aware of its brittleness and that it produces very sharp shards when it shatters. The brittleness and sharp are why toy swords are not made from PVC pipes. (The Society for Creative Anachronism experimented with PVC for their swords in the 1980s, when rattan was getting expensive, and decided that PVC was far too dangerous.)

There is a standard workaround that reduces (though does not completely eliminate) the hazard: wrapping all pressurized PVC components with a few layers of strapping tape.  The strapping tape does not prevent the PVC from shattering, but may contain the shrapnel or slow it down enough to reduce the danger zone. The strapping tape only has to hold the pieces together for long enough for the air to escape—the energy of the pressurized air is dissipated in tearing the tape rather than in propelling the shrapnel.

Even more common is keeping the pressure fairly low.  This is not as big a win as one might think, since even fairly modest pressures (like 40psi) can still store significant energy in a large pressure vessel, and PVC can shatter at low pressures if it is struck or if it has gotten UV damage from being left in the sun.  Keeping the vessel small and the pressure low is probably safe enough, so I have no concerns about the soda-bottle rocket launcher plans that I’ve published in the past—the friction-fit launching keeps 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 bar), so again is probably safe enough.

The one project we’ve done that potentially hazardous is the foam-dart launcher, which is very similar to the NASA design for the paper-rocket launcher (indeed, we have used it as a paper-rocket launcher).  We were pressurizing to 120psi (8 bar), and we had not wrapped the pressure vessels with strapping tape.  That project is on hold currently, not for safety reasons, but because the students in the robotics club got too busy to continue the robotics.  If we pick it up again, we’ll wrap the pressurized PVC pipes with a few layers of strapping tape, then add a layer of colored duct tape (for UV protection of the strapping tape and to make it look nice).  I think that would reduce the hazard to an acceptable level for experimenters, though not for a commercial product.

2012 June 12

Improved super pulley code

In Homemade super pulley, I described the “super pulley” that my son built and that I wrote some Arduino code for.

After our preliminary experiments using the super pulley to record the initial parts of soda-bottle rocket flights, my son was dissatisfied with the minimal code that I had written.  It was difficult to do the cut and paste from the Arduino serial monitor to save the data, and entering the necessary metadata was bit difficult with the bright sun washing out the laptop screen.

He decided it would be worth his time to improve the code.  He made two major changes to the code: first, he implemented a Python program on the laptop to capture the Arduino output and save it (with metadata) to a file.  Second, he rewrote the Arduino code to be interrupt driven and to have a queue of output, so that we could record far more than 400 points, while still being able to handle bursts of very rapid data.

These programs got him into three new coding concepts for him: circular buffers for queues, interrupts, and multi-threading.  He researched and chose these techniques, then implemented the program without my assistance, although I did go over it with him at the end, editing for style and clarity of documentation.  In fact, I tried to talk him out of using multi-threading—I favor lighter-weight techniques like the Unix select() or poll() functions, but portability won for him.  There is one place where portability was lost, though—he couldn’t find documentation on how to look for USB devices from Python on a Windows machine, and just used glob('/dev/tty.usb*') to look for the Arduino, which should work for Macs, Linux, and other Unix clones. On Windows, it looks like you have to open one of the COM ports, but I’m not sure how you figure out which one.

The only non-standard Python library used is PySerial, for talking to the USB port.

He decided to use curses rather than PyGUI or a plain terminal window, because he wanted to do some display of the number of samples as they were read, but did not want the hassle of setting up a full GUI.  (He has used both curses and PyGUI previously in science fair projects, and used PyGUI for the underwater ROV control, where curses would not have been adequate.)

Here is the revised Arduino code:


// Smart Pulley Rocket Timer
// By Abe Karplus and Kevin Karplus
//
// To use, connect the output from the photogate of a smart pulley to
//     digital pin 2 on the Arduino.
// On initialization, the Arduino will print "Arduino ready." to serial.
// For each tick of the photogate, the Arduino prints three tab-separated
//     values on a line. The first is the time since the start, in
//     microseconds. The second is the amount of string unwound, in meters.
//     The third is the instantaneous speed of unwinding, in meters per second.
// The ticks are buffered so that they may be recorded faster than the
//     serial can output them. The buffer holds 400 ticks, maximum.
// Serial baudrate is 115200.
// The time and distance may be reset by sending the character 'R'.
//
// License: CC-BY-NC http://creativecommons.org/licenses/by-nc/3.0/

#define BUFFER_SIZE (400) // size of tick buffer (limited by RAM)
#define TICKS_PER_REV (6) // number of holes in smart pulley
#define DIAMETER (0.044) // diameter of pulley sheave in meters

const float meter_per_tick = 3.14159265358979 * DIAMETER  / TICKS_PER_REV;
    // meters of string unwound per tick

// queue of times implemented as circular buffer
volatile unsigned long times[BUFFER_SIZE];
volatile unsigned int read_head  = 0; // position in buffer for removal
volatile unsigned int write_head = 0; // position in buffer for insertion
volatile bool wrote_last = 0;   // what was most recent operation on queue?
                                // 0 for read, 1 for write

// Next item to be read is at times[read_head].
// Next empty space to write into is at times[write_head].
// Number of entries in buffer is (write_head - read_head) mod BUFFER_SIZE.
//      Buffer is full if read_head==write_head and wrote_last.
//      Buffer is empty if  read_head==write_head and not wrote_last.

unsigned long count;        // how many ticks have been seen since reset
unsigned long old_time;     // if count>0, what was the time of the previous tick?
unsigned long first_time;   // what was the time for the first tick since reset?
// All times are in microseconds.

// print_time() prints one output line,
//      removing the corresponding time from the queue
//      and updating count, old_time, (and first_time when needed).
// It does nothing if the queue is empty.
void print_time(void)
{
    noInterrupts();     // don't add to the queue while reading from it
    if (read_head == write_head && !wrote_last)
    {   // buffer empty
        interrupts();
        return;
    }
    unsigned long datum = times[read_head++];   // time taken from queue.
    if (read_head>=BUFFER_SIZE)
    {   read_head = 0; // wrap around
    }
    wrote_last = 0;
    interrupts();       // finished getting the data, re-enable adding

    if (count==0)
    {   // On first call, do no output, just initialize the global times.
        old_time = first_time = datum;
        count++;
        return;
    }
    Serial.print(datum - first_time); // time since start (usec)
    Serial.print("\t");
    Serial.print(count++ * meter_per_tick, 3); // distance unwound (m)
    Serial.print("\t");
    Serial.println(meter_per_tick*1e6 / (datum - old_time), 3); // speed (m/s)
    old_time = datum;
}

// get_time() is the interrupt handler
// It adds the time returned by micros() to the queue,
// unless the queue is full, in which case nothing is done.
void get_time(void)
{
    // record time immediately to avoid
    // fluctuation due to variable amount of code executed
    unsigned long now = micros();

    if (read_head == write_head && wrote_last)
    {   // buffer full
        return;
    }
    times[write_head++] = now;
    if (write_head>=BUFFER_SIZE)
    {   // wrap_around
        write_head=0;
    }
    wrote_last = 1;
}

// do_reset() empties the queue and sents the count to 0
void do_reset(void)
{
    noInterrupts();     // don't take interrupts in the middle of resetting
    read_head = write_head = 0;
    wrote_last = 0;
    count = 0;
    // let the user know that the Arduino is ready
    Serial.print(F("\nArduino ready.\n"));
    interrupts();
}

// interrupt handler for data becoming available on the serial input
// Reset if an R is received
void serialEvent(void)
{
    if (Serial.read() == 'R')
    {
        do_reset();
    }
}

void setup(void)
{
    Serial.begin(115200);       // use the fastest serial connection available
    do_reset();
    // now start looking for ticks on pin 2
    attachInterrupt(0, get_time, FALLING); // pin 2 interrupt
}

void loop(void)
{   // keep trying to print out what is in the queue.
    print_time();
    delay(1);
}

And here is the Python code:

###
 #  Smart Pulley Rocket Timer
 #  By Abe Karplus
 #  12 June 2012
 #
 #  Interfaces with the Arduino program rockt_timer.ino for water rocket
 #   timing using a smart pulley. To use, upload rocket_timer.ino to the
 #   Arduino and run `python rocket_timer.py` on the command line.
 #  This program uses the PySerial library from pyserial.sourceforge.net
 #   in addition to the standard libraries glob, curses, threading, and time.
 #  To find the Arduino, this program searches in the /dev directory;
 #   this will need to be changed for Windows systems.
 #  If supported by the terminal (such as xterm-256color), uses colors.
 #
 #  License: CC-BY-NC http://creativecommons.org/licenses/by-nc/3.0/
###

from serial import Serial, SerialException
from glob import glob
import curses
from threading import Thread, Event
from time import sleep

def cursor_vis(n):
    """Change the visibility of the cursor (same as curs_set),
    but ignore errors.
    """
    try:
        curses.curs_set(n)
    except curses.error:
        pass

def send_reset(ard):
    """Send the Arduino a reset command character."""
    # discard any Arduino output lines that are waiting to be read
    ard.flushInput()

    ard.write('R')

def wait_for_ready(ard):
    """Discard input until a line with 'Arduino ready.' appears.
    Needed to avoid leftover trash from Arduino bootloader.
    """
    for line in ard:
        if 'Arduino ready.' in line:
            return

def prepare_arduino():
    """Create a serial connection to an Arduino
    (assuming it is the unique device with /dev/tty.usb*)
    and return the Serial file-like object.
    """
    # find Arduino port
    ports = []
    while len(ports) != 1:
        ports = glob('/dev/tty.usb*')
    port = ports[0]

    # create a Serial connection
    while True:
        try:
            arduino = Serial(port=port, baudrate=115200)
        except SerialException as e:
            # occurs upon physically connecting USB
            # so not really a problem.
            continue
        else:
            break

    wait_for_ready(arduino)
    return arduino

def add_readings(scr, ard, rds, cond):
    """Infinite loop run in a daemon thread.
    When Event cond is set,
        Read data from the arduino Serial object ard.
        Store the data in list rds.
        Display len(rds) on the curses window scr.
    """
    while True:
        cond.wait()
        line = next(ard).strip()
        if 'Arduino ready.' in line or not line:
            # Arduino has been reset
            rds[:] = []    # clear list in-place
            scr.move(3, 22)
            scr.clrtoeol()
        else:
            rds.append(line)

        # display running count of reads
        scr.addstr(3, 22, str(len(rds)), curses.color_pair(4))
        scr.refresh()

def prepare_screen(scr):
    """Initialize curses window scr.
    Set up color pairs and display initial message.
    """
    try:
        curses.init_pair(1, 0, 15) # Black on White
        curses.init_pair(2, 0, 10) # Black on Green
        curses.init_pair(3, 0, 11) # Black on Yellow
        curses.init_pair(4, 0,  9) # Black on Salmon
        curses.init_pair(5, 0, 14) # Black on Cyan
    except curses.error:
        for n in range(1, 6):
            curses.init_pair(n, curses.COLOR_BLACK, curses.COLOR_WHITE)

    # use black on white text by default
    scr.bkgdset(ord(' '), curses.color_pair(1))

    cursor_vis(0)   # hide cursor
    scr.clear()

    # setup permanent header
    scr.addstr(1, 1, 'Welcome to the Rocket Timer smart pulley interface.', curses.A_BOLD)

    # setup message line
    scr.addstr(3, 1, 'Waiting for the Arduino to connect.', curses.color_pair(3))

    scr.refresh()

def save_readings(scr, rds, d=dict(volume='1.3', pressure='3', water='0')):
    """Request metadata and file name, and save rds to file.
    Clear rds to become an empty list.
    d is a dictionary of default values for file metadata.
    """
    # clear message and prompts
    scr.move(3, 1)
    scr.clrtobot()

    # setup message
    scr.addstr(3, 1, 'Saving reads.', curses.color_pair(5))
    scr.refresh()
    cursor_vis(1)    # make cursor visible
    curses.echo()    # let user see what they type

    scr.addstr(4, 1, 'Save to file: ')
    fname = scr.getstr().strip()    # request file name
    scr.move(4, 1)
    scr.clrtoeol()

    scr.addstr(4, 1, 'Bottle volume in liters (default {}): '.format(d['volume']))
    # use default if blank, or changes default
    d['volume'] = scr.getstr().strip() or d['volume']
    scr.move(4, 1)
    scr.clrtoeol()

    scr.addstr(4, 1, 'Gauge pressure in bar (default {}): '.format(d['pressure']))
    # use default if blank, or changes default
    d['pressure'] = scr.getstr().strip() or d['pressure']
    scr.move(4, 1)
    scr.clrtoeol()

    scr.addstr(4, 1, 'Water mass in grams (default {}): '.format(d['water']))
    # use default if blank, or changes default
    d['water'] = scr.getstr().strip() or d['water']
    scr.move(4, 1)
    scr.clrtoeol()

    # return curses to display only (no character echo and hidden cursor)
    curses.noecho()
    cursor_vis(0)

    header = '\n'.join(l.strip() for l in """\
    # {fname}
    # Bottle volume: {d[volume]} liters
    # Gauge pressure: {d[pressure]} bar
    # Water mass: {d[water]} grams
    # Each line represents one tick of the photogate.
    # `time` is the time in microseconds since the beginning.
    # `distance` is the amount of string unwound, in meters.
    # `speed` is the instantaneous speed in meters per second of the unwinding.
    time\tdistance\tspeed
    N\tN\tN
    """.split('\n')).format(**locals())
    with open(fname, 'w') as f:
        f.write(header)
        for ln in rds:
            f.write(ln)
            f.write('\n')
    rds[:] = []    # clear list in-place

def do_interface(scr, ard, rds, cond):
    """Parse user interface commands, using a state machine.
    Arguments:
        scr is curses window
        ard is Arduino Serial object
        rds is list of reads
        cond is Event to control the add_readings daemon
    """
    # interpretations of keys typed in response to prompts:
    transitions_from_ready = {' ':'recording', 'q':'quitting', 'r':'reset'}
    transitions_from_save_query = {'y':'save readings', 'n':'reset'}

    state = 'reset'
    while True:
        scr.move(2, 0)
        scr.clrtobot()
        if state == 'reset':
            send_reset(ard)
            sleep(0.01)
            scr.addstr(3, 1, 'Waiting for Arduino to reset.', curses.color_pair(3))
            scr.refresh()
            sleep(0.5)
            wait_for_ready(ard)
            state = 'ready to record'
        elif state == 'ready to record':
            scr.addstr(3, 1, 'Ready to record.', curses.color_pair(2))
            scr.addstr(5, 2, 'Space to record, q to quit, r to reset:')
            scr.refresh()
            key = chr(scr.getch()).lower()
            state = transitions_from_ready.get(key,state)
        elif state == 'recording':
            scr.addstr(3, 1, 'Recording. Readings: 0', curses.color_pair(4))
            scr.addstr(5, 2, 'Space to finish.')
            scr.refresh()
            send_reset(ard)
            wait_for_ready(ard)
            cond.set()     # turn on recording in daemon thread
            while chr(scr.getch()) != ' ':
                pass
            cond.clear()    # turn off recording
            state = 'save query'
        elif state == 'save query':
            scr.addstr(3, 1, 'Save {} readings? (y/n)'.format(len(rds)), curses.color_pair(5))
            scr.refresh()
            key = chr(scr.getch()).lower()
            state = transitions_from_save_query.get(key,state)
        elif state == 'save readings':
            save_readings(scr, rds)
            state = 'reset'
        elif state == 'quitting':
            scr.addstr(3, 1, 'Quitting.', curses.color_pair(4))
            scr.refresh()
            sleep(0.5)
            return
        else:
            raise ValueError('Invalid state {}'.format(state))

def main(scr):
    """Set up screen and threads and start user interface.
    """
    prepare_screen(scr)

    arduino = prepare_arduino()

    readings = []

    cond = Event()
    cond.clear()
    reader_thread = Thread(target=add_readings, name='reader', args=(scr, arduino, readings, cond))
    reader_thread.daemon = True    # will not stop main from exiting
    reader_thread.start()

    do_interface(scr, arduino, readings, cond)

curses.wrapper(main)
print 'Done.'

The code here is released under Creative Commons license CC-BY-NC:

2012 June 10

Homemade super pulley

In Soda bottle rocket simulation take 2, I mentioned that the authors of  one paper used “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.”

Since I’m too cheap to buy a super pulley, especially since the vendors don’t even say how wide they are, so I can’t tell whether they’d work with the 1cm-wide photogates that I have, I decided to make my own.  Actually, since I’m also lazy, I had my son design and make the superpulley out of Lego.

Top view of the Lego super pulley, showing the Arduino, the photogate with the timing pulley, and the large hub used as the pulley sheave.

Closeup of the timing pulley and photogate. The timing pulley has 6 equally spaced holes, so the photogate provides a pulse every 60° of rotation. This is not as fine a resolution as the Vernier and Pasco pulleys, which provide a square wave, giving 10 or 20 timing events per revolution (36° or 18°). My sheave has a 4.4cm diameter, for a tick every 23.0mm, while the Vernier pulley could give a tick every 7.5mm.

To protect the photogate, there is a lid across the top that was removed for the previous photos. Here is a side view with the lid in place.

Because we will be unwinding a thread wrapped around the pulley, it is very convenient that the Lego wheel hub has a hole in it, which we can use for starting the winding of the thread.

The pulley was made last night, and we tested it by wrapping a thread around it and yanking as hard as we could.  We couldn’t get it to spin nearly as fast as the rocket was predicted to go, but it spun much more easily than I had expected. To get velocity, we took one-sixth of the circumference and divided by the time between ticks—this was smoother than I expected, so we did not try to do any computational smoothing. We also tried hooking up a Lego motor directly to the pulley to try to get a higher speed. The fastest speed we managed this way was still not as fast as we expected the rocket to go. (Trying a belt drive for more speed didn’t help—the Lego motor stalled.)

The first version of the Arduino software was very simple, and just reported each time interval (and equivalent velocity) back to the laptop over the USB line.  I noticed, though that the maximum baud rate of the Arduino serial port (115,200 baud) would only let us communicate reliably if the ticks were at least a millisecond apart.  At 23mm per tick, this would put an upper limit of 23 m/s on the measurements, but the predicted speed of the rocket was higher than that.  So I rewrote the code to record the times in an array, and print (with fancier formatting) the whole array at the end of the measurement.  To avoid having a complicated communication protocol, I have the program print the array and reset it to empty whenever there is a tick longer than a second.

One problem with the Arduino is that it has a really tiny RAM—only 2k bytes—and the time stamps in microseconds take 4 bytes each.  Because of some overhead, it is not possible to store 500 time stamps, but 400 was no problem.  This allows up to 9.2m of motion before the array fills up (at which point everything is printed).  Even if a very long string were not a serious drag problem for a rocket, we would not be able to use more than 9.2m, unless we traded off some resolution (say recording only alternate ticks or recording time with msec resolution instead of µsec).

We were not sure whether this pulley design would work with the rockets, so we did preliminary tests this afternoon, launching air-only rockets.  (That way we did not have to lug large bottles of water over to the field, and we did not have to figure out how to waterproof the photogate, Arduino, and laptop yet.)

The pump, launcher, pulley, and Arduino setup. We carried the concrete block the ¾ mile so that the pulley would not move when the string was pulled on.  Next time, we’ll take a bike trailer instead of walking—the block was a pain to carry.

The laptop was a bit hard to see in the sun, and the USB cable is definitely too short for launching rockets that spray water everywhere. Luckily I have USB-to-CAT5 converters, so I can use a CAT5 Ethernet cable as an extension cord, as we did for the tether on the underwater ROV.

A closer view of the connection between the empty water bottle used as a rocket and the pulley. Although the angle of this photo does not make it clear, the launcher and pulley were aligned so that the thread pulled straight out from the pulley.

For the first several launches, we ended up only saving the second burst of 400 data points, not the first one, so only saw the pulley spinning down, not the initial acceleration of the pulley.

When spinning down, the pulley shows a fairly constant deceleration, corresponding to -0.50 m/s^2. The sharp braking at the ends of the runs was from stopping the pulley by hand.  The jitter in the velocity seems to have a periodicity of about 6 ticks, and probably corresponds to the shaft wobbling in the “bearings”.  It is, after all, a cross-shaped Lego shaft in a loose round hole, which can hardly be expected to be a high-quality bearing.

We had another problem with the first few launches: the sewing thread that we were using snapped where it was tied onto the rocket. The first fix we tried was switching to buttonhole twist, the strongest thread we had brought with us. This was not sufficient. Two other fixes were tried. One was to reduce the launch pressure from 3bar (300,000 Pa) to 2bar and even 1bar, and the other was to use some Scotch tape as a strain relief. The combination of 2bar and the Scotch tape worked, but I’d like to try again without the tape, as I’m not sure it helped.  We’ll also try getting some stronger thread (perhaps some surf line from a fishing store), so that we can try higher pressures.

The string was tied around the neck of the bottle then taped to the side of the bottle. On launch, the string either pulled the tape off the bottle or cut through the tape. The hope was that this strain relief would allow the tension on the string to increase more gradually as the pulley spun up.

After clearing up the confusion about which data we were to cut and paste from the Arduino serial monitor and getting the string to stop snapping, we were finally able to record a few good launches.  For the next set of launches, we’ll use a Python interface on the laptop that my son wrote to capture the data directly from the Arduino and put it into files.

With the bottle pressurized to only 1 bar (gauge pressure, obviously, since that is about one atmosphere of pressure), the thread did not completely unwind from the pulley. Really only the initial acceleration pulse is tracking the rocket—after that the spinning of the pulley is unwinding the thread faster than the rocket is pulling it.  For one of the launches, the Lego supports for the shaft had been tweaked a bit without our noticing, and the drag was much higher.  Peak acceleration seems to be around 190 m/s^2, or about 19g.

With 2 bar of pressure, we got the full length of thread (about 6m) unspooled from the pulley, without breaking. The acceleration gets as high as 1270 m/s^2 (about 130g).
The spindown of the pulley on the three runs indicates different amounts of drag on the pulley. Note that the noise in the velocity measurement drops considerably once the thread detaches from the pulley.

I’m very confused by the data in the 2-bar launches. Why does the blue curve have so much lower initial acceleration than the rest? Was that the Scotch tape strain relief working? Why does the red curve have a curving velocity plot (decreasing deceleration) during spindown?

[UPDATE: 2012 June 11.  The predicted burnout height and velocity from the 2-bar simulation are 0.32m and 17.2 m/s.  The 1-bar simulation suggests 0.14m and 10.0m/s. The low-pressure launches peaked at about 7–8m/s and the high-pressure ones around 18 m/s. the burnout height looks to be about 0.3m in most of the launch measurements, but this may be more a function of how long it takes the pulley to get up to speed than how fast the rocket is really going. Given the low precision of our pressure gauge on the pump and the noise in the velocity measurements, this is pretty good agreement.]

I’m actually rather amazed that we managed to record accelerations of 130g without snapping the thread. Because the thread can’t really convey any information to the pulley once the rocket starts decelerating, I think we can get away with a shorter thread, at least until we try bottles that filled with so much water that they reach burnout just at the apex of their flight, which may be as high as 10m, based on our simulation.

The burnout takes longer with water in the bottle, so the pulley can get us more data on the flight, but the powered phase is still only about 1/10th of a second, and we won’t be able to observe much more than acceleration to peak speed and peak speed.

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]))
Next Page »

%d bloggers like this: