Прикладная математика и механика, 2020, T. 84, № 5, стр. 543-553

О линейных поперечных колебаниях троса космического лифта

Ю. А. Садов 1, А. Б. Нуралиева 1*

1 Институт прикладной математики им. М.В. Келдыша РАН
Москва, Россия

* E-mail: annanuralieva@yandex.ru

Поступила в редакцию 10.10.2019
После доработки 10.07.2020
Принята к публикации 18.07.2020

Полный текст (PDF)

Аннотация

Исследуются малые поперечные колебания сверхдлинного гибкого нерастяжимого троса космического лифта относительно равновесного вертикального положения. Показано, что в линейном приближении колебания в плоскости экватора и в меридиональной плоскости разделяются и имеют одинаковые наборы собственных функций, а их собственные частоты связаны простым соотношением. Описана постановка задачи Штурма–Лиувилля для определения собственных частот и собственных мод колебаний троса. Обсуждается ее преобразование с целью упрощения расчетов.

Ключевые слова: космический лифт, гибкий нерастяжимый трос, равнонапряженный трос, малые колебания, задача Штурма–Лиувилля, подстановка Прюфера

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

Первые данные о линейных колебаниях КЛ были получены еще в 1975 г. Пирсоном [1] с использованием простых динамических моделей. Так, на основе аналогии с центробежным математическим маятником была выполнена оценка периода маятниковых (с прямым тросом) колебаний, который оказался равным нескольким суткам, и в результате сделан вывод об отсутствии резонанса таких колебаний с основными лунно-солнечными возмущениями.

Более подробные вычисления были проведены на основе динамической модели Белецкого [2], в которой лифт моделируется точечной массой, связанной невесомым тросом с поверхностью Земли. Были получены [3] линеаризованные уравнения таких колебаний и выведена простая формула, связывающая частоты маятниковых колебаний в двух плоскостях – меридиональной (${{T}^{{{\text{mer}}}}}$) и экваториальной (${{T}^{{{\text{eq}}}}}$).

Из указанной формулы сразу следует [3], что в одномассовой модели ${{T}^{{{\text{mer}}}}}$ всегда меньше суток; вычисления показывают, что ${{T}^{{{\text{eq}}}}}$ больше суток. Оценка периодов изгибных колебаний не может быть сделана в рамках простой одномассовой модели КЛ, но такая же формула (см. (4.9)) ниже – справедлива и для изучаемой здесь модели троса с непрерывно распределенной массой; более того, она распространяется и на другие моды колебаний.

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

Уравнения колебаний троса КЛ с распределенной массой выведены в [4]; здесь вывод кратко воспроизводится. Работа [5] близка к данной статье, пересекается с ней в ряде результатов и тоже использует описанную [6] методологию анализа малых колебаний троса с грузом на конце.

С использованием дискретной, многомассовой модели троса были получены [7] численные данные о периодах и некоторых других характеристиках изгибных мод низкого порядка (учитывались также линейные упругие колебания на основе предполагаемых, малоизвестных пока данных об упругих свойствах материала троса). Обобщение этой модели, допускающее внеэкваториальную стартовую позицию троса КЛ, рассмотрено в статье [8]. Там же содержится достаточно полная библиография работ по динамике КЛ в рецензируемых изданиях.

Настоящая работа касается углубленной проработки важной, но частной задачи, связанной с тематикой космического лифта. Более широкая перспектива намечена в работе авторов [9]. Представляют интерес публикации американской организации ISEC (The International Space Elevator Consortium), где обсуждается комплекс научных и инженерных проблем, связанных с КЛ.

В ИПМ им. М.В. Келдыша РАН с конца 1960-х гг. ведутся работы по космическим тросовым системам (КТС). Изучение вопросов, связанных с космическим лифтом, было инициировано В.В. Белецким как ответвление от этого направления. Несколько страниц монографии [10] посвящено космическому лифту. Несмотря на схожесть основных уравнений динамики гибкого троса в моделях КТС и КЛ, количественные и качественные характеристики этих систем сильно различаются.

2. Уравнения движения троса. Рассматривается движение абсолютно гибкой, весомой, нерастяжимой нити, неподвижно закрепленной одним концом на экваторе Земли. На втором конце закреплено тяжелое твердое тело (материальная точка). Вся система движется в поле тяжести и инерционной центробежной силы. Для описания движения нити в непрерывной модели используем связанную с вращающейся Землей систему координат, центр которой $O$ совмещен с центром Земли, ось $x$ направлена из точки $O$ к точке закрепления нити на экваторе, ось $z$ параллельна оси вращения Земли и направлена на север, ось $y$ параллельна линии экватора в точке закрепления и направлена на восток.

