Исследование эволюции режима закрутки спутника в плоскости орбиты
( A study of the Evolution of Satellite's Rotations in the Orbit Plane
Preprint, Inst. Appl. Math., the Russian Academy of Science)

1. Введение. В [1] исследована закрутка в плоскости орбиты низколетящего искусственного спутника Земли. Спутник вращался вокруг своей главной центральной оси минимального момента инерции, совершавшей малые колебания относительно нормали к плоскости орбиты; угловая скорость вращения вокруг продольной оси в несколько раз превышала орбитальную угловую скорость. В уравнениях движения спутника относительно центра масс учитывались гравитационный и восстанавливающий аэродинамический моменты. Эти уравнения содержали малый параметр, характеризующий материальную асимметрию спутника и аэродинамический момент. Методом малого параметра и численно была исследована двумерная интегральная поверхность уравнений движения, описывающая невозмущенные движения спутника в режиме закрутки.

Существенным в [1] было предположение относительно формы внешней оболочки спутника. Она считалась эллипсоидом, главные оси которого получены из главных центральных осей инерции спутника сдвигом по оси минимального момента инерции. В результате выражения для компонент аэродинамического момента оказались специфическими четными и нечетными функциями угловых переменных так, что указанная интегральная поверхность состояла из периодических решений.

В данной работе рассматривается по существу тот же спутник, что и в [1], но при произвольном взаимном расположении главных центральных осей инерции и главных осей эллипсоида-оболочки. Интегральная поверхность невозмущенного движения спутника в режиме закрутки теперь уже не состоит из периодических решений. Более того, задающие ее ряды по степеням малого параметра являются формальными и могут служить лишь для приближенного описания представляющих интерес решений на ограниченных, хотя и продолжительных интервалах времени. Тем не менее, исследование интегральной поверхности удается свести к численному исследованию периодической краевой задачи для некоторой вспомогательной системы дифференциальных уравнений. Показана также возможность исследования этой поверхности посредством минимизации специального квадратического функционала.

2. Уравнения движения. Спутник будем считать твердым телом, центр масс которого - точка O - движется по неизменной круговой орбите вокруг Земли. Для записи уравнений движения спутника относительно центра масс введем три правые декартовы системы координат.

Орбитальная система OX₁X₂X₃ связана с плоскостью орбиты спутника и геоцентрическим радиусом-вектором точки O. Ось OX₃ направлена вдоль этого радиуса-вектора, ось OX₂ направлена по вектору кинетического момента орбитального движения спутника.

Система Ox₁x₂x₃ образована главными центральными осями инерции спутника. Оси Ox₁ и Ox₂ отвечают минимальному и максимальному моментам инерции спутника.

Система Oy₁y₂y₃ связана с внешней оболочкой спутника и используется при расчете действующего на него аэродинамического момента. Она занимает неизменное положение относительно системы Ox₁x₂x₃.

Ориентацию системы координат Ox₁x₂x₃ относительно орбитальной системы зададим с помощью углов y, q и j, которые определяются следующим образом. Система OX₁X₂X₃ может быть переведена в систему Ox₁x₂x₃ тремя последовательными поворотами: 1) на угол y вокруг оси OX₃, 2) на угол q вокруг новой оси OX₂, 3) на угол j вокруг новой оси OX₁, совпадающей с осью Ox₁. Введенные углы имеют следующий смысл: q - угол между осью Ox₁ и плоскостью орбиты (плоскостью OX₁X₃), y - угол между осью OX₁ и проекцией оси Ox₁ на плоскость орбиты, j - угол поворота спутника вокруг оси Ox₁. Матрицу перехода от системы Ox₁x₂x₃ к системе OX₁X₂X₃ обозначим ||a_ij||  ³_i,j=1 , где a_ij - косинус угла между осями OX_i и Ox_j. Элементы этой матрицы выражаются через углы y, q и j с помощью формул

a11=cosycosq , a12=-sinycosj+cosysinqsinj , a13 = sinysinj+cosysinqcosj ,

a21=sinycosq , a22=cosycosj+sinysinqsinj , a23=-cosysinj+sinysinqcosj ,        a31=-sinq , a32=cosqsinj , a33=cosqcosj.

Матрицу перехода от системы Ox₁x₂x₃ к системе Oy₁y₂y₃ обозначим ||b_ij ||  ³_i,j=1 , где b_ij - косинус угла между осями Oy_i и Ox_j. Элементы этой матрицы выражаются в функции углов g_c, a_c и b_c, на которые надо повернуть систему Oy₁y₂y₃ последовательно вокруг осей Oy₂, Oy₃ и Oy₁, чтобы перевести ее в систему Ox₁x₂x₃. Соответствующие выражения имеют вид

b11=cosaccosbc, b12=sinacsingc-cosacsinbccosgc, b13=sinaccosgc+cosacsinbcsingc,        b21=sinbc, b22=cosbccosgc, b23=-cosbcsingc,

b31=-sinaccosbc, b32=cosacsingc+sinacsinbccosgc, b33=cosaccosgc-sinacsinbcsingc.

В уравнениях движения спутника относительно центра масс будем учитывать гравитационный и восстанавливающий аэродинамический моменты. Компоненты гравитационного момента в системе Ox₁x₂x₃ имеют вид [2]

Mg1=3w02(I3-I2)a32a33 ,   Mg2=3w02(I1-I3)a33a31 ,

Mg3=3w02(I2-I1)a31a32 .

