The physics and mathematics for simulating a two dimensional rocket using lagrangian mechanics.
May 25, 2019
A two dimensional rocket is a much more complicated system contra a simple inverted pendulum. Though they are different in complexity they are similar regarding what they physically are. A two dimensional rocket, is simply, a inverted pendulum moving in free space. This also comes with a force acting on the rocket from down below, usually a rocket engine. This rocket engine also has an angle relative to the rocket if it is gimballed. This is what is going to get modelled in this post, a two dimensional rocket with a force acting on the rocket at an angle. As before, the equations of motion are the essential parts of the simulation of any system, but here, normal newtonian vector addition ect is not adequate. A new method and framework will be introduced and used, which is lagrangian mechanics. Besides this, some tweaks are done to the discretisation of the E.O.M in order to make the functions of motion working.
This section will firstly go through partial derivatives and time derivatives to build up to an understanding of lagrangian mechanics.
A partial derivative is a normal derivative but for just one variable where the others are constant. Consider the function
\[f(x, y) = x^2 + y^2\]
which would have the two partial derivatives
\[\frac{\partial f}{\partial x} = 2x\]
and
\[\frac{\partial f}{\partial y} = 2y\]
All rules of derivation apply, thinking of the other variables as an constant \(c\) and then deriving is helpful. If a functions contains \(x\) and \(\dot{x}\) the partial with respect to either \(x\) or \(\dot{x}\) still considers the other variable as constant.
A time derivative is as any other derivative but with respect to time but is treated as a function in the lagrangian equation. Consider the time dependent variable \(x\) which would have the time derivative of \(\dot{x}\) and the second order time derivative of \(\ddot{x}\). When taking the time derivative of a compound function the variables are treated as functions. For example given
\[f(x, y) = x + cos(y)\]
The time deriative of this function would be
\[\frac{d}{dt}\left(f(x, y)\right) = \dot{x} - sin(y)\dot{y}\]
Lagrangian mechanics is a very nice and easy framework one can use in order to find the equations of motion any system. Below is a list of the steps that need to be taken in order to find the equations of motion. That is, using the lagrangian method.
The linear kinetic energy is
\[E_k = \frac{1}{2}m\dot{x}^2 + \frac{1}{2}m\dot{y}^2 + \frac{1}{2}I\dot{\theta}^2\]
and the inertia is
\[I = \frac{1}{12}mL^2\]
together with the potential energy
\[E_p = mgy\]
forms the Lagrangian
\[\mathcal{L} = \frac{1}{2}m\dot{x}^2 + \frac{1}{2}m\dot{y}^2 + \frac{1}{2}I\dot{\theta}^2 - mgy\]
Taking the partial derivatives from the Lagrangian
\[ \begin{aligned} \frac{\partial \mathcal{L}}{\partial x} &= 0, \quad \frac{\partial \mathcal{L}}{\partial y} = -mg, \quad \frac{\partial \mathcal{L}}{\partial \theta} = 0, \\ \frac{\partial \mathcal{L}}{\partial \dot{x}} &= m\dot{x}, \quad \frac{\partial \mathcal{L}}{\partial \dot{y}} = m\dot{y}, \quad \frac{\partial \mathcal{L}}{\partial \dot{\theta}} = I\dot{\theta} \end{aligned} \]
and setting up the Euler-Lagrange equations gives us the equations of motion
\[\begin{aligned} \frac{d}{dt}\left(\frac{\partial \mathcal{L}}{\partial \dot{x}}\right) - \frac{\partial \mathcal{L}}{\partial x} &= m\ddot{x} = F_x \\ \frac{d}{dt}\left(\frac{\partial \mathcal{L}}{\partial \dot{y}}\right) - \frac{\partial \mathcal{L}}{\partial y} &= m\ddot{y} + mg= F_y \\ \frac{d}{dt}\left(\frac{\partial \mathcal{L}}{\partial \dot{\theta}}\right) - \frac{\partial \mathcal{L}}{\partial \theta} &= I\ddot{\theta} = \tau \end{aligned}\]
Solving these gives us
\[ \begin{aligned} \ddot{x} &= \frac{F_x}{m} \\ \\ \ddot{y} &= \frac{F_y - mg}{m} \\ \\ \ddot{\theta} &= \frac{\tau}{I} \end{aligned} \]
Lastly discretizing them
\[\begin{aligned} x_{k + 1} &= h^2\frac{F_x}{m} + 2x_k - x_{k-1}\\ \\ y_{i + 1} &= h^2\frac{F_y-mg}{m} + 2y_i - y_{i-1}\\ \\ \theta_{j + 1} &= h^2\frac{\tau}{I} + 2\theta_j - \theta_{j-1}\\ \\ \end{aligned} \]
From the figure above the following can be deduced
\[F_x = Fsin(\theta - \alpha)\]
\[F_y = Fcos(\theta - \alpha)\]
\[\tau = LFsin(\alpha)\]
The commented computer code for this simulation can be found on github.com/fullnitrous/two-dimensional-rocket-example and a in-browser simulation is shown below for the web version of this document.