Äîêóìåíò âçÿò èç êýøà ïîèñêîâîé ìàøèíû. Àäðåñ îðèãèíàëüíîãî äîêóìåíòà : http://www.ssau.ru/files/education/uch_posob/Attitude%20Dynamics-Balakin%20VL.pdf Äàòà èçìåíåíèÿ: Tue Dec 9 10:44:29 2014 Äàòà èíäåêñèðîâàíèÿ: Mon Apr 11 02:51:09 2016 Êîäèðîâêà: Ïîèñêîâûå ñëîâà: m 63

The Ministry of Education and Science of the Russian Federation Samara State Aerospace University (National Research University) Flight Dynamics and Control Systems Department

V. L. BALAKIN, A. V. DOROSHIN

ATTITUDE DYNAMICS AND CONTROL OF A DUALSPIN-SPACECRAFT, GYROSTAT-SATELLITES AND NANOSATELLITES WITH MULTYROTOR SYSTEMS

Electronic Textbook

SAMARA 2012




Authors: Balakin Victor L., Doroshin Anton V.
Translated by: Doroshin Anton V.

Balakin, V. L. Attitude Dynamics and Control of a Dual-Spin-Spacecraft, GyrostatSatellites and Nanosatellites with Multyrotor Systems = [Electronic resource] : Electronic Textbook / V. L. Balakin, A. V. Doroshin; The Ministry of Education and Science of the Russian Federation, Samara State Aerospace University. - Electronic text and graphic data (2.17 Mb). - Samara, 2012. - 1 CD-ROM.

This textbook gives a description of attitude dynamics and control of dual -spinspacecraft and gyrostat-satellites (also spacecraft with multyrotor systems). The textbook is a part of the Masters programmes for education directions 161100.68 «Motion Control Systems and Navigation», 160400.68 «Space-Rocket complexes and Cosmonautics», 151600.68 «Applied Mechanics». Material prepared at the Department of Flight Dynamics and Control Systems.

© Samara State Aerospace University, 2012


CONTENTS
Introduction ............................................................................................................4 Part I. Modeling of chaotic motion of gyrostats in resistant environment on the base of dynamical systems with strange attractors.......................................................................................5 Part II. Analysis of attitude motion evolutions of variable mass gyrostats and coaxial rigid bodies system..................................................................................................................22 Part III. Attitude control of spider-type multiple-rotor rigid bodies systems.................................52 Part IV. Example of MAPLE-program for attitude dynamics of the spacecraft modeling..................65

3


INTRODUCTION
Problem of rigid bodies motion and its practical engineering applications such as gyroscopes, gyrostats and dual-spin-spacecraft attitude dynamics and control are very important for modern science. Despite classical analytical research results and exact solutions this problem is still far from complete due to the existence of chaos phenomena. Among the basic directions of modern research within the framework of the indicated problem it is possible to highlight the following points: deriving exact and approximated analytical and asymptotic solutions, investigation into stability of motion, the analysis of motion under an influence of external regular and stochastic disturbance, research into dynamic chaos and study of non-autonomous systems with variable parameters.

This book gives a description of above mentioned scientific problems on the base of research results corresponded to the following scientific works (references are indicated in the original language): Part I: 1. A.V. Doroshin, Modeling of chaotic motion of gyrostats in resistant environment on the base of dynamical systems with strange attractors. Communications in Nonlinear Science and Numerical Simulation, Volume 16, Issue 8 (2011) 3188­3202. Part II: 1. A.V. Doroshin, Analysis of attitude motion evolutions of variable mass gyrostats and coaxial rigid bodies system, International Journal of Non-Linear Mechanics, Volume 45, Issue 2 (2010) 193­205. 2. .., .., .., .., .., .., .. - // . : , 1, 2003 . 3. A.V. Doroshin, Synthesis of attitude motion of variable mass coaxial bodies, WSEAS Transactions on Systems and Control, Issue 1, Volume 3 (2008) 50-61. 4. .., .., // . / , 2012.-440. //.54-55 Part III: 5. A.V. Doroshin, Attitude Control of Spider-type Multiple-rotor Rigid Bodies Systems. Proceedings of the World Congress on Engineering 2009, London, U.K. Vol II, pp.15441549. .. : . . . -.-:

4


PART I. MODELING OF CHAOTIC MOTION OF GYROSTATS IN RESISTANT ENVIRONMENT ON THE BASE OF DYNAMICAL SYSTEMS WITH 1 STRANGE ATTRACTORS
1 Introduction Problem of rigid bodies motion and its practical engineering applications such as gyroscopes, gyrostats and dual-spin-spacecraft are very important for modern science. Despite classical analytical research results and exact solutions this problem is still far from complete due to the existence of chaos phenomena [1-13]. Among the basic directions of modern research within the framework of the indicated problem it is possible to highlight the following points: deriving exact and approximated analytical and asymptotic solutions, investigation into stability of motion, the analysis of motion under an influence of external regular and stochastic disturbance, research into dynamic chaos and study of non-autonomous systems with variable parameters. Recently, chaotic dynamic has becomes one of the major part of nonlinear scien ce. Applications of dynamical systems with chaotic behavior and strange attractors are seen in many areas of science, including space-rocket systems [7-12]. E. N. Lorenz and O. E. RÆssler systems [1, 2] represent classical dynamical systems with strange attractors. R. B. Leipnik and T. A. Newton [3] found two strange attractors in rigid body motion. Since Leipnik and Newton 's work, the chaotic dynamics of rigid body motion investigates in many works. J. C. Sprott [4, 5] examined 19 systems of three-dimensional autonomous ordinary differential equations with strange attractors; also critical points, Lyapunov exponents and fractional dimensions of systems were found. Work [7] contains the analysis of chaotic behavior of a spacecraft with periodic time -dependent moments of inertia during its free motion. The equations of variable mass coaxial bodies system were developed in papers [10] where also the attitude motion of coaxial bodies system and double rotation spacecraft with time-dependent moments of inertia were analyzed on the base of special method of phase trajectory curvature analysis. The results [7-12] can be used for the analysis of attitude motion of a gyrostat-satellites and dual-spin spacecraft including motion with an active solid-propellant rocket engine. In this paragraph more attention is focused on chaotic attractors in phase space of angular velocity of gyrostat and on perturbed gyrostat motion in resistant environment with energy dissipation/excitation. Conditions of correspondence of mathematical models of gyrostats in resistant environment and dynamical systems with strange attractors (Lorenz, RÆssler, Newton-Leipnik and Sprott) are defined. To confirm the system chaotic behavior numerical computer simulations are used. These simulations are performed by means of numerical integration of the equations of motion with the help of several numerical tools: time history of phase coordinates, gyrostat longitudinal axis vector hodograph, PoincarÈ map, fast Fourier transform power spectrum. This characterizes the dynamical behavior of the gyrostat in resistant environment as regular or chaotic.
- Material corresponds to preprint version of article: A.V. Doroshin, Modeling of chaotic motion of gyrostats in resistant environment on the base of dynamical systems with strange attractors. Communications in Nonlinear Science and Numerical Simulation, Volume 16, Issue 8 (2011) 3188­3202.
1

5


2 Mathematical model Let us consider a gyrostat attitude motion about fixed point in resistant environment with energy dissipation/excitation (fig.1). Assume resistant environment effect corresponding to action of external forces moments that are constant M
e const



, linear M

e lin

velocity projections onto body frame axes x1x2x3 p, q, r





and nonlinear M

e quad



T





in main body angular

.

Fig.1 ­ Inertial







and gyrostat main body



x1 x2 x2 frames
(1)
3

The motion equations follow from angular moment's law: K + â K + R = Meonst + Mlein + Meuad c q where

K = I ; R = R1 , R2 , R3 ; M
T

e const

= d1 , d 2 , d



T

e ; M lin = A ;

M

= B p 2 , q 2 , r 2 ; A = aij ; B bij ; aij const; bij const; Ri const ; di const ; i, j 1..3
e quad

T

(2)

K ­ angular moment of gyrostat main body with "frozen" internal rotor; I ­ inertia tensor of main body with "frozen" internal rotor; R ­ constant angular moment of relative rotor motion (in body frame); A, B ­ constant matrixes. Matrix structure of external forces moments (2) can describe an action of viscous drag, V0 of main body with roughened hydro(aero)dynamic lift, nonuniform lift and friction in fluid flow surface and propeller elements. Assume coincidence of gyrostat center of mass, rotor center of mass and fixed point. Also let us I diag A, B, C consider case of spherical inertia tensor of rotor and gyrostat general inertia tensor In this case scalar form of eq. (1) can be write as follows 2 2 2 Ap ( B C )rq a11 p (a12 R3 )q (a13 R2 )r b11 p b12 q b13r d1 2 2 2 (3) Bq (C A) pr a22 q (a23 R1 )r (a21 R3 ) p b21 p b22q b23r d 2 Cr ( A B)qp a r (a R ) p (a R )q b p 2 b q 2 b r 2 d 33 31 2 32 1 31 32 33 3 Dynamical system (3) is supplemented with kinematical system for Euler type angles
6




about x1 about x2 about x3



:

p sin q cos ;

1 p cos q sin ; cos

sin r p cos q sin . cos In considered case gyrostat kinetic energy takes on form: 1 1 T Ap 2 Bq 2 Cr 2 pR1 qR2 rR3 R12 R22 R32 2 2J where J is axial inertia moment of the rotor.
3 Links between gyrostat chaotic motion and strange attractors

(4)



(5)

It is well known fact that unpredictable chaotic long-term solutions can exist for simple nonlinear deterministic systems. The study of nonlinear dynamics has brought new excitement to one of the oldest fields of science and, certainly, mechanics. So, many papers and, for example, works [3, 11, 12] describe chaotic motion of rigid body and gyrostats as modes corresponded to st range attractors in phase space. The paper [4] also contains several interesting and important chaotic dynamical systems with strange attractors. In this paragraph we will find conditions of reduction of gyrostats motion equations (3) to Lorenz, RÆssler, Newton-Leipnik and Sprott dynamical systems. General form of indicated dynamical systems of three autonomous first-order ordinary differential equations (ODE) can be write as:

x f x x, y, z ; y f

y



x, y, z ; z f z x, y, z



(6)

The system (6) has strange attractors in man y cases including classical dynamical systems, which presented in table 1 [4]. Cases A-S correspond to Sprott systems [4, 5], and LOR, ROS, NL ­ to Lorenz, RÆssler, Newton-Leipnik systems. It is possible to write condition of equivalence of dynamical systems (3) and A-NL (tabl.1), where variables change take place p x, q y, r z . First of all we take notice about signature (+/-) in table 1. Signature "+" means possibility of reduction of systems A-NL immediately to system (3): it implies definition of corresponded components values of vectors R, M eo n st and matrix A, B, I . Signature "-" means unrealizability of c





this reduction without presence of additional special control torque of gyroscopic type in right parts of systems (1) and (3): (7) G qr , pr , pq ; G gij ; i, j 1..3 This artificial forces moment (7) can be formed with the help of special technical actuators and thrusters.

M

gyro control

T

Now we can present following conditions of reductions of A-NL systems to system (3), and vice versa. These conditions establish connections between external and internal parameters (mass -inertia, gyrostat rotor angular moment, roughened surface and propeller elements properties, friction in fluid flow etc.).
7


Case A B C D E F G H I J K L M N O P Q R S LOR ROS NL LOR-case conditions:

fx y yz yz -y yz y+z 2x/5+z -y+z2 -y/5 2z xy-z y+3.9z -z -2y y 2.7y+z -z 0.9-y -x-4y -s(x-y) s=10 -y-z -kx+y+wyz k=0.4, w=10

fy -x+yz x-y x-y x+z x2-y -x+y/2 xz-y x+y/2 x+z -2y+z x-y 0.9x2-y -x2-y x+z2 x-z -x+y2 x-y 0.4+z x+z2 -y+wx-xz w=28 x+ky k=0.2 -x-my+5xz m=0.4

fz 1-y2 1-xy 1-x2 xz+3y2 1-4x x2-z -x+y x-z x+y2-z -x+y+y2 x+0.3z 1-x 1.7(1+x)+y2 1+y-2x x+xz+2.7y x+y 3.1+y2+0.5z xy-z 1+x -vz+xy v=8/3 v+(x-w)z v=0.2, w=5.7 vz-5xy v=0.175

Table 1 Signature + + + + + + + + + + + + -

A 2 C B B0 ; bij di gij 0 a11 2 B0 s; a12 2 B0 s R3 ; a13 R2 (8.LOR) a21 B0 w R3 ; a22 B0 ; a23 R1 a R ; a R ; a vC 2 32 1 33 31 If we use substitution of coefficient (8.LOR) into system (3), then we obtain classical Lorenz equations. It is need to note that for LOR-gyrostat ((3) with (8.LOR)) main body is dynamically symmetric (B=C) and third inertia moment is twice as large (A=2B).
ROS-case conditions:
A B C ; g32 C ; gij 0; bij d1,2 0; d3 vC i 3, j 2 a11 0; a12 A R3 ; a13 A R2 a21 B R3 ; a22 kB; a23 R1 a R ; a R ; a wC 2 32 1 33 31 For ROS-gyrostat spherical inertia-mass symmetry takes place (A=B=C).

(8.ROS)

8