Здесь w₀ - среднее движение спутника (орбитальная частота), I_i - моменты инерции спутника относительно осей Ox_i.

При выводе выражений для аэродинамического момента внешнюю оболочку спутника считаем эллипсоидом, который в системе координат Oy₁y₂y₃ задается уравнением

(y1-D1)2 L12 +  (y2-D2)2 L22 +  (y3-D3)2 L32 = 1 .

Атмосферу полагаем неподвижной в абсолютном пространстве, а ее плотность r вдоль орбиты спутника постоянной. Пусть спутник неподвижен относительно набегающего аэродинамического потока и молекулы воздуха при столкновении со спутником испытывают абсолютно неупругий удар. Тогда компоненты аэродинамического момента в системе Ox₁x₂x₃ имеют вид

Ma1=rv2S(a12d3-a13d2) ,    Ma2=rv2S(a13d1-a11d3) ,

Ma3=rv2S(a11d2-a12d1) ,    S=pL1L2L3   ж Ц  a12 L12 +  a22 L22 +  a32 L32    ,

di= 3 е k=1  Dkbki ,   ai= 3 е k=1  bika1k     (i=1,2,3) .

Здесь v - геоцентрическая скорость точки O. Эти формулы пригодны и в случае, когда угловая скорость спутника не очень велика по сравнению с w₀.

Исходные уравнения движения спутника относительно центра масс возьмем в виде динамических уравнений Эйлера для компонент W_i абсолютной угловой скорости спутника в системе координат Ox₁x₂x₃ и кинематических уравнений для углов y, q и j. Компоненты угловой скорости будем измерять в единицах w₀, в качестве единицы измерения времени t примем w₀^-1. Орбитальный период при этом равен 2p. Уравнения движения имеют вид

Ч j   =W1+(W2sinj+W3cosj)tg q-  siny cosq  ,

Ч q   =W2cosj-W3sinj-cosy ,

Ч y   =  W2sinj+W3cosj cosq -tg qsiny ,

Ч W   1  =m(W2W3-3a32a33)+eS(a12d3-a13d2) ,

Ч W   2  =  1-l 1+lm (W1W3-3a31a33)+  elS 1+lm (a13d1-a11d3) ,

Ч W   3  =-(1-l+lm)(W1W2-3a31a32)+elS(a11d2-a12d1) ,

l =  I1 I3  ,    m =  I2-I3 I1  ,   e =  rv2 I1w02  .

Здесь точкой обозначено дифференцирование по времени. По своему физическому смыслу e ≥ 0, параметры l и m должны удовлетворять неравенствам |m| < 1, 0 < l < 2/(1-m). Неравенства для l и m следуют из "неравенств треугольника"  для моментов инерции спутника. Правые части выписанных уравнений являются 2p-периодическими функциями j.

3. Закрутка спутника в плоскости орбиты. Интегральная поверхность ориентированного движения. Предположим, что спутник близок к осесимметричному с осью материальной симметрии Ox₁ и что влияние на него аэродинамического момента мало. Иными словами, |m| << 1, eE << 1, E=maxL_iL_j|d_k| (i,j,k=1,2,3). Для упрощения формул оставим только один малый параметр m, приняв e = me₁, e₁E ~ 1. Вместо переменных W₂,W₃ введем переменные

w2=W2cosj-W3sinj ,   w3=W2sinj+W3cosj .

Величины w₂,w₃ представляют собой проекции абсолютной угловой скорости спутника на оси Резаля, получающиеся из осей Ox₂ и Ox₃ при j = 0.

В новых переменных уравнения движения спутника принимают вид

Ч j   =W1+w3tg q-  siny cosq  ,

Ч W   1  =m[W2W3-3a32a33+e1S(a12d3-a13d2)] ,

Ч q   =w2-cosy ,    Ч y   =  w3 cosq -tg qsiny ,

(1)

Ч w   2  =- ж и lW1+w3tg q-  siny cosq ц ш w3+3(1-l)sinqcosq+mQq ,

Ч w   3  = ж и lW1+w3tg q-  siny cosq ц ш w2+mQy .

Здесь

Qq=Q2cosj-Q3sinj ,   Qy=Q2sinj+Q3cosj ,

Q2=  l 1+lm [-(1-l)(W1W3-3a31a33)+e1 S(a13d1-a11d3)] ,

Q3=-l[W1W2-3a31a32-e1S(a11d2-a12d1)] ,

причем переменные W₂, W₃ должны быть выражены через w₂, w₃. Правые части уравнений (1) периодически зависят от j с периодом 2p и аналитически зависят от m в окрестности точки m = 0.

При m = 0 (ось Ox₁ является осью материальной симметрии спутника, и аэродинамический момент отсутствует) уравнения (1) принимают вид

Ч j   =W1+w3tg q-  siny cosq  ,     Ч W   1  =0 ;

(2)

Ч q   =w2-cosy ,     Ч y   =  w3 cosq -tg qsiny ,

Ч w   2  =- ж и lW1+w3tg q-  siny cosq ц ш w3+3(1-l)sinqcosq ,

(3)

Ч w   3  = ж и lW1+w3tg q-  siny cosq ц ш w2 .

В силу второго уравнения (2) W₁=const, и не зависящие от j уравнения (3) образуют замкнутую систему, в которую W₁ входит как параметр. Система (3) описывает движение оси оси Ox₁ относительно орбитальной системы координат, уравнения (2) описывают движение спутника вокруг этой оси.

Система (3) допускает стационарные решения

