Как летает космическая ракета? (Продолжение) / Хабр

Здравствуйте, дорогие любители острых космических ощущений (хабровчане)!

В предыдущей своей публикации я посчитал траекторию космической ракеты "РН Союз", сравнив результаты с телеметрией из видеоролика на Ютуб. Расчёты были произведены без учёта силы сопротивления атмосферы, что в итоге привело к существенным расхождениям с реальностью (если конечно же верить той самой телеметрии из видеоролика). Разумеется, мне стало интересно, а что если учесть это сопротивление? Как оно повлияет на траекторию и другие параметры полёта космической ракеты? Попробуем по порядку в этом разобраться.

Благодарю пользователя в комментариях, это дополнительно послужило мотивацией к данному исследованию, спасибо!

Если коротко, то моё исследование можно описать так - увяз коготок, вся птичка пропала. Хотелось обойтись какими-нибудь упрощёнными вычислениями, но, как уж получилось.

Для тех, кто не прочитал начало, оно тут https://habr.com/ru/post/649961/

Постановка задачи

Физическая модель, системы координат и допущения, принятые в предыдущей публикации остаются справедливыми и для текущих расчётов, за исключением сопротивления атмосферы. Напомню, ракета имеет три ступени. Соответственно, полёт разделяется на три этапа: полёт с момента старта до отстыковки первой ступени, с момента отстыковки первой ступени до момента отстыковки второй ступени, и с момента отстыковки второй ступени до момента отстыковки третьей ступени. Изменения в вычислениях коснутся только первого этапа полёта, то есть от старта до момента отстыковки первой ступени. На этом участке полёта ракета преодолевает наиболее плотные слои атмосферы и испытывает вместе с этим наибольшее сопротивление трения. Забегая вперёд, из вычислений получилось, что в конце работы первой ступени сила сопротивления атмосферы, действующая на ракету (высота 45 км, скорость 1700 м/с), составляет около 5 тонн-сил!

Напишем уравнение динамики с учётом силы сопротивления:

$m\overrightarrow{w}=\overrightarrow{P_T}+ \overrightarrow{G} + \overrightarrow{R}$

где m - масса ракеты, $\overrightarrow{w}$ - вектор ускорения, $\overrightarrow{P_T}$ - вектор силы тяги двигателей, $\overrightarrow{G}$ - вектор силы тяжести, $\overrightarrow{R}$ - сила сопротивления атмосферы.

Разделив обе части на массу ракеты и сделав необходимые подстановки (см. первую публикацию), получим:

$\overrightarrow{w}=\frac{\overrightarrow{P_T}+\overrightarrow{R}}{M_0-kt}+\overrightarrow{g}$

Аэродинамическое сопротивление

Теперь давайте разберёмся, что такое .

Аэродинамическое сопротивление вычисляется по формуле:

$R=C_x\frac{\rho v^2}{2}S$

где C_x - коэффициент лобового аэродинамического сопротивления, $\rho$ - плотность атмосферы, - скорость движения в среде, - характерная площадь.

Сначала разберёмся с плотностью атмосферы.

Как известно, плотность атмосферы вслед за давлением убывает с высотой. Но не всё так просто. Плотность атмосферы также зависит и от температуры, которая тоже убывает с высотой. Но и это ещё н�� всё. Мы собираемся лететь так высоко, что будем пересекать такие слои атмосферы, где температура не изменяется или даже возрастает.

Теперь в правильных терминах.

Введём параметр L_i - градиент температуры. Не надо пугаться, в нашем случае это просто положительное или отрицательное число, которое характеризует быстроту и направление изменения температуры в i - том слое атмосферы. Нумерация слоёв начинается с самого нижнего слоя - тропосферы. Если градиент отрицательный, то температура атмосферы убывает, если положительный - возрастает. Атмосфера Земли хорошо изучена и градиенты температуры слоёв измерены и известны. Вот они:

Номер слоя	Диапазон высот, км	Градиент температуры,	Температура в начале слоя , K	Давление в начале слоя , гПа
1	0 - 11	-6,5	288	1030
2	11 - 20	0,0	216	229,8
3	20 - 32	+1,0	216	55,3
4	32 - 47	+2,8	227	8,7
5	47 - 51	0,0	270	1,1
6	51 - 71	-2,8	270	0,6
7	71 - 85	-2,0	216	0,03

Таким образом, зная эту таблицу, можем вычислить температуру на любой высоте по следующей формуле:

где H_i - высота начала слоя (например, для первого слоя H_1=0 ).

Изобразим эту зависимость графически:

Давление для каждого слоя соответственно вычисляется по формулам:

если градиент температуры $\neq0$

$P(h)=P_i(\frac{T_i}{T_i+L_i(h-H_i)})^\frac{34,163}{L_i}$

если градиент температуры

$P(h)=P_i\exp(\frac{-34,163(h-H_i)}{T_i})$

Плотность есть некоторая функция температуры и давления f(P,T) , где температура и давление в свою очередь являются функциями высоты. Плотность вычисляется по формуле:

$\rho=\frac{PM}{RT}$

где - молярная масса воздуха, - универсальная газовая постоянная.

Зависимость плотности от высоты будет выглядеть следующим образом:

Итак, с плотностью воздуха разобрались. Теперь вернёмся к формуле аэродинамического сопротивления и посмотрим на ещё один интересный параметр - C_x - аэродинамический коэффициент сопротивления. Наш полёт происходит на разной высоте, с разной скоростью. Поэтому этот коэффициент так же как и плотность воздуха не может считаться константой. Если рассматривать большой диапазон скорости летательного аппарата, например от близких к нулю значений, до нескольких Махов, а это как раз наш случай, то окажется, что коэффициент C_x значительно изменяется, и мы не можем этого не учитывать. В данном случае этот коэффициент будет зависеть от числа Маха, то есть от скорости полёта ракеты. Число Маха, в свою очередь, зависит от скорости звука, а скорость звука зависит от температуры среды, в которой он распостраняется. А, как мы выяснили раньше, температура среды изменяется с высотой. Давайте попробуем это записать:

Разберём по порядку все зависимости. Для начала займёмся функцией C_x(M) - зависимостью коэффициента сопротивления от числа Маха. После продолжительных исследований литературы на эту тему я решил найти готовый, наиболее подходящий под задачу по��ёта ракеты вариант, нежели самому проводить расчёты этой зависимости. Коэффициент сопротивления сильно зависит от формы обтекаемого газом тела, его геометрических параметров, плюс отдельно считаются боковые блоки, элементы аэродинамики и т. д. Методика таких расчётов довольно объёмна и муторна, приводить её здесь я посчитал излишним. Поэтому привожу то, что нашёл уже посчитанным для реальной ракеты. Вот оно:

Интересно то, что при приближении скорости к числу Маха M=1 и его пересечении коэффициент сопротивления резко возрастает. Происходит так называемый скачок уплотнения. После этого при дальнейшем возрастании скорости коэффициент несколько уменьшается.

Со следующей зависимостью - M(v) всё просто: число Маха есть отношение скорости движения в среде к скорости звука - $M(v)=\frac{v}{v_s}$ .

Зависимость скорости звука в воздухе от температуры v_s(T) тоже известна, её можно найти в любом справочнике по физике:

Теперь напишем формулу для вычисления силы сопротивления воздуха с учётом всех выше приведённых расчётов:

$R=C_x(M(v, v_s(T(h))))\frac{\rho(h)v^2}{2}S$

Подставим силу сопротивления в основное уравнение динамики и распишем его на оси координат:

$w_x=\frac{(P_T(t) - C_x(Mach(v,v_s(T(h))))\frac{\rho(h)v^2}{2}S)sin(\alpha(t))}{M(t)}$ $w_y=\frac{(P_T(t)-C_x(Mach(v, v_s(T(h))))\frac{\rho(h)v^2}{2}S)cos(\alpha(t))}{M(t)}-g$

Или в производных:

$\ddot{x}=\frac{(P_T(t)-C_x(Mach(\sqrt{\dot{x}^2+\dot{y}^2},v_s(T(y))))\frac{\rho(y)(\dot{x}^2+\dot{y}^2)}{2}S)sin(\alpha(t))}{M(t)}$ $\ddot{y}=\frac{(P_T(t)-C_x(Mach(\sqrt{\dot{x}^2+\dot{y}^2},v_s(T(y))))\frac{\rho(y)(\dot{x}^2+\dot{y}^2)}{2}S)cos(\alpha(t))}{M(t)}-g$

Таким образом, задача сводится к решению системы дифференциальных уравнений вида:

$\begin{cases}&\ddot{x}=f_1(t,\dot{x},y,\dot{y})\\&\ddot{y}=f_2(t,\dot{x},y,\dot{y})\end{cases}$

что мы и сделаем численным методом с помощью программы.

Входные данные

