АСТРОНОМИЧЕСКИЙ ВЕСТНИК, 2007, том 41, № 6, с. 548-554
УДК 523.682
ПРИБЛИЖЕНИЕ НАБЛЮДАЕМОГО ДВИЖЕНИЯ БОЛИДОВ АНАЛИТИЧЕСКИМ РЕШЕНИЕМ УРАВНЕНИЙ МЕТЕОРНОЙ ФИЗИКИ
© 2007 г. М. И. Грицевич
Институт механики МГУ им. М. В. Ломоносова, Москва Поступила в редакцию 01.02.2007 г.
К настоящему времени накоплен большой фактический материал по фоторегистрации траекторий метеорных тел в атмосфере Земли. Наибольшее количество снимков выполнено четырьмя болид-ными сетями, функционировавшими в разное время в США, Канаде, Европе и Испании. Аппроксимация реальных данных теоретическими зависимостями позволяет получить дополнительные оценки, не вытекающие непосредственно из наблюдений. В данной статье предложен алгоритм подбора параметров, при которых теоретическая зависимость высоты от скорости движения болида наилучшим образом аппроксимирует данные наблюдений на светящемся участке траектории. Основным отличием от предшествующих работ является приближение заданных точек аналитическим решением уравнений метеорной физики. Вычисления по представленной модели были проведены для ряда ярких метеоров Канадской болидной сети, Прерийной сети США, а также для болида Бене-шов, одного из крупнейших, зарегистрированных Европейской сетью. Правильное математическое моделирование метеорных явлений в атмосфере необходимо для последующей оценки ключевых параметров: внеатмосферной массы, коэффициента абляции, эффективной энтальпии испарения вторгающихся тел. В свою очередь, эти данные имеют значение для ряда приложений - исследования астероидно-кометной безопасности, разработки мер планетарной защиты, а также для поиска тел, способных достичь поверхности Земли.
РАС8: 96.30.Za, 94.20.Xa
ВВЕДЕНИЕ
В монографии (Стулов и др., 1995) рассмотрена задача о прямолинейном торможении болида в верхних слоях атмосферы Земли. Сконструирована динамическая система третьего порядка, фазовыми переменными которой являются масса болида М(1), его высота над поверхностью планеты И(1) и скорость У(1). Уравнения движения в проекции на касательную и нормаль к траектории имеют вид
чения тела, ра - плотность атмосферы, g - ускорение свободного падения, Я - радиус планеты, са, сЬ -коэффициенты сопротивления и подъемной силы.
Уравнения (1)-(3) дополняются уравнением для переменной массы тела:
Н*йМ = с р V35 Н й/ 2СьРа К 5'
(4)
а V
М=г- = - В + Р 8Ш у,
а/ '
й У MV MV —- = РС08у--— С08у - Ь,
а/ я
йИ ТГ ■ ----- = у,
й 1
(1)
(2)
(3)
где В = 2 СлР^З - лобовое сопротивление, Ь =
= сьРгУ25 - подъемная сила, Р = Mg - вес тела.
Здесь М, V- масса и скорость тела, 1 - время, И - высота над поверхностью планеты, у - местный угол траектории с горизонтом, V - площадь миделева се-
где Н* - теплота сублимации, сь - коэффициент теплообмена. Предполагается, что весь тепловой поток от окружающего газа расходуется на испарение материала поверхности тела.
Уравнение (3) позволяет перейти от 1 к новой независимой переменной И. Далее переходим к безразмерным величинам следующим образом: М = = Меш, V = Уеу, И = И<у, ра = р0р, V = (И0 - высота однородной атмосферы, р0 - плотность атмосферы у поверхности планеты, е - параметры входа в атмосферу).
В метеорной физике при расчете траекторий весом тела обычно пренебрегают. Поскольку форма метеорного тела, как правило, неизвестна, коэффициент подъемной силы еь считается неизвестным, так что уравнение (2) изменения угла у не рассмат-
ривается. С учетом изложенного, уравнения движения тела принимают простой вид:
т ¿т
¿У
¿У
1
= 2 С
= 1 Р 0 к0$е р V? 2^ Ме 31п у
2
Р 0 к0$е V2 р V ?
ь Ме Н * з1п у .
(5)
т = ехр
1- V2)
1-Ц
у = 1п а + в — 1п 2, А = Е1(Р) -Б1(ру2),
ёК х) = \ еА,
(6)
(7)
где а = 1 е^ Р° 0 85 - баллистический коэффици-2 ? ^
ент, в = (1 — ц)
сь V2 2^ Н *
- параметр уноса массы.
Для получения аналитического решения уравнений (5) принимается также предположение изотермической атмосферы. Также считается, что миде-лево сечение тела и его масса связаны следующей зависимостью: ? = тЦ, ц = сопй. Величина параметра Ц характеризует возможную роль вращения в полете: Ц = 0 - вращение отсутствует, ц = 2/3 - абляция тела за счет вращения происходит равномерно по всей поверхности, так что коэффициент формы тела А сохраняется.
Тогда решение уравнений (5) с начальным условием у = го, V = 1, т = 1 имеет вид:
В = 8
\Р = 4 | \ в = 6 III
В данной работе далее в качестве теоретической зависимости высоты от скорости будет использоваться аналитическое решение (7).
Необходимые начальные данные могут быть заимствованы из опубликованных таблиц, содержащих динамические и фотометрические данные реальных наблюдений на светящемся участке траектории. В частности, такие таблицы содержат значения высоты и скорости полета - к, и V, (/ = 1, ..., п) и, где V? = V1 - начальная скорость. Задача о приближении наблюдаемого движения болида сводится к поиску значений баллистического коэффициента а и параметра уноса массы в, при которых функция (7) наилучшим образом аппроксимирует данные наблюдений. Решение ищется методом наименьших квадратов.
Отметим, что при ограниченных значениях параметра уноса массы в аналитическое решение (7) обычно заменяется более простым выражением (Стулов и др., 1995; Кулаков, Стулов, 1992):
у = 1п а — 1п (—1п V) + 0.83 в( 1 — V).
(8)
0 0.2 0.4 0.6 0.8 1.0
V
Рис. 1. Аппроксимация траектории функцией (8); сплошные линии - аналитическое решение (7), штриховые линии - приближенное решение У^) (8).
В частности, для реализации метода наименьших квадратов, в качестве пробной функции ранее использовалось это приближенное решение. Вопрос об эффективности замены решения (7) на (8) был затронут ранее в работах (Грицевич, 2006; Грице-вич, Стулов, 2006). Было отмечено, что при любом в > е/0.83 = 3.275, функция (8) имеет на интервале 0 <
< V < 1 два локальных экстремума и корректно использовать ее, в качестве приближения монотонной функции (7), возможно лишь на отрезке ^^ < V <
< 1, где vrtln - точка локального минимума функции (8) (локальный максимум в данном случае всегда оказывается левее). На рис. 1 приведено графическое сопоставление функций (7) и (8) для в = 4, 6, 8 (построены графики функций У^, в) = у(^ а, в) -- 1па).
АЛГОРИТМ ПОИСКА ОСНОВНЫХ ПАРАМЕТРОВ
Баллистический коэффициент а характеризует интенсивность торможения, так как он пропорционален отношению массы столба атмосферы с поперечным сечением $ вдоль траектории к массе тела. Параметр уноса массы в пропорционален отношению доли кинетической энергии единицы массы тела, поступающей к телу в виде тепла, к эффективной энтальпии испарения. Поэтому при подборе этих параметров основной акцент следует делать на наилучшую аппроксимацию именно того участка траектории, где торможение и абляция достаточно четко выражены. Другими словами, особую ценность представляют собой наблюдения на заключительном этапе, предшествующем полному погасанию болида. Универсальным алгоритмом решения такой задачи является применение взвешенного метода наименьших квадратов, позволяющего регулировать вклад тех или иных данных в результат. С другой стороны, расстановка весов вручную затруднила бы исследование, сделала бы его недостаточ-
4
2
X
но объективным. Заметим, что по мере приближе ния метеорного тела к поверхности Земли значения следующем виде: высоты и скорости полета уменьшаются, поэтому к корректному результату удается придти, переходя от рассмотрения самих у к г~у, т.е. путем испольова-ния зависимости (7) в следующей форме:
приводит к необходимому условию экстремума
/ п
\
2аехр(-у) - Аехр(-в) = 0,
- - - Х
А = Е1(в) -Е1^2), Е1 (х) = | е-1-.
(9)
/(в) = £ (А £ ехр(-2у,)
1 = ^ г(1 =1 Л Л
ехр (-у,)
(12)
£Аг ехр ( -Уг)
г=1
/
(Аг - (Аг)в)
/
= 0.
Если бы существовали такие значения параметров а и в, при которых все пары у, = Иг/И0, V, = V в точности удовлетворяли (9), то, очевидно, теоретическая кривая (7) с этими параметрами проходила бы через все точки наблюдений. Будем ис-
Решение (12) ищем численно, например, методом Ньютона. Для реализации этого метода в качестве начального приближения в0 можно взять в, вычисленное методом Ш3(а, в) (Грицевич, Стулов, 2006).
Численный результат практически не меняется
кать оценки для а и в, основываясь на минимизации при вклютевш в ^ наблюдаемых точ<ж тра<жто-
рии с близкими к единице значениями V, эти точки
суммы квадратов значений левой части соотноше ния (9):
64 (а, в) = £( Р, (.уг, V,-, а,в))2
(10)
г = 1
где Р, (у,, v, а, в) = 2аехр(-у) - Агехр(-в), А,- =
= Е1(в) -Е1(в V2).
Поскольку параметр уноса массы - положительное действительное число, далее используется представление интегральной экспоненты в виде ряда:
Е1( х) = с + 1пх + £
п • п!
где с
вносят близкие к нулю слагаемые в саму сумму (10) и, как следствие, в соотношение (12). Это обстоятельство дополнительно свидетельствует о целесообразности применения пробной функции в форме (9), позволяющей включать в (10) весь имеющийся наблюдательный базис (у,, V), не отбрасывая заранее точки с малым торможением и абляцией.
АНАЛИЗ ПОЛУЧЕННОГО РЕШЕНИЯ Известно, что в точке минимума функции Q4, квадратичная форма й2Ш4(а, в) положительно опре-
Я Э^а, в-делена. Ясно, что --—-- =8 £ (
-2 у,
> 0. Зна-
Эа
г = 1
= 11т
п ^ га
£ к--1п п
и=1
= 0.5772 - постоянная Эйлера.
чит, при найденных оценках для а и в, для положительной определенности остается проверить выполнение следующего условия:
Параметр а, соответствующий нулю частной Эб4(а ,в)
Э2б4 (а,в)Э2б4(а,в) ^4 (а,вГ2
Эа
производной выражением:
Эа
, определяется следующим
эв2
ЭаЭв
> 0. (13)
£
- в-у д - • А,
а = <-=---
2£
г = 1
-2у
Подстановка в (13) выражений для вторых частных производных дает достаточное условие локального минимума функции Q4 в виде:
пп
£-"2у£(((А)в - Аг)2 +
= 1 = 1
(Ш + (А, - 2 а ехр (в - у, )))((Д,>; - 2 (Д,)в + А,) >
А = -21п V + £ & (I-V ).
>
£ ехр(-у,-)(А,- - (А,)в)
2
к = 1
V, = 1
Теперь
Ш( а, в ) Эв
п Если левая часть этого неравенства оказывается
= 2 £ р (у v а В)(А - меньше правой, то квадратичная форма знакопе-гКУг, 1, > н 1 ременная, т.е. исслед
Для дальнейшего прочтения статьи необходимо приобрести полный текст. Статьи высылаются в формате PDF на указанную при оплате почту. Время доставки составляет менее 10 минут. Стоимость одной статьи — 150 рублей.