q = 0 ,    y = ±  p 2  ,    w2=w3=0 ,

(4)

которые принимаются в качестве невозмущенных движений осесимметричного спутника в режиме закрутки в плоскости орбиты. Для практического использования стационарных решений (4) необходимо, чтобы они были устойчивы. Поскольку уравнения (3) инвариантны относительно замены переменных     y → -y, W₁ → -W₁, w₃ → - w₃, ограничимся описанием свойств устойчивости решения (4), в котором y = p/2.

Уравнения (3) допускают обобщенный интеграл энергии, используя который в качестве функции Ляпунова, можно найти достаточные условия устойчивости решения (4) при y = p/2 [2,3]:

lW1-1 > 0 ,    lW1-(4-3l) > 0 .

(5)

Необходимые условия устойчивости этого решения получаются из анализа соответствующей линеаризованной системы

Ч q   =w2+Dy ,    D Ч y   =w3-q ,

(6)

Ч w   2  =-(lW1-1) w3+3(1-l) q ,    Ч w   3  =(lW1-1) w2 .

Здесь Dy = y-p/2. Характеристическое уравнение системы (6) имеет вид

p4+d1(W1)p2+d2(W1)=0 ,

(7)

d1(W1)=l2W12-2lW1+3l-1 ,   d2(W1)=(lW1-1)(lW1+3l-4) .

Оно биквадратное, поэтому решения системы (6) ограниченны только в том случае, когда все корни уравнения (7) - чисто мнимые и простые, т. е. при d₁ > 0, d₂ > 0, d₁²-4d₂ > 0. Последние неравенства выражают необходимые условия устойчивости решения (4), в котором y = p/2. Эти условия удовлетворяются при выполнении неравенств (5) или неравенств [2]

lW1-1 < 0 ,    lW1+3l-4 < 0 ,   d12-4d2 > 0 .

(8)

Описываемые ниже числовые расчеты проводились при l = 0.7, что примерно соответствует орбитальной станции Мир в последние годы ее полета. Для такого l неравенства (5) выполнены при W₁ > 2.7143, неравенства (8) выполнены на интервалах W₁ < -1.8770 и 1.42462 < W₁ < 1.42857. Последний интервал практического значения не имеет, а движения с W₁ > 2.7143 и W₁ < -1.8770 были реализованы [4,5].

Решению (4) в случае y = p/2 отвечает двухпараметрическое семейство решений полной системы уравнений (2), (3)

j = j0+(W1-1)t ,    j0=const ,    W1=const ,

(9)

y =  p 2  ,    q = w2=w3=0

с параметрами j₀ и W₁. Если W₁ удовлетворяет неравенствам (5) (неравенствам (8)), то решения (9) устойчивы (устойчивы в первом приближении) по переменным W₁, q, y, w₂ и w₃.

Семейство (9) при m ≠ 0 порождает двумерную интегральную поверхность системы (1), которая может быть построена методом Н.Н.Боголюбова - Ю.А.Митропольского [6] в виде формальных рядов по целым неотрицательным степеням m:

j = x+ ∞ е k=1  mkjk(x,h) ,   W1=h+ ∞ е k=1  mkW1,k(x,h) ,

q = ∞ е k=1  mkqk(x,h) ,   y =  p 2 + ∞ е k=1  mkyk(x,h) ,

(10)

w2= ∞ е k=1  mkw2,k(x,h) ,   w3= ∞ е k=1  mkw3,k(x,h) .

Коэффициенты этих рядов суть 2p-периодические функции x. Изменение переменных x и h, т.е. движение по поверхности, определяется уравнениями

Ч x   =h-1+ ∞ е k=1  mkAk(h) ,    Ч h   = ∞ е k=1  mkBk(h) .

(11)

Построение рядов (10), (11) выполняется так. Подставим эти ряды в уравнения (1), разложим левые и правые части получившихся равенств в ряды, аналогичные рядам (10), и приравняем коэффициенты при одинаковых степенях m слева и справа. Получим цепочку линейных неоднородных дифференциальных уравнений относительно коэффициентов искомых рядов. Можно доказать, что при выполнении условия

[k(h-1)]4-d1(h)[k(h-1)]2+d2(h) ≠ 0   (k=0,1,2,...)

(12)

такая цепочка имеет единственное 2p-периодическое по x решение, удовлетворяющее соотношениям

2p у х 0  jk(x,h) dx = 0 ,    2p у х 0  W1,k(x,h) dx = 0   (k=1,2,...) .

(13)

Далее полагаем, что именно это решение использовано в качестве коэффициентов рядов (10), (11).

Условие (12) здесь является главным. Оно означает, что вращение спутника вокруг оси Ox₁ не должно находиться в резонансе с колебаниями этой оси по углам q и y, описываемыми уравнениями (6). Условие (12) выполнено, в частности, при l = 0.7, h < -1.8770 или h > 2.7143. Условия (13) нужны только для единственности решения цепочки, и хотя они использованы ниже, их, в принципе, можно опустить (с потерей единственности) или заменить какими-либо другими условиями.

Будем считать, что решения (10), (11) описывают невозмущенные движения спутника в режиме закрутки вокруг нормали к плоскости орбиты. По этой причине интегральную поверхность (10) назовем интегральной поверхностью ориентированного движения. Наибольший интерес в соотношениях (10), (11) представляет второе уравнение (11), описывающее вековое изменение угловой скорости вращения спутника вокруг оси Ox₁. Чтобы получить представление о таком изменении, достаточно найти первый отличный от тождественного нуля коэффициент B_k(h).

