Category Modeling and Simulation of Aerospace Vehicle Dynamics

10.1.2.3 Finally, we solve the basic kinematic problem, namely, given the inertial body rates in body coordinates calculate the direction cosine matrix [T]BI. We have dealt with this issue before in Chapter 4 and presented three possible approaches. For the flat Earth, six – DoF simulations of the preceding section, I recommended the quaternion solution for two reasons: 1) only four linear differential equations need be solved and 2) the Euler angles are readily computed from the direction cosine matrix [T]BL. Here, the question arises what are the Euler angles associated with [T]B/? They would be related to the inertial coordinate system and, therefore, void of any flight mechanical meaning.

Hence, let us calculate [T]BI directly from ojb/b, using the direction cosine differential equations (4.66)

– j T " BI

— = [0"]В[Г]В/ (10.47)

. df.

where [£2B/]B is the skew-symmetric form of [coBI]B, provided by Euler’s equation.

The initialization of [T]Bl is of immediate interest. Most likely, the simulation is initialized by the standard Euler angles, referred to as the geodetic coordinate system. Accordingly, the transformation matrix [T]BD can be calculated. Yet, we also need to compute [T]DI to complete the initialization process of [T]BI = [T BD[TDI. This is easily done: Given the initial geodetic longitude, latitude, and altitude, [T]DI is obtained from Eq. (10.39).

Once the vehicle is airborne and [T]B1 is being updated, we have to solve the inverse problem, namely, what are the Euler angles relative to geodetic coordinates? “No problem,” you say, let us just calculate [T][B from the current longitude, latitude, and altitude and solve for [T]BD = [T]B,[TDI. The Euler angles are then

9 — arcsin(—Гп)

where the ty are the elements of [T}BD. Correct, just make sure to program around cos 9 — 0 to avoid unwanted divide by zero run-time errors.

All would be well if we had a perfect computer to solve the differential equations. However, finite integration steps and limited word length cause the orthogonality of the direction cosine matrix to be slowly lost. Consequently, we have to reestablish orthogonality after every integration step. The following procedure, patterned after Savage,8 works quite well.

Simply orthogonalize [T]BI by the following recursive algorithm

[T(n + 1)]B/ = [T(n)]BI + {[E] – [T{n)]BI[T{n)]BI}[T(n)}Bl (10.49)

where [T(n +1 )]fl/ is the new orthogonalized matrix, based on the current direction cosine matrix [T(n)]BI attained by integration. Notice, the value in the parentheses is a zero matrix if [T(n)]BI is orthogonal. To show that indeed the orthogonality has improved, I invite you on a side tour, which you are free to bypass.

Let us expand the last term

[ST]BI = {[E] – [T(n)]B,[T(n)]BI}[T(n)]BI using the base vector format of the transformation matrix [see Eq. (3.5)]

where bi is the (3×1) column base vector and b, its transpose, a (1 x 3) row vector,

-bxb2b2 – Ьф3ь3 – b2bb – Ъ2ь3Ъ3 – ЬзЬфі – b3b2b2

Substituting into the recursive relationship and writing out the lines yields b(n + 1) = bi(n) – [b(n)b2(n)b2(n) + bx{n)b3(n)b3{n)}

b2(n + 1) = b2(n) – [b2(n)bx{n)b(n) + b2{n)b3{n)b3{n)]
b3{n + 1) = b3(n) – [b3{n)bi(n)bi{n) + b3(n)b2(n)b2(n)]

After these preliminaries we want to check if the orthogonality error decreases as a result of the recursive algorithm. Before the update the error of any two base vectors has the scalar value £,-* = Ьх(п)Ьфп). Let us say this error is of order О(є). Now we want to show that after the update the error Б,-(n + 1 jb^in + 1) has decreased by an order of magnitude to 0(e2). Multiply any two base vectors

bi(n + l)b2(n + 1) = b(n)b2(n) – bx{n)bx{n)bx{n)b2(n)

-bx{n)bx{n)bx{n)b2(n) + 0(є2) + 0(є3)

and recognize that the first three terms on the right-hand side cancel. Hence, the orthogonality error of the updated base vectors was reduced to 0(e2). Because this conclusion holds for base vectors one and two, it is true for any two combinations and thus also for the direction cosine matrix [T]BI.

You can use Eqs. (10.47-10.49) to solve the basic kinematic problem in your simulation. They will deliver the direction cosine matrix [ T]BI and the Euler angles xjr, в, and ф wrt the geodetic coordinates. Besides the inertial body rates [со81]8, you will also need the geodetic longitude, latitude, and altitude of the vehicle. For an example you may refer to CAD AC GHAME6 simulation, Module G3.

In summary, Fig. 10.6 combines the most important parts of a six-DoF simu­lation for a hypersonic vehicle over an elliptical Earth. The CADAC GHAME6 model is patterned after this diagram. Newton’s equations are given by Eqs. (10.42) and (10.43) and Euler’s equations by Eq. (10.44). The computation of the direc­tion cosine matrix follows Eqs. (10.47) and (10.48). For the representation of the aerodynamic forces, I chose the lift and drag coefficients, defined in stability coordinates.

The initialization of the state variables is also indicated. The geographic velocity is given in geodetic variables, speed |uf |, heading /, flight-path angle у, geodetic longitude l, latitude Xj, and altitude h. They are converted to the inertial velocity vector [ng]7 and the vehicle’s initial position [ід/]1. If not zero, the initial geo­graphic angular rate [coBE]D, given in geodetic coordinates, is converted to inertial

body rates [&>B/]/. Finally, the direction cosine matrix is initialized by the vehicle’s Euler angles.

10.1.2.2 In the pre­ceding Hat-Earth approach we had the choice of modeling the equations of mo­tion either in local-level, presumed inertial, or body coordinates. Now we have five possibilities: inertial, Earth, geographic, geodetic, or body coordinates. As bewildering as these choices can be, surveying the literature I find these op­tions: ENDOSIM6 integrates Newton’s equation in inertial coordinates, Euler’s equation in body coordinates and uses quaternions, requiring the introduction of inertial Euler angles. Another U. S. Army simulation7 solves Newton’s equa­tion in geocentric coordinates, while Euler’s equation is again given in body coordinates.

There is little doubt that Euler’s equation should be integrated in body coordi­nates because of the inherent invariance of the MOI tensor. However, for Newton’s equation we have a choice. I prefer the inertial approach, integrating both the ac­celeration and velocity equations in inertial coordinates. The other option is the integration of the geographic acceleration and velocity. We will describe it later. We also have to decide whether we use quaternions or the direction cosine matrix. Because quaternions are closely related to Euler angles, the introduction of inertial Euler angles becomes necessary. I propose that we use the direction cosine matrix and avoid confusing angle definitions.

The derivation of the equations of motion is easy and short. Newton’s law provides

mDlvlB – fa, p +mg

The left side is best expressed in inertial coordinates, but on the right-hand side the forces are formulated in body coordinates and the gravitational acceleration in geocentric coordinates:

m[D4]7 = [T]Br[fa, p]B +m[T]G,[g]G

This leads to the first set of differential equations:

(10.42)

The (T]BI transformation is the direction cosine matrix, and the other transforma­tion [f]GI is obtained from Eqs. (10.36) and (10.39):

Given the inertial velocity [vBY, the inertial displacement vector of the vehicle [яд/]7 is obtained from the second set of differential equations:

= КҐ (10-43)

Equations (10.42) and (10.43) are the translational differential equations, which couple with the rotational equations through the direction cosine matrix [T]BI and the forces [fa, pB ■ A precarious issue is the initialization of the state variables. We will deal with it after we have derived the attitude equation.

The attitude equation is obtained as before from Euler’s law

D’lg = mB

Express the angular momentum as the MOI tensor and the inertial body rates lBg = Ibblobi transform the rotational derivative to the body frame; and finally express the tensors in body axes

5 = + -тв]В) (10.44)

This differential equation is just like the one for the flat-Earth case [Eq. (10.10)], except that the Earth frame E is replaced by the inertial frame I.

Now I must address the initialization of the state variables. So far, we attained nine first-order differential equations, and therefore nine state variables must be initialized. They are divided into three pairs of three each: WBY from Newton’s equation (10.42), [.sg/J’ from the velocity integration Eq. (10.43), and [сошв from Euler’s equation (10.44). The initialization of the direction cosine matrix will be addressed in the next section.

When you build a simulation, you should make it user friendly. Particularly the input should be easy to visualize. Therefore, you want to replace WbY by the geographic velocity vector [uf ]D, expressed in geodetic coordinates. We generate such a relationship by first reviewing the definition of these velocities vB = D1 sBi and vf = DesBe, where sBi — sBe+ sei (E is any fixed point on the Earth and I is Earth’s center), and therefore

D1sBi = DesBe + Desei + EIeisB;

with Desei equals zero because sei is constant in the Earth frame. Thus

VB =VEB+ EIbesBi (10.45)

and expressed in the appropriate coordinate systems

КҐ = [t]D,[vEB]D + (10.46)

Furthermore, because [T]DI is only a function of longitude and latitude, we can summarize: Given the geodetic longitude, latitude, and altitude, and the geographic
velocity vector [i>f ]D, the inertial velocity vector [u^J7 can be initialized. A con­venient alternative to [i>|]D is the initialization by the geographic speed |nf | and the heading angle / and the flight-path angle у, both relative to the local tangent plane. All you need is a conversion routine as, for instance, the CAD AC subroutine MATCAR.

As far as the angular velocity [шві]в is concerned, it is usually initialized by zero values. Only in extreme cases, like close-in air-to-air combat, may the missile have to be given the angular rates of the launch aircraft. Then the aircraft’s initial geographic angular rate [coBE]D has to be converted to the inertial body rate in body coordinates

[шв1]в = [T]BD[coBE]D + [T]BI[coE,}!

where [coEI]! is the Earth’s rotation.

Now, there is an alternate way to initialize the inertial velocity. It is based on Euler and incidence angles. From these angles and the geographic speed the geographic velocity vector can be derived. Express the geographic velocity vector in body and relative wind coordinates ]w

[vE]B = [T]WB[vEB]W = [f]WB

With the transformation matrix Eq. (3.18) introduce the incidence angles

cos a cos /3
sin/6

sin a cos /6

and transform the velocity to geodetic coordinates

[vEB]D = mBD[v ir

The transformation matrix [T]BD contains the Euler angles referred to the geodetic coordinate axes and follows the format of Eq. (3.14). Once [nf ]D is computed, the determination of follows the preceding outline.

As we travel into the higher regions of the atmosphere, or leave it altogether for the wide expanse of space, our speed gets faster and our view of the globe broader. It becomes evident that the Earth is round and rotating with respect to the sun. As travelers, we experience such effects as centrifugal and Coriolis
accelerations and are looking for another immutable reference frame somewhere out there.

This situation motivates us to study the shape of the Earth in more detail. First, I shall give a brief overview of geodesy, the geodetic coordinate system, and the gravitational attraction of the Earth. Then we are prepared to formulate the six-DoF equations of motions in an inertial and geographic format. For the kinematic equa­tions we shall pursue the direction cosine approach rather than use quaternions.

Some of this material should be familiar to you from Chapters 8 and 9. Here, however, we expand our horizon to an elliptical Earth and deal with the mechanics of developing a simulation model for hypersonic vehicles.

10.1.2.1 Geodesy. Geodesy can be defined as the science concerned with the precise positioning of points on the surface of Earth and the determination of the exact size and shape of Earth. It also involves the study of the variations of Earth’s gravitational attraction and the application of these variations to the precise measurements of Earth’s shape.

The science of geodesy is as old as the early Greeks. Whereas Homer lived on a flat Earth, Pythagoras, the mathematician, envisioned the Earth to be fashioned according to a perfect sphere. Of the same opinion was Aristotle, the dominant scientific authority until the Middle Ages.

One of the great pursuits of the ages—and continuing today—is the accurate measurement of the Earth’s circumference. It amazes me to read that Eratosthenes, a Greek scholar experimenting in Egypt, deducted from measurements of sun an­gles that the circumference at the equator is 25,000 miles. Today we peg it at 24,899 miles. However, when Ptolemy came on the scene, he convinced the majority of the scholars that the number should be revised downward to 18,000 miles. It still stood at that value when Columbus made his plans to sail to India. Fortunately so, because otherwise Columbus may have been unwilling to travel such a great distance and America might have been discovered by Russian cosmonauts peeking down from space!

The Pythagorean spherical concept offers a simple surface, which is mathemati­cally easy to deal with. Many astronomical and navigational computations use it to represent Earth. Although the sphere is a close approximation to the true figure of Earth and satisfactory for many purposes, for the geodesist interested in measuring long distances—spanning continents and oceans—a more exact configuration is necessary.

Because Earth is in fact slightly flattened at the poles and bulges somewhat at the equator, the geometrical figure used in geodesy is an ellipsoid of revolution or spheroid. Two dimensions define a spheroid; the semimajor axis and the flattening parameter. The flattening parameter is defined as the fractional change between the semimajor axis a and semiminor axes b:

Some references also call it the ellipticity of the spheroid and invite some confusion with the term eccentricity, which is defined as

a2-b2

Thanks to the World Geodetic Systems (WGS) Committee, chartered in 1960, the numerical precision of the reference ellipsoid has been refined ever since. Today, we use the WGS 84 (Ref. 3) system with the values a = 6,378,137.0 m and e = 0.08181919.

For completeness’ sake, however, let me point out that the real shape of Earth is arbitrary and is only successively approximated by an alternating series of el­liptical and pear-shaped surfaces. So far, our classification is based on geometric considerations only (geometrical geodesy). The WGS 84 ellipsoid is that first ap­proximation beyond the sphere that fits the actual surface of the Earth best (at sea level). Alternately, the actual shape of the Earth is also referred to as geoid. It is defined as that surface to which the oceans would conform over the entire Earth if free to adjust to the combined effect of Earth’s mass attraction and the centrifugal force of Earth’s rotation. Notice the geoid is defined with the aid of physical forces. Therefore, this branch of geodesy is also called physical geodesy.

Do not be intimidated by geodesy. You will never learn all of the tricks of its trade unless you make it your life’s passion. If that is the case, I recommend you start with the book by Torge.4 Fortunately, for six-DoF simulations we are satisfied with the WGS 84 ellipsoidal approximation and do not distinguish between geometrical and physical geodesy.

Briefly review Secs. 3.1 and 3.2. There we define the inertial frame I and the associated inertial coordinate system ]7. The inertial frame’s attitude is fixed in the ecliptic, but the center of Earth is also one of its points. In effect, the inertial frame convects with Earth’s motion around the sun but does not follow Earth’s rotation. The inertial coordinate axes are aligned with the vernal equinox and Earth’s spin axis.

Earth’s frame E, the associated coordinate system ]£, and the geographic coor­dinate system ]G are also defined in Sec. 3.2. The transformation matrix [T]G£ is a function of longitude and latitude [Eq. (3.13)]. Because in Chapter 3 Earth was presumed to be spherical, latitude was measured from the equator to the geocentric displacement vector of the vehicle. However, for an elliptical Earth model we now have to refine the definition of latitude, while longitude remains the same.

We call the latitude of the WGS 84 model the geodetic latitude kj, and the earlier definition is relabeled the geocentric latitude kc. Figure 10.4 depicts the geometrical situation. To define the geodetic latitude, project the vehicle c. m. В onto the local tangent plane at Bq and continue the plumb line until it crosses the

plane through the equator. The angle between the plumb line and the equatorial plane is the geodetic latitude kd. This is the latitude that is used in cartography. If you want to go back to a round Earth, simply use the geocentric latitude kc.

Let us review and revise carefully the definition of the geographic coordinate system. We rename the former geographic system and call it the geocentric coor­dinate system while maintaining the same symbol ]G. Its 1G and 2G axes are still parallel to the plane normal to sbi and point north and east, respectively, while the 3G axis is parallel and points in the opposite direction of sBi­ts. new entity, the geodetic coordinate system, is required. Its lD and 2° axes are parallel to the plane normal to the plumb line and point north and east, respectively. The 3d axis is parallel to the plumb line, pointing downward. The geodetic and geocentric coordinate systems are related by the deflection angle <5, as indicated in Fig. 10.4, and defined by

8 = kd — kc

Its value depends on the geodetic latitude kd, the flattening parameter /, and the altitude of the vehicle h. Britting5 has given an approximation

5 = / sin(2krf) (l – L – (10.34)

where Ro is the distance between the points Bo and I, which can be calculated from the Earth’s major axis a, its flattening parameter, and the geodetic latitude of the vehicle

The accuracy of Ro is better than one meter, and 8 is within 4.5 arcsec at 30 km alti­tude. The transformation matrix of the geodetic wrt geocentric coordinate systems is governed by this deflection angle 5

(10.36)

For spherical Earth 5 = 0, the geocentric and geodetic coordinate systems coincide.

Given the longitude l (now also called geodetic longitude) and the geodetic latitude kd, the transformation matrix of the geodetic wrt the Earth coordinates is adapted from Eq. (3.13)

Now regard Fig. 10.5. We can extend Eq. (10.37) to inertial coordinates by re­placing the geodetic longitude with the celestial longitude /,, which is measured from the 11 axis. The celestial longitude is composed of the celestial longitude of Greenwich la and the geodetic longitude of the vehicle /, = lG + l. However,

because the Greenwich meridian rotates with the Earth rate co$ during the fly-out time of the vehicle, it increases from its initial location lGo by co$(t — to) becoming lG — lc0 + соф(г — t0). The celestial longitude of the vehicle at time t is then

U — /go + &>®(t — to) + l (10.38)

Substituting /, for / in Eq. (10.37) yields the transformation matrix of the geodetic wrt the inertial coordinate systems

Both transformation matrices will be useful when we build our six-DoF simula­tions.

Another important kinematic problem is the determination of the inertial posi­tion of the vehicle [лД/]/ in inertial coordinates, given the vehicle’s longitude /, latitude kd, altitude h, and initial Greenwich celestial longitude lGo – We start with the displacement vector of the vehicle in geocentric coordinates [sg)]G = [0 0 — (Ro +h), using the approximation Ro+h % , which is according to Bribing5

less than 1 m for endo-atmospheric flight. Ro is given by Eq. (10.35). Then we transform from geocentric to inertial coordinates and obtain the result

M7 = [f]DI[T]DG[sBI]G (10.40)

where the two transformation matrices are given by Eqs. (10.39) and (10.36), respectively.

The inverse problem is the determination of the vehicle’s longitude, latitude, and altitude, given the inertial displacement vector [.?£/]7. Unfortunately, a straight calculation is not possible because the deflection angle itself is a function of lati­tude. Therefore, we take recourse to an iterative approach that proceeds as follows: 1

2) Start an iterative loop in which you calculate the radius to the ellipse Ro with Eq. (10.35), the altitude h from the approximation h = |sB/| — Rq, and the deflection <5 using Eq. (10.34).

3) Now update the approximation кц Xj + 8 and stop the iterative loop if the improvement drops below a certain value, say 10-7 radians.

4) Finally, calculate the geodetic longitude

/ = arcsin {{sudi/l[(sm)[f + [Cv?/)^]2} – ho ~ coe(t – t0)

In most situations it takes fewer than 10 iterations for convergence. Yet, it is always advisable that you guard against excessive numbers of iterations.

Let us wrap up our adventure in geodesy with an excursion to the modeling of the gravitational field and the resulting gravitational acceleration that acts on the flight vehicle. As we have discussed in Sec. 8.2.2, the force of attraction F between two bodies M and m, separated by the r distance, is according to Newton

F = G-

where G is the universal gravitational constant. If m represents the flight vehicle and M the Earth, then the acceleration on the vehicle is

_F _GM

^ m r2

If the Earth were a perfect sphere of uniform density, g would sweep out a grav­itational field consisting of concentric spheres. Yet, Earth being a geoid, we have to call the mathematicians for help. They invented the spherical harmonics that successively improve the representation of the gravitational field. The first term represents the sphere, and, together with the second term, the desired model of an elliptical gravitational field is reached.

The WGS folks have worked since 1960 to improve this elliptical gravitational model. In 1984 they determined that the value of GM is 3.986005 x 1014 m3/s2 and that the second-degree zonal gravitational coefficient is C2.0 = —4.841668 x 10-4. The acceleration acting on a vehicle caused by the gravitational reference ellipsoid is given in geocentric coordinates as

with all variables already defined. The main contribution is in the 3G direction. If you neglect the ellipticity, i. e., Сг, о = 0, you recover the spherical value of g = GM/sbi2- Notice, because of symmetry, there is no gravitational component in the east direction.

Hitherto, we have labored through the preliminaries of modeling an elliptical Earth. You will use most of these equations to build the utility routines for your six-DoF simulation. In particular, you will have to program a subroutine that
calculates [T]GI, given /, Xd, h (called CADTGI84 in CADAC); another one that determines l, Xd, h, given [лд/]7 (called CADGE084 in CADAC); and the inverse that provides [ід/]7, given /, Xd, h (called CADINE84 in CADAC). You find these subroutines in CADAC utility file UTL3.FOR on the CADAC CD. Finally, the gravitational acceleration is given by Eq. (10.41) (in the CADAC GHAME6 simulation, Module G2, subroutine G2GRAV). Now we are prepared to derive the equations of motions.

10.1.1.4 If you liked the derivation of the spinning missile equations, you will really enjoy this excursion into the dynamics of bodies that spin about a horizontal axis. They are my favorite subject because I spent five years studying their characteristics, culminating in my dissertation.2

Spinning bodies were first investigated by Heinrich Gustav Magnus in 1852. He demonstrated experimentally that a body rotating in an airstream experiences a force that acts substantially normal to the airflow. An autorotating flight vehicle, designed to develop this Magnus force efficiently and to employ it as the major lift force in free flight, is called a Magnus rotor. It consists of a center body with driving vanes and endplates.

Figure 10.3 depicts a Magnus rotor spinning about a near horizontal axis and gliding with the velocity vf. The spin axes 2s and the velocity vector form a plane containing the roll axis Is normal to the spin axis. The yaw axis 3L is parallel to the gravitational acceleration.

The dynamics of Magnus rotors can best be explained in terms of a horizontally spinning gyroscope that is subjected to aerodynamic forces. The dynamic motions of the glide trajectory are rolling <j>, yawing ф, and sideslipping /8. They exhibit three modes: nutation, precession, and undulation. The nutation and precession
modes are those of a gyroscope with the modification that the nutation is aerody – namically dampened. New is the undulation mode. It arises from the interaction of aerodynamic moments—generated by /3—and gyroscopic precession. In general, this mode can be oscillatory or aperiodic.

You can understand the fun I had working on this project, developing the theory, building the computer simulation, and flight testing some specimens. For our current purpose we limit the discussion to the derivation of the equations of motion, the starting point of any six-DoF simulation of Magnus rotors.

As always, the translational equation is based on Newton’s second law

mDEv I = fa + fg

and the attitude equation derives from Euler’s law

DElf = mB

with aerodynamic forces and moments fa, mв and the gravitational force fg = mg.

Just as I used a nonspinning frame B’ for the spinning missile, in the same vein I introduce here a frame that is attached to the Magnus rotor but does not follow its spin. It is called the stability frame and is defined by the spin axis and the projection of the velocity vector on the symmetry plane of the rotor. The associated stability coordinate system has its 2s axis aligned and in the direction of the spin vector, its Is axes parallel and in the direction of the projection of the velocity vector on the symmetry plane. The third axis 3s completes the triad. The transformation matrix [T]5L of the stability wrt the local-level axes follows the conventional transformation of Euler angles: yaw j/sL, pitch 6sl. and roll ф%і (see Sec. 3.2.2.4).

The first step is to transform Newton’s equation to the stability frame and express the terms in stability coordinates

(10.25)

Focus on the meaning of [f25£]s. It is the skew-symmetric form of the angular velocity vector of the stability frame wrt Earth, which embodies the rolling and yawing rates, but excludes the spinning motion. The choice of the stability coor­dinates is dictated by the method of the aerodynamic force measurement. In the wind tunnel the sting supporting the model is centered at the spin axis, allowing the Magnus rotor to rotate freely. The force measurements, taken by the sting, are in stability axes. By the way, the term stability axes is borrowed from airplane dy­namics. There and here, the I s axis is also defined as the projection of the velocity vector on the symmetry plane of the aircraft.

From the velocity vector [uf ]s and the transformation matrix T]SL. we can get by integration the position of the Magnus rotor [_sbel, i. e., the displacement of the c. m. В with respect to an Earth reference point E, expressed in local-level coordinates.

Equations (10.25) and (10.26) are the translational equations of Magnus rotors. The angular velocity [f2S£]fl and the transformation matrix [ 7" ] s/ are still to be supplied by the rotational equations.

Similar to spinning missiles, the more exciting part is the derivation of the attitude equation. Applying Euler’s law and transforming the rotational derivative to the stability frame, we obtain

Dslf + nSElf – =mB

We divide the angular velocity шВЕ of the missile body В wrt the reference frame E into two parts:

, ,BE __ ,BS I . ,SE

UJ = LO — UJ

The first part is the spin rate, and the second one is the angular velocity of the nonspinning stability frame wrt Earth. Similarly, the angular momentum is split into two parts

If = IBBUJBE = IBBuBS + 1bujse (10.28)

We form the rotational derivative as required by Eq. (10.27)

Dslf = Ds(lBB u>BS) + Ds(Ibujse)

= IBBDsuBS + IbbDsuse (10.29)

where we made use of DSIB = 0 because the Magnus rotor’s MOI about the spin axis is independent of spin attitude. This simplifying assumption is true at least approximately because the effect of the driving vanes is insignificant. Substituting Eqs. (10.28) and (10.29) into Eq. (10.27) yields the atitude equation of motion

IbDsujbs + IbbDsuse + nSEIBuBS + tlSEIBBuSE = mB (10.30)

The total rate of change of angular momentum consists of the spin contribution— first term—and the stability frame rotations—second term. The primary gyroscopic coupling stems from the third term, which dominates the fourth term.

Just as the stability coordinates were used for the translational equations, they are also best suited for the attitude equation. The aerodynamic moment [mj]s is given in that form, the angular velocity [ws£]s is needed by the translational equations, and the rotational time derivatives become the ordinary derivatives. For these reasons we express the rotational equations in stability coordinates

+ inSE]S[lS]S([‘Js]s+ [**]*) = [mB]s

(10.31)

Another important fact is that the MOI tensor is constant in stability coordinates and has a particularly simple diagonal form (note because of symmetry /3 = /j)

with the angular velocities also included. Substituting the matrices into Eq. (10.31) yields

and the three scalar equations

h ( – hcor’ – q’rh ~ Л) = mBl

r(do> d q’

/2U+^J=W52

h ( J + htop’ + p’q’ih ~ h) = mB,

Although the pitch rate of the stability frame is not zero—magnitude of the un­dulating flight-path rate, it is much smaller than the spin rate, and therefore the terms containing q’ can be neglected. With this simplification and solving for the

derivatives yields the simple attitude equations of Magnus rotors:

dp’ _ (hcor1 + mBl) dt I

The first and last equations govern the gyroscopic coupling, which is a function of the spin MOI ratio h/h, the spin rate со, and the angular coupling rates r’ and p’. A positive yaw rate r’ induces a positive change in roll rate p’. whereas the roll rate has the opposite effect on the yaw rate. This purely gyroscopic behavior is modified by the aerodynamic moments mBl and m B}. The second equation simply describes the spin history. In a typical flight the Magnus rotor will reach a steady-state spin rate, indicating that the torque of the driving vanes equals the skin friction on the cylinder and the endplates.

To build a complete six-DoF simulation of Magnus rotors, you would have to program Eqs. (10.25), (10.26), and(10.32), supplemented by quaternion equations for the calculation of the transformation matrix [r]SL. Figure 10.2, with the proper substitutions for Newton’s and Euler’s equations, can serve as a blueprint.

You may have noticed already the great similarity between the equations of a spinning missile and a Magnus rotor. Both use a nonspinning frame that slips over the body, which we called either the nonspinning body frame B’ or the sta­bility frame S. Associated with the frames are the nonspinning body or stability coordinate systems. They serve as the axes for recording the aerodynamic forces and moments in wind tunnels. The gyroscopic coupling is evident in both cases and gives rise to nutation about the velocity vector of the missile and normal to the velocity vector for the Magnus rotor. In either case the development of a full six-DoF simulation would follow the same pattern.

The four models of missile, aircraft, spinning missile, and Magnus rotor cover every conceivable aerospace vehicle that you may have to simulate. You can pro­gram them in matrix form or, in some simple instances, as scalar equations. If guidance and control systems are part of the vehicle, you will have to supplement your simulation with the appropriate code.

One simplification is common to all: the equations of motion were derived with the Earth as inertial reference frame and its curvature unfurled into a flat surface. This assumption is acceptable as long as the vehicle travels at moderate speeds (less than Mach 5) and at low altitudes (less than 30 km). At higher speeds and altitudes, we have to abandon the flat-Earth assumption and consider Earth as a spheroid moving with respect to the sun’s ecliptic.

As you stroll through an armament museum and visit the tactical missile section, you will notice missiles with canard or tail control. Tail-controlled missiles, like the Advanced Medium-Range Air-to – Air Missile, have full three axes stabilization and tight control over the three body rates, roll, pitch, and yaw. Canard-controlled missiles on the other hand exhibit a nasty coupling of pitch and yaw control into the rolling moment. Because this rolling moment is difficult to suppress, the airframe is left to roll freely, or, like in the case of the Sidewinder air-to-air missile, the roll rate is dampened by so-called rolerons. A third type of missile is intentionally rolled for stability. A good example
is the Rolling Airframe Missile, which uses angular momentum for gyroscopic stability.

The first type of missile can be modeled by the equations summarized in Fig. 10.1. Spin-stabilized missiles, on the other hand, require a special set of six – DoF equations that decouple the fast roll rate from the pitch and yaw dynamics. A so-called nonspinning reference frame B’ is introduced that slips like a sheath over the spinning body. Its justification comes from the aerodynamic modeling of spinning symmetrical bodies that gives rise to such terms as the Magnus lift and moment coefficients.

The translational equations are transferable from the flat-Earth equations (10.1) and (10.4) by substituting frame B’ for B.

mDB v + mSlB’EvEB = fap + mg with its matrix version expressed in nonspinning body axes ]B

+ m[nB’E]B’[vE f = [fa, p]B’ + m[T]B’L[g]L (10.19)

The ]B coordinates, because they do not follow the rolling of the missile, are obtained from the local-level coordinates by only two rotations, namely, yaw and pitch, thus, the second axis of]5 always lies in the horizontal plane.

The position of the missile follows from the integration of the second set of differential equations (10.4), expressed in ]5 coordinates.

