Диссертация (Методика оценки долговременной эволюции техногенного засорения низких околоземных орбит при реализации активного удаления космического мусора), страница 8

Эволюционноеуравнениеиспользуетсядляпрогнозированияорбитальнойэволюциираспределения объектов КМ из подмножества D по высоте перигея. Принимаетсядопущение, что из учитываемых орбитальных параметров только высота перигеяоказывает существенное влияние на эволюцию распределения числа объектов КМпо высоте. Плотность распределения объектов КМ по высоте перигея изнекоторойгруппысопределеннымизначениямидругихпараметров(эксцентриситет, баллистический коэффициент и др.) в момент времени tобозначим как p(h,t).При расчете эволюции распределения объектов КМ по высоте перигеяучитываются следующие факторы: исходное распределение объектов КМ повысотеперигея;разбиениевсехобъектовнагруппыпопараметрам,отличающихся значениями эксцентриситета и баллистического коэффициента; ихторможение в атмосфере на высотах до 2000 км; ожидаемая интенсивностьобразования нового КМ.Высота перигея влияет на торможение объектов и изменяется под еговоздействием (высота перигея уменьшается, и КА сгорают в атмосфере навысотах около 100 км).

Баллистические коэффициенты и размеры КМ в процессеэволюции орбиты практически не изменяются. Значение эксцентриситетауменьшается под действием торможения в атмосфере, но эти изменения неоказывают существенного влияния на долгосрочную эволюцию КМ в областиНОО, так как большая часть объектов КМ имеет орбиты с малыми значениямиэксцентриситета.49Рассмотрим дискретное распределение p(h,t) на некотором заданноминтервале высот с шагом Δh.Рисунок 2.2.1 – Распределение объектов КМ p(h,t) по высоте перигеяНа рисунке 2.2.1 рассмотрены значения распределения p(h,t) для двух значенийаргумента h и h+Δh.

В этом интервале высот перигея количество объектов равноN (t )h , h + ∆ h = p (h , t ) ⋅ ∆ h.(2.1)Пусть V(h,t) и V(h+Δh,t)- скорость уменьшения высоты перигея объектов свысотой h и с высотой h+Δh.Через некоторое время, в момент времени t+Δt, распределение p(h,t)измениться. В рассматриваемом интервале высот количество объектов (2.1)изменится в результате следующих обстоятельств:1. в окрестности высоты h часть объектов снизится настолько, что их высотаперигея станет меньше высоты h, количество этих объектовN (t ,t + ∆t )(h1+)∆h = V (h ,t ) ⋅ ∆t ⋅ p (h ,t ) ,(2.2)2.

в окрестности высоты >h+Δh часть объектов снизится настолько, что их высотаперигея станет меньше высоты h+Δh, количество этих объектовN (t , t + ∆t )h+∆h = V (h + ∆h, t ) ⋅ ∆t ⋅ p(h + ∆h, t ) ,( 2)(2.3)3. добавляются вновь образованные объекты, количество этих объектовN (t , t + ∆t )h + ∆h = p(h, t )new ⋅ ∆h ⋅ ∆t .( 3)(2.4)В момент времени t+Δt суммарное количество объектов в рассматриваемоминтервале высот перигея будет равноN (t + ∆t )h ,h + ∆h = N (t )h ,h + ∆h − N (t ,t + ∆t )(h1+)∆h + N (t ,t + ∆t )(h2+ )∆h + N (t ,t + ∆t )(h3+)∆h .(2.5)50В момент времени t+Δt плотность распределения объектов по высотеперигея определяется на основе оценки (2.5)p(h, t + ∆t ) =N (t + ∆t )h,h+∆h∆h.(2.6)Соотношения (2.5) и (2.6) позволяют рассчитать при прогнозе на один шагпо времени изменение распределения числа объектов КМ по высоте перигея.Последовательное применение этих соотношений в цикле по высоте и по времениобеспечивает решение задачи прогнозирования распределения p(h,t) дляподмножества объектов D.Принимая дискретные значения Δh и Δt бесконечно малыми величинами (dhи dt) соотношения (2.5) и (2.6) можно преобразовать в дифференциальноеуравнение в частных производных, которое описывает эволюцию распределенияобъектов КМ по высоте перигея.∂p(h,t )∂p(h,t )∂V (h,t )= V (h,t ) ⋅+ p(h,t ) ⋅+ p (h,t )new .∂t∂h∂h(2.7)С учетом полного дифференциала распределение p(h,t) в точке (h,t)dp(h, t ) =∂p(h,t )∂p(h,t )⋅ dh +⋅ dt .∂h∂tС учетом того, что снижение высоты перигея V (h , t ) = − dh(2.8)dtподстановка в(2.8) выражения (2.7) приводит к уравнениюdp(h.t ) ∂V (h, t )=⋅ p(h, t ) + p(h, t )new = A(t ) ⋅ p(h, t ) + p(h, t )new .dt∂h(2.9)По своему содержанию - это эволюционное уравнение в дифференциальнойформе, а именно линейное неоднородное дифференциальное уравнение, котороепозволяет рассчитывать эволюцию во времени распределения объектов КМ повысоте перигея, а соотношение (2.5) - это решение эволюционного уравнения вдискретной форме.2.2.2 Определение скорости снижения высоты перигеяОпределения скорости уменьшения высоты перигея объектов V(h,t)необходимо для решения эволюционного уравнения и прогнозирования высотыперигея и эксцентриситета объектов при заданных других параметрах.51В области НОО одним из основных факторов, влияющих на эволюциюорбит КО, является торможение в верхних слоях атмосферы.

