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

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

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

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

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

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

Физическая модель, системы координат и допущения, принятые в предыдущей публикации остаются справедливыми и для текущих расчётов, за исключением сопротивления атмосферы. Напомню, ракета имеет три ступени. Соответственно, полёт разделяется на три этапа: полёт с момента старта до отстыковки первой ступени, с момента отстыковки первой ступени до момента отстыковки второй ступени, и с момента отстыковки второй ступени до момента отстыковки третьей ступени. Изменения в вычислениях коснутся только первого этапа полёта, то есть от старта до момента отстыковки первой ступени. На этом участке полёта ракета преодолевает наиболее плотные слои атмосферы и испытывает вместе с этим наибольшее сопротивление трения. Забегая вперёд, из вычислений получилось, что в конце работы первой ступени сила сопротивления атмосферы, действующая на ракету (высота 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

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

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

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

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

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

ссылка на оригинал статьи https://habr.com/ru/post/652407/