Category Modeling and Simulation of Aerospace Vehicle Dynamics

10.2.2.5

10.2.2.6 As mentioned in Sec. 9.2.3.3, there is little difference between an altitude hold autopilot of a pseudo-five-DoF or six-DoF implementation. Refer back to Fig. 9.19 and substitute the acceleration autopilot that we just designed. The altitude and rate measurements are provided by the INS. As before, the root locus analysis will also serve us here and will help us to determine the two gains Gv and Gh – However, I am not building an adaptive gain scheduler as before because altitude corridors are usually fixed, and a constant set of gains will suffice. The hypersonic vehicle GHAME is controlled by this autopilot. You can see the code implementation in the CADAC GHAME6 simulation, Module C2, Subroutine C2ALT.

10.2.2.7 Flight-path-angle controller. Finally, I want to present to you one more autopilot function, the flight-path-angle controller. The aircraft relies on it for climb and descent, especially if the flight-path angle is constrained by air safety rules. You will find this design eminently interesting because it is a perfect example for the pole placement technique of state variables, as embodied by Eqs. (10.77) and (10.78).

To be specific, we replace the general variable names of Fig. 10.16 with those of our example and present it inFig. 10.19. The commanded and achieved flight path

angle is yc and у, respectively, supplemented by two more state variables: pitch rate q and pitch angle в. All three states are provided by the INS. The three feedback gains c, multiplied by the state vector*, produce the scalar feedback variable r. The elevator deflection be is sent to the actuator. We proceed to calculate the feedback gains c and the feedforward gain k.

First, we copy the plant equation x = Fx + gSe from Eq. (7.73)

with the dimensional derivatives La, L^e, Ma. Mq, M$e. Substituting the plant into Eq. (10.76) yields on both sides third-order polynomials in s, the same order we encountered in the acceleration autopilot. We can, therefore, carry over Eqs. (10.82) and (10.83), but, unfortunately, we cannot solve for the gains by elimination. Instead we have to use Eq. (10.78) and invert the matrix P to calculate the feedback gains:

Once we have c. we proceed to calculate the feedforward gain к from Eq. (10.77). All variables are known, and h = [0 0 1] picks out the state variable y. Using к in the forward loop ensures that the steady-state gain is one, and у will track yc after the transients have died down.

10.2.2.7 Summary. We accumulated a whole assortment of autopilot de­signs for our six-DoF simulations. The rate dampers stabilize missiles off the launch rail and augment the directional stability of aircraft. Roll autopilots are used either to suppress the roll excursions for skid-to-turn missiles or to bank an aircraft into a turn. If this turn leads into a new direction, the heading tracker ensures that the new bearing is maintained. Acceleration autopilots are squarely in the domain of missiles. They execute guidance commands and exploit the full maneuvering capability of the missile airframe. Aircraft are less likely to employ acceleration autopilots, unless they are destined to fly low and hug the terrain while avoiding obstacles. Yet, we embedded an acceleration controller inside an altitude hold loop for air corridor flying. Finally, the climb and descent of an aircraft is controlled by the flight-path-angle controller.

The output of these autopilots are sent to actuators, which rotate control surfaces resulting in moments about the vehicle’s c. m. We will now turn to the dynamic modeling of these actuating devices.

Highly agile air-to-air missiles rely on tight autopilot control for successful target intercept. Because they zero miss dis­tance by lateral acceleration, you will find that all modern missiles are controlled by acceleration controllers. To improve performance, an inner rate loop is added for stability augmentation. Both acceleration and body-rate feedback are provided by the accelerometer and rate sensors of the INS.

Airplanes are less likely to be controlled by acceleration autopilots because passenger comfort is more important than maneuverability. However, for special
applications, like unmanned combat air vehicles or cruise missiles, the normal load factor plane may contain an acceleration feedback loop. Such guidance systems as terrain following and obstacle avoidance or target homing require rapid response that only an acceleration autopilot can provide.

Having used classical control for the rate and roll autopilots, I use modern pole placement techniques for designing the acceleration tracking loops. For best performance, proportional and integral (PI) techniques are applied—proportional control for quick response and integral control for zeroing steady-state errors. After a general introduction to state variable pole placement techniques, I will deal with the pitch plane acceleration autopilot and solve for the three feedback gains that satisfy the specified closed-loop response. For missiles the yaw plane implementation is the same.

Let us formulate the problem (see Fig. 10.16). Given is the linear time-variant plant* = F(t)x +g(t)u. What is the controller gain c in order for the closed-loop poles to take their places at the locations pt, і = 1,…,« for all times t; and what is the gain к so that the steady-state gain is one? We consider only SISO systems, with v as input and у as output.

The solution follows the same steps we used for the rate loop, except now we employ the state variable form. The closed-loop system is * = [F(t) — g(t)c]x + kg(t)v with the closed-loop fundamental matrix [F(t)—g{t) c], whose eigenvalues must be equal to the desired closed-loop poles. Therefore, the condition for pole placement is

П

DetjIs – F{t)+g(t)c = ]~[(s – Pi) (Ю.76)

;=i

with n the number of states.

The forward gain is calculated on the basis that у = v when steady state has been reached, which is the case when * = 0. Employing the steady-state closed-loop system equation and the output relationship у = hx, we can solve for the forward gain

1

h[F(t) – g(t)c]-lg(t)

Equation (10.76) is evaluated by equating terms of equal power in 5. This process leads to a system of linear equations, which can be solved for the controller gain c

c = P~’d (10.78)

where P is an n x n matrix and d an n x 1 vector.

For implementation, either in a simulation or onboard a flight vehicle, Eqs. (10.77) and (10.78) are calculated on-line, given the F(t) matrix and g(t) vec­tor. Note that in general two matrix inversions of order n must be carried out. Fortunately, for the acceleration autopilot we can solve Eq. (10.78) algebraically, whereas later on the flight-path-angle autopilot of Sec. 10.2.2.6 requires the full matrix inversion. To adapt the controller gains to highly variable flight conditions, the desired closed-loop poles pt can be made a function of some flight parameter, say, time, dynamic pressure, or altitude.

Let us get specific with the acceleration autopilot. Fig. 10.17 portrays the struc­ture. The proportional and integral feedforward branches with their respective gains Gp and G/ are clearly visible. The major feedback loop is via the rate and acceleration gains c — k], A second acceleration feedback with unit gain

