X-Git-Url: https://git.gag.com/?a=blobdiff_plain;f=core%2Fdoc%2Ftechdoc%2Fchapter-flight-simulation.tex;fp=core%2Fdoc%2Ftechdoc%2Fchapter-flight-simulation.tex;h=f452103bfa41ad63ac399ccd08f3c197a72acd3a;hb=84094f3e8b4e8d27310532c092ec9738156417d3;hp=0000000000000000000000000000000000000000;hpb=fc3d9fa48747ea14f7d77e573df16386c9fdf89e;p=debian%2Fopenrocket diff --git a/core/doc/techdoc/chapter-flight-simulation.tex b/core/doc/techdoc/chapter-flight-simulation.tex new file mode 100644 index 00000000..f452103b --- /dev/null +++ b/core/doc/techdoc/chapter-flight-simulation.tex @@ -0,0 +1,781 @@ + + +\chapter{Flight simulation} +\label{chap-simulation} + +In this chapter the actual flight simulation is analyzed. First in +Section~\ref{sec-atmospheric-properties} methods for simulating +atmospheric conditions and wind are presented. Then in +Section~\ref{sec-flight-modeling} the actual simulation procedure is +developed. + + + +\section{Atmospheric properties} +\label{sec-atmospheric-properties} + +In order to calculate the aerodynamic forces acting on the rocket it +is necessary to know the prevailing atmospheric conditions. Since the +atmosphere is not constant with altitude, a model must be developed to +account for the changes. Wind also plays an important role in the +flight of a rocket, and therefore it is important to have a realistic +wind model in use during the simulation. + + +\subsection{Atmospheric model} + +The atmospheric model is responsible to estimating the atmospheric +conditions at varying altitudes. The properties that are of most +interest are the density of air $\rho$ (which is a scaling parameter +to the aerodynamic coefficients via the dynamic pressure +$\frac{1}{2}\rho v^2$) and the speed of sound $c$ (which affects the +Mach number of the rocket, which in turn affects its aerodynamic +properties). These may in turn be calculated from the air pressure +$p$ and temperature $T$. + +Several models exist that define standard atmospheric conditions as a +function of altitude, including the Internaltional Standard +Atmosphere (ISA)~\cite{international-standard-atmosphere} and the +U.S. Standard Atmosphere~\cite{US-standard-atmosphere}. These two +models yield identical temperature and pressure profiles for altitudes +up to 32~km. + +The models are based on the assumption that air follows the ideal gas +law +% +\begin{equation} +\rho = \frac{Mp}{RT} +\end{equation} +% +where $M$ is the molecular mass of air and $R$ is the ideal gas +constant. From the equilibrium of hydrostatic forces the differential +equation for pressure as a function of altitude $z$ can be found as +% +\begin{equation} +\dif p = -g_0 \rho \dif z = -g_0 \frac{Mp}{RT} \dif z +\label{eq-pressure-altitude} +\end{equation} +% +where $g_0$ is the gravitational acceleration. If the temperature of +air were to be assumed to be constant, this would yield an exponential +diminishing of air pressure. + +The ISA and U.S. Standard Atmospheres further specity a standard +temperature and pressure at sea level and a temperature profile for +the atmosphere. The temperature profile is given as eight +temperatures for different altitudes, which are then linearly +interpolated. The temperature profile and base pressures for the ISA +model are presented in Table~\ref{table-ISA-model}. These values +along with equation~(\ref{eq-pressure-altitude}) define the +temperature/pressure profile as a function of altitude. + +\begin{table} +\caption{Layers defined in the International Standard + Atmosphere~\cite{wiki-ISA-layers}} +\label{table-ISA-model} +\begin{center} +\begin{tabular}{ccccl} +\hline +Layer & Altitude$^\dagger$ & Temperature & Lapse rate & + \multicolumn{1}{c}{Pressure} \\ + & m & $^\circ$C & $^\circ$C/km & \multicolumn{1}{c}{Pa} \\ +\hline +0 & 0 & $+15.0$ & $-6.5$ & 101\s325 \\ +1 & 11\s000 & $-56.5$ & $+0.0$ & \num22\s632 \\ +2 & 20\s000 & $-56.5$ & $+1.0$ & \num\num5\s474.9 \\ +3 & 32\s000 & $-44.5$ & $+2.8$ & \num\num\num\s868.02 \\ +4 & 47\s000 & \num$-2.5$ & $+0.0$ & \num\num\num\s110.91 \\ +5 & 51\s000 & \num$-2.5$ & $-2.8$ & \num\num\num\s\num66.939 \\ +6 & 71\s000 & $-58.5$ & $-2.0$ & \num\num\num\s\num\num3.9564 \\ +7 & 84\s852 & $-86.2$ & & \num\num\num\s\num\num0.3734 \\ +\hline +\end{tabular} +\end{center} +\vspace{-3mm} +{\footnotesize $^\dagger$ Altitude is the geopotential height which + does not account for the diminution of gravity at high altitudes.} +\vspace{3mm} +\end{table} + +These models are totally static and do not take into account any local +flight conditions. Many rocketeers may be interested in flight +differences during summer and winter and what kind of effect air pressure +has on the flight. These are also parameters that can easily be +measured on site when launching rockets. On the other hand, it is +generally hard to know a specific temperature profile for a specific +day. Therefore the atmospheric model was extended to allow the user +to specify the base conditions either at mean sea level or at the +altitude of the launch site. These values are simply assigned to the +first layer of the atmospheric model. Most model rockets do not +exceed altitudes of a few kilometers, and therefore the flight +conditions at the launch site will dominate the flight. + +One parameter that also has an effect on air density and the speed of +sound is humidity. The standard models do not include any definition +of humidity as a function of altitude. Furthermore, the effect of +humidity on air density and the speed of sound is marginal. The +difference in air density and the speed of sound between completely dry +air and saturated air at standard conditions are both less than 1\%. +Therefore the effect of humidity has been ignored. + + + + +\subsection{Wind modeling} + +Wind plays a critical role in the flight of model rockets. As has +been seen, large angles of attack may cause rockets to lose a +significant amount of stability and even go unstable. Over-stable +rockets may weathercock and turn into the wind. In a perfectly static +atmosphere a rocket would, in principle, fly its entire flight +directly upwards at zero angle of attack. Therefore, the effect of +wind must be taken into account in a full rocket simulation. + +Most model rocketeers, however, do not have access to a full wind +profile of the area they are launching in. Different layers of air +may have different wind velocities and directions. Modeling such +complex patterns is beyond the scope of this project. Therefore, the +goal is to produce a realistic wind model that can be specified with +only a few parameters understandable to the user and that covers +altitudes of most rocket flights. Extensions to allow for multiple +air layers may be added in the future. + +In addition to a constant average velocity, wind always has some +degree of turbulence in it. The effect of turbulence can be modeled +by summing the steady flow of air and a random, zero-mean turbulence +velocity. Two central aspects of the turbulence velocity are the +amplitude of the variation and the frequencies at which they occur. +Therefore a reasonable turbulence model is achieved by a random +process that produces a sequence with a similar distribution and +frequency spectrum as that of real wind. + +Several models of the spectrum of wind turbulence at specific +altitudes exist. Two commonly used such spectra are the {\it Kaimal} +and {\it von Kármán} wind turbulence +spectra~\cite[p.~23]{wind-energy-handbook}: +% +\begin{eqnarray} +\mbox{Kaimal:} & & \frac{S_u(f)}{\sigma_u^2} = + \frac{4 L_{1u} / U}{(1 + 6fL_{1u}/U)^{5/3}} \label{eq-kaimal-wind} \\ +% +\mbox{von Kármán:} & & \frac{S_u(f)}{\sigma_u^2} = + \frac{4 L_{2u} / U}{(1 + 70.8(fL_{2u}/U)^2)^{5/6}} \label{eq-karman-wind} +\end{eqnarray} + +Here $S_u(f)$ is the spectral density function of the turbulence +velocity and $f$ the turbulence frequency, $\sigma_u$ the standard +deviation of the turbulence velocity, $L_{1u}$ and $L_{2u}$ length +parameters and $U$ the average wind speed. + +Both models approach the asymptotic limit +$S_u(f)/\sigma_u^2 \sim f^{-5/3}$ quite fast. Above frequencies of +0.5~Hz the difference between equation~(\ref{eq-kaimal-wind}) and the +same equation without the term 1 in the denominator is less than 4\%. +Since the time scale of a model rocket's flight is quite short, the +effect of extremely low frequencies can be ignored. Therefore +turbulence may reasonably well be modelled by utilizing +{\it pink noise} that has a spectrum of $1/f^\alpha$ with $\alpha=5/3$. +True pink noise has the additional useful property of being +scale-invariant. This means that a stream of pink noise samples may +be generated and assumed to be at any sampling rate while maintaining +their spectral properties. + +Discerete samples of pink noise with spectrum $1/f^\alpha$ can be +generated by applying a suitable digital filter to {\it white noise}, +which is simply uncorrelated pseudorandom numbers. One such filter is +the infinite impulse response (IIR) filter presented by +Kasdin~\cite{pink-filter}: +% +\begin{equation} +x_n = w_n - a_1 x_{n-1} - a_2 x_{n-2} - a_3 x_{n-3} - \ldots +\label{eq-pink-generator} +\end{equation} +% +where $x_i$ are the generated samples, $w_n$ is a generated white +random number and the coefficients are computed using +% +\begin{equation} +\begin{array}{rl} +a_0 & = 1 \\ +a_k & = \del{k-1-\frac{\alpha}{2}} \frac{a_{k-1}}{k}. +\end{array} +\label{eq-pink-coefficients} +\end{equation} +% +The infinite sum may be truncated with a suitable number of terms. +In the context of IIR filters these terms are calles {\it poles}. +Experimentation showed that already 1--3 poles provides a reasonably +accurate frequency spectrum in the high frequency range. + +One problem in using pink noise as a turbulence velocity +model is that the power spectrum of pure pink noise goes to +infinity at very low frequencies. This means that a long sequence +of random values may deviate significantly from zero. However, when +using the truncated IIR filter of equation~(\ref{eq-pink-generator}), +the spectrum density becomes constant below a certain limiting +frequency, dependent on the number of poles used. By adjusting the +number of poles used, the limiting frequency can be adjusted to a value +suitable for model rocket flight. Specifically, the number of poles +must be selected such that the limiting frequency is suitable at the +chosen sampling rate. + +It is also desirable that the simulation resolution does not affect +the wind conditions. For example, a simulation with a time step of +10~ms should experience the same wind conditions as a simulation with +a time step of 5~ms. This is achieved by selecting a constant +turbulence generation frequency and interpolating between the +generated points when necessary. The fixed frequency was chosen at +20~Hz, which can still simulate fluctuations at a time scale of 0.1 +seconds. + +\begin{figure}[p] +\centering +\epsfig{file=figures/wind/pinktime, width=105mm} +\caption{The effect of the number of IIR filter poles on two 20 second + samples of generated turbulence, normalized so that the two-pole + sequence has standard deviation one.} +\label{fig-pink-poles} +\end{figure} + +\begin{figure}[p] +\centering +\epsfig{file=figures/wind/pinkfreq, width=95mm} +\caption{The average power spectrum of 100 turbulence + simulations using a two-pole IIR filter (solid) and the Kaimal + turbulence spectrum (dashed); vertical axis arbitrary.} +\label{fig-pink-spectrum} +\end{figure} + + + +The effect of the number of poles is depicted in +Figure~\ref{fig-pink-poles}, where two pink noise sequences were +generated from the same random number source with two-pole and +ten-pole IIR filters. A small number of poles generates values strongly +centered on zero, while a larger number of poles introduces more low +frequency variability. Since the free-flight time of a typical model +rocket is of the order of 5--30 seconds, it is desireable that the +maximum gust length during the flight is substantially shorter than +this. Therefore the pink noise generator used by the wind model was +chosen to contain only two poles, which has a limiting frequency of +approximately 0.3~Hz when sampled at 20~Hz. This means that gusts of +wind longer than 3--5 seconds will be rare in the simulted turbulence, +which is a suitable gust length for modeling typical model rocket +flight. Figure~\ref{fig-pink-spectrum} depicts the resulting pink +noise spectrum of the two-pole IIR filter and the Kaimal spectrum of +equation~(\ref{eq-kaimal-wind}) scaled to match each other. + + +%, which causes frequency +%components below approximately 0.3~Hz to be subdued. Therefore, gusts +%of wind longer than 3--5 seconds will be rare in the simulated wind, a +%suitable time scale for the flight of a model rocket. +%Figure~\ref{fig-turbulence}(a) shows a 20 second sample of the +%generated turbulence, normalized to have a standard deviation of one. +%Figure~\ref{fig-turbulence}(b) depicts the actual frequency spectrum +%of the generated turbulence and the Kaimal spectrum of +%equation~(\ref{eq-kaimal-wind}) scaled to match each other. + + + +To simplify the model, the average wind speed is assumed to be +constant with altitude and in a constant direction. This allows +specifying the model parameters using just the average wind speed and +its standard deviation. An alternative parameter for specifying the +turbulence amplitude is the {\it turbulence intensity}, which is the +percentage that the standard deviation is of the average wind +velocity, +% +\begin{equation} +I_u = \frac{\sigma_u}{U}. +\end{equation} +% +Wind farm load design standards typically specify turbulence +intensities around 10\ldots20\%~\cite[p.~22]{wind-energy-handbook}. +It is assumed that these intensities are at the top of the range of +conditions in which model rockets are typically flown. + +Overall, the process to generate the wind velocity as a function of +time from the average wind velocity $U$ and standard deviation +$\sigma_u$ can be summarized in the following steps: +% +\begin{enumerate} +%\item[Input:] Average wind velocity $U$ and standard deviation +% $\sigma_u$. +% +\item Generate a pink noise sample $x_n$ from a Gaussian white noise + sample $w_n$ using equations~(\ref{eq-pink-generator}) and + (\ref{eq-pink-coefficients}) with two memory terms included. + +\item Scale the sample to a standard deviation one. This is performed + by dividing the value by a previously calculated standard deviation + of a long, unscaled pink noise sequence (2.252 for the two-pole IIR + filter). + +\item The wind velocity at time $n\cdot\Delta t$ ($\Delta t = 0.05\rm~s$) + is $U_n = U + \sigma_u x_n$. Velocities in between are interpolated. +\end{enumerate} + + + + +\section{Modeling rocket flight} +\label{sec-flight-modeling} + + +Modeling of rocket flight is based on Newton's laws. The basic forces +acting upon a rocket are gravity, thrust from the motors and +aerodynamic forces and moments. These forces and moments are +calculated and integrated numerically to yield a simulation over a +full flight. + +Since most model rockets fly at a maximum a few kilometers high, the +curvature of the Earth is not taken into account. Assuming a flat +Earth allows us to use simple Cartesian coordinates to represent the +position and altitude of the rocket. As a consequence, the coriolis +effect when flying long distances north or south is not simulated +either. + + + +\subsection{Coordinates and orientation} + +During a rocket's flight many quantities, such as the aerodynamical +forces and thrust from the motors, are relative to the rocket itself, +while others, such as the position and gravitational force, are more +naturally described relative to the launch site. Therefore two sets +of coordinates are defined, the {\it rocket coordinates}, which are +the same as used in Chapter~\ref{chap-aerodynamics}, and +{\it world coordinates}, which is a fixed coordinate system with the +origin at the position of launch. + +The position and velocity of a rocket are most naturally maintained as +Cartesian world coordinates. Following normal convensions, the +$xy$-plane is selected to be parallel to the ground and the $z$-axis +is chosen to point upwards. In flight dynamics of aircraft the +$z$-axis often points towards the earth, but in the case of rockets it +is natural to have the rocket's altitude as the $z$-coordinate. + +Since the wind is assumed to be unidirectional and the Coriolis effect +is ignored, it may be assumed that the wind is directed along the +$x$-axis. The angle of the launch rod may then be positioned relative +to the direction of the wind without any loss of generality. + +Determining the orientation of a rocket is more complicated. A +natural choise for defining the orientation would be to use the +spherical coordinate zenith and azimuth angles $(\theta, \phi)$ and an +additional roll angle parameter. Another choise common in aviation is +to use {\it Euler angles}~\cite{wiki-euler-angles}. However, both of +these systems have notable shortcomings. Both systems have +singularity points, in which the value of some parameter is +ambiguous. With spherical coordinates, this is the direction of the +$z$-axis, in which case the azimuth angle $\phi$ has no effect on the +position. Rotations that occur near these points must often be +handled as special cases. Furthermore, rotations in spherical +coordinate systems contain complex trigonometric formulae which are +prone to programming errors. + +The solution to the singularity problem is to introduce an extra +parameter and an additional constraint to the system. For example, +the direction of a rocket could be defined by a three-dimensional unit +vector $(x,y,z)$ instead of just the zenith and azimuth angles. The +additional constraint is that the vector must be of unit length. This +kind of representation has no singularity points which would require +special consideration. + +Furthermore, Euler's rotation theorem states that a rigid body can be +rotated from any orientation to any other orientation by a single +rotation around a specific axis~\cite{wiki-euler-rotation-theorem}. +Therefore instead of defining quantities that define the orientation +of the rocket we can define a three-dimensional rotation that rotates +the rocket from a known reference orientation to the current +orientation. This has the additional advantage that the same rotation +and its inverse can be used to transform any vector between world +coordinates and rocket coordinates. + +A simple and efficient way of descibing the 3D rotation is by using +{\it unit quaternions}. Each unit quaternion corresponds to a unique +3D rotation, and they are remarkably simple to combine and use. The +following section will present a brief overview of the properties of +quaternions. + +The fixed reference orientation of the rocket defines the rocket +pointing towards the positive $z$-axis in world coordinates and an +arbitrary but fixed roll angle. The orientation of the rocket is then +stored as a unit quaternion that rotates the rocket from this +reference orientation to its current orientation. +This rotation can also be used to transform vectors from world +coordinates to rocket coordinates and its inverse from rocket +coordinates to world coordinates. (Note that the rocket's initial +orientation on the launch pad may already be different than its +reference orientation if the launch rod is not completely vertical.) + + + + +\subsection{Quaternions} + +{\it Quaternions} are an extension of complex numbers into four +dimensions. The usefulness of quaternions arises from their use in +spatial rotations. Similar to the way multiplication with a complex +number of unit length $e^{i\phi}$ corresponds to a rotation of angle +$\phi$ around the origin on the complex plane, multiplication with +unit quaternions correspond to specific 3D rotations around an axis. +A more thorough review of quaternions and their use in spatial +rotations is available in Wikipedia~\cite{wiki-quaternion-rotations}. + +The typical notation of quaternions resembles the addition of a scalar +and a vector: +% +\begin{equation} +q = w + x\vi + y\vj + z\vk = w + \vect v +\end{equation} +% +Addition of quaternions and multiplication with a scalar operate as +expected. However, the multiplication of two quaternions is +non-commutative (in general $ab \neq ba$) and follows the rules +% +\begin{equation} +\vi^2 = \vj^2 = \vk^2 = \vi\vj\vk = -1. +\end{equation} +% +As a corollary, the following equations hold: +% +\begin{equation} +\begin{array}{rl} +\vi\vj = \vk \hspace{15mm}& \vj\vi = -\vk \\ +\vj\vk = \vi \hspace{15mm}& \vk\vj = -\vi \\ +\vk\vi = \vj \hspace{15mm}& \vi\vk = -\vj +\end{array} +\end{equation} +% +The general multiplication of two quaternions becomes +% +\begin{equation} +\begin{array}{rl} +(a + b\vi + c\vj + d\vk)(w + x\vi + y\vj + z\vk)\;\; = + & (aw-bx-cy-dz) \\ + & + (ax+bw+cz-dy)\;\vi \\ + & + (ay-bz+cw+dx)\;\vj \\ + & + (az+by-cx+dw)\;\vk +\end{array} +\end{equation} +% +while the norm of a quaternion is defined in the normal manner +% +\begin{equation} +|q| = \sqrt{w^2+x^2+y^2+z^2}. +\end{equation} + +The usefulness of quaternions becomes evident when we consider a +rotation around a vector $\vect u$, $|\vect u|=1$ by an angle $\phi$. +Let +% +\begin{equation} +q = \cos\frac{\phi}{2} + \vect u \sin\frac{\phi}{2}. +\label{eq-rotation-quaternion} +\end{equation} +% +Now the previously mentioned rotation of a three-dimensional vector +$\vect v$ defined by $\vi$, $\vj$ and $\vk$ is equivalent to the +quaternion product +% +\begin{equation} +\vect v \mapsto q\vect v q^{-1}. +\end{equation} +% +Similarly, the inverse rotation is equivalent to the transformation +% +\begin{equation} +\vect v \mapsto q^{-1} \vect v q. +\end{equation} +% +The problem simplifies even further, since for unit quaternions +% +\begin{equation} +q^{-1} = (w + x\vi + y\vj + z\vk)^{-1} = w - x\vi - y\vj - z\vk. +\end{equation} +% +Vectors can therefore be considered quaternions with no scalar +component and their rotation is equivalent to the left- and right-sided +multiplication with unit quaternions, requiring a total of 24 +floating-point multiplications. Even if this does not make the +rotations more efficient, it simplifies the trigonometry considerably +and therefore helps reduce programming errors. + + +\subsection{Mass and moment of inertia calculations} +\label{sec-mass-inertia} + +Converting the forces and moments into linear and angluar acceleration +requires knowledge of the rocket's mass and moments of inertia. The +mass of a component can be easily calculated from its volume and +density. Due to the highly symmetrical nature of rockets, the rocket +centerline is commonly a principal axis for the moments of inertia. +Furthermore, the moments of inertia around the in the $y$- and +$z$-axes are very close to one another. Therefore as a simplification +only two moments of inertia are calculated, the longitudal and +rotational moment of inertia. These can be easily calculated for each +component using standard formulae~\cite{wiki-moments-of-inertia} and +combined to yield the moments of the entire rocket. + +This is a good way of calculating the mass, CG and inertia of a rocket +during the design phase. However, actual rocket components often have +a slightly different density or additional sources of mass such as +glue attached to them. These cannot be effectively modeled by the +simulator, since it would be extremely tedious to define all these +properties. Instead, some properties of the components can be +overridden to utilize measured values. + +Two properties that can very easily be measured are the mass and +CG position of a component. Measuring the moments of inertia is a +much harder task. Therefore the moments of inertia are still computed +automatically, but are scaled by the overridden measurement values. + +If the mass of a component is overridden by a measured value, the +moments of inertia are scaled linearly according to the mass. This +assumes that the extra weight is distributed evenly along the +component. If the CG position is overridden, there is no knowledge +where the extra weight is at. Therefore as a best guess the moments +of inertia are updated by shifting the moment axis according to the +parallel axis theorem. + +As the components are computed individually and then combined, the +overriding can take place either for individual components or larger +combinations. It is especially useful is to override the mass and/or CG +position of the entire rocket. This allows constructing a rocket from +components whose masses are not precisely known and afterwards scaling +the moments of inertia to closely match true values. + + + +\subsection{Flight simulation} + +The process of simulating rocket flight can be broken down into the +following steps: + +\begin{enumerate} +\setcounter{enumi}{-1} +\item Initialize the rocket in a known position and orientation at + time $t=0$. +\item Compute the local wind velocity and other atmospheric conditions. +\item Compute the current airspeed, angle of attack, lateral wind + direction and other flight parameters. +\item Compute the aerodynamic forces and moments affecting the rocket. +\item Compute the effect of motor thrust and gravity. +\item Compute the mass and moments of inertia of the rocket and from + these the linear and rotational acceleration of the rocket. +\item Numerically integrate the acceleration to the rocket's position + and orientation during a time step $\Delta t$ and update the current + time $t \mapsto t+\Delta t$. +\end{enumerate} + +Steps 1--6 are repeated until an end criteria is met, typically until +the rocket has landed. + +The computation of the atmospheric properties and instantaneous wind +velocity were discussed in Section~\ref{sec-atmospheric-properties}. +The local wind velocity is added to the rocket velocity to get the +airspeed velocity of the rocket. By inverse rotation this quantity is +obtained in rocket coordinates, from which the angle of attack and +other flight parameters can be computed. + +After the instantaneous flight parameters are known, the aerodynamic +forces can be computed as discussed in +Chapter~\ref{chap-aerodynamics}. The computed forces are in the +rocket coordinates, and can be converted to world coordinates by +applying the orientation rotation. The thrust from the motors is +similarly calculated from the thrust curves and converted to world +coordinates, while the direction of gravity is already in world +coordinates. When all of the the forces and moments acting upon the +rocket are known, the linear and rotational accelerations can be +calculated using the mass and moments of inertia discussed in +Section~\ref{sec-mass-inertia}. + +The numerical integration is performed using the Runge-Kutta~4 (RK4) +integration method. In order to simulate the differential equations +% +\begin{equation} +\begin{split} +x''(t) &= a(t) \\ +\phi''(t) &= \alpha(t) +\end{split} +\end{equation} +% +the equation is first divided into first-order equations using the +substitutions $v(t)=x'(t)$ and $\omega(t)=\phi'(t)$: +% +\begin{equation} +\begin{split} +v'(t) &= a(t) \\ +x'(t) &= v(t) \\ +\omega'(t) &= \alpha(t) \\ +\phi'(t) &= \omega(t) +\end{split} +\end{equation} +% +For brevity, this is presented in the first order representation +% +\begin{equation} +y' = f(y,\; t) +\end{equation} +% +where $y$ is a vector function containing the position and orientation +of the rocket. + +Next the right-hand side is evaluated at four positions, dependent on +the previous evaluations: +% +\begin{equation} +\begin{split} +k_1 &= f(y_0,\; t_0) \\ +k_2 &= f(y_0 + k_1\:\mbox{$\frac{\Delta t}{2}$},\; + t_0 + \mbox{$\frac{\Delta t}{2}$}) \\ +k_3 &= f(y_0 + k_2\:\mbox{$\frac{\Delta t}{2}$},\; + t_0 + \mbox{$\frac{\Delta t}{2}$}) \\ +k_4 &= f(y_0 + k_3\:\Delta t,\; t_0 + \Delta t) +\end{split} +\end{equation} +% +Finally, the result is a weighted sum of these values: +% +\begin{align} +y_1 &= y_0 + \frac{1}{6}\left(k_1+2k_2+2k_3+k_4\right)\,\Delta t \\ +t_1 &= t_0 + \Delta t +\end{align} + +Computing the values $k_1\ldots k_4$ involves performing steps~1--5 +four times per simulation iteration, but results in significantly +better simulation precision. The method is a fourth-order integration +method, meaning that the error incurred during one simulation step is +of the order $O(\Delta t^5)$ and of the total simulation +$O(\Delta t^4)$. This is a considerable improvement +over, for example, simple Euler integration, which has a total error +of the order $O(\Delta t)$. Halving the time step in an Euler +integration only halves the total error, but reduces the error of a +RK4 simulation 16-fold. + +The example above used a total rotation vector $\phi$ to contain the +orientation of the rocket. Instead, this is replaced by the rotation +quaternion, which can be utilized directly as a transformation between +world and rocket coordinates. Instead of updating the total rotation +vector, +% +\begin{equation} +\phi_1 = \phi_0 + \omega\,\Delta t, +\end{equation} +% +the orientation quaternion $o$ is updated by the same amount by +% +\begin{equation} +o_1 = \del{\cos\del{|\omega|\,\Delta t} + +\hat\omega\sin\del{|\omega|\,\Delta t}} \cdot o_0. +\end{equation} +% +The first term is simply the unit quaternion corresponding to the +3D rotation $\omega\,\Delta t$ as in +equation~(\ref{eq-rotation-quaternion}). It is applied to the +previous value $o_0$ by multiplying the quaternion from the left. +This update is performed both during the calculation of +$k_2\ldots k_4$ and when computing the final step result. Finally, in +order to improve numerical stability, the quaternion is normalized to +unit length. + +Since most of a rocket's flight occurs in a straight line, rather +large time steps can be utilized. However, the rocket may encounter +occasional oscillation, which may affect its flight notably. +Therefore the time step utilized is dynamically reduced in cases where +the angular velocity or angular acceleration exceeds a predefined +limit. This allows utilizing reasonably large time steps for most of +the flight, while maintaining the accuracy during oscillation. + + +\subsection{Recovery simulation} + +All model rockets must have some recovery system for safe landing. +This is typically done either using a parachute or a streamer. When a +parachute is deployed the rocket typically splits in half, and it is +no longer practical to compute the orientation of the rocket. +Therefore at this point the simulation changes to a simpler, three +degree of freedom simulation, where only the position of the rocket is +computed. + +The entire drag coefficient of the rocket is assumed to come from the +deployed recovery devices. For parachutes the drag coefficient is +by default 0.8~\cite[p.~13-23]{hoerner} with the reference area being the +area of the parachute. The user can also define their own drag +coefficient. + +The drag coefficient of streamers depend on the material, width and +length of the streamer. The drag coefficient and optimization of +streamers has been an item of much intrest within the rocketry +community, with competitions being held on streamer descent time +durations~\cite{streamer-optimization}. In order to estimate the drag +coefficient of streamers, a series of experiments were perfomed using +the $40\times40\times120$~cm wind tunnel of +Pollux~\cite{pollux-wind-tunnel}. The experiments were performed +using various materials, widths and lengths of streamers and at +different wind speeds. From these results an empirical formula was +devised that estimates the drag coefficient of streamers. The +experimental results and the derivation of the empirical formula are +presented in Appendix~\ref{app-streamers}. Validation performed with +an independent set of measurements indicates that the drag coefficient +is estimated with an accuracy of about 20\%, which translates to a +descent velocity accuracy within 15\% of the true value. + + + + +\subsection{Simulation events} + +Numerous different events may cause actions to be taken during a +rocket's flight. For example in high-power rockets the burnout or +ignition charge of the first stage's motor may trigger the ignition of +a second stage motor. Similarly a flight computer may deploy a small +drogue parachute when apogee is detected and the main parachute is +deployed later at a predefined lower altitude. To accomodate +different configurations a simulation event system is used, where +events may cause other events to be triggered. + +Table~\ref{tab-simulation-events} lists the available simulation +events and which of them can be used to trigger motor ignition or recovery +device deployment. Each trigger event may additionally include a +delay time. For example, one motor may be configured to ignite at +launch and a second motor to ignite using a timer at 5 seconds after +launch. Alternatively, a short delay of 0.5--1 seconds may be used to +simulate the delay of an ejection charge igniting the upper stage +motors. + +The flight events are also stored along with the simulated flight data +for later analysis. They are also available to the simulation +listeners, described in Section~\ref{sec-listeners}, to act upon +specific conditions. + +\begin{table} +\caption{Simulation events and the actions they may trigger (motor + ignition or recovery device deployment).} +\label{tab-simulation-events} +% +\begin{center} +\begin{tabular}{ll} +Event description & Triggers \\ +\hline +Rocket launch at $t=0$ & Ignition, recovery \\ +Motor ignition & None \\ +Motor burnout & Ignition \\ +Motor ejection charge & Ignition, recovery \\ +Launch rod cleared & None \\ +Apogee detected & Recovery \\ +Change in altitude & Recovery \\ +Touchdown after flight & None \\ +Deployment of a recovery device & None \\ +End of simulation & None \\ +\hline +\end{tabular} +\end{center} +\end{table} + + + + +