Simulating a rocket launch

One can make crude estimations of rocket trajectories analytically, usually neglecting the drag force and make assumptions of the gravitation. To do much better estimations of rocket trajectory and other parameters one can do numerical simulation of rocket launches. Let us first see the the equations of the forces on a launching rocket from the last section:

$m\,\dfrac{\text{d}v_{\parallel}}{\text{d}t} = T \cos\left( \alpha + \delta \right) - m g_0 \, \dfrac{r_0}{r} \, \sin \gamma - D \, ,$

$m\,\dfrac{\text{d}v_{\perp}}{\text{d}t} = T \sin\left( \alpha + \delta \right) - m g_0 \, \dfrac{r_0}{r} \, \cos \gamma - L + m \, \dfrac{v^2}{r} \, \cos \gamma \, ,$

These are differential equations that are perfect for numerical integration. We change the equations slightly to be in the horizontal and vertical directions instead:

$\dfrac{\text{d}\vec{v}}{\text{d}t} = \dfrac{1}{m} \left( T \hat{v} - m g_0 \, \dfrac{r_0}{r} \, \hat{k} - D \hat{v} \right),$

$\dfrac{\text{d}\vec{r}}{\text{d}t} = \vec{v} \, ,$

where $\hat{v}$ is the unit vector in the travelling direction, given by

$\hat{v} = \dfrac{\vec{v}}{|v|} \, ;$

while $\hat{k}$ is the unit vector from the center of the Earth to the rocket; $D$ is the drag, as defined on the previous page; and $\vec{r}$ is the position vector from the launch site to the rocket. The mass $m$ is a function of time, since it decreases as the rocket motor is burning.

Now we have a physical model of the rocket launch. We have not taken roll and other effects into account. Inducing spin takes some energy, and when leaving that out we overestimate the apogee altitude. We want to find the acceleration (which is given by the equation above), velocity and position as a function of time. Hence we need to integrate the equation twice.

So, how can we find the acceleration when it is dependent on the integrated result? The answer is that we solve the equations iteratively. We limit this case to two spatial dimensions. We start first with setting the initial values:

$\vec{v}_0 = \left[\begin{array}{c} 0 \\ 0 \end{array}\right]$

and

$\vec{r}_0 = \left[\begin{array}{c} 0 \\ 0 \end{array}\right].$

Then we find the acceleration at lift-off at $t = \Delta t$ from the acceleration equation above.

We need now to integrate the acceleration. The derivative is given as

$\vec{a} = \dfrac{\Delta \vec{v}}{\Delta t} \, ,$

which can be rewritten as

$\vec{v}_{i+1} = \vec{v}_i + \vec{a}_i \, \Delta t \, .$

Using this simple method we can do simple integrations. The method is called the Forward Euler integrator, and is an explicit solver (i.e. it only uses information from the previous timestep). This is not an accurate solver, and better approaches do exist, like the Runge–Kutta family.

When we have integrated the acceleration to velocity in the two spatial dimensions, we integrate the velocity to get the position. Then we have the acceleration, velocity and position at time $t = \Delta t$ . Then we start once again with updating the acceleration based on the velocity and position at the previous timestep and integrate to get the acceleration, velocity and position at time $t = 2 \, \Delta t$ . We do this until the vertical position of the rocket is negative, that is, the rocket has landed or after a given time. We store all acceleration, velocity and position points, and from this we can derive many parameters, e.g. the pitch angle.

The rocket is stable because of the velocity that keeps it pointing forward (small angle of attack). At the start of the launch before the rocket has gained enough velocity to be stable, we need to implement a launch rod to keep it pointing at the launching angle. This is done by fixing the unit vector to the launch angle while $r < 2\; \text{m}$ (or whatever the length of the rod is).

A few comments before we continue. We have assumed a flat Earth. As most people know, the Earth is indeed spherical (just don’t ask the flat Earth society if this is true), but for a small rocket this can be a good approximation. Also, we have used a constant time step. For numerically heavy simulations it can be smart to implement adaptive time steps to keep the simulation time to a minimum and a more accurate solution. And again, we have only assumed three degrees of freedom (DOF), which will overestimate the altitude, range and velocity. Such a simulation can be implemented in for example Python or MATLAB in about 20–30 lines of code.

Another approach is to use software designed to do this. For simple simulations of single stage rockets which are not intended to go into orbit, the open-source 6-DOF OpenRocket software can be recommended. For a professional solution, AGI System Tool Kit is recommended (but it has a steep learning curve).

<< Previous page – Content – Next page >>

This article is part of a pre-course program used by Andøya Space Education in Fly a Rocket! and similar programs.