is wrapped around the outside to improve performance. We derive from this schematic the relationship for the control variable

и — — cx + G[ j(v — h%x)dt + Gp(v — h^x)

Then we augment the states by introducing the scalar auxiliary variable § = f(v— h$x) dt with its state equation

I = v — x

Substituting и into the open-loop system yields the closed-loop system, augmented by the auxiliary variable §:

(10.79)

which we abbreviate by the equation

x’ = F'(t)x’ + g'(t)v

The eigenvalues of this closed-loop fundamental matrix F’ must be equal to the

desired closed-loop poles. The condition for pole placement is

Det[/s – F'(t)] = [](s – Pl)

i=i

and specifically

The evaluation of this equation leads to the on-line calculation of the gains c, G and Gp, given the desired closed-loop poles pi.

Now we are ready to design the acceleration autopilot for agile missiles. From Chapter 7 we carry over the linearized plant Eq. (7.69). It is similar to the rate equation, except for the angle of attack, which is replaced by the normal acceler­ation a:

To solve Eq. (10.80), we first establish the following correspondences with the plant equations:

и = S; = [ 0 1]; c = [k.2 k]

Then we apply them to Fig. 10.17 and display them in Fig. 10.18. Recognize the two acceleration feedback loops, the inner rate loop, and the two PI feedforward branches. The rate and acceleration signals are taken from the INS, and the <5

signal is sent to the pitch actuator. The gains k, k2, Gi are calculated from the pole placement condition Eq. (10.80)

Because of the simple structure, it was possible to solve for the gains without a matrix inversion. The position feedforward gain GP is not accessible to the pole placement technique. It must be determined based on a root locus analysis.

In a typical application the dimensioned aerodynamic derivatives are loaded into the missile computer as tables of Mach and angle of attack; the velocity is available from the onboard INS. If the missile experiences little change in dynamic pressure, it may be sufficient to specify fixed closed-loop poles. However, in the case of widely changing flight conditions, the closed-loop poles are scheduled as a function of dynamic pressure. As a rule, it is sufficient to schedule the natural frequency w„ only while keeping the damping £ and the real pole p location fixed. The gain GP is usually fixed and apportions additional feedforward for increased speed of response.

As you check out your design, make sure that the gains are not too high and cause the control to saturate prematurely. Lowering the bandwidth con can alleviate a potential problem. For tetragonal missiles you can use the same acceleration loop for both the pitch and yaw planes. A roll position controller will round out your autopilot design. You will find such an example in CADAC SRAAM6, Module C2 Subroutines C2PI and C2ROLL.

Heading changes of an aircraft are exe­cuted by roll control. As the lift vector is banked, a horizontal force component generates a lateral acceleration that turns the velocity vector horizontally. Direct sideslip control is ineffective because of the small lateral projected area of the aircraft; also, the ensuing adverse yaw-roll coupling would be undesirable.

We build the heading angle tracker by simply wrapping a heading loop around the roll position autopilot. The schematic is shown in Fig. 10.15. To derive the transfer function from the roll angle ф to the heading angle x, we set up the relationship between the lateral acceleration and the lateral force, caused by the lift vector, banked through the angle </>; and with the small angle assumption of ф, we can

simplify the expression

m Vj( = L sin ф = mg4>

and obtain the desired transfer function

X(s) _ J_

</>(5) Vs

shown in the last block of Fig. 10.15. The gain Kx translates the heading error Xc — X int0 the roll command фс after limits are applied.

We employ again the pole placement technique to calculate this gain. Clinching the open-loop transfer function from Fig. 10.15 and with Eq. (10.74)

t / „ ^ 8

*0S) — Л – V, – у у

fsJ + 2t;consl + co^s

and closing the feedback loop, while setting it equal to the generic third-order transfer function,

T(s) = Kx(8/VH = Kx(g/V)co2n

s3+ 2£wns2 + w2s + Kx(g/V)w2 (s + px)(s2+ 2$xconxs+co2x)

with the free parameters px,£x, and co„x, which determine the modal response of the heading tracker. Equating equal powers of s yields three equations, but, including Kx, there are four unknowns. Hence, one additional equation is required. After having studied the root locus for some time, I found that positioning px on the real axes, directly under the complex poles of the roll autopilot, grants a well – behaved heading response. Thus, we set

Px =

Eliminating and co„x, we calculate the gain from

Kx = -?<o„(l -£2) (10.75)

g

The heading gain depends on the data entry f and con of the roll autopilot and the flight speed V. If the optimal f = 0.707 is selected, the gain is simply Kx = 0.147VoL>n/g. Just as the roll autopilot, the heading loop is also a perfect tracker, thanks to the integrator in the forward loop. You can find the implementation of the heading tracker in the CADAC GHAME6 simulation, Module C2, subroutine C2HEAD. It operates in conjunction with subroutine C2ROLL.

The roll position autopilot has a dual purpose. Missiles, maneuvering independently in the pitch and yaw plane, use the roll autopilot to suppress any roll excursions that may be caused by aerodynamic cross coupling. Although earlier designs, like the Sidewinder air-to-air missile, incorporated only roll damping, all modern missiles possess roll position loops. If the airframe lacks tetragonal symmetry—like cruise missiles or airplanes, the roll loop is used to bank the vehicle for lateral maneuvers.

We build a dual feedback controller, utilizing the roll rate from the INS gyro and roll position from INS navigation computations. The inner rate loop augments the aerodynamic damping and the outer position loop executes the roll command. The transfer function is rather simple, consisting of the roll damping derivative LLp = (qSb/Iu)(b/2V)Cip and roll control derivative LLSa = (qSb/I\)Cis [see Eq. (7.57)]

function and set it equal to the desired form that introduces the parameters £ and con:

(pis’) KfpLL$a

<pc(s) s2 + (KPLLS„ – LLp)s + КфИца s2+ 2t;(0ns + (О2 Comparing terms of equal power of 5, we can calculate the two gains

The comments I made about rate loops apply here as well. To prevent aileron sat­uration in extreme flight conditions, you should schedule £ and con and calculate derivatives on-line from stored data. You will be glad to know that, thanks to the integrator in the forward loop, the steady-state error of the roll loop is zero, i. e., the roll position autopilot is a perfect tracker. An implementation of the roll autopilot is given in the CADAC GHAME6 simulation, Module C2, Subroutine C2ROLL.

The most elementary autopilot function is the improvement of the dynamic stability of the vehicle, both in the pitch and yaw planes. This is accomplished by augmenting the aerodynamic damping derivative by rate feedback loops. The feedback sensors are the rate gyros from the INS system. We shall have to deal with the pitch and yaw planes separately. The pitch – rate autopilot is the same for both aircraft and missiles, although we may have to use different symbols, like lift slope Cia and the normal slope C^. The first one is given in velocity axes and the second one in body axes. Missiles prefer Сца, whereas aircraft use the more traditional Cia. If the missile has tetragonal symmetry, the yaw aerodynamics behaves like those of the pitch plane, and therefore the same autopilot serves both planes. In the aircraft case, however, the lateral aerodynamics is distinctly different and requires a separate treatment.

Let us apply the pole placement technique to this simple example. We employ the open-loop transfer function of the pitch plane, close the loop with a gain, and calculate the closed-loop poles. Because the gain is a free variable, we can choose it such that the closed-loop damping achieves a desired value for all flight conditions.

From Sec. 7.4.2 we carry over the linear perturbation equations of the missile pitch plane in the two state variables x = [q a], pitch rate and angle of attack. (Note that in this autopilot section we adopt the bolded vector convention of modern control.)

where Na = (qS/m)C^a, Ns = (qS/m)CNs, Ma = (qSd/l2)Cma, Mq = [qSd2/ (2/2 W)JC,„;/, Ms = (qSd/l2)Cms are the dimensional derivatives for the normal

force N and the pitching moment M <5 is the pitch control surface deflection; and V the flight speed. These plant equations are in the general form of

x = Fx + gu

Now we build the transfer function of pitch rate q with respect to pitch control <5. Take the Laplace transformation of the plant equation and solve for x(s)

sx(s) = Fx(s) + gu(s)

(si – F)x(s) = gu(s)

x(s) = Ф (s)gu(s)

where Ф(5) = (si — F)~l is the resolvant matrix. The input u(s) is <5, whereas the output _y(5) could be either q or a. The vector h(s) = [10] will pick the output q from the general relationship

y(s) = h<b(s)gu(s)

Therefore, the desired transfer function is

Ф)=їіф( ) Ms[s +NJV -(Ma/MS)(NS/V)] S(s) {S)8 s2 + (Na/V – Mq)s – Ma – MqNa/V

You can verify this equation by comparing it with Roskam (Ref. 16, p. 426). The transfer function is positioned in a simple feedback loop (Fig. 10.11) and has two complex conjugate poles and a zero with the root locus pattern shown in Fig. 10.12. Selecting the closed-loop damping coefficient f determines the gain К.

To compute this relationship, we abbreviate the transfer function by

q(s) _ G(s + z)

S(s) s2 + as + b

and build the closed-loop transfer function

q(s) _ G(s + z) _ G(s + z)

q,(s) s2 + (a + KG)s + b + KGz, s2 + 2t;cons + со2

The last formulation introduces the desired closed-loop parameters £ and wn. Equating terms of equal powers in the denominator yields the formula for the rate loop gain

where the negative sign before the radical was discarded because it does not lead to a useful solution. The resulting natural frequency of the closed-loop poles, which follows from the selection of £, becomes

co„ — Vb + KGz

No matter what the flight conditions are, within the limitations of linear approx­imations the damping of the closed-loop response is guaranteed to be the value that you select. However, it could come at a penalty. A high gain K, causing a high natural frequency wn, may drive the controls prematurely into saturation. You need to check out the control deflections for some of the more extreme flight conditions before you commit to a desired f. If you cannot reconcile all conditions with one value £, you could make it dependent on dynamic pressure.

The calculation of A" in a simulation and in an actual flight control system is carried out on-line, making use of sensors and stored data. Tables provide the aerodynamic data and mass properties, while Mach number and dynamic pressure come from the air data computer. Combined, they establish the transfer function parameters. Together with the required damping f, the gain К can be computed.

This simple rate loop is used in the CADAC SRAAM6 simulation, Module C2, Subroutine C2RATE. It works quite well during the initial phase, as the missile leaves the launch rail, requiring stabilization. Because the SRAAM6 has tetragonal symmetry, both pitch and yaw channels have identical rate loops.

The GHAME6 simulation of the hypersonic SSTO vehicle requires a yaw sta­bility augmentation system (SAS). Our rate loop satisfies this need, yet must be adjusted to airplane yawing conventions. A positive /3 generates a negative side force Y, whereas a positive a causes a positive normal force N. Accounting for this difference, the linearized perturbation equations in the yaw plane are for airplanes

where Yp = (qS/m)CYl,,Ys = (qS/m)CYi, LNp = (qSb/hi)Cnil, LNr = [qSb2/ (2/33 V)]C„,, LNj = (qSb//зз)СЛі are the dimensioned derivatives for the side

force Y and the yawing moment LN <5 is the rudder deflection. According to the pole placement procedure, we obtain the open-loop transfer function, again confirmed by Ref. 16, p. 458:

Ф) = LNs[s — Yp/V + (LNp/LNs)(Ys/V)]

<5(s) s2-(Yp/V + LNr)s + LNp + LNrYp/V ( ‘

and with the abbreviated closed-loop transfer function

rjs) _ ________ Gjs + z)_______ _ GQ + z)

rc(s) s2 + (a + KG)s + b + KGz s2 + 2£ cons + w2

We can use Eq. (10.72) again to calculate the yaw loop gain K. You can find an example for a yaw loop autopilot in the code CADAC GHAME6, Module C2, Subroutine C2YSAS.

The aerodynamics, propulsion, and mass properties, activated by Newton’s and Euler’s equations, model the vehicle dynamics. They define the plant to be con­trolled. An aircraft needs rate damping for stability augmentation, roll control for coordinated turns and flight-path-angle tracking. For long duration flights, altitude and heading hold autopilots are indispensable for a weary crew. Missiles, flying in the atmosphere, are primarily controlled by rate and acceleration feedback loops. A roll controller maintains their orientation. For longer duration flights the trajec­tory can be shaped by a flight-path tracker or altitude hold autopilot. Outside the atmosphere missile attitude is controlled by angular feedback loops.

The feedback signals come from sensors that are located throughout the vehicle. The INS, an indispensable component in any modern aerospace vehicle, provides rate and acceleration feedback for all three axes, directly from its accelerometers and gyros, and attitude angles from its navigation computer. If wind is not a factor, it can also calculate vertical and horizontal flight-path angles. The air data computer provides Mach number and dynamic pressure. Using incidence angles as feedback signals is problematic because it is difficult to accurately measure or compute the angle of attack or sideslip angle. The alpha vanes that you see on the outside of airplanes are only used as warning devices and not as autopilot sensors.

The coverage in this section picks up from Sec. 9.2.3 where I discussed autopi­lots of pseudo-five-DoF simulations. As you may recall, in five-DoF simulations autopilots do not issue control commands to the control surfaces but model closed-loop response, delivering angle of attack to the trimmed aerodynamic ta­bles. In the real, nonpseudo world the autopilot, or the human operator, controls the aileron, elevator, and rudder directly, and only through the airframe and the sensors are the dynamic loops closed.

Most of the autopilot designs of Sec. 9.2.3 do not apply here, with the exception of the altitude hold autopilot. Its outer altitude position and rate loops can be maintained as long as the inner acceleration autopilot is replaced by its six-DoF equivalent. So let us get a new start and build control systems for aircraft and missiles that are suitable for our six-DoF simulations.

Before we begin, let me caution you, however, that my autopilot designs are simplified versions of what you actually find onboard these vehicles. Aircraft flight control systems can be very sophisticated and are beyond the scope of this discussion. If you are tasked to model an existing vehicle, you should get the autopilot specification and replicate it as faithfully as possible. Only then will the customer trust the results of your analysis. However, if your job is to study a generic system or to develop a new concept, you will find the following designs useful.

To make the autopilots serve a variety of vehicles and flight conditions, I employ the pole placement technique. You can find this method discussed by Stevens and Lewis.9 Yet if you have some familiarity with classic and modern linear control, you should be able to follow this self-contained presentation.

You may be surprised, but hypersonic vehi­cles are modeled just like airplanes. The framework of Lanchester10 applies even here. However, to gather the aerodynamic data is a gargantuan task. The greatest challenge is probably the SSTO vehicle. It takes off at low subsonic speeds, tran­sitions through the sonic barrier to the supersonic regime, and, after about Mach 5, enters the hypersonic region. Eventually, the air density becomes so small, at about 100 km, that the aerodynamic effects are essentially nil.

It would be nice to have the resources to gather the data in wind tunnels. It takes three different facilities (subsonic, supersonic, and hypersonic) to cover the flight regimes. Yet, even then one has to make compromises. It is difficult to match the wind-tunnel’s Reynolds number with free flight, and the Mach number tops out at about 10. You may take recourse to the highly touted computational fluid dynamics Mafia. However, their supercomputers devour your budget as well, and they will still send you back to the transonic tunnel to fill in the gaps.

There is a third option available to you. Since the dawn of the computer age and the maturing of the FORTRAN language, aerodynamicists have condensed their knowledge in semi-empirical codes. We already have met the Digital DATCOM

3B

Fig. І0.8 Hypersonic vehicle lift and drag coefficients.

for aircraft and the Missile DATCOM programs. Both operate from dimensioned vehicle data using analytical formulas and interpolated wind-tunnel databases. In the early 1970s, so-called panel methods were developed, which are particularly useful in the supersonic and hypersonic regimes. Best known is the U. S. Air Force Supersonic/Hypersonic Arbitrary Body Program (S/HABP),12 which is very popular with industry. It was upgraded in the late 1980s by McDonnell Douglas, St. Louis,13 and is based on the impact pressure calculation method, which is most accurate at hypersonic speeds. Its capabilities, including viscous effects, extend down into the supersonic regime using so-called embedded flowfields.

Another semi-empirical code, the so-called PAN AIR method,14 is NASA’s con­tribution to supersonic and hypersonic aerodynamic prediction. It solves the linear partial differential flow equations numerically by approximating the configura­tion’s surfaces by panels. PAN AIR is a higher-order panel method that is limited to inviscid flow.

NASA’s Generic Hypersonic Aerodynamics Model Example (GHAME) serves us as a test case.15 It is a hypersonic vehicle that takes off horizontally under turbojet propulsion. Once at Mach 2 the ramjet propels the vehicle to Mach 6, at which the scramjet takes over until the limits of the atmosphere are reached. The aerodynamic model was developed with the S/HABP code and from experimental data based on the space shuttle and the experimental vehicles X-24B and C.

Contrary to the FALCON6 model, GHAME aerodynamics uses lift and drag as force components, as depicted in Fig. 10.8. Not shown are the side force and the moment conventions because they remain unchanged. The NASP-like vehicle is controlled by rudder and two elevons: Svi, left and 8vr, right. The elevons function as ailerons 8a and elevators 8e according to the relationships

The aerodynamic forces are given in stability axes (see Sec. 3.2.2.5)

[T,S=qS[-CD CY -CL

whereas the moments are in body axes as before

[mB]B = qS[Ci x b Cm x c C„ x b]

NASA formulated the aerodynamics by simple Taylor-series expansion, includ­ing only the linear terms and making all derivatives tabular functions of Mach number and angle of attack:

Cl = CLll(M, a) + CLJM, a)a + CLJM, a)8e + CL (M, a)q

Cd — Cd0(M• «) + CDam, a)a

CY = CYn(M, a) + CYll(M, u)P + CyJM, a)8a + CYil (M, a)Sr

Yet, by no means do these equations represent a linear model. Look at the lift slope derivative C/.u(Af, a). It is a tabular function of a and is multiplied again by a. The same is true for the drag derivative CDa(M, a). The other derivatives are linear in terms of the remaining states, but still nonlinear in Mach and a. The effect of body rates on the lift and side forces, i. e., the Cl, CY, and CYr derivatives, can often be neglected.

The moment coefficients follow a similar pattern:

C, „ = a) + C„,a(M, a)a + CmJM, a)8e + Cmq(M,

C, = Ci„(M. a) + С/р(М, a)P + a)8a + Chr(M, a)8r

