1. INTRODUCTION
Solution of the six-degree-of-freedom flight equations for aircraft and missiles continues to represent one of the most important application areas for analog, hybrid, and digital computer systems. Important computer requirements such as accuracy and speed are very much dependent on the choice of axis system for the translation equations of motion. In this connection it is well known that the flight-path axis system (or wind-axis system) makes much lower accuracy and speed demands on the computer than does the body-axis system. /1/, /2/
Flight-path or wind axis (w), stability axis (s) and body axis (b)
Despite this, a number of current computer mechanizations continue to use body axes for solving the translational equations of motion. Because of this, unnecessary demands on accuracy or frequency response are placed on the computer, and many mechanizations which could be all-analog or all-digital have shifted to hybrid implementation. Even if the mechanization is hybrid from the outset, there is considerable advantage to be gained by using an efficient axis system. The purpose of this paper is to point out again the advantages of flight path axes and to summarize the overall equation requirements for solving the six-degree of freedom flight equations.
2. BODY-AXIS TRANSLATIONAL EQUATIONS
For comparison, we present first the body-axis translational equations. The body axes Xb, Yb, and Zb are defined as a right-hand set fixed to the vehicle with the Xb axis along the longitudinal axis and the Zb axis directed downward for normal level flight. The components of the total vehicle velocity vector Vp along the Xb, Yb , and Zb axes, respectively, are Ub, Vb, and Wb (see figure 2.1).
The components of the body-axis angular velocity vector Ω (and hence the vehicle angular velocity vector) along the Xb, Yb, and Zb axes are Pb, Qb, and Rb, i.e., roll rate, pitch rate, and yaw rate, respectively. Here we assume Vp and Ω represent vehicle translational and rotational velocity vectors as viewed from an inertial (nonaccelerating) frame of reference. If we denote the external forces along a set of coordinate axes by X, Y, and Z, respectively, then Euler’s translational equations of motion, which are obtained by summing forces along the coordinate axes, are the following:
where m is the mass of the vehicle.
The inefficiency of these equations in body axes is immediately apparent when one considers the approximate size of the various terms. Let the vehicle be, say, a Mach 2 aircraft with Vmax = 2000 feet per second. A reasonable upper limit on pitch-rate Q, might be 2 radians per second. Thus, the term UbQb in equation (2.3) might be as large as 4000 ft/s2 or 125 g’s. On the other hand Zb/m, the normal acceleration due to the external force (primarily gravity and aerodynamic lift) may have an upper limit of several g’s. Hence artificial accelerations which are perhaps 20 to 50 times greater than the actual accelerations are introduced because of the high rotation rates which the body-axes experience. This means much less favorable computer scaling and hence much poorer solution accuracy for a given computer precision. Furthermore, the high-speed dynamics of the rotational equations are coupled into the translational equations, thus placing severe dynamic response requirements on the computer.
The use of flight-path axes greatly alleviates these problems. As we shall see in the next section, the flight-path axes allow a more efficient calculation of the aerodynamic angle of attack α (alpha) and the aerodynamic angle of side slip β (beta) than body axes allow. Using body axes and assuming that the ambient air mass is not moving relative to the inertial frame used to define Vp, then the following formulas can be used to obtain α, β, and velocity magnitude V, from the body-axis velocity components Ub, Vb, and Wb:
3. FLIGHT-PATH AXIS TRANSLATIONAL EQUATIONS
Next consider the flight-path axes Xw, Yw, Zw, shown in figure 3.1.
These differ from the body axes Xb, Yb, Zb by the angle of attack α and the angle of side slip β, as sh own in the figure. To rotate from body axes to flight-path axes one first pitches the body axes about Yb through -α. This defines an intermediate axis system Xs, Ys, Zs, called stability axes. One then yaws about Zs, through β, which defines the flight-path axes Xw, Yw, and Zw. Note that Vp is the Xw component of vehicle velocity; the Yw and Zw components are zero by definition. Let us define the Xw, Yw, and Zw components of flight-path angular velocity relative to inertial space by Pw, Qw, and Rw, and the components of external force along Xw, Yw, and Zw by Xw, Yw, and Zw, respectively. Then, since Uw = Vp and Vw = Ww = 0, the translational equations (2.1), (2.2), and (2.3) referred to the flight-path axes become
Solution of these three equations results in total velocity Vp, flight-path axis yaw rate Rw, and flight-path axis pitch rate Qw.
Next consider the formulas for α and β. Reference to figure 3.1 shows that α_dot is directed along Ys, with a component α_dot cos β along Yw. Thus α_dot cos β is equal to the difference between body-axis and flight-path axis angular rates along Yw. Therefore, from equation (3.3) we can write
where Pb_s is the body-axis (not stability axis) angular rate along Xs and is given by (3.5)
Similarly, reference to figure 3.1 shows that β_dot is directed along Zw and is equal to the difference between flight-path axis and body axis angular rates along Zw. Thus from equation (3.2) we get (3.6) where Rb_s is the body-axis angular rate along Zs and is given by (3.7). Note that Qb_s = Qb since the Yb body axis and Ys stability axis are coincident (3.8).
Equations (3.1), (3.4), and (3.6) can be integrated to yield total velocity Vp, angle of attack α, and angle of side slip β. They do not present the scaling difficulty of equations (2.1), (2.2), and (2.3) in the body axes. The body-axis velocity components Ub, Vb, and Wb can be obtained from Vp, α, and β by the formulas (3.9) -3.11). Similarly, the flight-path-axis forces Xw, Yw, and Zw can be derived from body-axis force components Xb, Yb, and Zb by the formulas (3.12)-(3.14) where Xs, Ys, and Zs are the intermediate stability-axis force components.
Frequently the aerodynamic force components are computed along stability axes, in which case the power plant and gravity forces, computed in body axes, would be resolved into stability axes where the aerodynamic forces are added; then the total forces would be resolved into flight-path axes to allow use of equations (3.1), (3.4), and (3.6).
It should be noted that we have assumed throughout this section that the translational and rotational velocity vectors are velocities relative to an inertial reference frame. If the atmosphere through which the vehicle is flying can be considered to be fixed with respect to this inertial frame, then the velocity magnitude Vp is the vehicle velocity relative to the atmosphere, and α and β represent the aerodynamic angle of attack and side slip, respectively. Using the approximation that the earth is flat, and with constant surface winds, it is possible to define the inertial reference frame as a frame attached to the atmosphere. Then all the formulas, as presented, are correct, and α, β and Vp can be used for computation of aerodynamic forces and moments.
Unfortunately, for a rotating spherical earth, axes fixed in the atmosphere are not inertial. If we consider such a frame to be inertial, we will make acceleration errors in equations (3.1), (3.2), and (3.3) of the order of V²/r0 where V is the vehicle velocity relative to an inertial frame with origin at the center of the earth and r0 is the radius of the earth. For illustration, consider a vehicle flying eastward in still air with a velocity Va relative to the atmosphere. Then V ~ Va + r0 ωn cos L, where ωn is the earth spin rate and L is the latitude. If the vehicle is flying westward, V ~ Va - r0 ωn cos L. For other headings V lies between these values.
For example, if Va = 3000 ft/s and L = 0 degrees, V ranges between approximately 1600 and 4400 ft/s. The corresponding acceleration error in equations (3.1), (3.2), and (3.3), given by V²/r0, ranges between approximately 0.1 and 1 ft/s2. For many flight vehicles this is a negligible error. On the other hand, for a supersonic transport cruising eastward at 3000 ft,/s this lowers the required steady-state lift by about 3per cent, which could lower the drag significantly and hence make a noticeable difference in maximum range.
There appears to the authors to be no simple way to take these accelerations into account and still use a flight-path axis system referenced to the atmosphere. One could add an approximate correction acceleration in the vertical direction given by V²/r0 to the translatory forces in equations (3.1), (3.2), and (3.3). In fact, one could further simplify the computation by adding the term only to equation (3.3), based on the argument that most of the time a supersonic aircraft will be in near-level flight at cruise and that a moderate acceleration error during transient conditions can be accepted.
In this case equation (3.3) becomes
Note that this equation will exhibit acceleration errors of the order of V²/r0 for conditions far from level flight.
4. ROTATIONAL EQUATIONS
The only reasonable axes to use for the rotational equations of motion are the body axes. If we let the external moment components along Xb, Yb, and Zb, be Lb, Mb, and Nb, respectively, then summation of moments about the three body axes of a body symmetrical about the XbZb plane leads to the equations:
Here Ixx, Iyy, and Izz are the moments of inertia about Xb, Yb, and Zb respectively, and Ixz is the product of inertia of the symmetrical body. Note in equations (4.1), (4.2), and (4.3) that the second term in each equation represents a nonlinear inertial coupling term. For flight vehicles such as large transport aircraft which do not generate relatively high angular rates these terms often can be neglected. For many flight vehicles the roll rate Pb has a maximum value which is considerably higher than pitch-rate Qb or yaw-rate Rb. Hence the QbRb term in equation (4.1 ) often can be neglected in comparison with the RbPb and PbQb terms in equations (4.2) and (4.3), respectively.
The third term in each of equations (4.1), (4.2), and (4.3) represents the effect of the product of inertia Ixz. If the x, y, z body axes have been chosen to be almost coincident with the principal axes, this term may be negligible in all three equations, since Izz will be very small compared with the principal moments of inertia. For relatively low angular rates the nonlinear terms (PbQb, Pb² - Rb², and QbRb) usually can be neglected and in any event Rb² frequently can be neglected as small compared with Pb² in equation (4.2).
5. COMPUTATION OF EULER ANGLES FOR A FLAT EARTH
Solution of equations (4.1), (4.2), and (4.3) results in computation of the body-axis angular velocity components Pb Qb, and Rb. These must be integrated a second time to obtain the orientation of the vehicle body axes with respect to the desired references axes, typically Euler axes which point north, east, and toward the center of the earth (Xe, Ye, and Ze in figure 5.1 ).
This orientation usually is expressed in terms of the conventional aircraft Euler angles, i.e., heading-angle Ψ (psi), pitch angle Θ (theta), and bank angle Φ (phi). These angles usually are computed from Pb, Qb and Rb by the following well-known equations:
Note that, since Pb, Qb, Rb are the body-axis components of the vehicle angular velocity relative to an inertial reference system, there is a small error introduced by the angular velocity of axes which point north, east, and down relative to a spherical earth. This small error can be corrected, if necessary, using equations (7.1)-(7.3) in the next section.
The well known singularity of the Euler angle system at Θ = ±π/2 can be avoided, if necessary, by computing direction cosines or quaternions /4/ instead of Euler angles, or by introducing a fourth angle. /3/
The use of quaternions rather than direction cosines should be considered if a system free of singularities is needed, since there are only four quaternions with a single redundancy as compared with nine direction cosines with six redundancies.
5.1 Additional notes about Euler angles /6/
The most popular application is to describe aircraft attitudes, normally using a Tait–Bryan convention so that zero degrees elevation represents the horizontal attitude. Tait–Bryan angles represent the orientation of the aircraft respect a reference axis system (world frame) with three angles which in the context of an aircraft are normally called Heading Ψ (psi), Elevation Θ (theta) and Bank Φ (phi). When dealing with vehicles, different axes conventions are possible.
Tait–Bryan Euler angles in aircraft application (all zero degrees yield normal flight attitude north)
When studying rigid bodies in general, one calls the xyz system space coordinates, and the XYZ system body coordinates. The space coordinates are treated as unmoving, while the body coordinates are considered embedded in the moving body. Calculations involving acceleration, angular acceleration, angular velocity, angular momentum, and kinetic energy are often easiest in body coordinates, because then the moment of inertia tensor does not change in time. If one also diagonalizes the rigid body's moment of inertia tensor (with nine components, six of which are independent), then one has a set of coordinates (called the principal axes) in which the moment of inertia tensor has only three components.
The angular velocity of a rigid body takes a simple form using Euler angles in the moving frame. Also the Euler's rigid body equations are simpler because the inertia tensor is constant in that frame.
"It is important to appreciate the range of the Euler angles:
−π ≤ ψ ≤ +π
−π/2 ≤ θ ≤ +π/2
−π ≤ φ ≤ +πIt is not obvious that pitch should be constrained to the range ±90◦. Consider a loop; as the pilot pulls up into the loop, the pitch increases to 80◦, then85◦, 86◦, 87◦, 88◦ and 89◦. At this point, the aircraft is still pointing forwards, but almost vertically upwards. In the next two degrees of motion, the aircraft passes through 90◦ and subsequently is pointing backwards, but at 89◦ to the horizon. Notice also that during this maneuverer, as the pitch angle passes through 90◦, the aircraft heading changes from forwards to backwards (an 180◦ instantaneous change in yaw) and also, the aircraft changes from upright to inverted (an 180◦ instantaneous change in roll)." /8/
German words: "Die Lage des körperfesten Achsensystems wird mit Hilfe der drei Eulerwinkel Ψ (Azimut, Gierwinkel), Θ (Längsneigung, Nickwinkel) und Φ (Hängewinkel, Rollwinkel) beschrieben." /7/
6. COMPUTATION OF VEHICLE POSITION FOR A FLAT EARTH
Once the orientation of the flight vehicle with respect to the Euler axes has been established e.g., by means of the Euler angles, then it is possible to compute the velocity north Ue, the velocity east Ve, and the velocity downward We, or its negative, the rate of climb h_dot. Direct integration then yields the vehicle position.
Determination of Ue, Ve, We is complicated by the fact that the vehicle velocity vector Vp lies along the x wind axis and therefore must be resolved from wind to earth axes. Unfortunately, the complete orientation of the wind axes is known only relative to body axes; hence it is necessary to perform the resolution of Vp into earth axes by first resolving it into body axes,. then from body axes to earth axes. The resolution of Vr from wind axes to body axes is accomplished by the transformation given in equations (3.9), (3.10), and (3.11 ).
The resolution of vehicle velocity from body axes to earth axes can be accomplished by using direction cosines. Thus let l1, l2, and l3 be the projections of a unit vector along the x body axis onto the Xe, Ye, and Ze earth axes, respectively. Similarly, let m1, m2, and m3 be the projections of a unit vector along the y body axis onto the Xe, Ye, and Ze earth axes, respectively. In the same way let n1, n2, and n3 be the projections of a unit vector the z body axis along the Xe, Ye, and Ze earth axes, respectively. Then by definition
It is easy to show that the direction cosines are related to Euler angles by the following formulas:
(Equivalent formulas for direction cosines in terms of quaternions are given in reference 4.)
An alternative mechanization of equations (6.1 ) through (6.6) avoids computation of the direction cosines by instead performing successive resolutions of the velocity components Ub, Vb, and Wb through the angles -Φ, -Θ, and -Ψ. Consider the intermediate axis system x’, y’, z’ in figure 5.1. Clearly the velocity components U’, V’, and W’ along x’, y’, and z’ are given by
Defining U", V", and W" as the velocity components along the intermediate axes x", y", and z", we have (6.8). Finally, from figure 4.1 we see (6.9).
Successive application of equations (6.7), (6.8), and (6.9) for U, V, and W to obtain Ue, Ye, and We requires fewer mathematical operations than using equations (6.1) through (6.6). It therefore has computational advantages using either an analog or digital mechanization. In computing ground coordinates from Ue and Ve in equations (6.1 ) and (6.2), or equation (6.9), it is important to note that Ue and Ve represent airspeed components north and east, respectively. To convert them to groundspeed components, the north component of wind Wz must be subtracted from Ue, and the east component of wind Wy must be subtracted from Ve. Thus if Sz and Sy represent distance traveled north and east, respectively, then
Equations (6.10) and (6.11) are valid only for steady winds, since the implicit assumption was made that axes stationary relative to the atmosphere are inertial. Equations (3.1), (3.2), and (3.3) are referred to inertial space and correction terms must be added if the reference axes are not inertial.
7. COMPUTATION OF VEHICLE EULER ANGLES AND POSITION FOR A ROTATING SPHERICAL EARTH
In the previous section we presented the formulas for computing vehicle position over a flat earth with steady winds. We can use the same position formulas to obtain velocity north Ue and velocity east Ve over a rotating spherical earth with radius r0 and angular rate ωN. However, Ue and Ve will represent airspeed components and must be corrected to yield groundspeed components Sx_dot and Sy_dot respectively. Furthermore, when a spherical earth is considered, it may be necessary to correct the vehicle angular rates used to compute Euler angles in order to take into account the rotating reference frame, as pointed out earlier in section 6. It can be shown that the body-axis components of the vehicle angular velocity relative to the Euler reference frame, Pbe, Qbe, Rbe are given by the following formulas: /5/
Here r is the radial distance of the vehicle from the center of the earth and is given by (7.4) where r0 is the radius of the earth and h is the vehicle altitude. In many cases we can substitute r0 for r in equations (7.1), (7.2), and (7.3) and still obtain sufficient accuracy.
The values of Pbe, Qbe, and Rbe in these equations are then used to compute the Euler angle rates. Thus by analogy with equations (5.1), (5.2), and (5.4)
In many cases the computations involved in equations (7.1), (7.2), and (7.3) can be neglected i.e., we can assume that Pbe ~ Pb, Qbe ~ Qb, Rbe ~ Rb. This is particularly true if the overall six-degree-of-freedom computation involves a control system (automatic or human) which attempts to maintain Ψ, Θ, and Φ at specified values. In any case the correction rates are the order of Vp/r0. For example, if Vp = 2000 ft/s, the correction rate is equal to approximately 0.005 degree per second. On the other hand, if the flight-vehicle problem includes a stable platform, the rate corrections given by equations (7.1), (7.2), and (7.3) may be important.
It should be noted that Sz_dot and Sy_dot in equations (7.1 ), (7.2), and (7.3) represent vehicle velocity components north and east, respectively, over the surface of a nonrotating earth with steady winds. On the other hand equations (6.1 ) and (6.2), or alternatively, equation (6.9), gives us Ue and Ve, i.e., vehicle velocity components north and east, respectively, relative to the inertial reference frame for the translational equations of motion. We made the approximation in section 3 that this reference frame is fixed relative to the ambient atmosphere. We can account for the linear velocity of the atmospheric reference frame (but not the angular velocity) by noting that relative to the surface of a nonrotating earth it is moving northward with the northerly component of wind and eastward with the sum of the rate due to earth spin and that due to eastward component of wind. Thus we can write the following equations:
It can also be shown /5/ that the time rate of change of latitude and longitude are given by the formulas (7.10) and (7.11).
It should be noted that motion of a flight-vehicle over a rotating earth can be treated exactly, /5/ but that the exact translational equations referred to axes fixed relative to the ambient atmosphere are very complicated. Thus many of the computer mechanization advantages for flight-path axes are lost, and the computation of α, β, and Vp is much less elegant. We have attempted in this section to describe how one can utilize the flight-path axis system and still correct approximately for the fact that the vehicle is flying over a rotating earth with surface winds. This approach should be adequate for all but the most exacting requirements, unless the vehicle reaches hypersonic speeds. For subsonic vehicles or supersonic vehicles traveling over relatively short distances the flat-earth equations in sections 3 through 6 should be adequate.
8. COMPUTATION OF AERODYNAMIC FORCES AND MOMENTS
Aerodynamic forces and moments for flight vehicles normally are computed in stability axes. Thus let Xa, Ya, and Za be the aerodynamic forces along the Xs., Ys, and Zs stability axes shown in figure 3.1. Here -Xa corresponds to drag D, and -Za corresponds to lift L.
In addition to aerodynamic forces the flight vehicle experiences propulsion and gravity forces. These are most conveniently defined in body axes. Let us denote the propulsion force components along the x and z body axes by Xp and Zp, respectively. The gravity force components along x, y, and z will by definition be mgl3, mgm3, and mgn3 respectively. The sum of the propulsive forces and gravity forces along body axes can then be resolved to stability axes, where the aerodynamic forces are added to obtain the total force component Xs, Ys and Zs along the stability axes. Using equations (6.4), (6.5), and (6.6) to express the direction cosines l3, m3, n3 in terms of Euler angles, we obtain:
Finally, reference to figure 3.1 shows that the force components Xw, Yw, and Zw along the flight-path axes can be computed using the formulas (8.4) - (8.6).
These force components are used to mechanize the translational equations given in section 3. If the aerodynamic forces are given in body axes rather than stability axes, the modification of equations (8.1 ), (8.2), and (8.3) is apparent.
The external moments acting on the flight vehicle consist of aerodynamic moments La, Ma, and Na, normally given in stability axes, and power plant moments Lp, Mp, and Np, given in body axes. The total moments L, M, and N in body axes become the following:
Again, if the aerodynamic moments are given in body axes, the simplification of equations (8.7), (8.8), and (8.9) is obvious (set α = 0). Figure 8.1 shows a block diagram of the over-all six-degree-of-freedom equations for the case of a flat earth.
Figure 8.1
L. E. FOGARTY
R. M. HOWE,
University of Michigan
(REPRINT of the original with minor fixes and additions.)
REFERENCES
/1/ Howe R M
Coordinate systems for solving the three-dimensional flight equations
WADC Technical Note 55-747 June 1956
/2/ Howe R M
Coordinate systems and methods of coordinate transformation for three-dimensional flight equations
Proceedings of the First Flight Simulation Symposium WSPG
Special Report 9 September 1957
/3/ Greenwood D T
An extended Euler angle coordinate system for use with all-altitude aircraft simulators
WADC Technical Report 60-372 August 1960
/4/ Robinson A C
On the use of quaternions in simulation of rigid-body motion
WADC Technical Report 58-17 December 1958
/5/ Fogarty L E, Howe R M
Analog computer solution of the orbital flight equations
SIMULATION vol 1 no 1 Fall 1963 pp R-41 through R-47
/6/ http://en.wikipedia.org/wiki/Euler_angles
/7/ http://tumb1.biblio.tu-muenchen.de/publ/diss/mw/2004/holzapfel.pdf
/8/ David Allerton,
PRINCIPLES OF FLIGHT SIMULATION, Department of Automatic Control and Systems Engineering, The University of Sheffield
Posts
