﻿ Tensors of Rank Two in Tensor Flight Dynamics

Advances in Aerospace Science and Technology
Vol.03 No.02(2018), Article ID:85431,9 pages
10.4236/aast.2018.32002

Tensors of Rank Two in Tensor Flight Dynamics

Peter H. Zipfel

Modeling and Simulation Technologies, Shalimar, FL, USA    Received: April 12, 2018; Accepted: June 19, 2018; Published: June 22, 2018

ABSTRACT

Tensor flight dynamics solves flight dynamics problems using Cartesian tensors, which are invariant under coordinate transformations, rather than Gibbs’ vectors, which change under time-varying transformations. Three tensors of rank two play a prominent role and are the subject of this paper: moment of inertia, rotation, and angular velocity tensor. A new theorem is proven governing the shift of reference frames, which is used to derive the angular velocity tensor from the rotation tensor. As applications, the general strap-down INS equations are derived, and the effect of the time-rate-of-change of the moment of inertia tensor on missile dynamics is investigated.

Keywords:

Tensor Flight Dynamics, Covariance Principle, Einstein, Rotation Tensor, Angular Velocity Tensor, Moment of Inertia Tensor, Rotational Time Derivative, Euler Transformation, INS, Missile Dynamics 1. Introduction

Modeling flight dynamics with tensors enables building complex simulations by taking the two-step approach from tensor modeling to matrix programming. The new paradigm is called tensor flight dynamics. It generalizes the hitherto traditional vector flight dynamics as taught by Kane  , Likins  , Hughes  , and others. Vector analysis is based on Gibbs  , who simplified Hamilton’s  quaternions to three dimensional space, while tensor flight dynamics uses the Ricci and Levi-Civita  tensor formalism.

A tensor of rank two is in traditional dynamics a dyad, which is expressed in a particular set of orthogonal unit vectors. Special rules apply for the multiplication of a dyad with a vector  . Hughes  goes one step further and uses the concept of vectrix to express a dyad in basis vectors.

While vector flight dynamics uses only one dyad, i.e., the moment of inertia dyad, tensor flight dynamics uses three tensors of rank two, namely the moment of inertia, rotation, and angular velocity tensor. Simple tensor multiplication is used to combine them with tensors of rank one.

The moment of inertia tensor ${I}_{B}^{B}$ of body B (superscript) referred to the center of mass B (subscript) is a real, symmetrical tensor. (Tensors are expressed in bold face, upper case for rank two and lower case for rank one). When expressed in

body coordinates ${\right]}^{B}$ , the nine elements of the 3 × 3 matrix ${\left[{I}_{B}^{B}\right]}^{B}$ provide the numerical values of the axial and products of inertia.

The rotation tensor ${R}^{BA}$ describes the orientation of frame B with respect to (wrt) frame A and is an orthogonal tensor. Let B be an airplane and A the airport. Then ${R}^{BA}$ models the attitude of the airplane wrt the airport without reference to any coordinate system. If we want to calculate the attitude angles of the aircraft relative to the airport, we introduce two coordinates systems ${\right]}^{B}$ and ${\right]}^{A}$ associated with the airplane and the airport, respectively. Their transformation

matrix ${\left[T\right]}^{BA}$ , containing the attitude angles, is obtained from the rotation tensor expressed in either coordinate system ${\left[T\right]}^{BA}={\left[{\overline{R}}^{BA}\right]}^{B}={\left[{\overline{R}}^{BA}\right]}^{A}$ , where the overbar indicates the transpose of the matrix.

The angular velocity tensor ${\Omega }^{BA}$ models the angular rate of frame B wrt frame A. Because it is a skew-symmetric tensor of rank two it can be contracted to a tensor of rank one, ${\omega }^{BA}$ . (For brevity, tensors of rank two are called tensors, while tensors of rank one are called vectors). In vector analysis the angular velocity vector is called an axial vector, which has to abide by the right-hand rule, because, in actuality, it is contracted from a tensor of rank two.

Tensor flight dynamics elevates the modeling of flight dynamics to a coordinate-invariant form. It is based on Einstein’s Covariance Principle  , Page 107, which states, “All physical laws are invariant under all coordinate transformations”. While Einstein based his special and general theory of relativity on the covariance principle, tensor flight dynamics, a branch of classical mechanics, applies this principle to Newton’s laws of motion.