+ Clr(M, a)p^+Cli(M, a)rA (10.66)

c„ — C„0(M, a) + C„,(M, a)P + CnJM, a)8r + C„fo(M, a)Sa + C„p(M, a)p^ + Clh(M, a)r^

All four terms in the pitching moment coefficient Cm are important. Yet, the rolling and yawing moment coefficients С/ and Cn usually have negligible trim coefficients C/„ and C„„. In effect, Fig. 7.1 predicts that they are zero. Notice the strong cross coupling between the rolling and yawing moments. It occurs through the famous Dutch-roll derivative and the body rates p, q and the controls 8a, Sr.

We are fortunate that NASA has provided us with complete tables. If you want to get a feel for their numerical values, you can turn to the CADAC GHAME6

simulation and consult Module A1. Now, finished with these aircraft-type vehicles, we turn to missiles with tetragonal symmetry.

10.2.1.3 Missiles. Most missiles have simple geometrical shapes. In rock­etry the four fin, axis-symmetrical configuration, first fashioned by Wernher von Braun as the V2, are prevalent still today. Air-to-air missiles like ASRAAM and AMRAAM exhibit this tetragonal symmetry, as well as our CADAC SRAAM model.

The tetragonal symmetry fosters the expression of the aerodynamics in aerobal- listic axes. Refer back to Fig. 3.18 for the definition of the total angle of attack a’ and aerodynamic roll angle ф’. Both replace the standard incidence angles a and fi. The aerodynamic forces and moment are expressed in coefficient form in aeroballistic coordinates ]R