4. Численное исследование интегральной поверхности ориентированного движения. Для сокращения записи воспользуемся векторными обозначениями. Введем вектор¹ z=[q, y, w₂, w₃] и запишем систему (1) в виде

Ч j   =W1+g(z) ,    Ч W   1  =mh(j,z) ,

(1ў)

Ч z   =F0(W1,z)+mF1(j,W1,z) .

Векторная запись системы (3): [(z)\dot]=F₀(W₁,z). Ее стационарное решение (4), в котором y = p/2, обозначим z₀. В последующих выкладках учитываем равенство g(z₀)=-1.

Соотношения (10), (11) представим в виде

j = x+mu(x,h,m) ,   W1=h+mv(x,h,m) ,    z=z0+mf(x,h,m) ,

(10ў)

Ч x   =h-1+mA(h,m) ,    Ч h   =mB(h,m) .

(11ў)

Тот факт, что система (1 ′ ) допускает решения вида (10 ′), (11 ′) выражается тождествами

m(h-1+mA)   ∂ u ∂ x = 1+g(z0+mf)+m(v-A)-m2B   ∂ u ∂ h  ,

m(h-1+mA)   ∂ v ∂ x = m[ h(x+mu,z0+mf)-B ] -m2B   ∂ v ∂ h  ,

m(h-1+mA)   ∂ f ∂ x = F0(h+mv,z0+mf)+mF1(...)-m2B   ∂ f ∂ h  .

Здесь x и h - равноправные независимые переменные, и выписанные соотношения представляют собой уравнения в частных производных относительно u, v и f. В предыдущем разделе был описан способ построения формального решения этих уравнений в виде рядов по степеням m. Построенное решение 2p-периодически зависит от x и удовлетворяет соотношениям (ср. (13))

2p у х 0  u(x,h,m) dx = 0 ,    2p у х 0  v(x,h,m) dx = 0 .

(13ў)

Теперь пойдем другим путем. Отбросим в выписанных уравнениях члены m²B ∂ (...)/ ∂ h. В результате получим дифференциальные уравнения с одной независимой переменной x, т.е. обыкновенные, в которых h является параметром. По-прежнему будем искать решения этих уравнений, зависящие от x и h, причем от x - периодически с периодом 2p, и удовлетворяющее соотношениям (13 ′). Полученные усеченные уравнения и условия, определяющие их искомое решение, запишем используя исходные переменные системы (1) и обозначения a=mA, b=mB. В результате придем к краевой задаче

(h-1+a)  dj dx =W1+g(z) ,   (h-1+a)  dW1 dx =mh(j,z)-b ,

(h-1+a)  dz dx =F0(W1,z)+mF1(j,W1,z) ,

(14)

j(0)+2p = j(2p) ,    W1(0)=W1(2p) ,    z(0)=z(2p) ,

2p у х 0  (j-x) dx = 0 ,    2p у х 0  (W1-h) dx = 0 .

Здесь x - независимая переменная, h - параметр, j, W₁ и z - неизвестные функции, a и b - неизвестные постоянные.

Задачу (14) будем решать методом Пуанкаре. При m = 0 эта задача имеет решение j = x, W₁=h, z=z₀, a=b=0. Линеаризация соотношений (14) в окрестности указанного решения приводит к линейной краевой задаче

(h-1)  dDj dx =DW1-a ,   (h-1)  dDW1 dx =-b ,

(h-1)  dq dx =w2+Dy ,   (h-1)  dDy dx =w3-q ,

(h-1)  dw2 dx =-(lh-1) w3+3(1-l) q ,    (h-1)  dw3 dx =(lh-1) w2 ,

Dj(0)=Dj(2p) ,   DW1(0)=DW1(2p) ,   q(0)=q(2p) ,

Dy(0)=Dy(2p) ,   w2(0)=w2(2p) ,    w3(0)=w3(2p) ,

2p у х 0  Dj dx = 0 ,    2p у х 0  DW1 dx = 0

относительно функций Dj = j-x, DW₁=W₁-h, q, Dy, w₂, w₃ и постоянных a, b. Можно доказать, что при выполнении условия (12) выписанная задача имеет только тривиальное решение. Следовательно, в окрестности точки m = 0 краевая задача (14) имеет единственное решение j = j_*(x,h,m), W₁=W_1*(x,h,m), z=z_*(x,h,m), a=a_*(h,m), b=b_*(h,m), аналитически зависящее от m и удовлетворяющее соотношениям j_*(x,h,0)=x, W_1*(x,h,0)=h, z_*(x,h,0)=z₀, a_*(h,0)=b_*(h,0)=0.

Функции j_*(x,h,m), W_1*(x,h,m) и z_*(x,h,m), продолженные с отрезка 0 ≤ x ≤ 2p на всю действительную ось 2p-периодически по x, удовлетворяют вместе с постоянными a_*(h,m) и b_*(h,m) дифференциальным уравнениям в (14). Эти функции представляются рядами вида (10), коэффициенты которых удовлетворяют соотношениям (13), а постоянные a_*(h,m) и b_*(h,m) представляются рядами, аналогичными рядам, стоящим в правых частях уравнений (11). Если B_m(h) - первый отличный от тождественного нуля коэффициент второго ряда (11), то коэффициенты при первых m степенях m в рядах (10), (11) и в рядах, представляющих рассматриваемое решение задачи (14), совпадают. Более того, b_*(h,m)=m^m B_m(h)+m^m+1B_m+1(h)+O(m^m+2). Это соотношение позволяет использовать краевую задачу (14) для приближенного численного расчета правой части второго уравнения (11).