NL-case conditions:
: A, B, C ; bij di 0; G diag g11 wA B C ; g 22 5B A C ; g33 5C B A a11 kA; a12 A R3 ; a13 R2 a B R ; a mB; a R 3 22 23 1 21 a31 R2 ; a32 R1 ; a33 vC For NL-gyrostat general case of inertia-mass takes place (ABC). A-case conditions: B C A A0 ; g 21 A0 ; b32 A0 ; d 3 A0 a 0; a A R ; a R 11 12 0 3 13 2 a21 B R3 ; a22 0; a23 R1 a31 R2 ; a32 R1 ; a33 0

(8.NL)

(8.A)

where other components of A, B, G, B-case conditions:

e const

equal to zero.


F-case conditions:

A C ; B 2C ; gij 0; bij 0; d3 C a11 0; a12 R3 ; a13 R2 a21 B R3 ; a22 B; a23 R1 a31 R2 ; a32 R1 ; a33 0
(8.B)

A B C ; gij 0; b31 1; di 0 a11 0; a12 A R3 ; a13 A R2 (8.F) a21 B R3 ; a22 B / 2; a23 R1 a R ; a R ; a C 2 32 1 33 31 Other cases conditions can be write by analogy (by the way of equalization of corresponding coefficients of sys. (3) and A-NL). So we can conclude that dynamical s ystems with strange attractors A-NL correspond to gyrostats equation ((3) with conditions (8.A), (8.B), (8.NL)...), which allow chaotic modes of motion.
4 Perturbed motion examination 4.1 Inertia moments perturbation Haw we saw in previous paragraph dynamical system with strange attractors can correspond to system equations of gyrostat motion. Considered gyrostats possessed constant parameters (moments of inertia, relative rotor angular moment component, resistant environment and gyrostat outer surface properties, etc.). Now let us examine perturbed gyrostat motion with a time­dependent moments of inertia, motion of this gyrostat and influence of parameters variability on strange attractor change. It is need to note, that the inertia moment variability can describe small elastic vibrations in gyrostat construction [7]. LOR-gyrostat. Assume following time-dependencies of inertia moments in the case of LOR-gyrostat: B(t ) C (t ) B0 1 sin t ; A(t ) 2B0 1 sin t
9

(9)


where is small nondimensional parameter 0 1 ; other parameters in (8.LOR) are constant. Take into account conditions (8.LOR) and dependencies (9) we can write motion equations:
s x y x 1 sin t 1 (10) wx y 1 3 sin t xz y 1 sin t 1 z 1 sin t 1 3 sin t xy vz In order to examine of perturbed motion several numerical techniques are used. They are based on the numerical integration of the equations of motion (10) by means of a Runge­Kutta algorithm. So, we present perturbed strange attractor (fig.2-a) in phase space {x, y, z}, x(t) time-history (fig.2-b), power spectrum of x(t) fast Fourier transformation (fig.2-c), kinetic energy (5) time-history (fig.2-d), asymptotics of Lyapunov exponents (fig.2-e) and longitudinal axis vector hodograph (fig.2-f). Fig.2 was obtained at = 0.1 and =100 (1/s). Longitudinal axis vector hodograph e t was plotted with the help of numerical integration of

equations (3), (4) and matrix transformation of components of a unit vector of longitudinal z -axis of main body e
x1 x2 x3

0,0,1 into initial frame :
T

e
1

111

e

x1 x2 x3

(11)

cos 0 sin cos sin 0 , 1 0 , 1 sin cos 0 1 0 sin 0 cos 0 cos 0 1 All signs of chaotic motion are shown (fig.2): complexity and irregularity of phase coordinate, broadly distributed power spectrum, positive Lyapunov exponents. Lyapunov exponents for perturbed motion LOR-gyrostat was calculated on the base of Benettin algorithm [14] (with Gram­Schmidt process of orthogonalizaiton) and have following values (with accuracy 10-2): 0 1 0 cos 0 sin 0 sin

0:



1

0.89; 2 0; 3 14.56 ;

0.1 : 1 0.87; 2 0; 3 14.61 ; 0.5 : 1 1.04; 2 0; 3 16.73 ;

0.75 : 1 1.47; 2 0.14; 3 16.71 ;
0.90 : 1 3.66; 2 1.57; 3 13.51 .

10


(a)

(b )

(c)

(d )

(e) (f) Fig.2 ­ Numerical research results for LOR-gyrostat with inertia moment variability ( = 0.1) 11


The Kaplane-Yorke dimension of perturbed strange attractor increase as compared with classical Lorenz attractor: i 1D DKY D j ; D sup i : j 0; D 1 j 1 j 1 (12)
DKY
0

DKY

0.1

DKY

0.5

2.06;

DKY

0.75

2.08;

Calculation of divergence of perturbed system (10) phase flow F f x , f y , f
v 1 s div F s v 11 sin t s v 1 show that the perturbed system is dissipative div F 0 if



DKY

0.90
z

2.15.



T

(1 3 )

In the classical case of Lorenz s 10; v 8 / 3; w 28 from condition (14) follow limitation 41 19 2.16 , which guarantee the system dissipativity at 1 . Consequently, every finite (small) the system phase-space volume will reduce to zero value and every phase trajectory will attract to strange attractor. Comment about application of LOR-case. The Lorenz system, first of all, describes the convective motion of fluid [1]. This system also can be applied to the analysis of dynamos and lasers. In addition it is need to note that LOR-case can, for example, simulate attitude motion of the gyrostat A=mR2/2, B=C=mR2/4) at presence of propeller blades aij 0 and roughness of the body surface

1 v s

1

vs



(14)



R R1 , 0, 0



T



with inertia-mass parameters corresponded to a thin disk-shaped body (like a coin:

aii 0 . This makes it possible to apply the LOR-case investigation results to examination of vehicles special motion modes in resistant environments. Also these results can be used for the description of e gyrostat-spacecraft perturbed attitude motion with feedback control (interpreting the torques M lin as feedback control).
A-gyrostat. Assume following time-dependencies of inertia moments in the case of A-gyrostat: At Bt A0 1 sin t ; C t A0 1 sin t



(15)

Other parameters in (8.A) are constant. For numerical evaluation we take =100 (1/s). Take into account conditions (8.A) and dependencies (15) perturbed motion equations for A gyrostat can be write as follows: 2 sin t y x yz ; 1 sin t 1 sin t 2 sin t 1 (16) yz x ; y xz 1 sin t 1 sin t
z 1 1 y2 . 1 sin t





Lyapunov exponents for perturbed motion of A-gyrostat (with accuracy 10-2):

0:



1

0.01; 2 0; 3 0.01 ;

0.3 : 1 0.03; 2 0; 3 0.03
The Kaplane-Yorke dimension in this case always equals to 3; the system is conservative and phase 3 space volume conservation takes place i 0 . i 1
12


(a)

(b)

(c)

(d)

(e) (f) Fig.3 ­ Numerical research results for A-gyrostat with inertia moment variability ( = 0.3)
13


Integer (not fractional) dimension and presence of positive Lyapunov index means that this system has not strange attractor (like geometry objects with fractional dimension), but gyrostat motion is chaotic (positive -exponent mixes of phase trajectories). Numerical modeling results are presented at figures (fig.3 ­ fig.6). Fig.4-5 contain PoincarÈ sections (z=0) of the system phase space for unperturbed [4] (fig.4) and perturbed (fig.5, 6) cases. It is needed to note, that phase trajectory intersect the plane (z=0) in different region depending on direction of phase point motion along phase trajectory (fig.4-b): 1). Region y , 1 1, corresponds to intersection with direction z 0, z :
2). Region y 1,1 corresponds to intersection with direction z 0, z :

(a)

(b)

(c) (d) Fig.4 ­ PoincarÈ sections (z=0) in unperturbed A-gyrostat case ( = 0) [4]: a ­ general PoincarÈ section; b ­ with intersection direction control z : ; c, d - zoom

14


(a)

(b)

(a) (b) Fig.5 ­ PoincarÈ sections (z=0) in perturbed A-gyrostat ( = 0.3):
a ­ general PoincarÈ section; b ­ with initial condition from depicted rectangle; c, d - zoom

How can we see, perturbation generate heteroclinic loops and corresponding meander tori at the PoincarÈ sections (fig. 5). This circumstance even more complicates the system motion dynamics. Also it is need to note, that time history of kinetic energy T(t) show, on the one hand, gyrostat chaotic motion features and, on the other hand, nonregular characteristics of external environment and internal forces action. Kinetic energy change law imply T dW e dW i W (t ) const

where W(t) is total work of all external ("e") and internal ("i") forces. It corroborates the statement that deterministic chaos in dynamical system (and strange attractor like its geometrical image) can be explained on the base of mechanical description: presence of nonregular influence result in nonregular system behavior. Thus, we shall conclude that kinetic energy T(t) time history is also one of the primary technique for examine of chaotic motion.

15


(a) (b) Fig.6 ­ PoincarÈ sections (z=0) in perturbed A-gyrostat ( = 0.5):
a ­ general PoincarÈ section; b ­ zoom

Comment about application of A-gyrostat Sprott case. The Sprott system for A-gyrostat can be applied, for example, to the analysis of attitude motion of the gyrostat R 0, 0, R3





T



with inertia-

mass parameters of a spherical body (A=B=C), xy-propeller blades a12 a21 A R3 , smooth body surface aii 0 at presence of constant z-spin-up torque (d3=A) and special feedback control (g21=b32=A). This makes it possible to apply the A-case investigation results to examination of gyrostatvehicles special motion modes in resistant environments with feedback control.

4.2 Gyrostat internal rotor angular moment perturbation Let us investigate of gyrostat motion at presence of small harmonic perturbations in relative rotor angular moment R: T R R 1 sin t ; R R1 , R2 , R3 ; Ri const (17) This perturbation can be associated with existence of small harmonic disturbances in electric circuit of internal rotor-engine (simulation of simplest self-induction effects). Corresponding motion equations follow from angular moments law: K + â K + R = Meonst + Mlein + Me uad R (18) c q We conduct examination of perturbed motion on the base of NL-gyrostat. Other type of gyrostat (A-S, LOR, ROS) can be considered by analogy. Take into account conditions (8.NL) and (17) perturbed motion equations for NL-gyrostat will be write as follows: x kx y wyz Pert1 (19) y x my 5 xz Pert2 z vz 5 xy Pert 3 where Perti are components of vector

16


A R2 z R3 y sin t R1 cos t R x R z sin t R cos t (20) Pert 1 2 B 3 R1 y R2 x sin t R3 cos t C Let us note, that perturbation vector (20) will be the same also for other type (A-NL). Case 1. Firstly, consider main case of the NL system with w=10. Numerical research results are present at fig.7 and was obtained at following parameters and initial condition values: A=B=C=1; R1=1; R2=1.5; R3=2; =100; =0.01; x(0)=0.349; y(0)=0.0; z(0)=-0.16. In this case Lyapunov exponents and Kaplane-Yorke dimension for unperturbed and perturbed motion of NL-gyrostat (with accuracy 10-2) are equal:

0:



1

0.14; 2 0; 3 0.76 ; DKY 2.18 ;
.

0.01 : 1 0.12; 2 0.01; 3 0.74; DKY 2.18

Consequently, the system is dissipative (negative sum of all Lyapunov index) and has attractor; the system is chaotic (1>0); the system attractor is strange (fractional DKY). Case 2. Now consider case with w=1; other parameters are the same, like previous case. Numerical research results are present at fig.8. In this case Lyapunov exponents and Kaplane-Yorke dimension (with accuracy 10-2) are equal:

0:



1

0.01; 2 0.10; 3 0.53 ; DKY 1.1;

0.01 : 1 0.01; 2 0.11; 3 0.53; DKY 1.09.
The system also is dissipative, chaotic (1>0) and has strange attractor. But absolute value of positive 1-exponent is small (limiting close to zero with actual accuracy), therefore, trajectory mixing is weak. It allows conclude, that the system is quasichaotic. It is also supported by re gulation trend of time history of phase coordinate, kinetic energy, longitudinal axes hodograph, and by chaotic but degenerating power spectrum (fig.8). Case 3. Finally, let us consider case for w=10, v=0. In this case all Lyapunov exponents are negative and therefore motion is regular, system is dissipative, Kaplane-Yorke dimension equal to zero and attractor is stationary point (corresponded to permanent rotation of main body). The system regular motion represents transition to permanent rotation about body z-axis (x(t)0, y(t)0, z(t)z*=const). Numerical research results (fig.9) demonstrate signs of regular motion. Comment about application of NL-gyrostat case. The Newton-Leipnik system describes attitude motion of spacecraft with linear feedback control [3]. NL-gyrostat results can be applied to simulation of perturbed attitude nonregular motion of gyrostat-spacecraft.

17


(a)

(b)

(c)

(d)

(e) (f) Fig.7 ­ Numerical research results for NL-gyrostat chaotic motion with rotor relative angular moment variability ( = 0.01, w = 10)
18


(a)

(b)

(c)

(d)

(e) (f) Fig.8 ­ Numerical research results for NL-gyrostat quasichaotic motion with rotor relative angular moment variability ( = 0.01, w = 1)
19


(a)

(b)

(c)

(d)

(e) (f) Fig.9 ­ Numerical research results for NL-gyrostat regular motion ( = 0, w = 1, v = 0)
20