However, to achieve this invariance, the ordinary time derivative, must be replaced by a new time operator called the rotational time derivative (  , Section 4.2.1). The ordinary time derivative applied to a vector destroys its tensor characteristics, if the vector is subjected to a time-dependent coordinate transformation. By replacing the ordinary time operator in Newton’s Second Law by the rotational time derivative, the tensor characteristic of the law is maintained. With p the linear momentum, f the externally applied forces, and ${D}^{I}$ the rotational time derivative wrt the inertial frame I, Newton’s Second Law is expressed in tensor form

${D}^{I}p=f$

which is valid in all Cartesian coordinate systems even those that are related by time dependent coordinate transformations.

Besides the rotational time derivative, the Euler transformation is instrumental to enable tensor flight dynamics. It governs the shift of reference frames. Given frames I and B and the angular velocity tensor of frame I wrt B, ${\Omega }^{IB}$ , the rotational time derivative of tensors of rank one x transforms from frame I to frame B like

${D}^{B}x={D}^{I}x+{\Omega }^{IB}x$ (1)

This is a tensor relationship valid in all coordinate systems. Its equivalent in vector mechanics is for vector x

${\frac{\text{d}x}{\text{d}t}|}_{B}={\frac{\text{d}x}{\text{d}t}|}_{I}+{\omega }^{IB}×x$

The proof of the Euler transformation for tensors of rank one can be found in  , Appendix D. Here we present and prove the Euler transformation for tensors of rank two. It is the foundation for deriving the relationship between the rotation tensor and the angular velocity tensor, which hitherto was based on a definition, both in vector mechanics  and in tensor flight dynamics  , rather than on a strict mathematical derivation.

The Euler transformation of tensors of rank two has some precursors in vector mechanics. In  , Equation 8.93, we find for the dyad D and the coordinate systems a and b

$\frac{{}^{a}\text{d}}{\text{d}t}D=\frac{{}^{b}\text{d}}{\text{d}t}D+{\omega }^{ba}×D-D×{\omega }^{ba}$

However, because of the ordinary time derivative, this relationship is only valid for the two coordinate systems. Similarly,  , Appendix B.4, Equation 17 provides the relationship for the dyad D

$\underset{\to }{\stackrel{˙}{D}}={F}_{b}^{\text{T}}\left({\stackrel{˙}{D}}_{b}+{\omega }_{ba}^{×}{D}_{b}-{D}_{b}{\omega }_{ba}^{×}\right){F}_{b}$

where ${F}_{b}$ is the vectrix of coordinate system b. Again, the ordinary time derivative (over-dot) limits this expression to the two coordinate systems a and b.

As we will see, the Euler transformation of tensors of rank two has some similarities with these precursors from vector mechanics, but because it is valid for all coordinate systems it is a true tensor concept.

2. Euler Transformation for Tensors of Rank Two

Though the Euler transformation for tensors of rank one is sufficient for almost all flight dynamics situations, this extension to tensors of rank two makes possible a proof of an important kinematic relationship, which is given in Section 3. Additional examples are then provided in Sections 4 and 5.

We start with the formal statement of the new transformation, followed by the proof.

Let A and B be two arbitrary frames related by the angular velocity tensor ${\Omega }^{BA}$ . Then, for any tensor X of rank two, the following transformation of the rotational time derivative holds:

${D}^{A}X={D}^{B}X+{\Omega }^{BA}X+X{\overline{\Omega }}^{BA}$ (2)

The overbar indicates the transpose.

For the proof, we start with the Euler transformation for tensors of rank one, Equation (1)

${D}^{A}x={D}^{B}x+{\Omega }^{BA}x$

and substitute $x=Xy$ , where y is an arbitrary tensor of rank one

${D}^{A}\left(Xy\right)={D}^{B}\left(Xy\right)+{\Omega }^{BA}Xy$

Apply the chain rule

$\left({D}^{A}X\right)y+X{D}^{A}y=\left({D}^{B}X\right)y+X{D}^{B}y+{\Omega }^{BA}Xy$

Now the rotational time derivative of the second term on the left-hand side is transformed from frame A to frame B

$X{D}^{A}y=X{D}^{B}y+X{\Omega }^{BA}y$

and introduced into the previous equation

$\left({D}^{A}X\right)y+X{D}^{B}y+X{\Omega }^{BA}y=\left({D}^{B}X\right)y+X{D}^{B}y+{\Omega }^{BA}Xy$

