Effects of Payload Motion on the Attitude Behavior of a Spinning Spacecraft

Table of Contents

List of Figures

Abstract

Chapter 1 - Introduction

Chapter 2 - Mathematical Model
2.1 Model Description
2.2 Equations of motion
2.2.1 Kinematics
2.2.2 Kinetics

Chapter 3 - Case Studies
3.1 Numerical Integration
3.2 Stationary Payload (Case 0)
3.3 Single Payload
3.3.1 One oscillating payload (Case 1)
3.3.2 Varying Payload Inertia (Case 2)
3.3.3 Varying the frequency of payload oscillation (Case 3)
3.3.4 Varying the amplitude of payload oscillation (Case 4)
3.4 Motion of Multiple Payloads
3.4.1 Two Oscillating Payloads (Case 5)
3.4.2 Four Oscillating Payloads (Case 6)
3.5 Varying payload location
3.5.1 Single Oscillating Payload (Case 7)
3.5.2 Two Oscillating Payloads (Case 8)

Chapter 4 - Summary of Results
4.1 Characteristics of Nutation Angle
4.2 Effect of Payload Oscillation Frequency
4.3 Effect of Payload Oscillation Amplitude
4.4 Effect of Payload Location
4.5 Effect of Payload Inertia

Chapter 5 -Conclusions
5.1 Summary and Conclusions
5.2 Recommendations for Future Work

References

Appendix

A. Autolev code for obtaining equations of motion

B. Matlab code for simulation of nutation angle

C. Data used for plots in chapter

List of Figures

Figure 1.1 - Dual Spin Spacecraft

Figure 2.1 - Spacecraft Model

Figure 3.1 - Nutation Angle for Case 0

Figure 3.2 - Code Verification for u3 (0) = 0.5

Figure 3.3 - Code Verification for u3 (0) = 1

Figure 3.4 - Nutation angle for one rod oscillating

Figure 3.5 - Effect of inertia on nutation angle

Figure 3.6 - Effect of payload oscillation frequency on nutation angle

Figure 3.7 - Motion of Spacecraft

Figure 3.8 - Effect of payload oscillation amplitude on nutation angle

Figure 3.9 - Recommendations for reducing nutation

Figure 3.10 - Nutation Angle for Two Oscillating Opposite Rods

Figure 3.11 - Nutation Angle for Two Oscillating Adjacent Rods

Figure 3.12 - Nutation Angle for All Rods Oscillating

Figure 3.13 - Nutation Angle for P = 0.25 m (One Oscillating Rod)

Figure 3.14 - Nutation Angle for P = 0.5 m (One Oscillating Rod)

Figure 3.15 - Nutation Angle for P = 0.75 m (One Oscillating Rod)

Figure 3.16 - Nutation Angle for P = 1 m (One Oscillating Rod)

Figure 3.17 - Nutation Angle for P = 0.25 m (Two Opposite Oscillating Rods)

Figure 3.18 - Nutation Angle for P = 0.5 m (Two Opposite Oscillating Rods)

Figure 3.19 - Nutation Angle for P = 0.75 m (Two Opposite Oscillating Rods)

Figure 3.20 - Nutation Angle for P = 1 m (Two Opposite Oscillating Rods)

Figure 4.1 - Effect of Payload Oscillation Frequency on Nutation Amplitude

Figure 4.2 - Effect of Payload Oscillation Frequency on Mean Nutation Angle

Figure 4.3 - Effect of Payload Oscillation Amplitude on Nutation Amplitude

Figure 4.4 - Effect of Payload Oscillation Amplitude on Mean Nutation Angle

Figure 4.5 - Effect of Payload Location on Nutation Amplitude

Figure 4.6 - Effect of Payload Location on Mean Nutation Angle

Figure 4.7 - Effect of Payload Inertia on Nutation Amplitude

Figure 4.8 - Effect of Payload Inertia on Mean Nutation Angle

Title - Effects of Payload Motion on the Attitude Behavior of a Spinning Spacecraft Name - Srihari Gowri Shankar

Abstract

Spacecraft systems almost always include subsystems such as imaging or scanning instruments that are designed to move relative to the main spacecraft body during some part of the spacecraft’s mission. Depending on the relative size and location of such moving parts, their motion can have appreciable effects on the overall translational and/or rotational motion of the spacecraft. This research project investigates the impact of oscillatory motion of the spacecraft-mounted payload on the global attitude behavior of the spacecraft. The investigation is carried out by first deriving the equations of attitude motion of a spacecraft system that consists of a main cylindrical body to which are connected four slender rods each of which has a prescribed oscillatory motion relative to the main body. The main body is given initial spin and transverse angular velocities, and the resulting nutation angle is monitored over time.

The results obtained show that oscillatory motion of a payload relative to the main spacecraft body will, in general attenuate the spacecraft’s nutation. The degree of nutation attenuation depends on several factors, including the inertia of the payload compared to that of the spacecraft, the location of the payload on the spacecraft, the frequency of the payload’s motion, and the amplitude of such motion. The implication of these results is that payload motion can be exploited to assist a spacecraft’s nutation damper system. In fact, any movable payload can save as a backup nutation damper if the designated nutation damper should fail.

Chapter 1 - Introduction

1.1 Spinning Spacecraft

Developments in the study of the dynamics of spinning rigid bodies have led to the design of spacecraft that spin continuously while in motion. There are two main types of spinning spacecraft in use today - single spinners and dual-spinners. A single spinner consists essentially of one main body that spins in its entirety while in motion. On the other hand, a dual-spinner is made up of two main bodies connected together by a one- degree-of-freedom hinge. One of the bodies spins like a single spinner would, while the other body remains de-spun. The advantage of a dual-spinner over a single spinner is that it has a de-spun section, which can be used for experiments that call for a stable platform. The spin motion of part, or all of a spinner, imparts “spin rigidity” or “gyroscopic stiffness” to the system, thereby allowing it to maintain its desired direction of motion, even in the presence of small disturbances.

In the absence of applied torques, a spinner has the tendency to spin about one of its central principal axes of inertia. Transient disturbance torques will cause the spin axis to cone about its original orientation, so that the spin axis makes an angle θ, often referred to as nutation angle, with the original heading of the spinner. For a rigid single spinner, such spin motion is stable if the spin axis is either an axis of maximum inertia or an axis of minimum inertia¹. The term “stability” as used in this context is understood to mean that the angle θ eventually decays to zero, so that the body returns to a state of pure spin in which there is no coning. If there is some form of energy dissipation within the system, spin about the axis of minimum inertia ceases to be stable. The axis of maximum inertia then becomes the only axis of stable spin¹. Now, every real system is likely to have some type of energy dissipation present. This may be in the form of structural damping, fuel slosh, or similar phenomena. Hence, the only realistic requirement for stability (condition for the decay of θ) for a spinning spacecraft is spin about the central principal axis of maximum inertia. In other words, the spacecraft should be oblate like a wheel rather than prolate like a pencil in its overall shape for attitude stability to be guaranteed.

For a dual-spin spacecraft, the conditions for attitude stability, as established by Mingori², Likins³, and others⁴, are slightly more complex. If the rate of energy dissipation in the rotating portion (rotor) of the system is greater than that in the non- rotating part (stator), then the system will spin stably about the axis of maximum inertia of the rotor. On the other hand, if the energy dissipation rate in the stator is larger than that in the rotor, then spin about the rotor’s axis of minimum inertia is stable.

1.2 Nutation Control

Normally, a spinner is designed to spin continuously about one of its principal axes while this spin axis moves tangentially to the trajectory of the mass center. The trajectory is usually a simple curve. As stated above, slight transient disturbance torques will impart coning motion on the spinning spacecraft. Such coning motion is undesirable because it affects the pointing accuracy of devices such as communication antennas that are mounted on the spacecraft and used for communication with ground-based stations. It is thus important that the coning motion or nutation of the spacecraft be reduced to a pure spin or to as close to a pure spin as possible within a reasonable time period. Hence, every spinning spacecraft is equipped with some form of nutation control system or nutation damper that is designed to eliminate or drastically reduce spacecraft nutation. Nutation dampers can be active or passive. An active nutation damper uses one or more onboard actuators to generate a torque that is applied to the spacecraft in such a way as to counteract the coning motion. Such actuators would have to be powered, and would therefore need fuel to operate. Moreover, active nutation control would normally involve a feedback control system that includes a sensing component for measuring the nutation angle and comparing it to the desired value. Any resulting error is used to determine the required corrective torque, which is then implemented by the onboard actuators. Since external disturbances are present during most of the duration of a spacecraft’s mission, the nutation determination/correction cycle can be expected to continue indefinitely. The active approach to nutation damping can thus become costly in its fuel requirements, and the fuel cost can become prohibitive for space missions of very long duration. Passive nutation dampers exploit the dynamics of the system and knowledge of the stability criteria for the system. For example, if a spinner is oblate in shape and is designed to spin about its axis of maximum inertia, a mass-spring-damper system can be included in the spacecraft system, and arranged in such a way that any coning motion of the spacecraft immediately triggers vibratory motion of the mass spring damper system, thereby causing energy dissipation within the system. Since the original spinner (without the mass-spring-dashpot system) is already passively stable, energy dissipation due to the mass-spring-dashpot system causes the overall system to be “even more stable” - any nutation caused by external disturbances would be reduced or damped out faster than it would have been without energy dissipation. In other words, the coning motion will decay with reduced time constant, and the spacecraft will quickly return to the state of pure spin. A passive approach to nutation damping is advantageous because it simply exploits the inherent dynamics of the system and requires no fuel. Thus, a well-designed passive nutation damper has the same life span as the spacecraft hardware and is not limited by onboard consumables. Because of these advantages, most spinners do in fact use passive, rather than active devices to damp out nutation.

Passive nutation control systems are based on the fact that the nutational stability of a spinning spacecraft can be improved by deliberately building into the system some energy dissipation device that is triggered into action by the coning or nutational motion of the spacecraft. The assumption of course is that the inertia distribution of the spacecraft system is such that energy dissipation does in fact cause nutation to be damped. Various types of energy dissipating mechanisms have been used for the purpose of nutation damping. These include [5, 6, 7] mass-spring-dashpot systems, compound pendulum with damper, and viscous fluid in ring-shaped tubes.

1.3 Payload Motion

Spacecraft systems sometimes include instruments or payload that are mounted on the main spacecraft and that are allowed to move in some way relative to the main spacecraft body. Examples include antennas and cameras, or imaging instruments in general. For example, Fig 1.1 shows a dual-spin spacecraft in which A is the rotor and B is the stator. Ordinarily, the orientation of B remains fixed in space while A spins relative to B. If it should become necessary to rotate B back and forth about A’s axis of rotation, such motion will not affect the global attitude motion of the system if A’s axis of rotation is also a central principal axis of inertia for B⁸. Otherwise, back and forth motion of B will affect the overall motion of the spacecraft system.

Figure 1.1 - Dual Spin Spacecraft

Abbildung in dieser Leseprobe nicht enthalten

Several studies [9, 10, 11] have shown that the overall attitude motion of a spinner can sometimes be modified by allowing a payload mounted on the spinner to oscillate relative to the spinner about an axis parallel to the spinner’s spin axis. A case study on the Galileo Spacecraft [12, 13] showed that moving such payload about an axis perpendicular to the spacecraft spin axis can also influence nutation. Payload motion can drive or attenuate spacecraft nutation. The question then becomes whether such payload motion can drive nutation to the point of overpowering in-built spacecraft nutation control system; or whether, on the other hand, there are circumstances under which such payload motion reduce nutation so efficiently that a separate nutation control system becomes unnecessary.