5 Conclusion Links between mathematical models of gyrostats and dynamical systems with strange attractors (Lorenz, RÆssler, Newton-Leipnik and Sprott systems) were established. In order to examine of perturbed motion several numerical techniques was used: time-history of phase coordinate, kinetic energy, power spectrum of fast Fourier transformation, asymptotics of Lyapunov exponents and gyrostat longitudinal axis vector hodograph, and PoincarÈ sections. Mentioned numerical techniques showed chaotic and quasichaotic behavior of motion. Cases for perturbed gyrostat motion with variable periodical inertia moments and with periodical internal rotor relative angular moment were considered. References
[1] E. N. Lorenz, Deterministic nonperiodic flow, J. Atmos. Sci. 20 (1963) 130 ­141. [2] O. E. RÆssler, An equation for continuous chaos, Phys. Lett. A 57 (1976) 397 ­398. [3] R. B. Leipnik, T. A. Newton, Double strange attractors in rigid body motion with l inear feedback control, Phys. Lett A 86 (1981) 63-67. [4] J. C. Sprott, Some simple chaotic flows, Phys. Rev. B 50 (1994) R647 ­R650. [5] J. C. Sprott, Simplest dissipative chaotic flow, Phys. Lett. A 228 (1997) 271 -274. [6] A. V. Borisov, I. S. Mamaev, Dynamics of the Rigid Body. Izhevsk. RCD. 2001. P. 384. In Russian . [7] M. Inarrea, V. Lanchares, Chaotic pitch motion of an asymmetric non-rigid spacecraft with viscous drag in circular orbit, Int. J. Non-Linear Mech. 41 (2006) 86­100.

[8] X. Tong, B. Tabarrok, F.P.J. Rimrott, Chaotic motion of an asymmetric gyrostat in the gravitational field, Int. J. Non-Linear Mech. 30 (3) (1995) 191­203.
[9] J. Kuang, S. Tan, K. Arichandran, A.Y.T. Leung, Chaotic dyna mics of an asymmetrical gyrostat, Int. J. Non-Linear Mech. 36 (2001) 1213-1233 [10] A. V. Doroshin, Analysis of attitude motion evolutions of variable mass gyrostats and coaxial rigid bodies system, Int. J. Non-Linear Mech. 45 (2010) 193­205.

[11] A.P.M. Tsui, A.J. Jones, The control of higher dimensional chaos: comparative results for the chaotic satellite attitude control problem, Physica D 135 (2000) 41­62. [12] L. Zhou, Y. Chen, F. Chen, Stability and chaos of a damped satellite partially filled with liquid, Acta Astronautica 65 (2009) 1628-1638. [13] J. Baek, Adaptive fuzzy bilinear feedback control design for synchronization of TS fuzzy bilinear generalized Lorenz system with uncertain parameters, Physics Letters A 374 (2010) 1827­1834. [14] G. Benettin, L. Galgani, A. Giorgilli, J.-M. Strelcyn, Lyapunov characteristic exponents for smooth dynamical systems and for Hamiltonian systems: A method for computing all of them. P.I: Theory. P.II: Numerical application, Meccanica 15 (1980) 9-30.

21


PART II. ANALYSIS OF ATTITUDE MOTION EVOLUTIONS OF VARIABLE MASS 2 GYROSTATS AND COAXIAL RIGID BODIES SYSTEM
1. Introduction Research of attitude motion of a system of coaxial rigid bodies and gyrostats always was and still remains one of the important problems of theoretical and applied mechanics. The dynamics of the attitude motion of rotating rigid bodies and gyrostats is a classic mechanical research object. Basic aspects of such motion were studied by Euler, Lagrange, Kovalevskaya, Zhukovsky, Volterra, Wangerin, Wittenburg. The main results of the attitude motion research can be found in appropriate treatises [1-5]. However, the study of the dynamics of rotating bodies and gyrostats is still very important in modern science and engineering. Among the basic directions of modern research within the framework of the indicated problem it is possible to highlight the following points: deriving exact and approximated analytical and asymptotic solutions [1-5, 25, 28], research of a stability of motion conditions [6-16], the analysis of motion under the influence of external regular and stochastic disturbance, research of dynamic chaos [17-22], research of non-autonomous systems with variable parameters [23-30]. N. Ye. Zhukovsky studied the motion of a rigid body containing cavities filled with homogeneous capillary liquid. The research showed that the equations of motion in such case can be reduced to the equations of the attitude motion of a gyrostat. Also analytical solutions of some special modes of motion of a gyrostat were found. The ordinary differential equations of a gyrostat attitude motion with constant angular momentum were solved analytically by Volterra. Volterra solution has generalized a similar analytical solution for a rigid body in case of Euler. In the works of Wangerin and Wittenburg solution of Volt erra is reduced to the convenient parameterization expressed in elliptic integrals. The analytical solution for attitude motion of heavy dynamically symmetric gyrostat, colligating a classic solution for a heavy solid body in the case of Lagrange, is given in a paper [25]. In the indicated published works solutions for all Euler angles (precession, nutation, and intrinsic rotation) are found in elliptic functions and integrals. Also modes of motion with constant and variable relative angular momentum of gyrostat rotor are considered. The issues of the rotational motion dynamics of a gyrostat are very important for numerous applications such as the dynamics of satellite-gyrostats, spacecrafts and aircrafts. The attitude dynamics of gyrostat satellites and dual-spin (double-rotation) satellites has been studied by a number of scientists [6-22]. Most of these efforts were aimed on finding the equilibrium states in the presence of external disturbance torques [6-9], on analysis of the stability of spinning satellites under energy dissipation [10-16]. Some authors recently have investigated bifurcation and chaos in the gyrostat satellites [17-22].

- Material corresponds to preprint version of article: A.V. Doroshin, Anal ysis of attitude motion evolutions of variable mass gyrostats and coaxial rigid bodies system, International Journal of Non-Linear Mechanics, Volume 45, Issue 2 (2010) 193­205
2

22


Despite above-mentioned wide spectrum of research results the stated problem still remains actual, especially for the variable structure rigid bodies systems and variable mass dual-spin spacecrafts. The possibility of change in some mass and inertia parameters and the structure variability can be explained by the fact that a spacecraft (SC) performs active maneuvers with the use of the jet engine. Any SC in orbit is affected by external disturbances of different kind, e.g. the solar radiation pressure, the gravity gradient torque, the magnetic torque caused by the Earth's magnetic field, or the aerodynamic torque due to the action of a resisting medium like the Earth's atmosphere. However all these external disturbances are not large in comparison with the jet engine thrust of the SC on the active motion stage (e.g. inter-orbital transfer, orbit correction, attitude reorientation). Moreover, variability of mass parameters (mass and moments of inertia) has a considerable influence on attitude dynamics. The change of the moments of inertia entails change of angular momentum, which is the basic characteristic of attitude motion. Thereupon mass (structure) variation is one of the primary factors determining attitude motion of a SC. For the purposes of better understanding the essence of this problem it is important to give a brief overview of the main considered engineering peculiarities of SC's active motion. An SC in order to perform an active maneuver (e.g. inter-orbital transfer) should create a jet engine thrust and thus obtain acceleration or braking momentum V (reorbit/ deorbit burn). This momentum should be generated exactly in a pre-calculated direction. Engine thrust is usually focused along the SC's longitudinal axis, therefore it is necessary to stabilize the longitudinal axis in order to ensure the accurate momentum generation. Stabilization of the longitudinal axis can be carried out in a gyroscopic mode when SC spins around the longitudinal axis which is oriented in the calculated direction. Momentum generation is not instantaneous, it demands a continuous operation of the jet engine within several seconds (or minutes). During this period of time a SC performs two motions: trajectory motion of a center of mass and an angular motion around it. Such angular motion obviously changes the location of the longitudinal axis and, hence, a direction of thrust. The time history of thrust direction strongly affects the value and direction of a transfer momentum deviation. Consequently, the transfer is performed to the orbit different from the desired one. There is a "scattering" of thrust (Fig. 1). Therefore, it is very important to take SC angular motion into account during the analysis of the powered trajectory motion. It is necessary to obtain the angular motion which ensures that SC's longitudinal axis (and the thrust vector) performs precessional motion with monotonously decreasing nutation angle. Thus the longitudinal axis travels inside an initial cone of nutation and the thrust vector naturally comes nearer to an axis of a precession which is a desired direction of transitional momentum output ("is focused" along a necessary direction). When the angular motion does not provide a monotonous decrease in nutation angle the longitudinal axis moves in a complicated way. In such case the thrust vector also performs complicated motion and "scatters" the transitional momentum. A transfer orbit scatters as well. Among the works devoted to the rigid bodies systems with variable mass and inertia parameters it is possible to mark the following [17, 18, 23, 24, 28, 30]. The work [18] contains the analysis of chaotic behavior of a spacecraft with a double rotation and time-dependent moments of inertia during it's free motion. The main investigation results of variable mass system dynamics should be found in the monographies [23, 24]. These results include Ivan V. Meschersky theory of motion of bodies with variable mass, theory of "short-range interaction" and "solidification (freezing)". The equations of variable mass dynamically symmetrical coaxial bodies system were developed
23


in papers [28]. Also in [27] the attitude motion of coaxial bodies system and double rotation spacecraft with linear time-dependent moments of inertia were analyzed and conditions of motion with decreasing value of nutation were found. The results [27] can be used for the analysis of attitude mot ion of a dualspin spacecraft with an active solid-propellant rocket engine. Current material represents continuation of the research described in [27 -30] and is devoted to the dynamics of variable mass coaxial bodies systems, unbalanced gyrostats and dual-spin spacecrafts. The paragraph has the following structure: Section 1 ­ introduction of the primary theoretical and physical background, Section 2 ­ mathematical definition of the coaxial bodies attitude motion problem in terms of the angular momentum, Section 3 ­ main equations of attitude motion of two variable mass coaxial bodies system and unbalanced gyrostat, Section 4 ­ development of research method for the attitude motion of variable mass coaxial bodies and unbalanced gyrostat, Section 5 ­ examples of analysis of the attitude motion of variable mass unbalanced gyrostat, Section 6 ­ conclusion. 2. Problem Definition Below we will derive the equations of motion of a system of k coaxial rigid bodies of variable mass with respect to translating coordinate frame OXYZ. The motion of the system is analyzed with respect to the following coordinate systems (Fig. 2): P is a system of coordinates, fixed in absolute space, OXYZ is a moving coordinate system with origin at the point O, the axes of which remain collinear with the axes of the fixed system during the whole time of motion, and Oxi yi z i are systems of coordinates with a common origin, rigidly connected to the i-th body (i = 1, 2, ... , k), rotating with respect to the system OXYZ. OXYZ system has its origin in a point lying on the common axis of rotation of the bodies and matching with the initial position of the centre of mass t t0 : C O . Points in the different parts of the system are distinguished by the body they belong to, and in all expressions they are indicated by the subscript i (where i is a number of an appropriate body). To construct the obtain a relative velocity on it. In such case the th [23], written with respect equations of motion we use the "short-range" hypothesis ­ particles which when separated from the body no longer belong to the body and have no effect eorem on the change in the angular momentum of a system of variable mass to the fixed system of coordinates P, has the following form:
k dK P R M e M P Sie , P dt i 1 dm i Sie r i v i , dt i

(1)

where v i r i , m i is mass of appropriate point, M e is the principal moment of the external forces, P
M R is the principal moment of the reactive (jet) forces, and S ie is the sum of the angular momentum of P

the particles of body i, rejected in unit time in their translational motion with respect to the fixed system of coordinates. The angular momentum of a system of k bodies in coordinates system P (Fig. 2) is defined by the following expression:

24


KP


i 1 k
i

k

r i m i v

i


i 1 k
i

(r0 i ) m i ( v 0 i i )
i ,O

(2)


i 1

K

rO mi v Ci Ci mi v O ,



where

K

i ,O




i

m i i i i ,

mi


i

m i ,

vCi vO i Ci ,
is the radius vector of the centre of mass C i of body i in the OXYZ system and i is the absolute angular velocity of body i (and coordinates system Oxi yi z i ). In order to write the theorem of change in the angular momentum in the OXYZ coordinate system we need to implement some auxiliary expressions:
Ci



d

i

d dt dt w Ci w O

dt dv Ci

i i ,

d

Ci

dt

i Ci qCi ,

v O i Ci w Ci i qCi , i Ci i i Ci ,

(3)


i

dm dt

i



dm i i i i Ci mi qCi , dt



where q

Ci

is the relative velocity of the centre of mass C i due to a change in its position with respect to

the bodies, due to the variability of their masses; w Ci is an acceleration of the point of body i, that currently matches with its center of mass, i.e. it is an acceleration of translation for a center of mass Ci , wO is the acceleration of point O. Let's prove validity of the last expression from an expression group (3):

25



i

dm dt

i



i i i




i i

dm dt

i



i

d i dt


i

m i i
Ci



i

m i i i i

d i mi dt



m

i



dm i i Ci mi Ci i m i i i dt dm i i Ci mi q Ci i Ci i mi Ci dt dmi i C mi qCi . dt i





The relative motion of center mass with a velocity q

Ci

can be illustrated with Fig. 3, which

describes a burning process and corresponding geometric displacement of center of mass. Using expressions (2) and (3) it is possible to calculate the angular momentum derivative:
dK P dt


i 1