[‘I, f = qS[-cA c’y – c;v]

W~B]R=qSi[C, c; C’]

where q is the dynamic pressure, S the reference area (cross section), / the reference length (diameter), and the prime refers to the aeroballistic axes. Figure 10.9 depicts the positive sense of these coefficients. Pay particular attention to the positive direction of the aerodynamic roll angle ф’. It is taken positive from the aeroballistic axes to the body axes. Also notice that the axial force and the rolling moment coefficients are invariant under this transformation.

The functional dependency of a given configuration is reduced to the following form:

Ca = f(M. Sq’. Sr’. power) С/ — f(M. а’, ф’, p, Sp)

C’Y = ДМ, а’,ф’,8г’) C’m =/{М. а’.ф’. c. g., q’.Sq’)

C’v = f(M. а’, ф’. Sq’) C’n = f(M, а’. ф’, c. g., r’, Sr’)

where the Reynolds and incidence rate dependencies were neglected.

The body rates p. q’. r’ and the control commands Sp, Sq’. Sr’ are expressed in aeroballistic coordinates ]R and obtained from body coordinates by the

P

q’

r’

Sp

Sq’

Sr

The control commands require further explanations. The starting point is Fig. 10.10 with the convention of positive surface deflections. For reinforcement I have also included again the positive direction of the aerodynamic roll angle ф’. Because there are four fins, but only three attitude degrees of freedom, we combine them mathematically to form the three control commands