Положение точки нити, находящейся на расстоянии $s$ вдоль нити от точки закрепления, задается вектором ${\mathbf{r}} = (x(s,t),y(s,t),y(s,t))$, где $t$ – время. Пусть $\rho (s)$ – локальная линейная плотность нити в точке $s$ и $P(s,t)$ – сила натяжения, действующая вдоль нити на нижний участок (с меньшими значениями $s$) троса со стороны верхнего участка и направленная в точке $s$ в положительном направлении по $s$. Тогда уравнение движения для текущей точки троса с учетом силы тяжести, центробежной и кориолисовой сил во вращающейся вместе с Землей системе координат и силы натяжения троса записывается так:

(2.1)
${\mathbf{\ddot {r}}} = - 2\left[ {\omega ,{\mathbf{\dot {r}}}} \right] - \left[ {\omega ,\left[ {\omega ,{\mathbf{r}}} \right]} \right] - \frac{\mu }{{{{r}^{3}}}}{\mathbf{r}} + \frac{1}{\rho }\left( {P{\mathbf{r}}{\kern 1pt} '} \right),$
где $\mu $ – гравитационный параметр, $\omega $ – вектор сидерической угловой скорости Земли. Точкой обозначена производная по $t$, штрихом – по $s$.

К этому уравнению необходимо добавить условие для определения натяжения $P(s,t)$ (в данном случае – условие нерастяжимости нити)

(2.2)
${\text{|}}{\mathbf{r}}{\kern 1pt} '{\text{|}} = 1$

К системе уравнений с частными производными (2.1), (2.2) требуется добавить начальные и краевые условия. Начальными данными служат значения ${\mathbf{r}}(s,0)$ и ${\mathbf{\dot {r}}}(s,0)$.

Обозначим через $L$ длину троса, ${\mathbf{R}}$ – радиус-вектор точки его закрепления, $R = {\text{|}}{\mathbf{R}}{\text{|}}$ – экваториальный радиус Земли.

Краевые условия на нижнем конце троса выражают стационарность его закрепления в нижней точке. Условия на верхнем конце, где к тросу прикреплена массивная точка (противовес), выражают согласованность движения этой точки с движением верхнего конца троса. Уравнение движения точки на верхнем конце имеет вид

(2.3)
${{{\mathbf{\ddot {r}}}}_{L}} = - 2\left[ {\vec {\omega },{{{{\mathbf{\dot {r}}}}}_{L}}} \right] - \left[ {\vec {\omega },\left[ {\vec {\omega },{{{\mathbf{r}}}_{L}}} \right]} \right] - \frac{\mu }{{r_{L}^{3}}}{{{\mathbf{r}}}_{L}} + \frac{1}{M}{{P}_{L}}{{{\mathbf{r}}}_{{L'}}},$
где ${{{\mathbf{r}}}_{L}}$ – положение верхней точки, ${{P}_{L}}$ – сила натяжения на верхнем конце троса, $M$ – масса противовеса.

В результате краевые условия для уравнения (2.1) принимают вид

(2.4)
${\mathbf{r}}(0,t) = {\mathbf{R}},\quad {\mathbf{\dot {r}}}(0,t) = 0,\quad {\mathbf{r}}(L,t) = {{{\mathbf{r}}}_{L}}(t),\quad {\mathbf{\dot {r}}}(L,t) = {{{\mathbf{\dot {r}}}}_{L}}(t)$

Продифференцируем краевые условия в точке $L$ по времени, подставим вместо вторых производных их выражения из уравнений движения (2.1) и (2.3) и спроецируем найденные соотношения на касательный вектор к тросу ${\mathbf{r}}{\kern 1pt} '(L,t)$ в его конце. Получим краевое условие для величины натяжения троса в конечной точке

(2.5)
$P{\kern 1pt} '{\kern 1pt} (L,t) + \frac{{\rho (L){{P}_{L}}(t)}}{M} = 0$

Для завершения постановки задачи нужно задать функцию распределения линейной плотности $\rho (s)$ вдоль троса. Ее неоднородность существенна для динамики системы, но явный вид $\rho (s)$ не понадобится для описываемого ниже анализа. Для нескольких моделей КЛ он приведен в [9, уравнения (2), (8) и (11)].

Вертикальная равновесная конфигурация. Уравнения (2.1)(2.4) имеют частное стационарное решение, соответствующее равновесию троса, вытянутого вдоль радиус-вектора

(2.6)
$\begin{gathered} {{x}_{0}}(s) = r,\quad {{y}_{0}}(s,t) = 0,\quad {{z}_{0}}(s) = 0 \\ {{x}_{0}}(L) = \ell ,\quad {{y}_{0}}(L) = 0,\quad {{z}_{0}}(L) = 0 \\ \end{gathered} $

В соответствии с введенными ранее обозначениями, $r = r(s) = s + R$. Дополнительно введено обозначение $\ell = L + R$ (расстояние от центра Земли до конца троса).

Индексом 0 обозначаем величины, относящиеся к равновесному состоянию.

Натяжение троса ${{P}_{0}}(s)$ в равновесном состоянии, согласно уравнениям (2.1) и (2.5), определяется дифференциальным соотношением

(2.7)
${{P}_{{0'}}}(s) = - \rho (s)\left( {{{\omega }^{2}}r - \frac{\mu }{{{{r}^{2}}}}} \right)$
с граничным условием

(2.8)
${{P}_{{0L}}} = M\left( {{{\omega }^{2}}\ell - \frac{\mu }{{{{\ell }^{2}}}}} \right)$

Здесь $\omega = \left| \omega \right|$ – частота суточного вращения Земли. Заметим, что в выражение (2.8) не входит значение $\rho (L)$ плотности троса на конце.

Пусть ${{L}_{{gs}}}$ – высота геостационарной орбиты. Тогда ${{P}_{{0'}}}({{L}_{{gs}}}) = 0$, а условие ${{P}_{{0L}}} > 0$ равносильно тому, что $L > {{L}_{{gs}}}$. Это условие хорошо известно [2]: для обеспечения натянутого состояния троса противовес должен находиться за геостационаром.

3. Линеаризация уравнений. Рассматривая теперь равновесное состояние как невозмущенное, запишем уравнения близких к нему движений в линейном приближении. Для этого введем малые вариации и запишем переменные, входящие в исходные уравнения, в виде ${\mathbf{r}}(s,t) = {{{\mathbf{r}}}_{0}}(s) + \delta {\mathbf{r}}(s,t)$ и (учитывая, что ${{{\mathbf{\dot {r}}}}_{0}} = 0$) ${\mathbf{\dot {r}}}(s,t) = \dot {\delta }{\mathbf{r}}(s,t)$.

Из условия нерастяжимости троса (2.2), с учетом краевых условий (2.4) получим $\delta x(s,t) = 0$, т.е. в линейном приближении вариации как радиального смещения, так и радиальной скорости относительно вертикального положения равновесия отсутствуют. Введем вектор вариаций в трансверсальном к тросу направлении

$Y(s,t) = \left( {\delta y(s,t),\delta z(s,t)} \right)$
и положим

${{Y}_{L}}(t) = \left( {\delta {{y}_{L}}(t),\delta {{z}_{L}}(t)} \right),\quad {{Y}_{{L'}}}(t) = Y{\kern 1pt} '{\kern 1pt} (L,t)$

Введем также обозначения (напомним: $r = s + R$)

(3.1)
${{F}_{z}}(s) = - \frac{\mu }{{{{r}^{3}}}},\quad {{F}_{y}}(s) = {{F}_{z}}(s) + {{\omega }^{2}}$
и $F(s) = {\text{diag}}({{F}_{y}},{{F}_{z}})$ ($2 \times 2$ диагональная матрица).

Линеаризация уравнений (2.1), (2.3) дает систему

(3.2)
$\ddot {Y} = F(s)Y + \frac{1}{{\rho (s)}}\left( {{{P}_{0}}Y{\kern 1pt} '{\kern 1pt} } \right),\quad {{\ddot {Y}}_{L}} = {{F}_{L}}{{Y}_{L}} - \frac{{{{P}_{{0L}}}}}{M}{{Y}_{{L'}}}$
с краевыми условиями, соответствующими (2.4)

(3.3)
$Y(0,t) = \dot {Y}(0,t) = 0,\quad Y(L,t) = {{Y}_{L}}(t)$

Коэффициенты ${{P}_{0}}(s)$, ${{P}_{{0L}}}$ определены уравнениями (2.7), (2.8).

Итак, в линейной постановке смещения $\delta y(s,t)$ и $\delta z(s,t)$ разделяются. Скалярные линейные краевые задачи для $\delta y$ и $\delta z$ отличаются только постоянным сдвигом коэффициентов в силу (3.1). Для определенности мы подробно рассмотрим задачу для компоненты $\delta y$.

4. Задача Штурма–Лиувилля. Разделение переменных. Будем искать частные решения краевой задачи (3.2), (3.3) для экваториальной компоненты $\delta y$ в виде

(4.1)
$\delta y(s,t) = S(s)T(t)$

Посредством обычной процедуры разделения переменных находим

(4.2)
$\frac{{\ddot {T}}}{T} = {{F}_{y}}(s) + \frac{{({{P}_{0}}(s)S{\kern 1pt} '{\kern 1pt} ){\kern 1pt} '}}{{\rho (s)S}} = - \lambda ,\quad - {\kern 1pt} {{F}_{y}}(L) + \frac{{{{P}_{{0L}}}S{\kern 1pt} '{\kern 1pt} (L)}}{{MS(L)}} = \lambda ,$
где $\lambda $ – некоторая постоянная. Оказывается (см. ниже), что $\lambda > 0$, так что параметр $\lambda $ является квадратом частоты моды колебаний, соответствующей частному решению (4.1). Форма по $s$ такой моды описывается линейным дифференциальным уравнением второго порядка
(4.3)
$\left( {{{P}_{0}}(s)S{\kern 1pt} '} \right) + {{F}_{y}}(s)\rho (s)S = - \lambda \rho (s)S$
с однородными краевыми условиями

(4.4)
$S(0) = 0,\quad {{P}_{{0L}}}S{\kern 1pt} '(L) = M({{F}_{y}}(L) + \lambda )S(L)$

В результате для определения функции $S(s)$ получаем одномерную задачу Штурма–Лиувилля (4.3), (4.4) cо спектральным параметром $\lambda $. Отличие ее постановки от классической состоит в том, что здесь спектральный параметр входит также во второе краевое условие (4.4).

Классическая теории задачи Штурма–Лиувилля дает обоснование метода разделения переменных: доказывается, что при выполнении определенных условий coбственные функции, удовлетворяющие заданным однородным соотношения на границе, существуют для счетного дискретного множества значений параметра $\lambda $ и образуют ортогональный базис в соответствующем гильбертовом пространстве. В рассматриваемом случае дело обстоит так же.

Спектральные свойства задачи. Будем интерпретировать задачу Штурма–Лиувилля (4.3), (4.4) как поиск спектра дифференциального оператора $D$

(4.5)
$Du = - \frac{1}{{\rho (s)}}\left( {{{P}_{0}}(s)u{\kern 1pt} '} \right) - {{F}_{y}}(s)u,$
действующего на функции из ${{C}^{2}}[0,L]$, удовлетворяющие краевым условиям, ср. (3.3), (3.2):

(4.6)
$u(0) = 0,\quad \frac{1}{{\rho (L)}}\left( {{{P}_{0}}(s)u{\kern 1pt} '} \right)(L) = - \frac{1}{M}{{P}_{{0L}}}u{\kern 1pt} '(L)$

Множество (вещественное линейное пространство) таких функций обозначим $C_{0}^{2}$. (Для того, чтобы получить оператор на области определения, не зависящей от $\lambda $, приходится использовать краевое условие со второй производной.) Превратим $C_{0}^{2}$ в предгильбертово пространство, введя на нем скалярное произведение по аналогии с ([11], гл. II, приложение III)

(4.7)
$(u,v) = \int\limits_0^L {u(s)v(s){\kern 1pt} ds} + Mu(L)v(L)$

Интегрируя по частям с учетом (4.6), получаем

(4.8)
$(Du,v) = \int\limits_0^L {\left( {{{P}_{0}}u{\kern 1pt} 'v{\kern 1pt} '\; - {{F}_{y}}\rho uv} \right){\kern 1pt} ds} - M{{F}_{y}}(L)u(L)v(L)$

Отсюда видно, что оператор $D$ симметричен в $C_{0}^{2}$, так что собственные функции данной задачи с разными собственными значениями $\lambda $ ортогональны относительно введеного скалярного произведения.

Счетность множества собственных значений следует из описанного ниже алгоритма их вычисления.

При математически строгом анализе задачи необходимо доказать полноту системы собственных функций оператора $D$ в гильбертовом пространстве – пополнении пространства $C_{0}^{2}$ по норме, порожденной скалярным произведением (4.7). Для задач Штурма–Лиувилля такого типа это сделано, например, в [12].

Проверим, что все собственные значения оператора $D$ положительны. Для этого докажем, что $(Du,u) > 0$ для ненулевой функции $u$.

Полагая $v = u$ в (4.8) и используя соотношения (2.8) и (3.1), получаем

$(Du,u) = \int\limits_0^L {\left( {{{P}_{0}}u{\kern 1pt} {{'}^{2}}\; - {{F}_{y}}\rho {{u}^{2}} - ({{r}^{{ - 1}}}{{P}_{0}}{{u}^{2}}){\kern 1pt} '} \right)ds} $

Преобразуем подинтегральное выражение, используя (2.7), (3.1) и тождество $u{\kern 1pt} {{'}^{2}}\; - ({{u}^{2}}{\text{/}}r){\kern 1pt} ' = (u{\kern 1pt} '\; - u{\text{/}}r{{)}^{2}}$:

${{P}_{0}}u{\kern 1pt} {{'}^{2}}\; - {{F}_{y}}\rho {{u}^{2}} - \left( {\frac{{{{P}_{0}}{{u}^{2}}}}{r}} \right) = {{P}_{0}}u{\kern 1pt} {{'}^{2}}\; + \frac{{{{P}_{{0'}}}{{u}^{2}}}}{r} - \left( {\frac{{{{P}_{0}}{{u}^{2}}}}{r}} \right) = {{P}_{0}}{{\left( {u{\kern 1pt} '\; - \frac{u}{r}} \right)}^{2}}$

Поскольку ${{P}_{0}}(s) > 0$ при $0 < s \leqslant L$ (сила натяжения положительна), имеем $(Du,u) \geqslant 0$. Равенство $(Du,u) = 0$ могло бы иметь место только если $ru{\kern 1pt} '\; - u \equiv 0$, то есть при $u{\text{/}}r = {\text{const}}$, но тогда из граничных условий $u(0) = 0$, $r(0) = R$ следует $u \equiv 0$.

Сказанное выше относилось к задаче Штурма–Лиувилля для компоненты $\delta y$ малого отклонения троса в экваториальной плоскости. Для компоненты $\delta z$, в силу (3.1), справедливо следующее. Задачи Штурма–Лиувилля для обеих компонент имеют общий набор собственных функций, и если $\omega _{n}^{{{\text{eq}}}} = \sqrt {{{\lambda }_{n}}} $ – собственная частота $n$-й экваториальной моды, то собственная частота соответствующей меридиональной моды (пространственная форма которой такая же) равна

(4.9)
$\omega _{n}^{{{\text{mer}}}} = \sqrt {{{{(\omega _{n}^{{{\text{eq}}}})}}^{2}} + {{\omega }^{2}}} $

В частности, при любом распределении массы троса $\omega _{n}^{{{\text{mer}}}} > \omega $ и, значит, периоды всех собственных колебаний в меридиональной плоскости меньше суток.

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

1) большой геометрический размер системы и сильная неоднородность троса, что обусловливает сильное изменение коэффициентов ${{P}_{e}}(s)$, $\rho (s)$ и ${{F}_{Y}}(s)$;

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

Подстановка Прюфера. Опишем метод решения краевой задачи, использующий переход к полярным координатам на фазовой плоскости. Аналогичный метод для стандартной задачи Штурма–Лиувилля рассмотрен в [13] и назван тригонометрической прогонкой.

Для удобства дальнейших выкладок сделаем подстановку $S(s) = rV(s)$, предложенную А.В. Черновым. Введем обозначения $Q = {{P}_{0}}{{r}^{2}}$, $J = \rho {{r}^{2}}$. Физический смысл переменной $V$ – угол отклонения точки троса от вертикали (отсчитываемый от центра Земли). Функция $J(s)$ – момент инерции элемента троса, а $(QV{\kern 1pt} '){\kern 1pt} '$ – момент силы (относительно центра Земли), создаваемой натяжением троса. Задача Штурма–Лиувилля (4.3), (4.4) примет вид

(5.1)
$(Q(s)V{\kern 1pt} '){\kern 1pt} ' = - \lambda J(s)V,\quad V(0) = 0,\quad {{P}_{{0L}}}V{\kern 1pt} '(L) = \lambda MV(L)$

Введем полярные координаты $q$, $\varphi $ (подстановку Прюфера):

(5.2)
$qcos\varphi = f(s)V{\kern 1pt} ',\quad qsin\varphi = g(s)V$

Функции $f$ и $g$ пока не фиксируем. Дифференцируя по $s$ с учетом (5.1), имеем

$\begin{gathered} q{\kern 1pt} '\cos \varphi - q\varphi {\kern 1pt} '\sin \varphi = \left( {\frac{{f(s)}}{{Q(s)}}} \right)\frac{{Q(s)}}{{f(s)}}qcos\varphi - \frac{{f(s)}}{{Q(s)}}\frac{{J(s)}}{{g(s)}}\lambda qsin\varphi \\ q{\kern 1pt} '\sin \varphi + q\varphi {\kern 1pt} '\cos \varphi = \frac{{g{\kern 1pt} '{\kern 1pt} (s)}}{{g(s)}}qsin\varphi + \frac{{g(s)}}{{f(s)}}qcos\varphi , \\ \end{gathered} $
откуда следует
(5.3)
$\frac{{q{\kern 1pt} '}}{q} = \left( {log\frac{f}{Q}} \right){{\cos }^{{\text{2}}}}{\kern 1pt} {\varphi } + (logg){\kern 1pt} '{{\sin }^{{\text{2}}}}{\kern 1pt} {\varphi } + \left( {\frac{g}{f} - \lambda \frac{{fJ}}{{gQ}}} \right)cos\varphi sin\varphi $
и

(5.4)
$\varphi {\kern 1pt} ' = \lambda \frac{{fJ}}{{gQ}}{{\sin }^{2}}{\kern 1pt} {\varphi } + \frac{g}{f}{{\cos }^{2}}{\kern 1pt} {\varphi } + \left( {log\frac{{Qg}}{f}} \right)sin\varphi cos\varphi $

В силу (5.1) граничные условия для уравнения (5.4) имеют вид

(5.5)
$\varphi (0) = 0,\quad \varphi (L) = \frac{{{{P}_{{0L}}}g(L)}}{{\lambda Mf(L)}}$

Назовем уравнение (5.4) определяющим уравнением. Спектр задачи Штурма–Лиувилля состоит из таких значений $\lambda $, при которых краевая задача (5.4), (5.5) имеет решение. Когда $\lambda $ и $\varphi (s)$ уже вычислены, функция $q(s)$ находится интегрированием уравнения (5.3); после этого соответствующая собственная функция $V(s)$ определяется посредством второго уравнения (5.2).

Произвол в выборе функций $f$ и $g$ в подстановке Прюфера можно использовать для получения краевой задачи (5.4), (5.5) в наиболее удобной для определенных целей форме. Рассмотрим два конкретных варианта.

Вариант I. $f(s) = Q(s)$, $g = {\text{const}} = \lambda M{{\ell }^{2}}$. Определяющее уравнение имеет вид

(5.6)
$\varphi {\kern 1pt} ' = \lambda a(s){{\cos }^{2}}\varphi + b(s){{\sin }^{2}}\varphi ,$
где $a(s) = M{{\ell }^{2}}{\text{/}}Q(s)$, $b(s) = J(s){\text{/}}(M{{\ell }^{2}})$. Эта форма удобна тем, что правое граничное условие $\varphi (L) = 1$ не зависит от $\lambda $.

Для вычисления частоты главной моды можно разделить обе части уравнения (5.6) на ${{\cos }^{2}}\varphi $ и получить краевую задачу для уравнения Риккати относительно $\tau (s) = \varphi (s)$,

$\tau {\kern 1pt} ' = \lambda a(s) + b(s){{\tau }^{2}},\quad \tau (0) = 0,\quad \tau (L) = 1$

Отсюда, полагая $A = \int_0^L a(s){\kern 1pt} ds$, $B = \int_0^L b(s){\kern 1pt} ds$, получаем оценку

(5.7)
$\frac{{1 - B}}{A} < {{\lambda }_{0}} < \frac{1}{A}$

С вычислительной точки зрения, недостаток этого варианта преобразования Прюфера, особенно при больших $\lambda $, состоит в том, что функция $\varphi (s)$, определяемая уравнением (5.6), имеет характерную ступенчатую структуру, в которой чередуются участки с большим и малым наклоном к оси абсцисс. (Наглядное представление о поведении $\varphi (s)$ можно получить, рассмотрев решение уравнения (5.6) с постоянными коэффициентами $a$ и $b$.) Помимо этого, в нашей задаче коэффициенты изменяются на несколько порядков на интервале интегрирования.

Вариант II. $f = \sqrt {Q{\text{/}}\lambda J} $, $g = 1$. Определяющее уравнение таково:

(5.8)
$\varphi {\kern 1pt} ' = \sqrt {\frac{{\lambda J}}{Q}} + \frac{1}{4}(\log (JQ)){\kern 1pt} '\sin 2\varphi ,$
с правым граничным условием
(5.9)
$\varphi (L) = \eta {\kern 1pt} {\text{/}}\sqrt \lambda ,$
где $\eta = {{M}^{{ - 1}}}\sqrt {{{P}_{0}}(L)\rho (L)} $.

Можно еще упростить краевую задачу, заменив независимую переменную $s$ на $z = \int_0^s \sqrt {J{\text{/}}Q} {\kern 1pt} ds$ и обозначив $\varphi (z(s)) = \varphi (s)$, $Z = z(L)$. Получим

(5.10)
$\frac{{d\varphi }}{{dz}} = \sqrt \lambda + \frac{1}{4}\frac{d}{{dz}}ln(JQ)sin2\varphi ,\quad \varphi (Z) = \pi n + (\eta {\text{/}}\sqrt \lambda )$

Если $\lambda $ велико, то угловая переменная $\varphi $ изменяется (в нулевом приближении) равномерно по $z$, и на интервале периодичности правой части по $\varphi $ изменение коэффициента при $sin2\varphi $ – величина порядка ${{\lambda }^{{ - 1/2}}}$. Этим обстоятельством можно воспользоваться для построения полуаналитического алгоритма вычисления собственных мод высокого порядка. Из (5.10) нетрудно также вывести, что асимптотика собственных частот в экваториальной плоскости имеет вид

(5.11)
$\omega _{n}^{{{\text{eq}}}} = \frac{{\pi n}}{Z} + O\left( {\frac{1}{n}} \right)$

Из выражений (4.9) и (5.11) следует, что разность частот для меридиональной и экваториальной мод одного порядка имеет асимптотику

$\omega _{n}^{{{\text{mer}}}} - \omega _{n}^{{{\text{eq}}}}\sim \frac{{{{\omega }^{2}}Z}}{{2\pi n}},$
так что правая часть формулы (5.11) с такой же точностью описывает асимптотику собственных частот в меридиональной моде.

Случай равнонапряженного троса. Модель равнонапряженного троса [1] характеризуется соотношением

(5.12)
${{P}_{0}}(s) = \tau \rho (s),$
где $\tau = {\text{const}}$. Модель полностью определяется, помимо физических констант, двумя параметрами: длиной троса $L$ и напряжением троса $\tau $, имеющим размерность квадрата скорости.

При условии (5.12) уравнение (2.7) принимает вид

$\tau (log\rho ){\kern 1pt} ' = - {{\omega }^{2}}r + \mu {{r}^{{ - 2}}}$

Для коэффициентов в уравнении (5.8) получаем простые выражения: $J{\text{/}}Q = 1{\text{/}}\tau $ и

(5.13)
$\frac{1}{4}(log(JQ)){\kern 1pt} ' = \frac{1}{2}(log\rho {{r}^{2}}){\kern 1pt} ' = \frac{1}{r} + \frac{1}{{2\tau }}\left( {\frac{\mu }{{{{r}^{2}}}} - {{\omega }^{2}}r} \right)$

Из соотношения (2.8) находим выражение для параметра $\eta $ в (5.9):

(5.14)
$\eta = {{\tau }^{{ - 1/2}}}({{P}_{{0L}}}{\text{/}}M) = {{\tau }^{{ - 1/2}}}({{\omega }^{2}}\ell - \mu {\text{/}}{{\ell }^{2}})$

Для наглядности вместо параметра $\tau $ можно использовать параметр, имеющий размерность длины, а именно, эквивалентную разрывную длину ${{L}_{p}}$, связанную с $\tau $ формулой $\tau = {{g}_{0}}{{L}_{p}}$, где ${{g}_{0}} = 10$ м/c2 – верхняя оценка ускорения свободного падения вблизи поверхности Земли. Имеет место соотношение $k\tau = {{\tau }_{b}} = {{g}_{0}}{{L}_{b}}$, где ${{\tau }_{b}}$ – удельная прочность на разрыв материала троса, $k > 1$ – коэффициент, называемый запасом прочности, а ${{L}_{b}}$ – так называемая разрывная длина. За подробностями отсылаем читателя к статье [9].

Пример расчета собственных частот. Для примера взяты следующие значения, приведенные в табл. 1.

Таблица 1.
Физические параметры Земли
Гравитационный параметр Земли $\mu $ $3.986 \times {{10}^{{14}}}\;{{м}^{{\text{3}}}}{\text{/}}{{{\text{с}}}^{2}}$
Скорость вращения Земли $\omega $ $7.292 \times {{10}^{{ - 5}}}\;{{с}^{{ - 1}}}$
Радиус Земли $R$ $6.378 \times {{10}^{5}}\;{\kern 1pt} м$
Определяющие параметры модели
Длина троса $L$ $8.0 \times {{10}^{7}}\;{\kern 1pt} м$
Напряжение троса $\tau $ $3 \times {{10}^{7}}\;{{м}^{2}}{\text{/}}{{{\text{с}}}^{2}}$
Вспомогательные (зависимые) параметры
Эквивалентная разрывная длина ${{L}_{p}}$ $3 \times {{10}^{6}}\;м$
Отношение масс троса и баланса ${{M}_{{тр}}}{\text{/}}M$ 1.288
Параметр в гр. усл. (5.9) $\eta $ $7.41 \times {{10}^{{ - 5}}}$
Коэф. в асимптотике частот (5.11) $\pi {\text{/}}Z$ $2.15 \times {{10}^{{ - 4}}}\;{{с}^{{ - 1}}}$

В табл. 2 приведены периоды для первых 10 мод колебаний в экваториальной плоскости $T_{n}^{{{\text{eq}}}} = 2\pi {\text{/}}\omega _{n}^{{{\text{eq}}}}$ в часах и для некоторых мод высшего порядка – в секундах. Также приведены периоды первых 5 мод колебаний в меридиональной плоскости; далее разность $T_{n}^{{{\text{eq}}}} - T_{n}^{{{\text{mer}}}}$ становится очень малой.

Таблица 2.
$n$ $T_{n}^{{{\text{eq}}}}$, ч $T_{n}^{{{\text{mer}}}}$, ч $n$ $T_{n}^{{{\text{eq}}}}$, ч $n$ $T_{n}^{{{\text{eq}}}}$, с
0 138.25 23.58 5 1.615 20 1460
1 7.818 7.431 6 1.347 40 730.2
2 3.996 3.941 7 1.155 60 486.8
3 2.679 2.662 8 1.012 80 365.1
4 2.015 2.008 9 0.8996 100 292.1

Заметим, что оценка периода (в часах) главной моды по формуле (5.7) дает $114 < T_{0}^{{{\text{eq}}}} < 158$.

Доработка статьи осуществлена при участии С.Ю. Садова и Г.В. Калачева.

Список литературы

  1. Pearson J. The orbital tower: a spacecraft launcher using the Earth’s rotational energy // Acta Astron. 1975. V. 2. № 10. P. 785–799.

  2. Белецкий В.В., Иванов М.Б., Отставнов Е.И. Модельная задача о космическом лифте // Космич. исслед. 2005. Т. 43. № 3. С. 157–160.

  3. Нуралиева А.Б. О динамике троса космического лифта // Дисс. на соиск. уч. степ. канд. физ.-мат. наук, М., 2012. 103 с. https://keldysh.ru/council/1/nuralieva-diss.pdf

  4. Садов Ю.А., Нуралиева А.Б. Нелинейные поперечные колебания троса космического лифта // Матем. модел. 2011. Т. 23. № 12. С. 3–19.

  5. Калачев Г.В., Нуралиева А.Б., Чернов А.В. Малые колебания троса космического лифта // Тр. МФТИ. 2013. Т. 5. № 4. С. 26–36. https://mipt.ru/upload/medialibrary/090/26-36.pdf

  6. Петров А.Л., Садов С.Ю. Исследование близких к маятниковым движений тяжелой нити на круговой орбите // Препринт 11. М.: ИПМ им. М.В. Келдыша, 1989. 28 с. https://istina.msu.ru/publications/article/302515268/

  7. Williams P. Dynamic multibody modeling for tethered space elevators // Acta Astron. 2009. V. 65. P. 399–422.

  8. Wang Zhenkun, Cui Naigang, Fan Youhua, Liu Dun. Modal and dynamic analysis of a tether for a nonequatorial space elevator // IEEE Access. 2018. V. 6. P. 74940–74952. https://ieeexplore.ieee.org/document/8543785

  9. Садов Ю.А., Нуралиева А.Б. Нагруженный секционированный космический лифт // Космич. исслед. 2015. Т. 53. № 3. С. 246–252.

  10. Белецкий В.В., Левин Е.М. Динамика космических тросовых систем. М.: Наука, 1990. 329 с.

  11. Тихонов А.Н., Самарский А.А. Уравнения математической физики. М.: Наука, 1977. 735 с.

  12. Fulton C.T. Two-point boundary value problems with eigenvalue parameter contained in the boundary conditions // Proc. Royal Soc. Edinburgh. 1977. V. 77A. P. 293–308.https://doi.org/10.1017/S030821050002521X

  13. Федоренко Р.П. Введение в вычислительную физику. М.: Изд. МФТИ, 1994. 528 с.

Дополнительные материалы отсутствуют.