Численное интегрирование второго уравнения (11) можно выполнить двухцикловым методом [7]. Метод называется так потому, что включает два цикла интегрирования - внешний и внутренний. Внешний цикл состоит в интегрировании уравнения (11) каким-либо численным методом, внутренний заключается в расчете правой части этого уравнения посредством решения краевой задачи (14). В [8] описанный вариант двухциклового метода был применен в очень похожей на рассматриваемую задаче исследования режима одноосной гравитационной ориентации спутника.

Прежде чем перейти к более подробному описанию метода интегрирования упростим задачу (14). Заменим в ней краевое условие и интегральное соотношение для переменной j краевыми условиями j(0)=0, j(2p)=2p. Получившуюся краевую задачу назовем задачей (14 ′). Используя описанную выше схему метода Пуанкаре, можно доказать, что при выполнении условия (12) задача (14′) имеет в окрестности точки m = 0 единственное решение [^(j)](x,h,m), [^(W)]₁(x,h,m), [^(z)](x,h,m), [^(a)](h,m), [^(b)](h,m), аналитически зависящее от m и удовлетворяющее соотношениям [^(j)](x,h,0)=x, [^(W)]₁(x,h,0)=h, [^(z)](x,h,0)=z₀, [^(a)](h,0)=[^(b)](h,0)=0. Решения задач (14) и (14â) связаны между собой соотношениями

^ j   (x,h,m)=j*[x+x0(h,m),h,m] ,    ^ W   1  (x,h,m)=W1*[x+x0(h,m),h,m] ,

^ z   (x,h,m)=z*[x+x0(h,m),h,m] ,    ^ a   (h,m)=a*(h,m) ,    ^ b   (h,m)=b*(h,m) ,

где x₀(h,m) - корень уравнения j_*(x,h,m)=0, удовлетворяющий условию² x₀(h,0)=0 или, что то же самое,

x0(h,m)=  1 2p 2p у х 0  [ ^ j   (x,h,m)-x] dx .

В задаче (14′) перейдем к новой независимой переменной - времени t. Новую задачу будем решать на отрезке 0 ≤ t ≤ T, T=2p(h-1+a)^-1 при h > 1 и на отрезке T ≤ t ≤ 0 при h < 1. Эта задача имеет вид

Ч j   =W1+g(z) ,    Ч W   1  =mh(j,z)-b ,

Ч z   =F0(W1,z)+mF1(j,W1,z) ,

(15)

j(0)=0 ,    j(T)=2p ,    W1(0)=W1(T) ,    z(0)=z(T) ,

T у х 0  (W1-h) dt=0 .

Здесь j, W₁ и z - неизвестные функции, T и b - неизвестные постоянные. В окрестности точки m = 0 задача (15) имеет единственное решение

j = j*(t,h,m) ,    W1=W1*(t,h,m) ,    z=z*(t,h,m) ,

(16)

T=T*(h,m) ,    b=b*(h,m) ,

аналитически зависящее от m и удовлетворяющее соотношениям j_*(t,h,0)=(h-1)t, W_1*(t,h,0)=h, z_*(t,h,0)=z₀, T_*(h,0)=2p(h-1)^-1, b_*(h,0)=0.

Посредством численного решения задачи (15) решение (16) было построено в явном виде и исследовано в функции параметра h. Значения h менялись на равномерной сетке с шагом 0.01. В каждом узле сетки задача (15) решалась методом пристрелки. Неизвестными параметрами были W₁(0), z(0), T и b. Краевые условия в точке t=T и интегральное условие служили уравнениями для их определения. Система таких уравнений решалась методом Ньютона. Вычисление частных производных функций, задающих эти уравнения, по определяемым параметрам (кроме T) сводилось к интегрированию уравнений в вариациях для дифференциальных уравнений в (15). В процессе вычислений узлы сетки проходились последовательно, начиная с узла с максимальным значением |h|. Начальное приближение в этом узле бралось в виде W₁(0)=h, z(0)=z₀, T=2p(h-1)^-1, b=0. Во всех остальных узлах начальным приближением служило решение, найденное в предыдущем узле.

Расчеты проводились при l = 0.7, m = 0.1, g_c=0.01, a_c=-0.15, b_c=0.025, L₁=16 м, L₂=14 м, L₃=12 м, d₁=-0.5 м, d₂=d₃=1 м, e = 3·10^-4 м^-3. В этом случае рассматриваемый спутник может служить грубой моделью станции Мир на орбите с высотой около 330 км. Отвечающее такой высоте значение e получено с использованием статической модели атмосферы [9]. Найденные решения задачи (15) представлены на рис. 1, 2 графиками зависимости от h величин W₁(0), z(0), T и b. Кроме того, здесь представлен график функции d=d(h), характеризующей свойства устойчивости интегральной поверхности ориентированного движения (см. ниже). Найденные решения распадаются на две ветви. Одна ветвь отвечает значениям h > 1, другая - значениям h < 1. Эти ветви не были продолжены в область малых значений |h-1|, поскольку в этой области найденные решения оказались неустойчивы. Примеры решений краевой задачи (15) приведены на рис. 3, 4. Здесь эти решения представлены графиками функций Dj(t)=j(t)-2pT^-1t, q(t), y(t), W₁(t), w₂(t) и w₃(t), построенными на отрезках 0 ≤ t ≤ 2|T|. Как видно из графиков, отклонения оси Ox₁ от оси OX₂ в этих решениях весьма малы, т.е. интегральная поверхность, описываемая найденными решениями, обоснованно названа интегральной поверхностью ориентированного движения.