Sp = ±(-<51 -82 + 53 + <54)

Sq = ±(<51 +52 + 53 + 54) (10.68)

Sr = ±(-<51 +52-53+54)

A fourth combination is also possible, but it produces only drag and no control moment. Look at Fig. 10.10 and work out the three relationships. Four positive fin deflections create a positive pitch command; the first two fins up and the last two fins down generate a positive roll command; and alternating the fin deflections causes a positive yaw command. This convention is in agreement with aircraft practices and is summarized in the following equations.

Roll:

+8p —»■ +Д L rolling moment

Pitch:

+Sq —*■ +AN normal force

Yaw:

+<5r —>■ +ДТ side force

Now we are ready to expand the coefficients. Please refer to Eq. (10.69). The axial force coefficient CA consists of the base and skin friction drag СЛо and the power on/off effect AC a- Furthermore, it has a first-order term in a’ and a second – order three-dimensional table look-up in <5eff = (8q + <5r1 )/2 fthe effective fin deflection). You can see that great care has been taken to model the axial force because its fidelity determines the fly-out range of the missile. The side force and normal force coefficients C’Y and C’N are dependent on the roll orientation of the missile relative to the load factor plane, described by the angle 0′ (ф1 = 0 when the missile maneuvers in the cruciform configuration):

Ca — CAo(M) + AC,4(p0wer)(M) + САиЛМ)а’ + C/sj2rf (M, a’, <52ff)<52ff

Су — АСуф,(М, a’) sin4<£’ + C’Yg (M, a’)8r’ (10.69)

C’N = Сд,„(М, a’) + ДСд,_0.(М, a’) sin2 2ф’ + C’V^(M, a’)8q’

Inspecting Eq. (10.69), you see that this dependency is modeled by sine functions. Both terms, sin 40′ and sin2 0′, are periodic every 90 deg, corresponding to the tetragonal symmetry of the missile. However, the increment in side force reaches its maximum at 22.5 deg and its minimum at 67.5 deg. The normal force effect is always positive and reaches its peak every 45 deg. Realize that in aeroballistic axes the side force component is always small because, by definition, the missile responds in the load factor plane. Any perturbations caused by rotational asym­metries are of second order and are many times just neglected. The other terms in Exj. (10.69) are the control effectiveness derivatives Cy_ and C’N.

The moment coefficients are modeled in a similar fashion. They consist of primary terms, damping derivatives, control effectiveness, and c. g. adjustments. The primary roll term C;,^ ,(M. a’)a’2 sin 40′ of Eq. (10.70) is an attempt to model the roll coupling caused by vortices impinging on the tail surfaces at high angles of attack. It is dependent on the roll orientation of the missile, sin 40′, which we encountered already in the C’y coefficient:

Cl = CI<2(M, a)aa sin40′ + С/ДМ, or’)^ + Chp(M, a’)8p

+ a’)5q’- -f(*cg, fi – Mg)

This yaw-roll coupling is 90-deg periodic in <p’, causing a positive or negative rolling moment, depending on the orientation of the fins relative to the load factor plane. The other two terms of the rolling coefficient bear no surprises.

The pitching moment coefficient, similarly affected by the orientation of the load factor plane, possesses a perturbation term ЛС’г ф,(М, a’) sin2 2ф’. It vanishes at 90-deg increments between maxima. The yawing moment coefficient C’ repeats the pattern of С/.