The second terms on the left- and right-hand sides cancel; and since y is arbitrary, we have

${D}^{A}X+X{\Omega }^{BA}={D}^{B}X+{\Omega }^{BA}X$

Bringing the second term on the left-hand side to the right-hand side and using the fact that the angular velocity tensor is skew-symmetric, $-{\Omega }^{BA}={\overline{\Omega }}^{BA}$ , we have proven the transformation law of tensors of rank two

${D}^{A}X={D}^{B}X+{\Omega }^{BA}X+X{\overline{\Omega }}^{BA}$

3. Angular Velocity Tensor

The angular velocity tensor of rank two ${\Omega }^{BA}$ of frame B wrt to frame A is a skew-symmetric tensor and can therefore be contracted to a tensor of rank one ${\omega }^{BA}$ . Whether first or second rank, their derivations in the past were based on geometric consideration, without mathematical substantiation. Now, with the Euler transformation for tensors of rank two, we can provide such a mathematical derivation using the rotation tensor.

The rotation tensor ${R}^{BA}$ of rank two establishes the orientation of any frame B wrt frame A. The angular velocity tensor ${\Omega }^{BA}$ is then postulated (  , Equation 4.47) to be obtained by

${\Omega }^{BA}=\left({D}^{A}{R}^{BA}\right){\overline{R}}^{BA}$ (3)

To prove this relationship, we apply Equation (2) to the rotation tensor ${R}^{BA}$

${D}^{A}{R}^{BA}={D}^{B}{R}^{BA}+{\Omega }^{BA}{R}^{BA}+{R}^{BA}{\overline{\Omega }}^{BA}$

and post-multiply the right and left-hand sides by ${\overline{R}}^{BA}$

$\left({D}^{A}{R}^{BA}\right){\overline{R}}^{BA}=\left({D}^{B}{R}^{BA}\right){\overline{R}}^{BA}+{\Omega }^{BA}{R}^{BA}{\overline{R}}^{BA}+{R}^{BA}{\overline{\Omega }}^{BA}{\overline{R}}^{BA}$

The second term on the right-hand side is ${\Omega }^{BA}$ , because ${R}^{BA}$ is orthogonal and thus ${R}^{BA}{\overline{R}}^{BA}=E$ (unit tensor). Together with the left-hand term we get our relationship, Equation (3).

What remains to be shown is that

${R}^{BA}{\overline{\Omega }}^{BA}{\overline{R}}^{BA}=-\left({D}^{B}{R}^{BA}\right){\overline{R}}^{BA}$

is satisfied identically. We pre-multiply with ${\overline{R}}^{BA}$ and post-multiply with ${R}^{BA}$ . With ${\overline{R}}^{BA}{R}^{BA}={R}^{BA}{\overline{R}}^{BA}=E$ we get

${\overline{\Omega }}^{BA}=-{\overline{R}}^{BA}\left({D}^{B}{R}^{BA}\right)$ (4)

which is almost again in the form of Equation (3). To show that, we resort to a trick using the fact that the rotational time derivative of the unit tensor E is zero

${D}^{B}E={D}^{B}\left({R}^{AB}{\overline{R}}^{AB}\right)=0$

We apply the chain rule

${D}^{B}\left({R}^{AB}{\overline{R}}^{AB}\right)={D}^{B}\left({R}^{AB}\right){\overline{R}}^{AB}+{R}^{AB}{D}^{B}{\overline{R}}^{AB}=0$

and modify the last term, using ${R}^{AB}={\overline{R}}^{BA}$ (  , Section 4.1.1)

${D}^{B}\left({R}^{AB}\right){\overline{R}}^{AB}=-{R}^{AB}{D}^{B}{\overline{R}}^{AB}=-{\overline{R}}^{BA}{D}^{B}{R}^{BA}$

Substituting into Equation (4) we get

${\overline{\Omega }}^{BA}={D}^{B}\left({R}^{AB}\right){\overline{R}}^{AB}$

But because ${\overline{\Omega }}^{BA}={\Omega }^{AB}$ (  Section 4.2.4.1), we obtain again our Equation (3) with the frames A and B reversed

${\Omega }^{AB}={D}^{B}\left({R}^{AB}\right){\overline{R}}^{AB}$

