SIMULATION OF THE SPACE ELEVATOR'S OSCILLATIONS DURING PAYLOAD LIFTING INTO ORBIT

Full Text

На сегодняшний день проблема поиска альтернативных способов доставки грузов в космос является актуальной научной зада- чей. Одним из таких способов может стать космический лифт. Идея состоит в том, что- бы протянуть трос от поверхности Земли до космической станции за геостационарной орбитой и доставлять груз в космос на подъ- ёмнике. При этом энергия тратится только на увеличение высоты груза. Увеличение кине- тической энергии груза происходит автома- тически за счёт вращения всей конструкции вместе с Землей. Гравитационные и центро- бежные силы держат конструкцию в натяну- том состоянии [1; 2]. Существует несколько концепций космического лифта [3-5]. Наи- более проработанным является проект, пред- ложенный Bradley Edwards [5]. Данное ис- следование стало широко известным, и поч- image © Ледков А. С., Пикалов Р. С., 2016. Ледков Александр Сергеевич (ledkov@inbox.ru), доцент кафедры теоретической механики; Пикалов Руслан Сергеевич (pickalovrs@gmail.com), аспирант кафедры теоретической механики Самарского университета, 443086, Россия, г. Самара, Московское шоссе,34. ти все последующие работы по космическо- му лифту базировались на его концепции. Стоит отметить, что концепция Edwards имеет ряд не до конца проработанных во- просов, в частности почти не затрагиваются вопросы организации систем контроля за со- стоянием оборудования и управления дви- жениями лифта, обеспечивающие живучесть конструкции [6]. Однако эта концепция се- годня остается наиболее проработанной. В исследованиях [6-8] предложена концеп- ция нагруженного космического лифта, оп- ределяющая более реалистичные черты и характеристики лифта. Несмотря на большое количество работ по тематике космического лифта, до сих пор слабоизученным остаётся вопрос моделирова- ния и исследования движения космического лифта с учётом движения подъёмника. Данной проблемой, в разное время, занимался ряд учёных [4-17]. В работах [11-14] было иссле- довано воздействие, оказываемое подъёмни- ком на динамику космического лифта. Уста- новлено, что движение подъёмника приводит к возникновению колебаний лифта, которые не затухают после его остановки. В работах [11; 13] использовались математические моде- ли, в которых трос представлялся в виде стержня. В исследованиях [12; 14-16] по- строены более сложные многоточечные моде- ли. Для исследования поперечных колебаний троса в работах [6-8] использовалась контину- альная модель троса. Результаты [11-14] пока- зывают, что можно снизить или полностью убрать негативное влияние от движения подъ- ёмника. Были предложены разные схемы подъёма, позволяющие в разной степени уменьшить негативное влияние от движения подъёмника. Тем не менее, данная проблема требует дальнейшего изучения. Представленная работа является про- должением наших исследований, где для изучения динамики лифта была предложена двухзвенная математическая модель, в кото- рой трос рассматривался в виде двух нерас- Многоточечная модель лифта Рассмотрим механическую систему пространственного неэкваториального кос- мического лифта (рис. 1). Она состоит из троса, противовеса и подъёмника. Движение происходит в Ньютоновском гравитацион- ном поле. Влияние атмосферы, солнечного ветра, луны и других возмущаемых факторов отсутствует. Трос представим как набор ма- териальных точек, соединённых между со- бой невесомыми упругими стержнями, их число N +1. Масса точек троса, как и пло- щадь поперечного сечения стержней, зави- сит от расстояния до конца. Подъёмник рас- сматривается как материальная точка, дви- жущаяся вдоль троса. Введём системы координат: подвиж- ную вращающуюся вместе с Землёй систему тяжимых стержней переменной длины [18; Oxyz , неподвижную систему Ox0 y0 z0 . Ось 19]. Первый соединял точку закрепления троса и движущийся подъёмник, второй - подъёмник и противовес. Данная модель по- зволяла учитывать изгиб троса, но не учиты- Ox направим по направлению местной вер- тикали, Oz направлена в сторону северного полушария, Oy дополняет систему до пра- вала его упругость. Для получения более полной картины динамики лифта строится вой. Oz0 - направлена по оси вращения Зем- более сложная многоточечная модель кос- ли. Оси Ox0 и Oy0 направлены на непод- мического лифта, учитывающая упругость вижные звёзды и вместе с Oz0 составляют троса. За основу взята математическая мо- дель, предложенная Paul Williams [12]. По- строенная модель может быть использована при моделировании динамики космических тросовых систем при изучении динамики троса. В рамках данной статьи особое вни- мание уделено моделированию влияния подъёмника, вопрос влияния атмосферы и правую систему координат. Начало систем координат в центре Земли, точка O . Уравнения движения точек троса Определим положение j-ой точки троса в системе координат Oxyz с помощью радиус вектора: T нецентральности гравитационного поля ос- тавлен за рамками данной работы. где rj  x j , y j ,z j  j  0,...,N  2 . , (1) image Рис. 1. Многоточечная модель космического лифта Индекс 0 соответствует точке крепления T F g j - сила гравитационного воздействия троса к Земле A: r0  RE ,0,0 , где RE эк- Земли [12]. ваториальный радиус Земли. Индексы 1.N - тросу, N + 1 - противовесу C, N + 2 - подъ- ёмнику В. Движение подвижной системы коорди- нат Oxyz относительно Ox0 y0 z0 , определяется вектором угловой скорости вращения Земли, в системе Oxyz он имеет следующий вид: Распределение массы вдоль троса Конструкция космического лифта предполагает использование троса перемен- ного поперечного сечения, что позволяет минимизировать его массу, обеспечив необ- ходимую прочность. В модели это учитыва- ется тем, что масса сегментов троса разная. Для её определения запишем выражение для ω  E sin ,0,cos T , (2) массы, зависящий от расстояния до конца троса: где E угловая скорость вращения Земли; φ - угол, определяющий широту расположе- ния точки крепления лифта; E и  считаем постоянными. image m s    As  , image где  - средняя плотность материала троса; Абсолютная скорость j-ой точки в под- As функция площади поперечного сече- вижной системе координат будет опреде- ляться по следующей формуле [20]: ния троса, зависящая от расстояния до конца троса [3]. Она определяется по формуле: v j  rj  ω  rj . (3) image   g R  3 R s  R 2  As  Am exp  0 E    R  2 G s  R image  E  2R2  Абсолютное ускорение j-ой точки равно  G  E G  ,(6) a j  rj  ω  rj   2 ω  rj   ω  ω  rj  . (4) где RE экваториальный радиус Земли; Подставляя формулы (1) и (2) пооче- редно в (3) и (4), получим: RG - радиус геостационарной орбиты. Рассмотрим участок троса, представ- ленный на рис. 2. Массу j точки будем опре- делять как интеграл от функции (6), взятый  x j  E y j cos  по параметру s в пределах j  j E j E j  v   y   x cos   z sin   jl j , j  1l j  ,  j  1,N  1 :  z j  E y j sin  ,    x  2 x cos2   2 z cos sin   y cos  m   jl j m s ds,  j  1,N  1 j E j E j E j  a  y  2 y  2  x cos  z sin   j  j1l j ,  j  j E j E j j  z  2 z sin2   2 x cos sin   y sin   j E j E j E j  . где l j  l / N  1 начальная длина одного Используя общее уравнение динамики [20], уравнения движения троса можно запи- сать в следующем виде: сегмента троса после разбиения его на части. Предполагается, что масса рассчитывается для ненагруженного троса, когда его длина соответствует недеформированной длине r  image 1 Fs  Fg   2 ω  r   ω  ω  r  троса l . Масса последнего N -го сегмента m j j j j j j ,(5) троса добавляется к массе противовеса C . За площадь поперечного сечения Aj где s j  1,N , m j масса j-ого сегмента троса; j-ого сегмента троса будем принимать значе- ние площади, вычисленное в центре этого Fj - сила натяжения; сегмента (рис. 2). image Рис.2. Распределение массы троса То есть в формуле (6) будем использовать при деформации троса, и для полноты анали- значение параметра s  jl j  l j / 2 расстоя- за её следует учитывать. Демпфирование пропорционально скорости деформации, вы- ние от точки закрепления троса к Земле до се- редины j сегмента для ненагруженного троса. числяемой как производная по времени от величины деформации (7): Силы, действующие на точки троса Вектор, определяющий смещение j j image  q j  q j точки относительно j 1 , в подвижной сис- q j l j . теме координат Oxyz имеет вид: Суммарная сила натяжения, дейст- j q j  rj  rj1 . вующая на точку, складывается из силы натяжения, действующей на сегмент троса j Удлинение j-ого элемента троса: и j  1 , с учётом знаков будет иметь сле- image image  q j  l j  / l j , image image q j  l j , дующий вид [9]:    j   0, image image q j  l j , (7) F s j  Tj1 image q j1 q j1 q j image  T j q j ,  j  1,N  . (8) где j  1,N  1, а нумерация сегментов начи- нается с 1. Первый сегмент соединяет точку креп- В рамках модели рассматривается цен- тральное Ньютоновское гравитационное по- ле. В этом случае сила гравитационного воз- j ления троса к земле  x0 , y0 ,z0  и первую точ- действия Земли, действующая на троса, будет определяться формулой: точку ку троса  x1 , y1 ,z1  , последний - конечную точку троса  xN , yN ,zN  и противовес image Fg   m j r  xN 1 , yN 1 ,zN 1  . С учетом (6) сила упругости в j эле- j 3 j rj , (9) менте троса будет определяться как: Tj  EAj j  Cj , где первое слагаемое отвечает за упругость, а второе - за демпфирование; E - модуль Юнга для материала троса; C - коэффициент демпфирования. Трение между волокнами в тросе при- где  - гравитационный параметр Земли. Подставляя силы (8) и (9) в систему (5), получим систему дифференциальных урав- нений, описывающую движение троса. Уравнения движения противовеса Движение противовеса будет описы- ваться уравнением: водит к возникновению силы демпфирова- r image  1 Fs  Fg   2 ω  r   ω  ω  r  ния, которая стремится погасить продольные колебания троса. Она неизбежно возникает m C C C C C C , (10) s где FC сила натяжения, действующая со Пусть подъёмник находится на j сег- менте троса (рис. 3). При его движении стороны троса; g вдоль троса происходит деформация послед- FC - сила тяжести. него. Вектора p1 и p2 определяются на ос- нове текущего положения подъёмника в сис- Эти силы определяются формулами: Oxyz p1 j  1 F  T s C N 1 image qN 1 теме ; - соединяет p точку троса qN 1 , и подъёмник, 2 подъёмник и j точку. C Fg   image r mC 3 C rC , Данные вектора можно найти по формулам [12]: p1  rB  rj1. где rC радиус вектор противовеса во вра-  2 p  rj  rB . щающейся системе координат; mC - масса противовеса, которая складыва- ется из массы самого противовеса и массы последнего сегмента троса. Системы (5) и (10) составляют замкну- тую систему дифференциальных уравнений движения космического лифта без учёта подъёмника. Уравнения движения подъёмника Движение подъёмника будет описы- Для определения реакции, силы трения и движущей силы используется нормальная Bnb , спрямляющая Bτb и соприкасающаяся с Bτn плоскостью [12]. При определении данных сил для удобства будем проектиро- вать их на оси естественной системы коор- динат [12] Bτnb (рис. 3). Здесь b - единич- ный вектор бинормали, n - единичный век- тор нормали, τ - единичный вектор каса- тельной, определяются формулами [9]: ваться уравнением: p  p n 1  b  p b  p  b  n b  1 2 , n  B , n  1  2 , τ  image p  p n image image image image image B  b  p b  p  b  n . r  image 1 FN  F f  FP  Fg   2 ω  r   ω  ω  r  1 2 B 2  1 2  m B B B B B B B B ,(11) Подъёмник деформирует j сегмент тро- са, меняя тем самым значение силы натяже- где mB N масса подъёмника; ния на этом участке. Изменённое натяжение FB - сила реакции троса; F f B - сила трения; будет рассчитываться следующим образом: EA FB - движущая сила; image e image 1 image image 2 image j  l P F g B - сила гравитационного воздействия T  m p  p  l j . (12) Земли. Уравнение (11) строится по схеме, предложенной в работе [12]. Это значение предполагается постоян- ным по всему локальному сегменту лифта. image Рис.3. Сегмент троса с подъёмником Формула (12) должна применяться для определения силы натяжения используемых в уравнениях соседних с подъёмником точек троса. Суммарная реакция троса на подъём- ник будет определяться двумя компонентами вектора силы натяжения: ния, действующую на трос со стороны подъ- ёмника. Для этих точек сила натяжения троса будет определяться следующими соотноше- ниями: Fs  T image p2  T image q j  F f  Fg  Fk  Fc  p  n p  n  j e p j q j B B B N 1 2 2 j image image FB  Te  p image image   n p , (16)  1 2  . (13) F s j1  T q j 1 image q j1 j 1 p1 image p  Te 1  F f j1 B B B  Fg  Fk  Fc f f Реакция троса на подъёмник проециру- ется на направление нормального вектора. , (17) Когда напряжение в сегменте равно нулю, где Fj и Fj1 силы трения для j и j  1 реакция отсутствует. Кроме того, реакция точек, которые определяются следующими равна нулю, если вектора p1 и p2 парал- формулами: лельны, то есть сегмент не деформирован. При движении подъёмника возникает image F f   p1 image q j F f сила трения, деформирующая трос. Она дей- j p  p q B ствует по касательной к тросу по направле- нию вектора τ . Отметим тот факт, что в F f   1 2 j image p2 q j 1 , (18) image F f рамках данной работы подъёмник рассмат- ривается только во время движения вдоль троса. В этом случае сила трения определя- j1 p  p q B 1 2 j 1 ется формулой: Отметим, что в отличие от работы [9], в формулах (16) и (17) добавлены силы инер- Fmax  k FB , (14) ции Кориолиса B B  image image B  и центро- N Fk  2m   r с 2 бежная сила инерции FB  mB rB , дейст- где k коэффициент трения между тросом вующие на трос со стороны подъёмника. и подъёмником. Окончательный вид силы трения опре- Собирая 3N уравнений для точек троса и добавляя к ним три уравнения от подъёмни- ка и 3 от противовеса, получим замкнутую деляется формулой: систему из 3(N  2) дифференциальных урав- f FB  sign v B  τ Fmax τ , (15) нений, описывающих динамику космического лифта с учётом движения подъёмника. где vB - вектор относительной скорости Особенности моделирования движения подъёмника движения подъёмника вдоль троса. Сила трения направлена вдоль вектора τ . При подходе подъёмника к j  1 точке величина модуля вектора p1 уменьшается, а 1.6. Уравнения движения соседних с подъёмником точек троса Наличие подъёмника на j сегменте тро- са приводит к необходимости внесения из- менений в уравнения движения соседних с при прохождении через точку он обращается в ноль (рис. 4). Это создаёт проблему, поскольку из формул (13)-(19) видно, что модуль вектора p1 входит в знаменатель, и в этом случае воз- ним j и j  1 точек. Во-первых, следует из- никает неопределённость. Подобная ситуация происходит и после прохождения подъёмника менить выражения для сил натяжения, дей- ствующих на эти точки. Во-вторых, для дан- j  1 точки, в этом случае малой величиной ных точек необходимо добавить силу тре- становится модуль вектора p2 . image image Рис. 4. Подход/отход подъёмника к точкам троса Для решения данной проблемы будем отслеживать подход/отход подъёмника к точкам троса, а также момент выполнения одного из условий: Знак минус выбран из-за того, что вектор p1 направлен в противоположную движению подъёмника сторону. Учитывая введённый кинематический закон, движущая сила подъ- FP imagep1 image   , imagep2 image   , ёмника B может быть определена так: где  - величина, определяемая из требуе-  p   FP  m r (t) 1  FN  F f  Fg   2m ω  r   m ω  ω  r  мой точности расчётов. image p B B  B  .  1  tt B B B B B B B Поменяем вектора p1 и p2 на вспомо- Анализ влияния движения гательные вектора p1н и p2н (рис. 4). Такой подъёмника на динамику лифта подход позволяет избежать появления нулей в знаменателе, а следовательно, и появления ошибок при численном моделировании. Предложенный подход отличает данную мо- Будем рассматривать космический лифт, имеющий следующие параметры: средняя плотность материала троса image   1300 кг/м3, предел прочности дель от работы [12], в рамках которой пред-    ГПа, модуль Юнга E  630 ПА. Длилагалось удалять точки, фактически разбина недеформированного троса вать заново трос при подходе подъёмника к точке троса, что приводит к дополнительным l  1, 44 108 м, максимальная площадь повозмущениям системы за счёт перераспредеперечного mC  30000 сечения Am  10 мм2, ления массы. Это может сильно сказываться на результате при моделировании троса немасса противовеса 30 000 кг, масса подъёмбольшим числом точек. ника    м. mB  1000 кг. Примем величину Кинематический способ задания Для моделирования троса используем движения подъёмника вдоль троса N  10 точек. То есть трос будет разделён на Движение подъёмника в данной работе задается кинематическим законом. В этом случае необходимо заменить уравнения дви- жения подъёмника следующим уравнением: r  r (t) p1 11 сегментов. Начальным положением троса будет его стационарное состояние, которое определено уравнениями (5) и (10) путём подстановки в них первых и вторых произ- водных, которые равны нулю, а также реше- ния получившейся системы относительно image p B B x 1 , j  j  1, N  1 . Для заданных параметров полная масса где rB (t) заданный кинематический закон 3 движения подъёмника вдоль троса. троса составила 1197,110 кг. Полученное значение соответствует значениям, приво- будет равна 1,529 108 м. Временной интер- димым в работе [12]. Отметим, что к массе противовеса добавляется масса последнего сегмента, и при данных параметрах его масса вал моделирования движения составит в обоих случаях 34 дня. Будем рассматривать   0 mC  34790  30000  64790 кг. экваториальный лифт, . В качестве rB (t) будем использовать Для сравнения будет рассчитывать на основание положения подъёмника и проти- закон, использованный в работах [18, 19], согласно которому подъёмник движется с постоянной скоростью 50 м/с. В начале и вовеса в подвижной системе коорди- нат Oxyz , значения углов отклонения от ме- конце восхождения он разгоняется либо стной вертикали: 1 для подъёмника, 2 - тормозится до заданных значений скоростей. для противовеса (рис. 5). На рис. 6-7 приведены графики зави- Сравнение результатов симости углов 1 , 2 от времени t , для мно- с двухзвенной моделью Проведём сравнение результатов, по- лученных с помощью многоточечной модели и модели, описанной в работах [18, 19]. Па- раметры системы выберем одинаковыми, указанными выше. В случае многоточечной модели, стационарному положению лифта соответствует растянутое состояние троса, то есть начальная длина троса будет отличаться от длины используемой в двухзвенной моде- ли. Для сравнения будем использовать раз- готочечной - линия а и двухзвенной [18, 19] - б модели. На рис. 8 приведён график изменения относительной длины троса от времени. Из рис. 6-7 видно, что, с одной сторо- ны, амплитуды колебаний космического лифта в обоих случаях сопоставимы, а с дру- гой, - характер колебаний существенно раз- личается. Различия обусловлены влиянием силы Кориолиса на точки троса, масса кото- рых в рассматриваемом примере существен- биение на N  10 точек, длина троса в на- но превосходит массу противовеса. чальный момент времени для этого случая image Рис. 5. Углы 1 и 1 (пояснения в тексте) image Рис. 6. График изменения угла 1 (пояснения в тексте) image Рис. 7. График изменения угла 2 (пояснения в тексте) image Рис. 8. График изменения относительной длины троса (пояснения в тексте) image Рис.9. График зависимости i N , i  1, 2 от числа N (пояснения в тексте) Из рис. 8 видно, что движение подъём- ника приводит к возникновению продольных колебаний в тросе. Сила Кориолиса перево- дит эти колебания в поперечную плоскость. Двухзвенная модель [18; 19] не учитывает описанные эффекты. Исследование влияния числа точек на колебания космического лифта Исследуем влияния числа точек, ис- пользуемых для моделирования троса, на ко- лебания космического лифта. Для этого рас- смотрим космический лифт с параметрами приведенными выше, варьируя число точек Благодарности троса N  1..20 , с шагом 1. Для сравнения Работа выполнена при поддержке Рос- результатов, для углов 1 , 2 будем вычис- сийского фонда фундаментальных исследо- ваний (15-01-01456 А). лять абсолютную погрешность: t i N   max image1N t   1 N 1 t image, i  1, 2 . Литература Aslanov V. S., Ledkov A. S. Dynam- ics of the tethered satellite system. Cambridge: Wood-head Publishing Limited. 2012. 331 p. Графики зависимости i N , i  1, 2 Белецкий В. В., Левин Е. М. Дина- от величины N представлены на рис. 9, из которого видно, что с ростом числа точек N мика космических тросовых систем. М.: Наука, 1990. 329 с. Pearson J. The orbital tower: a space- величина i N , i  1, 2 убывает. Данный craft launcher using the Earth's potential ener- рисунок позволяет определить необходимое число точек троса, исходя из требуемой точ- ности расчёта. Неравномерность графиков при малых N указывает на то, что внесение дополнительной точки сильно сказывается на изменении геометрии всей системы. gy // Acta Astronautica. 2010. Vol. 2. № 9-19. P. 785-799. Smitherman D. V. Jr. Space elevator: an advanced Earth - space infrastructure for the millennium. NASA/CP-2000-210429. 2000. 48 p. Заключение В работе построена математическая модель, описывающая динамику космиче- ского лифта с учётом движения подъёмника. Трос был представлен как набор материаль- ных точек, соединенных невесомыми упру- гими стержнями. В модели учтена перемен- ная площадь поперечного сечения троса. Проведено исследование динамики космиче- ского лифта с учётом движения подъёмника. Результаты показывают, что равномер- ное восхождение подъёмника приводит к возникновению поперечных колебаний в тро- се. Было проведено сравнение полученных результатов с полученными нами ранее [18; 19]. Установлено, что амплитуда колебаний в обоих случаях сопоставима, но картина коле- баний отличается. Отличия связаны с упруго- стью троса, которая учитывается в многото- чечной модели. Движение подъёмника при- водит к возникновению продольных колеба- ний в тросе, которые, благодаря воздействию силы инерции Кориолиса, приводят к его рас- качке в поперечном направлении. Исследовано влияние числа точек, ис- пользуемых для моделирования троса, на точность полученных результатов. Получен график зависимости абсолютной погрешно-

About the authors

Alexander Sergeevich Ledkov

Samara University

Email: ledkov@inbox.ru
443086, Russia, Samara, Moskovskoye Shosse, 34

Ruslan Sergeevich Pikalov

Samara University

Email: pickalovrs@gmail.com
443086, Russia, Samara, Moskovskoye Shosse, 34

References

Supplementary files