k

dK i ,O dmi dt rO dt v

Ci

rO mi w Ci i q Ci



Ci



qCi mi v

O

(4)

dmi v O Ci mi w O . dt

Let's transform the terms on the right hand side of the equation (1):
e M e rO F e M O , P R R M P rO R M O ,

Sie rO

dmi v O rO mi i qCi dt dm dm rO i i Ci mi qCi i Ci v dt dt dm i i i i , dt i









(5)
O





where F e


i 1

k

Fie is resultant of system of external forces, R


i 1

k



R i

is resultant of reactive (jet)

R e forces, M O , M O are the principal moments of the external and reactive forces with respect to the point O. Using expressions (4) and (5), after like terms cancellation we can rewrite a theorem (1) in the following form:

26




dK i ,O dt rO mi w Ci Ci mi w O i 1 e R e R rO F rO M O M O
k

(6)




i 1
i

k

i

dm dt

i



i i .



From the definition of center of mass and from the theorem on the motion of the center of mass of a system of variable mass [23] the following expressions must hold:
mi w Ci Fie iR


j i R

Nij , Nij N ji ,


i 1

k

mi w Ci F ,
e


i 1

k

(7)
mi Ci mC ,

where m m(t )


i 1

k

mi (t ) is mass of the system, C is vector of center of mass of the system, N ij are

internal forces of interaction between bodies. Using expressions (6) and (7) we can write the theorem on the change in the angular momentum with respect to OXYZ system [28]:


i 1

k

dK

i ,O

dt

e MO M

R O




i 1
i

k

i

dm dt

i



i

i



(8)

C mw O .

Expression (8) corresponds to the assertion of the well-known theorem [23], taking into account the grouping of the terms according to the membership of the points of the body i (i = 1, ..., k). Using the idea of a local derivative for the angular momentum vector of each body in the system of coordinates Oxi yi z i connected with the body, rotating with respect to OXYZ with angular velocity i Eq. (8) can be rewritten as follows:


i 1

k



dK i ,O dt

i K
Oxi yi zi

i ,O

i i C mw O .

(9)

e R MO MO


i 1
i

k

i

dm dt

i





The subscript outside the brackets of the local derivatives indicates a coordinate system in which they were taken. Equation (9) expresses a vector-based form of the theorem of the change in the angular momentum of bodies of variable mass with respect to the translating axes. 3. Attitude Motion of Two Variable Mass Coaxial Bodies System
27


We will consider the motion of a system of two bodies, where only 1st has a variable mass. Body 2 does not change its inertial and mass characteristics, calculated in the system of coordinates Ox2 y2 z 2 connected to the body, and, consequently, produces no reactive forces. This mechanical model can be used for research of attitude motion of dual-spin SC with operating solid-propellant rocket engine (coaxial body 1). We will write the angular velocities and the angular momentum of the bodies in projections onto the axes of their connected systems of coordinates:
i pi ii qi ji rik i , K
i ,O

Ii i ,

(10)

where Ii are inertia tensors of body i; ii , ji , k i are the unit vectors of the system Oxi yi zi . If both tensors are general then angular momentum of the bodies in projections onto the axes of their connected systems of coordinates is defined by

K K

1,O 2 ,O

A1 (t ) p1i1 B1 (t )q1 j1 C1 (t )r1k1 , A2 p2i 2 B2 q2 j2 C2 r2k 2 ,

where Ai, Bi, Ci are general moments of inertia of body i, calculated in the corresponding system of coordinates connected to the body. The bodies of the system can only rotate with respect to one another in the direction of the common longitudinal axis, which coincides with Oz2 (and with Oz1 ). Here we will denote the angle and velocity of twisting of body 1 with respect to body 2 in the direction of the longitudinal axis Oz2 by



Ox1 , Ox2

and

respectively. The angles {, , } of spatial orientation of the coaxial

bodies with respect to the translating system of coordinates OXYZ are indicated in Fig. 4. The ratio between the angular velocities and the angular accelerations of two bodies in vector form are defined by

1 2 ,
where k
1

1 2 ,

(11)

is the vector of the relative angular velocity of the bodies, which has the only

projection ­ onto common axis of rotation Oz2 . The ratio between the components of the angular velocities for the two bodies is expressed by the following equations:
p1 p2 cos q2 sin , q1 q2 cos p2 sin , r1 r2 .

(12)

The theorem on the change in the angular momentum (9) in translating system of coordinates OXYZ can be rewritten in the form:

28




dK1, dt

O




Ox1 y1 z1


1

1

dm dt

1



1 1



dK 2, dt

O




Ox2 y2 z
2


i 1

2

i K

i ,O

(13)

e R M O M O C mw O ,

where

1 x1 i1 y1 j1 z1 k1.
By projecting the expression inside the square brackets (Eq. (13)) onto the axes of the system Ox1 y1 z1 and using expressions (10) we obtain:
L dK1,O dt
Ox1 y1 z1


1

1

dm dt

1



1

1



I1, xx (t ) p1 I1, xy (t )q1 I1, xz (t )r1 i1 I1, yx (t ) p1 I1, yy (t )q1 I1, yz (t )r1 j1 I1, zx (t ) p1 I1, zy (t )q1 I1, zz (t )r1 k1.



(14)

During the simplification of equation (14) terms containing the derivatives of time-varying moment of inertia cancel out with terms following from the sum in square brackets (vector L). This is vividly reflected in the projection of L onto the connected axis Ox1 :
d K d K
1,O x 1

Lx1



dt
1,O x 1




1

1 2 1 dm dt
1

1



1 1 m x1



1



dt q1

p1
1


1


1

y21 z21 dm dt
1


1


1

dm dt

x1 y1 r1



x1 z

d I1, xx (t ) p1 I1, xy (t )q1 I1, xz (t )r1 dt dI1, xy (t ) dI (t ) dI (t ) p1 1, xx q1 r1 1, xz dt dt dt I1, xx (t ) p1 I1, xy (t )q1 I1, xz (t )r1.

If tensors of inertia remain general for each moment of time I1,ij (t ) 0, i j , then vector L may be rewritten:
L A1 t p1i1 B1 (t )q1 j1 C1 t r1 k1.

(15)

Taking expressions (14) and (12) into account, we can write Eq. (13) in terms of projections onto the
29


axis of Ox2 y2 z 2 system, connected with body 2. When changing from system Ox1 y1 z1 to system Ox2 y2 z we will use an orthogonal matrix of rotation by angle . As a result we obtain:

2



L 1 K

1,O Ox y z 111



dK 2,O 2 K dt
e O R O

2 ,O



(16)
Ox2 y2 z
2

M M C mw O ,
where



cos sin 0

sin cos 0

. 1 0 0

From the theorem on the motion of the centre of mass of a system of variable mass [23] the following expressions must hold:
m1 (t )w C1 F1e 1R N12 , m2 w
e C2

F2e N 21 ,
C2

m1 (t )w C1 m2 w F F F ,
e 1 e 2

F e 1R ,

(17)

N12 N 21 ,

where F e is resultant of system of external forces, 1R is reactive (jet) force, N ij are internal forces of interaction between bodies (j, k=1, 2). The motion of center of mass of body 1 is easier to analyze as compound motion, where the motion of body 2 is translational. Considering the last remark, expressions for the acceleration of the center of mass will have the following form

w

C2

w O 2 C2 2 2 C2 ,
(18)

w C1 w e w r w c , w e w O 2 C1 2 2 C1 , w r C1 C1 , w c 2 2



C1



,

where we ­ acceleration of translation, wr ­ relative acceleration, wc ­ Coriolis acceleration. Expressions (18) imply:

30


m1 (t )w C1 m2 w

C2

mw O 2 m

C

2 2 mC m1 (t ) C1 2 2





C1



C1

.

From last relation and (17) expression for acceleration wO follow:
wO 1e F 1R 2 mC 2 2 m m
C

m1 (t ) C1 C1 2 2





C1



(19)
.

The C mw

O



vector is represented using Eq. (19):
C mw O C F e 1R 2 m(t )C 2 2 m(t ) m1 (t ) C1 2
2



C

(20)

C1



C1



.

If the change of mass of body 1, which has general tensors of inertia, is uniform along the whole volume, then tensors of inertia remain general tensors of inertia and the centre of mass of body 1 remains on a common axis of rotation Oz2 . Thus we will consider, that the body 2 also has general tensors of inertia. The following expressions will take place in this case in terms of projections onto the axes of Ox2 y2 z 2 :
T 0, 0, , C1 0, 0, C1 (t ) ,
T

qC 0, 0, qC



T

,
T

C 0, 0, C (t ) ,



C mw

O Ox2 y2 z



2

C F e C 1R

(21)
Ox2 y2 z
2

2 m(t ) C p2 r2 q2 i 2 q2 r2 p2 j2 0 k 2 .

Let's transform the moments of external and reactive (jet) forces in expression (16):
e e M O C F e M C , R R R MO C 1 M C .

Taking the expressions (21) into account, we will write Eqs. (16) in the matrix form:

31




A t p C1 t B1 (t ) q1r1 1 1 B1 (t )q1 A1 t C1 (t ) p1r1 C t r B1 t A1 (t ) q1 p1 1 1 A2 p2 C2 B2 q2 r2 B2 q2 A2 C2 p2 r2 C2 r2 B2 A2 q2 p2 p2 r2 e R 2 M C M C m(t ) C q2 r2 0 q2 p2 .





(22)

Components of p1 , q1 , r1 in Eq. (22) must be expressed via p2 , q2 , r2 using (12). We will add an equation describing the relative motion of the bodies to the Eq. (16). A theorem on the change in the angular momentum projected onto the axis of rotation for the first body will have the following form:



L 1 K

1,O Oz 1



M

e 1,Oz

M zR M



C1 m1w O

Oz1

,

(23)

where M



is the moment of the internal interaction of the bodies (e.g. action of internal engine or

bearing friction), M ie,Oz is the moment of external forces acting only on body i. If tensors of inertia of body 1 remain the general ones for every moment of time and the centre of mass of body 1 remains on a common axis of rotation Oz2 , then Eq. (23) can be rewritten in the follow form:
C1 t r1 B1 t A1 (t ) q1 p1 M
e 1,Oz

M zR M .

(24)

We will supplement the dynamic equations (16) and (23) (or their simplified analogs (22) and (24)) by the following kinematic equations (Fig.4):

;


p2 sin q2 cos ,
(25)

1 p2 cos q2 sin , cos sin r2 p2 cos q2 sin . cos

Let's analyze the motion of a system of two dynamically symmetrical bodies, equations (22) and (24) will be written in the following form:

32


A(t ) p2 C (t ) A(t ) q2 r2 C1 (t ) q2 M
e C,x

M xR ,
R My ,

A(t ) q2 C (t ) A(t ) p2 r2 C1 (t ) p2 M C (t ) r2 C1 (t ) M
e C,z e C,y

(26)

M zR ,
e 1,Oz

C1 (t ) r M M zR M
2 where A(t ) A1 (t ) A2 m(t ) C (t ), C (t ) C1 (t ) C2 .

,

Systems (26) and (25) together form a complete dynamic system for the research of attitude motion of dynamically symmetrical unbalanced gyrostat with variable mass. 4. Research method of attitude motion of variable mass coaxial bodies and unbalanced gyrostat Let's refer to a motion of coaxial bodies (unbalanced gyrostat) of variable mass under an action of dissipative and boosting external moments depending on components of angular velocities. Let the gyrostat consists of dynamically symmetrical main body (coaxial body 2) of a constant mass and a rotor (coaxial body 1) of the variable mass, which remains dynamically symmetrical during modific ation of a mass (Fig. 4). The fixed point O coincides with an initial geometrical position of a system's center of mass. The unbalanced gyrostat has a varying relative angular velocity of rotor rotation around the main body. It is possible in connection with the existence of internal moment M acting between coaxial bodies. Let's assume there is a moment of jet forces only around a longitudinal axis Oz
1



R M xR M y 0 .

Let's implement the new variables corresponding to the magnitude G of a vector of transversal angular velocity and angle F between this vector and axis Oy2:
p2 G (t ) sin F (t ), q2 G (t ) cos F (t ).

(27)

Equations (26) will be rewritten in new variables as follows:
1 C (t ) A(t ) r2 C1 (t ) f F A(t ) M e M fG G, F G , r2 2,Oz , A(t ) C2 R e e C (t ) M M z M 1,Oz M 2,Oz . C1 (t )C2 C1 (t ) C2

F



G, F ,

(28)

In equations (28) the following disturbing functions describing exposures take place:

33


f f

G



G, F M G, F

e C,x

sin F M
e C,x

e C,y

cos F ,
e C,y

F

1 M G

cos F M

sin F .

We will consider a case when the module of a transversal angular velocity of main body is small in comparison to relative longitudinal rotation rate of the rotor:
2 2 p2 q2 1.

(29)

From spherical geometry the formula for a nutation angel (an angle between axes OZ and Oz2) follow
cos cos cos .

We will assume angles and to be small O , O . Then the nutation angle will be defined by the following approximated formula:

2 2 2 .

(30)

Using the expressions (22) and kinematic equations (20) we can write (second order infinitesimal terms are omitted):
G cos (t ), G sin (t ), r2 , , (t ) F (t ) (t ).

(31)

