Every equation behind the paper rocket simulator, from pneumatic launch to landing
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.
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.
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:
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:
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:
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.
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:
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:
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.
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:
Equivalently, using the circumference of the mandrel to find the paper width per wrap:
Tape adds significant mass. Common types:
For a single wrap of tape around the circumference: Lapplied = π · d.
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:
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:
Total fin system mass for n fins:
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 key insight: drag scales with velocity squared. Doubling speed quadruples drag force. This is why paper rockets decelerate so rapidly after launch.
The total Cd is a sum of contributions from different rocket components:
The Reynolds number determines whether the boundary layer is laminar (smooth, less drag) or turbulent (rough, more drag):
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.
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).
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.
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.
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):
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.
The simulator uses the Euler method: step forward in small time increments Δt, updating position and velocity at each step.
At each time step:
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:
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.
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.
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.
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.
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.
CG is the mass-weighted average position of all components, measured from the nose tip:
Component centroid positions (measured from nose tip):
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:
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.
For n identical trapezoidal fins, the normal force coefficient and CP location are:
The CP of the fins, measured from the fin leading edge at the root:
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.
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.
Converting to revolutions per second:
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.
For a thin-walled cylinder (our paper tube) of mass m and radius r:
Each canted fin produces a lift force perpendicular to the airflow. The component of this force in the tangential direction (creating torque) is:
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.
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.
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.
The sensitivity of altitude to angle measurement error is found by differentiating:
The altitude error from a small angle measurement error δθ (in radians):
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.
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:
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.
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.
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.
During descent, gravity pulls the rocket down while drag pushes up. The rocket accelerates until these forces balance:
Solving for terminal velocity:
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.
If you can estimate the descent rate (from altitude and descent time), you can invert the terminal velocity equation to find Cd:
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.
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.
| 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) |
| 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) |