L

= [tfL[vE}B’ (10.20)

These translational equations of spinning missiles have the same appearance as those of nonspinning missiles. The only difference consists in applying the nonspinning reference frame B’ and its associated coordinate system ]5 . Major differences exist, however, in the attitude equations.

The derivation of the attitude equations starts with Eq. (10.7), applied to the nonspinning body frame B’

DB’lBBE + ПВ’Е1ВВЕ = mB (10.21)

We divide the angular velocity шВЕ of the missile body В wrt the reference frame E into two parts:

The first part is the spin rate, and the second part is the angular velocity of the nonspinning shell wrt Earth. Similarly, the angular momentum is split into two parts:

jBE _ jB. .BE _ jB, ,BB’ , jB. ,B’E lB —1ВШ =lBu> +ІВШ

As called for by Eq. (10.21), we apply the rotational derivative wrt to B’

DB'(lBBuBB’) + Db'(Ibbu:b’e) IBBDB’uBB’ +IBBDBшв’Е

where we made use of DB IB = 0, based on the assumption that the missile has rotational symmetry. Thus, the attitude equation of motion is

lBDB’uBB’ + IbDb’u:b’e + Пв’ЕІвшвв’ + ClB’E Івшв’Е = тв

The first term models the change in angular momentum caused by the variable spin rate of the missile. For a missile with constant spin rate, this term is zero. The second term represents the change of angular momentum of the nonspinning body shell as it is subjected to yaw and pitch rate changes. The gyroscopic coupling moment is expressed by the third term, which usually dominates the weaker coupling of the fourth term.

