NONLINEAR STABILITY ANALYSIS OF RELATIVE EQUILIBRIA OF A SOLID CARRYING A MOVABLE POINT MASS IN THE CENTRAL GRAVITATIONAL FIELD Текст научной статьи по специальности «Физика»

Russian Journal of Nonlinear Dynamics, 2019, vol. 15, no. 4, pp. 505-512. Full-texts are available at http://nd.ics.org.ru DOI: 10.20537/nd190409

MSC 2010: 70H14, 70H15, 70M20

Nonlinear Stability Analysis of Relative Equilibria of a Solid Carrying a Movable Point Mass in the Central Gravitational Field

O. V. Kholostova

The motion of a solid (satellite) carrying a moving point mass in the central Newtonian gravitational field in an elliptical orbit of arbitrary eccentricity is considered. The law of motion of a point mass is assumed to allow for the existence of relative equilibria of the "body-point" system in the orbital coordinate system. A nonlinear stability analysis of these equilibria is carried out, based on the construction and normalization of the area-preserving mapping generated by the motions of the system.

Keywords: solid carrying a point mass, elliptical orbit, relative equilibrium, stability, resonance

1. Introduction

Problems of stabilization and motion control of a dynamically symmetric body (satellite) with a point mass, which moves according to a chosen law, are studied in [1]. In [2], for different laws of motion of a point mass, dynamics of the carrying solid (satellite) relative to the center of mass in a circular and elliptical orbit is considered. In particular, the law of motion of a point mass is found, which allows for relative equilibria of the satellite in an elliptical orbit; an investigation of their linear stability is carried out. In [3], approximate autonomous equations of motion of a satellite in a circular orbit are obtained for the cases of fast motions of a carried point; and a stability analysis of relative equilibria of the body in the orbital coordinate system is carried out. The dynamics of a spacecraft with a variable mass geometry in an elliptical orbit is studied in [4, 5], under the assumption that the change of the three principal central moments of inertia allows for the existence of 24 relative equilibria of the spacecraft in the orbital coordinate system. A linear analysis of their stability is carried out.

Received May 29, 2019 Accepted August 01, 2019

This research was supported by the Russian Foundation for Basic Research (project No. 17-01-00123).

Olga V. Kholostova [email protected]

Moscow Aviation Institute (National Research University) Volokolamskoe sh. 4, GSP-3, A-80, Moscow, 125993 Russia

In this paper, we continue the study of relative equilibria of a solid (satellite) carrying a point mass in an elliptical orbit, which was begun in [2]. A nonlinear analysis of their stability in the linear stability regions as well as on the boundary curves of stability and instability regions is carried out.

2. Equations of motion of a body-point system

Consider the motion of a solid body (satellite) of mass M in the central Newtonian gravitational field with an attracting center O* of mass M*. We assume that the body contains a material point P of mass m, moving in a given way along a curve fixed in the body. Let Oxyz be the coordinate system with the origin at the center of mass O of the satellite and the axes directed along its principal central inertia axes. The projections of the vector OP on the axis of the coordinate system Oxyz are denoted by x(t), y(t), z(t), these functions being considered to be twice continuously differentiable in t.

We introduce the orbital coordinate system CXCYCZC with the origin at the center of mass C of the body-point system and the axes CXC, CYC and CZC directed, respectively, along the transversal to the orbit of C, binormal to it, and the radius vector of the point C relative to the attracting center. We also introduce the coordinate system OXOYOZO obtained from the orbital one by translation. The orientation of the body-fixed coordinate system Oxyz relative to OXOYOZO will be specified using the Euler angles d, defined in a usual way.

If the characteristic size of the system is small compared to the orbital size, then the usual [2, 6] assumptions can be accepted about the independence of the motion of the system relative to its center of mass from the motion of the center of mass itself. The orbit of the latter is assumed to be elliptic with eccentricity e (0 < e < 1).

The motion of the system relative to the center of mass can be described by differential equations of dynamic Euler equation type [2]

[Jx + n(y2 + z2)]wx — fixyuy — ¡ixzwz + (Jz — Jy)wywz + (2.1)