1.4 Project Description

The purpose of this research work is to investigate the effect of payload motion on the attitude behavior of spinning rigid spacecraft. As stated above, it is already a known fact that payload relative motion can influence the attitude motion of a spacecraft. However, the studies conducted so far relate to special spacecraft systems and to specific payload types. In some of the results, the effects observed can be attributed to imbalance due to mass center offset and / or existence of products of inertia⁹. The goal here is to focus on the case where the axis of motion of the payload is perpendicular to the spin axis of the spacecraft, and to monitor effects that are due solely to motion of the payload relative to the spacecraft bus, and not those effects that are attributable to static or dynamic imbalance.

Motivation for this work comes from the potential benefits that can accrue from the results. For example, it may turn out that by slightly modifying the placement of a payload and permitting it to move, the same payload can serve as a backup nutation damper for a spinner. It is also possible to increase the efficiency of an existing nutation damper system by slightly modifying the inertia properties and / or motion sequence of some payload on a spacecraft.

For this study, a spacecraft system containing movable payloads located at various points on the spacecraft body is assumed. Equations of motion for such a system are derived, and these equations are integrated numerically for specified relative motion of the payload. The main output from the study is the spacecraft nutation angle. The effect on the time history of the nutation angle of various payload motion scenarios, and various payload/spacecraft mass properties are studied. The ultimate goal is to arrive at some recommendations on how payload motions can be exploited in the attenuation of spacecraft nutation.

Chapter 2 - Mathematical Model

2.1 Model Description

The model used for this study is shown in Fig 2.1. It consists of a uniform, homogeneous, cylindrical rigid body A to which is attached four uniform slender rods B, C, D, and E. A has diameter R, height H and mass MA, and represents the main bus of a spacecraft. Each of the rods B, C, D, and E has length L and mass M and is connected to the cylinder through a one degree-of-freedom hinge, whose axis is perpendicular to the axis of the cylinder. The rods act as payloads and are arranged in such a way that their axes of relative rotation all lie on a plane which is at a distance P from the diametral plane containing the mass center A* of the cylinder. The mass center of the entire system is located at S*.

For this study, the four rods B, C, D, and E are assumed to be distributed symmetrically around the cylinder as shown, and each can oscillate about its axis of relative rotation, which is perpendicular to the axis of the cylinder. The symmetrical arrangement of the rods guarantees that the mass center of the entire system always lies on the axis of the cylinder, so that mass center offset is not an issue in the system's dynamics. The investigation then consists essentially of the following steps: (a) In a torque free space environment, the cylinder is given an initial spin as well as transverse angular velocity, with the spin angular speed being much larger than the transverse angular speed, resulting in coning or nutation of the cylinder; (b) one or more of the four rods is then made to oscillate relative to the cylinder at a certain frequency; (c) the nutation angle of the cylinder is monitored over time to find out whether the oscillatory motion of the rod(s) has any noticeable effect on the time history of the cylinder's nutation angle.

Figure 2.1 - Spacecraft Model

Abbildung in dieser Leseprobe nicht enthalten

2.2 Equations of motion

To begin this analysis, the equations of motion of the system for the general configuration shown in Fig. 2.1 are derived. The method used is that of Kane. Kane's equation gives the dynamical equations of a system of rigid bodies as¹⁴:

Abbildung in dieser Leseprobe nicht enthalten

where Fr is the generalized active force corresponding to the rth generalized speed, Fr* is the generalized inertia force, and n is the number of degrees of freedom of the system.

2.2.1 Kinematics

The interest in this study is in the attitude behavior (not translational motion) of the system of Fig. 2.1. Since the oscillatory motion of the rod(s) is prescribed, the system has three degrees of freedom. The three generalized speeds that characterize the system’s attitude motion are taken to be:

Abbildung in dieser Leseprobe nicht enthalten

where ω is the inertial angular velocity of body A, and the unit vectors a1, a2 , a3 are fixed in A and directed as shown in Fig 2.1. Thus,

Abbildung in dieser Leseprobe nicht enthalten

Hence the generalized speeds are really the a1, a2, a3 components of the inertial angular velocity of the cylinder A. The mass center S* of the system is assumed to be fixed in space and the oscillations of the payloads relative to body A are assumed to be defined by the function:

Abbildung in dieser Leseprobe nicht enthalten

where X is the amplitude of oscillation and Ω is the frequency of oscillation.

The inertial angular velocities ω B, ω C, ω D, and ω E of the payloads B, C, D, and E are given by:

Abbildung in dieser Leseprobe nicht enthalten

Relevant linear velocities are those of the mass center A* of body A, and the mass centers B*, C*, D*, and E* of payloads B, C, D, and E respectively. These are given by:

Abbildung in dieser Leseprobe nicht enthalten

where rS*A*, rS*B*, rS*C*, rS*D*, and rS*E* are the position vectors from the point S* to A*, B*, C*, D*, and E* respectively.

Angular acceleration for each of the rigid bodies in Fig 2.1 is obtained by taking the time derivative of each of the equations (2.3), (2.5) and (2.6) in an inertial reference frame. Similarly, the acceleration of each of the mass centers is obtained by taking the time derivative of each of the equations (2.7) - (2.11) in an inertial reference frame.

2.2.2 Kinetics

Assuming that no external torque is applied to the system as it moves in space, the generalized active forces for the system are given by¹⁴:

Abbildung in dieser Leseprobe nicht enthalten

where T B, T C, T D, and T E are the torques that drive the oscillatory motion of B, C, D, and E respectively, and the subscripted angular velocities represent the corresponding partial angular velocities. The generalized inertia forces are given by¹⁴:

Abbildung in dieser Leseprobe nicht enthalten

where the superscripted a vectors are accelerations of the various mass centers, the

subscripted velocity/angular velocity vectors are the partial velocities/angular velocities, H represents the central angular momentum whose derivative is taken in an inertial reference frame as signified by the dot over the H vector.

Equations (2.5) - (2.13) are assembled with the aid of the software AUTOLEV, and the equations of motion were then generated by implementing equation (2.1) with the same software. The AUTOLEV code used to generate the equations of motion is shown in appendix A. The resulting equations of motion for the system of Fig 2.1 in its most general configuration are given by:

Abbildung in dieser Leseprobe nicht enthalten

Chapter 3 - Case Studies

3.1 Numerical Integration

The equations of motion derived in chapter 2 and shown as equations (2.14) - (2.16), are clearly too complex to be dealt with analytically. A numerical approach is therefore adopted. Since the goal of this study is to assess the impact of payload motion on the nutational behavior of a spinning spacecraft, the solution strategy is as follows:

a. For a set of parameters representative of a real spacecraft, and for specified motion of one or more of the payloads, the equations of motion are integrated numerically using MATLAB;
b. The values obtained for the generalized speeds from the numerical integration are used to determine the value of the nutation angle for each integration interval;
c. The time history of the nutation angle is recorded in graphical form using MATLAB;
d. A parametric study is conducted by varying the mass properties, location, and motion of each payload to determine how these influence the nutation angle.

3.2 Stationary Payload (Case 0)

To begin the numerical study, the following parameter values are assigned to the cylinder A and the slender rods B, C, D, and E:

Abbildung in dieser Leseprobe nicht enthalten

The mass and dimensions chosen for the cylinder are such as to result in spin and transverse moments of inertia values of a typical axisymmetric spinner that has been flown in the past. It is recalled that, because of the assumption of axisymmetry, the spin moment of inertia of the cylinder is given by the expression¹⁵:

Abbildung in dieser Leseprobe nicht enthalten

and the central transverse moment of inertia is given by:

Abbildung in dieser Leseprobe nicht enthalten

For the slender rods, the axial moment of inertia is: [Abbildung in dieser Leseprobe nicht enthalten]

and the transverse moment of inertia is given by¹⁵:

Abbildung in dieser Leseprobe nicht enthalten

First, the case is considered (case 0) in which P = 0 and the rods do not move relative to the cylinder, so that β = 0. The initial conditions used for the purposes of numerical integration for this case are u1 (0) = u2 (0) = 0.05 rad/s and u3 (0) = 0.5 rad/s. Case 0 thus corresponds to the situation where S* coincides with A*, and the cylinder A spins about its axis while its axis of spin cones about a fixed axis in space as shown in Fig 3.1. In this well known case, the nutation angle shown as angle θ in the figure, remains constant and is given by the expression¹:

Abbildung in dieser Leseprobe nicht enthalten

where I1 is the combined central transverse moment of inertia of the cylinder A and the four rods B, C, D, and E all viewed as one rigid body, with the rods stationary with

Figure 3.1 - Nutation Angle for Case 0

Abbildung in dieser Leseprobe nicht enthalten

respect to the cylinder and aligned with the cylinder axis. I3 is the axial moment of inertia of the combined system, ω 3 = u3 (0) and ω 12 = [u1 (0)[[2]] + u2 (0)[[2]]][[1]]/[[2]]. Using parallel axes theorem,

Abbildung in dieser Leseprobe nicht enthalten

If u1 (0) and u2 (0) are kept at 0.05 rad/s and the initial spin rate is increased to u3 (0) = 1.0 rad/s, then θ becomes 2.518 [[0]].

The MATLAB code that is used to integrate the equations of motion (2.14) - (2.16) and to obtain the value of the nutation angle is given in appendix B. To verify the accuracy of this code as well as that of the equations of motion, the code is used to generate the time history of the nutation angle for case 0 described above. Plots of the nutation angle produced by this case are shown below in Fig 3.2 for a spin rate of 0.5 rad/s and in Fig 3.3 for a spin rate of 1.0 rad/s.

It is clear from these plots that the results obtained from the numerical integration code are identical to those obtained analytically above. Therefore these results give us confidence that the equations of motion derived in chapter 2 as well as the numerical integration code are correct. Having now established the reliability of the two critical tools (equations of motion and numerical integration code) needed for the completion of this study, the investigation outlined in section 3.1 can proceed.

Figure 3.2 - Code Verification for u 3 (0) = 0.5

Abbildung in dieser Leseprobe nicht enthalten

Figure 3.3 - Code Verification for u 3 (0) = 1

Abbildung in dieser Leseprobe nicht enthalten

3.3 Single Payload

In this section, the effect on spacecraft nutation of moving only one of the slender rods (payload) relative to the main spacecraft body (cylinder A) is examined. Hence, one of the rods B, C, D, or E is given a prescribed motion relative to body A while the others remain stationary and serve simply to maintain the axisymmetric nature of the entire system. Thus, any effect that may be observed will indeed be due to the motion imparted by the rod and not because of any static or dynamic imbalance in the system.

3.3.1 One oscillating payload (Case 1)

In this first case involving the motion of the payloads, the parameter values and initial conditions are kept the same as in Case 0 (see section 3.2 above), but one of the rods is made to oscillate back and forth relative to the main spacecraft body following the relation:

Abbildung in dieser Leseprobe nicht enthalten

Only one rod oscillates, the other three remain fixed to the cylinder. Fig 3.4 shows the time history of the nutation angle as obtained by the numerical procedure described in section 3.1. The light solid line is for the case where only rod D is allowed to oscillate (about a1 axis), while the dotted line corresponds to the situation where only rod B is allowed to oscillate (about a2 axis). The dark solid line represents the case when all payloads are stationary (Case 0).

Figure 3.4 - Nutation angle for one rod oscillating

Abbildung in dieser Leseprobe nicht enthalten

The following observations are made from the plot:

1. Except for a phase shift, the effect of moving one rod is pretty much identical to that of moving another rod. This is obviously to be expected due to the symmetry of the system.
2. In the case of stationary payload, the nutation angle was constant. Now, with an oscillating payload, the nutation angle varies sinusoidally with time, though the amplitude of oscillation is quite small (~ 0.2°).
3. The mean value of the nutation angle shown by the ‘Y’ in the data box for the case of oscillating payload (4.908°) is slightly less than the constant value of the nutation angle for the case of stationary payloads (5.027°).
4. It appears that moving the payload did some good by bringing down the average value of the nutation angle. However, the improvement (decrease) in the nutation angle is small.

The above result indicates that it might be worthwhile to investigate other motion scenarios for the payload to check for the possibility of further and more substantial decrease in the nutation angle. For example, increasing the mass of the moving payload might further reduce nutation. Varying the frequency and/or amplitude of payload oscillation could have noticeable effects on the spacecraft nutation.

3.3.2 Varying Payload Inertia (Case 2)

In this subsection, the inertia of the moving payload is increased, and the effect on nutation angle is examined. To increase payload inertia, the mass of the moving payload is simply doubled. To maintain the axisymmetric nature of the system, the mass of each of the slender rods (B, C, D, and E) is doubled while the mass of the main spacecraft body (cylinder) is adjusted so that the entire system maintains the same overall mass as in Case 1. Hence the assigned mass values are:

Abbildung in dieser Leseprobe nicht enthalten

the other parameter values as well as the initial conditions are kept the same as in Case 1. Numerical integration of the equations of motion yields the time histories shown in figure 3.5 for the spacecraft nutation angle. The light solid line is for a payload mass of 100 kg and the dotted line is for a payload mass of 200 kg. It is clear from the figure that an increase in the moment of inertia of the moving payload does in fact bring about a further decrease in the nutation angle. In this specific case, changing the mass of the oscillating payload from 100 kg to 200 kg reduces the average nutation angle from 4.909° to 4.633°, and the nutation angle always remains below the constant value of 5.027° obtained for the case of stationary payload.

Figure 3.5 - Effect of inertia on nutation angle

Abbildung in dieser Leseprobe nicht enthalten

3.3.3 Varying the frequency of payload oscillation (Case 3)

So far, the motion prescribed for the moving rod relative to the main cylinder A has been an oscillatory motion governed by the relationship:

Abbildung in dieser Leseprobe nicht enthalten

with Ω = 0.25 rad/s. To check for a possible effect of rod oscillation frequency on spacecraft nutation, other values - Ω = 0.1, 0.5, - are now assigned to payload oscillation frequency. The mass of each payload is returned to M = 100 kg, and all other parameters and initial conditions (except for the new values of Ω) are the same as in Case 1 (subsection 3.3.1).

Fig 3.6 shows the time history of the nutation angle for the above values of payload oscillation frequency. The dotted line in each plot corresponds to the mean nutation angle and the mean values are shown by ‘Y’ in the data boxes.

Figure 3.6 does not appear to show any clear pattern. For example, increasing the rod’s oscillation frequency from 0.1 rad/s to 0.25 rad/s results in (top and middle graph) a decrease in the frequency of oscillation of the nutation angle as well as a decrease in the average nutation angle. A further increase in Ω from 0.25 rad/s to 0.5 rad/s (middle and bottom graph) now leads to increase in both the mean value and the frequency of the nutation angle. So based on figure 3.6 alone, it is clear that payload oscillation frequency has some effect on nutation angle, but the nature of its effect is not predictable from the figure.

Figure 3.6 - Effect of payload oscillation frequency on nutation angle

Abbildung in dieser Leseprobe nicht enthalten

There are a number of frequencies that are of interest in the study of the dynamics of spinning rigid bodies. Consider an axisymmetric rigid body that is moving in a torquefree environment. Imagine that to start with, the body is spun up and given some small transverse angular velocity. It is known that the body will maintain a constant inertial spin rate n, and its attitude motion will be governed by:

Abbildung in dieser Leseprobe nicht enthalten

where ω 1 and ω 2 are the body’s transverse angular velocity components and λ is given by ¹:

Abbildung in dieser Leseprobe nicht enthalten

with I3 being the spin inertia of the body and I1 the body’s transverse moment of inertia. Geometrically, the motion of the body may be visualized as shown in figure 3.7.

Figure 3.7 - Motion of Spacecraft

Abbildung in dieser Leseprobe nicht enthalten

The angular momentum vector H remains fixed in inertial space while the precession plane P rotates about H at a rate ρ. At the same time, the body spins relative to P at the rate λ. The angular momentum vector H, the symmetry axis X3 of the body, and the inertial angular velocity vector of the body all remain in plane P at all times. θ is the nutation angle, and ρ is given by:

Abbildung in dieser Leseprobe nicht enthalten

There are at least three frequencies that are of interest in the motion described above: (a) the body’s spin frequency λ relative to the precession plane P, (b) the spin frequency ρ of the precession plane P in an inertial frame, and (c) the inertial spin rate n of the body. If the cylinder A and the four rods B, C, D, and E are considered as one axisymmetric spinning body, and using the same mass/inertia/geometric properties given for the system in the study of Case 0 above, I3 = 5600 kg-m[[2]] and I1 = 3482.7 kg-m[[2]]. Hence, for a value of ݊ = 0.5 rad/s,

Abbildung in dieser Leseprobe nicht enthalten

Looking back at figure 3.6, it appears that when Ω is close in value to and less than λ, the nutation angle is reduced, but when Ω is far removed from λ, there is hardly any effect on the nutation angle.

3.3.4 Varying the amplitude of payload oscillation (Case 4)

Finally, in this subsection, the effect of payload oscillation amplitude on spacecraft nutation is studied. Keeping the system parameters and the initial conditions the same as in Case 1 above, the amplitude X of payload motion is varied as follows: X = 0.25, 0.75, 1 rad

Fig 3.8 shows time plots for the spacecraft’s nutation angle for the various values of X.

The plots indicate that the mean value of the nutation angle decreases as payload oscillation amplitude is increased. However the excursions of the nutation angle relative to its mean value do also increase. So it is not necessarily a good idea to make X as high as possible.

Figure 3.8 - Effect of payload oscillation amplitude on nutation angle

Abbildung in dieser Leseprobe nicht enthalten

This section is now concluded by examining a case that combines all the above recommendations. The parameters are the same as in Case 2, but Ω = λ = 0.31 rad/s and X = 1 rad. Figure 3.9 shows the results, from which it is clear that the mean nutation angle of 4.293° is much smaller than what was obtained for the case of stationary payloads (5.027°).

Figure 3.9 - Recommendations for reducing nutation

Abbildung in dieser Leseprobe nicht enthalten

3.4 Motion of Multiple Payloads

3.4.1 Two Oscillating Payloads (Case 5)

First the effect on the nutation angle, of having two opposite rods oscillate is considered. The frequency and amplitude of rod oscillation are taken to be 0.25 rad/s and 0.5 rad respectively. All other parameters are the same as in Case 1. Fig 3.10 shows the nutation angle for two opposite rods oscillating both in phase and 180° out of phase. The solid line is for the case of rods oscillating in phase and the dotted line is for the case of rods oscillating 180° out of phase.

Figure 3.10 - Nutation Angle for Two Oscillating Opposite Rods

Abbildung in dieser Leseprobe nicht enthalten

Compared to the case of stationary rods (Case 0), the figure shows a decrease in the mean nutation angle when two opposite rods oscillate in phase. However when two opposite rods oscillate 180° out of phase, the nutation angle is approximately constant at the same value obtained for stationary payloads indicating that the effect of moving one of the rods is effectively cancelled by the motion of the second rod.

It is useful to compare the case where the two rods oscillate to the one presented in Case 2 (subsection 3.3.2). The two cases can be viewed as being equivalent, in the sense that the total mass as well as the inertia of the oscillating rod(s) are the same in both cases. In one case, only one rod oscillates and it does so at one end of the axis of oscillation. In the other case, it is as if the rod is split into two equal rods with one oscillating at each end of the axis of oscillation. A comparison of fig 3.5 for case 2 with fig 3.10 above shows that the mean value of the nutation angle is approximately the same in bothe cases. Case 2, where a single large payload moves does show lower maximum nutation angle that the case of fig 3.10. It is therefore advantageous to maintain one payload with increased inertia in motion instead of moving multiple payloads. Next the effect of the oscillation of two adjacent rods on nutation angle is considered. Fig 3.11 shows the nutation angle for the case where two adjacent rods are made to oscillate both in phase and 180° out of phase. The solid line is for the case of rods oscillating in phase while the dotted line is for the case of rods oscillating 180° out of phase oscillation. It is seen that then two adjacent rods oscillate 180° out of phase, the nutation angle increases compared to the case of stationary payloads. Hence this case is not useful to us. However there is a decrease in the nutation angle when two adjacent rods oscillate in phase. It is also found that the decrease in mean nutation angle for the case of two adjacent payloads oscillating in phase is slightly smaller than the case for two opposite payloads oscillating in phase.

Figure 3.11 - Nutation Angle for Two Oscillating Adjacent Rods

Abbildung in dieser Leseprobe nicht enthalten

3.4.2 Four Oscillating Payloads (Case 6)

Finally, the effect on spacecraft nutation angle of having all the payloads oscillate is assessed. The frequency and amplitude selected for the rods oscillation are 0.25 rad/s and 0.5 rad respectively. All other parameters are the same as in Case 1.

Three types of rod oscillations are considered - all rods oscillate in phase, opposite rods oscillate in phase with a 180° phase difference between the two opposite pairs, and finally adjacent rods oscillate in phase with a 180° phase difference between the two adjacent pairs (B,D oscillate in phase and C,E oscillate in phase with a phase difference of 180° between the two pairs).

Figure 3.12 - Nutation Angle for All Rods Oscillating

Abbildung in dieser Leseprobe nicht enthalten

Fig 3.12 shows the comparison of the nutation angle between the above three types of rod oscillations. The solid light line is for the case of all rods oscillating in phase, the dashed line is for the case of opposite rods oscillating in phase with a 180° phase difference between the two opposite pairs, and the solid dark line is for the case of adjacent rods oscillating in phase with a 180° phase difference between the two adjacent pairs. It is seen from the graph that when adjacent rods oscillate in phase with a 180° phase difference between the two adjacent pairs as shown by the dark solid line, the nutation angle is constant and is similar to fig 3.2 for stationary payloads. It is also seen that the mean nutation angle for the case of all rods oscillating in phase as shown by the light solid line produces the best value of the mean nutation angle of all the cases involving the oscillation of multiple payloads.

3.5 Varying payload location

All of the case studies considered so far assume that the mass center of the main spacecraft body (the cylinder) and the mass center of each of the payloads (rods) all lie in the same plane. Hence, P = 0 for all the cases considered so far. In this section, the effect on spacecraft nutation of having P ≠ 0 is studied. In other words, each payload is shifted axially a distance P (≠ 0) away from the mass center of the cylinder, and the payload axis is then kept stationary at this new position. The values chosen for P are, in meters:

Abbildung in dieser Leseprobe nicht enthalten

First, the case where all the payloads remain stationary is considered. Since the entire system remains axisymmetric, the nutation angle is constant and is once more given by¹:

Abbildung in dieser Leseprobe nicht enthalten

Since P ≠ 0, the axial location of the mass center of the entire system from A* is denoted as Z and is given by:

Abbildung in dieser Leseprobe nicht enthalten