The attitude equation is expressed in nonspinning body coordinates:

Carefully distinguish between the letters В and B’. The point В is always the c. m. of the missile, and the coordinate system is always the nonrotating ]s ; but the frame may be the body fixed frame B, as in /f, or the nonrolling frame B’. as in uB’E.

Rolling missile airframes usually possess rotational symmetry relative to the Is axis. The MOI tensor therefore has a particular simple diagonal form, which is the same in the ]fl and ]fl coordinate systems. With the two cross-principal moments of inertia being equal /3 = /2, we have

I also included the components of the spin angular velocity [coBB]B and the angular velocity of the nonrolling frame wrt Earth [coDEB. Substituting these matrices into Eq. (10.22) yields

Matrix multiplications deliver the three scalar differential equations

(dco dp’

71U+dг)=тв>

h ( + hcor’ – p’r’ih – h) = тв2

h ( – hm’ + p’q’ih – h) = тВз

By definition, the rolling motion of the nonrolling airframe is zero; therefore dp’/dt = p’ = 0, and the equations simplify to

= h l{-h<or’ + mBl) dr’

— = + тВз)

These are the Euler equations for spinning missiles. Compare them with the missile equations (10.16). They differ in the gyroscopic terms. Because of an overriding spin angular velocity со, the gyroscopic moments lcor’ and lcoq’ dominate the cross coupling. The yaw rate r’ couples into the pitch degrees of freedom and, inversely, the pitch rate affects the yaw rate equation. This coupling gives rise to a nutational mode, which effects the missile dynamics profoundly. The programming of these equations follows the same layout as provided in Fig. 10.1, except that the body rates p, q, and r are replaced by the primed variables q’ and r’, while p’ and ф are zero.

