CC BY
294
УЧЕНЫЕ ЗАПИСКИ ЦАГИ
Том XI 20 09 № 3
УДК 551.510
629.78.015.076.8:525.7
АППРОКСИМАЦИЯ МОДЕЛИ СТАНДАРТНОЙ АТМОСФЕРЫ
В. А. ЯРОШЕВСКИЙ
Предлагается модель зависимости плотности атмосферы от высоты, удобная для выполнения математического моделирования процесса входа космического аппарата в атмосферу Земли. Обсуждаются вопросы, связанные с возможностью замены этой модели упрощенной экспоненциальной моделью.
Ключевые слова: модель плотности атмосферы, математическое моделирование, вход космического аппарата в атмосферу земли, упрощенная экспоненциальная модель.
При проведении предварительных параметрических расчетов траекторий входа космического аппарата (КА) в атмосферу Земли (да и других планет) обычно используют так называемую экспоненциальную модель плотности атмосферы
р = р0ехр (-Ш ). (1)
В действительности эта модель является приближенной, главной причиной отличия этой
модели от истинной является неизотермичность земной атмосферы. Рассмотрим уравнение гид-
ростатики
1И(2)
вместе с уравнением газового состояния
Р = РЯ/. (3)
Здесь р — давление, Т — температура, Я. — газовая постоянная.
Тогда нетрудно представить истинную модель плотности атмосферы в виде
Р = Ро ехр
( н \
-| Х(Н )сіН1
(4)
V 0
где местный логарифмический градиент плотности атмосферы по высоте определяется формулой
х = -л.р=_^+_^
дИ Я§Т ёИу ’
Таким образом, этот градиент в основном обратно пропорционален температуре (первое слагаемое в сумме (5)), которая в диапазоне высот от нуля до 100 км варьируется в пределах ±20% от среднего значения. Неудобство представления зависимости А,(И) в виде (5) проявляется
в том, что зависимость Т(Н) для принятых моделей стандартной атмосферы обычно аппроксимируется кусочно-линейной функцией высоты, т. е. имеет изломы, благодаря чему зависимость А,(Н) претерпевает разрывы [1], которые достигают 20 — 25%, что не отражает реальности. Поэтому целесообразно предложить другие формулы для описания этой зависимости.
Возьмем за основу стандартную модель плотности атмосферы 1981 г., вычислим значения ln р в зависимости от высоты H (в километрах) и поместим их во второй столбец табл. 1. Аппроксимируем эту зависимость квадратным трехчленом с добавлением различного количества гармонических слагаемых:
ln р = а0 + a1 (Я/ 50 -1) + а2 (Я/ 50 -1)2 +
+b1 cos[п(Н/50 -1)] + b2 cos [2п(Н/50 -1)] +
+b3 cos [3п (Я/ 50 -1)] + b4 cos [4п (Я/ 50 -1)] +
+b5 cos [5п (Я/ 50 -1)] + c1 sin [п (Н/ 50 -1)] +
+c2 sin [2п (Я/ 50 -1)] + c3 sin [3п(Н/ 50 -1)] +
+c4 sin [ 4п(Н/ 50 -1)] + c5 sin [5п(Я/ 50 -1)]. (6)
Определяя коэффициенты методом наименьших квадратов, получим приближенные формулы для расчета значений ln р в диапазоне от нуля до 100 км с разной степенью точности. Результаты расчетов приведены в табл. 1 для различного количества гармонических слагаемых в сумме (6), (N — номер столбца в табл. 1).
Т аблица 1
Значения 1прЯН) для различных вариантов коэффициентов аппроксимации
Н, км Стандартная атмосфера N и 3 N II 4 N II 5 N II 6 N II 7
0. 0.2029 0.2859 0.2213 0.2047 0.1992 0.2000
5. -0.3060 -0.3115 -0.3135 -0.2961 -0.2923 -0.2957
10. -0.8831 -0.9740 -0.9287 0.9114 -0.9006 -0.8978
15. -1.6358 -1.6881 -1.6250 -1.6265 -1.6266 -1.6253
20. -2.4201 -2.4378 -2.3889 -2.4062 -2.4193 -2.4224
25. -3.2169 -3.2061 -3.1943 -3.2104 -3.2180 -3.2187
30. -3.9949 -3.9765 -4.0078 -4.0083 -3.9988 -3.9955
35. -4.7720 -4.7349 -4.7965 -4.7817 -4.7681 -4.7677
40. -5.5225 -5.4708 -5.5373 -5.5208 -5.5225 -5.5258
45. -6.2318 -6.1783 -6.2223 -6.2184 -6.2336 -6.2337
50. -6.8811 -6.8567 -6.8600 -6.8715 -6.8789 -6.8756
55. -7.4732 -7.5108 -7.4718 -7.4882 -7.4770 -7.4771
60. -8.0799 -8.1500 -8.0845 -8.0915 -8.0772 -8.0805
65. -8.7205 -8.7878 -8.7227 -8.7148 -8.7178 -8.7175
70. -9.3987 -9.4402 -9.4027 -9.3875 -9.4041 -9.4008
75. -10.1286 -10.1242 -10.1302 -10.1216 -10.1281 -10.1288
80. -10.8999 -10.8563 -10.9033 -10.9081 -10.8946 -10.8978
85. -11.7089 -11.6505 -11.7170 -11.7289 -11.7147 -11.7134
90. -12.5865 -12.5175 -12.5694 -12.5753 -12.5825 -12.5797
95. -13.4755 -13.4631 -13.4656 -13.4601 -13.4783 -13.4812
100. -14.4043 -14.4878 -14.4172 -14.4145 -14.4032 -14.4024
В табл. 2 приведены коэффициенты аппроксимации для расчета приближенных значений 1п р по формуле (6).
Значения коэффициентов аппроксимации
N = 3 N и 4 N и 5 N и 6 N и 7
а0 -6.3515 -6.3643 -6.3068 -6.3515 -6.3759
-7.3869 -7.3193 -7.3096 -7.3012 -7.3012
а2 -1.2547 -1.2168 -1.3876 -1.2545 -1.1817
Ь\ -0.5052 -0.4894 -0.5599 -0.5052 -0.4754
Ь2 0 -0.0062 0.0123 -0.0019 -0.0096
Ь3 0 0 -0.0172 -0.0104 -0.0068
Ь4 0 0 0 -0.0098 -0.0120
Ь5 0 0 0 0 0.0042
С1 0.2344 0.1917 0.1856 0.1803 0.1803
с2 0 0.0816 0.0846 0.0872 0.0872
с3 0 0 -0.0136 -0.0153 -0.0153
с4 0 0 0 0.0145 0.0145
с5 0 0 0 0 0
ІЛрІ , % 1 ' Ішах 9 4.5 2.8 1.7 1.5
В последней строке табл. 2 приведена максимальная по модулю относительная погрешность вычисления плотности в диапазоне от нуля до 100 км в процентах. Представляется, что учет шести гармонических слагаемых в сумме (6) (пятый столбец) позволяет получить приемлемую точность. Тогда функцию А,(Я) можно определить, вычисляя производную от функции (6):
ЦИ)« 0.1462 + 0.0555(И/50-1)-0.0117ео8[п(И/50-1)]--0.0106ео8 [ 2я(Я/ 50 -1)] + 0.0026ео8 [3п(И/ 50 -1)] --0.03528ІП [я(Я/ 50 -1)] + 0.00168ІИ [2я(Я/ 50 -1)] -
-0.00328ІИ [3л(Я/ 50 -1)]. (7)
Зависимость А,(Я), определяемая формулой (7), изображена на рис. 1. На рис. 2 иллюстрируется точность аппроксимации стандартной атмосферы (точки) суммой (6) для этого случая.
1пр
2 0
-4
8 10 -12 14 16
-
N
Рис. 1. Зависимость логарифмического градиента плотности от высоты 'к(Я) (формула (5))
О 10 20 30 40 50 60 70 80 90 100
Н, км
Рис. 2. Зависимость логарифма плотности от высоты 1п р (И) (формула (6)): точки — стандартная атмосфера
Обсудим возможности использования простейшей модели плотности, описываемой формулой (1), для исследования траекторий входа в атмосферу. Рассмотрим гипотетический спускаемый КА с массой т = 10 т, характерной площадью £ = 8 м2, аэродинамическим качеством К = 0.7 (в диапазоне гиперзвуковых скоростей), коэффициентом сопротивления сх = 2. Характеристики
этого аппарата близки к рассмотренному в [2]. Чтобы иметь возможность анализировать пространственную структуру траекторий, примем, что КА совершает спуск в атмосферу Земли с углом крена, равным 45°. Начальные условия при входе КА в атмосферу невращающейся Земли зададим в виде:
V0 = 7845 м/с,
H0 = 100 км,
0О = -0.5°
р0 = 0.000000555 кг/м3.
Добавим к уравнениям движения в местной вертикальной плоскости уравнения бокового движения в виде:
V cos 0
dt R + H
Y 2,
d Y 2 dt
V cos 0 R + H
Yr
cyaSpV siny 2mV cos 0
(8)
dy3 cyaSPV sinY
dt 2mV cos 0
Y 2.
Здесь углы У!, у2, У3, связанные соотношением у1 + у2 +Уз = 1, определяют смещение КА относительно исходной плоскости орбиты, фиксируемой при входе в атмосферу. Если считать эту плоскость экваториальной, то значения этих углов определяются формулами:
Y1 = sin ф, y2 = sin у cos ф, y3 = cos у cos ф,
(9)
где ф — широта, у — угол наклона горизонтальной проекции скорости к местной параллели. Тогда начальные условия для этих углов можно записать в виде: yj (0) = уг (0) = 0, у3 (0) = 1.
Проинтегрируем уравнения продольного и бокового движения, задавая модель плотности атмосферы в простейшем виде (1) (А, = 1/7000 м) и в виде (4), где функция ’к(И) определяется формулой (7). Конечные условия определим моментом времени, когда скорость уменьшается до 785 м/c (одна десятая от скорости входа в атмосферу). На рис. 3 сравниваются полученные
численные решения в виде зависимостей продольной перегрузки от скорости пш (V). Как видно,
различие между представленными кривыми небольшое. Причина близости этих кривых объясняется следующими соображениями. Подсчитанные траектории, соответствующие достаточно «большому» эффективному аэродинамическому качеству K cos у, в целом близки к траекториям квази-стационарного планирования, для которых в нулевом приближении выполняется условие
1 - V 2 K cos y
(10)
Рис. 3. Зависимость продольной перегрузки от скорости
Пха (У) :
К = 0.7; 0о = -0.5°;-------экспоненциальная модель плотности
Как видно, значение логарифмического градиента плотности Х(И) в эту формулу не входит и при одинаковых значениях эффективного аэродинамического качества K cos у зависимости оказываются близкими. Рассмотрим траектории КА
n
xa
с меньшими значениями K cos у. Для этого, при прочих равных условиях, уменьшим аэродинамическое качество до значений K = 0.5, 0.35, 0.25 и до 0. Соответствующие кривые изображены на рис. 4 — 7. Различие становится более ощутимым. Действительно, согласно результатам аналитических исследований [3, 1] при K = 0 и малом угле входа в атмосферу перегрузка пропорциональна , поэтому различие между кривыми более заметно.
Рис. 4. То же, что на рис. 3, но K = 0.5
п
4
3.5 3
2.5 2
1.5 I
0.5
сч; г ч\ - х 1
/ к : V :
/ * Г""Ж1 [ дк
/.. I \
1 1 \ : V
: \
т ? ; : \
Рис. 5. То же, что на рис. 3, но K = 0.35
О
1000 2000 3000 4000 5000 6000 7000 8000
V, м/с
Рис. 6. То же, что на рис. 3, но K = 0.25
Рис. 7. То же, что на рис. 3, но K = 0
Обратимся к другим параметрам траекторий и выделим «интегральные» и «локальные» параметры. К первым отнесем, например, боковую геоцентрическую дальность к, накапливаемую к моменту достижения скорости 785 м/с, ко вторым — конечное значение угла наклона траектории 0. Сопоставим результаты расчетов с использованием простейшей модели плотности атмосферы (формула (1)) и реальной (формула (4)). Некоторые результаты расчетов приведены в табл. 3.
Т аблица 3
Сравнение результатов расчетов
K 0.7 0.5 0.35 0.25 0
h°, ф.(1) 5.807 3.216 1.740 1.001
h°, ф.(4) 5.797 3.206 1.733 0.996
0°, ф.(1) -11.97 -13.65 -15.16 -16.38 -20.53
0°, ф.(4) -11.19 -12.78 -14.26 -15.36 -19.56
H f, км 27.6 26.5 25.5 24.7 21.7
Как видно, интегральные параметры вычисляются с использованием простейшей модели очень точно, а локальные параметры вычисляются с заметными погрешностями: они варьируются в зависимости от конечного значения высоты Н ^.
Приведем еще один пример, относящийся к вычислению конечного значения боковой дальности для КА с «большим» аэродинамическим качеством. Примем, что коэффициенты сопротивления и подъемной силы остаются неизменными при скорости, превышающей 3000 м/с, и далее коэффициент сопротивления уменьшается от 2 до 1, а коэффициент подъемной силы увеличивается от 1.4 до 4 при уменьшении скорости до нуля. Используем вначале модель (1) и зададим параболическую зависимость угла крена от скорости, отнесенной к скорости входа в атмосферу, определяемую тремя значениями у(1), у(0.5), у(0). Подберем эту зависимость из условия максимизации конечной боковой дальности. В работе [5] показано, что такой подход позволяет практически точно оценить маневренные возможности КА. Далее повторяем описанную процедуру уже при использовании модели (4). Результаты поместим в табл. 4.
Таблица 4 Вычисление конечного значения боковой дальности
Номер формулы У°(1) У°(0.5) У°(0) к° тах
ф. (1) 66.5 35.5 24.5 8.019
ф. (4) 66 38 19 8.039
Отсюда следует вывод, что применение простейшей модели к решению данной задачи вполне допустимо.
Различие между моделями более очевидно при рассмотрении другого класса траекторий, а именно — для траекторий крутого баллистического спуска. Увеличим (по модулю) угол входа в атмосферу до -5° и рассмотрим три варианта КА с различными значениями характерной площади S, уменьшенными и увеличенными в 10 раз, при К = 0, не придавая значения физической природе этого КА. Результаты расчетов изображены на рис. 8 —10. Как видно, для модели, описываемой формулой (4), максимальная перегрузка варьируется в достаточно больших пределах. В работе [6] показано, что основное слагаемое при вычислении этой перегрузки оказывается обратно пропорциональным температуре атмосферы на высоте, соответствующей достижению максимальной перегрузки. В свою очередь, можно усмотреть из формулы (5), что величина, обратная температуре, почти пропорциональна значению Х( Нп +1/ X), где Нп — высота, на которой достигается максимум перегрузки. Действительно, если считать функцию Т(Н) достаточно гладкой и учесть два первых члена разложения в ряд Тейлора, то формулу (5) можно переписать в виде:
/ > / У ■к
/ / / \ Л
// / / // ч \\
// Л V
\
\
\
0 1000 2000 3000 4000 5000
6000 7000 8000
V, м/с
Рис. 8. Зависимость продольной перегрузки от скорости пха {V) :
К = 0; 00 = -0.5°; 8 = 0.8 м2:
- экспоненциальная
модель плотности
Х = ^— +------------(1п Т )
ЯЯТ дНу ’
Я
1
Я
Я
Т(н-ят/я) )Т( -Vх)'
¿г
// // / / \ N \ \ Ч Л
// \ \ V \ \ V \\
// у \\
// / / // \ л \\
.. /, \і \\
\
\
У / І,''”''1' 1.
/ / "ч" \ \ \
// \ \ \ \
// // / / \ч \\ \л
О 1000 2000 3000 4000 5000 6000 7000 8000 о 1000 2000 3000 4000 5000 6000 7000 8000
У, м/с у м/с
Рис. 9. То же, что на рис. 8, но 5 = 8 м2 Рис. 10. То же, что на рис. 8, но 5 = 80 м2
В рассматриваемом случае сумма Нп +1/X составляет 26, 38 и 55 км при 5 = 0.8, 8 и 80 м2 соответственно. Легко убедиться, используя рис. 1, что максимальная перегрузка практически пропорциональна Х(Нп +1/X). В то же время, при наличии разрывов в этой функции подобные
оценки сделать затруднительно.
Работа выполнена при поддержке РФФИ, проект № 06-01-08053 офи.
ЛИТЕРАТУРА
1. Ярошевский В. А. Вход в атмосферу космических летательных аппаратов. — М.:
Наука, 1988.
2. БобылевА. В., Дядькин А. А., КобзевВ. И., Поединок В. М., РешетинА. Г., Супруненко С. Н., ЯрошевскийВ. А. Проблемы управления возвращаемым аппаратом с умеренным аэродинамическим качеством на этапе входа в атмосферу // Космические исследования. 2007. Т. 45, № 5.
3. Копнин Ю. М. Об уравнениях пространственного движения спутника, использующего аэродинамическую подъемную силу // Космические исследования. 1968. Т. 6, № 3.
4. Чепмен Д. Р. Приближенный аналитический метод исследования входа тел в атмосферу планет. — М.: Изд. иностр. лит., 1962.
5. Ярошевский В. А. Приближенное вычисление потребной и располагаемой боковой дальности, реализуемой при спуске космического аппарата в заданную точку на поверхности Земли // Космические исследования. 2005. Т. 43, № 6.
6. БобылевА. В., ЯрошевскийВ. А. Аналитическая оценка максимальной перегрузки, достигаемой при входе в атмосферу космического аппарата // Космические исследования. 1999. Т. 37, № 2.
Рукопись поступила 21/ХІ2007 г.