Because this part of the proof is satisfied identically, we have succeeded in deriving the mathematical relationship of the angular velocity tensor ${\Omega }^{BA}$ with the rotation tensor ${R}^{BA}$ as shown in Equation (3)

${\Omega }^{BA}=\left({D}^{A}{R}^{BA}\right){\overline{R}}^{BA}$

4. Derivation of the Strap-Down INS Equations

Any kind of tactical missile today has a strap-down INS. Its sensors consist of gyros and accelerometers mounted on the missile body frame B. Integrating the accelerometer measurements twice yields the missile position, provided they are first converted to the inertial frame I with the help of the gyro measurements. We will derive these equations that need be programmed for the navigation computer. The gyros measure the angular velocity ${\Omega }^{BI}$ of the missile body B wrt to the inertial frame I. This measurement is related to the attitude of the missile by the rotation tensor ${R}^{BI}$ , as expressed by Equation (3)

${\Omega }^{BI}=\left({D}^{I}{R}^{BI}\right){\overline{R}}^{BI}$

Post-multiplying by ${R}^{BI}$ yields the time differential equations that govern the attitude of the missile, given the gyro measurements

${D}^{I}{R}^{BI}={\Omega }^{BI}{R}^{BI}$ (5)

To implement this relationship in the navigation processor we must convert the tensors to matrices by introducing coordinate systems. Though ${\Omega }^{BI}$ is measured in body coordinates, we start by applying inertial coordinates to Equation (5)

${\left[{D}^{I}{R}^{BI}\right]}^{I}={\left[{\Omega }^{BI}\right]}^{I}{\left[{R}^{BI}\right]}^{I}$

The rotational derivative wrt to frame I becomes the ordinary time derivative, as the coordinate system ${\right]}^{I}$ is associated (fixed) with frame I

${\left[\frac{\text{d}}{\text{d}t}{R}^{BI}\right]}^{I}={\left[{\Omega }^{BI}\right]}^{I}{\left[{R}^{BI}\right]}^{I}$

Because the gyro measurements are made in missile body coordinates,

${\left[{\Omega }^{BI}\right]}^{B}$ , we need to introduce the coordinate transformation matrix ${\left[T\right]}^{BI}$ of body coordinates ${\right]}^{B}$ wrt inertial coordinates ${\right]}^{I}$ , ${\left[{\Omega }^{BI}\right]}^{I}={\left[\overline{T}\right]}^{BI}{\left[{\Omega }^{BI}\right]}^{B}{\left[T\right]}^{BI}$

${\left[\frac{\text{d}}{\text{d}t}{R}^{BI}\right]}^{I}={\left[\overline{T}\right]}^{BI}{\left[{\Omega }^{BI}\right]}^{B}{\left[T\right]}^{BI}{\left[{R}^{BI}\right]}^{I}$

Now, because the coordinate systems ${\right]}^{B}$ and ${\right]}^{I}$ are associated with frames B and I, respectively, there exists a relationship between the rotation tensor

${\left[{R}^{BI}\right]}^{I}$ and the transformation matrix ${\left[T\right]}^{BI}$ , namely ${\left[{R}^{BI}\right]}^{I}={\left[\overline{T}\right]}^{BI}$ (  , Section 4.1.1), which we now introduce

${\overline{\left[\frac{\text{d}}{\text{d}t}T\right]}}^{BI}={\left[\overline{T}\right]}^{BI}{\left[{\Omega }^{BI}\right]}^{B}{\left[T\right]}^{BI}{\left[\overline{T}\right]}^{BI}$

with ${\left[T\right]}^{BI}{\left[\overline{T}\right]}^{BI}=\left[E\right]$ and ${\left[\overline{T}\right]}^{BI}={\left[T\right]}^{IB}$ (  , Section 3.2.1.3)

${\left[\frac{\text{d}}{\text{d}t}T\right]}^{IB}={\left[T\right]}^{IB}{\left[{\Omega }^{BI}\right]}^{B}$ (6)

here we have the nine linear differential equations to be programmed for the navigation processor. The transfer alignment will provide the initial conditions.

Applying ${\left[T\right]}^{IB}$ to the measured accelerations ${\left[{a}_{B}^{I}\right]}^{B}$

${\left[{a}_{B}^{I}\right]}^{I}={\left[T\right]}^{IB}{\left[{a}_{B}^{I}\right]}^{B}$