+M(x^x + yUy + z^z )(yWz — z^y) + 2(yy + zz )Wx — —2x(yUy + zUz) + yz — zy] = Mx,

3k 3

MX = ~ Jy)a32a33 + + 0-322/ + a33z)(a33y - a32z)}

C

(x, y, z; a3!, (I32, a33),

where [7]

n n2 k n2

mM

^ = -77) K = 7 -M*.

r m + M' ' *

(2.2)

Here, Jx, Jy, Jz

are the principal central moments of inertia of the body; wx, wy, wz; a21, a22, a23; a31, a32, a33 and Mx, My, Mz are projections on the Ox, Oy, Oz axes of the angular velocity vector of the body, the unit vector of the axis OYO, the unit vector of the axis OZO, and the gravitational moment vector, respectively. The symbol 7 denotes the universal gravitational constant, v the true anomaly of the mass center, n = 2n/T where T is the orbital period of the center of mass of the system.

3. Equations of planar motion. Relative equilibria

Let the carried point P move along the principal body inertia axis Oy (that is, x(t) = z(t) = 0). In this case, Eqs. (2.1) have the particular solution ^ = n, d = n/2, and the change in the angle y is determined by equation [2]

jt\(Jz +M2)((P + v)\ + ~ Jy + ^y2)sinycosy = 0.

C

(3-1)

This equation describes the planar motion of the solid, for which its principal inertia axis Oz is perpendicular to the orbital plane.

If a moving point is absent, then in a circular orbit there are relative equilibrium positions of the body in the orbital coordinate system, corresponding to the particular solutions y = 0 and y = n/2 of Eq. (3.1). In an elliptical orbit, the relative equilibria of the satellite do not exist in the absence of the point P. When a moving point is added, such equilibria can appear only when a special law of its motion is chosen. It is easy to verify that Eq. (3.1) admits the particular solutions y = 0 and y = n/2 when the relation y = y(v) of the point motion law written as a function of v is given by the expression [2]

2

Jz + ¡iy2(v) = J + fiy (0)]

1 + e

1 + e cos v

(3-2)

Note that, for the first solution (y = 0), the body principal inertia axes Ox and Oy are parallel to the transversal to the orbit of mass center C of the system and the vector O^C, respectively. For the second one, (y = n/2), these axes are parallel to the vector O^C and transversal to the orbit of the point C, respectively, and the carried point mass P moves in a straight line perpendicular to the vector O^C.

Let us introduce the new independent variable v in Eq. (3.1) using the relations (2.2) and taking into account the expression (3.2). The equation of planar oscillations of the satellite with the chosen law of motion of the point P has the form (the prime indicates the differentiation with respect to v)

J' +

(s — 3)(1 + e cos v )

(1 + e)2

+

3

1 + e cos v

sin y cos y = 0,

s = 3