An aircraft, in contrast to a missile, does not have all of its principle axes aligned with the body axes. The 1B axis may be parallel to what is called the water line (an expression borrowed from naval architecture), aligned with the zero-lift attitude, or just poke through the nose of the radome. Whatever the case may be, it is not, in general, the principal MOI axes, nor is the 3B axis. However, thanks to the planar symmetry wrt the Ій and 3B axes, the 2B axis is indeed a principal axis. Such a configuration gives rise to a MOI tensor in body coordinates of the form

In 0 /13

0 /2 0

/31 0 /33

with product of inertia components at the 1,3 and 3,1 locations and /в = /31. On the diagonal we have the axial moments of inertia /ц and /33 and the principal moment of inertia /2.

We can program Euler’s equation (10.10) directly, tolerating the inversion of the moment of inertia matrix, or, because there are still four zero elements in [IB]B, it is possible to solve for the body rate components explicitly and write down the scalar equations.

Before proceeding, I want to expand the angular momentum vector with addi­tional rotary devices like propellers or turbojets, which are an important feature of airplanes. Let the rotary frame be BR and its c. m. BR. The angular momentum is designated Ibbrre and the angular momentum of the remaining rigid aircraft cell is as previously lBBE = ІвшВЕ. Combining both according to Example 6.6, the total