Исследуем вопрос об устойчивости какого-либо решения j ^° (t), W₁ ^° (t), z°(t), T^°, b^° краевой задачи (15). Под устойчивостью будем понимать устойчивость неподвижной точки [W₁ ^°(0),z ^°(0)] преобразования пространства R⁵[W₁,z], задаваемого решениями дифференциальных уравнений в (15) при b=b ^°. А именно, решение j(t), W₁(t), z(t) с начальными условиями j(0)=0, W₁(0)=W ′ ₁, z(0)=z ′ задает преобразование точки [W ′ ₁,z ′  ] в точку [W"₁,z" ], если в некоторый момент времени T выполнены соотношения W"₁=W₁(T), z"=z(T), j(T)=2p. В окрестности указанной выше неподвижной точки такое преобразование, называемое преобразованием Пуанкаре, определено.

Устойчивость исследовалась следующим образом. Пусть dj(t), dW₁(t) dz(t) - решение системы уравнений в вариациях для решения j^°(t), W₁ ^° (t), z^°(t), T^°, b^° краевой задачи (15). Построим вектор

x= й л  dW1(T°)- Ч W   ° 1  (T°) Ч j   °   (T°) dj(T°),  dz(T°)- Ч z   °   (T°) Ч j   °   (T°) dj(T°)  щ ы  .

Обозначим через x_i (i=1,2,...,5) такие векторы x, вычисленные при dj(0)=0 и [dW₁(0),dz(0)]=[1,0,...,0], [0,1,0,...,0], ... [0,...,0,1]. Составим 5×5 - матрицу X=(x₁,x₂,...,x₅). В случае T ^° > 0 необходимое условие устойчивости неподвижной точки [W₁ ^° (0),z ^° (0)] состоит в том, что модули всех собственных чисел матрицы X должны быть ≤ 1. В случае T ^° < 0 для устойчивости этой неподвижной точки необходимо, чтобы модули всех собственных чисел матрицы X были ≥ 1. Представленная на рис. 1, 2 величина d определяется следующим образом. Пусть r₁ и r₂ - минимум и максимум модулей собственных чисел матрицы X. Тогда d=r₁-1 при T ^° > 0, d=r₂^-1-1 при T ^° < 0. В терминах величины d необходимое условие устойчивости выражается неравенством d ≤ 0.

Как показало вычисление функции d(h), все найденные решения задачи (15) неустойчивы. Однако при h > 3.93 (рис. 1) и h < -2.26 (рис. 2) эта неустойчивость весьма слаба - выполнено неравенство d < 10^-3. На графиках функции d(h) участки слабой неустойчивости выглядят как горизонтальные прямые. Такая слабая неустойчивость встречается во многих задачах вращательного движения спутников, в которых учитывается аэродинамический момент, и обусловлена непотенциальностью этого момента [10]. При h < 3.93 и h > -2.26 неустойчивость решений задачи (15) ярко выражена (см. рис. 1, 2). В обоих семействах решений функция d(h) реализуется действительными собственными числами матрицы X. К тому же, эта неустойчивость возникает при значениях h, лежащих достаточно далеко от границ интервалов устойчивости решения (4) в случае l = 0.7 (см. раздел 3), поэтому ее трудно связать с этими границами.

Вернемся к двухцикловому методу численного интегрирования уравнения [(h)\dot]=b_*(h,m). Заменим его эквивалентным уравнением

dt dh =  1 b*(h,m)  ,

(17)

которое будем интегрировать методом трапеций на сетке по параметру h в процессе продолжения по этому параметру решений краевой задачи (15). Результаты такого интегрирования будем сравнивать с результатами непосредственного интегрирования уравнений (1). Начальные условия сравниваемых решений: для уравнения (17) h(0)=h₀, для системы (1) j(0)=0, W₁(0)=W_1*(0,h₀,m), z(0)=z_*(0,h₀,m).

Сравнение проведем так. Для решения системы (1) определим последовательности чисел

qўN= min t О SN  q(t) ,   q"N= max t О SN  q(t) ,

DyўN= min t О SN  й л y(t)-  p 2 щ ы  ,   Dy"N= max t О SN  й л y(t)-  p 2 щ ы  ,

Wў1N= min t О SN  W1(t) ,   W"1N= max t О SN  W1(t) ,

wў2N= min t О SN  w2(t) ,   w"2N= max t О SN  w2(t) ,

wў3N= min t О SN  w3(t) ,   w"3N= max t О SN  w3(t) ,

SN={t:  2p(N-1) ≤ t ≤ 2pN} ,    N=1,2,3,...