Jx-Jy + Aty2(0) Jz+(iy2( 0) ;

(3-3)

-3 < s < 3-

Equation (3.3) is equivalent to a system of two first-order differential equations written in the form of Hamiltonian canonical equations, in which the generalized momentum is introduced using the formula pv = y', and the Hamiltonian function has the form

H

1 2

r2Pv +

(s — 3)(1 + e cos v)

(TTW2

+

3

1 + e cos v

sin2 y.

(3-4)

Assuming y = y0 + x, = y (y0 = 0 or n/2) in (3.4), we consider the Hamiltonian of perturbed motion as follows:

H'2 2V ^ 2

H = H2 + H4 + --(s — 3)(1 + e cos v)

H4 = T

(1 + e)2 (s — 3)(1 + e cos v)

(TTë)2

+

3

1 + e cos v

x2,

(3-5)

+

3

1 + e cos v

x4.

Here, the ellipsis means the set of terms of even degrees (the sixth and above) with respect to perturbations.

A linear stability analysis of particular solutions p = po, = 0 of the system with the Hamiltonian (3.5) has been carried out numerically and analytically in [2]. The results are illustrated in Fig. 1 and Fig. 2 in the plane of the parameters e, s for the cases po = 0 and po = n/2, respectively. In these figures, the instability regions are shaded; the linear stability regions are white. Note that the solution po = n/2 is always unstable if 0 ^ s ^ 3 and 0 < e < 1; this part of the parameter plane in Fig. 2 is not shown. In this region and in the shaded regions in Figs. 1,2, instability takes place not only in the linear problem but in the complete nonlinear system.

The boundary curves of stability and instability regions in Fig. 1 originate from the points of the axis Oe with the coordinates s = 0, 1/4, 9/4, as well as from points with s = 4, 25/4 lying outside the range of admissible values of s. The generating points for the boundaries of the regions in Fig. 2 are the points s = 0, —1/4, —1, —9/4 of the axis Oe. From each of the points indicated (except for the point (0, 0) in both Fig. 1 and Fig. 2), two boundary curves are initiated; from the points (0, 0), one curve. Note that for small values of e, the equations of the boundary curves going out from the point s = —9/4 coincide in terms up to the order of e2 inclusive [2]; in Fig. 2 these curves are indistinguishable. On all boundary curves, the solutions considered are linearly unstable.

-3-2-10 1 2 s

Fig. 1. The case p0 = 0.

1

0.8-

0.6- \ e \ 0.4

0.2-

Fig. 2. The case p0 = n/2.

The region above the upper boundary curve in Fig. 1 consists of a countable set of stability and instability regions, in [2] and in this paper this part of the plane is not considered.

The purpose of this paper is to carry out a nonlinear stability analysis of the two considered relative equilibria of the system (with respect to p and ) in the linear stability regions, as well as on the boundary curves separating the stability and instability regions.

4. The algorithm of nonlinear analysis

To draw conclusions on stability in a complete nonlinear formulation, one should normalize the Hamiltonian function (3.5) in terms up to the fourth order inclusive with respect to perturbations, and then apply the known criteria for stability analysis of nonlinear, time-periodic one-degree-of-freedom Hamiltonian systems (see, for example, [8]). From the computational point of view, it is more convenient to normalize not the Hamiltonian function itself, but the

area-preserving mapping generated by the motions of the Hamiltonian system under study over a time interval equal to period 2n. We have normalized this mapping using the algorithm developed in the paper [9]. Here is a brief summary of the relevant theoretical results outlined in this article.

Let X(v) = Wxj(v)||2j=1 be the fundamental matrix of solutions of a linear system with the Hamiltonian H2 from Eq. (3.5), with initial conditions X (0) = E where E is a second-order unit matrix. Consider the characteristic equation of the matrix X (2n):

p2 — 2Ap + 1 = 0, 2A = x11(2n)+ x22(2n). (4.1)

At the first stage, a linear normalization is carried out. It differs for linear stability regions and boundary curves (see Sections 4.1 and 4.2 below).

Then a near-identity canonical transformation x, y — x,,y, is made. It is represented as

x, = x + Xs(x, y) + ..., y, = y + Y3(x, y) + ... (4.2)

and given by the relations

dF dF _

x = —, y* = T—, F = x#y + FA(x„ y) + ..., F4= frsxlys. (4.3)

°x* r+s=4

Here F(x, X) is the generating function of this transformation. Its coefficients are chosen in such a way as to obtain the simplest form in the nonlinear part of the mapping. For our purposes, it is enough to obtain the normal form of the function F4 with the fourth-degree terms. The criteria of stability and instability can be expressed in terms of the coefficients of this function.

4.1. Linear stability regions

Consider the regions of linear stability in Figs. 1 and 2. In these regions, |A| < 1, the roots P1 ,2 = cos a ± i sin a (cos a = A) of Eq. (4.1) are complex conjugate. The matrix X (2n) has the form

(a b\ \c d).

A canonical transformation (with valence a = sign b) of the form

/ Ibl _ a — cos(a) _ /sin(a)_

x = —a\/ —t^x, V = — -x ~ -i / —m—V

y sin(a) ^J\b\ sin(a) V \b\

reduces the matrix of the linear part of the mapping to the normal form

(cos a sin a\ cos a sin a .

— sin a cos a J

After that, we normalize the function F4 from the relation (4.3) and restore the form of the normalized Hamiltonian in terms to the fourth degree inclusive with respect to perturbations. One should distinguish between the nonresonant case (0 < |A| < 1) and the fourth-order resonance case (A = 0).

In the nonresonant case, the normalized Hamiltonian written in the symplectic polar coordinates tp, r (x* = \/2rsint£>, y* = \/2rcosf£>) has the form

H = Xr + c2 r2 + O(r5/2).

Here the coefficients X and c2 are defined from the formulas

cos 2ttX = A, c2 = -¿(/22 + 3/40 + 3/04)- (4.4)

If c2 = 0, then, according to the Arnold-Moser theorem, the solution under study is stable. If c2 = 0 (the degeneracy case), then to solve the question of stability, the Hamiltonian should be normalized in terms of higher than the fourth order.

In the fourth-order resonance case, the normalized Hamiltonian has the form

H = Xr + [c2 + a2 sin(4^ - 4Xt) + p2 cos(4^ - 4Xt)]r2 + O(r5/2), (4.5)

in which X = (2n + 1)/4 (n is an integer), the coefficient c2 is given in (4.4), and

«2 = -¿(/13 - /31), 02 = -¿(-/22 + /40 + /04)- (4.6)

If the condition holds

c2 >a2 + &2, (4.7)

then the solution under study is stable. If the inequality is satisfied with the opposite sign, the solution is unstable.

4.2. Stability criteria on boundary curves

On the boundary curves, the condition |A| = 1 is fulfilled and p1 = p2 = 1 (first-order resonance) or p1 = p2 = —1 (second-order resonance). The matrix X(2n) is represented by one of the following types:

/±1 b\ /±1 0\

(0 ±J or (c ±1), (48)

where the upper and lower signs correspond to the cases A = 1 and A = —1.

If the matrix X (2n) is represented by the first option from (4.8), then we make a linear canonical change (with valence sign b) of the form

/—sign b _ x = \/\b\x, y = ^=y.

If the matrix is represented by the second option, then the canonical change (with valence — sign c) has the form

x = —j=, y = signed. |c|

As a result, the matrix of the linear part of the mapping is reduced to the normal form

Co1 ±1), (49)

where the upper and lower signs correspond to the upper and lower signs in (4.8).

At the stage of nonlinear normalization, we obtain the normal form of the function F4 of the fourth degree terms in the generating function. The conclusion on stability of the solution considered on the boundary curves depends on the sign of the coefficient f40 [9]. If A = 1 (upper signs in matrix (4.9)), then the solution is stable for f40 < 0, and unstable for f40 > 0. If A = —1 (lower signs in (4.9)), then, on the contrary, we have stability for f40 < 0, and instability for f40 < 0.

5. Results of nonlinear analysis

The examination of the criteria presented in the previous section has been carried out numerically and analytically using the computational formulas given in [9]. The calculations were performed in the system of analytical calculations MAPLE. Let us describe the results obtained.

5.1. Conclusions on stability in the linear stability regions

In the absence of the fourth-order resonance, a linear and then nonlinear normalization of the mapping has been performed in the considered linear stability regions as it is described in Section 4.1. The values of the coefficient c2 were examined in these regions. It is revealed that in each region the sign of this coefficient is preserved, that is, there are no degeneracy cases. Thus, in nonresonant cases, the solutions under study are stable.

Consider the resonance cases. In each of the investigated linear stability regions, there is one fourth-order resonance curve. For the regions in Fig. 1, these curves originate at the points s = 1/16, 9/16, 25/16 of the axis e = 0; two more resonance curves originate at points s = 49/16 and 81/16 of this axis, which are outside the permissible range of parameter values. For the regions in Fig. 2, the generating points of the resonance curves have the coordinates s = —1/16, -9/16, —25/16.

For small values of e, analytical equations of resonant curves are obtained in the form of series in e:

1

S = l6" _9_

16 25

16

2 9601 3 160292701 4 e - -TTT—e + e4 +

s = — —

192 6609 3

s =---

s = — — —

s = ——--

1

Î6 9_

Î6 25

8~

47 10729 —c —-

8 ' 384

39 5049 2

—e H--e H--e" —

8 640 320

23 2855 2 1009 3

—e--e H--eà -

8 2688 1344

49 7657 2 8833 3

—e H--e2 H--e

8 384 192

57 2853 2 2397 —e--e--

8 128 64

122880 1508149383 4

7168000 179005699 4

e4 +

e3-

227598336

44891257 4 -p

24576 44891257 4

e4 + ... ; e4 + ...

e4 + ....

73

—p —-

4 1344

30937 2 18673

e -

672

24576

3 893279741 4

e--e +

113799168

For arbitrary values of e, the resonance curves were constructed numerically. They are shown in Figs. 1 and 2 by dotted lines. The coefficients c2, a2 and (32 of the normalized Hamiltonian (4.5) were calculated using formulas (4.4), (4.6), and criterion (4.7) was examined. It is revealed that on all fourth-order resonance curves this criterion is satisfied, and, thus, in the resonance cases, the solutions under study are stable.

5.2. Conclusions on stability on the boundary curves

On the boundary curves going out from the points s = 0, ±1 (Figs. 1 and 2) and the (external) point s = 4 (Fig. 1), we have A = 1, on the other curves (when s = ±1, ±9/4), the equality A = —1 holds. For the points of the boundary curves, the corresponding coefficients f40 of the generating function were calculated, and the criterion described in Section 4.2 was applied.

The result is as follows. For both curves (in Figs. 1 and 2) initiated at the point (0, 0), the instability of the corresponding solutions holds. For the solution = 0, on all left boundaries of the instability regions (see Fig. 1) we have stability, on all right boundaries, instability. For the solution = n/2, on the contrary, on all left and all right boundaries of the instability regions (see Fig. 2) there is instability and stability, respectively.

References

[1] Aslanov, V. S. and Bezglasnyi, S. P., Gravitational Stabilization of a Satellite Using a Movable Mass, J. Appl. Math. Mech, 2012, vol.76, no. 4, pp. 405-412; see also: Prikl. Mat. Mekh, 2012, vol.76, no. 4, pp. 563-573.

[2] Markeev, A. P., Dynamics of a Satellite Carrying a Point Mass Moving about It, Mech. Solids, 2015, vol. 50, no. 6, pp. 603-614; see also: Izv. Akad. Nauk. Mekh. Tverd. Tela, 2015, no. 6, pp. 3-16.

[3] Markeev, A. P., Equations of Motion of a Solid in a Circular Orbit at Fast Relative Motion of Its Carried Material Point, Dokl. Phys, 2016, vol.61, no.4, pp. 184-187; see also: Dokl. Akad. Nauk, 2016, vol.467, no.4, pp.414-417.

[4] Burov, A. A. and Kosenko, 1.1., Planar Vibrations of a Solid with Variable Mass Distribution in an Elliptic Orbit, Dokl. Phys., 2011, vol.56, no. 10, pp. 548-552; see also: Dokl. Akad.. Nauk, 2011, vol.440, no. 6, pp. 760-764.

[5] Burov, A., Guerman, A., and Kosenko, I., Satellite with Periodical Mass Redistribution: Relative Equilibria and Their Stability, Celestial Mech. Dynam. Astronom., 2019, vol. 131, no. 1. Art. 1, 12pp.

[6] Beletskii, V. V., Motion of an Artificial Satellite about Its Center of Mass, Jerusalem: Israel Program for Scientific Translations, 1966.

[7] Duboshin, G.N., Celestial Mechanics: Basic Problems and Methods, Dayton, Ohio: Foreign Technology Division, 1969.

[8] Markeev, A. P., Libration Points in Celestial Mechanics and Space Dynamics, Moscow: Nauka, 1978 (Russian).

[9] Markeyev, A. P., A Method for Analytically Representing Area-Preserving Mappings, J. Appl. Math. Mech, 2014, vol.78, no. 5, pp. 435-444; see also: Prikl. Mat. Mekh., 2014, vol.78, no. 5, pp. 612-624.