This article, received April 21 1964, gives the results of an analytical and numerical study of a two-gyro, gravity-oriented communications satellite. The principal purpose of the study was to uncover and solve ihe analytical problems arising in the design of passive gravity-gradient attitude control systems. Although the study was directed at satellite orientation, it is felt that many of the techniques developed have general use in the investigation of dynamical systems.
We consider both small and large motions about the desired earth-pointing orientation. In the small-motion study, the goal is simultaneous optimization of the transient response and the forced response to perturbations caused by orbital eccentricity, magnetic torques, solar torques, thermal rod bending, and micrometeorite impact. In the large-motion study, we enumerate all possible equilibrium positions of the satellite and then consider initial despin after injection into obit, inversion of the satellite from one stable equilibrium position to another by switching of gyro bias torques, and the decay of transient motions resulting from large initial angular rates.
As a specific numerical example, we hae treated a 300-lb satellite in a 6000-nm orbit, stabilizcd by a 60-ft extensible rod with a 20-lb tip mass, and by two single-degree-of-freedom gyros, each with an angular momentum of 106 cgs units. Without a detailed discussion of hardware, it is concluded that such a system, having a total weight of 50 to 75 pounds including povwer supply, will provide a settling time for small disturbances of less thau one orbit and will hold the antenna pointing error within a few degrees.
A computer animated movie showing the system is available.
It has been known for over two hundred years that the variation in the gravitational field over the length of an earth satellite generates torques which tend to keep the axis of minimum inertia of the satellite pointing toward the earth. In particular, this mechanism keeps one face of the moon earth-pointing.
Such gravity-gradient orientation of communications satellites is very attractive because the simplicity of the effect leads to the possibility of simple attitude control and hence high reliability and long life. On the other hand, the tiny size of the gravity-gradient torques means formidable mechanization problems, and although Pierce suggested its use as early as 1955 [1] gravity-gradient stabilization has been widely held to be impractical.
However, several recent analytical and hardware studies have resulted in proposals for practical, gravity-gradient controlled satellites. All the proposed schemes work on the same principle. Steady-state perturbations, due, for example, to magnetic and solar torques, are kept within tolerable limits by making the satellite inertia sufficiently large, usually with some sort of extensible rod-tip mass combination. Damping of transient perturbations is provided by connecting the satellite through a dissipative joint to an anchor, that is, to some object that will allow energy dissipation by virtue of relative motion between itself and the satellite proper. The anchor may be one or more gyros, as in the schemes discussed by Ogletree, et al [2][3][4], by Burt [5], and by Scott [6]; a second rigid body, either hinged to the satellite, as proposed by Kamm [7], by Paul, West, Yu, et al [8][9]; or a second rigid body at the end of a compliant dumbbell as discussed by Paul [10], by Newton [11], and by Fischell and Mobley [12]; or a second, fluid body, as considered by Lewis [13].
In this article we examine a gravity-gradient system anchored by two gyros. A schematic of the system is shown in Fig. 1, where also is indicated the standard nomenclature for axes: the pitch axis is normal to the orbit plane, the yaw axis is along the local vertical, and the roll axis is along the orbital track. Each gyro rotor is contained in a gimbal can (not shown in the schematic), mounted on bearings, and immersed in a fluid bath. Thus, fluid shear produces the required energy dissipation. The gyros are single-axis gyros: that is, the spin vectors are constrained by the gimbal bearings to lie in a single plane within the satellite. In the position shown, this is the pitch-yaw plane.
Because of the small physical dimensions of the gyro anchor, this system has the virtue that the dissipative joints can be sealed within the satellite; the joints are not exposed to the space environment when the satellite is operating. Also, the required inertia augmentation is particularly simple: only a single extended rod-tip mass.
A simple explanation of how single-axis gyros damp out an arbitrary motion can be given in terms of the rate or torque-seeking property of gyros. A torque applied to a gyro will cause it to precess. By conservation of angular momentum, the precession will try to line up the gyro spin vector with the applied torque or angular rate vector.
Bearing this in mind, assume that the satellite is in orbit in its earth-pointing orientation. It then is rotating at the rate of one revolution per orbit about the pitch axis. If the gyro gimbals were free, this pitch rate would cause the spin vectors to align themselves in the direction of the pitch axis. However, in order to obtain three-axis damping, the spin vectors are held in a vee position by equal and opposite constant torques (see Fig. 1), applied to the gimbals.
Now, if the satellite is disturbed about the pitch axis, both gyros seek the disturbance, resulting in a scissoring motion of the gimbals relative to the satellite, damping out the pitch disturbance. A yaw disturbance causes an in-phase motion of the gyros and again energy is dissipated. Since the gyro spin vectors are constrained to move in the pitch-yaw plane, they are constrained from moving toward a disturbance about the roll axis. However, the roll and yaw motions are coupled. Hence in this case the gyros again try to line up with the yaw axis. Thus three-axis damping is obtained.
Our work continues a study carried out by the Instrumentation Laboratory [2][3] of the Massachusetts Institute of Technology, under the sponsorship of Bell Telephone Laboratories, in which the particular two-gyro configuration studied here was shown to be the most promising of several possible gyro-anchored systems. Our primary objective, however, was not to design a specific attitude control system, which in any case would have to he integrated with the design of a specific satellite, but rather to develop general guiding principles and analytical and numerical techniques useful in such a design problem. Thus, we consider only the broad hardware questions that affect the analysis for example, the design of extensible rods necessary to augment satellite inertias - but we do not go into the detail of specific gyro hardware, as would be required in a complete design.
The organization of the article is as follows. Section II, for the general reader, summarizes the results of our study in some detail in non-mathematical terms. Following some general remarks about inertia levels, applicable to all gravity-gradient systems, we more fully describe the two-gyro system studied. We then summarize the system's small-angle performance, stressing, in particular, the performance obtained when the inertia of the satellite is augmented by the erection of a single rod. Next we discuss the effects of the main small-angle perturbations: orbital eccentricity, magnetic torques, solar radiation pressure, micrometeorite bombardment and thermal rod bending. Finally, we consider large-angle motions, starting with initial despin upon orbital injection by a combination of rod erection and unengaging of the gyros. In the discussions of large-angle motions, we indicate that there may exist equilibrium positions far removed from the desired, earth-pointing position; we also show how these may be avoided.
Gravity-gradient systems are bistable: that is, associated with a stable, earth-pointing orientation is a second, equally stable orientation obtained by a 180° rotation about the pitch axis. In the concluding section of Section II we describe how the satellite can be flipped from one stable orientation to the other by means of a torque pulse applied to the gimbals.
The results pertaining to the two-gyro system given in Section II serve as an outline of the analysis required for the design of any gravity-gradient attitude control system. They also serve as an introduction to the theory in Section III and Section IV. In these parts we present several results and methods that we feel apply generally to the design of many-parameter, linear dynamical systems (see Section III) and to large-angle motions of a satellite (see Section IV).
Specifically, in Section III we develop various bounds on system settling time, and then show how series expansions in terms of system parameters can be used to explore the behavior of a linear system as a function of its parameters. We next describe a computer program based on the Routh criteria, which allows very rapid computation of system response as a function of system parameters. By these means, we are able to survey system behavior over the entire range of six system parameters.
In Section IV, we develop the equations of large-angle motion, including the case of variable inertia, occurring during rod erection. Here we stress the superiority of direction cosines or Euler parameters as compared to Euler angles in satellite kinematics, both from the point of view of computing speed and of ease in visualizing satellite motions. We then give the analysis of equilibrium positions, despin, and flipping or inversion.
All gravity-gradient systems have one feature in common, namely the low magnitude of the gravity-gradient restoring torque, of the order of IΩ2, where I is a typical satellite moment of inertia and Ω is the orbital rate. The level of this torque is the main factor determining the steady-state response to constant and periodic disturbing torques. In particular, in the case of a typical communications satellite at an altitude of 6000 nm, the magnitude of the torque exerted by the geomagnetic field on the residual magnetic moment of the satellite is such that the satellite inertia must be increased by a factor of about forty to reduce the steady-state response to an acceptable level.
The low level of the gravity-gradient restoring torque also implies low system natural frequencies, of the order of the orbital rate Ω. Corresponding to this low natural frequency is a minimum I/e settling time of the order of a fraction of an orbit. Zajac [14] has shown that all the systems mentioned above have pitch settling times no less than about one tenth of an orbit. This, of course, is a lower bound on minimum settling time for three-axis motion.
Based on these simple considerations, we would expect that all well designed gravity-gradient attitude control systems would have about the same transient and steady-state performance, that they would all have settling times of a fraction of an orbit, and that they would all require some form of inertia augmentation to obtain acceptable steady-state response. Thus the choice of a particular gravity-gradient attitude control system should be based mainly on ease of mechanization and longtime reliability, rather than system performance.
In the present case, at least, requirements on the large-angle performance of the system (despin, satellite inversion, etc.) preclude choosing the system parameters to give minimum settling time, although the settling time is not greatly increased by meeting the other requirements. It is likely that such a compromise would be necessary for optimum overall performance of any gravity-gradient system, so that the minimum settling time is of academic interest only. Of more importance is the variation of system performance with variation in system parameters. We have thus taken the view that a broad survey of performance as a function of system parameters is of more interest than an optimization based on a single measure of system performance, e.g., settling time.
In the following sections we describe the configuration and performace of the two-gyro system in detail. The interested reader may find the corresponding theoretical analyses in Section III and Section IV.
Fig. 2 shows the more important features of the typical single-axis gyro indicated schematically in Fig. 1. The basic element of the gyro is a rotor which spins rapidly about the spin axis and generates a certain angular momentum vector.
The spinning rotor element is enclosed in a sealed gimbal can, mounted on bearings so that it can rotate about a single axis, the gimbal or output axis. A fluid-filled gap between gimbal can and gyro case provides damping as the gimbal rotates.
In the system considered, the two gyros have their gimbal axes along the satellite roll axis. The gyro spin axes are disposed in a vee configuration around the satellite pitch axis, which is also the axis about which the satellite rotates to remain aligned with the local vertical as it traverses its orbit. To distinguish this arrangement from other possible two-gyro configurations [2] [3] [4] [5] [6] , it will be called a roll-vee configuration.
In the vee arrangement, torques must be supplied constantly to change the direction of the gyro angular momentum vectors, as the satellite traverses its orbit. These torques, constant in magnitude and exerted about the gimbal axes, are provided by a constant electrical signal into electromechanical torquers on the gimbal axes.
It is also possible to inject a signal into the torquers on ground command. This can be used to invert the satellite if it should get into an undesirable equilibrium position. This possibility is discussed in the sequel.
In order to spin the gyro motor, current must he brought into the gimbal can. This is done by means of highly compliant flex leads. In the present application, the flex-lead spring constants, exerting a small restoring torque around the gimbal axis, can be neglected. However, for a typical communications satellite without inertia augmentation, the flex-lead torques can be of the same order as the gravity-gradient torques.
In any case the gimbal excursions must be limited by suitably placed stops. The location and nature of these stops is an important design consideration. In the first place, undesired equilibrium positions, with the gyro gimbals against the stops, may occur if the stop positions are not carefully chosen. In the second place, large tumbling rates may force the gimbals against the stops, wbere they are capable of only limited relative motion, depending on the stop elasticity. In botb cases the available damping may be greatly reduced. The equilibrium positions may be dealt with analytically, while the large motion may be studied numerically witb the stops simulated by hardening springs.
For the attitude control of a typical communications satellite in a 6000-nm altitude orbit, we require two single-axis gyros, each with a rotor angular momentum of about 106 cgs units, weighing about 10 pounds and requiring from 7 to 10 watts power to drive the rotor motor. In addition we require some sort of inertia augmentation which we shall assume is supplied by a single extensible 60-foot rod of the STEM (self-storing tubular extensible member) type, designed and developed by DeHavilland Aircraft of Canada, Ltd., and described in detail in Ref. 8, together with a 20-pound tip mass, which also serves as the tape storage drum. We then have the attitude control system weight breakdown given below:
2 106 cm-gm-sec gyros 20 lbs
1 tip mass 20 lbs
1 extensible rod 4 lbs
gyro power supply
(2 lbs of solar cells/watt) 40 lbs
total 84 lbs
We bave assumed that the satellite proper is a four-foot diameter sphere, weighing 300 pounds, with a moment of inertia of 20 slug-ft2.
It is believed that the above estimates are quite conservative and subject to considerable reduction. The power is used to maintain the gyro rotor speed constant mainly against bearing drag. In a zero-g enviroument the bearing drag might he substantially reduced. In any case there is probably a trade-off between gyro life, requiring heavily lubricated bearings, and minimum rotor power, requiring light lubrication.
The rod length is chosen to increase the satellite pitch and roll inertias from 20 to 2000 slug-ft2. This inertia augmentation sufficiently despins the satellite from currently estimated injection rates of 0.5-1.0 rpm to cause capture by the gravity-gradient field. The required inertia augmentation varies roughly linearly with initial injection rates (see Section 2.4.1). With a sufficiently small injection rate, the augmented inertia could be reduced to the 700 slug-ft2 level required to counter magnetic torques (see Section 2.3.4.1). Such a reduction in inertia would mean smaller gyros, and, again, less power.
In order to study the small-angle transient and steady-state response of the roll-vee gyro attitude control system, extensive tables giving decay rates, response to orbital eccentricity, and response to periodic torques at zero, one, and two times orbital frequency Ω as functions of the system parameters were produced by an IBM 7090 computer in a running time of 0.04 hour by a procedure outlined in Section III. Figs. 3 through 15 summarize this broad survey. For each pair of inertia ratios, B/A, C/A, where A,B,C, are the satellite pitch, roll, and yaw moments of inertia, satisfying the inequalities
A ≥ B ≥ C, B + C ≥ A
values of gyro parameters were chosen from these tables to minimize the asymptotic settling time, i.e., the time in which the most lightly damped mode of motion is reduced by 1/e. Fig. 3 shows the corresponding settling times, while Fig. 4 and Fig. 5 give the gyro dimensionless parameters
h = (H/AΩ) cos α, h' = (H/CD) cos α
where α is the vee half-opening angle, H the gyro angular momentum, and CD the gyro damping constant for both gyros. Since the small roll-yaw motion depends only on H in the form H cos α, the above is a convenient choice of parameterization. In all cases, except those indicated, the best value of α was 60°, at least over the relatively coarse grid of Δh = 0.25, Δh' = 0.25 and Δα = 20° used in the tables.
Figs. 6 through 15 give the steady-state response to an orbit eccentricity ε= 0.01 and to periodic torques of amplitude 0.01 AΩ2 for the same values of gyro parameters. Note that the eccentricity response when ε = 0.01 is of the order of 1° over the entire range of inertia ratios, having a maximum value of less than 3°. Both the pitch offset, due to a constant pitch torque, and the roll amplitude, due to a periodic roll torque at orbital frenquency, depend only on the satellite inertias, being given in radians by the simple relations
|φx| 0 = M/3(B-C)Ω2 |φy| 1 = M/3(A-C)Ω2
for a torque of amplitude M. Similarly, for torques at frequency ω ≫ Ω, the pitch, roll, and yaw amplitudes tend to the values M/Aω2, M/Bω2, M/Cω2, again independent of the gyro parameters.
These plots do not show the values of inertia ratios and gyro parameters which yield the smallest settling time. A search over a finer grid of parameter values gives a minimum value of settling time of 0.332 orbits, attained for B/A = 0.925, C/A = 0.175, h = 0.260, h' = 0.688, α= 64°. To attain this value, a slightly negative gimbal spring K = -0.15 HΩ cos α must be used. A negative spring constant may be realized by a simple feedback circuit between gimbal pickoff and gimbal torquer, This minimum settling time is useful as a lower bound, but of more practical interest is the broad range of system parameters over which settling times of less than one orbit can be obtained.
The figures also do not give performance values for a dumbbell or spindle, i.e., a body for which A = B ≫ C. This case is of particular interest, since it may be realized by the erection of a single rod-tip mass combination. To describe the spindle, as well as to give a sample of the tables made by the computer, we reproduce the computer output for B/A=1, C/A=0.01 in Table I. Since the IBM printer has only a limited range of symbols, the following replacements were used:
BB = b = B/A CC = c = C/A KAPPA = κ =1+|K/(HΩ cos α)|
where K is the gimbal spring constant, so that κ = 1 means zero gimbal spring constant,
HH = h HP = h' ALPHA = α
The remaining quantities give the transient and steady-state responses. In particular, P0, P1, P2, R0, R1, R2, Y0, Y1, Y2 are the pitch, roll, and yaw amplitudes in degrees for pitch, roll, and yaw torques of amplitude 0.01AΩ2 at zero, one, and two times orbital frequency. Note that P0 and R1 are constant, since they depend only on b and c, while R0,Y0 are fixed for fixed values of h. The quantity E is the pitch amplitude in degrees at orbital frequency for an orbit eccentricity ε = 0.01. Finally the quantities labeled QUINT and C give the real parts of smallest magnitude of the characteristic roots of the roll-yaw quintic and the pitch cubic in terms of the orbital rate Ω. The smallest of these values D (say) determines the asymptotic settling time Ts= 1/(2 πD). Inspection of the table reveals that, for h=0.750, h'=1.25, α=40°, we have the smallest settling time, for QUINT=-0.340(2), -0.657(2),i.e., two roots with real parts -0.340 and two roots with real part -0.657, and one negative real root (not listed) of larger magnitude. Similarly, in this ease C=-0.279(2), i.e., two roots of the cubic with real part -0.279 and one unlisted negative real root of larger magnitude. The asymptotic settling time is then given by 1/[(2π) (0.279)]=0.57 orbits. Despite the coarseness of the table, this is very close to the minimum value of 0.50 orbits for a spindle, attained for h=0.77, h'= 1.29, α=38°. This minimum value may be calculated by an asymptotic expansion in the large quantity h/c = (H cos α)/CΩ.
Let us now attempt a specific design. This design must be regarded as illustrative, rather than definitive, since a real design must take into account the fine details of gyro hardware as well as requirements imposed by the use of the satellite in an actual communications system. For example, it is not at all clear what limits on maximum settling time would be imposed by system requirements. We have tentatively set this maximum settling time at one orbit.
Fig. 16 shows the asymptotic settling time in orbits as a funetion of the dimensionless gyro angular momentum H/AΩ for α=40°, h'=1.25 and for α=60°, h'=1.00. The former gives a minimum settling time very near the optimum value for a spindle for H/AΩ nearly unity, but varies more rapidly with H/AΩ than does the other system. Also, we are particularly interested in large values of H/AΩ i.e., H/AΩ>2 - since we propose to use the gyros as inertia wheels in the initial despin of the satellite after injection into orbit. ln this case the second system gives a considernbly smaller settling time (0.85 orbits, compared with 1 orbit, at H/AΩ=2). We may actually increase H/AΩ to about 2.4 in this case and stay within the maximum settling time of 1 orbit. Undoubtedly, by trimming the values of α and h, we may increase H/AΩ even more, but, since this is intended to be an illustrative design, we do not considcr these questions further here; instead, we simply take as our design α-60° h'=h= 1.00(H/AΩ=2.00). In the illustrative examples of the sequel, these parameter values will be assumed. From the table, they yield
QUINT = -0.189(2), -1.318(2), C = -0.190(2)
Given the real parts of 4 roots of the roll-yaw quintic and 2 roots of the pitch cubic, it is a simple matter to calculate all the characteristic roots completely, especially for a spindle. In the present case we find solutions of the form
e-0.19Ωt sin (0.64Ωt) e-0.19Ωt cos (0.64Ωt) e-6.62Ωt
for the pitch motion, and
e-0.19Ωt sin (1.4Ωt) e-0.19Ωt cos (1.4Ωt) e-1.32Ωt sin (0.53Ωt) e-1.32Ωt cos (0.53Ωt) e-200Ωt
for the roll/yaw motion. Note that the period of the oscillatory solutions is comparable to the orbital period, as we would expect, and that both motions have rapidly damped exponential solutions. The latter feature, typical of spindle-shaped bodies, causes difficulties in numerical integration of the differential equations, both for small and large motion, for it implies that derivatives may be very much larger than the dependent variable itself.
One source of transient disturbance is the angular momentum imparted by micrometeorite impact. It was estimated in Ref. 8 that for a satellite of comparable inertia level, impacts producing offsets greater than 5° would occur every two years and impacts large enough to tumble the satellite every 23 years, on the average. A more recent study [24] of the present two-gyro system indicates similar times if Whipple's 1958 micrometeorite data are used. For Whipple's 1963 data, the corresponding 5° and tumbling times are 40 years and 1000 years. From a systems point of view, the low frequeney of occurrence of these disturbances suggests that a settling time of 1 orbit is quite adequate.
Disturbances producing constant or periodic pointing errors may be classified as either kinematic or dynamic. The response to the former, of which orbit eccentricity is a typical example, is essentially independent of satellite inertia; the response to the latter type of disturbance may he reduced simply by increasing the level of satellite inertia to a suitable value. The kinematic response limits the minimum attainable pointing error; the dynamic response to given disturbing torques sets the level of satellite inertia.
In the case of the spindle, the table yields an eccentricity response amplitude at orbital frequency of E = 1.81° for an orbit eccentricity ε = 0.01. We also find the steady-state response amplitudes due to torque amplitudes of 0.01AΩ2 given by
P0 = 0.19° P1 = 0.18° P2 = 0.05° R0 = 0.14° R1 = 0.19° R2 = 0.17° Y0 = 0.29° Y1 = 0.56° Y2 = 0.09°
We are particularly interested in the pitch offset P0, due to a constant pitch torque, and the roll amplitude R1, due to a roll torque at orbital frequency. Both of these are independent of gyro parameters and equal in the case of a spindle. Together with a given disturbing torque, they serve to set the satellite inertia level.
In the case of a communications satellite, one of the principal disturbing torques is the torque exerted by the geomagnetic field on the residual magnetic moment of the satellite. It has been estimated in Ref. 8 that this torque might be as large as 5×10-6 ft-lb for a satellite like the Telstar satellite at an altitude of 6000mm. At this altitude Ω=2.73×10-4 rad/sec (∼ 1 rad/hr). Because of the steady rotation of the earth-pointing satellite, this torque does not have a constant pitch component, but it will have a roll component at orbital frequency. Thus R1 is the response amplitude of interest. To make R1 equal to the eccentricity response of 1.81° requires a satellite pitch moment of inertia A such that AΩ2 is ten times the above torque, yielding A = 670 slug-ft2. Since a typical moment of inertia for a satellite somewhat larger than the Telstar satellite is 20 slug-ft2, this calculation indicates that some sort of inertia augmentation is required. We shall assume that the satellite proper has equal moments of inertia A0 = B0 = C0 = 20 slug-ft2 and that the pitch and roll moments of inertia are increased to 2000 slug-ft2 by the erection of a single 60-ft extensible rod and a tip mass of 20 pounds. As indicated in Section 2.2.1, we make the inertia somewhat larger than the minimum required to counter magnetic torque in order that rod erection may he used for satellite despin.
When rods are erected to augment the satellite inertia, the perturbing torque due to incident solar radiation is in general increased, but not at the same rate as the inertia, if a thin rod and a dense tip mass are used. To give some idea of the order of magnitude of this torque, let us consider a 2-foot radius, 300-pound satellite, joined to a 0.2-foot radius, 20-pound tip mass by a 58-foot-long, 0.04-foot diameter rod. This yields a maximum solar torque of about 2×10-6 ft-lb (0.0075 AΩ2 for A = 2000 slug-ft2) around the center of mass of the system, whieh lies about 4 feet from the center of the satellite proper. This low torque is the result of a partial balance between the resultant force on the satellite and the resultant force on the rod, both yielding torques of the order of 5×10-6ft-lb (0.019 AΩ2). Even using this figure, the deflection due to solar pressure will be no larger than that due to magnetic torque. Thus such a satellite need not be especially designed to balance out solar torques.
In Ref. 8 the bending of an extensible rod of the STEM type, due to differential solar heating, was analyzed. Further unpublished work by P. Hrycak and by J. G. Engstrom at Bell Telephone Laboratories leads to the formula
d/L = (L/4r)π α TB / [κ + 4 + 16β/3]
for the deflection d of a rod of length L, radius r, and expansion coefficient α, where
T0 = (a0 S / π e0 σ)¼ κ = π k h T0 / r2 a0 S β = ei / e0
with a0 and e0 the rod external absorptivity and emissivity, ei the rod internal emissivity, σ the Stefan-Boltzmann eonstant, S the flux of solar radiation through unit area in unit time, k the rod thermal conductivity, and h the rod wall thickness. The dimensionless quantity gives the ratio between heat transferred by conduction and by radiation. Typical values of the above quantities are:
L = 60 ft r = 0.02 ft h = 2×10-4ft S = 442 Btu/ft2-hr σ = 171×10-11 Btu/ft2-hr-(°R)4 k = 65 Btu/ft-hr-°F α = 10-5/°F, a0 = 0.67 e0 = ei = 0.33,
the last five values being appropriate for beryllium copper. In this case T0 = 635°R, k = 219, β = 1, and d/L = 0.0656. Note that in this case κ ≫ 1, so that conduction effects dominate. Unless β ≫ 1, i.e., the outside of the rod is highly reflecting and the inside black, we may use the simpler formula
d/L = (L/4r) (α r2 a0 S / k h)
obtained by neglecting the remaining terms in the denominator of the previous formula in comparison with κ. In the present case this yields d/L = 0.0684, compared with the more exact value of 0.0656.
The above displacement d is the displacement of the tip mass at the end of the rod and hence produces a corresponding rotation of the yaw principal axis of inertia through an angle of order d/L = 0.0656 = 3.7° and an antenna pointing error of the same size. Note that this angle inereases linearly with L, so that thermal bending sets an upper limit on the length of a given type of rod which may be used for inertia augmentation. In observations of the Applied Physics Laboratory 1963 22A satellite [12] thermal bending manifested itself apparently as a high-frequency oscillation of the satellite's attitude attributed to the rapid heating of an extensible rod on passage from shadow into sunlight.
We shall discuss various large-angle motions of the gravity-oriented, gyro-stabilized satellite in order of their occurence. First we consider the injection, despin, and capture of the satellite in orbit and the equilibrium positions into which it may settle. Next we discuss the use of the gyros to flip the satellite in case it settles into the inverted equilibrium position. Finally we report the results of computer studies of various large motions.
We assume that the satellite is injected into a nearly circular, 6000-nm orbit with an initial spin rate of less than 1 rpm around an arbitrary axis. After injection, erection of a single 60-foot rod with a 20-pound tip mass then increases the moment of inertia around axes normal to the rod from 20 slug-ft2 to 2000 slug-ft2 and decreases the spin rate around these axes by a corresponding factor, e.g., from 250 rpo (revolutions per orbit) to 2 rpo. The component of spin around the rod axis is, of course, unaffccted by rod erection. This component of spin is removed by uncaging the gyros from their nominal equilibrium position, in which they have a zero net component of angular momentum around the rod axis (the body yaw axis) and allowing them to precess toward the spin rate vector. The net change of yaw angular momentum due to this precession is of the order of the gyro angular momentum H = 2AΩ, where A is the final moment of inertia (2000 slug-ft2) around the body pitch and roll axes normal to the rod, and the angular momentum due to the initial spin around the rod axis is 250 CΩ = 250 A0Ω = 2.5 AΩ, of the same order of magnitude. Note that this latter despin is in proportion to the difference of angular momenta, rather than their ratio, so that we might expect difficulties with the small differences of large numbers, leaving us with a sizeable angular velocity around the body yaw axis. However, the yaw component of angular velocity rapidly settles out; it is the yaw angular momentum, rather than the yaw angular velocity, which is of importance.
This is shown in Figure 17 and Figure 18 where, for the design of Section 2.3.2, yaw, pitch, and roll rate, obtained by a digital computer, are plotted against time. At t=0, the satellite was assumed injected into the desired orientation, with gyros at the null position, and with a yaw rate of 250 rpo (approximately ⅔ rpm). The elapsed time of Fig. 17 is two minutes, corresponding to boom erection. In this time, the yaw rate decreases to 20 rpo, while pitch and roll first peak at -12 rpo and -6 rpo respectively and then decay to -2.5 rpo at the end of boom erection. Subsequently, as shown in Fig. 18, all three rates decrease to less than 1 rpo at the end of 1 orbit.
Four equilibrium positions, in which the satellite is stationary with respect to the rotating local vertical, may be found by inspection. Two of these, shown schematically in Fig. 19, are the stable roll-vee positions with the gyro angular momentum vectors making a symmetrical vee with the orbit pitch axis (normal to the orbit plane), the gyro gimbal axes along the orbit roll axis (tangent to the orbit track), and the rod along the local vertical. The satellite antenna in this case is either directed toward the earth or away from the earth. We discuss the inversion of the satellite from the latter position in the sequel. Two other equilibrium positions are yaw-vee positions, Fig. 20, with tho gimbal axes along the orbit yaw axis, i.e., the local vertical, and the gyro angular momentum vectors again making a symmetrical vee with the orbit pitch axis, and the rod along the orbit roll axis. These two positions are unstable, however, just as they would be without the gyros.
Other equilibrium positions occur because of the presence of the gyro gimbal stops. Suppose, for example, that the satellite is rotated around the local vertical through 180° from its normal operating position. The gyro gimbal torquers which normally hold the gyro vee open against the 1 rpo steady precession of the satellite in orbit, now act with the precession to force the gyro gimbals against stops located near the body yaw axis. The resulting symmetrical reverse vee configuration (see Fig. 21a) is a possible satellite equilibrium position. Although the satellite antenna is still directed toward the earth in this position, it is an undesirable equilibrium position, because, when the gyro gimbals are against stops, their damping capability is severely reduced. This reversed equilibrium position can be made unstable by moving the gimbal stops in from the yaw axis and by choosing the gyro angular momentum to have a suitable value, as discussed in Section IV. The inversion of the position shown in Fig. 21 (a) and the corresponding reversed yaw-vee positions are unstable as before.
Finally we note the possibility of skewed equilibrium positions, in which the body principal axes do not coincide with the orbit axes and both gyro gimbals are against stops (see Fig. 21b). Examples are discussed in Section 4.3. Such unsymmetrical equilibrium positions may be easily eliminated by appropriate choices of stop positions and gyro angular momentum, but their occurrence suggests the necessity of a thorough investigation of equilibrium positions for any attitude control system, especially one in which constraints due to stops are present. The investigation of equilibrium positions also may serve as a guide in singling out lightly damped modes of large motion.
As we have already noted, the satellite may be captured, after injection into orbit, in inverted position with its antenna pointing away from the earth. With sufficiently large gyros it may be flipped from this position by changing the net gyro angular momentum by means of a simple signal injected by ground command into the gyro gimbal torquers. We simply reverse the polarity of the bias signals into the gyro torquers for a preset short time interval. The resulting change in angular momentum is just enough to cause the satellite to tumble, so that it is captured again in its normal operating position. Fig. 22 shows the result of such an inversion procedure, where the polarity is switched for ½ orbit. Here, we have plotted the cosine of the error angle, i.e., the angle between the body yaw axis and the local vertical.
Since only a limited amount of energy may be imparted by initial displacement of the satellite, computer studies were directed at the effects of high initial angular velocities. In Figs. 23-26 are shown some sample results of computer runs for the response of the satellite design of Section 2.3.2 to high initial rates, applied to the satellite in the stable roll-vee orientation. These may be regarded as responses to micrometeorite impacts, or as representative of initial transients following inadequate despin.
To save space, we again plot, as a function of time, only the cosine of the error angle between the yaw axis and the local vertical. However, the orientation of the satellite and of the gyro spin vectors are shown every half orbit in computer-made perspective drawings of a rectangular parallelepiped representing the satellite. The view is along the orbital track in the rotating, earth-pointing reference frame, so that the local vertical and normal to the orbital plane are in the plane of the paper. Plus signs are placed on the faces of the parallelepiped to avoid optical illusions. The gyro stops are indicated by dots. The reader may find more details about these drawings, as well as a description of computer-made movies showing large motions of the two-gyro satellite, in Ref. 15.
As is seen from Figs. 23-26, rates of the order of 4 rpo about pitch and roll damp out in about 10 orbits, whereas yaw rates of even 100 rpo settle out in about 5 orbits. In roll and yaw, the settling time and motion are similar if negative rather than positive rates are applied. The responses to positive and negative pitch rates are, however, different in character. A high positive pitch rate collapses the gyros toward the pitch axes, and a slowly decaying, essentially single-axis spin ensues. A high negative pitch rate opens up the gyros and drives them against the yaw stops. This sends the satellite into a complicated tumbling which eventually settles out.
If a micrometeorite of linear momentum m strikes the satellite at a lever arm L from a principal axis with moment of inertia I, the angular velocity ω imparted around that axis will be ω = mL/I. This velocity varies directly with L and inversely with I. For the design of Section 2.3.2, the yaw and pitch or roll lever arms are in the ratio 2/60, while the inertias are in the ratio 2000/20. A micrometeorite which imparts a pitch or roll rate of 4Ω will impart a yaw rate of (2/60) (2000/20) 4Ω = 13.3Ω. Therefore we see from Figs. 23-26 that the two-gyro spindle satellite is considerably more resistive to micrometeorite impact about yaw than about pitch and roll.
It is well known that gravity-gradient satellites will tumble if placed in a sufficiently eccentric orbit. Computer experiments showed that for the design of Section 2.3.2 this occurred at an eccentricity of about 0.2. Computer results for ε = 0.22 are shown in Fig. 27.
To settle the vexing questions of nomenclature and sign convention once and for all, we commence with a brief description of the quantities characterizing gravity-oriented satellite moving in a circular orbit (including the effect of small eccentricity later on) at the orbital angular velocity Ω. A convenient number to remember, to fix the magnitude of Ω, is that an orbit at an altitude of about 6800 statute miles coresponds to an orbital rate Ω = 1 radian/hr and an orbital period of 2π or about 6 hours, 15 minutes.
For our purposes the satellite is described by its principal moments A≥B≥C about principal axes x',y',z', along which the principal unit vectors i',j',k' lie. These principal axes form the body pitch, roll, and yaw axes, respectively.
When the satellite is in a position of stable equilibrium, the x',y',z' body axes coincide with orbit pitch, roll, and yaw axes x,y,z as in Fig. 19, with corresponding unit vectors i,j,k, normal to the orbital plane, along the orbit track, and along the local vertical toward the center of the earth. These orbit reference axes rotate at the orbital rate Ω (1 rpo) about the orbit pitch axis. It should be noted that, although a spindle-shaped body, formed by the extension of a single rod and tip mass, is shown, for which B ≈ A and C ≪ B, the small-angle analysis which follows covers the whole range of inertias, given by the inequalities
A > B > C
required for stability, and
B + C > A
imposed by rigid-body geometry
For small perturbations from equilibrium the satellite orientation is specified by the small pitch, roll, and yaw angles φx, φy,φz through which the body axes x',y',z' arc rotated from the orbit axes x,y,z, as in Fig. 28.
The corresponding satellite angular velocity vector ωs, with respect to inertial space, is given by
ωs = i(Ω + φx) + j φy + k φz
with respect to orbit axes, or
ωs = i'(Ω + φx) + j'(-Ω φz + φy) + k'(Ωφy +φz)
with respect to body axes.
By neglecting second-order terms in the dynamical equations of Section IV, as indicated in Section 4.2.2, we obtain the satellite equations of motion
The mathematics in the original document is only just legible and there is likely to be some mistakes in the following mathematical text.
Here Hc = H cos α and Hg = H sin α. This is an eighth-order linear system of equations for φx, φy, φz, ψg, φg, which splits immediately into a cubic pitch system for φx and ψg and a quintic roll-yaw system for φy, φz, φg, since the pitch motion depends only on the out-of-phase, or scissoring motion of the gyro gimbals, given by the difference angle ψg, and the roll-yaw motion depends only on the in-phase gimbal motion, given by the sum angle φg.
These equations can be reduced to dimensionless form by setting
p = (1/Ω)d/dt b = B/A, c = C/A h = (H/AΩ) cos α h' = (H/CD cos α κ = 1 + [K/(HΩ cos α)]
yielding the two sets of equations:
Pitch: [p2 + 3(b-c)]φx + (2h tan α)pψg = 0 - (h' tan α)pφx + (p+κh')ψg = 0 Roll-Yaw: [bp2 + 4(1-c) + 2h]φy + (1-b-c+2h)pφz + 2hφg = 0 -(1-b-c+2h)pφy + [cp2 + (1-b) +2h]φz - 2hpφg = 0 h'φy + h'pφz + (p+κh')φg = 0
If we insert appropriate terms on the right-hand sides of these equations to include the effect of given initial conditions and perturbing torques, we may regard the above systems as the Laplace transforms of the original set of differential equations, with transform variable p. The solution is then found by solving this set of linear, algebraic equations for φx, ψg, etc now interpreted as Laplace transforms, and calculating the residues at the poles of these functions of p. The transient response is entirely determined by the location of these poles and by the specific initial conditions. The steady-state response to a periodic perturbing torque at frequency NΩ may be determined by inserting constant right-hand sides, in general complex, setting p = iN, and solving for the amplitudes |φx|, | ψg| etc.
For given initial conditions and zero perturbing torques, the transforms are rational functions of p, with the characteristic pitch and roll-yaw polynomials as denominators, given by
Pitch Cubic:
f3(p) = p3 + c1p2 + c2p + c3 .............(2)
Roll-Yaw Quintic:
f5(p) = p5 + a1p4 + a2p3 + a3p2 + a4p + a5 ..........(3)
where
For stability, it is necessary that all the roots of the above polynomials have negative real parts. In particular, the magnitude of the real part of the root nearest the imaginary axis determines the rate of decay of the most lightly damped mode of motion. If this real part is -D, we can define a 1/e asymptotic settling time Ts = 1/2πD (in orbit periods) and use Ts as a measure of transient response, particularly suited for use with a digital computer. In Section 3.6 we discuss the determination of D as a function of b, c, α, h, h', and κ. Oncce it is reduced to a suitably small value by some choice of system parameters, the short-time transient response can be determined by solution of the differential equations with specific initial conditions and the system parameters readjusted, if necessary. Actually systems chosen on the basis of minimum asymptotic settling time seem to have quite adequate short-time, as well as steady-state, response.
The steady-state response to periodic perturbing torques at frequency NΩ, determined as previously outlined, is given in various cases by the following relations:
These amplitudes, independent of the gyro parameters, limit the minimum permissible satellite inertias for given perturbing torques.
Finally, an elliptic orbit of small eccentricity induces forced pitch vibrations at the orbital frequeney Ω. with amplitude
By straightforward differentiation, it is easily shown that the eccentrieity response |φx|ε has a single maximum as a function of the gyro opening angle α. For a larger than the value at which the maximum is attained, the eccentricity response decreases monotonically, approaching 2ε as α approaches 90°.
As mentioned in Section 3.3, a convenient single measure of transient response is the parameter D, the distance from the imaginary axis of the right-most root of the characteristic equations. One would like to know D as a function of the system parameters, D=D(b,c,a,h,h',κ). In general, this function is impossible to determine analytically and must be computed numerically. In order to limit such computations to ranges of the parameters b,c,a.,h,h',κ that give reasonable values of D, it is convenient to have bounds on D.
One set of bounds is given by the following theorem [16]:
We note that in both f3(p) and f5(p), in (2) and (3), q0= 1.Hence by the above theorem with l=0, if the system parameters are such that any of the coefficients c1, c2, c3, a1, ...a5, is small, then D will be small.
Likewise, if any coefficient in f3(p) or f5(p) is large compared to a subsequent coefficient, then the theorem tells us that D will again be small.
We note also that b- c < 1, so that the theorem applied to c2 gives, in pitch,
D2 ≤ b-c ≤1
i.e., the asymptotic settling time Ts, in pitch for the roll-vee system is bounded by Ts ≥ 1/2π = 0.159 orbit. (This is slightly larger than the corresponding bound Ts ≥ 5¼/2 √3π = 0.137 orbit obtained in Ref. 17 for a two-body satellite.)
From these bounds we conclude immediately that at best D can be of order unity, and to get a D of this order of magnitude the coefficients and ratios of coefficients in f3(p) and f5(p) must be at least of order unity.
Another useful bound is obtained by shifting the origin in the p plane to p = -D and applying the Routh criteria (see Section 3.7), to the shifted f3(p) polynomial. It is then found that one of the terms in the Routh array is
r = -2Dx2 + [8D2 + 3(b-c)]x - 8D3 - 6(b-c)D - 3(b-c)κh' where x = κh' + 2hh' tan2α
To have all roots to the left of the line Rep = - D, r must be positive. But it is easily verified that r > 0 only if
D ≤ ⅜ (b-c) / κh'
which gives an additional bound on D.
When the coefficients in f3(p) or f5(p) are either large or small, D can sometimes be expanded in a power series around a known root. This again restricts the parameter ranges over which D must be determined numerically. For example, suppose h' is small. The roots of f3(p) at h'= 0 are p=0, p = ±i√[3(b - c)]. However, it is well known [18] that each of the three branches of the triple-valued function p = p(h') is analytic in h'. Expanding around h'=0, say for p(0)=i√[3(b-c)], we have
p = i √[3(b-c)] + h' (dp/dh')h'=0 + ... p = i √[3(b-c)] - h'h tan2 α +...
with similar expressions easily obtained for the other two branches of p = p(h').
A particular case of interest is that of a spindle-shaped body. In this case, c→0,b→1, and the coefficients a11,...a5 of f5(p) all become large. One can then consider the equation ef5(p), in which the leading coefficient is small. However at c = 0, this equation is singular because it is reduced in degree from a quintic to a quartic. The quartic, with coefficients ca1...ca5 gives only four of the limiting roots as c→0. The fifth limiting root is however easily found by setting p = ρ/c, yielding f5*(ρ/c). Application of the expansion theorem to f5*(ρ) then yields the fifth limiting root p→-2hh'/c as c → 0.
This root gives a highly damped mode, and has a real part far in the left half-plane. The roots of interest in finding D are thus those of the limiting quartic:
f4(p) = p4 + (β1/h')p3 + β2p2 + (β3/h')p + β4 = 0 where β1 = 2h + 1 β2 = κ(2h + 1) + 4 - 2h β3 = 2h + 4 β4 = κ(2h + 4) - 2h
The limiting quartic f4(p) is a function of the three parameters κ,h,h', whereas, in this limiting case, the cubie f3(p) is a function of κ,h,h', and α. It turns out that D = D(h,h',κ,α) can be obtained graphically. Further, Dm, the maximum D for all possible h,h',κ,α can be found and has the value Dm = 0.317, attained at the values
0.77 < h < 0.78, h' = 1.29, κ = 0.92, α = 38°
(The value Dm corresponds to an asymptotic settling time
Ts = 1/2πDm = 0.502
orbits.) However, the description of the graphical technique and the derivation of D are too lengthy for inclusion here.
The over-all small-angle performance of the satellite attitude control system is characterized by its steady-state response to constant and periodic disturbances (solar torques, magnetic torques, orbital eccentricity) and by its transient response to sporadic disturbances (initial injection, micrometeorite impact). In proper design, one wants to diminish the response to all disturbances to below a suitably small level.
The steady disturbances have their main components at zero, orbital, and twice orbital frequency. As indicated earlier, their amplitudes may be diminished by inertia augmentation with extensible rods. Fortunately, it is easy to write down the formulas, (4)-(8), for satellite response to steady disturbances, and also easy to program these formulas for digital computation.
The computation of the transient response is not so straightforward, even in terms of the single measure D. An interesting theoretical problem is to find the maximum D as a function of all system parameters. Gradient or steepest descent methods, which first come to mind for the solution to this problem, seem to be difficult to apply, since the maximum D usually occurs in the neighborhood of multiple roots where the function D = D(b,c,h,h',κ, α) is singular.
However, although this is a theoretically interesting problem, its solution is not of great practical importance, as indicated in Section 2.1. It is more important to have a cheap method of computing D. A method that we have found useful involves the Routh criteria as follows:
Write the polynomial f(p) as
f(p) = a0pn + a1pn-1...+ ampn-m + ... + an = 0
and form the Routh array
a0, a2, a4, ...
b0, b2, b4, ...
c0, c2, c4, ...
d0, d2, d4, ...
where
b0 = a1, b2 = a3 ..., b2i = a2i+1
and
c2i = a2i+2 - (b2i+2a0 / b0)
d2i = b2i+2 - (c2i+2b0 / c0) etc i=0,1,2,...
Then the number of sign changes in the sequence a0, b0, c0, d0,etc (providing no term is zero) is the number of roots in the right-half plane. Because of its recurrence structure, this scheme is easily programmed on a digital computer.
To determine the real parts of the roots of f(p), one applies the scheme to a succession of translated half-planes as follows. If p = -D + ζ, then
f(-D+ζ)
= f(-D) + f'(-D)ζ + f0(-D)ζ2/2! +... + f(n)(-D)ζn/2! = 0
= q0ζn + q1ζn-1+ ... + qn = 0
where it is easily verified that
The Routh array applied to the coefficients q0, q1... qn then indicates the number of roots to the right of the line Re p = - D (Re ζ = 0). In order to locate the real parts of the roots to an arbitrary degree of accuracy, one applies this array on a sequence of nested intervals. For example, start with some large D = D* such that the Routh array applied on Re p = - D* indicates roots to its right. Take as the initial interval -D* < Re p < 0. In a stable system there will be roots between the right boundary (Re p = 0), and the left boundary (Re p = -D*). Next apply the Routh eriteria on Re p = - D*/2. There are two possibilites: (a) if there are no roots to the right of Re p = - D*/2, make this the new right boundary; the interval -D* < Re p < -D*/2 now has the same properties as the initial interval, (b) if there are roots to the right of Re p = -D*/2, make this line the left boundary of the new interval -D*/2 < Re p < 0, which again has the same properties as the initial interval. By applying this process n times, one ends up with an interval of width D*/2n, which contains roots but has no roots to its right. The accuracy of the loeation of the real parts of the roots closest to the imaginary axis can be set by prescribing the width of "the final interval. Since the widths of the successive intervals go down as 1/2n, the process converges rapidly.
After the real parts of the roots closest to the imaginary axis are found within some interval of desired width, say Interval 1, the same procedure can be used to find the next closest roots to the imaginary axis. One starts again at some sufficiently large D*, such that some roots fall to the right of Re p = -D* and to the left of Re p = - D L1, the left boundary of Interval 1. One makes these the left and right boundaries respectively, of an initial Interval 2, and applies the nested interval iteration again. The right boundary of Interval 2 in each iteration is characterized by having m roots to its right, where m is the number of roots contained in Interval 1.
The starting value D* can be chosen in various ways. If one is interested only in the roots closest to the imaginary axis, he can pick D* as
for then (see Section 3.5) there will be at least one root to the right of Re p = -D*. If it is desired to find the real parts of all the roots, D* can be chosen as
since it is well known [19] that this value of D* is a bound on the modulus of the maximum root and hence all the roots will be to the right of Re p= -D*.
We remark that this procedure may be easily extended to a method for finding both the real and imaginary parts of the roots of a real polynomial. It is only necessary to use well-known relations between the imaginary parts of the roots and certain members of the Routh array.
The above scheme goes rapidly on the IBM 7090 computer. For example, if the widths of Interval 1, Interval 2, ete. are set at 0.005, the running time is about 1000 cases a minute to find the real parts of all the roots of both the quintic, f5(p), and the cubic, f3(p). Tables calculated by this process were used in making the parameter survey whose results are summarized in Section 2.3.
The large-angle motion of the satellite is of course governed by non-linear differential equations, which in general must be integrated numerically. Nevertheless, a few analytical and intuitive insights are available. These are pointed out in the sections which follow.
We begin with a discussion of the pertinent dynamical and kinematic equations, including the effect of variable inertia, due to rod erection. Then we enumerate the equilibrium positions of the satellite, in which it is at rest with respect to the orbiting reference frame in a circular orbit, and show that certain restrictions must be placed on gyro angular momentum to eliminate undesired positions. This is followed by a discussion of satellite despin by the erection of a single rod and tip mass. Finally we show how the satellite may be inverted by ground command to the gyro torquers.
In the following, we make an explicit distinction between dynamical equations, valid in any coordinate system when written in proper vector form, and kinematic relations between various specific coordinate systems, This allows us to introduce a minimum number of different coordinate systems and to avoid a good deal of irrelevant algebraic complexity.
The rate of change of angular moment L about the satellite center of mass, with respect to a reference frame rotating at the satellite angular velocity Ωs, is governed by the equation
where M is the resultant torque around the center of mass, the sum of the gravity-gradient torque Ma, the total gyro precession torque MH, and the external disturbing torque ME. For a rigid body
where I is the inertia dyadic, given in terms of the principal moments of inertia A > B > C and corresponding principal vectors i',j',k', by
If ω is the satellite angular velocity relative to orbit reference axes,
where i is a unit vector normal to the orbit plane and ψ(t) is the polar angle of the satellite center of mass, measured from orbit perigee in earth-centered coordinates and satisfying the orbit equation
where ε is the orbit eccentricity and Ω = 2π/T0 being the orbit period. The gravity-gradient torque MG is given by
where k is a unit vector directed along the local vertical toward the center of the earth. Here and in the following, we consider only what Beletski [20] calls the restricted problem for which the motion of the center of mass is given by (13) and is unaffected by the motion around the center of mass. Finally the resultant gyro torque for a two-gyro, roll-vee configuration is
where Hi's are the gyro angular momenta, of fixed magnitude H, and the ωgi's are the gyro gimbal angular velocities. In terms of the gimbal angles φgi and the nominal roll-vee half-opening angle α, we have
The set of dynamical equations is completed by the gimbal equations of motion. If the gyro gimbals are not against stops, these are
where CD is the gyro damping constant, K the gimbal spring constant, and the constant bias torques Mgi are given by Mg2 = -Mg1 = HΩsin α. When tho gimbals are against stops, the reaction torques from the stops on the gimbals must be added to (19).
The orientation of the satellite body axes x',y',z', or the corresponding unit vectors i',j',k', with respect to the orbit axes x,y,z or corresponding unit vectors i,j,k (see Section 3.1), may be specifed in a number of ways. In classical dynamics, Euler angles have been traditionally used. They specify a rigid body's orientation with a minimum set of three numbers, and, in some of tho soluble problems of rigid body dynamics, lead to straightforward analytical manipulations.
From a computing point of view, Euler angles, however, have three serious disadvantages: (1) they involve trigonometric functions, which are expensive to compute, (2) they are singular when the rotation angle is zero, and (3) they are difficult to use in the visualization of complicated motions. We have chosen to use the so-called Euler parameters, rather than the Euler angles. A set of variables, perhaps even more suitable for the matrix algebra typical of modern computer programming, might be the direction cosines α, β, γ, etc., satisfying the relations
i = i'α + j'α' + k'α" j = i'β + j'β' + k'β" (20) k = i'γ + j'γ' + k'γ"
These α's should not be confused with the nominal vee opening angle, which we shall distinguish from the direction cosines, whenever they are used together, with a subscript. By using identities of the form
ḱ + ω × k = 0 (21)
satisfied by i, j, and k, we nay obtain 9 equations giving the rates of change of the direction cosines in terms of the direction cosines and the components of ω. We would then have 15 equations, including 3 satellite equations of angular motion, 1 equation of motion for the satellite's mass center, 2 gimbal equations of motion, and 9 equations for the direction cosine rates, yielding the 3 components of ω, the satellite polar angle ψ, the 2 gimbal angles, and the 9 direction cosines. The identities
α2 + β2 + γ2 = α2 + α'2 + α"2 = 1 etc αβ + α'β' + α2β2 = αα' + ββ' + γγ' = 0 etc
which must be satisfied initially, would then serve as checks on the numerical solution. Incidentally, it should be noted that the cosine of γ", the antenna pointing error angle is given by the direction cosine, between the local vertical and the body yaw axis.
We shall use the direction cosines to study equilibrium positions, but Euler parameters in the study of general satellite motion, since they are simply related to the deflection angles for small motion. If we assume that the (x',y',z' )-axes are formed by rotation of the (x,y,z)-axes through the angle θ around an axis with direction cosines mx, my, mz, the Euler parameters ψx, ψy, ψz, are defined by the relations
(ζx, ζy, ζz) = (mx, my, mz) sin(θ/2), t = cos(θ/2)
we now have [21]
i = i'( ξx2-ξy2-ξx2+χ2 ) + 2j'( ξx ξy - χξz) +2k'(ξxξz + χξz) j = 2i'( ξxξy+χξz )+j'( -ξx2+ξy2-ξz2+χ2 ) + 2k'( ξyξz - χξx ) (22) k = 2i'( ξxξy - χξy ) + 2j'( ξyξz + χξz) + k' ( - ξx2 - ξy2 + ξz2 + χ2 )
giving the direction cosines, and
completing a set of 10 equations for the three components of angular velocity ω, two ginbal angles, for Euler parameters, and the polar angle, ψ. Just as in the case or the direction cosines, the single identity
ξx2 + ξy2 + ξz2 + χ2 = 1
serves as a check on the numerical solution, Also note that, for small rotation θ ,
where φx, φy, φz are the small pitch, roll, and yaw angles. If these relations are inserted into the dynamical equations and second-order terms neglected, the linear equations for the small motion, (1), are obtained.
In coding the differential equations for the digital computer, it was found convenient to define cross-product and dot-product subroutines:
A × B = (-A3B2 + A2B3, A3B1 - A1B3, -A2B1 + A1B2) A.B = A1b1 + A2B2 + A3B3
This allowed the coding to follow closely the vector form of (9)-(11), which was useful from the standpoint of both coding simplicity and debugging.
For non-dumbbell satellite shapes, say b = 0.9, c = 0.5, the five-point predictor-corrector, with ⅔ rule as given by Hamming (Ref. 22, Chapter 15) was used. However, in the spindle case, b = 1.0, c = 0.01, the differential equations become singular because the small number c multiplies a derivative, and the five-point, ⅔ rule scheme was found to be very slow. Following a suggestion of R. W. Hamming, a simple three-point predictor-corrector scheme (Ref. 22, p. 186) was then tried. It turned out to be three to four times faster than the five-point scheme and to give about the same accuracy.
In the computer runs, the gyro stops were simulated by hardening springs. For example, for φg1 > β, the normal spring restoring torque of Kφg1, was replaced by
Kφg1 + B (φg1 - β) + C / (θ - φg1)
where B, C, β, and θ are constants, to simulate the pitch stop of the first gyro. The same expression but with different constants was used for the yaw stop. Specifically, in all the computer runs for the spindle-shaped satellite reported here, K = 0 and the pitch-stop values for Gyro 1 were
B = 59AΩ2, C = 0.01AΩ2, β = 58°, θ = 60°
For the yaw stop, {β and θ above were replaced by β = -20° and θ = -30°. Corresponding, symmetrical stop constants were used for Gyro 2.
For a circular orbit we can easily obtain a useful expression for the rate of change of kinetic and potential energy, relative to orbit axes. We take the scalar product of the satellite equation of motion (9) with the relative angular velocity ω and combine it with the two gyro equations, (11), multiplied by φgi. After some routine algebra, we obtain the relation
This expression is useful in the estimation of various quantities, in partioular the velocity required to tumble the satellite, and conditions necessary for capture [23].
When the satellite moves in a circular orbit in such a way as to remain stationary with respect to the local vertical, its motion satisfies the equilibrium equation
as defined by the table of direction cosines, (20). Thus we have the following general result:
I. The equilibrium positions of any symmetrical, gravity-oriented, gyro-stabilized body in a circular orbit must be such that the principal axis of least inertia and the resultant angular momentum are perpendicular either to the orbit roll axis (β" = Hy= 0) or to the orbit yaw axis (γ" = Hz= 0).
In the case of a roll gyro system, with all gimbal axes parallel to the body roll axis, the resultant gyro angular momentum must have the form
Thus, β"=Hy=0 implies βHz=0 and γ"=Hz implies γHx'=0. If we now assume that the motion of the gyro gimbals is restricted by stops along the body yaw axis, so that Hx'>0, an assumption appropriate for the case of the two-gyro roll-vee, we have the following result:
II. The equilibrium positions of any symmetrical, gravity-oriented, roll-gyro-stabilized body in a circular orbit must be such that the body roll axis is either parallel to the orbit roll axis (β'= ±1) or parallel to the orbit yaw axis (γ'= ±1).
So far we have not made use of the gyro gimbal equilibrium equations. If the gimbals are not agaiust stops, from (19) for the roll-vee, these take the form
where we now denote the nominal gimbal angle by α0 to avoid confusion with the direction cosines. If the flex-lead constraint is negligible, i.e., K = 0, these two equations imply that γ'Hy-β'Hz=0 which in either case in II (β'=±1,Hy=0 or γ'=±1,Hz=0) implies that Hy= Hz= 0, so that no torque is exerted on the satellite by the gyros. We then have:
III. The equilibrium positions of a symmetrical, gravity-oriented, free roll-vee-gyro-stabilized body must be sutch that the resultant gyro torque vanishes (i × H = 0) and either the body pitch, roll, and yaw axes are parallel to the orbit pitch, roll, and yaw axes (i'=±i, j = ±j, k'±k) or the body pitch, roll, and yaw axes are parallel to the orbit pitch, yaw, and roll axes, respectively (i'= ±i, j'=±k, k'=±j).
The signs of course must be chosen so that the above represents a proper rotation. Note that the above applies to any roll gyro system for which the resultant torque around the body roll axis exerted by the gyros on the satellite vanishes. If i'=+i, the second set of equilibrium positions gives Fig. 20, with the gyro gimbal axes along the orbit yaw axis in a yaw-vee configuration. Small pitching motion around these equilibrium positions is governed by a characteristic equation of the same form as that for the roll-vee, ( 1), except that the coefficient 3(B-C)/A>0 is replaced by (C-B)/A<0). Thus these equilibrium positions are unstable.
The equilibrium position i'=i,j'=j,k'=k of the first set is shown in Fig.19a. It corresponds to the normal operating position with the body yaw axis, on which the antenna is situated, directed toward the earth. The inverted position (see Fig. 19b) i'=i, j'=-j, k'=-k is also stable, since it merely corresponds to an interchange of the two gyros. This bistability is characteristic of gravity-oriented bodies and a gravity-oriented communications satellite must either use two antennas, with associated switching, or incorporate some means of flipping the satellite in response to ground command. The latter possibility is discussed in some detail in the sequel.
Tbe reversed roll-vee equilibrium positions Fig. 21(a), with i' = -i, remain to be investigated. The corresponding yaw-vee positions are still unstable. If the gyro gimbals were completely free, satellite precession and the bias torques, now acting together, would rotate the gyro gimbals from the reversed roll-vee until they formed a normal roll-vee around the orbit pitch axis. But with gimbal stops, making the angles ±α* with the body pitch axis, the gimbals rotate until the stop reaction torque and the bias torque sustain the 1 rpo steady precession of the satellite in orbit. The stability of this reversed roll-vee position can be investigated by using the characteristic equations for the normal roll-vee, with H cos α0 replaced by-H cos α* and a large spring constant K*, introduced to take the stop compliance into account. In particular the coefficient
a4 = [ (A-B)Ω + 2H cos α0] [4(A-C)Ω + 2H cos α0] / BC Ω2
in the roll-yaw characteristic equation, is replaced by
a4* = [ (A-B)Ω - 2H cos α*] [4(A-C)Ω - 2H cos α*] / BC Ω2
If
(A-B)Ω < 2H cos α* < 4(A-C)Ω
this is negative und the equilibrium position is unstable.
The instability of the reversed roll-vee when α* satisfies the above inequalities is shown in Fig. 29. In this case, the system parameters are the same as those of the sample design in Section 2.3.2 with α*=80°. Initially the gyros are against the stops and the satellite has rates of 0.05 rpo about all three axes. It is seen that the satellite turns around the yaw axis and settles down to rest in the desired orientation in less than five orbits.
When the gyro gimbals are against stops, the gyros exert a torque on the body and in general there are other, skewed equilibrium positions. To investigate these positions without getting involved in the details of stop compliances, etc., which depend on the specific gyros used, we consider only two idealized cases, the first with stops along the positive and negative body yaw axes but with no stops along the body-pitch axis, and the second with stops along the pitch axis as well as along the yaw axis.
In both cases the gyro spin axes may be back-to-back along the yaw axis, but this is a case of zero net gimbal torque already treated and is easily eliminated by moving the gimbal stops in slightly. In the first case both spin axes may lie along the body yaw axis against stops, so that
H1 = H2=k'H, H = 2k'H
From (24) and (25) we again have two cases to consider: (a) β" = Hy = 0, and (b) γ" = Hz = 0. In case (a) we have
Hz = 2γ"H = 4(A-C)Ωα"γ"
The subcase β" = Hy = 0 ; γ"=0 is easily shown to be unstable, so there remains only the position given by
α" = H/[2(A-C)Ω]
Unless H > 2(A - C)Ω, this yields an equilibrium position which can be maintained by the stop reaction torques. These torques are of course one-sided, since the stop can only push and not pull. This undesirable skew equilibrium position can be eliminated by making
H > 2(A-C)Ω
A similar position for case (b) (A similar position for case (b) γ" = Hz=0 can be eliminated by satisfying the less restrictive condition H > (A - C)Ω/2.
The corresponding situation with yaw and pitch stops finds one gyro against the yaw stop and the other against the pitch stop (see Fig 21b). For example, suppose
H1=k'H, H2=i'H, H=H(i'+k')
Now in case (a), β" = Hy = 0, we have
Hz=H(γ+γ") = 4(A-C)Ωα"γ"
or, in terms of the angle θ between the x and x' axes,
sin 2θ = -[H/2½(A-C)Ω] sin[θ+π/4)]
The two roots of this equation in the interval -π/4 < θ < π are excluded, because they require stops which pull on the gimbals. On the other hand, the two roots in the interval π < 0 < 3θ/2 yield possible equilibrium positions. These roots exist only if H <2½(A - C)Ω. Again the case (b) γ" = Hz = 0, yields no equilibrium positions of this type under the less restrietive condition H > 2-½(A - C)Ω. Since an increase in gyro angular momenturn tends to degrade the transient performance of the system, we shall assume in the following that the gyro gimbals arc limited in excursion by both yaw and pitch stops, so that only the restriction H > 2½(A - C)Ω need to be satisfied.
In the case of an unsymmetrie satellite, a similar but more complicated analysis of the equilibrium positions can be earried out.
We have already indicated the necessity of augmenting the satellite inertia to increase the gravity-gradient restoring torques to required levels. If this inertia augmentation is done after injection into orbit, it also reduces the satellite angular velocity to a level where the gravity-gradient torques may become effective in aligning the satellite with the local vertical. One method of inertia augmentation is the extension of so-called STEM rods deseribed in detail in Ref. 8. These are beryllium copper tapes which form straight, tubular rods when unwound from a drum. If they are used together with dense tip masses, the satellite inertia may be increased by several hundred-fold without a proportional increase in solar torque. In the following sections we first consider the effect a variable inertia has on the general form of the satellite equation of motion and then discuss satellite despin using a single extensible rod in combination with two gyros.
We may derive all of the dynamical equations for the motion of a gravity-oriented body by integration of the general equations of motion for a continuous medium. In fact this is perhaps the most direct way of calculating the gravity-gradient torque, which is due to the variable gravitational body force acting on each mass element of the body. The resulting equation of motion, (9), applies to rigid and flexible bodies alike, provided that the angular momentum L is calculated correctly. L is given in general by the integral
L = ∫B r × (∂r/∂t + ωs × r)dm (26)
where r is the radius vector from the center of mass of the body B to the mass element dm. For a rigid body, r differs from its initial value rg only by a rotation and ∂r/∂t = 0, yielding the usual form
L = ∫B r × (ωs × r)dm = I.ws
but in general r depends both on r0 and t, so that
L = I.ωs + ∫B r × (∂r/∂t)dm
where the inertia dyadic I depends on t.
Let us now consider the extension of a single massless rod with tip mass ma . If a(t) is the radius vector from the center of mass of the satellite proper to the tip mass, (26) yields
for satellite mass ms, and the inertia dyadic for the satellite around its center of mass
I0 = A0i'i' + B0j'j' + C0k'k'
When the rod is erected parallel to itself, as would normally be the case,
and the effect of rod extension is entirely taken into account by the time-dependent inertia I. If the rod is extended along the axis of minimum moment of inertia, a=k'a(t) and
I = (A0+ma2)i'i' + (B0+ ma2)j'j' + C0k'k'
When the moments of inertia of a torque-free spinning body are increased by a factor of N, conservation of angular momentum requires that the angular velocity of the body decrease by a factor of 1/N. On the other hand, if the spinning body contains a spinning rotor, an increase in the angular momentum of the rotor produces a corresponding decrease in the angular momentum of the body and hence a reduction in the angular velocity of the body. Both elements exist in the gravity-oriented, gyro-stabilized satellite. Inertia augmentation is required to obtain gravity-gradient torques of the proper level; rotation of the gyro gimbals provides a change in angular momentum.
A single extensible rod-tip mass combination provides adequate gravity-gradient torques, if erected along the satellite yaw axis. However, erection of such a rod reduces only pitch and roll injection angular velocities; the yaw component is unaffected. This may be removed by using the gyros as reaction wheels.
Suppose the gyros are caged at their null position during the rod erection phase. Then, neglecting gravity-gradient and external torques during the short erection time, we have a freely spinning body. If the initial components of angular momentum are all of the order of magnitude A0NΩ, where A0 is the moment of inertia about all three axes, we may reduce the pitch and roll components of angular velocity to order Ω (1 rpo), with respect to inertial space, by extending a single rod that increases the pitch and roll moments of inertia from A0 to A = NA0. The yaw angular velocity remains equal to NΩ.
The yaw angular momentum, A0NΩ, however, is of the order AΩ, the same order of magnitude as the angular momentum H of each gyro. The gyros then are large enough to absorb the residual angular momentum. If the gyro gimbals are now released, the spin axes will tend to line up with the residual angular yelocity around the yaw axis. One gyro spin axis rotates until constrained by the yaw axis gimbal stop; the other rotates until constrained by the pitch axis gimbal stop. There is thus a net change in yaw gyro angular momentum of order H and, furthermore, because of the rate-seeking property of the gyro, it always occurs in the correct sense to reduce the satellite angular momentum.
This qualitative argument has been supported by computer runs for given parameter values and initial rates. In particular one may determine the gyro size needed to despin from a given initial yaw rate. Actually it is not necessary to separate the erection and uncaging phases of injection, so long as the gyros are not uncaged before erection.
The result of such a computer run was given in Fig. 17 and Fig.18. In one orbit the angular rates are reduced to less than 1Ω and capture occurs.
As we have already noted, gravity-oriented bodies are bistable, i.e., they are in stable equilibrium with the axis of least inertia, on which the antenna would be mounted, both directed toward the earth and away from the earth. In this section we discuss a method of flipping the satellite by means of a simple ground command injected into the gyro gimbal torquers.
When the satellite is in either of the above stable equilibrium positions, its total angular momentum is AΩ + 2H cos α0 around the pitch axis. If we could somehow rotate the two gyro gimbals instantaneously, so that both spin axes pointed along the pitch axis, the total angular mometum would become A(Ω + ω) + 2H, where ω is the pitch angular velocity with respect to the orbit frame. Since the gimbal rotation is assumed instantaneous, the total angular momentum is conserved, i.e.,
AΩ + 2H cos α0 = A(Ω + ω) +2H or ω = -2H(1 - cos α0)/A
Thereafter, the single-axis, pitching motion is governed by an equation of the form
For a gyro angular momentum satisfying this condition we can excite a tumbling motion by collapsing the gyro spin axes toward the pitch axis. Actually, since we only want to rotate the satellite through a half revolution, the gyro angular momentum need barely exceed this minimum value. Furthermore, we may collapse the gyro spin axes toward the pitch axis by simply reversing the bias torques applied to the gimbals for a suitable length of time. For a spindle satellite with (B- C) / A = 0.99, computer runs (see Fig. 22) show that the satellite may be inverted by applying this reversed bias for about a half orbit. For a satellite with (B - C)/A = 0.4 it turns out that it is only necessary to remove the bias torques for a fraction of an orbit. For any satellite a suitable combination of bias torque and time can always be found to flip the satellite into its desired operating position, providing the relation (27) is satisfed.
We remark that bias torques could also be used to rotate the gyro gimbals against the yaw stops. A similar, single-axis argument then gives an expression like (27). However, with the gyros back-to-back against yaw stops, the satellite has negligible yaw stiffness, and is vulnerable to yaw disturbances. This possibility of inversion was therefore not pursued.
The authors would like to acknowledge the contributions of Dr. J.E. DeLisle, B. M. Hildebrant, and E.G. Ogletree, of the MIT Instrumentation Laboratory, not only in the initial study but also in subsequent consultations. They would also like to thank J. Sautter and Mrs. C. Kimme for assistance with computer programming.
1. Pierce, J R, Orbital Radio Relays, Jet Propulsion, 25, 1955, pp. 153-157.
2. Ogletree, E G, Sklar, S J, Mangan, J G, Satellite Attitude Control Study, Part I, Report No. R-308, Instrumentation Laboratory, MIT, July, 1961.
3. Ogletree, E G, Mangan, J G, Petrie, T D, Sklar, S J, Merkle, R F, Hutchinson, R C, Satellite Attitude Control Study, Part II, Report No. R-308, Instrumentation Laboratory, MIT, Feb., 1962.
4. DeLisle, JE, Ogletree, E G, Hildebrant, B M, Applications of Gyro Stabilizers to Satellite Attitude Control, Progress in Astronautics and Aeronautics, Vol. 13, Academic Press, 1964, pp. 149-175.
5. Burt, E G C, On the Attitude Control of Earth Satellites, Eighth Anglo-American Aeronautical Conference, London, Sept., 1961.
6. Scott, E D, Control Moment Gyro Gravity Stabilization, Progress in Astronautics and Aeronautics, Vol. 13, Academic Press, 1964, pp. 103-149.
7. Kamm, L J, Vertistat: An Improved Satellite Orientation Device, ARS Journal. 32, 1902, pp. 911-913.
8. Paul, B, West, J W, Yn, E Y, A Passive Gravitational Attitude Control System for Satellites, BSTJ, 42, Sept 1963, pp. 2195 2238.
9. Fletcher, H J, Rongved, L, Yu, E Y, Dynamics Analysis of a Two-Body Gravitationally Oriented Satellite, BSTJ, 42, Sept 1963, pp. 2239 2266.
10. Paul, B, Planar Librations of an Extensible Dumbbell Satellite, AIAA Journal, 1, Feb, 1963, pp. 411 418.
11. Newton, R R, Damping of a Gravitationally Stabilized Satellite, AIAA Journl, 2, Jan., 1064, pp. 20-25.
12. Fischell, R E, Mobley, F F, A System for Passive Gravity-Gradient Stabilization of Earth Satellites, Progress in Astronautics and Aeronautics, Vol 13, Academic Press, 1964, pp. 37-71.,
13. Lewis, J A, Viscous Damping of Gravitationally Stabilized Satellites, Proc. Fourth U.S. Nat. Congr, Appl. Mech., Berkeley, Calif., June, 1962,
14, Zajac, E E, Limits on the Damping of Two-Body Gravitationally Oriented Satellites, AIAA Journal, 1, Feb, 1963, pp. 498-499.
15. Zajac, E E, Computer-Made Perspective Movies as a Scientific and Communication Tool, CACM, 7, March, 1964, pp. 169-170.
16. Zajac, E E, Bounds on the Decay Rate of Damped Linear Systems, Quart. Appl. Maths, 20, Jan 1963, pp 383-384
17. Zajac, E E, Damping of a Gravitationnlly Oriented Two-Body Satellite, ARS Journal, 32, Dec, 1962, pp. 1871 -1875.
18. Knopp, K, Theory of Functions, Vol. II, Dover, N, Y., 1947, p. 121.
19. Uspensky, J V, Theory of Equations, McGraw-Hill, N, Y., 1948,p. 75.
20. Beletskii, V V, The Librations of a Satellite, in Artificial Earth Satellites, Vol. 3 (translated from the Russian), Plenum Press, N. Y,, 1961, pp, 18-45.
21. Whittaker, E T, Analytical Dynamics of Particles and Rigid Bodies, Cambridge Univ. Press, 1961, p. 8.
22. Hamming, R W, Numerical Methods for Scientists and Engineers, McGraw-Hill, 1962.
23. Zajac, E E, Capture Problem in Gravitational Attitude Control of Satellites, ARS Journal, 31, Oct 1961, pp. 1464-1466.
24. Morgan, S P, Yu, E Y, Expected Frequencies of Meteoritic Disturbances of Gravitationally Oriented Satellites, to appear.