Math & Physics Reference
Every equation behind the paper rocket simulator, from pneumatic launch to landing
Contents
- Pneumatic Launch Physics
- Rocket Construction & Mass
- Aerodynamic Drag
- Trajectory & Ballistics
- Stability (Barrowman Equations)
- Spin & Gyroscopic Effects
- Altitude Measurement
- Drag Analysis from Flight Timing
These equations have been cross-validated against the NASA Rockets Educator Guide (EG-2020-11-46-MSFC). See the validation report on the Physics page.
1. Pneumatic Launch Physics
A pneumatic launcher stores energy as compressed air in a chamber. When the valve opens, the gas expands rapidly behind the rocket, accelerating it down the barrel. Understanding this process lets us predict muzzle velocity from chamber pressure, volume, and barrel dimensions.
Adiabatic Expansion
The release happens so quickly that very little heat transfers to or from the gas. We model it as an adiabatic (no heat exchange) process. For an ideal gas expanding adiabatically:
- P1 = initial absolute pressure in the chamber (Pa)
- V1 = initial chamber volume (m3)
- P2 = pressure after expansion
- V2 = expanded volume (chamber + barrel volume behind rocket)
- γ = ratio of specific heats = 1.4 for air (diatomic gas)
As the rocket travels down the barrel, V2 grows and P2 drops. The pressure behind the rocket at any position x along the barrel is:
- Abarrel = cross-sectional area of the barrel (m2)
- x = distance the rocket has traveled down the barrel (m)
Work Done by Expanding Gas
The net force on the rocket at position x equals the pressure difference times the barrel area: F(x) = [P(x) − Patm] · Abarrel. The total work done as the rocket traverses the barrel of length L is the sum of this force over the barrel distance.
Algebra approach — divide the barrel into N small steps of width Δx = L / N and sum:
Exact Integral for Work
Replacing the sum with a continuous integral over the barrel length:
Substituting the adiabatic pressure expression:
Let u = Vc + Abx, so du = Ab dx. The pressure integral becomes:
Evaluating:
With γ = 1.4, the exponent 1 − γ = −0.4. The total net work also subtracts the atmospheric back-pressure contribution: Watm = Patm · Ab · L, giving W = Wgas − Watm.
Conversion to Kinetic Energy
By the work-energy theorem, the net work on the rocket equals its kinetic energy at the muzzle (neglecting friction in the barrel):
Solving for muzzle velocity:
Valve Flow: Cv Rating
In reality, the valve restricts airflow. The valve's Cv (flow coefficient) determines how fast air can pass through. A larger Cv means less restriction. For a fast-opening valve with Cv much larger than needed, the ideal adiabatic model is a good approximation. For smaller valves, flow choking limits the effective pressure behind the rocket.
The mass flow rate through the valve under choked conditions (when upstream pressure > ~1.9× downstream) is:
Effect of Chamber Volume, Pressure, and Barrel Length
- Increasing chamber volume — the gas has more stored energy and maintains higher pressure longer as it expands, increasing muzzle velocity.
- Increasing pressure — directly increases force on the rocket. Velocity scales roughly as √P for a given chamber/barrel geometry.
- Increasing barrel length — gives the gas more distance to do work. But the pressure drops as the volume behind the rocket grows. Beyond a certain length, the gas pressure drops to atmospheric and extra barrel length adds friction without benefit.
2. Rocket Construction & Mass
Accurate mass estimation is critical. Even a gram matters for a paper rocket that typically weighs 5-20 grams. Each component is computed from material properties and geometry.
Paper Tube (Body)
The body tube is formed by wrapping paper around a mandrel (the launch tube). The mass depends on paper density, tube dimensions, and number of wraps:
- ρpaper = areal density of paper, typically 75-90 g/m2 for copy paper
- Ltube = length of the tube (m)
- wsheet = width of each paper sheet used for wrapping (m)
- nlayers = number of paper layers (wraps)
Equivalently, using the circumference of the mandrel to find the paper width per wrap:
Tape Mass
Tape adds significant mass. Common types:
- Packing tape: ~35 g/m2 (clear, thick)
- Electrical tape (e-tape): ~120 g/m2 (heavy, stretchy)
- Masking tape: ~55 g/m2
For a single wrap of tape around the circumference: Lapplied = π · d.
Nose Cone Mass
A paper nose cone is a section of a circle folded into a cone. The lateral (side) surface area of the cone determines paper usage:
- r = base radius of the cone (matches tube outer radius)
- s = slant height = √(r2 + hcone2)
- hcone = height of the nose cone
Fin Mass with Attachment Tape
Fins are cut from card stock or folded paper. Each fin has a trapezoidal or triangular shape with area Afin. The attachment tape formula accounts for both the tape covering the fin surface and the fillet strips along the root chord:
- 2 · Afin = tape on both sides of the fin
- 4 · sfin · croot = fillet strips: each fin has two fillets (one each side), each strip is sfin (fin span) by croot (root chord), times 2 sides × 2 fillets per fin... simplified as 4 · span · root
Total fin system mass for n fins:
Total Mass Budget
3. Aerodynamic Drag
Once the rocket leaves the barrel, drag is the dominant force (besides gravity) acting on it. Understanding drag is essential for predicting flight performance.
The Drag Equation
- ρair = air density, approximately 1.225 kg/m3 at sea level, 15 °C
- v = velocity of the rocket relative to the air (m/s)
- Cd = drag coefficient (dimensionless)
- Aref = reference area, usually the frontal cross-section = πr2
The key insight: drag scales with velocity squared. Doubling speed quadruples drag force. This is why paper rockets decelerate so rapidly after launch.
Drag Coefficient Breakdown
The total Cd is a sum of contributions from different rocket components:
- Nose cone (Cd ~ 0.1-0.3): A pointed cone has lower drag; blunt noses are higher. A 3:1 length-to-diameter ogive is nearly optimal.
- Body friction (Cd ~ 0.2-0.5): Skin friction along the body tube. Depends on Reynolds number and surface roughness. Tape seams and bumps increase this.
- Fins (Cd ~ 0.1-0.3): Both skin friction on the fin surfaces and form drag from fin thickness.
- Base drag (Cd ~ 0.1-0.2): The low-pressure wake behind the flat base of the rocket. A boat-tail reduces this but paper rockets rarely have one.
- Interference & roughness: Additional drag from tape edges, glue bumps, fin-body junctions, and imperfections.
Reynolds Number
The Reynolds number determines whether the boundary layer is laminar (smooth, less drag) or turbulent (rough, more drag):
- L = characteristic length (body length for body drag, chord for fin drag)
- μ = dynamic viscosity of air ≈ 1.81 × 10−5 Pa·s
For a typical paper rocket at 30 m/s with a 25 cm body length, Re ≈ 500,000. This is in the transitional regime, meaning parts of the boundary layer may be laminar and parts turbulent. The simulator uses empirical correlations for skin friction that account for this transition.
Drag Integral in Trajectory
Because drag depends on velocity, and velocity changes continuously, the energy lost to drag over a flight segment requires integration (see Section 4 for the full trajectory treatment).
Energy Lost to Drag
The work done by drag over a path from position s1 to s2:
Since v(s) is not generally known in closed form (it depends on drag, which depends on v), this integral must typically be evaluated numerically — exactly what the simulator does with its step-by-step Euler method.
4. Trajectory & Ballistics
After the rocket leaves the barrel at muzzle velocity, its path is governed by gravity and aerodynamic drag. Unlike the simple parabolic trajectories of introductory physics, drag makes this problem nonlinear — there is generally no closed-form solution.
2D Equations of Motion
We decompose the motion into horizontal (x) and vertical (y) components. The velocity vector has magnitude v at angle θ from horizontal:
The drag force always opposes the velocity vector. During ascent with velocity components (vx, vy):
Differential Equations of Motion
Newton's second law in component form, with drag opposing velocity:
During descent, the drag opposes downward motion (so drag acts upward):
More compactly, the sign handling is automatic if we write:
These coupled, nonlinear ODEs have no closed-form solution. We must solve them numerically.
Euler Method (Numerical Integration)
The simulator uses the Euler method: step forward in small time increments Δt, updating position and velocity at each step.
At each time step:
- Compute the total speed: v = √(vx2 + vy2)
- Compute the drag acceleration magnitude: adrag = ½ ρ v2 Cd A / m
- Apply drag opposite to velocity components and subtract gravity from vertical component
- Update velocities and positions
Maximum Height Estimation
For a vertical launch (simplest case), max height occurs when vy = 0. Without drag:
With drag, height is always less than this. A useful approximation for a vertical shot defines the ballistic coefficient:
Analytical Max Height with Drag (Vertical Launch)
For vertical motion with quadratic drag during ascent:
Define k = ρ Cd A / (2m). Using the substitution v dv = a dy (chain rule: dv/dt = (dv/dy)(dy/dt) = v dv/dy):
This is separable. Rearranging and integrating from v0 to 0 and y = 0 to hmax:
The left integral evaluates using the substitution u = g + kv2:
Substituting back k = ρ Cd A / (2m):
Notice that when drag is negligible (k → 0), ln(1 + x) ≈ x for small x, and this reduces to v02 / (2g), matching the vacuum case.
Time of Flight
Without drag, the total flight time for a launch at angle θ from height y = 0:
With drag, the ascent takes less time (drag + gravity decelerate it) and the descent also takes less time (drag limits fall speed). The simulator computes flight time directly from the numerical integration by detecting when y ≤ 0.
Range Calculation
Without drag, range is maximized at 45 degrees:
With drag, the optimal angle shifts lower (typically 30-40 degrees for paper rockets). The simulator finds the exact range by tracking horizontal position during the Euler integration.
Launch Angle Optimization
For maximum range with drag, the optimal launch angle depends on the ratio of muzzle velocity to terminal velocity. When drag is significant (v0 much greater than vt), the optimal angle drops well below 45 degrees. The simulator sweeps angles in 1-degree increments to find the maximum.
5. Stability (Barrowman Equations)
A rocket is stable if it naturally corrects itself when disturbed (e.g., by a gust of wind). Stability requires that the center of pressure (CP) be behind the center of gravity (CG). The Barrowman equations (1966) provide a straightforward method to locate the CP for a subsonic, axially symmetric rocket at small angles of attack.
Center of Gravity (CG)
CG is the mass-weighted average position of all components, measured from the nose tip:
- mi = mass of each component (nose, body, fins, clay, tape, etc.)
- xi = distance from nose tip to each component's centroid
Component centroid positions (measured from nose tip):
- Nose cone: xnose = 2/3 · Lnose (for a conical nose, the centroid is 2/3 of the way back)
- Body tube: xbody = Lnose + Lbody / 2
- Fins: xfins = distance from nose to fin centroid (depends on fin placement)
- Nose weight (clay): xclay ≈ Lnose / 3 (packed into the nose tip)
Center of Pressure (CP) via Barrowman Method
Each component generates a normal force coefficient CN when the rocket flies at a small angle of attack α. The CP is the force-weighted average of each component's CP location:
Nose Cone Contribution
For a conical nose:
For an ogive nose shape, the normal force coefficient is still 2, but the CP shifts slightly forward to about 0.466 · Lnose.
Fin Contribution (Full Barrowman)
For n identical trapezoidal fins, the normal force coefficient and CP location are:
- n = number of fins (typically 3 or 4)
- s = fin semi-span (height from body surface to tip)
- d = body diameter (reference)
- cr = root chord (fin length at body)
- ct = tip chord (fin length at tip, 0 for triangular fins)
- r = body radius = d / 2
The CP of the fins, measured from the fin leading edge at the root:
- xLE = distance from nose tip to the fin root leading edge
- msweep = sweep distance (horizontal offset of the tip leading edge from the root leading edge)
Stability Margin
- SM > 1.0 caliber: Stable. The rocket will weathercock into the wind and fly straight.
- SM = 1.0-2.0 calibers: Ideal range for paper rockets.
- SM > 3.0 calibers: Overstable. The rocket will weathercock aggressively and may not fly in the intended direction in wind.
- SM < 1.0 caliber: Marginally stable to unstable. The rocket may tumble.
- SM < 0: Unstable. CG is behind CP. The rocket will tumble and loop.
6. Spin & Gyroscopic Effects
Canting (angling) the fins slightly induces a roll that spin-stabilizes the rocket, much like a rifle bullet. Spin stabilization supplements fin-based aerodynamic stability and helps maintain a consistent flight path.
Fin Cant Angle and Spin Rate
When fins are mounted at a cant angle δ from the airflow, they generate a side force component that creates a rolling torque. The equilibrium spin rate occurs when the aerodynamic roll moment equals the roll damping moment.
- ω = angular velocity of roll (rad/s)
- CL,fin = lift coefficient of each fin (typically 2πα for small angles)
- δ = fin cant angle (typically 1-5 degrees)
- v = airspeed
- reff = effective moment arm (roughly the average radial position of fin area, ≈ r + s/2)
Converting to revolutions per second:
Gyroscopic Stability
A spinning body resists changes to its spin axis. The angular momentum vector points along the spin axis, and any disturbance torque causes precession rather than tumbling.
- I = moment of inertia about the roll axis (kg·m2)
- ω = spin rate (rad/s)
For a thin-walled cylinder (our paper tube) of mass m and radius r:
Roll Moment from Canted Fins
Each canted fin produces a lift force perpendicular to the airflow. The component of this force in the tangential direction (creating torque) is:
Spin-Up Dynamics
The roll acceleration is governed by Newton's second law for rotation:
The damping torque increases linearly with spin rate: τdamping = Cdamp · ω. This gives a first-order ODE:
The solution is an exponential approach to equilibrium:
The time constant is τ = I / Cd. For a typical paper rocket, spin-up to 63% of max spin rate happens within a few rocket-lengths of travel.
Why Spin Stabilization Works
A disturbance (like a gust of wind) applies a pitch/yaw torque. For a non-spinning rocket, this torque directly rotates the rocket off-axis. For a spinning rocket, the gyroscopic effect converts this torque into precession — a slow, cone-shaped wobble rather than a tumble. The higher the spin rate relative to the disturbance torque, the smaller the precession angle.
7. Altitude Measurement
Without electronic altimeters, we measure altitude using basic trigonometry. An observer stands a known distance from the launcher and measures the elevation angle to the rocket at apogee.
Single Observer Method
- H = altitude of rocket at apogee (above ground)
- d = horizontal distance from observer to launcher (m)
- θ = elevation angle measured by observer (from horizontal)
- heye = observer's eye height above ground (m)
Error Analysis
The sensitivity of altitude to angle measurement error is found by differentiating:
The altitude error from a small angle measurement error δθ (in radians):
Optimal Observer Distance
A good rule of thumb: stand at a distance approximately equal to the expected altitude. This gives angles near 45 degrees where the tangent function is well-behaved and errors are moderate.
Multiple Observer Triangulation
Two observers stationed at known positions can triangulate the rocket's 3D position. If observer A is at distance dA measuring angle θA and azimuth φA, and observer B likewise:
With two observers, you can also account for wind drift. If the rocket drifts laterally, a single observer gets the wrong distance. Two observers at 90 degrees to each other can resolve the actual position.
The best estimate of altitude combines both measurements, weighted by their estimated uncertainties:
Full Error Propagation
Using standard error propagation for H = d · tan(θ), the variance in H from independent errors in d and θ:
Evaluating the partial derivatives:
The second term dominates at high angles, confirming why angle accuracy matters more than distance accuracy for steep elevation angles.
8. Drag Analysis from Flight Timing
By timing the ascent and descent of a vertical launch, we can work backwards to estimate the rocket's drag coefficient. This is one of the most powerful field techniques available without electronic instrumentation.
Ascent Phase: Deceleration
During ascent, both gravity and drag act downward (opposing upward motion). The deceleration is:
The ascent time is always shorter than it would be without drag. The ratio of actual ascent time to vacuum ascent time (tvac = v0/g) tells us how much drag matters.
Descent Phase: Terminal Velocity
During descent, gravity pulls the rocket down while drag pushes up. The rocket accelerates until these forces balance:
Solving for terminal velocity:
Deriving Terminal Velocity from the ODE
During descent (taking downward as positive), the equation of motion is:
Setting dv/dt = 0 (no more acceleration) gives the terminal velocity equation above. We can also solve the full ODE. Let k = ρ Cd A / (2m) and note Vt = √(g/k). The equation becomes:
This is separable. Integrating with initial condition v(0) = 0 (starting from rest at apogee):
The left side is Vt · arctanh(v/Vt). Solving for v:
Integrating once more gives the distance fallen:
The tanh function approaches 1 asymptotically, meaning the rocket approaches terminal velocity exponentially. It reaches 99% of Vt after falling approximately 5.3 · Vt2 / g meters.
Estimating Cd from Measured Terminal Velocity
If you can estimate the descent rate (from altitude and descent time), you can invert the terminal velocity equation to find Cd:
Practical Procedure
- 1 Launch vertically and measure total flight time ttotal and max altitude H (via tracking).
- 2 Time the ascent tup (launch to apogee) and descent tdown (apogee to landing) separately.
- 3 Estimate average descent velocity: Vavg ≈ H / tdown.
- 4 For a long descent (rocket falls > 5 · Vt2/g), the average speed approximates Vt.
- 5 Plug Vt into the Cd equation above with known mass, air density, and reference area.
Ascent Timing Method
An alternative uses ascent time alone. If you know the muzzle velocity v0 (from the simulator) and measure the ascent time tup, the effective drag can be estimated. Without drag, the ascent time would be v0 / g. The ratio of actual to vacuum ascent time constrains Cd.
Ascent Time with Drag (Vertical Launch)
During ascent (upward positive), with k = ρ Cd A / (2m):
Separating variables and integrating from v0 to 0:
If you measure tup and know v0, you can numerically solve for Vt, and then find Cd.
Consistency Checks
- Asymmetry ratio: tup / tdown should be less than 1 (ascent is always faster than descent when drag is present). The further below 1, the more drag matters.
- Cross-check: Cd estimated from descent should match Cd from ascent timing. Large discrepancies suggest the rocket's attitude changed (tumbling on descent, for example).
- Multiple flights: Average Cd over several launches to account for measurement noise and wind effects.
Quick Reference: Notation & Constants
Physical Constants
| Symbol | Quantity | Value |
|---|---|---|
| g | Gravitational acceleration | 9.81 m/s2 |
| ρair | Air density (sea level, 15 °C) | 1.225 kg/m3 |
| μ | Dynamic viscosity of air | 1.81 × 10−5 Pa·s |
| γ | Ratio of specific heats (air) | 1.4 |
| Patm | Standard atmospheric pressure | 101,325 Pa |
| R | Specific gas constant (air) | 287 J/(kg·K) |
Common Symbols
| Symbol | Meaning | Typical Units |
|---|---|---|
| m | Rocket mass | kg (or g) |
| v | Velocity | m/s |
| Cd | Drag coefficient | dimensionless |
| A | Reference area (frontal) | m2 |
| CN | Normal force coefficient | dimensionless |
| Re | Reynolds number | dimensionless |
| β | Ballistic coefficient | kg/m2 |
| ω | Angular velocity (roll) | rad/s |
| SM | Stability margin | calibers (body diameters) |