- Email: [email protected]
Control system design of flexible-body launch vehicles
Control Engineering Practice 7 (1999) 1163}1175
Control system design of #exible-body launch vehicles Hidehiko Mori* Sumitomo Heavy Industries, Ltd., Research and Development Centre, 63-30 Yuhigaoka, Hiratsuka-shi, 254 Japan Received 4 November 1998; accepted 30 April 1999
Abstract Based on a series of H-I launches in Japan, a review of control system designs of expendable launch vehicles is presented in terms of their #exibility. The equations of motion are "rst formulated for lateral and longitudinal dynamics of a #exible vehicle, after which the aerodynamic loads caused by in-#ight atmospheric wind are described in relation to trajectory planning. Included are mode analyses for launch vehicles and time responses to the #ight environment. Finally, stability analyses for lateral and longitudinal vibrations of launch vehicles are presented with special emphasis on POGOs. 1999 Elsevier Science Ltd. All rights reserved. Keywords: Attitude control; Mode analysis; Oscillation; POGO; Stability analysis; Trajectory planning; Vehicle dynamics
1. Introduction The launching of expendable vehicles began a rich history following von Braun's successful #ight of the prototypical V-1 rocket some 60 years ago. Earlier, basic principles of control designs were established, and a computation procedure incorporating empirical methods became available for the practical work of trajectory planning and control system design. The procedure, however, is considered too detailed and is limited with regards to providing an overall picture. This situation led to the present paper which uses simpli"ed analytical methods to provide a perspective view of the treatment of #exible launch vehicles. As numerical examples, launches of Japan's H-I vehicle from 1986 to 1992 are considered. Due to signi"cant aims in reducing the structure masses of launch vehicles such that the payload mass ratio is raised, they become highly #exible. In-#ight winds cause bending moments that negatively a!ect vehicle strength, while bending motion coupled with sloshing of propellants produces stability problems. Much research in 1960s was directed at the control of launch vehicles for atmospheric #ights (e.g., Frosh & Valley, 1967; Johnston & Johnson, 1967). Using control designs based on conventional theories, Geisler (1970) went on to edit analytical treatments for a highly #exible launch
* Corresponding author. Fax: #82-463-21-8454. E-mail address: hdh}[email protected] (H. Mori)
vehicle intended for Saturn-V; while much more recently, new approaches using H theory have been applied by Morita, Kawaguti, Nakatani, Goto and Ohtsuka (1994) for the M-V launch vehicle and Mau!rey and Schoeller (1998) for the Ariane type. However, since the basic concept of control systems remains more or less shaped by conventional ideas, here their applicability is reviewed. During satellite launching, the environmental conditions are mainly determined by longitudinal vibrations; and in fact, resonant oscillations called POGO have appeared during a number of launchings. Early studies of this phenomenon were carried out by Mckenna (1964) and Rubin (1966), and subsequently led to the development of a practical computational program. Although Pain and Rubin (1974) applied the method numerically to the Delta vehicle, and a similar method was applied to the H-I, an over abundance of data was required to accurately compute all features of POGO phenomena. Here, such features are elucidated by application of a simpli"ed method using H-I #ight data.
2. Lateral and longitudinal dynamics of a 6exible vehicle The #exible characteristics of a launch vehicle can be expressed using a beam model representing the lateral and longitudinal motions, respectively. Fig. 1 shows notations used for the analyses. Since the treatment of both motions is similar, analogous descriptions are adopted in subsequent sections.
0967-0661/99/$ - see front matter 1999 Elsevier Science Ltd. All rights reserved. PII: S 0 9 6 7 - 0 6 6 1 ( 9 9 ) 0 0 0 8 7 - 8
1164
H. Mori / Control Engineering Practice 7 (1999) 1163}1175
satisfy the orthogonal relations
Fig. 1. Employed notations for analysis of lateral and longitudinal dynamics.
The moment M(x, t) on the cross section of a beam at position x is represented by
"
#
*
*z *x
(m!x) !m(m)
V
*z(m, t) dm, *t
(1)
where z(x, t) is the displacement of the cross section at x in the z-direction, EI(x) the bending sti!ness, q(x, t) the distributed force in the z-direction, and m(x) the mass of the cross section per unit length. Twice di!erentiating Eq. (1) gives the equation of motion for the lateral dynamics of the beam, i.e.,
* *z *z EI #m(x) "q(x, t). *x *x *t
(2)
Boundary conditions are
while initial conditions are
The homogeneous equation is de"ned as (3)
where the set of eigenvalues and the eigenfunctions, i.e., u , u , u ,2. Z (x), Z (x), Z (x),2,
q (x, t) dx, * I $ (t)" (x!x )q(x, t) dx, !% * M (gK #ug )" Z q(x, t) dx. (4) G G G G G Now let us consider that the distributed aerodynamic force is proportional to the angle of attack a and that the z-component of the thrust is proportional to the gimbal de#ection angle b. Accordingly, q(x, t) is q(x, t)"q SC (x)a#¹ b*(x!x ). (5) 8? A @ Substituting Eq. (5) into the "rst and the second of expressions of Eq. (4) subsequently gives zK "Z$ a#Z$ b, ? @
$ "U$ a#U$ b, ? @ where
(6)
* C (x) dx/M, 8? Z$ "¹ /M, @ A * U$ "q S (x!x )C (x) dx/I, ? !% 8? U$ "l ¹ /I. @ @ A Eq. (6), respectively, represents the translation and rotation of the vehicle's rigid body motion, while the third expression in Eq. (4) gives the equation of motion for the ith bending mode, i.e.,
z(x, 0)"0, z (¸, 0)"0.
*
Z$ "q S ?
z (0, t)"0, z (¸, t)"0,
dZ d EI !um(x)Z"0, dx dx
MzK (t)"
z (0, t)"0, z (¸, t)"0
(m!x)q(m, t) dm
V *
z(x, t)"Z (x)z (t)#Z (x) (t)# Z (x)g (t), G G G and substituting it into Eq. (2) yields
2.1. Lateral motion
M(x, t)"EI
* Z m(x)Z dx"M d , G H G GH * d dZ H dx"M ud (i"1,2,32), EI Z G G GH G dx dx with the eigenfunctions for u "0 being Z (x)"1, Z (x)"x!x . !% Assuming the solution of Eq. (2) to be
gK #ug " q S G G G
*
Z (x) ¹ Z (x )b C G a# A G @ b. 8? M M G G
(7)
H. Mori / Control Engineering Practice 7 (1999) 1163}1175
Substituting Eq. (5) into Eq. (1) gives the bending moment as M(x, t)"M (x, t) a(t)#M (x) b(t)# M (x)gK (t), ? @ EG G G where
M (x)"q S ?
* V
!U$
?
(m!x)C (m) dm!Z$ X? ?
*
V
!U$ @
"!
*
V
V
m(m)(m!x) dm
V
*
m(m)(m!x ) dmM K G(x) !% E
m(m)(m!x)Z (m) dm. G
* * *u *u " f (m, t) dm# !m(m) dm, (9) *x *t V V where u(x, t) is the axial de#ection of the cross section at x, AE(x) the elasticity of the beam, and f (x, t) the distributed axial force on the cross section at x. Di!erentiating Eq. (9) by x gives the equation of motion for the longitudinal dynamics of the beam, i.e.,
(10)
Boundary conditions are
u(x, 0)"0, u (x, 0)"0. The homogeneous equation is de"ned as
(11)
where the set of eigenvalues and the eigenfunctions, i.e., u , u , u ,2. ; (x), ; (x), ; (x),2, satisfy the orthogonal relations
* ; m(x); dx"M d , G H G GH * d d; H dx"!M ud , ; AE G dx G G GH dx
uK (t)"a !g cos h, V where
(14)
a "(¹!X )/M. V ? Eq. (14) represents the longitudinal acceleration of the vehicle's rigid body motion, while the second expression of Eq. (12) gives the equation of motion for the ith longitudinal #exible mode, i.e., ; (¸) ; (x ) (15) uK #uu " G @ ¹! G X . ? G G G M M G G Substituting Eq. (13) into (9) gives the axial force as
while intial conditions are
* m(m) dm V * ! m(m); (m) dm uK (t) for x)x H H @ H V * "!X !a m(m) dm ? V V * ! m(m); (m) dm uK (t) for x*x . (16) H H @ H V
P(x, t)"a M! V
u (0, t)"0, u (¸, t)"0,
d d; AE #um(x);"0, dx dx
f (x, t) dx,
f (x, t)"¹*(x!x )!X *(x!¸)#m(x)g cos h. (13) @ ? Substituting Eq. (13) into the "rst expression of Eq. (12) gives
The axial force on the cross section of the beam at x is represented by
* ; (x) f (x, t) dx. (12) G Now, let us consider the axial thrusting force ¹ at the gimbal point, the aerodynamic drag X at the top of the ? vehicle, and the gravity force distributed along the vehicle. Accordingly, f (x, t) is
V
*u *u * AE !m(x) "!f (t). *x *t *x
*
M (uK #uu )" G G G G
m(m)(m!x) dm
u(x, t)"; (x)u (t)# ; (x)u (t), G G G and substituting it into Eq. (10) yields
2.2. Longitudinal motion
P(x, t)"AE
with the eigenfunctions for u "0 being ; (x)"1. Assuming the solution of Eq. (10) to be
MuK (t)"
m(m)(m!x)(m!x ) dm, !%
M (x)"(x !x)¹ !Z$ @ @ A @ *
*
(8)
1165
3. Trajectory planning One of the most important tasks of the #ight system design of a launch vehicle is to reduce the aerodynamic lateral loads during atmospheric #ight. Eq. (16) shows that the bending moment is proportional to the angle of attack a, the gimbal de#ection angle b, and the secondorder time derivatives of general coordinates g s. G
1166
H. Mori / Control Engineering Practice 7 (1999) 1163}1175
Geisler (1970) reported that the e!ect of vibration modes is dominant at the upper part of a huge launch vehicle such as Saturn-V; and in fact, even for an ordinary vehicle size, such as that for Japanese rockets, lateral vibration has a substantial e!ect on the structure of a launched satellite. However, if the system is stable, the vibration modes do not greatly a!ect the launch vehicle structure such that the #ight system design is executed mainly to suppress the magnitudes of the attack angle a and gimbal de#ection angle b, respectively. The e!ect of lateral vibration is discussed in Section 4. In-#ight wind is the main factor causing large values of a and b, and an e!ective method to suppress them is to set the vehicle attitude command such that it is aligned with the relative velocity so as to anticipate the e!ect of wind. Since the discussion on pitching dynamics is applicable to yawing with minor changes, only pitching motion is considered using the #ight notations shown in Fig. 2, where the direction of the vehicle velocity relative to earth rotation is indicated by #ight path angle c. In order to realize a zero angle of attack, the direction of vehicle axis h must satisfy condition h"c!a , (17) U where a is the contribution of in-#ight wind to the angle U of attack. Under this condition, the gimbal de#ection angle b is also maintained at zero, with the resultant #ight mode being referred to as `zero-lift turna. Fig. 3 shows the average monthly wind speeds and directions for Kagoshima prefecture where Japan's launch sites are located. Note that the wind direction at altitudes higher than 3 km is almost exactly west. Rocket launchings in Japan are allowed only during two seasons: the windless summer and windy winter. Since the angle of attack caused by wind a is nearly zero U in summer, but reaches nearly 103 in winter, the respect-
Fig. 2. Notations for pitching dynamics.
Fig. 3. Average monthly wind speeds and directions for Kagoshima prefecture.
ive attitude angle programs required for zero-lift turn are naturally quite di!erent. During launch, the attitude command is given in a series of step functions of attitude rate. An example of the actual pitch command rate hQ used for winter launching of the H-I vehicle is shown in A Fig. 4. The H-I lifts o! vertically and 8}12 s after lift-o! takes a kick turn at a constant pitch rate. Then the zero-lift turn is executed for about 90 s until aerodynamic pressure su$ciently decreases. In the task of #ight analysis, the attitude commands are determined via numerical calculations to minimize fuel consumption from lift-o! to orbit injection under given constraints. Fig. 4 also shows #ight pro"les for c and a , obtained as the result of the U #ight simulation with the attitude rate command hQ . It is A apparent that the zero-lift turn rule is observed well. Similar results were obtained for the yaw rate command by replacing #ight path angle c with #ight direction angle t. If the in-#ight wind coincides with the average wind used to calculate the attitude rate programs, then the aerodynamic loads on the vehicle are maintained at close to zero. However, this is not the case for all #ight altitudes, and actual in-#ight wind has turbulence with high spatial frequencies which are smoothed out in the average wind. It is therefore necessary to con"rm that the loads remain within allowed limits during all possible wind conditions. The vehicle's response to an unexpected wind condition varies depending on the design of the attitude control system. Such a condition can be described using
H. Mori / Control Engineering Practice 7 (1999) 1163}1175
1167
Fig. 5. Signal #ow diagram for the pitching control system.
Fig. 4. Pitch command, #ight path angle, and angle of attack pro"les for a winter launching of the H-I.
perturbational equations derived around the nominal equations of motion. The zero-lift turn trajectory is taken as nominal motion, u, *c, and *a are the perturbations U of h, c, and a from nominal trajectory, and a and b are U themselves perturbed quantities. Under these conditions, Eq. (6) represents the perturbational equations for such pitching motion. The lateral displacement z is related to *c by z"<*c,
(18)
with Eq. (17) yielding
"*c!*a . U
(19)
The sensors used for the pitch attitude control system are a position gyro, rate gyro, and accelerometer. Fig. 5 shows the signal #ow diagram of the control system for pitching motion, where a and a s are the feedback gains of the position and rate gyros, respectively. Broken lines show the feedback paths of the accelerometer and g is the feedback gain, while other lines show the pitching dynamics represented by Eqs. (6), (18) and (19). Note that the line of Z$ , which represents coupling ( between the lateral and longitudinal motions, has been added to Eq. (6). Accelerometer feedback is not typically used, although severe #ight conditions make such feedback necessary, i.e., it can give the control system the characteristics of drift minimum and load relief for the wind; concepts developed by Hoelker (1961).
Under the signal #ow system (Fig. 5), the transfer functions from *a to zK /< are given by U zK /< A s#A s#B " . (20) *a B s#B s#B s#B U For a step function of *a , the lateral acceleration zK /< U reaches the input level unless B "0; a case in which the stationary value of the transfer function becomes A /B which is 1. The relation for B "0 is (a !g Z$ )(Z$ U$ !Z$ U$ )!Z$ U$ "0, (21) ( ? @ @ ? ( ? being called drift minimum condition, which is only realized when accelerometer feedback is present. Fig. 6 shows the typical gain schedule for an H-II vehicle with accelerometer feedback. Transfer functions from *a to a and b are both U represented as A s#A s#A . (22) B s#B s#B s#B Thus, for a step function of *a , both a and b normally U have stationary values. If the position feedback gain a "0, A s for both parameters become zero; and in such case, their stationary values are reduced to zero. The relation a "0 is called the load relief condition. Even if a is not exactly zero, the e!ect of load relief can be expected to occur. The application of the drift minimum condition leads to a smaller value of a , which is useful for the purpose of drift minimum and load relief as well. Eq. (22) indicates that a and b have stationary values for a step function of *a , although they become zero for U a step function of *a . In other words, a severe wind U for launch vehicles has large step functions of *a . U A synthetic wind is used to con"rm that the #ight trajectory does not incur aerodynamic loads beyond the vehicle structure limit. Fig. 7 shows an example of synthetic wind which has the peak at 9 km in altitude. The wind in the upper atmosphere is measured using balloon devices, and the measured data is considered independent along 500-m altitude intervals. At each interval, a 95% level wind is de"ned as one whose magnitude is on the 2.45p probability ellipse in the #ight
1168
H. Mori / Control Engineering Practice 7 (1999) 1163}1175
Fig. 6. Gain schedule for the H-II vehicle.
Fig. 8. a and a for a winter launching of the H-I. U U
and its time derivative a and a , where it is apparent U U that a for the synthetic wind shows a large step function U at an altitude of 9-km in comparison with a of the U average wind. Thus the synthetic wind causes large values of a and b. The angle of attack a and the gimbal de#ection angle b generated by the synthetic wind are used to calculate the bending moment using Eq. (8), with Fig. 9 showing the resultant bending moment and axial force acting on the H-I vehicle. Several synthetic winds are formed with peaks distributed from 2 to 15 km, being used to con"rm the loads imposed on the vehicle are below allowed limits. The total loads are described in terms of the equivalent axial force, i.e.,
Fig. 7. Synthetic wind with peak velocity at an altitude of 9 km.
direction. There are two 95% levels at each altitude: high and low. The 95% pro"les shown in Fig. 7 are obtained by connecting measured 95% level winds along successive 500-m altitude intervals. The wind shear for each measurement is taken as the di!erentiated value of wind with respect to altitude, with 95% wind shear being similarly de"ned. The synthetic wind that peaks at 9-km is formed by connecting 95% wind shears along successive 500-m altitude intervals such that the peak wind occurs at that altitude. Fig. 8 shows the angle of attack due to the wind
4 P (x, t)"P(x, t)$ M(x, t), CO D
(23)
where P(x, t) is the axial force, M(x, t) the bending moment, and D the diameter of the beam. A negative sign indicates compression.
4. Mode analysis and time response Three-dimensional (3D) computation of the mode frequencies and mode shapes of a launch vehicle is carried out using a model composed of a number of series of beam elements, each of which has 5 degrees of freedom, i.e., x, y, z, h and h . Torsion is not considered in the W X
H. Mori / Control Engineering Practice 7 (1999) 1163}1175
1169
Fig. 9. H-I bending moment and axial force for a 9-km peak synthetic wind.
model. The axial displacement x is not coupled with lateral displacements or rotations, while lateral displacement y is only coupled with h , and z only with h . X W Accordingly an element has three non-coupled members whose degrees of freedom are in x, y!h , and z!h , X W respectively. Members are represented by an equivalent mass/inertia-spring component, and are connected to similar members of adjacent elements. Under this approach, the vehicle model is basically represented by three series of members in axial, pitch, and yaw axes. These three series are connected to each other at the engine sections, guidance section, and satellite model, while also being coupled as a whole system. About 500 degrees of freedom are considered "rst. Application of the constrained normal mode method reduces the degrees of freedom of the model to about 100, which are subsequently used in the computation of mode frequencies and shapes. Fig. 10 shows the distribution of the mode frequencies of the H-I vehicle through a #ight, where lateral modes appear in the lower frequency range and axial in the higher. Each mode shape is given in 3D coordinates of pitch, yaw, and axial. The dimension in which a mode mainly belongs can be identi"ed by the dominant coordinate, although more than a few modes do not have a speci"c dominant coordinate. The "rst axial mode appears in the 19th order of frequency, after the 6th pitch and 6th yaw modes, respectively, appear in the 17th and 18th orders, while the second and third axial modes appear in the 24th and seemingly in the 44th, respectively. It should be noted, however, that (i) this model does not show good accuracy in the higher orders and (ii) changes in the trend of modes greatly depends on the time after lift-o!. Since much work is required to construct models for di!erent #ight times, this makes it di$cult to predict an accurate
Fig. 10. Distribution of H-I mode frequencies.
trend of modes for the high-frequency range; a problem discussed later in relation with vehicle longitudinal stability. The model composed of beam elements is also used to compute the response of the satellite to shocks and oscillations generated by the launch vehicle. The greatest load imposed on the satellite occurs during launching, being the load which prescribes basic requirements for the satellite structure. Naturally then, it is important to precisely predict the response of the satellite to such shocks and oscillations. Eq. (16) indicates that the axial force P(x, t) is the sum of the weight of the upper part of the structure * produced by static acceleration of the vehicle * and inertia forces of the distributed mass which are related to the general coordinates due to modal oscillations. At the time of main engine cut-o! (MECO), static acceleration ultimately increases. Moreover, during H-I launches, longitudinal oscillations called MECO POGO occurred due to resonance of the propellant feed system and engine chamber. Hence, the load at MECO prescribed the main
1170
H. Mori / Control Engineering Practice 7 (1999) 1163}1175
design requirement. Fig. 11 shows the axial acceleration at the "rst engine section of the H-I vehicle. Another resonant oscillation (20}33 Hz) called preMECO POGO also appeared during H-I launches. The satellite launched by the H-I vehicle was required to have its own 1st longitudinal mode frequency higher than 35 Hz and the 1st lateral mode frequency higher than 15 Hz. The pre-MECO pogo, whose frequency was close enough to the satellite's lowest longitudinal frequency, imposed a serious burden on the structure design. The dynamics of POGOs are treated in Section 6. If the satellite is considered to have a simple beam structure } normally being more complex } then the oscillation produces an inner force to the satellite represented by the last terms of Eqs. (8) and (16). Thus subsequent time response computations using the employed element model were necessary to examine loads on each part of the satellite. The results were re#ected in vibration tests of the structure models and actual satellites. Extreme shocks act on the vehicle at the time of the lift-o! and at max q a, i.e., the maximum value of the product of the dynamic pressure and angle of attack. For the H-I, max q a occurred at an altitude of 9 km. These shocks had lateral components and stimulated the bending modes shown in Fig. 10. Those transmitted to the satellite had a high spectrum in the lower frequency
range ((15 Hz), being close to the satellite's lowest lateral mode frequency. This is another reason making it essential to carry out time response investigations using the developed model.
5. Stabilization of lateral motion The basic equation of motion for lateral vibration of a vehicle in each bending mode is represented by Eq. (7), where through control feedback, the motions of bending modes are coupled with each other and with the rigid body. Moreover, the control system considers the sloshing motions of propellants and engine-actuator dynamics, as these characteristics must be considered together for the stability analysis of the vehicle's lateral motion. The total equations reached almost 30th-order, and included at least "ve bending modes and two sloshing modes. The rigorous equations of motion for a highly #exible vehicle can be derived using Lagrange's equations, and an example of such derivation can be seen in the work of Luckens, Schmit and Broucek (1961). The derived equations of motion can be written as (Ms#Bs#K) x"k b#k w, (24) @ U where M, B and K are matrices, x is a vector whose component takes one degree of freedom, b is the gimbal de#ection angle, w is the perturbation of the wind and k and k are scalar coe$cients. Fig. 12 shows a typical @ U feedback system used to represent the lateral motion of a vehicle, in this case that of the H-I. Feedback loops can be represented mathematically as b"G(s)x,
(25)
where G(s) is a column vector whose components are rational functions of s. Stability is examined in time-slice pattern along the trajectory using either Nyquist diagrams or Bode plots. Feedback gains are determined as described in Section 3, and they should not be altered for stabilizing purposes. That is, since the natural frequency of rigid body motion is low enough compared to those of other modes, it is not di$cult to stabilize the system by shaping frequency characteristics; and a combination of lead and lag compensations is generally used.
Fig. 11. Longitudinal oscillations measured at the engine section during H-I launchings.
Fig. 12. A block diagram of the H-I's guidance and control system.
H. Mori / Control Engineering Practice 7 (1999) 1163}1175
As the compensating "lters in the feedback loops are associated with the inertial guidance unit (Fig. 12), the "lter transfer functions must be changed into the digital algorithms via bilinear transformations, i.e., 2 1!z\ D(z)"D(s)" , f (z)" , QD X ¹ 1#z\
(26)
where z"e2Q and ¹ is the sampling period for digital computation. The digital algorithms obtained via Eq. (26) are accurate for a frequency range below ¹. Given that the sampling period for the H-I was 25 ms, the digital "lters were accurate for a frequency range below 5 Hz, although they might also be e!ective for stabilization at a higher range. For example, Goto and Suzuki (1981) designed the H-Is digital "lters and con"rmed their e!ectiveness below 25 Hz. After passing through the atmosphere, the guidance loop becomes active, and the stability of the system during this phase of #ight, which includes both guidance and control loop's, must be examined since the computation period for guidance is longer than that for control. Fig. 13 shows a block diagram of the sampled-data system incorporating two di!erent sampling periods, i.e., ¹ and ¹ are, respectively, the control and guidance Q E loop sampling periods * z "e2Q and z "e2E * while Q E D and D are the compensating "lters and H and H the Q E zero-order holds for each sampling. When the guidance loop sampling period is an integer multiple of the control loop's, the characteristic equation of the sampled-data system is represented by
H FD (z )[H ]H2Q H2Q Q Q E Q(z )"0. (27) E 1#D (z )[H F]H2Q E Q System stability can be examined using the Nyquist diagram of the open-loop transfer function of Eq. (27), with this stability analysis allowing the gains of the guidance to be con"rmed. In Japan, most of solid motor rockets use tail "ns for aerodynamic stability. Since strong control forces are not present, such rockets can #y without an attitude control system, and do not necessarily fail if they happen to be unstable during #ight. In fact, several #ights were suc1#
Fig. 13. Sampled-data system incorporating samplers for two di!erent periods.
1171
cessful even though the control system was found to be unstable. In comparison with these conventional vehicles, the solid-motor M-V vehicle, newly developed by the Institute of Space and Astronautical Science (ISAS), has no tail "ns and is inherently aerodynamically unstable. It is obliged to be equipped with a control system for stability, which works in the presence of uncertain system parameters such as the bending rigidity, the aerodynamic coe$cients and the actuator bandwidth. Morita et al. (1994) applied the H control theory in the construction of its feedback system; and under a framework of the mixed sensitivity problem, they achieved better control performance for a relatively low-frequency region and preferable robust stability in a higher one. A similar application of the theory is investigated by Mau!rey and Schoeller (1998) for a launch vehicle of the Arian type with liquid propellants. 6. Stability analysis of londitudinal vibration As described in Section 4, longitudinal vibration is the main factor determining satellite environmental conditions during launch. Three types of POGO resonant oscillations have appeared in N-II and H-I rocket #ights. MECO POGOs appear just before main engine cuto!, being caused by #uctuations in LOX tank bottom pressure, and two kinds of pre-MECO POGOs sometimes appearing in the midst of #ight. Although the latter POGOs have not been well analyzed, possible causes were investigated by Mori (1991) using a simpli"ed analysis method in conjunction with #ight data. Here an outline of the analysis is described. Eq. (15) represents the basic equation of longitudinal motion of a vehicle, where a single mode is excited by resonance with thrust variation. By adding structure damping, the equation of motion becomes ; (x )C A uK #2fu u #uu " G @ 2 RFp , (28) G G G G G A M G where C is the thrusting coe$cient, A the nozzle throat 2 RF area of the engine, and p the chamber pressure. During A resonant oscillation, p is a!ected by structure vibrations. A Fig. 14 shows signal #ow diagrams for such dynamics,
Fig. 14. Signal #ow diagram representing Eq. (28).
1172
H. Mori / Control Engineering Practice 7 (1999) 1163}1175
where [ p /uK ] is the transfer function from uK to p , while A G G A K is the gain from p to uK . 2 A G The open-loop transfer function of this system is given by
u u p A . O¸(s)"2f G# G !K G s 2 uK s G
(29)
Letting s"ju, Eq. (29) takes the form
O¸( ju)"!
u p G #K Re A 2 u uK G
u p !j 2f G#K Im A Gu 2 uK G
.
(30)
The Nyquist diagram of Eq. (30) is represented conceptually by Fig. 15. That is, when the structure is not excited by chamber pressure, i.e., when K 0, the Nyquist diagram is in2 dicated by the broken line. In this case, for u"u , G Re(O¸)"!1, Im(O¸)"!2f , G and the systm is stable. On the other hand, when K 0, 2 there is a possibility that the last term of Eq. (30) becomes large such that, as shown by the solid line, the diagram rises up over (!1, 0) depending on the value of Im[ p /uK ]. It indicates the system is unstable. The transA G fer function [ p /uK ], whose notations are indicated in Fig. A G 16 and applied to the fuel and the oxygen feed system, respectively, is determined as follows. A pressure perturbation at the bottom of the tank is given by p "oh; (x)su (s), R G G
(31)
where o is the density of propellant and h the height of propellant in the tank. The characteristics of #ow and pressure transmitted through the pipe line is represented by ¸ cos uq ) p # J sin uq ) s Q "p , Q uq Q R
(32)
Fig. 16. Notations for pogo stability analysis.
where p and Q are, respectively, perturbations of presQ Q sure and #ow at the outlet of the feedline due to oscillation; and ¸ is the inertance of propellant in the feed line J and q the time required for the sound wave to pass through it. The #ow perturbation at the outlet of the accumulator, Q is represented by Q 1 Q !Q " sp , (33) Q Q ¸ s#R s#K Q ? ? ? where R is the resistance, ¸ the inertia, and K the ? ? ? volume elasticity of the accumulator, respectively. The #ow to the discharge line is represented by Q "Q !C sp #A su(x , s), (34) B Q @ Q J N where C is the cavitation compliance and A the cross@ J sectional area of the feedline. The perturbation in chamber pressure is expressed as p "ap !(¸ #¸ )sQ #(R #R )Q , (35) A Q N B B N B B where a is the ampli"cation coe$cient of the pressure variation from the pump inlet to outlet, and ¸ , ¸ , and N B R and R are the inertance and resistance of the pump N B and discharge line, respectively. The propellant #ow raises the chamber pressure in accordance with
Fig. 15. Nyquist diagrams of the open loop transfer functions of Eq. (30).
cHo (1#q s)p " Q, (36) A A A B RF where q is the time constant of the combustion and A cH the characteristic exhaust velocity. The transfer function [ p /uK ] is obtained eliminating "ve variablesA G p , p , Q , Q and Q * from Eqs. (32)}(36). The bubble R Q Q Q B
H. Mori / Control Engineering Practice 7 (1999) 1163}1175
1173
angular velocity u , i.e., @ u "1/(C ¸ (37) @ @ R is considered an important parameter related to POGO phenomena. This is also expressed in the term of bubble frequency by f "u /(2n). If we consider a simpli"ed case @ @ in which no accumulator is installed, and further assume that R"R #R "0 and q "0, then Eq. (30) becomes N B A O¸( ju)
"!
u cd G #K Ka 2 N u *
u c¸ #a¸ tan uq/(uq) J !j 2f G#K Ka NB du , Gu 2 N * (38) where *"(cK)#+c¸ #a¸ tan uq/(uq),u, N NB R u tan uq , c(u)"1! u uq @ o sin uq d(u)" h; (x )# l; (x ) . G R G N cos uq uq
The unstable condition shown in Fig. 15 occurs when c(u ) is nearly zero, while the condition c(u )"0 is satisG G "ed when u equals `the resonant bubble angular velo@ citya - , i.e., @ tan u q @ " G . (39) u uq G G - is a good parameter for indicating which bubble @ frequency produces vehicle POGO, although some simpli"cation is applied. This is also expressed in the term of resonant bubble frequency fM by @
tan 2nf q G . (40) 2nf q G Fig. 17 shows the #ight time periods and frequencies of typical POGO phenomena appearing during three H-I launches (*, 䉭, and ;). Two pre-MECO POGOs, whose frequencies are indicated by f and f , are named PMP-1 and PMP-2, G G respectively, while MECO POGO related with the "rst longitudinal frequency f is abbreviated to MP. The resonant bubble frequency fM corresponding to each @ POGO frequency f was determined using Eq. (40), with G Table 1 summarizing values of fM corresponding to the @ central frequency of each POGO. The value of q for the LOX and fuel system is 0.26;10\and 1.72;10\ s, respectively. Since Eq. (40) stands for 2nf q n in the fuel system, the fM for PMP-2 G @ can be markedly smaller than f . While POGO oscillaG fM "f @ G
Fig. 17. Pogo frequency/bubble frequencies of the H-I.
tions during PMP-1 and MP can occur only when 2nf q 0 for the LOX system. G Since the fuel system is not equipped with an accumulator, Eq. (40) is expected to have a good accuracy and f likely coincides with fM . In the accumulator-equipped @ @ LOX system, however, fM is not expected to have a close @ value to f . Use of Eq. (37) to directly determine bubble @ frequencies is problematic because it is di$cult to determine cavitation compliance accurately; hence, bubble frequencies were estimated such that stability conditions match the #ight status. To compute Eq. (25) using Eqs. (26)}(31), the oscillated (POGO) mode shapes and their generalized masses must be known. Using measured #ight data for each abovementioned H-I launch, Fig. 18 shows the longitudinal mode shapes corresponding to POGOs, i.e., PMP-1-1, PMP-1-2, PMP-2, and MP. PMP-1 are divided into two sections of time. Longitudinal accelerations were measured at indicated positions (AF,GS,2GP, and IES). Also shown is the mass distribution at lift-o!. The mode shapes of POGOs in MP and PMP-2 were calculated using the model described in Section 3. The mode shapes for those in PMP-I, however, were not calculated as these oscillations did not critically impact the satellite; instead they were estimated based on measured data. With mode shapes and the mass distribution in hand (Fig. 18), time-varying values of the generalized masses can be determined for each POGO, with results being shown in Fig. 19. Next, bubble frequencies are estimated in order to explain the pogo phenomena. Table 2 shows bubble frequencies corresponding to POGO frequencies of Table 1.
1174
H. Mori / Control Engineering Practice 7 (1999) 1163}1175
Fig. 19. Generalized masses corresponding to longitudinal mode shapes in Fig. 18.
Table 1 Resonant bubble frequency corresponding to respective POGO frequency Fig. 18. Longitudinal mode shapes of POGOs during H-I launches. PMP-1 PMP-2 MP
Resultant values were determined by trial and error. The broken lines in Fig. 17 indicate time-varying values of bubble frequencies for the LOX and fuel systems. Finally, the estimated bubble frequencies help in determining respective Nyquist diagrams at various #ight times. Fig. 20 shows the result, where the stability characteristics represent the time-varying POGO phenomena of Fig. 11 well. No POGO-induced longitudinal vibrations have occurred during launch of the H-II, which unlike the H-I has its lox tank located above the LH2 tank. While this change in con"guration is thought to prevent MECO POGOs, the possibility of Pre-MECOs still exists, although none have been observed thus far. The design of the POGO suppressor (i.e., accumulator) installed in the H-II may be e!ective.
7. Concluding remarks The results obtained using the employed analytical model indicate that the design of a control system for a #exible-body launch vehicle is highly interrelated with trajectory planning and structure analysis consideration.
f (Hz) G
LOX fM (Hz) @
Fuel fM (Hz) @
22.0 32.0 17.5
22.5 * 17.7
* 9.85 *
Table 2 Bubble frequency corresponding to respective POGO frequency
PMP-1 PMP-2 MP
f (Hz) G
LOX f (Hz) @
Fuel f (Hz) @
22.0 32.0 17.5
33.5 * 23.8
* 9.9 *
In addition, through the use of analytical methods, critical #ight characteristics are elucidated, i.e., aerodynamic loads acting on the launch vehicle, lateral and longitudinal stability, and launching environments of the housed satellite. Based on a stability analysis of POGOs, it is shown that their occurrence is estimated by examining resonant bubble frequencies, around which bubble frequencies of the feed systems are determined such that time-varying POGO phenomena of the H-I can be represented by Nyquist diagrams.
H. Mori / Control Engineering Practice 7 (1999) 1163}1175
1175
References
Fig. 20. Nyquist diagrams for indicated POGO periods during H-I launches.
Frosh, J. A., & Valley, D. P. (1967). Saturn AS501/S-IC #ight control system design. Journal of Spacecraft, 4(8), 1003}1009. Geisler, E.D. (1970). Wind ewects on launch vehicle. AGUARDograpph No. 115. Goto, S., Suzuki, H. (1981). A design of digital yight control system of a launch vehicle. NASDA TR-10, National Space Development Agency of Japan. Hoelker, R.F. (1961). Theory of artixcial stabilization of missiles and space vehicle with exposition of four control principles. NASA TND-555. Johnston, D.E., Johnson, W.A., (1967). Feasibility of conventional control techniques for large highly coupled elastic boost vehicles. STI Technical Report, No.67-592. Lukens, D.R., Schmit, A.F., Broucek, G.T. (1961). Approximate transfer functions for yexible-booster-and-autopilot analysis. WADD-TR-6193. Mau!rey, S., Schoeller, M. (1998). Non-stationary H control law for launcher with bending modes. Preprint of 14th IFAC symposium on automatic control in aerospace. (pp. 438}443). McKenna, K.J. et al. (1964). A method for studying the coupled engine airframe longitudinal instability of liquid rocket systems. AIAA paper 64}81. Mori, H. (1991). POGO analysis based on N-II/H-I vehicle #ight data. Journal of the Japan Society for Aeronautical and Space Sciences, 39, 71}80. Morita, Y., Kawaguti, J., Nakatani, I., Goto, S., Ohtsuka, H. (1994). Attitude control algorithm of the M-V launching rocket, Proceeding of 4th workshop on astrodynamics and yight mechanics. (pp. 169}174). Kanagawa, Japan. Pain, J.G., Rubin, S. (1974). Pogo suppression on the Delta vehicle. Report SAMSO-TR-74-180. Rubin, S. (1966). Longitudinal instability of liquid rockets due to propulsion feedback (POGO). Journal of Spacecraft, 3, 1185}1195.
Control system design of #exible-body launch vehicles Hidehiko Mori* Sumitomo Heavy Industries, Ltd., Research and Development Centre, 63-30 Yuhigaoka, Hiratsuka-shi, 254 Japan Received 4 November 1998; accepted 30 April 1999
Abstract Based on a series of H-I launches in Japan, a review of control system designs of expendable launch vehicles is presented in terms of their #exibility. The equations of motion are "rst formulated for lateral and longitudinal dynamics of a #exible vehicle, after which the aerodynamic loads caused by in-#ight atmospheric wind are described in relation to trajectory planning. Included are mode analyses for launch vehicles and time responses to the #ight environment. Finally, stability analyses for lateral and longitudinal vibrations of launch vehicles are presented with special emphasis on POGOs. 1999 Elsevier Science Ltd. All rights reserved. Keywords: Attitude control; Mode analysis; Oscillation; POGO; Stability analysis; Trajectory planning; Vehicle dynamics
1. Introduction The launching of expendable vehicles began a rich history following von Braun's successful #ight of the prototypical V-1 rocket some 60 years ago. Earlier, basic principles of control designs were established, and a computation procedure incorporating empirical methods became available for the practical work of trajectory planning and control system design. The procedure, however, is considered too detailed and is limited with regards to providing an overall picture. This situation led to the present paper which uses simpli"ed analytical methods to provide a perspective view of the treatment of #exible launch vehicles. As numerical examples, launches of Japan's H-I vehicle from 1986 to 1992 are considered. Due to signi"cant aims in reducing the structure masses of launch vehicles such that the payload mass ratio is raised, they become highly #exible. In-#ight winds cause bending moments that negatively a!ect vehicle strength, while bending motion coupled with sloshing of propellants produces stability problems. Much research in 1960s was directed at the control of launch vehicles for atmospheric #ights (e.g., Frosh & Valley, 1967; Johnston & Johnson, 1967). Using control designs based on conventional theories, Geisler (1970) went on to edit analytical treatments for a highly #exible launch
* Corresponding author. Fax: #82-463-21-8454. E-mail address: hdh}[email protected] (H. Mori)
vehicle intended for Saturn-V; while much more recently, new approaches using H theory have been applied by Morita, Kawaguti, Nakatani, Goto and Ohtsuka (1994) for the M-V launch vehicle and Mau!rey and Schoeller (1998) for the Ariane type. However, since the basic concept of control systems remains more or less shaped by conventional ideas, here their applicability is reviewed. During satellite launching, the environmental conditions are mainly determined by longitudinal vibrations; and in fact, resonant oscillations called POGO have appeared during a number of launchings. Early studies of this phenomenon were carried out by Mckenna (1964) and Rubin (1966), and subsequently led to the development of a practical computational program. Although Pain and Rubin (1974) applied the method numerically to the Delta vehicle, and a similar method was applied to the H-I, an over abundance of data was required to accurately compute all features of POGO phenomena. Here, such features are elucidated by application of a simpli"ed method using H-I #ight data.
2. Lateral and longitudinal dynamics of a 6exible vehicle The #exible characteristics of a launch vehicle can be expressed using a beam model representing the lateral and longitudinal motions, respectively. Fig. 1 shows notations used for the analyses. Since the treatment of both motions is similar, analogous descriptions are adopted in subsequent sections.
0967-0661/99/$ - see front matter 1999 Elsevier Science Ltd. All rights reserved. PII: S 0 9 6 7 - 0 6 6 1 ( 9 9 ) 0 0 0 8 7 - 8
1164
H. Mori / Control Engineering Practice 7 (1999) 1163}1175
satisfy the orthogonal relations
Fig. 1. Employed notations for analysis of lateral and longitudinal dynamics.
The moment M(x, t) on the cross section of a beam at position x is represented by
"
#
*
*z *x
(m!x) !m(m)
V
*z(m, t) dm, *t
(1)
where z(x, t) is the displacement of the cross section at x in the z-direction, EI(x) the bending sti!ness, q(x, t) the distributed force in the z-direction, and m(x) the mass of the cross section per unit length. Twice di!erentiating Eq. (1) gives the equation of motion for the lateral dynamics of the beam, i.e.,
* *z *z EI #m(x) "q(x, t). *x *x *t
(2)
Boundary conditions are
while initial conditions are
The homogeneous equation is de"ned as (3)
where the set of eigenvalues and the eigenfunctions, i.e., u , u , u ,2. Z (x), Z (x), Z (x),2,
q (x, t) dx, * I $ (t)" (x!x )q(x, t) dx, !% * M (gK #ug )" Z q(x, t) dx. (4) G G G G G Now let us consider that the distributed aerodynamic force is proportional to the angle of attack a and that the z-component of the thrust is proportional to the gimbal de#ection angle b. Accordingly, q(x, t) is q(x, t)"q SC (x)a#¹ b*(x!x ). (5) 8? A @ Substituting Eq. (5) into the "rst and the second of expressions of Eq. (4) subsequently gives zK "Z$ a#Z$ b, ? @
$ "U$ a#U$ b, ? @ where
(6)
* C (x) dx/M, 8? Z$ "¹ /M, @ A * U$ "q S (x!x )C (x) dx/I, ? !% 8? U$ "l ¹ /I. @ @ A Eq. (6), respectively, represents the translation and rotation of the vehicle's rigid body motion, while the third expression in Eq. (4) gives the equation of motion for the ith bending mode, i.e.,
z(x, 0)"0, z (¸, 0)"0.
*
Z$ "q S ?
z (0, t)"0, z (¸, t)"0,
dZ d EI !um(x)Z"0, dx dx
MzK (t)"
z (0, t)"0, z (¸, t)"0
(m!x)q(m, t) dm
V *
z(x, t)"Z (x)z (t)#Z (x) (t)# Z (x)g (t), G G G and substituting it into Eq. (2) yields
2.1. Lateral motion
M(x, t)"EI
* Z m(x)Z dx"M d , G H G GH * d dZ H dx"M ud (i"1,2,32), EI Z G G GH G dx dx with the eigenfunctions for u "0 being Z (x)"1, Z (x)"x!x . !% Assuming the solution of Eq. (2) to be
gK #ug " q S G G G
*
Z (x) ¹ Z (x )b C G a# A G @ b. 8? M M G G
(7)
H. Mori / Control Engineering Practice 7 (1999) 1163}1175
Substituting Eq. (5) into Eq. (1) gives the bending moment as M(x, t)"M (x, t) a(t)#M (x) b(t)# M (x)gK (t), ? @ EG G G where
M (x)"q S ?
* V
!U$
?
(m!x)C (m) dm!Z$ X? ?
*
V
!U$ @
"!
*
V
V
m(m)(m!x) dm
V
*
m(m)(m!x ) dmM K G(x) !% E
m(m)(m!x)Z (m) dm. G
* * *u *u " f (m, t) dm# !m(m) dm, (9) *x *t V V where u(x, t) is the axial de#ection of the cross section at x, AE(x) the elasticity of the beam, and f (x, t) the distributed axial force on the cross section at x. Di!erentiating Eq. (9) by x gives the equation of motion for the longitudinal dynamics of the beam, i.e.,
(10)
Boundary conditions are
u(x, 0)"0, u (x, 0)"0. The homogeneous equation is de"ned as
(11)
where the set of eigenvalues and the eigenfunctions, i.e., u , u , u ,2. ; (x), ; (x), ; (x),2, satisfy the orthogonal relations
* ; m(x); dx"M d , G H G GH * d d; H dx"!M ud , ; AE G dx G G GH dx
uK (t)"a !g cos h, V where
(14)
a "(¹!X )/M. V ? Eq. (14) represents the longitudinal acceleration of the vehicle's rigid body motion, while the second expression of Eq. (12) gives the equation of motion for the ith longitudinal #exible mode, i.e., ; (¸) ; (x ) (15) uK #uu " G @ ¹! G X . ? G G G M M G G Substituting Eq. (13) into (9) gives the axial force as
while intial conditions are
* m(m) dm V * ! m(m); (m) dm uK (t) for x)x H H @ H V * "!X !a m(m) dm ? V V * ! m(m); (m) dm uK (t) for x*x . (16) H H @ H V
P(x, t)"a M! V
u (0, t)"0, u (¸, t)"0,
d d; AE #um(x);"0, dx dx
f (x, t) dx,
f (x, t)"¹*(x!x )!X *(x!¸)#m(x)g cos h. (13) @ ? Substituting Eq. (13) into the "rst expression of Eq. (12) gives
The axial force on the cross section of the beam at x is represented by
* ; (x) f (x, t) dx. (12) G Now, let us consider the axial thrusting force ¹ at the gimbal point, the aerodynamic drag X at the top of the ? vehicle, and the gravity force distributed along the vehicle. Accordingly, f (x, t) is
V
*u *u * AE !m(x) "!f (t). *x *t *x
*
M (uK #uu )" G G G G
m(m)(m!x) dm
u(x, t)"; (x)u (t)# ; (x)u (t), G G G and substituting it into Eq. (10) yields
2.2. Longitudinal motion
P(x, t)"AE
with the eigenfunctions for u "0 being ; (x)"1. Assuming the solution of Eq. (10) to be
MuK (t)"
m(m)(m!x)(m!x ) dm, !%
M (x)"(x !x)¹ !Z$ @ @ A @ *
*
(8)
1165
3. Trajectory planning One of the most important tasks of the #ight system design of a launch vehicle is to reduce the aerodynamic lateral loads during atmospheric #ight. Eq. (16) shows that the bending moment is proportional to the angle of attack a, the gimbal de#ection angle b, and the secondorder time derivatives of general coordinates g s. G
1166
H. Mori / Control Engineering Practice 7 (1999) 1163}1175
Geisler (1970) reported that the e!ect of vibration modes is dominant at the upper part of a huge launch vehicle such as Saturn-V; and in fact, even for an ordinary vehicle size, such as that for Japanese rockets, lateral vibration has a substantial e!ect on the structure of a launched satellite. However, if the system is stable, the vibration modes do not greatly a!ect the launch vehicle structure such that the #ight system design is executed mainly to suppress the magnitudes of the attack angle a and gimbal de#ection angle b, respectively. The e!ect of lateral vibration is discussed in Section 4. In-#ight wind is the main factor causing large values of a and b, and an e!ective method to suppress them is to set the vehicle attitude command such that it is aligned with the relative velocity so as to anticipate the e!ect of wind. Since the discussion on pitching dynamics is applicable to yawing with minor changes, only pitching motion is considered using the #ight notations shown in Fig. 2, where the direction of the vehicle velocity relative to earth rotation is indicated by #ight path angle c. In order to realize a zero angle of attack, the direction of vehicle axis h must satisfy condition h"c!a , (17) U where a is the contribution of in-#ight wind to the angle U of attack. Under this condition, the gimbal de#ection angle b is also maintained at zero, with the resultant #ight mode being referred to as `zero-lift turna. Fig. 3 shows the average monthly wind speeds and directions for Kagoshima prefecture where Japan's launch sites are located. Note that the wind direction at altitudes higher than 3 km is almost exactly west. Rocket launchings in Japan are allowed only during two seasons: the windless summer and windy winter. Since the angle of attack caused by wind a is nearly zero U in summer, but reaches nearly 103 in winter, the respect-
Fig. 2. Notations for pitching dynamics.
Fig. 3. Average monthly wind speeds and directions for Kagoshima prefecture.
ive attitude angle programs required for zero-lift turn are naturally quite di!erent. During launch, the attitude command is given in a series of step functions of attitude rate. An example of the actual pitch command rate hQ used for winter launching of the H-I vehicle is shown in A Fig. 4. The H-I lifts o! vertically and 8}12 s after lift-o! takes a kick turn at a constant pitch rate. Then the zero-lift turn is executed for about 90 s until aerodynamic pressure su$ciently decreases. In the task of #ight analysis, the attitude commands are determined via numerical calculations to minimize fuel consumption from lift-o! to orbit injection under given constraints. Fig. 4 also shows #ight pro"les for c and a , obtained as the result of the U #ight simulation with the attitude rate command hQ . It is A apparent that the zero-lift turn rule is observed well. Similar results were obtained for the yaw rate command by replacing #ight path angle c with #ight direction angle t. If the in-#ight wind coincides with the average wind used to calculate the attitude rate programs, then the aerodynamic loads on the vehicle are maintained at close to zero. However, this is not the case for all #ight altitudes, and actual in-#ight wind has turbulence with high spatial frequencies which are smoothed out in the average wind. It is therefore necessary to con"rm that the loads remain within allowed limits during all possible wind conditions. The vehicle's response to an unexpected wind condition varies depending on the design of the attitude control system. Such a condition can be described using
H. Mori / Control Engineering Practice 7 (1999) 1163}1175
1167
Fig. 5. Signal #ow diagram for the pitching control system.
Fig. 4. Pitch command, #ight path angle, and angle of attack pro"les for a winter launching of the H-I.
perturbational equations derived around the nominal equations of motion. The zero-lift turn trajectory is taken as nominal motion, u, *c, and *a are the perturbations U of h, c, and a from nominal trajectory, and a and b are U themselves perturbed quantities. Under these conditions, Eq. (6) represents the perturbational equations for such pitching motion. The lateral displacement z is related to *c by z"<*c,
(18)
with Eq. (17) yielding
"*c!*a . U
(19)
The sensors used for the pitch attitude control system are a position gyro, rate gyro, and accelerometer. Fig. 5 shows the signal #ow diagram of the control system for pitching motion, where a and a s are the feedback gains of the position and rate gyros, respectively. Broken lines show the feedback paths of the accelerometer and g is the feedback gain, while other lines show the pitching dynamics represented by Eqs. (6), (18) and (19). Note that the line of Z$ , which represents coupling ( between the lateral and longitudinal motions, has been added to Eq. (6). Accelerometer feedback is not typically used, although severe #ight conditions make such feedback necessary, i.e., it can give the control system the characteristics of drift minimum and load relief for the wind; concepts developed by Hoelker (1961).
Under the signal #ow system (Fig. 5), the transfer functions from *a to zK /< are given by U zK /< A s#A s#B " . (20) *a B s#B s#B s#B U For a step function of *a , the lateral acceleration zK /< U reaches the input level unless B "0; a case in which the stationary value of the transfer function becomes A /B which is 1. The relation for B "0 is (a !g Z$ )(Z$ U$ !Z$ U$ )!Z$ U$ "0, (21) ( ? @ @ ? ( ? being called drift minimum condition, which is only realized when accelerometer feedback is present. Fig. 6 shows the typical gain schedule for an H-II vehicle with accelerometer feedback. Transfer functions from *a to a and b are both U represented as A s#A s#A . (22) B s#B s#B s#B Thus, for a step function of *a , both a and b normally U have stationary values. If the position feedback gain a "0, A s for both parameters become zero; and in such case, their stationary values are reduced to zero. The relation a "0 is called the load relief condition. Even if a is not exactly zero, the e!ect of load relief can be expected to occur. The application of the drift minimum condition leads to a smaller value of a , which is useful for the purpose of drift minimum and load relief as well. Eq. (22) indicates that a and b have stationary values for a step function of *a , although they become zero for U a step function of *a . In other words, a severe wind U for launch vehicles has large step functions of *a . U A synthetic wind is used to con"rm that the #ight trajectory does not incur aerodynamic loads beyond the vehicle structure limit. Fig. 7 shows an example of synthetic wind which has the peak at 9 km in altitude. The wind in the upper atmosphere is measured using balloon devices, and the measured data is considered independent along 500-m altitude intervals. At each interval, a 95% level wind is de"ned as one whose magnitude is on the 2.45p probability ellipse in the #ight
1168
H. Mori / Control Engineering Practice 7 (1999) 1163}1175
Fig. 6. Gain schedule for the H-II vehicle.
Fig. 8. a and a for a winter launching of the H-I. U U
and its time derivative a and a , where it is apparent U U that a for the synthetic wind shows a large step function U at an altitude of 9-km in comparison with a of the U average wind. Thus the synthetic wind causes large values of a and b. The angle of attack a and the gimbal de#ection angle b generated by the synthetic wind are used to calculate the bending moment using Eq. (8), with Fig. 9 showing the resultant bending moment and axial force acting on the H-I vehicle. Several synthetic winds are formed with peaks distributed from 2 to 15 km, being used to con"rm the loads imposed on the vehicle are below allowed limits. The total loads are described in terms of the equivalent axial force, i.e.,
Fig. 7. Synthetic wind with peak velocity at an altitude of 9 km.
direction. There are two 95% levels at each altitude: high and low. The 95% pro"les shown in Fig. 7 are obtained by connecting measured 95% level winds along successive 500-m altitude intervals. The wind shear for each measurement is taken as the di!erentiated value of wind with respect to altitude, with 95% wind shear being similarly de"ned. The synthetic wind that peaks at 9-km is formed by connecting 95% wind shears along successive 500-m altitude intervals such that the peak wind occurs at that altitude. Fig. 8 shows the angle of attack due to the wind
4 P (x, t)"P(x, t)$ M(x, t), CO D
(23)
where P(x, t) is the axial force, M(x, t) the bending moment, and D the diameter of the beam. A negative sign indicates compression.
4. Mode analysis and time response Three-dimensional (3D) computation of the mode frequencies and mode shapes of a launch vehicle is carried out using a model composed of a number of series of beam elements, each of which has 5 degrees of freedom, i.e., x, y, z, h and h . Torsion is not considered in the W X
H. Mori / Control Engineering Practice 7 (1999) 1163}1175
1169
Fig. 9. H-I bending moment and axial force for a 9-km peak synthetic wind.
model. The axial displacement x is not coupled with lateral displacements or rotations, while lateral displacement y is only coupled with h , and z only with h . X W Accordingly an element has three non-coupled members whose degrees of freedom are in x, y!h , and z!h , X W respectively. Members are represented by an equivalent mass/inertia-spring component, and are connected to similar members of adjacent elements. Under this approach, the vehicle model is basically represented by three series of members in axial, pitch, and yaw axes. These three series are connected to each other at the engine sections, guidance section, and satellite model, while also being coupled as a whole system. About 500 degrees of freedom are considered "rst. Application of the constrained normal mode method reduces the degrees of freedom of the model to about 100, which are subsequently used in the computation of mode frequencies and shapes. Fig. 10 shows the distribution of the mode frequencies of the H-I vehicle through a #ight, where lateral modes appear in the lower frequency range and axial in the higher. Each mode shape is given in 3D coordinates of pitch, yaw, and axial. The dimension in which a mode mainly belongs can be identi"ed by the dominant coordinate, although more than a few modes do not have a speci"c dominant coordinate. The "rst axial mode appears in the 19th order of frequency, after the 6th pitch and 6th yaw modes, respectively, appear in the 17th and 18th orders, while the second and third axial modes appear in the 24th and seemingly in the 44th, respectively. It should be noted, however, that (i) this model does not show good accuracy in the higher orders and (ii) changes in the trend of modes greatly depends on the time after lift-o!. Since much work is required to construct models for di!erent #ight times, this makes it di$cult to predict an accurate
Fig. 10. Distribution of H-I mode frequencies.
trend of modes for the high-frequency range; a problem discussed later in relation with vehicle longitudinal stability. The model composed of beam elements is also used to compute the response of the satellite to shocks and oscillations generated by the launch vehicle. The greatest load imposed on the satellite occurs during launching, being the load which prescribes basic requirements for the satellite structure. Naturally then, it is important to precisely predict the response of the satellite to such shocks and oscillations. Eq. (16) indicates that the axial force P(x, t) is the sum of the weight of the upper part of the structure * produced by static acceleration of the vehicle * and inertia forces of the distributed mass which are related to the general coordinates due to modal oscillations. At the time of main engine cut-o! (MECO), static acceleration ultimately increases. Moreover, during H-I launches, longitudinal oscillations called MECO POGO occurred due to resonance of the propellant feed system and engine chamber. Hence, the load at MECO prescribed the main
1170
H. Mori / Control Engineering Practice 7 (1999) 1163}1175
design requirement. Fig. 11 shows the axial acceleration at the "rst engine section of the H-I vehicle. Another resonant oscillation (20}33 Hz) called preMECO POGO also appeared during H-I launches. The satellite launched by the H-I vehicle was required to have its own 1st longitudinal mode frequency higher than 35 Hz and the 1st lateral mode frequency higher than 15 Hz. The pre-MECO pogo, whose frequency was close enough to the satellite's lowest longitudinal frequency, imposed a serious burden on the structure design. The dynamics of POGOs are treated in Section 6. If the satellite is considered to have a simple beam structure } normally being more complex } then the oscillation produces an inner force to the satellite represented by the last terms of Eqs. (8) and (16). Thus subsequent time response computations using the employed element model were necessary to examine loads on each part of the satellite. The results were re#ected in vibration tests of the structure models and actual satellites. Extreme shocks act on the vehicle at the time of the lift-o! and at max q a, i.e., the maximum value of the product of the dynamic pressure and angle of attack. For the H-I, max q a occurred at an altitude of 9 km. These shocks had lateral components and stimulated the bending modes shown in Fig. 10. Those transmitted to the satellite had a high spectrum in the lower frequency
range ((15 Hz), being close to the satellite's lowest lateral mode frequency. This is another reason making it essential to carry out time response investigations using the developed model.
5. Stabilization of lateral motion The basic equation of motion for lateral vibration of a vehicle in each bending mode is represented by Eq. (7), where through control feedback, the motions of bending modes are coupled with each other and with the rigid body. Moreover, the control system considers the sloshing motions of propellants and engine-actuator dynamics, as these characteristics must be considered together for the stability analysis of the vehicle's lateral motion. The total equations reached almost 30th-order, and included at least "ve bending modes and two sloshing modes. The rigorous equations of motion for a highly #exible vehicle can be derived using Lagrange's equations, and an example of such derivation can be seen in the work of Luckens, Schmit and Broucek (1961). The derived equations of motion can be written as (Ms#Bs#K) x"k b#k w, (24) @ U where M, B and K are matrices, x is a vector whose component takes one degree of freedom, b is the gimbal de#ection angle, w is the perturbation of the wind and k and k are scalar coe$cients. Fig. 12 shows a typical @ U feedback system used to represent the lateral motion of a vehicle, in this case that of the H-I. Feedback loops can be represented mathematically as b"G(s)x,
(25)
where G(s) is a column vector whose components are rational functions of s. Stability is examined in time-slice pattern along the trajectory using either Nyquist diagrams or Bode plots. Feedback gains are determined as described in Section 3, and they should not be altered for stabilizing purposes. That is, since the natural frequency of rigid body motion is low enough compared to those of other modes, it is not di$cult to stabilize the system by shaping frequency characteristics; and a combination of lead and lag compensations is generally used.
Fig. 11. Longitudinal oscillations measured at the engine section during H-I launchings.
Fig. 12. A block diagram of the H-I's guidance and control system.
H. Mori / Control Engineering Practice 7 (1999) 1163}1175
As the compensating "lters in the feedback loops are associated with the inertial guidance unit (Fig. 12), the "lter transfer functions must be changed into the digital algorithms via bilinear transformations, i.e., 2 1!z\ D(z)"D(s)" , f (z)" , QD X ¹ 1#z\
(26)
where z"e2Q and ¹ is the sampling period for digital computation. The digital algorithms obtained via Eq. (26) are accurate for a frequency range below ¹. Given that the sampling period for the H-I was 25 ms, the digital "lters were accurate for a frequency range below 5 Hz, although they might also be e!ective for stabilization at a higher range. For example, Goto and Suzuki (1981) designed the H-Is digital "lters and con"rmed their e!ectiveness below 25 Hz. After passing through the atmosphere, the guidance loop becomes active, and the stability of the system during this phase of #ight, which includes both guidance and control loop's, must be examined since the computation period for guidance is longer than that for control. Fig. 13 shows a block diagram of the sampled-data system incorporating two di!erent sampling periods, i.e., ¹ and ¹ are, respectively, the control and guidance Q E loop sampling periods * z "e2Q and z "e2E * while Q E D and D are the compensating "lters and H and H the Q E zero-order holds for each sampling. When the guidance loop sampling period is an integer multiple of the control loop's, the characteristic equation of the sampled-data system is represented by
H FD (z )[H ]H2Q H2Q Q Q E Q(z )"0. (27) E 1#D (z )[H F]H2Q E Q System stability can be examined using the Nyquist diagram of the open-loop transfer function of Eq. (27), with this stability analysis allowing the gains of the guidance to be con"rmed. In Japan, most of solid motor rockets use tail "ns for aerodynamic stability. Since strong control forces are not present, such rockets can #y without an attitude control system, and do not necessarily fail if they happen to be unstable during #ight. In fact, several #ights were suc1#
Fig. 13. Sampled-data system incorporating samplers for two di!erent periods.
1171
cessful even though the control system was found to be unstable. In comparison with these conventional vehicles, the solid-motor M-V vehicle, newly developed by the Institute of Space and Astronautical Science (ISAS), has no tail "ns and is inherently aerodynamically unstable. It is obliged to be equipped with a control system for stability, which works in the presence of uncertain system parameters such as the bending rigidity, the aerodynamic coe$cients and the actuator bandwidth. Morita et al. (1994) applied the H control theory in the construction of its feedback system; and under a framework of the mixed sensitivity problem, they achieved better control performance for a relatively low-frequency region and preferable robust stability in a higher one. A similar application of the theory is investigated by Mau!rey and Schoeller (1998) for a launch vehicle of the Arian type with liquid propellants. 6. Stability analysis of londitudinal vibration As described in Section 4, longitudinal vibration is the main factor determining satellite environmental conditions during launch. Three types of POGO resonant oscillations have appeared in N-II and H-I rocket #ights. MECO POGOs appear just before main engine cuto!, being caused by #uctuations in LOX tank bottom pressure, and two kinds of pre-MECO POGOs sometimes appearing in the midst of #ight. Although the latter POGOs have not been well analyzed, possible causes were investigated by Mori (1991) using a simpli"ed analysis method in conjunction with #ight data. Here an outline of the analysis is described. Eq. (15) represents the basic equation of longitudinal motion of a vehicle, where a single mode is excited by resonance with thrust variation. By adding structure damping, the equation of motion becomes ; (x )C A uK #2fu u #uu " G @ 2 RFp , (28) G G G G G A M G where C is the thrusting coe$cient, A the nozzle throat 2 RF area of the engine, and p the chamber pressure. During A resonant oscillation, p is a!ected by structure vibrations. A Fig. 14 shows signal #ow diagrams for such dynamics,
Fig. 14. Signal #ow diagram representing Eq. (28).
1172
H. Mori / Control Engineering Practice 7 (1999) 1163}1175
where [ p /uK ] is the transfer function from uK to p , while A G G A K is the gain from p to uK . 2 A G The open-loop transfer function of this system is given by
u u p A . O¸(s)"2f G# G !K G s 2 uK s G
(29)
Letting s"ju, Eq. (29) takes the form
O¸( ju)"!
u p G #K Re A 2 u uK G
u p !j 2f G#K Im A Gu 2 uK G
.
(30)
The Nyquist diagram of Eq. (30) is represented conceptually by Fig. 15. That is, when the structure is not excited by chamber pressure, i.e., when K 0, the Nyquist diagram is in2 dicated by the broken line. In this case, for u"u , G Re(O¸)"!1, Im(O¸)"!2f , G and the systm is stable. On the other hand, when K 0, 2 there is a possibility that the last term of Eq. (30) becomes large such that, as shown by the solid line, the diagram rises up over (!1, 0) depending on the value of Im[ p /uK ]. It indicates the system is unstable. The transA G fer function [ p /uK ], whose notations are indicated in Fig. A G 16 and applied to the fuel and the oxygen feed system, respectively, is determined as follows. A pressure perturbation at the bottom of the tank is given by p "oh; (x)su (s), R G G
(31)
where o is the density of propellant and h the height of propellant in the tank. The characteristics of #ow and pressure transmitted through the pipe line is represented by ¸ cos uq ) p # J sin uq ) s Q "p , Q uq Q R
(32)
Fig. 16. Notations for pogo stability analysis.
where p and Q are, respectively, perturbations of presQ Q sure and #ow at the outlet of the feedline due to oscillation; and ¸ is the inertance of propellant in the feed line J and q the time required for the sound wave to pass through it. The #ow perturbation at the outlet of the accumulator, Q is represented by Q 1 Q !Q " sp , (33) Q Q ¸ s#R s#K Q ? ? ? where R is the resistance, ¸ the inertia, and K the ? ? ? volume elasticity of the accumulator, respectively. The #ow to the discharge line is represented by Q "Q !C sp #A su(x , s), (34) B Q @ Q J N where C is the cavitation compliance and A the cross@ J sectional area of the feedline. The perturbation in chamber pressure is expressed as p "ap !(¸ #¸ )sQ #(R #R )Q , (35) A Q N B B N B B where a is the ampli"cation coe$cient of the pressure variation from the pump inlet to outlet, and ¸ , ¸ , and N B R and R are the inertance and resistance of the pump N B and discharge line, respectively. The propellant #ow raises the chamber pressure in accordance with
Fig. 15. Nyquist diagrams of the open loop transfer functions of Eq. (30).
cHo (1#q s)p " Q, (36) A A A B RF where q is the time constant of the combustion and A cH the characteristic exhaust velocity. The transfer function [ p /uK ] is obtained eliminating "ve variablesA G p , p , Q , Q and Q * from Eqs. (32)}(36). The bubble R Q Q Q B
H. Mori / Control Engineering Practice 7 (1999) 1163}1175
1173
angular velocity u , i.e., @ u "1/(C ¸ (37) @ @ R is considered an important parameter related to POGO phenomena. This is also expressed in the term of bubble frequency by f "u /(2n). If we consider a simpli"ed case @ @ in which no accumulator is installed, and further assume that R"R #R "0 and q "0, then Eq. (30) becomes N B A O¸( ju)
"!
u cd G #K Ka 2 N u *
u c¸ #a¸ tan uq/(uq) J !j 2f G#K Ka NB du , Gu 2 N * (38) where *"(cK)#+c¸ #a¸ tan uq/(uq),u, N NB R u tan uq , c(u)"1! u uq @ o sin uq d(u)" h; (x )# l; (x ) . G R G N cos uq uq
The unstable condition shown in Fig. 15 occurs when c(u ) is nearly zero, while the condition c(u )"0 is satisG G "ed when u equals `the resonant bubble angular velo@ citya - , i.e., @ tan u q @ " G . (39) u uq G G - is a good parameter for indicating which bubble @ frequency produces vehicle POGO, although some simpli"cation is applied. This is also expressed in the term of resonant bubble frequency fM by @
tan 2nf q G . (40) 2nf q G Fig. 17 shows the #ight time periods and frequencies of typical POGO phenomena appearing during three H-I launches (*, 䉭, and ;). Two pre-MECO POGOs, whose frequencies are indicated by f and f , are named PMP-1 and PMP-2, G G respectively, while MECO POGO related with the "rst longitudinal frequency f is abbreviated to MP. The resonant bubble frequency fM corresponding to each @ POGO frequency f was determined using Eq. (40), with G Table 1 summarizing values of fM corresponding to the @ central frequency of each POGO. The value of q for the LOX and fuel system is 0.26;10\and 1.72;10\ s, respectively. Since Eq. (40) stands for 2nf q n in the fuel system, the fM for PMP-2 G @ can be markedly smaller than f . While POGO oscillaG fM "f @ G
Fig. 17. Pogo frequency/bubble frequencies of the H-I.
tions during PMP-1 and MP can occur only when 2nf q 0 for the LOX system. G Since the fuel system is not equipped with an accumulator, Eq. (40) is expected to have a good accuracy and f likely coincides with fM . In the accumulator-equipped @ @ LOX system, however, fM is not expected to have a close @ value to f . Use of Eq. (37) to directly determine bubble @ frequencies is problematic because it is di$cult to determine cavitation compliance accurately; hence, bubble frequencies were estimated such that stability conditions match the #ight status. To compute Eq. (25) using Eqs. (26)}(31), the oscillated (POGO) mode shapes and their generalized masses must be known. Using measured #ight data for each abovementioned H-I launch, Fig. 18 shows the longitudinal mode shapes corresponding to POGOs, i.e., PMP-1-1, PMP-1-2, PMP-2, and MP. PMP-1 are divided into two sections of time. Longitudinal accelerations were measured at indicated positions (AF,GS,2GP, and IES). Also shown is the mass distribution at lift-o!. The mode shapes of POGOs in MP and PMP-2 were calculated using the model described in Section 3. The mode shapes for those in PMP-I, however, were not calculated as these oscillations did not critically impact the satellite; instead they were estimated based on measured data. With mode shapes and the mass distribution in hand (Fig. 18), time-varying values of the generalized masses can be determined for each POGO, with results being shown in Fig. 19. Next, bubble frequencies are estimated in order to explain the pogo phenomena. Table 2 shows bubble frequencies corresponding to POGO frequencies of Table 1.
1174
H. Mori / Control Engineering Practice 7 (1999) 1163}1175
Fig. 19. Generalized masses corresponding to longitudinal mode shapes in Fig. 18.
Table 1 Resonant bubble frequency corresponding to respective POGO frequency Fig. 18. Longitudinal mode shapes of POGOs during H-I launches. PMP-1 PMP-2 MP
Resultant values were determined by trial and error. The broken lines in Fig. 17 indicate time-varying values of bubble frequencies for the LOX and fuel systems. Finally, the estimated bubble frequencies help in determining respective Nyquist diagrams at various #ight times. Fig. 20 shows the result, where the stability characteristics represent the time-varying POGO phenomena of Fig. 11 well. No POGO-induced longitudinal vibrations have occurred during launch of the H-II, which unlike the H-I has its lox tank located above the LH2 tank. While this change in con"guration is thought to prevent MECO POGOs, the possibility of Pre-MECOs still exists, although none have been observed thus far. The design of the POGO suppressor (i.e., accumulator) installed in the H-II may be e!ective.
7. Concluding remarks The results obtained using the employed analytical model indicate that the design of a control system for a #exible-body launch vehicle is highly interrelated with trajectory planning and structure analysis consideration.
f (Hz) G
LOX fM (Hz) @
Fuel fM (Hz) @
22.0 32.0 17.5
22.5 * 17.7
* 9.85 *
Table 2 Bubble frequency corresponding to respective POGO frequency
PMP-1 PMP-2 MP
f (Hz) G
LOX f (Hz) @
Fuel f (Hz) @
22.0 32.0 17.5
33.5 * 23.8
* 9.9 *
In addition, through the use of analytical methods, critical #ight characteristics are elucidated, i.e., aerodynamic loads acting on the launch vehicle, lateral and longitudinal stability, and launching environments of the housed satellite. Based on a stability analysis of POGOs, it is shown that their occurrence is estimated by examining resonant bubble frequencies, around which bubble frequencies of the feed systems are determined such that time-varying POGO phenomena of the H-I can be represented by Nyquist diagrams.
H. Mori / Control Engineering Practice 7 (1999) 1163}1175
1175
References
Fig. 20. Nyquist diagrams for indicated POGO periods during H-I launches.
Frosh, J. A., & Valley, D. P. (1967). Saturn AS501/S-IC #ight control system design. Journal of Spacecraft, 4(8), 1003}1009. Geisler, E.D. (1970). Wind ewects on launch vehicle. AGUARDograpph No. 115. Goto, S., Suzuki, H. (1981). A design of digital yight control system of a launch vehicle. NASDA TR-10, National Space Development Agency of Japan. Hoelker, R.F. (1961). Theory of artixcial stabilization of missiles and space vehicle with exposition of four control principles. NASA TND-555. Johnston, D.E., Johnson, W.A., (1967). Feasibility of conventional control techniques for large highly coupled elastic boost vehicles. STI Technical Report, No.67-592. Lukens, D.R., Schmit, A.F., Broucek, G.T. (1961). Approximate transfer functions for yexible-booster-and-autopilot analysis. WADD-TR-6193. Mau!rey, S., Schoeller, M. (1998). Non-stationary H control law for launcher with bending modes. Preprint of 14th IFAC symposium on automatic control in aerospace. (pp. 438}443). McKenna, K.J. et al. (1964). A method for studying the coupled engine airframe longitudinal instability of liquid rocket systems. AIAA paper 64}81. Mori, H. (1991). POGO analysis based on N-II/H-I vehicle #ight data. Journal of the Japan Society for Aeronautical and Space Sciences, 39, 71}80. Morita, Y., Kawaguti, J., Nakatani, I., Goto, S., Ohtsuka, H. (1994). Attitude control algorithm of the M-V launching rocket, Proceeding of 4th workshop on astrodynamics and yight mechanics. (pp. 169}174). Kanagawa, Japan. Pain, J.G., Rubin, S. (1974). Pogo suppression on the Delta vehicle. Report SAMSO-TR-74-180. Rubin, S. (1966). Longitudinal instability of liquid rockets due to propulsion feedback (POGO). Journal of Spacecraft, 3, 1185}1195.