Assuming that mass values are kept at

Abbildung in dieser Leseprobe nicht enthalten

Finally, the values of the combined central transverse moment of inertia of the system I1 and the combined axial moment of inertia of the system I3 are calculated. The other terms in equation 3.17 are given by ω 3 = u3 (0) and ω 12 = [u1 (0)[[2]] + u2 (0)[[2]]][[1]]/[[2]]. Using parallel axes theorem,

Abbildung in dieser Leseprobe nicht enthalten

where IA is the central transverse moment of inertia of cylinder A as given by equation 3.2, JA is the spin moment of inertia of the cylinder A as given by equation 3.1 and I is the transverse moment of inertia of each payload as given by equation 3.4. Hence on substitution of all values into equation 3.20 and 3.21,

Abbildung in dieser Leseprobe nicht enthalten

All the above values are substituted in equation 3.17 to get,

Abbildung in dieser Leseprobe nicht enthalten

Finally, substituting the values of P as selected above in equation 3.24, the value of the constant nutation angles are obtained as:

θ = 5.056° (for P = 0.25 m) θ = 5.145° (for P = 0.5 m)

θ = 5.295° (for P = 0.75 m) θ = 5.504° (for P = 1 m)

Hence as P is increased from zero, the nutation angle also increases.

The new values of frequency λ as discussed in chapter 3 are also calculated for each value of P using equation 3.13 and are as follows:

λ = 0.2992 rad/s (for P = 0.25 m) λ = 0.2852 rad/s (for P = 0.5 m)

λ = 0.2630 rad/s (for P = 0.75 m) λ = 0.2338 rad/s (for P = 1 m)

Due the conclusion from Case 3 (subsection 3.3.3), the selected value of the frequency of oscillation of the payloads for the simulations is kept below the value of λ shown above.

3.5.1 Single Oscillating Payload (Case 7)

One rod is now allowed to oscillate while the others remain stationary. The parameters and initial conditions are the same as in case 1, and figures 3.13 - 3.16 show the nutation angle for P = 0.25, 0.5, 0.75 and 1 m respectively. In each figure the solid dark line is for the case where all the payloads are held stationary in their displaced positions and the solid light line is for the case of one oscillating payload. The mean values of nutation angle are also shown by the Y values in the data boxes. The following observations are made from the plots:

1. For all the above values of P, the nutation angle for the case of stationary payloads does in fact increase with P, as was established analytically above.

2. For P = 0.25, 0.5, 0.75 m, the mean nutation angle for the case of one oscillating payload is found to be less than that obtained for P = 0 from Case 1(4.908°). However for P = 1 m, the mean nutation angle when one payload is made to oscillate is found to be higher than that obtained for P = 0.

3. The mean nutation angle for one oscillating payload is seen to increase as P increases and overall, the smallest mean nutation angle occurs when P = 0.25 m.

Figure 3.13 - Nutation Angle for P = 0.25 m (One Oscillating Rod)

Abbildung in dieser Leseprobe nicht enthalten

Figure 3.14 - Nutation Angle for P = 0.5 m (One Oscillating Rod)

Abbildung in dieser Leseprobe nicht enthalten

Figure 3.15 - Nutation Angle for P = 0.75 m (One Oscillating Rod)

Abbildung in dieser Leseprobe nicht enthalten

Figure 3.16 - Nutation Angle for P = 1 m (One Oscillating Rod)

Abbildung in dieser Leseprobe nicht enthalten

4. Therefore as seen above keeping P as small as possible but not equal to zero results in the biggest decrease in the nutation angle. This suggests that the payloads should be axially located as close as possible to the mass center A* of the cylinder.

3.5.2 Two Oscillating Payloads (Case 8)

In this case, two opposite payloads are made to oscillate in phase and the other two are kept stationary while the plane of location of the payloads remains offset from the mass center A*. System parameters and initial conditions are the same as in Case 5 for two oscillating payloads. The values considered for P are the same as in Case 7 (0.25, 0.5, 0.75 and 1 m respectively). Figures 3.17 - 3.20 show the nutation angle for the above values of P. In each figure the solid dark line is for the case of stationary payloads and the solid light line is for the case of two oscillating payloads. The mean values of nutation angle are also shown by the Y values in the data boxes. The observations that result from the above plots are the same as those for single oscillating payload (case 7).

Abbildung in dieser Leseprobe nicht enthalten

Figure 3.17 - Nutation Angle for P = 0.25 m (Two Opposite Oscillating Rods)

Abbildung in dieser Leseprobe nicht enthalten

Figure 3.18 - Nutation Angle for P = 0.5 m (Two Opposite Oscillating Rods)

Abbildung in dieser Leseprobe nicht enthalten

Figure 3.19 - Nutation Angle for P = 0.75 m (Two Opposite Oscillating Rods)

Abbildung in dieser Leseprobe nicht enthalten

Figure 3.20 - Nutation Angle for P = 1 m (Two Opposite Oscillating Rods)

Chapter 4 - Summary of Results

4.1 Characteristics of Nutation Angle

The plots shown by the various cases (Case 0 - Case 9) in chapter 3 highlighted the benefits of impacting oscillatory motion to the payloads on spacecraft nutation angle and hence stability. The two important characteristics of the spacecraft’s nutation angle are: (a) its value, when nutation angle is constant and (b) when the nutation angle fluctuates with time, its mean value and amplitude. Therefore to have a better understanding of the potential benefits of payload motion on spacecraft nutation, plots of the nutation amplitude and mean value against various payload motion parameters such as frequency of payload oscillation, amplitude of payload oscillation, payload location and payload inertia are generated and displayed in figures 4.1 - 4.8. Results are given only for the case of a single oscillating payload. The parameter values and system initial conditions are kept the same as in Case 1. The values of the nutation amplitude and mean value used for all the plots are shown in appendix C.

4.2 Effect of Payload Oscillation Frequency

The payload oscillation amplitude is taken to be X = 0.5 rad and the payload location is defined by P = 0. Figures 4.1 and 4.2 show the effect of payload oscillation frequency on the amplitude and mean value of nutation angle. It is recalled here that 0.3 is the value of the frequency λ discussed in Case 3 (subsection 3.3.3).

Abbildung in dieser Leseprobe nicht enthalten

Figure 4.1 - Effect of Payload Oscillation Frequency on Nutation Amplitude

Abbildung in dieser Leseprobe nicht enthalten

Figure 4.2 - Effect of Payload Oscillation Frequency on Mean Nutation Angle

The following observations are made from the above plots:

1. From figure 4.1, it is seen that the maximum nutation amplitude occurs when Ω = λ = 0.3 rad/s. This indicates some kind of resonance is occurring around this region. The amplitude at all other regions remains fairly small.

2. From figure 4.2, it is seen that the smallest mean value and hence the biggest reduction in nutation angle also occurs when Ω = λ = 0.3 rad/s which seems to be the cutoff point for reduction in nutation angle. In addition, for values of Ω > λ, the mean value is higher than the value obtained for the case where all the payloads are fixed to the spacecraft (shown in the graph at Ω = 0).

3. The preferred scenario is of course to have small values of both the mean value and the nutation amplitude. Therefore from the above two figures, it is clear that selecting a value of Ω very close to λ but less than λ, will achieve this desirable combination of low mean value and low nutation amplitude.

4.3 Effect of Payload Oscillation Amplitude

The payload oscillation frequency is taken to be Ω = 0.25 rad/s and the payload location is defined by P = 0. Figure 4.3 and 4.4 show the effect of payload oscillation amplitude on the amplitude and mean value of nutation angle. The following observations are made from the plots:

1. From figure 4.3, a clear pattern emerges where it is seen that increasing the payload oscillation amplitude increases the nutation amplitude.

Abbildung in dieser Leseprobe nicht enthalten

Figure 4.3 - Effect of Payload Oscillation Amplitude on Nutation Amplitude

Abbildung in dieser Leseprobe nicht enthalten

Figure 4.4 - Effect of Payload Oscillation Amplitude on Mean Nutation Angle

2. From figure 4.4 once again a clear pattern is seen where increasing the payload oscillation amplitude decreases the mean nutation angle.

3. These two figures can be used to select a payload oscillation amplitude that will keep both the mean nutation angle and nutation amplitude small.

4.4 Effect of Payload Location

The payload oscillation frequency Ω is kept smaller than the corresponding value of the frequency λ as calculated in chapter 3 (section 3.5) and the payload oscillation amplitude is defined by X = 0.5 rad. Figure 4.5 and 4.6 show the effect of location of the payloads on the amplitude and mean value of nutation angle. The following observations are made from the plots below:

1. From figure 4.5, it is seen that maximum nutation amplitude occurs when P = 0.25m.

2. From figure 4.6, once again the smallest mean nutation angle occurs when P = 0.25m, which appears to be a cutoff point for reduction in nutation angle. Any value of P above this causes an increase in the mean nutation angle.

3. For good attitude behavior and to keep both the amplitude and the mean angle of the spacecraft nutation, it is recommended to select a value of P close to and smaller than 0.25 m as seen from both the figures below.

Abbildung in dieser Leseprobe nicht enthalten

Figure 4.5 - Effect of Payload Location on Nutation Amplitude

Abbildung in dieser Leseprobe nicht enthalten

Figure 4.6 - Effect of Payload Location on Mean Nutation Angle

4.5 Effect of Payload Inertia

The payload oscillation frequency is taken to be Ω = 0.25 rad/s, the payload oscillation amplitude is taken as X = 0.5 rad and the payload location is defined by P = 0. Figure 4.7 and 4.8 show the effect of payload inertia on the amplitude and mean value of nutation angle. The following observations are made from these plots:

1. From figure 4.7, it is seen that an increase in the payload inertia increases the nutation amplitude.
2. From figure 4.8, it is seen that the mean nutation angle decreases upon increasing the payload inertia.
3. Here again, these two figures can be used to choose an optimum value of payload inertia which results in a small value of both the amplitude and the mean value of nutation angle.

Abbildung in dieser Leseprobe nicht enthalten

Figure 4.7 - Effect of Payload Inertia on Nutation Amplitude

Abbildung in dieser Leseprobe nicht enthalten

Figure 4.8 - Effect of Payload Inertia on Mean Nutation Angle

Chapter 5 -Conclusions

5.1 Summary and Conclusions

Spacecraft systems almost always include subsystems such as imaging or scanning instruments that are designed to move relative to the main spacecraft body during some part of the spacecraft’s mission. Depending on the location and the relative mass or inertia of such moving parts, their motion can have noticeable influence on the overall translational and rotational motion of the main body. The main goal of the research project is to investigate the effect of oscillatory motion of spacecraft-mounted payload on the attitude behavior of the spacecraft, especially the spacecraft’s nutational motion. The study is designed to focus on the actual effect of the motion of the moving parts and not on accompanying effects that may be due to the static or dynamic imbalance of the moving parts. The results obtained show that oscillating motion of a payload relative to the main spacecraft will, in general, attenuate the spacecraft’s nutational motion. The degree of nutation attenuation depends on several factors. The higher the ratio of the payload’s moment of inertia to that of the spacecraft, the larger is the reduction in spacecraft mean nutation angle. The amplitude of oscillation of the payload is also found to influence spacecraft nutation. As this amplitude is increased, while maintaining the same frequency, spacecraft nutation (mean nutation angle) is found to decrease. The frequency of the payloads motion is found to play an important role in how the spacecraft nutation is affected. When the payload’s motion frequency is equal to or close to the frequency λ at which the spacecraft spins relative to the precession plane, spacecraft nutation is reduced to its lowest value.

