Rocket Maths

Rocket Mechanics


Referring to the diagram and using a earth fixed polar coordinate system. The vehicle is assumed to be a point mass m=m(t).
\begin{align} {\ddot r} & = {r{\dot \theta^2}} + \frac{F_r}{m}-\frac{\mu}{r^2} \\ {\ddot \theta} & = -2{\dot r}{\dot \theta} + \frac{F_\theta}{m} \end{align} These are the Newtonian state equations for r and θ. The first term RHS of the r equation represents the centrifugal acceleration due to a curved flight path - balances the gravitational acceleration (third term) on a circular orbit. \begin{align} {\mu} = {G}{M}\\ \end{align} where G is the universal gravitational constant (6.6743 x 10-11 m3 kg-1 s-2), M is the mass of the earth (5.9722 x 1024) kg. The first term RHS of the θ equation is a Coriolis acceleration due to r and θ velocity coupling.
As fuel is burned, the mass decreases with time: \begin{align} {\dot m} = - {\dot m_f} \\ \end{align} where mf is the remaining fuel mass. Thrust is assumed constant during the burn.
The radial and tangential forces are given by: \begin{align} {F_\theta} & = ({T-D})cos({\gamma}) \\ {F_r} & = ({T-D})sin({\gamma}) \end{align} where γ is the angle between the vehicle and the local horizontal, \begin{align} {T} = {I_{sp}}{g}{\dot {m_f}} \\ \end{align} where Isp is the specific impulse of the fuel and g is the sea level acceleration of gravity g = 9.81 m/s2.
The drag force is defined by \begin{align} {D} = {C_D}\frac{\rho}{2}{V^2}{A} \\ \end{align}
where \begin{align} {C_D} = f({M}) \\ \end{align} with the Mach number \begin{align} {M} = \frac{V}{c} \\ \end{align}
where c is the acoustic speed. The acoustic speed as well as air pressure, air density, and air temperature vary with altitude. These properties are calculated from an algorithm based on the 1976 US Standard Atmosphere, up to 600 km altitude. \begin{align} {\rho},{p},{T_s},{c} = f({r}) \\ \end{align} The drag coefficient is constant up to about M=0.8 and then increases sharply to a peak at M=1, after which it falls off to a value less than the subsonic value. The dynamic pressure is defined by \begin{align} {q} = \frac{\rho}{2}{V^2} \\ \end{align} The dynamic pressure is used to calculate aerodynamic forces (D) and is an indicator of stresses acting on the vehicle.
The vehicle fixed Cartesian coordinates and velocities are: \begin{align} {V_{\theta}} & = {\dot x} = {r}{\dot \theta} = {r}{\omega} \\ {V_{r}} & = {\dot r} \\ \end{align} with \begin{align} {V}^2 & = {V_{\theta}}^2 + {V_r}^2 \\ \end{align} \begin{align} {x} & = {r}{\theta} \\ {h} & = {r}-{R_E} \\ \end{align} The target orbital velocity (circular orbit) is \begin{align} {V_T} & = {V_{\theta}} =\sqrt{\frac{\mu}{r}} \\ \end{align} There is a orbital velocity boost due to the earths rotation: \begin{align} {V_{B}} = {R_E}{\omega_{E}}cos({\phi})cos({i}) \\ \end{align} where ωE is the earth rotation rate, φ is the launch latitude, and i is the launch inclination with 90° being a polar orbit. All velocities reported are relative to launch site, but orbital velocity, VA, is absolute. \begin{align} {V_{A}} = {V_{B}}+{V_{R}} \\ \end{align} The aerodynamic heat rate due to air friction is \begin{align} {P} = {D}{V} \\ \end{align} When the heat rate per unit mass of the vehicle exceeds a certain threshold, melting occurs and the mass of the vehicle is decreased accordingly. Typically this will happen at high Mach number (>5) at altitudes less than 80 km.

Orbit Geometry and Properties

The specific energy, ε, (J/kg) remains constant while in orbit. \begin{align} {\epsilon} = \frac{V^2}{2}-\frac{\mu}{r} \end{align} where V and r are the instantaneous velocity and position of the spacecraft. In general, an orbit will be be elliptic in shape. The semi-major axis of the orbit, a, is given by \begin{align} {a} = -\frac{\mu}{2\epsilon} \end{align} To calculate additional properties, the specific angular momentum, ℏ, is needed. Calculation require some vector algebra (not shown here), but starting with \begin{align} \vec{\hbar} = \vec{r} \times \vec{V} \end{align} From this the orbit ecentricity, e, is \begin{align} {e} = \sqrt{1 + \frac{2\epsilon\hbar^2}{\mu^2}} \end{align} A circular orbit has an eccentricity of zero, and an elliptic orbit e < 1; From this, the closest approach (perigee) is \begin{align} {r_p} = a(1-e) \end{align} and the farthest point (apogee) is \begin{align} {r_a} = a(1+e) \end{align} In terms of altitude \begin{align} {h_p} = {r_p}-{R_E} \end{align} and \begin{align} {h_a} = {r_a}-{R_E} \end{align} The period of the orbit is a function of the semi-major axis: \begin{align} {T} = 2\pi \sqrt{\frac{a^3}{\mu}} \end{align} Finally the True Anomaly, ν, the position angle measured from the perigee of the orbit, is found using the relationship \begin{align} \cos(\nu) = \frac{a(1-e^2)-r}{re} \end{align} This is used to determine the tilt of the orbit relative to the launch point.

...and of course, all programmed in Javascript.