CC BY
17
УДК 531
Г. В. Калачев, А. Б. Нуралиева, А. В. Чернов
Московский государственный университет им. М.В.Ломоносова Институт прикладной математики им. М.В. Келдыша РАН Московский физико-технический институт (государственный университет)
Малые колебания троса космического лифта
Рассматриваются малые колебания троса космического лифта относительно устойчивого вертикального положения равновесия. Трос гибкий, нерастяжимый, переменной линейной плотности. Нижний конец закреплен на Земле в районе экватора, а верхний уходит за геостационар. Система держится в равновесии гравитационными и центробежными силами. Для линеаризованной по отношению к поперечным отклонениям системы решается задача Штурма-Лиувилля и найдены собственные частоты и собственные формы линейных колебаний. Приведены асимптотические (для высоких мод) формулы и описан алгоритм их численного нахождения. Найденные формы использовались для вычисления движений в численной модели нелинейных колебаний троса. Периодичность полученных движений показывает соответствие численного и аналитического подходов. Проанализировано распределение натяжения вдоль троса при малых колебаниях.
Ключевые слова: тросовые системы, космический лифт, малые колебания, задача Штурма^ Лиу вилля.
1. Введение
Рассматриваются малые колебания в плоскости экватора гибкого нерастяжимого троса, закрепленного на экваторе, в гравитационно-центробежном поле. Такой трос является основным элементом космического лифта [1], и его динамика определяет характеристики всего сооружения. Исследование выполнено на основе непрерывной математической модели, изложенной в [2, 3]. Рассматривается динамика свободного троса (без учета движения кабины). Разработана методика вычисления собственных частот и форм линейных колебаний троса, найдены асимптотические по номеру моды формулы для этих частот и мод. Полученные результаты проверены и до некоторой степени развиты численными расчетами для нелинейной модели. Найденные разными методами результаты достаточно хорошо совпадают.
2. Уравнения движения
Рассмотрим движение гибкого троса в плоскости экватора, один конец которого закреплен на планете, а другой свободен, но на нем закреплена концевая масса М. Введем систему координат, вращающуюся вместе с планетой с угловой скоростью ш,ис центром, совпадающим с центром планеты О. Ось Ох проходит через точку прикрепления троса, ось О г направлена вдоль оси вращения, а ось Оу дополняет их до правой тройки векторов. Положение точки на тросе задается линейной координатой в, отсчитываемой вдоль троса от точки крепления троса к планете. Для компактности записи будем обозначать точками производные от переменных по £ и штрихами - производные по в:
(1)
В этом случае уравнения динамики такого троса записываются в виде (1) с краевыми условиями (2). Здесь Ь — длина троса, р = р(,в) — линейная плотность троса в точке в, ^ — гравитационный параметр, Т = Т(в) — сила натяжения троса в точке в, £ — время, г — расстояние текущей точки троса от центра. К уравнениям (1) и (2) необходимо добавить условие нерастяжимости троса (3):
(X )2 + (у')2 = 1. (3)
Определение. Рабочим состоянием троса назовем его вертикальную равновесную конфигурацию (4):
х = К + в; у = 0. (4)
Будем считать, что трос состоит из двух компонент: несущей и дополнительной. Площадь поперечного сечения несущей компоненты обозначим 5 = Б (в), ее объемная плотность ру, линейная плотность дополнительной компоненты постоянна по всей длине троса и равна ра. Тогда будем иметь
р = Бру + ра. (5)
Для равнонапряженного троса Т(в) = Б(в)а, где константа а — допустимое по условиям прочности напряжение троса. В этом случае, известно из [2, 3], условие (6)
Б = (5о + ^) ехр ру (и (КЕ + з) -и (КЕ)) - —; и (г) = - ^ - \ш2г2. (6)
ру а ру г 2
II (г) — гравитационно-центробежный потенциал, Бо — площадь поперечного сечения троса в нижнем сечении. Соотношения (4) - (6) полностью определяют равновесное вертикальное (рабочее) состояние троса.
Для изучения динамики троса в окрестности этого состояния выполним линеаризацию задачи (1) - (3). Введем вариации х\(,з, Ь), у\(з, Ь), Т\(в, Ь) переменных х, у и Т, малые вместе со своими первыми и вторыми производными, так что
х = хо +Xi, у = Уо + УъТ = То + Ti,
где xo(s, t) = Re + s, yo(s, t) = 0, To(s, t) = S(s)a.
В дальнейшем будем пренебрегать членами выше первого порядка относительно Xi, yi, Ti. № условия нерастяжимости (3) следует x^s, t) = 0, г = х0 = Re + s, а изменение сечения вдоль троса дается формулой (6). В результате линеаризации и замены yi = (Re + s )yi получаем краевую задачу (7) на вариации yi и Ti(s, t):
T i
0 = MRe + s№ + -j, ii = ^(j+REy-2 (t"(s + Re )W,
2u,(Re + L)ii(L) - Ш=0. ш + T§)vi =0.
3. Разделение переменных
Из (7) выделяется краевая задача на вариацию yi(s, t). В дальнейшем индекс «1» при обозначении вариаций будем опускать.
Воспользуемся методом разделения переменных и будем искать решения в виде ряда y(s, t) = (s)@k(t), где Yk(s)@k(t) — частные решения системы.
Ykвk = (To(s + Re)2Yl)'@k ^ в = (Tf lRE)ikУ = -\k = const.
p(s + Re )2 вк p(s + Re )2Yk
Далее покажем, что Хк > 0, поэтому Ок = Ак соъ(у/\ЦЬ+фк) — решение вспомогательной задачи.
(Т0 (в + КЕ )2Ук)' + Хк р(з + КЕ )2 Гк = 0; ¥к(0) = 0; Т0(Ь)¥к(Ь) - ХкМУк(Ь) =0.
Для Ук находим краевую задачу (8), где Хк — спектральный параметр. Это задача Штурма-Лиувилля, в решении которой и состоит основная сложность. Здесь Хк входит в краевое условие, что не позволяет использовать в данном случае стандартные результаты теории Штурма-Лиувилля.
В дальнейших выкладках индексы к опущены.
Утверждение 1. Собственные значения вещественны и положительны, а собственные функции ортогональны.
Доказательство. Интерпретируем задачу (8) как поиск спектра дифференциального оператора Б:
ш = (То (8 + Пе )2У')' р(з + Пе )2 .
Данный оператор действует на множестве В функций из С2[0,Ь], удовлетворяющих краевым условиям из системы (8). В — линейное пространство в силу того, что оба ограничения линейны. Введем на множестве В скалярное произведение:
ь ь
(и,ь)= ! и(8)у(8)йV = ! и^Ь^р^ + КЕ)2йв + Ми(Ь)ь(Ь)(Ь + КЕ)2. о о
Оператор И является симметричным оператором относительно введенного скалярного произведения. Действительно, с учетом формул (8) имеем
МОи(Ь)у(Ь)(Ь + Ке )2 = -Т0 (Ь + Ке )2и' (Ь)у(Ь).
Поэтому
ь ь
(Ии, V) = !(Т0(,в + Ке)2и')'ьйв - Т0(Ь + КЕ)2и'(Ь)ь(Ь) = - ^ Т0(,в + КЕ)2и'^й8 = (и, Бь). о о
Отсюда следует вещественность собственных значений оператора и их ортогональность. Из полученного равенства легко видеть, что оператор И является отрицательно определенным, а равенство (Ии, и) = 0 достигается только при и = 0. Утверждение доказано.
4. Асимптотика собственных значений и собственных функций
Для упрощения дальнейших действий введем следующие обозначения:
, 2 То (в + КЕ )2 ( 1 V' , , р Х = Х; Р1(8) =-2-(го{в + КЕ)2) ; Р2(8) = Т0;
ОД = у/рф); ВД= [ П(т)йт.
Jо
Заменой У = + КЕ)^Т0) 1 У система (8) сводится к краевой задаче (9):
У" + Р1У + Х2Р2У = 0;
П0) = * У'(ь) - ущ (Ж + + ££) = (9)
2То(Ь) Ь + Ке То(Ь))
Для этой задачи выполняются следующие условия:
1) функции Р\ £ С\ Р2 £ С2.
2) Р2($) > 0, поэтому корни Ф^ Ф2 характеристического уравнения Ф2 + '2(5) = 0 различны между собой: Ф^) = г°(,в), Ф2(в) = — )■
В работе [5] показано (см. также [6], [7] и [8]), что фундаментальная система решений для задачи (9) представима в виде (10) при выполнении условий 1) - 2).
ьк(*,Х)=ехр ((-фхЩв)) (Е^, к = 0,1. (Ю)
Здесь Ик^(з) — некоторые непрерывно дифференцируемые функции. Подставляя эти решения в дифференциальное уравнение в системе (9) и заменяя и^^ = и^1, находим
= 0; ик>] = (-1)4 1 (и^—О-2)" + Ц. Для упрощения дальнейших выкладок введем новое обозначение
г8 '1 1 ( 1 ,0 ^2( 1°-') ■ 2
ик,
^=К1 (/2 )"+1
ик,0 = Ск,о; ик,з = Ск,у + (-1)кг Аи,^-1.
Здесь си,] — это константы интегрирования. Отсюда находим следующее выражение для функций
^ = J2((-i)ki Уа1сы-1. 1=0
Подставляя выражение, полученное для ик,}, в (10), находим функции, определяющие фундаментальное решение краевой задачи (9):
(œ j
T,T,((-1)ki ?А1сЫ-1 Х-
3=0 1=0
Переставляя порядок суммирования в получившемся ряде, находим пару независимых решений дифференциального уравнения (9).
exp ((-1)кгХВД) ~
Vk(s, Х) =-^^-¿ Y,((-1)H) А Х- • (И)
1=0
Решения краевой задачи (9) будем искать в виде линейной комбинации функций
vo(s, х) = vo(s,Х)+ vi(s, Х), vi(s, Х) = M*,Х) - П(s,Х)-
А
Ш) U=o=0; + А) 'А = А{^ + А) 1 •
Используя Q(0) = 0, неодим, что v1(0,x) = 0, a г>0(0,х) = 2Q.-1/2(0) = 0. Следовательно, собственные функции краевой задачи (9) имеют вид
Ф,х) = У-1/2 (вш (ХУ)( 1+ 1 (1)+сов (ХУ)(1 + Ц) .
Возвращаясь к представлению обратного оператора в виде ряда и подставляя полученный результат в краевое условие при 8 = Ь, находим
= 0;
в=Ь
№) - ^Х2) + сов [ХУ) +(п - Щ Х) =( ) у2У + <18 2То(8 + Ке)2 + То )
V( & ( й А2+ + аАШ1 йА*+2 (Т0(8 + ЕЕ)2У А23+1 + МА23+Л
=(-1) ^- йА - 2То(8 + Ке)2 А + %А )
,,кр = ( ^423+1 + (1А23+1 ПА.23+2 _ СТо(в + Ке)2) А2]+1 + ^А2]+3\ у-2]-1
^2У + йз 2То(в+ " Л2 + ^ ,А
те
Собственные значения ищем в виде ряда х(п) = Хо(п) + Т1, Хгп-ггДе Хг ~ некоторые
1=1
константы. Тогда
_1_ _ у/р(Ь)То(Ь) -МР(Ь)_
Хо(п) = жпУ(Ь) ; Х1 =
жМ
р <* = \ I ((?)1 ((?)1)' + ^ (^)"(?)21-
Подставляя полученную асимптотику в формулу для собственной функции, находим
Уп(8,х) =й-1/2 +
1 I й(8) \ ( у/Щ
+--сов жп ~ | | „ ,
жп
Ущ) ( /Г | ь й(*) - Р(Ь)у(8) + у(Ь)Р(в)^ .
Возвращаясь к \п и Уп(в), получаем следующие выражения для собственных значений и функций:
Ап = ж2п2 ( + 2 (- Р (^ (¡\ild8] ; (12)
То \ М
оо
„ , \ (рТо)-1/Ч . ( й(8) \ п ' 8 + Ке \ ^ УУЩ
+уо) сов (жп м) (мщу)1» ш - р +р
жп ^ УЩ)\ м ущ к ;У(Ь)
(13)
Выписанные формулы дают асимптотические приближенные для собственных значений п п
описанный ниже.
5. Алгоритм вычисления собственных значений и собственных функций
Следуя методике, изложенной в [9], сделаем замену
То(8 + Ке )2У = й(8)соъф);
(14)
ХМ (Ь + НЕ )2У (з) = а^втф).
Поскольку Т(в) > 0 и У (з) и У' (в) одновременно не равны нулю, то эта замена взаимно однозначна и 1(в) = 0. Подставляя в (8) соотношения (14), получаем
$ и р(НЕ + в)2 1 = * ** — М(Не + Ь)2 8ш *(0) = 0; 1(Ь) сов *(Ь) — 1(Ь) 8т *(Ь) = 0.
Без потери общности можно считать, что *(0) = 0. Из второго краевого условия получаем, что tg*(Ь) = 1. Дифференцированием второго равенства (14) вместе с полученным выше соотношением получаем независимую краевую задачу (15) для
' ХМ(НЕ + Ь)2 2 , Р(Не + 8)2 . 2 * = Т(КЕ + 8)2 СОВ * + М(Не + Ь)2 (15)
*(0) = 0; %Е*(Ь) = 1.
Хк
которых краевая задача (15) имеет решение.
Лемма 1. Пусть у(х, а) - решение задачи Коши у'х = д(х, у, а), у(0) = 0, причем д(х, у, а) строго возрастает по а при х € [0,Ь),у > 0. Тогда Vх0 € [0, Ь) у(х0, а) возрастает по а. Если х0 < Ь, то возрастание строгое.
Доказательство. Допустим а < Ъ и у(х0,а) > у(х0, Ь). Пусть I := тах{х € [0; х0] : у(х,а)<у(х, Ь)}. Тогда
д(1, у(1, а), а) — д(1, у(1, Ь), Ь) = у'(1, а) — у'(1, Ь) =
= Иш У(1 + Ъ,а) — у(1, а) — у(1 + Ъ, Ь) + у(1, Ь) > 0 ь^+о Ъ ~ '
так как сумма нечетных слагаемых в числителе по предположению положительна, а чет-аргументу.
Следствие. Если *\(8) - решение задачи (15) без второго краевого условия, то *\(Ь) возрастает по Х. (Вытекает из доказанной леммы, т.к. правая часть уравнения для * из
Х
Хк
ровать уравнение (15) и, пользуясь монотонностью *(Ь) по Х, находить последовательные значения Хк, при которых *(Ь) принимает значения ж/4 + жк(к = 0,1, 2,...), на которых удовлетворяется второе краевое условие задачи (15).
После вычисления собственного значения соответствующая ему собственная функция находится совместным интегрированием уравнений в переменных 1 и
' ХМ(НЕ + Ь)2 2 , Р(НЕ + 8)2 . 2 * = Т{Не+ 8)2 СОВ * + М(Не + Ь)2
..ХМ (Не + Ь)2 Р(Не + в)2 . 1 = 1( Т(НЕ + 8)2 — М(НЕ + Ь)2СО8*,
и затем от этих переменных нужно вернуться к исходным по формулам (14).
32
Аэрокосмические исследования
ТРУДЫ МФТИ. — 2013. — Том 5, № 4
6. Численное изучение малых колебаний космического лифта
Результаты описанной выше работы использовались для расчетов в программах, численно моделирующих движение лифта в нелинейной модели. Во-первых, это служит обоюдной проверкой численных и аналитически полученных результатов. Во-вторых, численные результаты дополняют аналитические, например, исследованием влияния движения троса на напряжение в нем. Численные результаты также позволяют определить границы применимости линейной модели.
При вычислениях собственные формы линеаризованной задачи задавались как начальные формы троса. Результаты расчетов показаны на рис. 1 на примере первых 8 мод колебаний. Для наглядности на каждой фигуре рисунка изображены положения троса через некоторые интервалы времени для одного движения. Рисунок растянут по горизонтальной оси, т.к. отклонение составляет порядка 0.001 от длины. (На дальнейших рисунках изображения троса тоже растянуты по горизонтали.) При достаточно малой амплитуде колебания получаются периодическими. Узлы остаются на месте на протяжении десятков периодов. Это показывает соответствие численных и аналитических расчетов.
Рис. 1. Первые собственные моды колебаний
Здесь можно видеть собственные формы колебаний, п-я собственная форма имеет п пучностей п п +1 узел и близка к синусоиде, у которой на длине троса располагаются п/2 периодов (хотя синусоидой, конечно, не является). В случае четной собственной моды п п/2
Фурье, в случае нечетной собственной моды п две ((п - 1)/2, (п + 1)/2).
Поведение напряжения в линейных колебаниях. Вариация 5Т(в, I), то есть поправка 1-го порядка к силе натяжения троса, раньше не вычислялась, так как с точностью д0 линейных членов не входит в уравнения для вариаций отклонений. Из формул для из (1), (2) следуют соотношения
о гт
-+ 2^Ерёу = 0, 8Т(Ь, г) - 2мЕМ6У(Ь, I) = 0,
дв
которые представляют собой дифференциальное уравнение для 5Т(в, I) и краевое условие. Отсюда 5Т(в, I) находится квадратурой
5Т( 8, г) = 5Т(ь, г) + 2ше1 Р(£)Н£, № = 2^Е И , № + М6У(Ь, г)) . (16)
б \ в
Подставляем вариацию скорости ёу(в, I) как производную но времени 5у:
5Т = 2ш^Ак вк №к (з), к
(ь \
где Ак — амплитуда колебаний к-й моды, 2к($) = / Ук(0р(0$^ + Ук(Ь)М . Для одно-модовых колебаний
6Т = 2шЕАк 9к №к (в). (18)
Если при Ь = 0 трос неподвижен и форма его соответствует к-ш моде колебаний, то
5у(в, г) = АкУк (в) сов(\/Х^, 5 у(в, г) = —АкУк (в)\/Х ^(^/Хг). (19)
Формула (17) принимает вид
6Т(8, г) = —2шЕАклД~к вш(^Хс*)Ек(в). (20)
Напомним, что Т(в, ¿) = То (в) + 5Т(в, ¿), где То — сила натяжения в вертикальном неподвижном тросе. Из (17) видно, что сила натяжения изменяется периодически по времени и достигает экстремумов со сдвигом по фазе на четверть периода по сравнению с поперечным отклонением §у. Т совпадает с То, когда фаз а = жп, то есть когда отклонения троса принимают амплитудные значения п € 2 или когда 2к (я) = 0. Заметим, что вариация силы натяжения в точке синфазна или противофазна скорости 5у(,в, Ь) в (16), в зависимости от знака выражения 2к(в). Вообще, знак этого выражения может меняться
Ак
определяется из начального отклонения.
Для вычисления напряжения а(в, ¿) разделим натяжение в точке на площадь сечения:
2. (8)
а(з, г) = <го + $<г = ао + 2шЕвк(Ь)Ак^^-. (21)
Или, в другой форме:
ь
Мб Ук (Ь, $ + / р(£)5 Ук (£, Ш
а(в, Ь) = <о + = <о + 2ше--. (22)
где <о """"" номинальное напряжение, то есть напряжение в вертикальном неподвижном тросе. Для равнонапряженного троса <о > 0 — константа.
Трос может порваться, если напряжение больше разрывного напряжения <ьг, которое является характеристикой материала. Но при малых колебаниях напряжение не приближается к разрывному, поскольку номинальное напряжение <о выбирается с запасом прочности. Напряжение в рассматриваемых примерах меняется мало — на 0.1-0.01% от начальной величины, но ведет себя достаточно интересно.
Например, на рис. 2 представлены парами профиль троса и напряжение в один и тот же момент времени.
Начинается движение с гладкой кривой напряжения, в которой напряжение монотонно убывает (2а), потом быстро образует волнообразную кривую, у которой столько экстремумов, сколько узлов у собственной формы троса (26). Такая примерно форма кривой напряжения держится до половины периода. Амплитуда колебаний напряжения увеличивается по мере движения от крайнего положения до среднего, т.е. при фазе от 2пТ до (2п + 1)Т (2в — четверть периода), и уменьшается при движении от среднего положения к крайнему, т.е. при фазе от (2п + 1) 2 до (2п + 2) п € 2. Потом на половине периода кривая становится монотонной (2г), затем перегибается в обратную сторону, т.е. образует столько же экстремумов (по числу узлов исходной формы), но минимумы и максимумы поменялись местами. И держится в таком виде еще полпериода.
34
Аэрокосмические исследования
ТРУДЫ МФТИ. — 2013. — Том 5, 1
Рис. 2. Пары профиль троса и напряжение в один момент времени
На следующих рисунках для наглядности изображены положения троса через некоторые интервалы времени, как на рис. 1. Чтобы отразить на том же рисунке изменение напряжения, участки с разным напряжением показаны разным оттенком. Кривые, соответствующие разным моментам времени, могут накладываться, более поздние перекрывать более ранние. Наблюдаются несколько интересных эффектов.
На рис. 3 показаны колебания для 14-й моды. Участки троса с напряжением больше номинального сто закрашивались более светлым цветом, меньше - более темным. Напряжение больше при более прямом тросе и меньше при амплитудных значениях. Этот эффект наблюдается на более высоких модах, где напряжение сильнее зависит от фазы движения, от того, находится ли трос в самом изогнутом положении (время кратно половине периода) или распрямлен (время составляет нечетное число четвертей периода). Скорее всего, большее напряжение при прямом тросе объясняется большей скоростью (22), которая в свою очередь математически объясняется большим множителем (19), а физически - тем, что колебания высоких мод более высокочастотные.
Иногда на протяжении всего движения напряжение на одном участке троса почти не меняется. На рис. 4а изображена первая половина периода, на рис. 46 вторая, дальше все продолжается периодически. Рисунок выглядит сплошным, потому что положения троса отрисовывались часто. Окраска следующая: номинальное напряжение обозначено серым,
Рис. 3. 14-я мода
Рис. 4. 1-я и 2-я половины периода. Стрелками показаны участки с почти постоянным по времени напряжением
чем напряжение больше, тем цвет светлее, чем меньше, тем темнее. Показаны стрелками два узких участка троса (две серые полосы), на которых напряжение практически сохраняется. Остальные участки троса кажутся черными или светло-серыми потому, что вариация напряжения на них, хотя в абсолютном значении мала, значительно больше вариации напряжения на участках постоянства напряжения, которая почти равна нулю. Видно, что 1-я и 2-я половины периода «инвертированы» по напряжению. Этот эффект наблюдается на низких модах.
Существование же участков, на которых напряжение почти не меняется, связано с тем, что для этой формы троса находятся координаты s , при которых L
^■k(s) = fYk+ (L)M = 0 то есть Sa обращается в 0 (см. формулу (21)).
S
Границы линейности. Вычисления показывают, что периодичность сохраняется при отклонении порядка десятков километров (для более низких мод больше, для более высоких - меньше), т.е. для отклонения порядка 0.001 от длины. Данная работа продолжает и развивает результаты, изложенные в |2|, [3|. Она позволяет связать аналитические результаты для линейных колебаний с численными результатами, изложенными в [2|. Возможно дальнейшее развитие, и полученные данные (в частности, собственные формы) могут дать новые инструменты для изучения нелинейных колебаний.
Авторы выражают благодарность Ю.А. С адову за внимание к работе и полезные сове-
Литература
1. Садов Ю.А., Нуралиева А.Б. О концепции нагруженного секционированного космического лифта: препринт / Институт прикладной математики им. М.В. Келдыша РАН. - 2011. - У* 39. - 24 с.
2. Садов Ю.А., Нуралиева А.Б. Нелинейные поперечные колебания троса космического лифта // Математическое моделирование. — 2011. — Т. 23, .Т2 12. — С. 3-19.
3. Чернов A.B. О динамике космического лифта в окрестности рабочего состояния // Информационные технологии: модели и методы. — 2010. — М.: МФТИ. — С. 3-12.
4. Поляков Г.Г. Собрание сочинений. Т. 1. Привязные спутники, космические лифты и кольца. — Астрахань, 1967-1974.
5. Тамарьшн Я. Д. О некоторых общих задачах теории обыкновенных дифференциальных уравнений и о разложении произвольной функции в ряды. — Петроград, 1917.
6. Наймарк М.А. Линейные дифференциальные операторы. — М.: Наука, 1969.
7. Садовничий В.А., Султанаев Я. Т., Ахтямов А.М. Обратные задачи Штурма-Лиувилля с нераспадающимися краевыми условиями. — М.: Изд-во Моск. ун-та, 2009.
8. Tamarkin J.D. The notion of the Green's function in the theory of integro-differential equations 11 Trans. Amer. Math. Soc. - 1927. - V. 29. - P. 755-800.
9. Петров А.Л., Садов С.Ю. Исследование близких к маятниковым движений тяжелой нити на круговой орбите: препринт / 1111 \! им. М.В. Келдыша АН СССР. — 1989. — № 11.
Поступила в редакцию 12.07.2012