The results obtained in this study lead to the conclusion that relative payload motion can indeed be exploited to supplement the nutation damping actions of a spacecraft’s nutation damper system. Payload motion can also serve as a backup for a spacecraft’s nutation damper system in case of partial or total malfunction of such a system.

5.2 Recommendations for Future Work

The study presented here can be expanded on many fronts. First, this study assumed that the oscillatory motion of each payload is without any damping. Real systems always have non-zero damping. So, it is suggested that future work include damping for the rotational oscillatory motion of the payload. Our suspicion is that adding damping may improve the nutation damping effects observed. The payload motion studied here is restricted to simple rotational oscillations. It may be worthwhile to study a combination of rotational motion about payload axis and translational oscillations of the payloads axis parallel to the cylinder axis. Finally, this work is restricted to single spinners. Future work should include the effect of payload motion on dual-spin spacecraft.

References

1. Marshall H. Kaplan., “Modern Spacecraft Dynamics & Control”. 1976: John Wiley & Sons.

2. Mingori, D.L., “Effects of Energy Dissipation on the Attitude Stability of Dual- Spin Satellites”, AIAA Journal, Vol.7, No.1, 1969.

3. Linkins, P.W., “Attitude Stability Criteria for Dual-Spin Spacecraft”, Journal of Spacecraft and Rockets, Vol.4, No.12, 1967.

4. Bainum, P.M., Fuechsel, P.G., and Mackison, D.L., “Motion and Stability of a Dual-Spin Satellite with Nutation Damping”, Journal of Spacecraft and Rockets, Vol.7, No.6, 1970.

5. Auelmann, R.R., and Lane, P.T., “Design and Analysis of Ball-In-Tube Nutation Dampers”, Proceedings of the Symposium on Attitude Stabilization and Control of Dual-Spin Spacecraft, 1968.

6. Taylor, R.S., and Conway, J.J., “Viscous Ring Precession Damper for Dual-Spin Spacecraft”, Proceedings of the Symposium on Attitude Stabilization and Control of Dual-Spin Spacecraft, 1968.

7. Eke, F.O., “Effectiveness of Large Booms as Nutation Dampers for Spin Stabilized Spacecraft”, Proceedings, Flight Mechanics/Estimation Theory Symposium, Greenbelt MD, 1991.

8. Michael D. Griffin and James R. French., “Space Vehicle Design”. 2004: Reston, VA.

9. Smay, J.W., and Slafer, L.I., “Dual-Spin Spacecraft Stabilization using Nutation Feedback and Inertia Coupling”, Journal of Spacecraft and Rockets, Vol.13, 1976.

10. McIntyre, J.E., and Gianelli, M.J., “The Effect of an Autotracking Antenna on the Nutational Stability of a Dual-Spin Spacecraft”, AIAA Paper 72-571, April, 1972.

11. Farrenkopf, R.L., “Stability of Dual-Spin Spacecraft Containing Tracking Payloads”, AIAA Paper 72-890, 1972.

12. Man, G.K., and Eke, F.O., “Effects of Payload Motions on the Nutational Stability of the Galileo Spacecraft”, Journal of Guidance, Control, and Dynamics, Vol.8, No.6, 1985.

13. Eke, F.O., and Eke, E.M., “Backup Nutation Damping Strategy for the Galileo Spacecraft”, ASME Journal of Dynamic Systems, Measurement, and Control, Vol.113, June 1991.

14. Kane, T.R., and Levinson, D.A., Dynamics: Theory and Applications, McGraw- Hill, 1985.

15. Meriam, J.L., and Kraige, L.G., “Engineering Mechanics - Dynamics”. 2008: John Wiley & Sons.

Appendix

A. Autolev code for obtaining equations of motion

(1) AUTOZ ON

(2) NEWTONIAN N

(3) POINTS SO

(4) BODIES A,B,C,D,E

(5) MASS A=Ma,B=M,C=M,D=M,E=M

(6) CONSTANTS H,R,P,L,X,OMEGA

(7) SPECIFIED OB'',OC'',OD'',OE''

(8) VARIABLES THETA

(9) MOTIONVARIABLES' U{3}'

(10) INERTIA A,Ma*R*R/16+Ma*H*H/12,Ma*R*R/16+Ma*H*H/12,Ma*R*R/8

-> (11) Z1 = Ma*(R^2+1.333333*H^2)

-> (12) Z2 = Ma*R^2

-> (13) I_A_AO>> = 0.0625*Z1*A1>*A1> + 0.0625*Z1*A2>*A2> + 0.125*Z2*A3>*A3>

(14) INERTIA B,M*L*L/12,M*L*L/12,0

-> (15) Z3 = M*L^2

-> (16) I_B_BO>> = 0.08333333*Z3*B1>*B1> + 0.08333333*Z3*B2>*B2>

(17) INERTIA C,M*L*L/12,M*L*L/12,0

-> (18) I_C_CO>> = 0.08333333*Z3*C1>*C1> + 0.08333333*Z3*C2>*C2>

(19) INERTIA D,M*L*L/12,M*L*L/12,0

-> (20) I_D_DO>> = 0.08333333*Z3*D1>*D1> + 0.08333333*Z3*D2>*D2>

(21) INERTIA E,M*L*L/12,M*L*L/12,0

-> (22) I_E_EO>> = 0.08333333*Z3*E1>*E1> + 0.08333333*Z3*E2>*E2>

(23) SIMPROT(A,B,2,OB)

-> (24) Z4 = COS(OB)

-> (25) Z5 = SIN(OB)

-> (26) A_B = [Z4, 0, Z5; 0, 1, 0; -Z5, 0, Z4]

(27) SIMPROT(A,C,2,OC)

-> (28) Z6 = COS(OC)

-> (29) Z7 = SIN(OC)

-> (30) A_C = [Z6, 0, Z7; 0, 1, 0; -Z7, 0, Z6]

(31) SIMPROT(A,D,1,OD)

-> (32) Z8 = COS(OD)

-> (33) Z9 = SIN(OD)

-> (34) A_D = [1, 0, 0; 0, Z8, -Z9; 0, Z9, Z8]

(35) SIMPROT(A,E,1,OE)

-> (36) Z10 = COS(OE)

-> (37) Z11 = SIN(OE)

-> (38) A_E = [1, 0, 0; 0, Z10, -Z11; 0, Z11, Z10]

(39) W_A_N>=U1*A1>+U2*A2>+U3*A3>

-> (40) W_A_N> = U1*A1> + U2*A2> + U3*A3>

(41) W_B_A>=OB'*A2>

-> (42) W_B_A> = OB'*A2>

(43) W_C_A>=OC'*A2>

-> (44) W_C_A> = OC'*A2>

(45) W_D_A>=OD'*A1>

-> (46) W_D_A> = OD'*A1>

(47) W_E_A>=OE'*A1>

-> (48) W_E_A> = OE'*A1>

(49) P_AO_BO>=-R/2*A2>+P*A3>

-> (50) P_AO_BO> = -0.5*R*A2> + P*A3>

(51) P_AO_CO>=R/2*A2>+P*A3>

-> (52) P_AO_CO> = 0.5*R*A2> + P*A3>

(53) P_AO_DO>=-R/2*A1>+P*A3>

-> (54) P_AO_DO> = -0.5*R*A1> + P*A3>

(55) P_AO_EO>=R/2*A1>+P*A3>

-> (56) P_AO_EO> = 0.5*R*A1> + P*A3>

(57) P_AO_SO>=CM(AO)

-> (58) Z12 = Ma + 4*M

-> (59) Z13 = M*P/Z12

-> (60) P_AO_SO> = 4*Z13*A3>

(61) V_SO_N>=0>

-> (62) V_SO_N> = 0>

(63) V2PTS(N,A,SO,AO)

-> (64) V_AO_N> = -4*Z13*U2*A1> + 4*Z13*U1*A2>

(65) V2PTS(N,A,SO,BO)

-> (66) Z14 = P - 4*Z13

-> (67) V_BO_N> = (Z14*U2+0.5*R*U3)*A1> - Z14*U1*A2> - 0.5*R*U1*A3>

(68) V2PTS(N,A,SO,CO)

-> (69) V_CO_N> = (Z14*U2-0.5*R*U3)*A1> - Z14*U1*A2> + 0.5*R*U1*A3>

(70) V2PTS(N,A,SO,DO)

-> (71) V_DO_N> = Z14*U2*A1> + (-Z14*U1-0.5*R*U3)*A2> + 0.5*R*U2*A3>

(72) V2PTS(N,A,SO,EO)

-> (73) V_EO_N> = Z14*U2*A1> + (0.5*R*U3-Z14*U1)*A2> - 0.5*R*U2*A3>

(74) H>=MOMENTUM(ANGULAR,SO)

-> (75) Z15 = Z1*U1

-> (76) Z16 = Z1*U2

-> (77) Z17 = Z2*U3

-> (78) Z18 = Z4*U1 - Z5*U3

-> (79) Z19 = Z4*U3 + Z5*U1

-> (80) Z20 = OB' + U2

-> (81) Z21 = Z3*Z18

-> (82) Z22 = Z3*Z20

-> (83) Z23 = Z6*U1 - Z7*U3

-> (84) Z24 = Z6*U3 + Z7*U1

-> (85) Z25 = OC' + U2

-> (86) Z26 = Z3*Z23

-> (87) Z27 = Z3*Z25

-> (88) Z28 = OD' + U1

-> (89) Z29 = Z8*U2 + Z9*U3

-> (90) Z30 = Z8*U3 - Z9*U2

-> (91) Z31 = Z3*Z28

-> (92) Z32 = Z3*Z29

-> (93) Z33 = OE' + U1

-> (94) Z34 = Z10*U2 + Z11*U3

-> (95) Z35 = Z10*U3 - Z11*U2

-> (96) Z36 = Z3*Z33

-> (97) Z37 = Z3*Z34

-> (98) Z38 = Z13*U2

-> (99) Z39 = Z13*U1

-> (100) Z40 = Z14*U2 + 0.5*R*U3

-> (101) Z41 = Z14*U1

-> (102) Z42 = R*U1

-> (103) Z43 = Z14*U2 - 0.5*R*U3

-> (104) Z44 = Z14*U2

-> (105) Z45 = -Z14*U1 - 0.5*R*U3

-> (106) Z46 = R*U2

-> (107) Z47 = 0.5*R*U3 - Z14*U1

-> (108) Z48 = Ma*Z13

-> (109) Z49 = M*Z14

-> (110) Z50 = M*R

-> (111) H> = (0.0625*Z15+0.5*Z50*Z42+2*Z49*Z41+16*Z48*Z39-Z49*Z45-Z49*Z47)*A1> + (0.0625*Z16+Z49*Z40+Z49*Z43+0.5*Z50*Z46+2*Z49*Z44+16*Z48*Z38)*A2> + (0.125*Z17+0.5*Z50*Z40+0.5*Z50*Z47-0.5*Z50*Z43-0.5*Z50*Z45)*A3> + 0.0 8333333*Z21*B1> + 0.08333333*Z22*B2> + 0.08333333*Z26*C1> + 0.08333333

*Z27*C2> + 0.08333333*Z31*D1> + 0.08333333*Z32*D2> + 0.08333333*Z36* E1> + 0.08333333*Z37*E2>

(112) THETA=ACOS(DOT(UNITVEC(H>),A3>))

-> (113) Z51 = 0.0625*Z15 + 0.5*Z50*Z42 + 2*Z49*Z41 + 16*Z48*Z39 - Z49*Z45 - Z49*Z47