Function (t ) is a phase of spatial oscillations. Precession motion of the gyrostat with small nutation angles is obviously described by a phase space of variables , . The phase trajectory in this space completely characterizes motion of the longitudinal axis Oz2 (an apex of the longitudinal axis). Therefore our further researches will be connected to the analysis of this phase space and chances of behaviors of phase curves in this space. We can develop a special qualitative method of the analysis of a phase space. Main idea of the method is the evaluation of a phase trajectory curvature in the phase plane , . On the indicated plane the phase point will have following components of a velocity and acceleration:

V , V , W , W .
With the help of expressions (26) the curvature of a phase trajectory (k) is evaluated as follows:

k 2



2



2

23



2 G2 .

(32)

If curvature magnitude increase, there will be a motion on a twisted spiral trajectory similar to a steady focal point (Fig. 5, case "a") and if decreases - on untwisted. On twisted spiral trajectory motion condition can be noted as:

k



kk 0 G G2 0.
34

(33)


For the analysis of the condition realization it is necessary to study a disposition of zero points (roots) of a following function:

P(t ) G G2 .
Function (34) will be defined as a Different qualitative cases zero points of (Fig.5). In the first considered slice of time t 0, T ,

(34)

function of phase trajectory evolutions. of phase trajectory behaviors are possible depending on function P(t) case (Fig. 5, case "a") the function is positive and has no zero on a thus the phase trajectory is spirally twisting. In the second case (Fig.

5, case "b") there exists one zero point and there is one modification in a monotonicity of the trajectory curvature. The Cornu spiral, also known as clothoid, take place in case "b". The third case (Fig. 5, case "c") represents a number of zero points and the trajectory has alternation of untwisted and twi sted segments of motion; also there are some points of self-intersection. 5. Research of attitude motion. 5.1. Example 1. As an first example we will refer to a motion of coaxial bodies of variable mass under the influence of constant internal moment ( M const ) and constant moment of jet forces ( M zR const ). The analysis of a phase space is conducted using a developed method of curvature evaluation. We will suppose that the mass and moments of inertia are linear functions of time:

m1 (t ) mr kt , A1 (t ) m1 (t ), C1 (t ) m1 (t ), A2 Am const , C2 Cm const ,
(35)

where mr is initial mass of rotor, k is rate of mass change, and , are constants. Dependencies (35) are valid for a dual-spin SC when one of coaxial bodies is a solid-propellant rocket engine with packed and roll shaped grains. A linear law of mass change provides a constant thrust. In the rocket engines of a described type the grain usually burns uniformly over the whole volume, the grain density changes uniformly as well (Fig.3-c). The center of mass of an engine-body will show no displacement relative to the bod y, therefore C lr const. Center of mass of the main
1

body does not move as well, because the body doesn't change its mass, therefore C lm const. Let's
2

mark that constants lr and lm are in fact the distances between the bodies' center of masses and the point O (Fig.4): ln OC2 , lr OC1 . For example, if body1 is solid cylinder, then the constants and are

H 2 12 R2 4 lr2 , R2 2 ,
where H is the height of a cylinder and R is the radius. Magnitude (t) will be defined by a linear-fractional form:

C (t )

lm m2 lr m1 (t ) . m2 m1 (t )

(36)

At t=0 the system's center of mass C and the point are matching, therefore
35


C (0) 0, lr m1 (0) lm m2 ,
lm 0, lr 0.

On base (36) and (35) we will write time-dependences for A(t) and C(t):
A(t ) A at C (t ) C ct , k 2lr2t 2 , m kt

(37)

where
A Ar Am , Ar A1 (0), C Cr Cm , Cr C1 (0), a m1 (0), c m1 (0).

In a considered case equations (28) will obtain the following form:

G 0, 1 F C (t ) A(t ) r (Cr ct ) , A(t ) r2 M , Cm C (t ) M M zR . (Cr ct )Cm Cr ct
(38)

Analytical solutions for angular velocity r2(t) and t are derived from equations (38):
r2 r0 M t, Cm

0 s1t s2 ln 1 c1t ,
R z

1 s1 M Cm , s2 M M c



(39)
, c1 c / Cr .

Kinematic equations (31) can be used to receive a solution for the angle

:

0 r0t

M 2 t. 2Cm

Expression for a time derivative of a spatial oscillations phase can be obtained using (31), (39) and (38):

1 C (t ) A(t ) s1t r0 A(t )


(40)

(Cr ct ) 0 s1t s2 ln 1 c1t s1t r0 .

Formula (40) make it possible to receive an explicit expansion for evolution function (34), but
36


this expression is difficult to analyze. It is reasonable to expand this expression into a series:

P(t )


i 0



fi t i .

(41)

Writing formula (41) on the basis of (34) we don't take into account a constant multiplier G0. Further we will investigate the simplest case when the expansion for P(t) has only linear part (other terms of expansion are not taken into account). On the basis of e xpressions (40) we can get a polynomial of the first degree for the phase trajectory evolutions function (34):

P(t ) f 0 f1t ,
where
f0 Cr 2 aCr cA 0 AM zR 0 , 3 A

(42)

f1

2 3a 2Cr2 0 Cr 0 3 4aM zR A4 A C k 2l 2 2 0 r r 2ca m 1 2 2 c 2 0 c 0 M 3M A

R z



M

R2 z



.

There is a unique zero point of function P(t): t1 f0 / f1. For implementation of a condition (33) of twisted spiral motion it is necessary for the polynomial to be steady ( t1 0 ) and positive for all t 0 . It is possible only in case the following conditions fulfills:

f0 0, f1 0.
We will consider a case when following contingencies are correct:

(43)

r0 0, 0 0, M zR 0.
In this case value f0 will be positive if following condition is true:

c Cr a A ,

0 cA aCr AM zR .

(44)

In order f1 to be also positive the following conditions must be satisfied:

0 c
3M zR



Cr k 2lr2 2am

R Mz , c 0 M .



M

R2 z 0



(45)

c

37


Also f1 > 0 if following conditions hold true:

Cr k 2lr2 2acm,

M 3M zR 0.

(46)

Figure 6 illustrates the results of evolution function and appropriate phase trajectories numeric calculations. Figures Fig.6-a and Fig.6-b demonstrate the situation when (44) and (45) are satisfied, fig.6-c show the opposite case. Point indicates the beginning of phase trajectories. In case "a" evolution function has two roots and phase trajectory has three evolutions: twisting-untwisting-twisting. In case "b" evolution function has no root and single evolution of phase trajectory takes place ­ this evolution is spiral twisting. System parameters and initial conditions for obtained solutions are listed in table 1. Table 1 Value of the system parameters for figure 6 Quantity M, Nm MRz, Nm 0, radian/sec G0, radian/sec Am, kgm2 Ar, kgm2 Cm, kgm2 Cr, kgm2 , kgm2/sec c, kgm2/sec lr, m k, kg/sec m1(0), kg m2, kg a 1 15 10 0.2 2.5 2.5 1 1.5 0.08 0.08 0.5 1 35 35 b -10 10 1 0.2 2.5 1.5 1 1.5 0.05 0.08 0.4 1 25 35 c 200 0.35 16 0.2 2.5 2 1 2 0.1 0.08 0.6 1.2 45 35

Conditions (44) and (45) can be used for the synthesis of dual-spin spacecraft parameters. In order to enhance the accuracy of SC's longitudinal axis positioning it is necessary to realize preces sion motion with a decreasing nutation angle. This motion is realized when the conditions (44) and (45) are satisfied. For realization of more accurate researches, certainly, it is necessary to take higher degrees polynomials P(t) (34) into account. However it was shown that an implemented analysis provides an adequate description of the precessional motion evolutions of variable mass coaxial bodies. Examined above case of investigation does not take into account many important aspects of variable mass coaxial bodies motion. However the introduced example has illustrated the approach to a research of non-autonomous dynamical systems of indicated type. 5.2. Example 2. Let's refer to the other mode of motion with the following external and internal dissipative force moments:

38


M M

e C,x e 2 ,Oz

p, r ,

M M

e C, y

q , r .

e 1,Oz

(47)

Constants , , describe the influence of resisting substance on a gyrostat and the dissipation of energy. Let the mass and inertial parameters be described by polynomial functions of time:
2 A(t ) A2 A1 (t ) m(t ) C (t )


i 0

n

ait i ,
(48)

C1 t


i 0

m

cit ,

i

Dynamical equations of motion will have the following form (26):
A t p2 C (t ) A(t ) q2 r2 C1 (t )q2 p2 , A t q2 C (t ) A(t ) p2 r2 C1 (t ) p2 q2 , C (t )r2 C1 (t ) M zR r2 r2 , C1 (t ) r2 M M zR r2 ,

(49)

where C (t ) C2 C1 (t ). The last equation in (49) provides an equation and a general solution for an absolute longitudinal angular velocity of rotor :
1 M M zR , C1 (t ) 1

t






M M

R z



(50)

1





M M zR 0 exp J C (t ) ,

where
r2 , J C (t )

C t
0 1

t

dt

.

(51)

Using expressions (50) we can notice that the third equation in (49) gives an equation and a general solution for a longitudinal angular velocity of the main body:
C2 r2 M r2 , M r2 t r0 exp t M . C2

(52)

In the case involved disturbing functions will have the form:

39


f

G



G, F G, f

F



G, F 0,

(53)

thereby first two equations can be rewritten the following way:

G G , A(t ) 1 F C2 r2 C1 (t ) A(t )r2 . A(t )
This implies a solution for the amplitude of the transverse angular velocity:
G t G0 exp J A (t ) ,

(54)

(55)

where
J A (t )

A t
0

t

1

,

G0 0.

(56)

The value of integrals contained in expressions (50) and (55) can be calculated analytically:
J

A

t
i 1 t

n

ln t i A
i

t

,
0

J

C

t C1 t i 1 0

dt

m

ln t i , C1 i 0

t

(57)

where i , i are the roots of A(t ), C1 (t ) polynomials, which contain no real roots within a considered interval t 0, T , because the moments of inertia are strictly positive quantities. Let's mark that formulas can be checked by differentiation. Expressions (50), (52) and (55) provide final analytic solutions for the longitudinal angular velocities and the amplitude of transverse angular velocity of the gyrostat. Using the first equation in (54) we can obtain an expression for the evolution function:
1 d 2 2 P(t ) G . A(t ) 2 dt

(58)

Let's consider a case when there are no moments of internal interaction and jet forces and the main body has no initial longitudinal angular velocity [26 ­ 29]:
M M zR 0, r0 0.

(59)

Let's perform some supplementary transformations of quantities presented in expression (58). Replacing derivatives , with corresponding right hand sides of equations (54), (50) and considering that r (t ) 0 , we can write the following expressions:

40


2 2 C1 (t ) 2 , A2 (t ) d 2 2 2 3 C1 C1 A AC12 C1 A . dt A

(60)





Now function (58) will have the following form:
P(t ) G2C1 C1 A AC1 A C1 . A3





(61)

Solutions (55) and (50) imply that G 0, 0 , therefore a multiplier, placed before brakets in (61), is strictly positive over the whole interval t 0, T . Considering the last remark, function P(t) may be replaced by the following expression:
P(t ) C1 A AC1 A C1 .

(62)

Expressions (62) and (48) imply, that function P(t) ­ is a polynomial with a degree of N=m+n-1, that is why a theoretically valid number of phase trajectory evolutions can not exceed N+1. Using function (62) we can obtain the following limitations for the moment of inertia functions, providing the twisted-in phase trajectory and therefore the decreasing amplitude of nutation angle:
C (t ) A(t ) C (t ) A(t ) , C (t ) A(t ) .

(63)

This limitations particularly result in the condition for the linear functions of moments of inertia in the absence of dissipative moments ( = =0), given in articles [27, 28]:

a A c Cr .
In found. It trajectory have been general case, when constraints (59) are not satisfied, conditions similar (63) have not been is possible in this case to receive numerical results for the evolution function and phase (Fig. 7). The following parameters and polynomial time-dependences of inertia moments used for computation:
A(t ) 0.0001t 4 0.0064t C (t ) 0.0001t 4 0.0027t
3 3

0.1053t 2 0.4491t 3, 0.0344t 2 0.1997t 1.92, C2 1.5 (kgâ m 2 ); T 24 (sec);

0.02, 0.01, 0.03 (kgâ m 2 / sec);
M 11, M zR 21 (kgâ m 2 / sec 2 ) ; G0 0.1, r0 15, 0 10 (rad / sec).

It can be noticed (Fig. 7) that function P(t) has five roots on examined time-interval and consequently phase trajectory has six evolutions.
41


5.3. Example 3. Omitting the solution details we consider the numerical simulation of previous example when inertia moments are simplest harmonic function:

A(t ) a1 sin t a0 ,

a0 a1 0,

C1 (t ) c1 cos t c0 , c0 c1 0.

(64)

For this case solutions (50) and (55) remain valid, but expressions for integrals (51) and (56) take on a value:
J C (t ) arctan 2 c0 c12 2 c c tan t 2 0 1 2 c0 c12

t div , 2 2 J A (t ) arctan 2 a0 a12 2 a tan t 2 a 1 0 2 2 a0 a1

(65)

t 2 2 where operation (x div y) corresponds to evaluation Evolution function take on form:

div , of integer part of division x/y.

K z2 (t ) K z (t ) C (t ) P(t ) G (t ) 3 A 2 A (t ) A (t )