MOI referred to the common c. m. C becomes

‘YaBE ____ jBE і jBr,

c — lB lB

and Euler’s law

DElf + DElBBf= mc

As before, we transform the rotational derivative to the reference frame associated with the aircraft cell В:

DBlf + ПВЕ1ВВЕ + DBlBBf + nBElB/RE = mc

Now, let us simplify the solution by neglecting changes in rotary motion DBlBBE = 0, i. e., assuming that the propellers or turbojets spin at constant speed. With the additional stipulation that the aerodynamic and propulsive moments can be referred to the cell’s c. m. В (rather than C), we obtain the tensor equation for the attitude motions, without incurring any significant errors

IBDBu, BE + Пве(івшве+і%Е) = mB

Expressed in aircraft cell coordinates ]B

and solved for the body-rate derivative

Compare this equation with Eq. (10.10) and note the coupling of the rotary angular momentum [1Brr e]B with the aircraft body rates [Q. BE]B. This vector product models the gyroscopic effect that the rotary engine imparts on the aircraft.

Before we express the component equations, let us make the reasonable assump­tion that the spin axes of all rotary parts are parallel to the aircraft longitudinal axes; therefore, [lBRRE]B = Ur 0 0]. Now with In = hi, we express Eq. (10.17) in coordinates

derivatives is somewhat convoluted. I will spare you the details and just provide the result:

As you can see, the scalar equations of motions of an aircraft are much more complicated than those for missiles, Eq. (10.16). You have a choice now, either to program these scalar equations or the matrix equations (10.17). With today’s computer power at you fingertips, I recommend you use the simpler, less error prone matrix formulation.

Comparing Eqs. (10.16) and (10.18) provides us with some insight into the differences of missiles and aircraft dynamics. Let us first contrast the roll equations. The product of inertia Icouples the pitch and yaw rates to the roll degrees of freedom of the aircraft. The greater its value, the stronger the coupling. If 713 = 0, the missile roll equation is recovered. The pitch and yaw equations exhibit similar trends. Other phenomena caused by A 3 are the coupling of the aerodynamic yawing moment m Bj into the roll axis and the aerodynamic rolling moment m Bl into the yaw axis. Finally, the rotary angular momentum Ir couples through the pitch rate into both the roll and yaw derivative equations, while the yaw rate connects it with the pitch derivative.

Could we have drawn the same conclusions from the matrix equation (10.17)? The answer is yes, in a general sense; it shows the coupling of the body rates IcoBE]B through the MOI tensor [IB]B into the body-rate derivatives, as well as the coupling of the aerodynamic moments [mB]B through the inverse ([Ig ]B)_1. Furthermore, the vector product of [ QDE |B couples the rotary angular momentum [lBlk]B. However, by comparing Eqs. (10.17) with (10.10) we would not have discovered the much simpler missile equations, unless we had scrutinized the simplifying effect of the diagonal moment of inertia matrix.

To sum up, I have assembled the important aircraft equations in Fig. 10.2, just as I did with the missile equations in Fig. 10.1. Euler’s equations are copied from Eq. (10.18), and Newton’s equations are adopted unchanged from the missile model. The incidence angles are now angle of attack a and sideslip angle /3, and the aerodynamic forces are modeled in body axes directly.

The missile and aircraft six-DoF equations are the bread and butter of simulation engineers. They were derived with the assumption that the Earth is the inertial frame and that the curved Earth can be unfurled into a flat plane. How these equations are integrated with other subsystems, such as autopilots, actuators and navigation systems will be described later. Before we leave this section, however, I would like

to derive the especially interesting equations of motions for spinning missiles and Magnus rotors.

Missiles have simple MOI tensors. Because of symmetry, the body axes coincide with the principal axes, and the tensor exhibits only diagonal elements. A particular important case is the missile with tetrago­nal symmetry, which manifests two equal principal MOI. (We exclude here the cruise missiles with their aircraft-like symmetries and attribute them to the aircraft equations in the next section.)

Equation (10.10) with a diagonal MOI tensor [IB]B is in component form

(10.15)

and it is easily expressed in three scalar equations:

~ = A-1 [(A ~ h)qr + mBl]

= 72_1[(7з – h)pr +mBl] (10.16)

^ = A_‘[(A – h)pq+mB3]

These first-order differential equations are nonlinear and coupled only through the aerodynamic moments mд, mj2, and mB, to the translational equations. Their integration yields the body rates p, q, and r. Let me point out again that the term body rates in flight mechanics refers to the angular velocity of the body wrt the inertial frame expressed in body axes, in short [шВЕ]в — [p q r].

Many missiles exhibit tetragonal symmetry, i. e., /3 = /2, and therefore possess rotational equations that are even simpler:

We are now in a position to summarize the equations that form the core of a missile’s six-DoF simulation. They are displayed in Fig. 10.1. The translational degrees of freedom, represented by the velocity components u, v, and w, are solved by Newton’s equation (10.3); and the rotational DoF, expressed in body rates p, q, and r, are governed by Euler’s equation (10.16). The kinematic equations block calculates body attitudes in the form of Euler angles ф, 6, and ф and the direction cosine matrix [T]BL from body rates using the quaternion methodology [Eqs. (10.11), (10.13), and (10.14)]. For aerodynamic table look-ups the incidence angles a1 and ф’, peculiar to missiles, are displayed [see Eqs. (3.22) and (3.23)]. Finally, the aerodynamic coefficients and propulsive forces are combined in the forces and moments block. Notice the special missile features of expressing the force coefficients as Сд, Су, and Сц in body axes. Refer to Sec. 10.2.1.3, in which I treat the modeling of missile aerodynamics in greater depth.

As I mentioned, these blocks form just the basic loop of a six-DoF simulation. In a full-up missile model they are joined by an autopilot and actuator to form the inner loop. An outer loop, supplying the navigation and guidance functions, completes the missile simulation. Modeling details of these additional components will be presented in Secs. 10.2.2, 10.2.4, and 10.2.5.

Unless you work exclusively with space applications, you will encounter the flat-Earth equations of motion more often in six-DoF simulations. We first derive the general form of the equations of motion and then specialize them for missiles
and aircraft based on the characteristic of the MOI tensor. Spinning missiles and Magnus rotors will each require a fresh approach.

Let us first derive the translational equations of a flight vehicle subjected to the aerodynamic and thrust forces fa and the gravitational acceleration g. Newton’s second law wrt the inertial frame I states that the time rate of change of linear momentum equals the externally applied forces:

= f a p + mg

where m is the vehicle mass and vB is the velocity of the missile c. m. В wrt the inertial reference frame I. The flat-Earth assumption allows us to declare the Earth frame E as the inertial frame. Therefore, Newton’s law becomes

mDEvEB = fa’P+mg

To model the aerodynamic forces, we need access to the incidence angles. They are calculated from the velocity vector as perceived from the vehicle as reference frame. Its time rate of change DBvf is therefore referred to the vehicle B. This desired shift of reference frame is accomplished through the Euler transformation

mDBvB + mflBEvB = fa p + mg (10.1)

These are the translational equations of a vehicle over a flat Earth. The second term on the left-hand side is identified as the tangential acceleration term. The equation is valid in any coordinate system.

Now comes the important event of picking the coordinate system. Actually, I have already preempted the decision by transforming the rotational derivative to the body frame. To generate the ordinary time derivative, the rotational derivative must be expressed in the coordinate system ]B associated with the frame B. Expressing all terms in ]B yields