-> (114) Z52 = 0.0625*Z16 + Z49*Z40 + Z49*Z43 + 0.5*Z50*Z46 + 2*Z49*Z44 + 16* Z48*Z38

-> (115) Z53 = 0.125*Z17 + 0.5*Z50*Z40 + 0.5*Z50*Z47 - 0.5*Z50*Z43 - 0.5*Z50* Z45

-> (116) Z54 = Z4*Z6 + Z5*Z7

-> (117) Z55 = Z4*Z7 - Z5*Z6

-> (118) Z56 = Z5*Z6 - Z4*Z7

-> (119) B_C = [Z54, 0, Z55; 0, 1, 0; Z56, 0, Z54]

-> (120) Z57 = Z5*Z9

-> (121) Z58 = Z5*Z8

-> (122) Z59 = Z4*Z9

-> (123) Z60 = Z4*Z8

-> (124) B_D = [Z4, -Z57, -Z58; 0, Z8, -Z9; Z5, Z59, Z60]

-> (125) Z61 = Z5*Z11

-> (126) Z62 = Z5*Z10

-> (127) Z63 = Z4*Z11

-> (128) Z64 = Z4*Z10

-> (129) B_E = [Z4, -Z61, -Z62; 0, Z10, -Z11; Z5, Z63, Z64]

-> (130) Z65 = Z7*Z9

-> (131) Z66 = Z7*Z8

-> (132) Z67 = Z6*Z9

-> (133) Z68 = Z6*Z8

-> (134) C_D = [Z6, -Z65, -Z66; 0, Z8, -Z9; Z7, Z67, Z68]

-> (135) Z69 = Z7*Z11

-> (136) Z70 = Z7*Z10

-> (137) Z71 = Z6*Z11

-> (138) Z72 = Z6*Z10

-> (139) C_E = [Z6, -Z69, -Z70; 0, Z10, -Z11; Z7, Z71, Z72]

-> (140) Z73 = Z8*Z10 + Z9*Z11

-> (141) Z74 = Z9*Z10 - Z8*Z11

-> (142) Z75 = Z8*Z11 - Z9*Z10

-> (143) D_E = [1, 0, 0; 0, Z73, Z74; 0, Z75, Z73]

-> (144) Z76 = Z51^2 + Z52^2 + Z53^2 + 0.006944444*Z21^2 + 0.006944444*Z22^2 +

0.006944444*Z26^2 + 0.006944444*Z27^2 + 0.006944444*Z31^2 + 0.0069444 44*Z32^2 + 0.006944444*Z36^2 + 0.006944444*Z37^2 + 0.01388889*Z22*Z27 + 0.01388889*Z31*Z36 + 0.1666667*Z22*Z52 + 0.1666667*Z27*Z52 + 0.1666 667*Z31*Z51 + 0.1666667*Z36*Z51 + 0.01388889*Z4*Z21*Z31 + 0.01388889* Z4*Z21*Z36 + 0.01388889*Z6*Z26*Z31 + 0.01388889*Z6*Z26*Z36 + 0.013888 89*Z8*Z22*Z32 + 0.01388889*Z8*Z27*Z32 + 0.01388889*Z10*Z22*Z37 + 0.01 388889*Z10*Z27*Z37 + 0.01388889*Z54*Z21*Z26 + 0.01388889*Z73*Z32*Z37

+ 0.1666667*Z4*Z21*Z51 + 0.1666667*Z6*Z26*Z51 + 0.1666667*Z8*Z32*Z52 + 0.1666667*Z9*Z32*Z53 + 0.1666667*Z10*Z37*Z52 + 0.1666667*Z11*Z37*Z53

- 0.1666667*Z5*Z21*Z53 - 0.1666667*Z7*Z26*Z53 - 0.01388889*Z57*Z21*Z32

- 0.01388889*Z61*Z21*Z37 - 0.01388889*Z65*Z26*Z32 - 0.01388889*Z69*Z26

*Z37

-> (145) Z77 = Z76^0.5

-> (146) Z78 = Z51/Z77

-> (147) Z79 = Z52/Z77

-> (148) Z80 = Z53/Z77

-> (149) Z81 = Z21/Z77

-> (150) Z82 = Z22/Z77

-> (151) Z83 = Z26/Z77

-> (152) Z84 = Z27/Z77

-> (153) Z85 = Z31/Z77

-> (154) Z86 = Z32/Z77

-> (155) Z87 = Z36/Z77

-> (156) Z88 = Z37/Z77

-> (157) THETA = ACOS(Z80+0.08333333*Z9*Z86+0.08333333*Z11*Z88-0.08333333*Z5* Z81-0.08333333*Z7*Z83)

(158) ALF_A_N>=DT(W_A_N>,N)

-> (159) ALF_A_N> = U1'*A1> + U2'*A2> + U3'*A3>

(160) A_SO_N>=DT(V_SO_N>,N)

-> (161) A_SO_N> = 0>

(162) A2PTS(N,A,SO,AO)

-> (163) Z89 = Z13*U2*U3

-> (164) Z90 = Z13*U1*U3

-> (165) Z91 = Z13*(U1^2+U2^2)

-> (166) A_AO_N> = (-4*Z13*U2'-4*Z90)*A1> + (4*Z13*U1'-4*Z89)*A2> + 4*Z91*A3>

(167) A2PTS(N,A,SO,BO)

-> (168) Z92 = 0.5*R*U1^2 + 0.5*U3*(R*U3+2*Z14*U2)

-> (169) Z93 = -Z14*U1^2 - 0.5*U2*(R*U3+2*Z14*U2)

-> (170) Z94 = U1*(R*U2-2*Z14*U3)