(66)

M zR r2 (t )

,

where
K z (t ) C2 r2 C (t ) .

Though quite simple analytical description of evolution function (66), the phase trajectory behavior in considered case can be complex. Phase trajectory can be regular (Fig.8-a), or can demonstrate unpredictable forms, which typical in chaotic d ynamics (Fig.8-b, c). On Fig.8 three cases of phase trajectories are shown. These trajectories are cal culated for parameters from table 2; constants for inertia moments dependences (64) are equal a0=c0=2, a1=c1=1 (kgâm2), = 3 (1/sec). Case "a" correspond to motion without reactive and internal forces moments M M zR 0 . In this case quasi-periodic evolution function (with slow damped amplitude) take place and phase trajectory is also quasi-periodic and regular. In presence reactive and internal forces moments (cases "b" and "c") evolution function become nonperiodic with complex changing amplitude rate. In these cases phase trajectories become nonregular

42


and similar to chaotic. Cases "b" and "c" correspond to a positive , , 0 and negative , , 0 dissipation. It is necessary to note, what all calculations were conduct in MAPLE 11 [31] with use of numerical solution of stiff initial value problem (absolute error tolerance is equal 0.0001). Table 2 Value of the system parameters for figure 8 Quantity M, Nm MRz, Nm r0, radian/sec 0, radian/sec G0, radian/sec C2, kgm2 , kgâm2/sec , kgâm2/sec , kgâm2/sec T, sec 6. Conclusion The material described a research of the phase space of non-autonomous dynamical system of coaxial bodies of variable mass using a new method of phase trajectory curvature analysis. Developed method allows to estimate the phase trajectory form. System motion can be both simple regular and very complicated nonregular (chaotic). Regular motions are realized at evolution functions with finite number of roots (polynomial) or at periodic evolution functions. Complicated non-regular motions arise at nonperiodic alternating evolution functions with infinite number of roots. Results of the research have an important applied value for the problems of space flight mechanics and especially for coaxial spacecrafts. References [1] J. Wittenburg, Dynamics of Systems of Rigid Bodies, B.G.Teubner, Stuttgart, Germany, 1977. [2] S.V. Kovalevskaya, Sur le probleme de la rotation d'un corps solide autor d'un point fixe. Acta. Math. 12 (1889). [3] V. Volterra, Sur la theories des variations des latitudes, Acta Math. 22 (1899). [4] N. Ye. Zhukovsky, The motion of a rigid body having cavities filled with a homogeneous capillary liquid. In Collected Papers, Vol. 2. Gostekhizdat, Moscow (in russian), 1949, pp. 152-309. [5] A. Wangerin, ýber die bewegung miteinander verbundener kÆrper, UniversitÄts-Schrift Halle, 1889. [6] T.R. Kane, D.L., Mingori, Effect of a rotor on the attitude stability of a satellite in circular orbit, AIAA J. 3 (1965) 938.
43

a 0 0 15 -18 0.01 2.5 0.00001 0.00001 0.00001 50

b 0.1 0.1 -1 -10 0.01 2.5 0.001 0.001 0.001 470

c 0.5 0.03 1 -10 0.01 2.5 -0.001 -0.002 -0.003 500


[7] P.W. Likins, Spacecraft Attitude Dynamics and Control - A Personal Perspective on Early Developments, J. Guidance Control Dyn. Vol. 9, No. 2 (1986) 129-134. [8] P.W. Likins, Attitude Stability Criteria for Dual Spin Spacecraft, Journal of Spacecraft and Rockets, Vol. 4, No. 12 (1967) 1638-1643. [9] G.J. Cloutier, Stable Rotation States of Dual-Spin Spacecraft, Journal of Spacecraft and Rockets, Vol. 5, No. 4 (1968) 490-492. [10] D.L. Mingori, Effects of Energy Dissipation on the Attitude Stability of Dual-Spin Satellites, AIAA Journal, Vol. 7, No. 1 (1969) 20-27. [11] P.M. Bainum, P.G. Fuechsel, D.L. Mackison, Motion and Stability of a Dual-Spin Satellite With Nutation Damping, Journal of Spacecraft and Rockets, Vol. 7, No. 6 (1970) 690-696. [12] K.J. Kinsey, D.L. Mingori, R.H. Rand, Non-linear control of dual-spin spacecraft during despin through precession phase lock, J. Guidance Control Dyn. 19 (1996) 60-67. [13] C.D. Hall, Escape from gyrostat trap states, J. Guidance Control Dyn. 21 (1998) 421-426. [14] C.D. Hall, Momentum Transfer Dynamics of a Gyrostat with a Discrete Damper, J. Guidance Control Dyn., Vol. 20, No. 6 (1997) 1072-1075. [15] A.E. Chinnery, C.D. Hall, Motion of a Rigid Body with an Attached Spring-Mass Damper, J. Guidance Control Dyn. Vol. 18, No. 6 (1995) 1404-1409. [16] A.I. Neishtadt, M.L. Pivovarov, Separatrix crossing in the dynamics of a dual-spin satellite, J. Appl. Math. Mech., V. 64 (5) (2000) 709-714. [17] M. Inarrea, V. Lanchares, Chaotic pitch motion of an asymmetric non -rigid spacecraft with viscous drag in circular orbit, Int. J. Non-Linear Mech. 41 (2006) 86­100 [18] V. Lanchares, M. Iarrea, J.P. Salas, Spin rotor stabilization of a dual -spin spacecraft with time dependent moments of inertia, Int. J. Bifurcation Chaos 8 (1998) 609­617. [19] T.S. Parker, L.O. Chua, Practical Numerical Algorithms for Chaotic System s, Springer, Berlin, 1989. [20] A. Guran, Chaotic motion of a Kelvin type g yrostat in a circular orbit, Acta Mech. 98 (1993) 51­ 61. [21] X. Tong, B. Tabarrok, F.P.J. Rimrott, Chaotic motion of an asymmetric gyrostat in the gravitational field, Int. J. Non-Linear Mech. 30 (3) (1995) 191­203. [22] J. Kuang, S. Tan, K. Arichandran, A.Y.T. Leung, Chaotic dynamics of an asymmetrical gyrostat, Int. J. Non-Linear Mech. 36 (2001) 1213-1233 [23] A.A. Kosmodem'yanskii, Course in Theoretical Mechanics, Part 2, Prosveshcheniye, Moscow, 1966 (in russian). [24] F.R. Gantmaher, L.M. Levin, The theory of rocket uncontrolled flight, Fizmatlit, Moscow, 1959 (in russian).
44


[25] V.S. Aslanov, A.V. Doroshin, About two cases of motion of unbalanced gyrostat. Izv. Ross. Akad. Nauk. MTT. V.4 (2006). (Proceedings of the Russian Academy of Sciences. Solid bodies mechanics. In russian) [26] V.S. Aslanov, A.V. Doroshin, Stabilization of the descent apparatus by partial twisting when carrying out uncontrolled descent. Cosmic Research, Vol. 40, No. 2 (2002). 193-200. [27] V.S. Aslanov, A.V. Doroshin, G.E. Kruglov, The Motion of Coaxial Bodies of Varying Composition on the Active Leg of Descent, Cosmic Research, Vol.43, No. 3 (2005) 213-221. [28] V. S. Aslanov, A. V. Doroshin, The Motion of a System of Coaxial Bodies of Variable Mass, PMM. J. Appl. Math. Mech. Vol.68 (2004) 899-908. [29] A.V. Doroshin, Stabilization of a descent space vehicle with double rotation at presence of small dynamic asymmetry, The Technological collection on space-rocket engineering. "Calculation, designing and tests of space systems", Space-Rocket Corporation (Volga branch, Samara) "RKK Energia" (2001) 133-150 (in russian). [30] A.V. Doroshin, Phase Space Research of One Non-autonomous Dynamic System, Proceedings of the 3rd WSEAS/IASME International Conference on DYMAMICAL SYSTEM and CONTROL. Arcachon, France (2007) 161-165. [31] M.B. Monagan, K.O. Geddes, K.M. Heal, G. Labahn, S.M. Vorkoetter, J. McCarron, P. DeMarco, MAPLE. User Manual. Waterloo Maple Inc. 2007.

45


Fig.1. Scattering of thrust and transfer orbit

46


Fig.2. Coaxial bodies system and coordinates system

Fig. 3. Change of center of mass position with respect to the bod ies, due to the variability of their masses (cases of body burn): a ­ right to left burn, b ­ left to right burn, c ­ burn with uniform reduction of body density

47


Fig.4. Two coaxial bodies system, coordinates systems and Euler's type angels

Fig. 5. Cases of phase trajectory behaviors

48


Fig.6. Evolution function and cases of phase trajectory evolutions depends on conditions (44) and (45) fulfillment: In cases "a" and "b" conditions are satisfied and first evolution of phase trajectory is spiral twis ting; In case "c" condition are unsatisfied and first evolution is untwisting

49


Fig.7. Numerical simulation results for the evolution function and phase trajectory

50


Fig.8. Evolution functions and phase trajectories in case with harmon ic inertia moments 51


PART III. ATTITUDE CONTROL OF SPIDER-TYPE MULTIPLE-ROTOR 3 RIGID BODIES SYSTEMS

1. Introduction Most research into attitude motions of rigid bodies systems always has been and still remains one of the most important problems of theoretical and applied mechanics. Dynamics of the attitude motion of such systems is a classical mechanical research topic. Basic aspects of this motion were studied by Euler, Lagrange, Kovalevskaya, Zhukovsky, Volterra, Wangerin, Wittenburg. However, the study of the dynamics of rigid bodies is very important in modern science and engineering. Among the basic directions of modern research into the framework of the indicated problem it is possible to highlight the following points: mathematical modeling and analysis of multibody systems motion [1], multibody spacecraft (SC) attitude dynamics and control [2]-[18], multibody systems approach to vehicle dynamics and computer-based technique [19], simulation of multibody systems motion [21], multibody dynamics in computational mechanics [20]. If we speak about practical use of system of rigid bodies dynamics research results we have to note first of all SC with momentum wheels, reaction wheels and control moment gyroscopes (dual-spin satellites, gyrostats, space stations, space telescopes, etc.) [1-18]. Let us briefly describe current rotortype systems for attitude control of SC. Reaction wheel is a spinning wheel that can be moved to change the orientation of SC to which the wheel is attached. Reaction wheels are used in many satellites, including the Hubble Space Telescope, to allow precise pointing. By attaching an electric motor to a heavy wheel, and spinning the wheel one way quickly, a satellite rotates the other way slowly by conservation of angular momentum. This method can provide very precise orientation, and does not require any fuel . The problem is that if there is some small continuous torque on a satellite, e.g., radiation pressure, the wheels end up spinning all the time, faster and faster, to counteract it. The solution of the problem is to have some way of dumping momentum. The usual technique is utilizing a set of electromagnets that can be used to exert a weak torque against the Earth's magnetic field. Momentum wheels are a different type of actuator, mainly used for gyroscopic stabilization of SC: momentum wheels have high rotation speeds and mass, while reaction wheels work around a nominal zero rotation speed. A control moment gyroscope (CMG) is an attitude control device which consists of a spinning rotor and one or more motorized gimbals that tilt the rotor's angular m omentum. As the rotor tilts, the changing angular momentum causes a gyroscopic torque that rotates the spacecraft. Certainly, a CMG is a very powerful device which differs from reaction wheels. The latter applies torque simply by changing rotor spin speed, but the former tilts the rotor's spin axis without necessarily changing its spin
- Material corresponds to preprint version of article: A.V. Doroshin, Attitude Control of Spider-type Multiplerotor Rigid Bodies Systems. Proceedings of the World Congress on Engineering 2009, London, U.K. Vol II, pp.1544-1549.
3

52


speed. CMGs are also far more efficient. For a few hundred watts and about 100 kg of mass, large CMGs have produced thousands of newton meters of torque. CMG has got "gyroscopic power increase". A reaction wheel of similar capability would require megawatts of power. Thus, we have a wide choice of devices for SC attitude control. Nevertheless, it is still possible to evolve new equipment for changing SC orientation. This paper sets out to develop a new multiple-rotor system which also can be used for attitude control of a SC. The new multirotor system differs from typical attitude control devices. It has no disadvantages of a reaction wheel, and on the contrary, it possesses the advantages of a control moment gyroscope. Due to the large number of rays with rotors, we called new systems a "Rotor-type hedgehog" and a "Rotor-type spider". The paragraph has the following structure: Section 1 is an introduction of the problem background, Section 2 comprises mechanical and mathematical models of the new multiple-rotor system, Section 3 is devoted to development of a new reorientation method of the multiple-rotor system, Section 4 contains analytical solution of multiple-rotor system attitude motion, Section 5 include numerical simulation of the method implementation, Section 6 is conclusion. 2. Mechanical and Mathematical Models of Multiple-Rotor Systems We shall investigate an attitude motion about fixed point O of multirotor systems which are depicted in Fig.1. Fig.1 contains mechanical models of multirotor systems: cases (a) and (b) correspond to "Rotortype spider" systems; case (c) corresponds to "Rotor-type hedgehog". Firstly we consider rotor-type spider with six rotors which spin about general orthogonal axis of main (central) body (Fig.1-a). Let's assume symmetry of rotors disposition with respect to point O and equivalence of their mass-inertia parameters. Angular momentum of the system in projections onto the axes of frame Oxyz connected with main body is defined by
K Km K
r