To adjust for weight shifts during flight, the pitching and yawing moment must be corrected for the shift in c. g. Just like the aircraft equations (10.62), the adjustment is made with the moment arm (xcg д> — xcg), where лС1, д> is the reference location (at launch or from wind-tunnel tests.) and xcg the true location.

I have presented to you just one missile model, which is quite useful for pre­liminary performance studies. It has been adopted for many air-to-ground and air-to-air simulations and linked with the Digital DATCOM aerodynamic predic­tion program. If you want to see it in action, you should refer to the CADAC SRAAM6 simulation, study the Module Al, and run the test cases.

10.2.1.4 Summary. The three aerodynamic models, although limited in sophistication, should give you an appreciation of the six-DoF treatment of aero­dynamics. There is a distinct difference between aircraft and missiles, correspond­ing to their respective symmetries. Planar vehicles are modeled with the Cartesian incidence angles a and /3, whereas tetragonal missiles lend themselves to the treat­ment with polar incidence angles a’ and ф’. However, you will find many missile simulations that use also the Cartesian incidence angle scheme. Next time you get your hands on a six-DoF simulation, see if you recognize some of the modeling features of this section.

In the preceding section we concentrated on the left side of Newton’s and Euler’s equations. Now we move over to the right-hand sides and discuss the aerodynamic forces and moments. The modeling of aerodynamics is more of an art than a science because the physics of the airflow over a vehicle is so complicated that even today’s supercomputers have to bow to the supremacy of the wind tunnel. However, the art of modeling aerodynamics is well understood, and we make use of the extensive background material available to us.

Before we proceed, I encourage you to review the five-DoF aerodynamics sec­tion 9.2.1 and the Taylor-series development in Sec. 7.3. For six-DoF modeling the book by Stevens and Lewis9 provides additional up-to-date insight into the linearization process.

Let us revisit the functional dependencies of the aerodynamic forces and mo­ments of Sec. 9.2.1:

aero forces and moments

M, Re, a, p, a, p, p, q, r, 8p, 8q, 8r, c. g., power, shape, scale

^ flow incidence angles body rates control surface

characteristics and rates deflections

where the flow characteristics are determined by Mach number M = velocity/sonic speed and Reynolds number Re = inertia forces/frictional forces. The incidence angles, here in Cartesian form, are the main dependencies. Their derivatives are usually of secondary importance and sometimes combined with the body rates, although their physical origins are quite distinct. The control surfaces are the main moment generators and are crucial for controlled flight. If the vehicle bums much fuel, the c. g. is bound to shift and has to be accounted for. Equally important is the effect of the exhaust plume on the vehicle drag. Hence, power-on/off effects must be included. Finally, the geometric descriptors are shape and scale, but are constant for a given vehicle.

You can imagine there are many techniques to simplify the functional form. The two primary methods are 1) the brute force modeling by tables and 2) the expansion in Taylor series. Diversity arises from the mix of the two approaches. Rather than being all inclusive, I will select three types of vehicles and furnish their aerodynamic models. The FALCON6 represents the aircraft model, NASA’s GHAME6 vehicle exemplifies the hypersonic flight regime, and tactical missiles are embodied by the SRAAM6 concept.

10.2.1.1 Aircraft. Aircraft aerodynamicists, like Lanchester,10 have pione­ered the modeling of aerodynamics. His system of aerodynamic coefficients has served us well over the years. Even rocketeers, after first using ballistic coefficients, have joined up and have used since the early 1960s11 the same framework.

The six-DoF aerodynamic model of an aircraft can be as simple as the linear terms of the Taylor-series expansion or as complicated as tables with five inde­pendent variables. You will find the simple representations in simulations that emphasize autopilot design studies, whereas flight simulators require high-fidelity aerodynamic response for pilot training. We will take the middle road and essen­tially follow Stevens and Lewis.9

Today, the formatting of force coefficients is mostly in body coordinates with the positive sense following the direction of the body axes

[7a?=qS[Cx CY Cz]

where q is the dynamic pressure and 5 the reference area. Gone are the times of using lift and drag coefficients, although we shall encounter this archaic form shortly at the hypersonic vehicle. The format of the moments, referenced to the c. g. B, is noncontroversial and has always been in body coordinates

ImB]B — qS[C/ x b C„, x c C„ x b]

with the reference lengths b as span and c as chord.

Figure 10.7 should help you visualize the positive directions of the force and moment coefficients. The controls of an aircraft are aileron 8a, elevator 8e, and rud­der 8r. Their positive deflections are defined by the following equations. Whereas a positive aileron deflection causes a positive rolling moment, elevator and rudder deflections are defined as positive if they produce positive increments of lift and side force.

Aileron:

+ALL rolling moment

+<5r -* +ДТ

Now let us expand the aerodynamic force coefficients into a Taylor series but only with respect to the body rates p, q, r and the control deflections aileron, elevator, and rudder, Sa, Se, Sr. The coefficients remain in general a function of Mach M, angle of attack a, and sideslip angle /3. With some variations the following model has been used for the F16 (Ref. 9). Please refer to Eq. (10.61). The axial force coefficient Cx essentially represents drag. Its dependency on Mach number and angle of attack is in tabular form. In addition, the pitch-rate elevator deflection is also tabular because of its nonlinear behavior. Only one Taylor-series expansion term relative to pitch rate is included. It models the drag effect, caused by the increased local angle of attack on the tail, as a result of the pitch rate. Because the side force Су remains usually small, it contains only tables of Mach numbers with the exception of the damping derivatives CYr and CYp, which are also a function of a; otherwise, only the linear terms in the Taylor series of /3, Sa, Sr are taken into account. The normal force coefficient Cz is similar to Cx, except that its dependency on Se can be linearized:

Cx = Cx„(M, a,8e)+^Cx,(M, a)q CY = CYf(M)P + CYJM)Sa + CYJM)Sr

(10.61)