-> (171) A_BO_N> = (Z14*U2'+0.5*R*U3'-0.5*Z94)*A1> + (-Z14*U1'+Z92)*A2> + (-0.5

*R*U1'+Z93)*A3>

(172) A2PTS(N,A,SO,CO)

-> (173) Z95 = -0.5*R*U1^2 - 0.5*U3*(R*U3-2*Z14*U2)

-> (174) Z96 = 0.5*U2*(R*U3-2*Z14*U2) - Z14*U1^2

-> (175) Z97 = U1*(R*U2+2*Z14*U3)

-> (176) A_CO_N> = (Z14*U2'-0.5*R*U3'+0.5*Z97)*A1> + (-Z14*U1'+Z95)*A2> + (0.5* R*U1'+Z96)*A3>

(177) A2PTS(N,A,SO,DO)

-> (178) Z98 = U2*(R*U1-2*Z14*U3)

-> (179) Z99 = 0.5*R*U2^2 + 0.5*U3*(R*U3+2*Z14*U1)

-> (180) Z100 = -Z14*U2^2 - 0.5*U1*(R*U3+2*Z14*U1)

-> (181) A_DO_N> = (Z14*U2'+Z99)*A1> + (-Z14*U1'-0.5*R*U3'-0.5*Z98)*A2> + (0.5* R*U2'+Z100)*A3>

(182) A2PTS(N,A,SO,EO)

-> (183) Z101 = U2*(R*U1+2*Z14*U3)

-> (184) Z102 = -0.5*R*U2^2 - 0.5*U3*(R*U3-2*Z14*U1)

-> (185) Z103 = 0.5*U1*(R*U3-2*Z14*U1) - Z14*U2^2

-> (186) A_EO_N> = (Z14*U2'+Z102)*A1> + (0.5*R*U3'-Z14*U1'+0.5*Z101)*A2> + (-0 .5*R*U2'+Z103)*A3>

(187) ZERO=FR()+FRSTAR()

-> (188) Z104 = 0.0625*U1*Z16 - 0.0625*U2*Z15

-> (189) Z105 = 0.0625*U3*Z15 - 0.125*U1*Z17

-> (190) Z106 = 0.125*U2*Z17 - 0.0625*U3*Z16

-> (191) Z107 = OB'*U3

-> (192) Z108 = OB'*U1

-> (193) Z109 = -Z4*Z107 - Z5*Z108

-> (194) Z110 = Z4*Z108 - Z5*Z107

-> (195) Z111 = Z3*Z5

-> (196) Z112 = Z3*Z4

-> (197) Z113 = Z3*Z109

-> (198) Z114 = Z3*OB''

-> (199) Z115 = 0.08333333*Z18*Z22 - 0.08333333*Z20*Z21

-> (200) Z116 = Z19*Z21

-> (201) Z117 = Z19*Z22

-> (202) Z118 = OC'*U3

-> (203) Z119 = OC'*U1

-> (204) Z120 = -Z6*Z118 - Z7*Z119

-> (205) Z121 = Z6*Z119 - Z7*Z118

-> (206) Z122 = Z3*Z7

-> (207) Z123 = Z3*Z6

-> (208) Z124 = Z3*Z120

-> (209) Z125 = Z3*OC''

-> (210) Z126 = 0.08333333*Z23*Z27 - 0.08333333*Z25*Z26

-> (211) Z127 = Z24*Z26

-> (212) Z128 = Z24*Z27

-> (213) Z129 = OD'*U3

-> (214) Z130 = OD'*U2

-> (215) Z131 = Z8*Z129 - Z9*Z130

-> (216) Z132 = -Z8*Z130 - Z9*Z129

-> (217) Z133 = Z3*OD''

-> (218) Z134 = Z3*Z8

-> (219) Z135 = Z3*Z9

-> (220) Z136 = Z3*Z131

-> (221) Z137 = 0.08333333*Z28*Z32 - 0.08333333*Z29*Z31

-> (222) Z138 = Z30*Z31

-> (223) Z139 = Z30*Z32

-> (224) Z140 = OE'*U3

-> (225) Z141 = OE'*U2

-> (226) Z142 = Z10*Z140 - Z11*Z141

-> (227) Z143 = -Z10*Z141 - Z11*Z140

-> (228) Z144 = Z3*OE''

-> (229) Z145 = Z3*Z10

-> (230) Z146 = Z3*Z11

-> (231) Z147 = Z3*Z142

-> (232) Z148 = 0.08333333*Z33*Z37 - 0.08333333*Z34*Z36

-> (233) Z149 = Z35*Z36

-> (234) Z150 = Z35*Z37

-> (235) Z151 = 0.5*M*R^2 + 0.0625*Z1 + 0.1666667*Z3 + 4*M*Z14^2 + 16*Ma*Z13^2

-> (236) Z152 = Z151 + 0.08333333*Z4*Z112 + 0.08333333*Z6*Z123

-> (237) Z153 = -0.08333333*Z4*Z111 - 0.08333333*Z6*Z122

-> (238) Z154 = 0.08333333*Z133 + 0.08333333*Z144 + Z106 + Z5*Z115 + Z7*Z126 +

0.08333333*Z4*Z113 + 0.08333333*Z6*Z124 + 0.5*Z49*Z98 + 0.5*M*(R*Z96-2

*Z14*Z95) - 0.08333333*Z139 - 0.08333333*Z150 - 16*Z48*Z89 - 0.5*Z49* Z101 - 0.08333333*Z4*Z117 - 0.08333333*Z6*Z128 - 0.5*M*(R*Z93+2*Z14* Z92)

-> (239) Z155 = Z151 + 0.08333333*Z8*Z134 + 0.08333333*Z10*Z145

-> (240) Z156 = 0.08333333*Z8*Z135 + 0.08333333*Z10*Z146

-> (241) Z157 = 0.08333333*Z114 + 0.08333333*Z125 + Z105 + 0.08333333*Z116 +

0.08333333*Z127 + 0.08333333*Z8*Z136 + 0.08333333*Z8*Z138 + 0.08333333

*Z10*Z147 + 0.08333333*Z10*Z149 + 0.5*Z49*Z97 + 16*Z48*Z90 + 0.5*M*(R* Z100+2*Z14*Z99) - Z9*Z137 - Z11*Z148 - 0.5*Z49*Z94 - 0.5*M*(R*Z103-2* Z14*Z102)

-> (242) Z158 = M*R^2 + 0.125*Z2

-> (243) Z159 = Z158 + 0.08333333*Z5*Z111 + 0.08333333*Z7*Z122 + 0.08333333*Z9*

Z135 + 0.08333333*Z11*Z146

-> (244) Z160 = -0.08333333*Z5*Z112 - 0.08333333*Z7*Z123

-> (245) Z161 = 0.08333333*Z9*Z134 + 0.08333333*Z11*Z145

-> (246) Z162 = Z104 + Z4*Z115 + Z6*Z126 + Z8*Z137 + Z10*Z148 + 0.08333333*Z5* Z117 + 0.08333333*Z7*Z128 + 0.08333333*Z9*Z136 + 0.08333333*Z9*Z138 +

0.08333333*Z11*Z147 + 0.08333333*Z11*Z149 + 0.25*Z50*Z98 + 0.25*Z50* Z101 - 0.25*Z50*Z94 - 0.25*Z50*Z97 - 0.08333333*Z5*Z113 - 0.08333333* Z7*Z124

-> (247) ZERO¹ = -Z154 - Z152*U1' - Z153*U3'

-> (248) ZERO² = -Z157 - Z155*U2' - Z156*U3'

-> (249) ZERO³ = -Z162 - Z159*U3' - Z160*U1' - Z161*U2'

(250) OUTPUT T,U1,U2,U3,THETA

(251) CODE DYNAMICS() THESIS1.M, SUBS

B. Matlab code for simulation of nutation angle

function THESIS1

SolveOrdinaryDifferentialEquations

Abbildung in dieser Leseprobe nicht enthalten

% Unit conversions

Pi = 3.141592653589793;

DEGtoRAD = Pi/180.0; RADtoDEG = 180.0/Pi;

% Reserve space and initialize matrices

z = zeros(162,1);

COEF = zeros(3,3); RHS = zeros(1,3);

% Evaluate constants

z(2) = Ma*R^2;

z(1) = Ma*(R^2+1.333333333333333*H^2); z(3) = M*L^2;

z(12) = Ma + 4*M; z(13) = M*P/z(12); z(14) = P - 4*z(13); z(48) = Ma*z(13); z(49) = M*z(14); z(50) = M*R;

z(151) = 0.5*M*R^2 + 0.0625*z(1) + 0.1666666666666667*z(3) + 4*M*z(14)^2 + 16*Ma*z(13)^2; z(158) = M*R^2 + 0.125*z(2);

% Set the initial values of the states

VAR(1) = U1;

VAR(2) = U2; VAR(3) = U3;

Abbildung in dieser Leseprobe nicht enthalten

FileIdentifier = fopen('THESIS1.1', 'wt'); if( FileIdentifier == -1 ) error('Error: unable to open file THESIS1.1'); end

fprintf( 1, '%% T U1 U2 U3

fprintf( 1, '%% (UNITS) (UNITS) (UNITS) fprintf(FileIdentifier, '%% FILE: THESIS1.1\n%%\n' );

fprintf(FileIdentifier, '%% T U1 U2 U3

fprintf(FileIdentifier, '%% (UNITS) (UNITS) (UNITS)

THETA\n' );

(UNITS) (UNITS)\n\n' );

THETA\n' ); (UNITS) (UNITS)\n\n' );

Abbildung in dieser Leseprobe nicht enthalten

global H L M Ma P R;

global OB OC OD OE OBp OCp ODp OEp OBpp OCpp ODpp OEpp; global U1 U2 U3;

global THETA U1p U2p U3p;

global DEGtoRAD RADtoDEG z COEF RHS SolutionToAxEqualsB; global TINITIAL TFINAL INTEGSTP PRINTINT ABSERR RELERR;

OpenOutputFilesAndWriteHeadings

VAR = ReadUserInput;

OdeMatlabOptions = odeset( 'RelTol',RELERR, 'AbsTol',ABSERR, 'MaxStep',INTEGSTP ); T = TINITIAL;

PrintCounter = 0;

mdlDerivatives(T,VAR,0); while 1,

if( TFINAL>=TINITIAL & T+0.01*INTEGSTP>=TFINAL ) PrintCounter = -1; end if( TFINAL<=TINITIAL & T+0.01*INTEGSTP<=TFINAL ) PrintCounter = -1; end if( PrintCounter <= 0.01 ),

mdlOutputs(T,VAR,0);

if( PrintCounter == -1 ) break; end PrintCounter = PRINTINT;

end

[TimeOdeArray,VarOdeArray] = ode45( @mdlDerivatives, [T T+INTEGSTP], VAR, OdeMatlabOptions, 0 );

TimeAtEndOfArray = TimeOdeArray( length(TimeOdeArray) );

if( abs(TimeAtEndOfArray - (T+INTEGSTP) ) >= abs(0.001*INTEGSTP) ) warning('numerical integration failed'); break; end

T = TimeAtEndOfArray;

VAR = VarOdeArray( length(TimeOdeArray), : ); PrintCounter = PrintCounter - 1;

end

mdlTerminate(T,VAR,0);

Abbildung in dieser Leseprobe nicht enthalten

global H L M Ma P R;

global OB OC OD OE OBp OCp ODp OEp OBpp OCpp ODpp OEpp; global U1 U2 U3;

global THETA U1p U2p U3p;

global DEGtoRAD RADtoDEG z COEF RHS SolutionToAxEqualsB; global TINITIAL TFINAL INTEGSTP PRINTINT ABSERR RELERR;

% Update variables after integration step

U1 = VAR(1);

U2 = VAR(2);

U3 = VAR(3);

z(15) = z(1)*U1;

z(16) = z(1)*U2;

z(17) = z(2)*U3;

z(89) = z(13)*U2*U3; z(90) = z(13)*U1*U3;

z(92) = 0.5*R*U1^2 + 0.5*U3*(R*U3+2*z(14)*U2); z(93) = -z(14)*U1^2 - 0.5*U2*(R*U3+2*z(14)*U2); z(94) = U1*(R*U2-2*z(14)*U3);

z(95) = -0.5*R*U1^2 - 0.5*U3*(R*U3-2*z(14)*U2); z(96) = 0.5*U2*(R*U3-2*z(14)*U2) - z(14)*U1^2; z(97) = U1*(R*U2+2*z(14)*U3);

z(98) = U2*(R*U1-2*z(14)*U3);

z(99) = 0.5*R*U2^2 + 0.5*U3*(R*U3+2*z(14)*U1); z(100) = -z(14)*U2^2 - 0.5*U1*(R*U3+2*z(14)*U1); z(101) = U2*(R*U1+2*z(14)*U3);

z(102) = -0.5*R*U2^2 - 0.5*U3*(R*U3-2*z(14)*U1); z(103) = 0.5*U1*(R*U3-2*z(14)*U1) - z(14)*U2^2; z(104) = 0.0625*U1*z(16) - 0.0625*U2*z(15);

z(105) = 0.0625*U3*z(15) - 0.125*U1*z(17); z(106) = 0.125*U2*z(17) - 0.0625*U3*z(16);

% Quantities to be specified

X = 0.5;

OMEGA = 0.25;

OB = X*sin(OMEGA*T);

OBp = OMEGA*X*cos(OMEGA*T);

OBpp = -X*OMEGA^2*sin(OMEGA*T);

OC = 0;

OCp = 0; OCpp = 0;

OD = 0; ODp = 0; ODpp = 0;

OE = 0;

OEp = 0; OEpp = 0;

% OB = X*sin(OMEGA*T);

% OBp = OMEGA*X*cos(OMEGA*T);

% OBpp = -X*OMEGA^2*sin(OMEGA*T); %

% OC = X*sin(OMEGA*T);

% OCp = OMEGA*X*cos(OMEGA*T);

% OCpp = -X*OMEGA^2*sin(OMEGA*T); %

% OD = 0; % ODp = 0; % ODpp = 0; %

% OE = 0; % OEp = 0; % OEpp = 0;

% OB = X*sin(OMEGA*T);

% OBp = OMEGA*X*cos(OMEGA*T);

% OBpp = -X*OMEGA^2*sin(OMEGA*T); %

% OC = X*sin(OMEGA*T);

% OCp = OMEGA*X*cos(OMEGA*T);

% OCpp = -X*OMEGA^2*sin(OMEGA*T); %

% OD = X*sin(OMEGA*T);

% ODp = OMEGA*X*cos(OMEGA*T);

% ODpp = -X*OMEGA^2*sin(OMEGA*T); %

% OE = X*sin(OMEGA*T);

% OEp = OMEGA*X*cos(OMEGA*T);

% OEpp = -X*OMEGA^2*sin(OMEGA*T);

z(4) = cos(OB);

z(5) = sin(OB); z(6) = cos(OC); z(7) = sin(OC); z(8) = cos(OD); z(9) = sin(OD); z(10) = cos(OE); z(11) = sin(OE);

z(18) = z(4)*U1 - z(5)*U3; z(19) = z(4)*U3 + z(5)*U1; z(20) = OBp + U2;

z(21) = z(3)*z(18); z(22) = z(3)*z(20);

z(23) = z(6)*U1 - z(7)*U3; z(24) = z(6)*U3 + z(7)*U1; z(25) = OCp + U2;

z(26) = z(3)*z(23); z(27) = z(3)*z(25); z(28) = ODp + U1;

z(29) = z(8)*U2 + z(9)*U3; z(30) = z(8)*U3 - z(9)*U2; z(31) = z(3)*z(28);

z(32) = z(3)*z(29); z(33) = OEp + U1;

z(34) = z(10)*U2 + z(11)*U3; z(35) = z(10)*U3 - z(11)*U2; z(36) = z(3)*z(33);

z(37) = z(3)*z(34); z(107) = OBp*U3; z(108) = OBp*U1;

z(109) = -z(4)*z(107) - z(5)*z(108); z(111) = z(3)*z(5);

z(112) = z(3)*z(4); z(113) = z(3)*z(109); z(114) = z(3)*OBpp;

z(115) = 0.08333333333333333*z(18)*z(22) - 0.08333333333333333*z(20)*z(21); z(116) = z(19)*z(21);

z(117) = z(19)*z(22); z(118) = OCp*U3; z(119) = OCp*U1;

z(120) = -z(6)*z(118) - z(7)*z(119); z(122) = z(3)*z(7);

z(123) = z(3)*z(6); z(124) = z(3)*z(120); z(125) = z(3)*OCpp;

z(126) = 0.08333333333333333*z(23)*z(27) - 0.08333333333333333*z(25)*z(26); z(127) = z(24)*z(26);

z(128) = z(24)*z(27); z(129) = ODp*U3; z(130) = ODp*U2;

z(131) = z(8)*z(129) - z(9)*z(130); z(133) = z(3)*ODpp;

z(134) = z(3)*z(8); z(135) = z(3)*z(9); z(136) = z(3)*z(131);

z(137) = 0.08333333333333333*z(28)*z(32) - 0.08333333333333333*z(29)*z(31); z(138) = z(30)*z(31);

z(139) = z(30)*z(32); z(140) = OEp*U3; z(141) = OEp*U2;

z(142) = z(10)*z(140) - z(11)*z(141); z(144) = z(3)*OEpp;

z(145) = z(3)*z(10); z(146) = z(3)*z(11); z(147) = z(3)*z(142);

z(148) = 0.08333333333333333*z(33)*z(37) - 0.08333333333333333*z(34)*z(36); z(149) = z(35)*z(36);

z(150) = z(35)*z(37);

z(152) = z(151) + 0.08333333333333333*z(4)*z(112) + 0.08333333333333333*z(6)*z(123);

z(153) = -0.08333333333333333*z(4)*z(111) - 0.08333333333333333*z(6)*z(122);

z(154) = 0.08333333333333333*z(133) + 0.08333333333333333*z(144) + z(106) + z(5)*z(115) + z(7)*z(126) + 0.08333333333333333*z(4)*z(113) + 0.08333333333333333*z(6)*z(124) +

0.5*z(49)*z(98) + 0.5*M*(R*z(96)-2*z(14)*z(95)) - 0.08333333333333333*z(139) -

0.08333333333333333*z(150) - 16*z(48)*z(89) - 0.5*z(49)*z(101) - 0.08333333333333333*z(4)*z(117) -

0.08333333333333333*z(6)*z(128) - 0.5*M*(R*z(93)+2*z(14)*z(92));

z(155) = z(151) + 0.08333333333333333*z(8)*z(134) + 0.08333333333333333*z(10)*z(145); z(156) = 0.08333333333333333*z(8)*z(135) + 0.08333333333333333*z(10)*z(146); z(157) = 0.08333333333333333*z(114) + 0.08333333333333333*z(125) + z(105) +

0.08333333333333333*z(116) + 0.08333333333333333*z(127) + 0.08333333333333333*z(8)*z(136) +

0.08333333333333333*z(8)*z(138) + 0.08333333333333333*z(10)*z(147) +

0.08333333333333333*z(10)*z(149) + 0.5*z(49)*z(97) + 16*z(48)*z(90) +

0.5*M*(R*z(100)+2*z(14)*z(99)) - z(9)*z(137) - z(11)*z(148) - 0.5*z(49)*z(94) - 0.5*M*(R*z(103)- 2*z(14)*z(102));

z(159) = z(158) + 0.08333333333333333*z(5)*z(111) + 0.08333333333333333*z(7)*z(122) +

0.08333333333333333*z(9)*z(135) + 0.08333333333333333*z(11)*z(146);

z(160) = -0.08333333333333333*z(5)*z(112) - 0.08333333333333333*z(7)*z(123); z(161) = 0.08333333333333333*z(9)*z(134) + 0.08333333333333333*z(11)*z(145); z(162) = z(104) + z(4)*z(115) + z(6)*z(126) + z(8)*z(137) + z(10)*z(148) +

0.08333333333333333*z(5)*z(117) + 0.08333333333333333*z(7)*z(128) +

0.08333333333333333*z(9)*z(136) + 0.08333333333333333*z(9)*z(138) +

0.08333333333333333*z(11)*z(147) + 0.08333333333333333*z(11)*z(149) + 0.25*z(50)*z(98) +

0.25*z(50)*z(101) - 0.25*z(50)*z(94) - 0.25*z(50)*z(97) - 0.08333333333333333*z(5)*z(113) -

0.08333333333333333*z(7)*z(124);

COEF(1,1) = -z(152);

COEF(1,2) = 0;

COEF(1,3) = -z(153); COEF(2,1) = 0;

COEF(2,2) = -z(155); COEF(2,3) = -z(156); COEF(3,1) = -z(160); COEF(3,2) = -z(161); COEF(3,3) = -z(159);

RHS(1) = z(154); RHS(2) = z(157); RHS(3) = z(162);

SolutionToAxEqualsB = COEF\RHS';

% Update variables after uncoupling equations

U1p = SolutionToAxEqualsB(1);

U2p = SolutionToAxEqualsB(2); U3p = SolutionToAxEqualsB(3);

% Update derivative array prior to integration step

VARp(1) = U1p;

VARp(2) = U2p; VARp(3) = U3p;

sys = VARp';

Abbildung in dieser Leseprobe nicht enthalten

global H L M Ma P R;

global OB OC OD OE OBp OCp ODp OEp OBpp OCpp ODpp OEpp; global U1 U2 U3;

global THETA U1p U2p U3p;

global DEGtoRAD RADtoDEG z COEF RHS SolutionToAxEqualsB; global TINITIAL TFINAL INTEGSTP PRINTINT ABSERR RELERR;

% Evaluate output quantities

z(40) = z(14)*U2 + 0.5*R*U3; z(47) = 0.5*R*U3 - z(14)*U1; z(43) = z(14)*U2 - 0.5*R*U3; z(45) = -z(14)*U1 - 0.5*R*U3;

z(53) = 0.125*z(17) + 0.5*z(50)*z(40) + 0.5*z(50)*z(47) - 0.5*z(50)*z(43) - 0.5*z(50)*z(45); z(42) = R*U1;

z(41) = z(14)*U1; z(39) = z(13)*U1;

z(51) = 0.0625*z(15) + 0.5*z(50)*z(42) + 2*z(49)*z(41) + 16*z(48)*z(39) - z(49)*z(45) - z(49)*z(47); z(46) = R*U2;

z(44) = z(14)*U2; z(38) = z(13)*U2;

z(52) = 0.0625*z(16) + z(49)*z(40) + z(49)*z(43) + 0.5*z(50)*z(46) + 2*z(49)*z(44) + 16*z(48)*z(38); z(54) = z(4)*z(6) + z(5)*z(7);

z(73) = z(8)*z(10) + z(9)*z(11); z(57) = z(5)*z(9);

z(61) = z(5)*z(11); z(65) = z(7)*z(9); z(69) = z(7)*z(11);

z(76) = z(51)^2 + z(52)^2 + z(53)^2 + 0.006944444444444444*z(21)^2 +

0.006944444444444444*z(22)^2 + 0.006944444444444444*z(26)^2 + 0.006944444444444444*z(27)^2 +

0.006944444444444444*z(31)^2 + 0.006944444444444444*z(32)^2 + 0.006944444444444444*z(36)^2 +

0.006944444444444444*z(37)^2 + 0.01388888888888889*z(22)*z(27) +

0.01388888888888889*z(31)*z(36) + 0.1666666666666667*z(22)*z(52) +

0.1666666666666667*z(27)*z(52) + 0.1666666666666667*z(31)*z(51) +

0.1666666666666667*z(36)*z(51) + 0.01388888888888889*z(4)*z(21)*z(31) +

0.01388888888888889*z(4)*z(21)*z(36) + 0.01388888888888889*z(6)*z(26)*z(31) +

0.01388888888888889*z(6)*z(26)*z(36) + 0.01388888888888889*z(8)*z(22)*z(32) +

0.01388888888888889*z(8)*z(27)*z(32) + 0.01388888888888889*z(10)*z(22)*z(37) +

0.01388888888888889*z(10)*z(27)*z(37) + 0.01388888888888889*z(54)*z(21)*z(26) +

0.01388888888888889*z(73)*z(32)*z(37) + 0.1666666666666667*z(4)*z(21)*z(51) +

0.1666666666666667*z(6)*z(26)*z(51) + 0.1666666666666667*z(8)*z(32)*z(52) +

0.1666666666666667*z(9)*z(32)*z(53) + 0.1666666666666667*z(10)*z(37)*z(52) +

0.1666666666666667*z(11)*z(37)*z(53) - 0.1666666666666667*z(5)*z(21)*z(53) -

0.1666666666666667*z(7)*z(26)*z(53) - 0.01388888888888889*z(57)*z(21)*z(32) -

0.01388888888888889*z(61)*z(21)*z(37) - 0.01388888888888889*z(65)*z(26)*z(32) -

0.01388888888888889*z(69)*z(26)*z(37);

z(77) = z(76)^0.5;

z(80) = z(53)/z(77); z(86) = z(32)/z(77); z(88) = z(37)/z(77); z(81) = z(21)/z(77); z(83) = z(26)/z(77);

THETA = acos(z(80)+0.08333333333333333*z(9)*z(86)+0.08333333333333333*z(11)*z(88)-

0.08333333333333333*z(5)*z(81)-0.08333333333333333*z(7)*z(83));

Output(1)=T; Output(2)=U1; Output(3)=U2; Output(4)=U3; Output(5)=THETA; FileIdentifier = fopen('all');

WriteOutput( 1, Output(1:5) );

WriteOutput( FileIdentifier(1), Output(1:5) );

%=========================================================================== function WriteOutput( fileIdentifier, Output )

numberOfOutputQuantities = length( Output ); if numberOfOutputQuantities > 0,

for i=1:numberOfOutputQuantities,

fprintf( fileIdentifier, ' %- 14.6E', Output(i) ); end

fprintf( fileIdentifier, '\n' ); end

Abbildung in dieser Leseprobe nicht enthalten

FileIdentifier = fopen('all');

fclose( FileIdentifier(1) );

fprintf( 1, '\n Output is in the file THESIS1.1\n' );

fprintf( 1, '\n To load and plot columns 1 and 2 with a solid line and columns 1 and 3 with a dashed line, enter:\n' );

fprintf( 1, ' someName = load( ''THESIS1.1'' );\n' );

fprintf( 1, ' plot( someName(:,1), someName(:,2), ''-'', someName(:,1), someName(:,3), ''--'' )\n\n' ); sys = [];

Abbildung in dieser Leseprobe nicht enthalten

function [sys,x0,str,ts] = Sfunction(t,x,u,flag) switch flag,

case 0, [sys,x0,str,ts] = mdlInitializeSizes; and sample times ts

case 1, sys = mdlDerivatives(t,x,u); in sys

case 2, sys = mdlUpdate(t,x,u); case 3, sys = mdlOutputs(t,x,u);

% Initialization of sys, initial state x0, state ordering string str, % Calculate the derivatives of continuous states and store them

% Update discrete states x(n+1) in sys % Calculate outputs in sys

case 4, sys = mdlGetTimeOfNextVarHit(t,x,u); % Return next sample time for variable-step in sys

case 9, sys = mdlTerminate(t,x,u); % Perform end of simulation tasks and set sys=[]

otherwise error(['Unhandled flag = ',num2str(flag)]); end

Abbildung in dieser Leseprobe nicht enthalten

sizes = simsizes; % Call simsizes to create a sizes structure

sizes.NumContStates = 3; sizes.NumDiscStates = 0; sizes.NumOutputs = 5;

sizes.NumInputs = 0; sizes.DirFeedthrough = 1;

% sys(1) is the number of continuous states % sys(2) is the number of discrete states % sys(3) is the number of outputs

% sys(4) is the number of inputs

% sys(6) is 1, and allows for the output to be a function of the input

sizes.NumSampleTimes = 1; % sys(7) is the number of samples times (the number of rows in ts)

sys = simsizes(sizes); % Convert it to a sizes array

stateOrderingStrings = [];

timeSampling = [0 0]; % m-by-2 matrix containing the sample times OpenOutputFilesAndWriteHeadings

VAR = ReadUserInput

C. Data used for plots in chapter 4

Frequency of payload oscillation ‘ Ω ’

Abbildung in dieser Leseprobe nicht enthalten

Amplitude of payload oscillation ‘ X ’

Abbildung in dieser Leseprobe nicht enthalten

Payload Location ‘ P ’

Abbildung in dieser Leseprobe nicht enthalten

PAYLOAD INERTIA

Abbildung in dieser Leseprobe nicht enthalten