(1) (2)

Km

Ap Bq Cr

p 1 2 q , K I r 4 4 J 2I 3 r 5 6

where K m is angular moment of a main rigid body with resting ("frozen") rotors; K r is relative angular moment of rotors; =[p, q, r]T is vector of absolute angular velocity of main body; i is relative
angular velocity of i-th rotor with respect to main body; A, B, C are general moments of inertia of main

body; I is longitudinal moments of inertia of single rotor; J is equatorial moments of inertia of single rotor calculated about point O. Motion equations of the multirotor system can be obtained with help of law of changing of angular momentum in frame Oxyz
dK K M dt
e

(3)

53


(a)

(b)

(c) Fig.1. Multirotor systems
54


where M e is principal moment of the external forces. Eq. (3) can be rewritten as
e Ap I 12 C B qr I q 56 r 34 M x 12 56 e 34 Bq I A C pr I r p M y 34 12 e 56 Cr I B A pq I p q M z

(4)

In the last equations following terms take place
ij i j ,
A A 4 J 2I

B B 4 J 2I , C C 4 J 2I

(5)

We need to add equations for rotors relative motion. These equations can also be written on the base of the law of the change in the angular momentum
i I p 1 M 1i M 1ex ; I p 2 M 2 M i e i I q 3 M 3 M 3 y ; I q 4 M 4 M i e i I r 5 M 5 M 5z ; I r 6 M 6 M e 2x e 4y e 6z

(6)

where M

i j

is a principal moment of the internal forces acting between main body and j-th rotor;

M ex , M ey , M ez are principal moments of external forces acting only at j-th rotor. j j j

Equation systems (4) and (6) together completely describe the attitude dynamics of the rotortype spider (Fig.1-a). Motion equations (4) and (6) corresponding to the spider with six rotors can be generalized for description of attitude dynamics of rotor-type spider with 6N rotors (Fig.1-b). As presented in Fig.1-b multirotor system has got N rotors on every ray - N rotor layers (levels). Similarly to previous case we can obtain the same equation system (4) for attitude motion of the multirotor system with N rotor layers (levels), but equations (4) will receive the following new terms
ij


l 1

N

il jl ,

A A 4


l 1 N

N

J l 2 NI

B B4


l 1

N

(7)

J l 2 NI , C C 4


l 1

J l 2 NI

where kl is the relative angular velocity of the kl-th rotor (Fig.1-b) with respect to main body; Jl is equatorial moments of inertia of the kl-th rotor (correspond to the l-th layer of rotors) calculated about point O. Equations of rotors' relative motion are given by
I p 1l M 1il M 1elx ; I p i e I q 3l M 3l M 3ly ; I q i e I r 5l M 5l M 5lz ; I r
2l 4l 6l



i M 2l M i M 4l M i M 6l M

e 2 lx e 4 ly e 6 lz

(8)

where l 1..N and therefore we have got N systems like (8) for each kl-th rotor.
55


Equation system (4) with terms (7) and N systems like (8) completely describe of attitude dynamics of the rotor-type spider with 6N rotors (Fig.1-b). Thus, we have dynamic equations of attitude motion. Let's define kinematic parameters and corresponding kinematic equations. We will use well-know [23] Euler parameters 0 , 1 , 2 , 3 describing a finite rotation of main body by an angle
e cos , cos , cos
T

about an arbitrary unit vector

in inertial fixed frame OXYZ which coincides with the initial position of Oxyz

(Fig.2).

Fig.2. Finite rotation The Euler parameters are defined by
1 cos sin 0 cos 2 ; 2 cos sin ; cos sin 3 2 2 2

(9)

Following system of kinematical equation takes place for Euler parameters
2

(10)

where


0 1 2 3

,



0 p q r

p 0 r q

r r q 0 p p 0

q

(11)

Set of dynamic and kinematic equations completely describe of attitude motion of multirotor systems. 3. A New Method of Attitude Reorientation of Multiple-Rotor System First of all, we will give a series of definitions. Def.1. Conjugate rotors are pair rotors located in the same layer on the opposite rays. For example, rotor 3N and rotor 4N (Fig.1-b) are conjugate rotors (also rotor 12 and rotor 22, etc.). Def.2. Conjugate spinup mean a process of spinning up conjugate rotors in opposite directions up to a desired value of relative angular velocity with help of internal forces moments from main body.
56


Velocities of conjugate rotors will equal in absolute value and opposite in sign. Def.3. Rotor capture is an immediate deceleration of rotor relative angular velocity with help of internal forces moment from the main body. So, rotor capture means an "instantaneous freezing" of rotor with respect to the main body. The capture can be performed with help of gear meshing, large friction creation or other methods. Now we provide an explanation of essence of an attitude reorientation method. Let's consider conjugate spinup of conjugate rotors 1 and 2 (Fig.1-a) in the absence of external e forces moments M xe M y M ze 0 assuming initial rest of main body and all rotors and mass-inertia symmetry of the system
p 0 q 0 r 0 0; i : i 0 0

(12) (13)

A BC D In simplest case we can use following piecewise constant spinup internal forces moments
M 12 , if t 0, t1s2 M otherwise 0,
i 1

M 12 , if t 0, t1s2 i M2 otherwise 0,

(14)

s where t12 is the time point of spinup stopping of rotors 1 and 2; M12 const 0 . After the conjugate spinup rotors 1 and 2 will reach an absolute value S12 M12 t1s2 of relative angular

velocity 1 S12 , 2 S12 but the main body will remain in rest; angular momentum of system will be equal null as before. After the conjugate spinup we will capture rotor 1 at time point t1c t1c t1s2 and then the relative angular velocity of the rotor 1 will be null 1 0 , but main body will take absolute
angular velocity p and rotor 2 will change relative angular velocity up to 2 . Conservation of angular

momentum of full system makes it possible to write Ap I 2 0 Similarly, conservation of angular momentum of rotor 2 makes it possible to write I p 2 IS12

(15) (16)

Numerical values for angular velocities after caption of rotor 1 are obtained from expressions (15) and (16)
p

At the time t

c 2



c t2 t1c we will capture of rotor 2 and then all of system's bodies (main body and

IS12 AS ; 2 12 A I A I

(17)

both conjugate rotors) will return to the absolute rest. Thus we can conclude that conjugate spinup and two serial captures of conjugate rotors bring to piecewise constant angular velocity of main body
c 0, t 0, t1c t2 , p IS cc P 12 , t t1 , t2 A I

(18)

It can be used for main body angular reorientation about corresponding axis. In our case the main body performed the rotation about Ox axis by finite angle
57


x

c IS12 t2 t1c



A I

(19)

In the next section we will put forward this method of attitude reorientation and solve the equation of motion. 4. Analytical Solutions The research results presented in previous section showed that angular velocity of main body is piecewise constant value at realization of conjugate spinup and serial capture of one pair of conjugate rotors (rotors 1 and 2). Similarly we can explain the same process at realization of conjugate spinup and serial capture of another one pair of conjugate rotors. For example, for the rotors 3 and 4 we will write
c 0, t 0, t3 q IS Q 34 , t BI BS 4 34 ; S34 BI



c t4 ,



cc t3 , t4

(20)

M 34 t

s 34

s where tij is during (ending time point) of spinup of conjugate rotors i and j; tic is time point of rotor i

capture; M ij is absolute value of principle momentum of internal forces for spinup of conjugate rotors i and j; Sij is absolute value of relative angular velocity of conjugate rotors i and j after spinup; i is relative angular velocity of rotor i after capture of its conjugate rotor i*
i
i

t tic* , t

c i



c c c Assuming t1c t3 , t2 t4 with help of (20), (17), and (13) we can prove following auxiliary expression
c at t t1c , t2

p 34 q

12

p 4 q 2 0

(21)

Expression (21) show that gyroscopic term in the third equation (4) equals null at concurrent execution of spinups and captures of conjugate rotors {1, 2} and {3, 4}. Therefore angular velocity r will remain constant. Also we can prove equality to zero of all gyroscopic terms and equality to constant of all angular velocity components at every set of concurrent spinups and captures. Thus, angular velocity components are constant at concurrent execution of spinups and between conjugate captures of rotors (22) p P; q Q; r R Let's assume synchronical capture of all conjugate rotors {1, 2}, {3, 4}, {5, 6} and coincidence of frame Oxyz initial position and fixed frame OXYZ (Fig.2).
c c t1c t3 t5 t c start c c c ; t 2 t 4 t6 t c finish c start

0 t

c start



1; 1 t

c start



2 t

c start



3 t



0

(23)

Now we can solve kinematic equation (10)

58


0 t cos 2 t

t t

c start



Q sin
2 2

2 t t 2
2

;

c start



c t tstart P sin 2 c t tstart R ; 3 t sin 2

1 t


(24)

P Q R ; t t

c start

,t

c finish

; T t

c finish

t

c start

Solutions (24) show that at the time point t by angle (Fig.2)

c finish

the main body performed finite rotation about vector e

P Q R e cos , cos , cos T

T

(25)

So, solutions (24) illustrate an opportunity of practical use of the multirotor system for attitude reorientation of spacecraft. How we can see, attitude reorientation of SC by desired finite rotation was performed with help of conjugate spinup and serial captures of rotors into the SC. Above we have considered the case of motion with zero initial value (12) angular momentum of system with equal inertia moments (13). It is easy to show that without external forces at zero initial value angular momentum of system the main body with different inertia moments also has piecewise constant components of angular velocity in frame Oxyz. First of all, in this case vector of angular momentum identically equals zero in every coordinate frame including moving frame Oxyz connected with the main body
Me 0 K K K
OXYZ Oxyz OXYZ OXYZ Z

const

K X , KY , K


T

T

Kx , K y , Kz 0 0 = - â K

(26)

K X KY K Z K x K y K z 0
Oxyz Oxyz

- â 0 0

Kx K y Kz 0



A, B, C



Then we can write
p Kx 0 q Ky 0 r Kz 0

c if t 0, t1c t2 , c p P const, 1 0, 2 2 if t t1c , t2

0, 1 2 0 S12





c c if t 0, t3 t4 , cc p Q const, 3 0, 4 4 if t t3 , t4

0, 3 4 0 S

34



c c if t 0, t5 t6 , cc r R const, 5 0, 6 6 if t t5 , t6

0, 5 6 0 S

56





Last expression proves that main body has piecewise constant components of angular velocity in case zero value of system angular momentum. The case of main body attitude reorientation at zero value of system angular momentum is inertialess and reaction-free method. This method has no disadvantages of reorientation method with help of reaction wheels system like a negative affect of nonlinearity of internal spinup engines, use of large value torque of internal spinup engines and, therefore, large energy consumption. If initial angular momentum of system equals nonzero vector, then at rotors captures main body
59


angular velocity will depend on time. In this case the effect of "gyroscopic power increase" will take place, which means gyroscopic torque initiation. We can illustrate gyroscopic torque initiation and variability of main body angular velocity
K K
Oxyz OXYZ

const 0

0

0, tic t tic* 0
Oxyz

K

(27)

Oxyz

= - â K K
x, y , z

0

var t



Gyroscopic power increase allow to use of small value of initial main body angular veloc ity (at conjugate rotors capture) to appear large value of torque which accelerate reorientation process. This effect is used for attitude control of spacecraft with control moment gyroscope. In addition, nonzero initial angular momentum can be transferred from main body on several (or one) rotors with help of internal forces moments (by engines): main body comes to a stop, but rotors comes to a spins. We need to note that initiation of nonzero angular momentum can be performed by rocket jet engines. On the contrary, if this initial momentum is nonzero, then it can be reduced to null also by jet engines and internal dampers. So, it is possible to use both modes (with/without initial value of system angular momentum) of realization of attitude method. Systems like multirotor-spider (Fig.1-b) and multirotor-hedghog (Fig.1-c) allow using sequences of serial rotor spinups and captures to perform compound attitude motion. In these cases we can write a program for realization of complex serial spinups and captures of set conjugate rotors. For rotor-type hedgehog this research and engineering problem is highly interesting.

5. Numerical Simulation of Reorientation Process Let's calculate two sets of results of numerical simulation of multirotor system attitude reorientation. We used the following internal forces moment
M ij M
jj *

M ij* M M
jj *

H t H t t s j* H t t c* j const 0; 1; j 1, 3, 5; j* 2, 4, 6
jj*

H t H t t s j H t t

c j


(28)

where H(t) is Heaviside function. In expressions (28) the first term corresponds to piecewise constant spinup moment and the second ­ to capture moment of viscous friction type. System parameters for calculations are presented in table I. Also the following parameters are common for all calculations: = 300 Nms, t s 3 s. Simulation results in Fig.3 correspond to the case of reorientation process with zero value of angular momentum (12) and different inertia moments of main body. These results demonstrate constant components of angular velocity of the body between conjugate captures. Fig.4 corresponds to reorientation process at nonzero system angular momentum
(0) 1 (0) ... 4 (0) 0; 5 (0) 100 1 / s

for main body with the same inertia moments. In this mode of motion we also can see that main body comes to permanent rotation, which illustrate redistribution of angular momentum. Certainly, it is possible to perform the next series of spinups and captures to new reorientation of main body and transfer of angular momentum back to rotors.

60


Table I