Параметры ракеты (в основном, массовые характеристики) были уточнены. Поэтому есть расхождения между траекториями без учёта сопротивления воздуха в предыдущей публикации и в этой. В данный момент они являются более точными. Ссылки на источники будут ниже.

Результаты вычислений

Результаты весьма интересны. Честно говоря, они меня впечатлили. Я не думал что атмосфера настолько сильно влияет на траекторию полёта и конечные орбитальные параметры. Разницу траекторий без учёта силы сопротивления и с учётом этой силы Вы можете видеть на этом изображении:

Давайте сравним полученные данные.

В момент перед отстыковкой первой ступени:

	Параметры телеметрии	Расчёты программы с учётом R	Расчёты программы без учёта R
Высота, км	45	44	51
Дальность, км	48	47	51
Скорость, км/ч	6312	6198	6785
Перегрузка, g	4	3,95	3,99

В момент перед отстыковкой второй ступени:

	Параметры телеметрии	Расчёты программы с учётом R	Расчёты программы без учёта R
Высота, км	154	153	185
Дальность, км	452	459	480
Скорость, км/ч	13732	13864	14266
Перегрузка, g	2,3	2,3	2,3

В момент перед отстыковкой третьей ступени:

	Параметры телеметрии	Расчёты программы с учётом R	Расчёты программы без учёта R
Высота, км	202	204	281
Дальность, км	1675	1725	1770
Скорость, км/ч	26737	27120	27386
Перегрузка, g	2,9	2,8	2,8

Хотел бы привести ещё один график, который мы немного проанализируем:

Это зависимость сопротивления атмосферы от высоты.

Ну во-первых, сразу бросается в глаза значение максимума - 740 кН, это 75 тонн-сил! Да, уже на высоте чуть больше 10 км ракета набирает такую скорость, что сила сопротивления воздуха составляет такую большую величину, даже с учётом того, что атмосфера на этой высоте значительно разреженная. Для сравнения, когда ракета стартует, избыток тяги (разница между тягой двигателей и весом ракеты) составляет 1130 кН. То есть сила сопротивления на максимуме составляет две трети от тяги на старте!

Также интересно, насколько быстро нарастает сила сопротивления, но это и не удивительно. Ракета - тело переменной массы. Ракета теряет массу, ускорение стремительно возрастает. Эффекта добавляет здесь ещё тот факт, что двигатели существенно прибавляют мощности с ростом высоты (тяга в ваккууме больше, чем на уровне моря).

Ещё один интересный результат - сопротивление атмосферы в момент отстыковки первой ступени. Казалось бы, высота уже 45 км, атмосфера крайне разреженная. Но не тут то было, получите: 46 кН (4,7 тонн-сил)! Неожиданно, правда? Но если учесть, что в этот момент ракета летит со скоростью 1722 м/с, что уже является даже не сверхзвуковой, а гиперзвуковой скоростью (> 5 Маха), то можно в это поверить. К тому же если сравнить с тягой двигателя в этот момент, а осталась у нас только вторая ступень, вполне приемлемо:

$\frac{990200кН-46000кН}{990200кН}=0.95$

95% тяги остаётся, потери на сопротивление всего 5%, и оно продолжает уменьшаться, мы же взлетаем.

Ну и в завершении обратим внимание на то место, которое обозначено красным овалом. Там явно прослеживается излом. Давайте посмотрим, что в этот момент происходит со скоростью:

Зелёный график - зависимость скорости от высоты, чёрный - зависимость силы сопротивления от высоты. Ось абсцисс выдержана в одном масштабе. а ось ординат теперь это значение скорости. Отсюда видно, что в момент излома скорость составляет почти 400 м/с. Что это за скорость? Вычислим число Маха для данной высоты. На высоте, соответствующей излому (~8 км) скорость звука составляет примерно 308 м/с

$M=\frac{400 м/с}{308м/с}=1,3$

Теперь обратимся к графику зависимости аэродинамического коэффициента сопротивления от числа Маха:

Данное зачение числа Маха соответствует резкому прекращению возрастания коэффициента сопротивления. Физически это означает, что ракета в данный момент закончила преодолевать трансзвуковой барьер (0,8 < M < 1,2).

На этом всё, спасибо за внимание!

Ссылка на программу здесь, бранч soyz

Использованные источники:

Параметры РН Союз

Аэродинамический коэффициент сопротивления

Параметры атмосферы

Зависимость скорости звука от температуры: справочник по физике.