yields the acceleration in inertial coordinates, ready to be integrated for the missile’s velocity and position after a gravity compensation.

5. Effect of Variable Moment of Inertia on Missile Dynamics

It is common practice in flight dynamics to neglect the time-rate-of-change of mass except in rocketry, where the expulsion of fuel is called thrust, but considered an external force. The same simplification is made in the attitude equations of motion, where the time-rate-of-change of the moment of inertia tensor is neglected.

Thomson  discusses this effect as early as 1966. He derives the equations of motion with vector notation, and uses a fixed body point as reference. Quadrelli  addresses the variable mass effects of a supersonic inflatable decelerator using a coordinate-free approach based on Reynolds’ transport theorem from continuum mechanics. Here we use our tensor formulation to investigate the effect of the varying moment of inertia on a short-range air-to-air missile and refer the equations of motion to the center of mass.

The attitude equations of motion are governed by Euler’s law: The inertial time-rate-of-change ${D}^{I}\left(*\right)$ of the angular momentum ${l}_{B}^{BI}$ of the body B wrt the inertial frame I, referred to the c.m. B equals the eternal moments ${m}_{B}$ applied at the c.m. B

${D}^{I}\left({l}_{B}^{BI}\right)={m}_{B}$

where the angular momentum is the product of the moment of inertia tensor ${I}_{B}^{B}$ of body B referred to the c.m. B and the angular velocity ${\omega }^{BI}$ of body B wrt the inertial frame I

${l}_{B}^{BI}={I}_{B}^{B}{\omega }^{BI}$

Substitution into the previous equation yields

${D}^{I}\left({I}_{B}^{B}{\omega }^{BI}\right)={m}_{B}$

here the moment of inertia ${I}_{B}^{B}$ is a tensor of rank two, ${\omega }^{BI}$ and ${m}_{B}$ are tensors of rank one. We apply the Euler transformation for tensors of rank one, Equation (1), in order to shift the rotational time derivative to frame B

${D}^{B}\left({I}_{B}^{B}{\omega }^{BI}\right)+{\Omega }^{BI}{I}_{B}^{B}{\omega }^{BI}={m}_{B}$

and apply the chain rule to the argument of the rotational time derivative

${D}^{B}\left({I}_{B}^{B}\right){\omega }^{BI}+{I}_{B}^{B}{D}^{B}\left({\omega }^{BI}\right)+{\Omega }^{BI}{I}_{B}^{B}{\omega }^{BI}={m}_{B}$ (7)

In flight dynamics it is usually assumed that ${D}^{B}\left({I}_{B}^{B}\right)$ is negligibly small, so

that the attitude equations of motion assume the simpler form

${I}_{B}^{B}{D}^{B}{\omega }^{BI}+{\Omega }^{BI}{I}_{B}^{B}{\omega }^{BI}={m}_{B}$

This assumption may be acceptable for aircraft. But is it also justified for missiles, particularly for air-to-air missiles, whose fuel may be 40% of the launch mass and expended in less than 3 seconds? To pursue this investigation, we use a typical short range air-to-air missile, called AIM6  , and use the complete attitude equations of motion, Equation (7). Introducing body coordinates ${\right]}^{B}$ and multiplying out the matrices to obtain the scalar equations in the state variables

${\left[{\overline{\omega }}^{BI}\right]}^{B}=\left[\begin{array}{ccc}p& q& r\end{array}\right]$ , yields the attitude equations of motion

$\begin{array}{l}\stackrel{˙}{p}=\left({m}_{B1}-{\stackrel{˙}{I}}_{1}p\right)/{I}_{1}\\ \stackrel{˙}{q}=\left\{\left({I}_{2}-{I}_{1}\right)pr+{m}_{B2}-{\stackrel{˙}{I}}_{2}q\right\}/{I}_{2}\\ \stackrel{˙}{r}=\left\{-\left({I}_{2}-{I}_{1}\right)pq+{m}_{B3}-{\stackrel{˙}{I}}_{2}r\right\}/{I}_{2}\end{array}$

where the over-dot indicates the ordinary time derivative, and where the moment of inertia for a missile with tetragonal symmetry ( ${I}_{3}={I}_{2}$ ) is given in body coordinates

${\left[{I}_{B}^{B}\right]}^{B}=\left[\begin{array}{ccc}{I}_{1}& 0& 0\\ 0& {I}_{2}& 0\\ 0& 0& {I}_{2}\end{array}\right]$