A, kg m
Fig.3 Fig.4 60 100
2

B, kg m
80 100
2

C, kg m
100 100
2

I, kg m
10 10
2

t1c , s
4 .0 4 .0

c t2 , s

c t3 , s

c t4 , s

c t5 , s

c t6 , s

M 12 , N m
10 10

M 34 , N m
20 20

M 56 , N m
30 0

5 .5 4 .7 5

4 .5 4 .0

5 .0 4 .7 5

4 .7 5 4 .0

6 .0 4 .7 5

61


Fig.3

Fig.4

6. Conclusion In the paragraph new multirotor systems and their attitude reorientation method have been developed. Systems like multirotor-spider and multirotor-hedghog allow using programs for
62


serial rotor spinups and captures to perform compound attitude motion. The analysis of such motion is a subject of further research. RE
FERENCES

[1] J. Wittenburg, Dynamics of Systems of Rigid Bodies, Teubner, Stuttgart Germany, 1977. [2] P. W. Likins, Spacecraft Attitude Dynamics and Control - A Personal Perspective on Early Developments, J. Guidance Control Dyn. Vol. 9, No. 2 (1986) 129-134. [3] D. L. Mingori, Effects of Energy Dissipation on the Attitude Stability of Dual-Spin Satellites, AIAA Journal, Vol. 7, No. 1 (1969) 20-27. [4] P. M. Bainum, P. G. Fuechsel, D. L. Mackison, Motion and Stability of a Dual-Spin Satellite With Nutation Damping, Journal of Spacecraft and Rockets, Vol. 7, No. 6 (1970) 690-696. [5] K. J. Kinsey, D. L. Mingori, R. H. Rand, Non-linear control of dual-spin spacecraft during despin through precession phase lock, J. Guidance Control Dyn. 19 (1996) 60-67. [6] C. D. Hall, Escape from gyrostat trap states, J. Guidance Control Dyn. 21 (1998) 421-426. [7] C. D. Hall, Momentum Transfer Dynamics of a Gyrostat with a Discrete Damper, J. Guidance Control Dyn., Vol. 20, No. 6 (1997) 1072-1075. [8] A. E. Chinnery, C. D. Hall, Motion of a Rigid Body with an Attached Spring-Mass Damper, J. Guidance Control Dyn. Vol. 18, No. 6 (1995) 1404-1409. [9] W. H. Steyn, A dual-wheel multi-mode spacecraft actuator for near-minimum-time large angle slew maneuvers, Aerospace Science and Technology, Vol. 12, Issue 7 (2008) 545-554. [10] D. Izzo, L. Pettazzi, Command shaping for a flexible satellite platform controlled by advanced fly-wheels systems, Acta Astronautica, Vol. 60, Issues 10-11 (2007) 820-827 [11] V.S. Aslanov, A.V. Doroshin, Two cases of motion of an unbalanced gyrostat, Mechanics of solids, Allerton Press, Inc., Vol. 41, No. 4, (2006). 29-39. [12] V.S. Aslanov, A.V. Doroshin, Stabilization of the descent apparatus by partial twisting when carrying out uncontrolled descent. Cosmic Research, Vol. 40, No. 2 (2002). 193-200. [13] V.S. Aslanov, A.V. Doroshin, G.E. Kruglov, The Motion of Coaxial Bodies of Varying Composition on the Active Leg of Descent, Cosmic Research, Vol.43, No. 3 (2005) 213-221. [14] V. S. Aslanov, A. V. Doroshin, The Motion of a System of Coaxial Bodies of Variable Mass, J. Appl. Math. Mech. Vol.68 (2004) 899-908. [15] V.S. Aslanov, A.V. Doroshin, Influence of Disturbances on the Angular Motion of a Spacecraft in the Powered Section of Its Descent, Cosmic Research, Vol. 46, No. 2 (2008) 166­171. [16] A.V. Doroshin, Evolution of the Precessional Motion of Unbalanced Gyrostats of Variable Structure. Journal of Applied Mathematics and Mechanics 72 (2008) 259­269. [17] A.V. Doroshin, Synthesis of Attitude Motion of Variable Mass Coaxial Bodies. WSEAS TRANSACTIONS on SYSTEMS and CONTROL, Issue 1, Volume 3 (2008) 50-61. Available: http://www.wseas.us/e-library/transactions/control/2008/25-273.pdf [18] A. Pothiawala, M. A. Dahleh, H Optimal Control For The Attitude Control And Momentum Management Of The Space Station, MIT, Cambridge, MA 02139. Available: http://dspace.mit.edu/bitstream/1721.1/3208/1/P-1985-22200134.pdf

63


[19] M. Blundell, D. Harty, Multibody Systems Approach to Vehicle Dynamics, Elsevier, 2004, 288 p. [20] W. Schiehlen, N. Guse, R. Seifried, Multibody dynamics in computational mechanics and engineering applications, Computer Methods in Applied Mechanics and Engineering, Vol. 195, Issues 41-43 (2006) 5509-5522. [21] T. M. Wasfy, A. K. Noor, Multibody dynamic simulation of the next generation space telescope using finite elements and fuzzy sets, Computer Methods in Applied Mechanics and Engineering, Vol.190, Issues 5-7 (2000) 803-824. [22] P. Santini, P. Gasbarri, Dynamics of multibody systems in space environment; Lagrangian vs. Eulerian approach, Acta Astronautica, Vol. 54, Issue 1 (2004) 1-24. [23] H. Goldstein, Classical Mechanics, 2-nd ed. Reading, MA: Addison-Wesley, 1980.

64


PART IV. EXAMPLE OF MAPLE PROGRAM FOR ATTITUDE DYNAMICS OF THE DUAL-SPIN-SPACECRAFT MODELING
The following listing of MAPLE-program can be used for attitude dynamics of dual-spinspacecraft. > restart; > with(plots): Warning, the name changecoords has been redefined > ur1:=A*diff(p(t),t)+(C2-B)*q(t)*r(t)+C1*(r(t)+sigm(t))*q(t)=0; d ur1 := A p( t ) ( C2B ) q( t ) r( t )C1 ( r( t )sigm( t ) ) q( t )0 dt > ur2:=B*diff(q(t),t)+(A-C2)*p(t)*r(t)-C1*(r(t)+sigm(t))*p(t)=0; d ur2 := B q( t ) ( AC2 ) p( t ) r( t )C1 ( r( t )sigm( t ) ) p( t )0 dt > ur3:=C2*diff(r(t),t)+C1*diff((r(t)+sigm(t)),t)+(B-A)*p(t)*q(t)=0; d d d ur3 := C2 r( t ) C1 r( t ) sigm( t ) ( BA ) p( t ) q( t )0 dt dt dt > ur4:=C1*diff((r(t)+sigm(t)),t)=eps*sin(nu*t); d d ur4 := C1 r( t ) sigm( t ) eps sin( t ) dt dt > #A2:=4;B2:=6; C2:=8;A1:=1;C1:=2;# - OBLATE GYROSTAT > A2:=15;B2:=8; C2:=6;A1:=5;C1:=4;# - PROLATE GYROSTAT A2 := 15

B2 := 8 C2 := 6
A1 := 5

C1 := 4
> eps:=0;nu:=1;

eps := 0
:= 1

> A:=(A1+A2);C:=C2+C1;B:=B2+A1;

A := 20 C := 10
B := 13

> DELTA:=10;

DELTA := 10

> p0:=0.15;q0:=0.15;r0:=0.1;sigm0:=DELTA/C1-r0; #evalf(sqrt(p0^2*((B2+A1)*(A1+A2)(A1+A2)^2)/(C2*(C2-B2-A1))))+0.5;sigm0:=-r0;
65


p0 := 0.15 q0 := 0.15

r0 := 0.1
sigm0 := 2.400000000
> resh:=dsolve({ur1,ur2,ur3,ur4, p(0)=p0,q(0)=q0,r(0)=r0,sigm(0)=sigm0},{p(t),q(t),r(t),sigm(t)},type=numeric,output=listproced ure); resh := [ t( p ro c(t) ... e nd p ro c , p( t )( p ro c(t) ... e nd p ro c , ) ) q( t )( p ro c(t) ... e nd p ro c , r( t )( p ro c(t) ... e nd p ro c , ) ) sigm( t )( p ro c(t) ... e nd p ro c ] ) > TT:=30; > > > > >

TT := 30

pic_ch_1:=odeplot(resh,[t,p(t)],0..TT,numpoints=1500): pic_ch_2:=odeplot(resh,[t,q(t)],0..TT,numpoints=500): pic_ch_3:=odeplot(resh,[t,r(t)],0..TT,numpoints=500): b:=evalf(sqrt(p0^2*(C2-A2-A1)*(A1+A2)/(C2-B2-A1)/(B2+A1))); b := 0.2631174058

> lambda:=evalf(sqrt(p0^2*(B2-A2)*(C2-A1-A2)/(C2)/(B2+A1))); := 0.1681345615 > pp:=t->p0/cosh(lambda*t);
pp := t p0 cosh( t )

> qq:=t->b*tanh(lambda*t); > rr:=t->r0/cosh(lambda*t);

qq := tb t anh( t )
rr := t r0 cosh( t )

> > > > >

pic_an_1:=plot(pp(t),t=0..TT,style=point,color=blue): pic_an_2:=plot(qq(t),t=0..TT,style=point,color=blue): pic_an_3:=plot(rr(t),t=0..TT,style=point,color=blue): display([pic_ch_1,pic_an_1]);

66


> display([pic_ch_2,pic_an_2]);

> display([pic_ch_3,pic_an_3]);

>> > urfi:=diff(fi_(t),t)-r(t)+(cos(tetta(t))/sin(tetta(t)))*(p(t)*sin(fi_(t))+q(t)*cos(fi_(t))); d cos( t et t a( t ) ) ( p( t ) sin( fi_ ( t ) )q( t ) cos( fi_ ( t ) ) ) urfi := fi_ ( t ) r( t ) dt sin( t et t a( t ) ) > urteta:=diff(tetta(t),t)-p(t)*cos(fi_(t))+q(t)*sin(fi_(t)); d urteta := t et t a( t ) p( t ) cos( fi_ ( t ) )q( t ) sin( fi_ ( t ) ) dt > urpsie:=diff(psie(t),t)-(1/sin(tetta(t)))*(p(t)*sin(fi_(t))+q(t)*cos(fi_(t))); d p( t ) sin( fi_ ( t ) )q( t ) cos( fi_ ( t ) ) urpsie := p sie( t ) dt sin( t et t a( t ) ) > > urdel:=diff(del(t),t)=sigm(t);
urdel :=

d del( t )sigm( t ) dt

> # Initial conditions ­ case of coincidence of vector K and axis OZ. > K:=sqrt(A^2*p0^2+B^2*q0^2+(C*r0+C1*sigm0)^2); K := 11.18760475 > Kz:=C*r0+C1*sigm0; >

Kz := 10.60000000
67

> tetta0:=arccos(Kz/K);fi0:=evalf(Pi/2);psi0:=0;


tetta0 := 0.3255431246 fi0 := 1.570796327

:= 0
F:=dsolve({ur1,ur2,ur3,ur4,urteta,urfi,urpsie,urdel,p(0)=p0,q(0)=q0,r(0)=r0,sigm(0)=sigm0,tetta( 0)=tetta0,fi_(0)=fi0,psie(0)=psi0,del(0)=0},{p(t),q(t),r(t),sigm(t),psie(t),tetta(t),fi_(t),del(t)},type =numeric,output=listprocedure); F := [ t( p ro c(t) ... e nd p ro c , del( t )( p ro c(t) ... e nd p ro c , ) ) fi_ ( t )( p ro c(t) ... e nd p ro c , p( t )( p ro c(t) ... e nd p ro c , ) ) p sie( t )( p ro c(t) ... e nd p ro c , q( t )( p ro c(t) ... e nd p ro c , ) ) r( t )( p ro c(t) ... e nd p ro c , sigm( t )( p ro c(t) ... e nd p ro c , ) ) t et t a( t )( p ro c(t) ... e nd p ro c ] ) pp:=subs(F,p(t));qq:=subs(F,q(t));rr:=subs(F,r(t));SIGM:=subs(F,sigm(t));FIe:=subs(F,fi_(t));PPs ie:=subs(F,psie(t));TET:=subs(F,tetta(t));DEL:=subs(F,del(t)); pp := p ro c t) ... end p ro c (

qq := p ro c t) ... end p ro c (
rr := p ro c t) ... end p ro c (

SIGM := p ro c(t) ... end p ro c FIe := p ro c t) ... end p ro c (
PPsie := p ro c(t) ... end p ro c

TET := p ro c(t) ... end p ro c DEL := p ro c(t) ... end p ro c
> plot([pp(t),qq(t),rr(t),SIGM(t)],t=-0..TT);

> plot([PPsie(t),FIe(t),DEL(t)],t=0..TT);
68


> plot(TET(t),t=0..2*TT,numpoints=4000);

> fw:=fopen(`angle_.txt`,APPEND);

fw := 0

> h_time:=0.01;NNN:=round(2*TT/h_time); h_time := 0.01

NNN := 6000
> for i from 0 by 1 to NNN do t:=i*h_time: fprintf(fw,"%f %f %f %f %f\n",t,PPsie(t),TET(t),FIe(t),DEL(t)): od: > close(fw);
69