Аэродинамическиесилы, возникающие при движении КО, невелики в сравнении с их значениями внижних слоях атмосферы, однако большое время полета КО и диссипативныйхарактер влияния приводит к тому, что в области НОО торможение оказываетсущественноевлияниенаэволюциюорбитыКО.Действиесилыаэродинамического сопротивления на КО с массой m вызывает ускорениеa = kb ⋅ ρ ⋅Vrel2(2.10)где величинаC ⋅Skb = x[м2/кг]2m(2.11)это баллистический коэффициент, ρ - плотность атмосферы, Vrel - скоростьнабегающего потока газа, равная скорости полета КО относительно воздуха, C x безразмерный коэффициент аэродинамического сопротивления, S - характернаяплощадь КО.

Коэффициент C x зависит от многих величин: свойств материала егоповерхности, геометрической формы и ориентации КО, состава атмосферы и еепараметров. Для верхних слоев атмосферы в большинстве случаев этоткоэффициент находится в пределах 2.00 – 2.50, вследствие чего для определенияkb можно использовать просто отношение площади к массе КО. Площадь S - этоплощадь максимального сечения КО, нормального к вектору скорости Vrel . ДляориентированныхнеориентированныхКАSКОстрогоонаопределяетсяоказываетсяегопеременной,геометрией;чтодляприводиткнеобходимости использования ее средних значений.Дляоценкидолговременнойэволюциицелесообразнымявляетсяиспользование формул для возмущения большой полуоси ( δ a ) и эксцентриситета( δ e ) за виток [13], опубликованных в [27, 28].a23δ a = −4π ⋅ (k b ⋅ ρ ⋅ p ) ⋅ exp( − z ) ⋅  I 0 ( z ) + 2eI 1 ( z ) + e 2 ⋅ [I 0 ( z ) + I 2 ( z )] + ...p4δ e = −4π ⋅ (k b ⋅ ρ ⋅ p ) ⋅ exp(− z ) ⋅ I 1 ( z) + 0.5e[ I 0 ( z ) + I 2 ( z )] + e 2 ⋅ [3I 1 ( z ) + I 3 ( z)] + ...18(2.12)52δ h = (1 − e ) ⋅ δ a − a ⋅ δ e = −4π ⋅ (kb ⋅ ρ ⋅ p )a⋅ exp( − z ) ⋅1+ e⋅ {I 0 ( z ) − I1( z ) + e ⋅ [I1( z ) − 0.5 ⋅ I 0 ( z ) − 0.5 ⋅ I 2 ( z )] + ...}Где ρ - плотность атмосферы в перигее, а и p –большая полуось ифокальный параметр орбиты, z = a ⋅ e / H , I j ( z ) - функции Бесселя мнимогоаргумента порядка j .Множитель (k b ⋅ ρ ⋅ p ) характеризует уровень атмосферных возмущений иявляется безразмерным.