Cz = CZ„{M, a, (3) + CZse(M, a)8e + ^CZq(M, a)q

The rolling and yawing moment coefficients С/ and C„ of Eq. (10.62) are mod­eled in the same fashion with respect to their nonlinear behavior in Their de­pendencies on aileron and rudder deflections have been linearized. In contrast, the
pitching moment Cm is not a function of /3. but the elevator deflection is included in the Cm„ table. An important effect has the center of gravity (c. g.) location on the pitching and yawing moments. Just consider the shift of fuel or cargo during flight. Furthermore, in wind-tunnel testing the moment center of the wind-tunnel model was probably not placed at the c. g. location of the full-up aircraft. In either case, be it dynamic or static, a moment arm xcgR — xcg appears, which couples Cz into C„, and Су into C„. The fixed reference point is at xcgR, whereas the actual c. g is at xcg.

C, = C,„(M, a. P) + QJM, a. P)Sa + Chi (M. a, P)Sr + a)r + Clp(M, a)p

C,„ = a. Se) + ^С,„ч(М, a)q + ^-(xcgR – xcg) (10.62)

C„ — C, h (M, a, P) + C„Sa(M, a. P)Sa + C„Sr(M, a, P)Sr Су b r

– —Ucgr ~ Xcg) + — [C„r(Af, a)r + C„p(M, a)p

Stevens and Lewis9 provide numerical values for the coefficients at the single Mach number of 0.6. You can refer to the aerodynamic tables and mass properties in Appendix A of their book or can go to the CADAC FALCON6 simulation, Module A1.

I had difficulties locating six-DoF aerodynamics of modern aircraft in the open literature. Many of the models are either classified or considered proprietary by the manufacturer. If you have access to flight simulators, ask the software engineer to show you the aerodynamic model. You will be surprised by its sophistication. Here, we are contented with a middle-of-the-road approach.

The equations of motions are the most important part of a simulation, but by no means the most complicated ones. Their models are already in mathematical form and easily programmed. However, an aerospace vehicle consists of many more subsystems encompassing numerous technical disciplines.

The aerodynamicist provides the mathematical model of the aerodynamic forces and moments that act on the vehicle; thrust tables are generated by the propulsion expert; the control engineer designs the flight control system; and INS, sensors, and seekers are modeled by the appropriate experts.

We first discuss the modeling of the aerodynamics of aircraft, hypersonic vehi­cles, and missiles. The approach differs from the five-DoF implementation because now the aerodynamic forces and moments are untrimmed and depend on the con­trol surface deflections. I shall acquaint you with the aerodynamic models of a typical aircraft, like the Fighting Falcon, a NASA hypersonic vehicle, and an air- to-air missile. Each model represents a different approach. The first two, with their planar airframes, use the Cartesian incidence angles of angle of attack and sideslip angle. The missile is modeled with the polar incidence angles of total angle of attack and aerodynamic roll angle.

For propulsion you can reach back to Secs. 8.2.4 and 9.2.2, where I discussed rocket motors, subsonic turbojets, ramjets, and scramjets. Those models can also be used in six-DoF simulations. If you want to build a detailed engine model however, you need to consult with propulsion experts and the references that they suggest.

In six-DoF simulations autopilots and flight control systems are represented by their full-up designs. No simplifications are required as in the pseudo-five-DoF simulations. When introducing you to six-DoF autopilots, I provide fairly simple
structures. They are based on pole placement techniques that have the advantage of adapting to varying flight regimes.

The actuators execute the autopilot commands and turn them into control surface deflections. For aircraft we deal with aileron, elevator, and rudder actuators and in some cases with elevons and rudder. Missiles most likely have four control surfaces with four actuators. Their autopilot commands roll, pitch, and yaw, and the signals are divided up for the four surface actuators. Besides aerodynamic control, I also give you a model of thrust vector control for agile missiles.

Then we take on the subsystems of navigation and guidance. The INS is an inte­gral part of any modern aerospace vehicle. I derive the equations of space stabilized and local-level systems and provide you with a nine-state error formulation.

Earlier you were introduced to proportional navigation. I expand on that discus­sion and derive the compensated PN and advanced guidance law (AGL). Whereas PN was formulated in inertial LOS rates, AGL uses differential velocity and dis­placement, expressed in inertial coordinates, as navigation input.

Finally, I augment the seeker discussion of Sec. 9.2.5 with an imaging IR seeker. We will investigate the intricate modeling of focal plane arrays, and mechanical and virtual gimbals. We also address noise and biases that corrupt the signals.

These subsystems are intended to serve as examples for modeling important components of aerospace vehicles. Although complete in themselves, they do not represent the most sophisticated models that you will encounter in six-DoF simulations. Regard them as instructional material that will get you started, guiding you from concepts to mathematical models and code implementation. You will find most of the subsystems, discussed in this section, coded up in the six-DoF simulations FALCON6, GHAME6, and SRAAM6, available on the CADAC CD.

10.1.2.4 The geo­graphic modeling of the equations of motion over an elliptical Earth focuses on the integration of the geographic acceleration DEvER relative to the Earth rather than the inertial acceleration D’v’B. Make no mistake, the J2000 is still the inertial frame for Newton’s and Euler’s laws. We just introduce the Earth as a convenient reference frame. I will derive the relationship between the two velocities and their accelerations before presenting the translational equations. As before, the attitude equation holds no surprises.

There are two reasons for focusing on the geographic velocity vf: 1) the at­mosphere is rotating with the Earth, and the aerodynamic forces are a function of the geographic velocity; 2) the vehicle’s position and velocity are traced more conveniently over the Earth than in the inertial frame.

By definition, the inertial velocity is v B = D’sbi, where s bi is the location of the vehicle c. m. В wrt the center of the Earth /. To introduce the geographic velocity, we change the reference frame from I to E

D! sbi — DeSbi +

Now, introduce a reference point on the Earth E (any) as an intermediary between В and /:

вві = sbe + sei (10.50)

and substitute it into the preceding equation

D1sbi = Desbe + Desei + Qeisbi

DesEi = 0 because sei is constant in the Earth frame. Using the definition of geographic velocity vB = Desbe, we get the desired relationship between the inertial and geographic velocities:

v’B =vb+^EIsbi (10.51)

The translational differential equation is derived from Newton’s second law

mDlv’B = f a p + mg

It must be modified to bring out the geographic velocity. We first deal with the rotational derivative and transform it to the frame E

D’v’g = DEv’B + nE’v’B (10.52)