Эти последовательности будем представлять графически ломаными. Например, последовательности {W_1N′} и {W_1N ^′′ } будем изображать ломаными, проходящими соответственно через точки (N,W_1N ′) и (N,W_1N ^′′ )   (N=1,2,...) в плоскости (N,W₁). Две такие ломаные в единых координатных осях дают представление о пределах изменения переменной W₁ на длительном интервале времени. В тех же координатных осях представим график решения h(t) уравнения (17), причем связь между N и t зададим формулой t=2pN, а ось h отождествим с осью W₁. Аналогичным образом последовательности {Dy′_N} и {Dy"_N} будем представлять в единых координатных осях со второй компонентой векторной функции z_*[t,h(t),m]-z₀ и т.д. Такие тройки графиков дают наглядное представление о близости решений системы (1) к интегральной поверхности ориентированного движения.

Сравнение решений системы (1) и уравнения (17) представлено на рис. 5, 6. Помимо указанных выше графиков на этих рисунках приведены еще графики функции d[h(t)], которые характеризуют свойства устойчивости решений, лежащих на интегральной поверхности ориентированного движения. Графики функций h(t) и z_*[t,h(t),m]-z₀ отмечены маркерами. Как видно из рисунков, двухцикловой метод правильно описывает эволюцию решений системы (1), лежащих в окрестности интегральной поверхности ориентированного движения. Точность описания сохраняется до тех пор, пока неустойчивость решений краевой задачи (15) при h = h(t) остается слабой, т.е. до тех пор, пока значения функции d[h(t)] малы.

Нарушение точности после возникновения сильной неустойчивости показано на рис. 7, 8. На этих рисунках представлены те же решения, что и на рис. 5, 6, но на более продолжительном интервале времени. Способ представления несколько изменен. Во-первых, в плоскостях (N,w₂), (N,w₂), (N,q) и (N,Dy) приведены только ломаные, отвечающие последовательностям

^ w   2N  = max t О SN  |w2(t)| ,    ^ w   2N  = max t О SN  |w2(t)| ,

^ q   N  = max t О SN  |q(t)| ,   D ^ y   N  = max t О SN  к к y(t)-  p 2 к к  ,

вычисленным вдоль решения системы (1). Во-вторых, добавлена новая ломаная, характеризующая последовательность

LN= max t О SN  L(t) ,   L = arccos a21 .

Здесь L - угол между осями Ox₁ и OX₂. Этот угол также вычисляется для решения системы (1). В третьих, для удобства анализа график функции d[h(t)] приведен дважды. Как видно из рисунков, резкое возрастание функции d[h(t)] сопровождается резким возрастанием последовательностей {[^(w)]_2N}, {[^(w)]_3N}, {[^(q)]_N}, {D[^(y)]_N} и {[^(L)]_N}, т. е. развитием сильной неустойчивости. При этом резко ухудшается точность описания эволюции последовательностей {W′_1N}, {W"_1N} решением уравнения (17). На рис. 7, 8 маркеры отсутствуют, но сравнение этих рисунков с рис. 5, 6 позволяет легко установить происхождение кривых в плоскости (N,W₁).

На рис. 7 резкое возрастание последовательностей {[^(w)]_2N}, {[^(w)]_3N} и т.п. заметно запаздывает по сравнению с резким подъемом функции d[h(t)]. Это известное явление задержки развития неустойчивости, объясняемое весьма слабыми возмущениями ориентированного движения спутника. После возникновения неустойчивости эти возмущения нарастают экспоненциально, но должно пройти некоторое время прежде, чем они станут заметны в привычном масштабе изучения данной задачи. Для сравнения на рис. 9 приведены такие же графики, полученные для решения системы (1) с возмущенными начальными условиями - к значению w₂(0), задаваемому третьей компонентой вектора z_*(0,h₀,m), было добавлено возмущение 0.1. Теперь запаздывание практически исчезло, а точность описания эволюции последовательностей {W′_1N}, {W"_1N} функцией h(t) практически не изменилась.

Исследуем влияние параметров спутника на эволюцию решений, лежащих на интегральной поверхности ориентированного движения. Для этого повторим описанные выше расчеты решений системы (1) и уравнения (17), изменив знак одного из параметров d₁, d₂, d₃, g_c, a_c и b_c. Результаты четырех таких расчетов, отвечающих сменам знака d₁, d₂, d₃ и a_c, приведены на рис. 10, 11. Каждого расчет представлен графиками функций h(t), d[h(t)] и первых двух компонент векторной функции z_*[t,h(t),m]-z₀, а также ломаными для последовательностей {W′_1N}, {W"_1N}, {q′_N}, {q"_N}, {Dy′_N} и {Dy"_N}. Это сокращенный вариант графиков, приведенных на рис. 5, 6. Как видно из рисунков, изменение знака одного из параметров d₁, d₂, a_c меняет направление эволюции - торможение вращения спутника сменяется его раскруткой (рис. 10 и 11б). Изменение знака одного из параметров d₃, g_c и b_c к изменению направления эволюции не приводит. В случае параметра d₃ этот факт иллюстрируется рис. 11а.

5. Еще один способ отыскания решений, лежащих на интегральной поверхности ориентированного движения. Во время выполнения работы [7] уровень развития вычислительной техники был таким, что только применение двухциклового метода (или какой-либо его альтернативы) позволяло исследовать эволюцию движения спутников на длительных интервалах времени. В настоящее время ситуация существенно иная, и краевая задача (15) вместе с уравнением (17) интересны главным образом тем, что позволяют лучше понять устройство решений системы (1), лежащих на интегральной поверхности ориентированного движения. Кроме того, начальные условия решений краевой задачи (15) дают весьма точное приближение начальных условий ориентированного движения спутника.

Недостаток использованного варианта двухциклового метода состоит в том, что он пригоден только для систем с одной быстрой переменной. В частности, им нельзя исследовать интегральную поверхность ориентированного движения уравнений (1), если в последних учесть эллиптичность орбиты спутника и периодическое изменение плотности атмосферы вдоль орбиты. Существует однако простой метод, который позволяет справиться с трудностью такого рода. Согласно этому методу поиск начальных условий ориентированного движения сводится к минимизации интеграла

J= t у х 0  (z-z0)T(z-z0) dt

(18)

по начальным условиям z(0)=p=[q₀,y₀,...] при j(0)=0 и заданном значении W₁(0). Здесь t - достаточно большое положительное число. Решения системы (1), доставляющие при таких условиях минимум интегралу (18), будем называть минималями.

Минимизация функции J=J(p) сводилась к решению уравнения ∂J(p)/ ∂p=0 методом Гаусса-Ньютона. Рсчет производных ∂J(p)/ ∂p и ∂² J(p)/∂p² выполнялся по формулам вида

∂ J ∂ q0 = t у х 0  (z-z0)T  ∂ z ∂ q0  dt ,     ∂ 2 J ∂q0 ∂y0 ≈ t у х 0  ж и   ∂z ∂q0 ц ш T     ∂z ∂y0  dt .

Производные ∂z/ ∂q₀ и т.п.вычислялись посредством интегрирования уравнений в вариациях для системы (1).

Начальные условия минималей в случае t = 30p приведены на рис. 12. Здесь приведены также графики зависимости от W₁(0) величин d и e, которые характеризуют преобразование Пуанкаре пространства R⁵[W₁,z], задаваемое решениями системы (1) из окрестности найденных минималей. А именно, решение j(t), W₁(t), z(t) с начальными условиями j(0)=0, W₁(0)=W′₁, z(0)=z′ задает преобразование точки [W′₁,z ′ ] в точку [W"₁,z" ], если в некоторый момент времени T выполнены соотношения W"₁=W₁(T), z"=z(T), j(T)=2p. Пусть точка [W₁ ^° ,z ^° ] задает начальные условия минимали и переходит при таком преобразовании в точку [W₁^*,z^*]. Тогда e= √ {(W₁^*-W₁ ^° )²+||z^*-z ^° ||²}. Величина e характеризует смещение точки [W₁ ^° ,z ^° ] под действием введенного преобразования. При e=0 эта точка была бы неподвижной. Величина d, характеризует устойчивость преобразования Пуанкаре в окрестности точки [W₁ ^° ,z ^° ]. Она определяется по линеаризации этого преобразования тем же способом, каким одноименная величина была определена для решений краевой задачи (15). В целом графики на рис. 12 похожи на графики, приведенные на рис. 1. Функция d[W₁(0)]) ведет себя очень похоже на функцию d(h).

Пример минимали приведен на рис. 13. На рис. 13а изображены графики функций W₁(t), q(t) и y(t) (0 ≤ t ≤ 2p), из которых видно, что на минимали отклонения оси Ox₁ от оси OX₂ весьма малы, т.е. минималь лежит в окрестности интегральной поверхности ориентированного движения. На рис. 13б для той же минимали показаны ломаные, описывающие последовательности {W′_1N}, {W"_1N}, {[^(q)]_N} и {D[^(y)]_N}. Как видно из рисунка, при достижении величинами W′_1N и W"_1N окрестности точки резкого возрастания функции d[W₁(0)] (примерно при W₁=3.92) происходит быстрое разрушение ориентированного движения (ср. рис. 7).

Данная работа выполнена при поддержке РФФИ (проект 02-01-00323).

Литература

[1]

Бозюков А.Ю., Сазонов В.В. Об одном способе гравитационной ориентации вращающегося спутника. Препринт Института прикладной математики РАН, 2004, N 24.

[2]

Белецкий В.В. Движение искусственного спутника относительно центра масс. М., Наука, 1965.

[3]

Черноусько Ф.Л. Об устойчивости регулярной прецессии спутника. Прикладная математика и механика, 1964, т. 28, вып. 1, с. 155-157.

[4]

Бабкин Е.В., Беляев М.Ю., Ефимов Н.И., Сазонов В.В., Стажков В.М. Неуправляемое вращательное движение орбитальной станции Мир. Космические исследования, 2001, т. 39, N 1, с. 27-42.

[5]

Бабкин Е.В., Беляев М.Ю., Ефимов Н.И., Сазонов В.В., Стажков В.М. Неуправляемое вращательное движение орбитальной станции Мир в последние месяцы ее полета. Космические исследования, 2003, т. 41, N 3, с. 285-294.

[6]

Боголюбов Н.Н., Митропольский Ю.А. Асимптотические методы в теории нелинейных колебаний. М., Физматгиз, 1963.

[7]

Охоцимский Д.Е., Энеев Т.М., Таратынова Г.П. Определение времени существования искусственного спутника Земли и исследование вековых возмущений его орбиты. Успехи физических наук, 1957, т. 63, N 1а, с. 33-50.

[8]

Сазонов В.В., Петров А.Л. Эволюция режима гравитационной ориентации вращающегося спутника под действием непотенциального аэродинамического момента. Космические исследования, 1987, т. 25, N 4, с. 508-522.

[9]

Навигационное обеспечение полета орбитального комплекса Салют-6 - Союз - Прогресс. М., Наука, 1985.

[10]

Сарычев В.А. Вопросы ориентации искусственных спутников. Итоги науки и техники. Исследование космического пространства, т. 11. М., ВИНИТИ, 1978.

Footnotes:

¹  Запись [x₁,...,x_n] означает вектор-столбец с соответствующими компонентами.

² Поскольку j_*(0,h,0)=0 и dj_*(0,h,0)/dx = 1, по теореме о неявной функции в окрестности точки m = 0 такой корень существует и единствен.

Рисунки