with the external moments ${\left[{\overline{m}}_{B}\right]}^{B}=\left[\begin{array}{ccc}{m}_{B1}& {m}_{B2}& {m}_{B3}\end{array}\right]$ .

It appears that the changing moment of inertia terms introduce a lag into the three channels. Let’s assume the missile executes a pure yaw maneuver, then $p=q=\stackrel{˙}{p}=\stackrel{˙}{q}=0$ and what remains is

$\stackrel{˙}{r}+\frac{{\stackrel{˙}{I}}_{2}}{{I}_{2}}r=\frac{{m}_{B3}}{{I}_{2}}$

This is a simple first order lag system with the time constant $T={I}_{2}/{\stackrel{˙}{I}}_{2}$ . However, because ${\stackrel{˙}{I}}_{2}$ is negative, the time constant is negative and we are faced with destabilization. How serious is this effect?

I ran a horizontal engagement of the prototype missile AIM6 using its 6 DoF simulation  and recorded ${I}_{2},{\stackrel{˙}{I}}_{2}$ at 1.345 sec, Mach=2, and lateral acceleration 33 g’s. The time constant at that point is $T=51.195/-5.745=-8.911$ sec. Since the thrust phase lasts only 2.5 sec, the destabilizing effect should not be too severe.

To see the actual evidence, I used the 6 DoF simulation without autopilot, which is possible because the missile is aerodynamically stable. Figure 1 displays the results of two runs with and without the unsteady moment of inertia, taken around the time of our test point at 1.345 sec.

The unsteady effect lowers the pitch angle and increases the amplitude of the pitch rate. The discrepancies are small and can be neglected. With the autopilot engaged they are not noticeable.

6. Conclusion

With the proof of the Euler transformation for tensors of rank two, all elements of tensor flight dynamics are now in place. This transformation provides us with a mathematically sound derivation of the relationship between angular velocity

Figure 1. Effect of the time-rate-of-change of moment of inertia.

and rotation tensor, which we used for the straight-forward derivation of the navigation equations of a strap-down INS. By applying the transformation to the moment of inertia tensor, we could discuss how its rapidly changing values affect the dynamics of an air-to-air missile. Because the effect is negligible, the common practice of neglecting the time-rate-of-change of moment of inertia is justified.

Cite this paper

Zipfel, P.H. (2018) Tensors of Rank Two in Tensor Flight Dynamics. Advances in Aerospace Science and Technology, 3, 11-19. https://doi.org/10.4236/aast.2018.32002

References

1. 1. Kane, T.R. (1968) Dynamics. Holt Rinchart and Winston, New York.

2. 2. Likins, P.W. (1974) Elements of Engineering Mechanics. McGraw-Hill, New York.

3. 3. Hughes, P.C. (2004) Spacecraft Attitude Dynamics, Dover. John Wiley & Sons, New York.

4. 4. Gibbs, W.J. (1961) The Scientific Papers of J Willard Gibbs, Vol II Dynamics. Dover, New York.

5. 5. Hamilton, W.R. (1843) On a New Species of Imaginary Quantities Connected with a Theory of Quaternions. Dublin Proceedings, 2, 424-434.

6. 6. Ricci, G. and Levi-Civita, T. (1901) Methodes de Calcul Differentiel Absolu et Leurs Applications. Mathematische Annalen, 54.

7. 7. Adler, R., Bazin, M. and Schiffer, M. (1965) Introduction to General Relativity. McGraw Hill, New York.

8. 8. Zipfel, P.H. (2014) Modeling and Simulation of Aerospace Vehicle Dynamics. 3rd Edition, AIAA Education Series, AIAA, Washington DC.

9. 9. Thompson, W.T. (1966) Equations of Motion for the Variable Mass System. AIAA Journal, 4, 766-768. https://doi.org/10.2514/3.3544

10. 10. Quadrelli, M.B., Cameron, J., Balaram, B., Branwal, M. and Bruno, A. (2014) Modeling and Simulation of Flight Dynamics of Variable Mass Systems. AIAA/AAS Aerodynamic Specialist Conference, San Diego, 4-7 August 2014.

11. 11. Zipfel, P.H. (2014) Fundamentals of Six Degrees of Freedom Aerospace Simulation and Analysis in C++. 2nd Edition, AIAA, Washington DC.