then we modify the first term on the right-hand side by substituting Eq. (10.51)

DEv’g=DEvE + DE(nE’sBI)

and subject the last term to the chain rule

De(EIe’sBi) = DEnEIsH, + EIeiDesBi

The first term on the right-hand side is zero because the Earth’s angular velocity is constant. The second term can be expanded:

nEIDEsBi = ftEI(DEsBE + Desei) = CIeiDesBe = ПЕІуЕв

Collecting terms

Dev’b = DEvEB + SlE’vEB

substituting them into Eq. (10.52) yields

D! v!b = Deveb + SlEIvEB + nE! v!B

and finally replacing the last v! B by Eq. (10.51) we get a relationship between the inertial and gravitational accelerations

D‘v’B = Deveb + 2flEIvEB + flEInE! sBI (10.53)

We identify the Coriolis term 2flEIvB and the centrifugal term flEIflEIsBi that we first encountered in Sec. 5.3.1. Replacing the inertial acceleration in Newton’s equation by Eq. (10.53), we obtain the translational equation

DEvE = – f +g – 2SlE, vEB – ПЕ! ПЕІ*ВІ (10.54)

m

It is exclusively a function of the geographic velocity and acceleration, although it governs the translational equations of motion relative to the inertial frame. It is also still in the invariant tensor form and therefore holds in any coordinate system. In particular, if we chose a coordinate system that is associated with the Earth’s reference frame E, the rotational time derivative will turn into the ordinary time derivative. We have three choices; the Earth ]£, the geocentric ]G, or the geodetic ]° coordinate systems. Most likely, we want the geographic velocity computed in geodetic coordinates. In this case we program

^11 = i[f]BD[f f + [Г]ОС[ ]G

at J m

-2[T]d,[£2ei]![T]di[uf]D – [7′]°/[^£/]/[^£/]/[sb/]/ (10.55)

Each term on the right-hand side is expressed in the most suitable coordinate sys­tem: the forces in body coordinates, Earth’s gravitational acceleration in geocentric coordinates, Earth’s angular velocity in inertial coordinates, and the geographic velocity in geodetic coordinates. The vehicle’s displacement is referenced to the center of Earth [.vb/]1 and is shown in inertial coordinates. It could as easily, with the help of Eq. (10.50), be expressed in geodetic Earth-referenced position [sbe]1* (point E can be any reference point fixed in the Earth frame)

[SB/]7 = [T]d,[sbe]D +

The transformation matrices [T]DG and [T]DI have already been given, where [T]BD is obtained from the direction cosine matrix calculations.

An alternate form of Eq. (10.55) is also useful. It stems from the require­ment of calculating the incidence angles a and /8. They are derived from the geographic velocity expressed in body axes [uf]B. We could of course obtain it from [u|]B = [T]BDvf ]°, but some aerospace professionals prefer to formulate the geographic acceleration in body coordinates proper. In that case we take the perspective of the geographic velocity from the vehicle’s body frame and make an Euler transformation to body frame В

Substitution into Eq. (10.54) yields

Dbveb + nBEvEB = – f +g- 2nE! vEB – nE! nE1sBI

m

Expressed in body coordinates the rotational derivative will convert to the ordinary time derivative, and Newton’s equations are programmed as

= -[£2B£]B[uf ]B + -[fa, p]B + [T]BG[gf L J m

– 2[rf[n0];[f ]BI[vEB]B – [Г]в/([^£/]/[^£/]/[^/]/) (10.56)

where [QBE]B = [^й/]й – [Т]Ш[ПЕІУ[Т]В!, [QBI]B is obtained from Euler’s equation, and [T]BI from the direction cosine equation. Equation (10.56) is also found in this form in the book by Stevens and Lewis [Eqs. (1.3-1.6)].9

Having obtained the geographic velocity vf defined as vf = Desbi from inte­grating Newton’s equation, a second differential equation solves for the inertial position sB1:

D! sBi = DesB! + flEIsBI

Note that the inertial point I in DEsB; is located at the center of Earth and belongs also to the rotating Earth frame E. For this reason the geographic velocity can be defined two ways: vf = DesBe — DEsBj. However, let me caution you that the point E in this context can only be the center of the Earth and not just any point of the Earth frame, as was the case earlier. With this stipulation we formulate the second translational equation:

£>7sB/ = vf + flEIsBI (10.57)

For programming purposes we have to distinguish whether the geographic velocity is given in geodetic or body coordinates. Accordingly, we have two options:

= [Г]°/[иІ]° + [П£/]/М/ (10.58)

= [L]B/[uf]5 + [f2£/]/M/ (10.59)

Equations (10.55) and (10.58) constitute the translational equations with emphasis on geodetic coordinates, whereas Eqs. (10.56) and (10.59) focus on the body coordinate system. Both give the same results. It is up to your customer to express her preference. The initialization of the geographic velocity v f may be the deciding factor. [vB]D is easier to visualize than [uf ]B. Yet in any case, you will need to calculate both: [vB]B to derive the incidence angles and [uf ]D to plot the velocity components over the Earth.

The derivation of the attitude equation harbors no surprises. Euler’s law, referred to the inertial frame 7 and expressed in body axes, leads to the same rotational

equations (10.44), already used and repeated here for completeness:

= ([/l]Br1Hn"]B[/2]lV/]B + -тв}в) (10-6°)

The direction cosine matrix is as already calculated from the inertial body rates [coBI]B, Eqs. (10.47-10.49).

Let us list the state variables that are to be integrated: the geographic velocity [uf ]D or [uf]B; the inertial position [sBiY: the inertial body rates [coB,]B – and the direction cosine matrix. Now, the second formulation of the six-DoF equations of motions over an elliptical Earth is complete.

Here are your choices. You can model your hypersonic vehicle relative to the inertial or the geographic frames over an elliptical rotating Earth and exercise several methods of initialization. I recommend the inertial approach because of its simplicity and I would combine it with the direction cosine matrix computations. If your customer requests a spherical Earth instead of its elliptical shape, you can accommodate him easily by setting the flattening parameter to zero and the semimajor axis of the ellipsoid to the mean Earth radius.

As you build your simulation, let the CADAC GHAME6 code be your guide. There you will also find other interesting modeling applications for aerodynamics and autopilots, which are the subject of the following sections.