+ m[QBE]B[vE]B = [fa, p]B+m[g]B

This selection suits us well because it models [vB]B, the linear velocity of vehicle wrt Earth, as needed for incidence angle calculations; [Q. BE]B, the angular veloc­ity of vehicle wrt Earth; and [fa, P]B, the aerodynamic and propulsive forces, all expressed in the preferred coordinate system ]B. Only the last term is better for­mulated in the local-level coordinate system ]L so that g takes on the simple form [g]L = [0 0 g], Alas, the direction cosine matrix [T]BL appears, which we will derive later from Euler’s equations. The translational equations in matrix form, suitable for programming, are then

+ m[0.BE]B[ vE]B = [fa, p]B + m[T]BL[g]L (10.2)

Much labor is involved in modeling the forces [fa, p]B. In the next section I will give you a taste for the complexity of six-DoF aerodynamics, both for missiles and aircraft. Here we proceed blissfully and code the equations directly with program languages like FORTRAN or C. You can also choose the simulation environment

CADAC, which provides the utility subroutines that enable vector state variable integration and matrix manipulations. Written out in coordinate form, we get

These equations are simple enough to be expressed in scalar form:

The integration of these differential equations yields the velocity vector that must be integrated once more to obtain the location of the missile c. m. В wrt an Earth reference point E

[Ле] = [uf] (10.4)

The integration is best carried out in the local-level coordinate system. Therefore, given [uf ]B, we program the second set of differential equations as

and in coordinate form

(10.6)

Here we have the six first-order differential equations that govern the transla­tional motions of a vehicle with the Earth as the inertial reference frame. Equa­tions (10.3) are nonlinear and coupled by the body rates p, q, and r and the direction cosine matrix [T]BL with the attitude equations. The nonlinearity enters through the incidence angles in the aerodynamic force and moment calculations [for instance, seeEqs. (10.61) and (10.62)]. Equations (10.6) are linear differential equations, again coupled by [T]BL with the attitude equations.

The rotational degrees of freedom are governed by Euler’s law that states that the time rate of change of angular momentum equals the externally applied moments.

To conform to the translational equations, we pick E as the inertial frame:

DElBE = mB (10.7)

where lBBE = ІвшВЕ is the angular momentum of body В wrt frame E referred to the c. m. B, IBB is the MOI of missile body В referred to the c. m., and mB the aerodynamic and thrust moments referred to the c. m.

Just as in the case of the translational equations, we take the perspective of the vehicle body frame В and use therefore the Euler’s transformation to transfer the rotational derivative to the body frame

DBlf + ПВЕ1ВВЕ = mB

Let us focus on the rotational derivative, expanding the angular momentum vector, and applying the chain rule

DBlBBE = Db(Ibbube) = DbIbbube + 1BBDB шВЕ = IbDbujbe

For a rigid body DBIB is zero, a simplification that motivates the transformation to the body frame. The equation reduces to

IbDbube + ПВЕІвиВЕ = mB (10.8)

These are the rotational equations of a vehicle with respect to a flat Earth. The second term on the left-hand side is the apparent gyroscopic effect, arising from the choice of the body as reference frame.

As a coordinate system, we pick again the body coordinates because they express the MOI tensor in constant form, yield the ordinary time derivative of the body rates, and express the aerodynamic moments in their preferred form:

+ [Qbe]b[Ib]B[wbe]b = [mB]B (10.9)

For an arbitrary [IB]B we premultiply the equation by its inverse and solve for the time derivative

B = {[Ibb}B)–[Q. be]b[Ib]BWbe}b + [mB]B) (10.10)

These three first-order nonlinear differential equations couple with the translational equations (10.2) only through the aerodynamic moments [mB]B. The form of the MOI tensor [IB]B establishes a major demarcation. Missiles exhibit pure diagonal forms, resulting in simple equations, whereas aircraft have one off-diagonal ele­ment that leads to more complex solutions. We reserve therefore in the following sections different treatments for missiles and aircraft. Here we complete that part of the treatment of the attitude kinematics, which is common to both.

Given the body rates [coBE]B = [p q /’] from Eq. (10.10), a second set of differ­ential equations provides the body attitudes. We discussed the three possibilities in Chapter 4: Euler angle, quaternion, or direction cosine equations. Here we pursue the quaternion approach. Go back to Sec. 4.3.3 and review quaternions, or, like many other impatient engineers, just program the following equations.