Произведение exp( − z ) ⋅ {...} учитывает влияние формыорбиты. Для круговой орбиты (е=0) оно равно 1. При увеличении эксцентриситетаего значение уменьшается.Формула (2.12) является приближенной [13]. В фигурных скобках нетслагаемых, пропорциональных e2, e3 и т.д. Формула (2.12) не учитывает такжевлияние "вздутия" атмосферы и отличие орбиты КО от эллипса. Суммарнаявеличинаэтихпогрешностейимеетпорядок10%.Такойжепорядокпогрешностей имеют расчетные значения плотности атмосферы.Основные погрешности формулы (2.12) для расчета эволюции космическогомусора связаны с большим разбросом возможных значений баллистическихкоэффициентов.

Этот разброс достигает 4-го порядка и более. Однако, даннымиформулами можно пользоваться для решения задачи долгосрочного прогнозаснижения высоты перигея [13].532.2.3 Модель верхней атмосферы и солнечной активностиДля определения скорости уменьшения высоты перигея и эволюциибольшой полуоси и эксцентриситета КО необходимо рассчитывать значенияплотности атмосферы с учетом 11-летнего цикла солнечной активности, которыйоказывает на нее наибольшее влияние.

Плотность атмосферы ρ являетсяпространственно-временной функциейρ = f (h, α , δ , F10.7 , a p , Ci , i = 1,2 ..) .(2.13)Модели верхней атмосферы задают конкретный вид этой функции [29].Основными аргументами динамических моделей верхней атмосферы являютсяследующие величины: h - высота точки над поверхностью Земли, F1 0 .7 - индекссолнечной активности, равный интенсивности радиоизлучения Солнца (1 SolarFlux Units (SFU) = 10-22 Вт/(м2∙Гц) на волне 10.7 см), α , δ - сферическиекоординаты точки в геоцентрической инерциальной системе координат, a p (илиK p ) - индекс, характеризующий геомагнитную активность; C i - параметрымодели.

Высота h зависит от радиуса-вектора (r) и широты точки (φ)(h = r − RE ⋅ 1 − ε ⋅ sin 2 ϕ).(2.14)Наиболее существенным аргументом в модели (2.13) является высота. Вотносительно небольшом диапазоне высот зависимость плотности от высотыаппроксимируется выражением h − h0  ,H ρ ( h ) = ρ ( h0 ) ⋅ exp  −(2.15)где H - так называемая высота однородной атмосферы (шкала высот).Влияние координатα и δ на плотность атмосферы связано с суточнымэффектом.

Происхождение этого названия объясняется различной степеньюразогрева верхней атмосферы в дневное и ночное время. Для решения задачидолгосрочного прогнозирования данными эффектами можно пренебречь, т.к. ихвлияние на порядок меньше основных факторов, таких как высота и солнечнаяактивность.54Для расчета плотности верхней атмосферы используется модель ГОСТ Р25645.166-2004 [29]:=exp (=+(1 +h++h ++h ++)h +h +h)(2.16)Коэффициенты модели (2.16) выбираются в зависимости от индексовсолнечной активности, высоты и периода времени, для которого рассчитываетсяплотность верхней атмосферы. При прогнозировании на большие интервалывремени используются среднегодовые значениягод."– индекса солнечнойактивности, равного плотности потока радиоизлучения Солнца на длине волны10.7 см (на частоте 2800), 10-22 Вт(м2∙Гц) и средние значения коэффициента a p ,т.к. он слабо влияет на значение плотности.Для описания циклов солнечной активности используется модель ГОСТ25645.302-83 [30].

Среднегодовое значение индексаформулегде, a=0.895 10)[Вт/( мгод."=& год + ',Гц)] , b= 61.17 10)год."рассчитывается по(2.17)[Вт/( мГц)], & год –среднегодовое значение чисел Вольфа. Методика расчета среднегодовыхзначений чисел Вольфа состоит из расчета чисел внутри 11-ти летнего цикласолнечной активности и расчета чисел каждого следующего цикла солнечнойактивности.Для расчета среднегодовых чисел Вольфа внутри текущего цикла, каждоеследующее число & год вычисляют через предыдущее по уравнениям линейныхрегрессий, приведенных в [30].