The four-quaternion elements qa, qi, qi, and q- are related to the three body-rate components by four linear differential equations:

(10.11)

To solve the equations, they must be initialized. Unfortunately, quaternions are hard to visualize; therefore, it is more likely that the initial Euler angles are given. The conversions are

For example, an aircraft’s Euler angles are zero if it starts flying north, straight, and level. Its quaternions are also zero with the exception of the scalar component Яо = 1.

If the differential equations (10.11) could be solved exactly, we would not have to worry about maintaining the orthonormality of the quaternion. However, the use of numerical integration schemes introduces errors caused by finite word length and discretization of time. We can maintain the unit norm of the quaternion by a proven trick.1 Let the orthonormality error X be

* = 1 – (Яо + V? + ЯІ + Яз)

and add to the right side of Eq. (10.11) the factor kX{q}, where {q} is the four – element quaternion vector. The constant к is chosen such that к At < 1, with At the integration interval. Its value is not crucial. I have tried several к and found the effect quite acceptable. For the CADAC SRAAM6 simulation, after some testing I picked к = 0.5 at At = 0.001.

Given the quaternions, the Euler angles are derived from the following relation­ships:

The first equation has singularities at [r = ±90 deg, and the last equation at ф = ±90 deg. These singularities are not serious because they occur during off-line calculations and not inside the differential equations. They can easily be bypassed by programming around them.

The direction cosine matrix [T]BL could be calculated from the Euler angles. However, because of the singularities, it is better to use the quaternion relationship directly:

(10.14)

In summary, the three translational degrees of freedom are governed by the six first-order differential equations (10.2) and (10.6). The three rotational degrees of freedom are calculated from the three first-order differential equations (10.10) and the four kinematic differential equations (10.11).

With the general equations of motion of a vehicle over a flat Earth in place, we now turn to special expressions of Eq. (10.10). The form of the MOI tensor will determine whether we call them missile or aircraft equations.

Before we embark on our journey, it would be to your advantage to stop by at Chapter 5 and review Newton’s law and Chapter 6 on Euler’s law. Otherwise, starting from these general principles the derivation here is self-contained.

As always, I shall formulate the equations in an invariant tensor form first, followed by the all-important decision of coordinate systems. The form of the MOI tensor will lead us into two different directions. For missiles with tetragonal sym­metry, the equations are so simple that we write them in scalar form. For aircraft with planar symmetry however, we have to live with an inverted moment of inertia matrix.

The selection of the inertial reference frame is always an important decision that you have to make before you begin your modeling task. For near-Earth aerospace vehicles, including low-Earth orbiting satellites, the Earth-centered, ecliptic-oriented frame serves as inertial reference (J2000). It can be combined with a spherical or elliptical Earth shape. However, frequently the trajectory of a vehicle is so slow compared to orbital speed and so close to the Earth that the cen­trifugal and Coriolis accelerations can be neglected. In these cases we choose the Earth as the inertial frame and unwrap the Earth’s surface into a plane tangential to the launch point. We speak now of the flat-Earth approach.

The choice is yours to make. I will first discuss the easier flat-Earth case, appli­cable to tactical missiles and aircraft and add two special cases: spinning missile and Magnus rotor. Then I will lure you over to the much more difficult modeling task of the elliptical Earth. We will make excursions into geodesy and study the shape and gravitational field of the Earth before we derive the equations of mo­tions. The inertial frame is the Earth-centered, ecliptic-oriented frame. To provide you with the full spectrum of formulations, I will not only use the inertial frame as a reference frame but also the Earth frame. You will encounter both options in six-DoF simulations.

We have assembled all of the gear we need to embark on the ultimate trip of full six degrees of freedom. In Chapter 5 you learned to model the three translational degrees of freedom, and Chapter 6 taught you the secrets of three rotational degrees of freedom. Newton and Euler should have become your favorite companions.

You also have become familiar with simple simulations. Chapter 8 instructed you how to build a three-DoF simulation and under what circumstances it can satisfy your needs. Then you moved up one notch and added two more degrees of freedom in Chapter 9 to arrive at that peculiar form of pseudo-five-DoF simulations. Both have important roles to play during the development cycle of an aerospace vehicle. As an engineer, you have to adjust your tools to the task at hand and the resources available. Start simple! It is better to have a reliable three-DoF simulation ready in time than a nonvalidated and tardy six-DoF simulation.

Eventually, the time will come when you have to strike out and go for maximum fidelity. But count your cost, and let your sponsor know the extensive resources he has to set aside. You will need lots of data: full aerodynamic and thrust tables— trimmed data will no longer suffice; a complete flight control design; mass pa­rameters that include moments of inertia; and other subsystems, like sensors and guidance logic. You also will need lots of time. I have yet to finish a simulation, completely verified and validated, within the projected period. From experience. I know that once a simulation runs reasonably well it takes another third of the development time to produce a reliable product.

The challenge is great, but very rewarding. I will help you with the details. We develop together two types of equations of motion. For in-atmosphere, near-Earth flight we assume the flat Earth to be the inertial frame. Tactical missiles and aircraft fall into this category, and I will add spinning missiles and Magnus rotors to the collection. Quaternions will serve us to compute the attitude angles.

Other types of equations are required for hypersonic and orbital vehicles. They are based on an inertial frame aligned with the solar ecliptic and convecting with the center of the Earth—the inertial frame J2000 you encountered in Chapter 3. We shall study the elliptical geometry of Earth, the gravitational field, and introduce the geodetic coordinate system. Then I shall formulate the equations of motion relative to the inertial and geographic frames. The attitude angles will be computed from the direction cosine matrix equations.

Next, we have to flesh out the right-hand side of the equations of motion. I present aerodynamic models for airplanes, hypersonic vehicles, and missiles. The autopilot for six-DoF simulations requires an entirely new treatment with such options as rate damping loop, roll autopilot, heading angle, and flight-path-angle controller, acceleration autopilot, and altitude hold autopilot. Actuators are now also part of the control loop.

All of the seeker and guidance options of five-DoF simulations can directly be integrated into your six-DoF. I shall just expand your horizons. The modeling of INS errors becomes now important; extensions of proportional navigation and the advanced guidance law will improve the intercept performance of your missile, and a gimbaled IR sensor provides you with a modem missile seeker.

This chapter will not satisfy all of your needs. Its purpose is mainly to furnish enough examples so that you can develop the specialized components for your own simulation. It should also help you decipher six-DoF simulations when you are asked to evaluate or to modify them. For that reason, I refer you also to the simula­tion examples on the CADAC CD, which enrich this chapter; theFALCON6 aircraft simulation, the GHAME6 hypersonic vehicle, and the SRAAM6 short – range air-to-air missile.