Радиус падения ступеней ракет

ПОВАЛЯЕВ П.П., Московский авиационный институт (национально исследовательский университет), инженер

КАЛЁНОВА Н.В., Московский авиационный институт (национально исследовательский университет), доцент

РОМАНЕНКОВ А.М., Московский авиационный институт (национально исследовательский университет), доцент

ШАРИКОВ Д.В., Московский авиационный институт (национально исследовательский университет), инженер

Движения ступеней ракет-носителей
для уменьшения зон «отчуждения»

В данной работе, рассмотрен подход, который позволит исследовать движение ступеней ракет-носителей. Для примера был рассмотрен запуск ракеты-носителя «Союз», а точнее критическое падение второй ступени. Чтобы выработать рекомендации для уменьшения зоны «отчуждения», реализована методика численного моделирования свободных колебаний модели. Была доработана базовая программа, основанная на методе Годунова – Родионова интегрирования нестационарных уравнений Эйлера, путем дополнения блоком интегрирования обыкновенного дифференциального уравнения плоского движения тела (включающий в себя плоское угловое движение и движение центра масс) с учетом изменений параметров набегающего потока за счет изменения высоты. Для решения газодинамической задачи, разработана программа построения трехмерной структурированной вычислительной сетки, в прямоугольной системе координат. Обработаны полученные данные и визуализированы в специализированных программных пакетах такие как Tecplot и Grapher 2.0

При запуске ракеты-носителя образуются зоны падения отработанных ступеней и их обломков в различных районах площадью в сотни квадратных километров. Это приводит к тому, что возникают зоны «отчуждения», в которых невозможна хозяйственная деятельность и размещение населенных пунктов.

Рассмотрим запуск ракеты-носителя «Союз».

• Первая ступень совершает отделение в непосредственной близости к месту старта ракеты-носителя, где отсутствуют населенные территории.
• Вторая ступень, при падении с огромной скоростью разрушается в плотных слоях атмосферы, и обломками разносится на большие территории, тем самым создавая зону «отчуждения».

Задача данной работы – рассмотреть подход, который позволит исследовать движение ступеней ракет-носителей с целью выработок рекомендаций для уменьшения зоны «отчуждения».

Известны способы расчетов динамики падения ступеней ракет-носителей, в ходе которых заранее формировались массивы аэродинамических коэффициентов ступеней, в зависимости от числа Маха, высоты, углов атаки и набегающего потока. Далее с использованием интерполированных данных этих массивов решалась задача динамики.

Нужно реализовать методику комплексного расчета, решить одновременно уравнение газовой динамики и уравнение движения центра масс, а также и углового движения относительно центра масс.

Для интегрирования уравнения движения будет использован простейший метод Эйлера. Поскольку идет разработка методики, для ускорения времени счета используется простейший метод Эйлера, который позволяет интегрировать уравнение динамики, используя неравномерный шаг по времени.

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

Исследование движения газообразной среды, т.е. определение в каждой точке пространства параметров, характеризующих это движение, состоит в решении соответствующих уравнений, которые связывают между собой эти параметры. Все эти уравнения независимы и составляют систему уравнений газодинамики. Число независимых уравнений системы определяем количеством отыскиваемых независимых параметров газа.

Рассмотрим движения идеального газа. Если скорости потока невелики, то можно пренебречь изменением удельных теплоемкостей от температуры и не учитывать излучения. В этом случае газовое течение представляет собой термодинамическую изолированную систему и является адиабатическим. Система уравнений газодинамики должна включать шесть независимых уравнений.

Запишем систему уравнений по осям прямоугольной системы координат X , Y , Z , и проекции вектора скорости u , v , w на данные оси.

Система уравнений Эйлера, записанная в дивергентной форме, имеет вид:

(1)

В ряде случаев для определения характера движения системы требуется знать закон движения ее центра масс. Чтобы найти этот закон просуммируем левую и правую часть системы дифференциальных уравнений движения, система состоит из n материальных точек, точка системы с массой m , равнодействующая всех приложенных к точке внешних сил , и равнодействующая всех внутренних сил ,точка имеет ускорение a_n .

Преобразуем левую часть равенства. Из формулы (2) для радиуса-вектора центра масс получаем

Беря от обеих частей этого равенства вторую производную по времени и замечая, что производная от суммы равна сумме производных, найдем

Так как по свойству внутренних сил системы , получим окончательно из равенства (2) учитывая (4)

Уравнение (5) и выражает теорему о движении центра масс системы. Проектируя обе части равенства на координатные оси, получим

Эти уравнения представляют собой дифференциальные уравнения движения центра масс в проекции на оси декартовой системы координат.

Одним из наиболее эффективных методов численного решения уравнений газодинамики является метод С.К. Годунова. Годунов С.К. в 1957 году предложил численный метод решения уравнения газовой динамики, который основывается на использовании точного или приближенного решения задачи Римана.

Рассмотрим равномерную пространственную сетку с шагом ∆х . Значения сеточной функции будем обозначать, как U k _i и U k _i±1/2 . Нижний целый индекс i=1,2,… обозначает значения функции, отнесенные к центру масс i -й дискретной ячейке. Полуцелые нижние индексы i±1/2 обозначают значения сеточной функции на границе между ячейками с номерами i и i±1 . Верхний целый индекс k=0,1,2. обозначает номер слоя (шага) по времени.

Положим, что все сеточные функции являются постоянными внутри каждой из дискретных ячеек. Тогда на границе с номером i±1/2 на каждом шаге по времени будем решать задачу Римана со следующими начальными данными:

Пусть U k _i+1/2 – это решение такой задачи в точке x=x_i+1/2 . Аналогично, U k _i–1/2 – это решение такой задачи в точке x=x_i–1/2 . Тогда явная конечно-объемная схема Годунова имеет вид:

Эта схема обладает первым порядком точности по времени и по пространству.

Схему (7) можно обобщить на неравномерные пространственные и подвижные сетки. Она также может применяться в расчетах с выделением разрывов. Рассмотрим, например, двумерную гиперболическую систему уравнений:

Конечно-объемная схема для этой системы на равномерной декартовой сетке записывается в виде:

Таким образом, для построения схемы следует решать задачу Римана на каждой из границ дискретных ячеек.

Схему Годунова (9) можно обобщить на произвольную пространственную сетку. Перепишем систему (10) в интегральной форме:

Построим явную конечно-объемную схему Годунова для уравнений в интегральной форме (13). Покроем для этого вычислительную область дискретной сеткой, состоящей из произвольных выпуклых конечных многоугольников с площадями G_i , где i=1,2,… , и с числом сторон m = m ( i ), каждая из которых имеет длину S_j и внешнюю нормаль n_j , где j=1,2,…,m(i) . Аппроксимируем интегральные уравнения на каждом из многоугольников следующим образом:

Явная конечно-объемная схема Годунова на произвольной движущейся сетке записывается в виде:

Уравнение (16) является уравнением изменения объема дискретной ячейки G_i . Оно играет роль дискретного условия совместности аппроксимаций для G и S . При заданной аппроксимации площадей S уравнение (16) можно использовать для определения объемов G_i (k+1) . В ряде случаев это дает поправку к величине объема, вычисленного исходя только из геометрических соображений.

Выполнения уравнения (16) можно также добиться выбором подходящих величин D_j . Важнейшее свойство аппроксимаций G и S , которые удовлетворят соотношению (16), заключается в следующем. Решение системы уравнений (8) в виде состояния U=U₀=const является также решением дискретных уравнений (15) и (16). Соблюдения этого условия следует добиваться при использовании движущихся и криволинейных систем координат.

Схема (14) устойчива на равномерной прямоугольной сетке при условии:

где C_x и C_y являются числами Куранта, соответственно, для осей х и у и определяются ∂F/∂U и ∂E/∂U . Аналогично строится, конечно-объемный метод Годунова для многомерного случая.

Рассмотрим метод Эйлера рассмотрим на примере простейшего дифференциального уравнения y’ = f ( x , y ) с начальным условием y(x₀)=y₀ . Подставив x₀ , y₀ , в уравнение, получим значение производной в точке x₀ :

При малом ∆x₀ имеет место:

Обозначим f(x₀,y₀)=f₀ , перепишем последнее равенство в виде:

Принимая x₁ , y₁ теперь за новую исходную точку, точно также получим:

В общем случае будем иметь:

Рассмотрим пример задачи.

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

Предполагается, что с высоты 80 км начинается интенсивное воздействие атмосферы на ступень ракеты-носителя со скорость ϑ₀ и углом атаки α . Дальше движущая система координат двигается со скорость ϑ₀ равномерно.

Строим местную систему координат, фиксируем под каким углом θ по оси О_х она войдет в атмосферу. Связываем первый момент времени с центром масс и строим систему координат x , y , так что вектор скорости ϑ направлен по оси О_х . Интегрируем систему уравнений в движущейся системе координат.

Система координат центра масс

Положение в земной системе координат

Высота будет равна , по которой определяем параметры атмосферы на данной высоте.

Для интегрирования уравнения движения используем простейший метод Эйлера. При этом выполняется шаг интегрирования уравнения газовой динамики, а затем на этом же временном шаге производится расчет движения ступени. Определяется пространственное и угловое движение.

Расчетная область представляет собой сферу, в центре располагается исследуемое тело. Сфера выбрала с точки зрения, чтобы шел нормальный расчет, и формировалась поверхность ударной волны. Вокруг ступени выделяем вторую подобласть, которая со сферой образует единую область. Такая расчетная область позволит считать возле поверхности тела более точные данные, так как сетка более мелкая, и сгущена к поверхности ступени ракеты-носителя.

Для численной реализации программы используется модификация метода С.К. Годунова, предложенная в работах. Для сокращения выкладок его описание приводится на примере расчета двумерных течений. Обобщение метода на трехмерный случай, очевидно.

Численное решение системы

(20)

осуществляется методом установления, при последовательном переходе от одного слоя t = const к другому.

Предположим, что задан алгоритм построения регулярной четырехугольной конечно-разностной сетки на каждом слое t=t l ( l=1,2,3,… ). Обозначим координаты узлов разностной сетки через ( j=1. J1 ; k=1. K1 ), а осредненные по элементарным плоским ячейкам параметры через ω l _j–1/2 , k–1/2 ( ω = u , v , w , ρ , P , γ , h₀ ). Проинтегрируем систему дифференциальных уравнений (31) по одной из элементарных ячеек от слоя х = x l до х = х l+1 . В качестве примера выпишем конечно-разностный аналог уравнения:

Здесь ∆S – площадь элементарной ячейки в плоскости t = const , ∆t=t l+1 –t l ∙∏ – поток энтальпии через боковую грань ячейки, определяемый соотношением:

Где и – осредненные параметры на боковых гранях ячейки:

– площадь проекции четырехугольника с вершинами на плоскость х = const .

– площадь проекции четырехугольника с вершинами на плоскость х = const .

Предположим, что известны осредненные параметры на l -ом слое и на четырех боковых гранях ячейки , , и . Тогда разностные уравнения переписываются относительно неизвестных параметров :

Уравнения сразу позволяют определить величины:

Для совершенного газа оставшиеся величины определяются по соотношениям:

Таким образом, определена процедура расчета параметров : по известным параметрам на l -том слое и на четырех боковых гранях ячейки – , , и .

Прежде, чем перейти к описанию полной процедуры перехода от слоя t = t l к слою t = t l+1 , сформулируем одно из основных предположений.

Будем полагать, что на каждом l -оме распределение параметров ω является кусочно-линейным с разрывами на линиях конечно-разностной сетки.

Для задания кусочно-линейного распределения кроме параметров, приписываемых центрам ячеек ω l _{j–1/2,k–1/2} , будем использовать их приращения по двум направлениям:

Итак, полная процедура расчета параметров по известным ω l _{j–1/2,k–1/2} ( переход от l -го слоя к l+1 -му) состоит из следующих этапов.

Этап I. Для каждой ячейки определяются приращения параметров [∆jῶ] l _{j–1/2,k–1/2} и [∆kῶ] l _{j–1/2,k–1/2} по известным параметрам в центрах четырех соседних ячеек и их координатам. Затем эти приращения корректируются для обеспечения монотонности кусочно-линейного распределения параметров:

Здесь функция minmod определяется следующим образом:

Этап II.(Шаг предиктор). Предварительно рассчитываются параметры по описанной выше процедуре, в которой:

Этап III. Определяются параметры на боковых гранях ячеек из решения задачи Римана (распад разрыва) с начальными параметрами

для определения и

Этап IV (Шаг корректор). Окончательный расчет .

Описанная разностная схема имеет второй порядок аппроксимации в областях в равномерной конечно-разностной сеткой и первый порядок – в местах сильной неравномерности сетки. Она обладает свойством монотонности при выполнении условий

где – локально максимально допустимый расчетный шаг по условию Куранта-Фридрихса-Леви.

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

Начальные условия определяются значениями параметров газа для некоторого момента времени и имеют смысл, очевидно, для неустановившегося движения. M_∞=10 , α=30° .

Граничное условие определяется характером течения газа на обтекаемой поверхности. Если газ невязкий и не проникает сквозь такую поверхность, то течение характеризуется условием безотрывного обтекания.

В соответствии с этим условием в каждой точке поверхности составляющая скорости, нормальная к ней, равна нулю, а вектор полно скорости совпадает с направлением касательной к поверхности.

При неустановившемся движении газа для некоторого момента t должны быть заданы значения всех искомых параметров движения v_x , v_y , v_z , ρ , р , Т . Если движение газа установившееся, то начальные условия отпадают.

Граничные условия идеального газа. Если течение ограничено неподвижной границей, сквозь которую газ не проникает, то условием безотрывного обтекания границы газом будет равенство нулю нормальной составляющей скорости газа к поверхности границы. Если F(x,y,z)=0 – уравнение поверхности границы, то это условие запишется так:

Если же граница является подвижной поверхностью F(x,y,z,t)=0 и течение газа безотрывное, то нормальная составляющая скорости частицы в любой точке поверхности должна равняться нормальной скорости D точки самой поверхности. (20), выражается формулой:

нормальная составляющая скорости к стенке

Из равенства этих величин (поскольку G≠0) получаем:

Вычислительная сетка внутри расчетной области строится с помощью простого алгебраического генератора.

Программа для интегрирования систем уравнения, а также построение вычислительной сетки была написана на языке FORTRAN. Рассмотрим блок-схему программы (см. рис. 1).

Рисунок 1. Блок-схема программы расчета

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

Далее следуя по смоделированной блок-схеме, идет вызов подпрограммы PRC. PRC является основной подпрограммой в расчетном модуле, а для ввода данных и запись результатов используются такие подпрограммы, как PRCINP и PRCOUT.

Все три подпрограммы отвечают за ввод исходных данных и их преобразование к требуемому виду, построение конечно-разностной сетки, задание граничных условий и последовательное обращение к программам пошагового интегрирования уравнений газовой динамики, обработку и вывод результатов расчета. Решаемая задача – трехмерное обтекание осесимметричного тела совершенным газом.

PRС является основной в модуле.

PRСINP определяет начальные данные.

PRСOUT выводит результаты расчета.

Файлы входных данных имеют фиксированную структуру, т.е. не допускается перестановка строчек и переменных в файле. Границы задания параметров ограничены угловыми скобками. Формат ввода целых и действительных чисел внутри скобок произвольный.

Координаты узлов исходной вычислительной сетки должны быть записаны в виде:

Пример заполнения файла приведен ниже:

Выходными данными программы являются файлы двух типов:

• графические файлы;
• файлы со значениями основных аэродинамических характеристик и координат узлов сетки.

Здесь значения компонент скорости отнесены к модулю скорости набегающего потока, плотность – к плотности набегающего потока, давления – к удвоенному скоростному напору набегающего потока.

Пример содержания файла для визуализации PPROD:

Рассмотрим результаты программы.

Наши исходные данные для вычислений: α=30° , M=10 , начальная высота 80 км .

На рис. 2-3 в качестве иллюстративного материала, представлены области распределения давления вблизи ступени при различных положениях относительно системы координат.

Рисунок 2. Распределение давление в расчетной области при числе Маха набегающего потока равное 10, при δ =5° и δ =55°

Рисунок 3. Распределение давление в расчетной области при числе Маха набегающего потока равное 10, при δ =120° и δ =180°

Результаты обработки t , δ , ρ , P , h , x_c , y_c представлены на рис. 4-6 в виде графиков зависимостей угла атаки от времени при колебании ступени при различных положениях ступени ракеты-носителя относительно системы координат (угла δ ). Колебания несут затухающий характер.

Рисунок 4. Зависимости δ от времени для угла δ =5° и δ =55°

Рисунок 5. Зависимости δ от времени для угла δ =120° δ =180°

Рисунок 6. Зависимости угла отклонения δ от времени при колебаниях для различных углов δ

На рис. 7 ступень двигается относительно движущейся системы координат. Система координат движется с постоянной скоростью, а тело в этой системе обувается потоком. По оси Х видно, на какое расстояние относительно начало координат отдалилась ступень.

Рисунок 7. Движение относительно движущейся системы координат на высоте h =40 км, t =26.4сек и h =31,7 км, t =31сек

Ниже на рис. 8-9 показаны зависимости движение ступени по координате X_c от времени. Можно наблюдать как отклоняется ступень во время падения с высоты 80 км от начальной точки.

Рисунок 8. График движения зависимость X_c от t при начальном угле отклонение δ =5°и δ =55°

Рисунок 9. График движения зависимость X_c от t при начальном угле отклонение δ =120° и δ =180°

2-ая ступень в отличие от первой успевает за свое время работы набрать достаточно скорости и высоты, чтобы начать при падении разваливаться на десятки объектов. C помощью полученных результатов можно оценить размеры территории падения обломков ступени. Зона рассеивания обломков зависит от движения ступени и ее положении (угла отклонения δ ).

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

Рассматривать ступень под разными углами ориентации ( δ ). Рассчитывать такие параметры, как t , δ , ρ , P , h , x_c , y_c с помощью которых, можно определить координаты падения и построить область разброса обломков.

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

1. Годунов С.К. и др. Численное решение многомерных задач газовой динамики. М.:»Hаука».1976.
2. Родионов А.В. Численный метод решения уравнений Эйлера с сохранением аппроксимации на деформируемой сетке. ЖВМ и МФ. 1996. 36. №3. сс.117-129.
3. Рахматулин Х. А., Сагомонян А. Я., Бунимович А.И., Зверев И. Н. Газовая динамика. Изд. «Высшая школа» 1965.
4. J. Blazek. Computational Fluid Dynamics: Principles and Applications. Elsevier. 2005.
5. Г. В. Липман, А. Рошко. Элементы газовой динамики. Перевод с английского В. П. Шидловского. Издательство иностранной литературы. Москва 1960.-505с.
6. Липницкий Ю.М., Красильников А.В., Покровский А.Н., Шманенков В.Н. «Нестационарная аэродинамика баллистического полета», Москва, ФИЗМАТЛИТ, 2003.
7. Дэниел Д. Мак-Кракен, Уильям С. Дорн. Численные методы и программирование на ФОРТРАНе. Ориг. название – Numerical Methods and Fortran Programming with Application in Engineering and Science, 1977.

Ключевые слова: Числа Маха, метода Эйлера, уравнение динамики, движение идеального газа, задачи Римана, схема Годунова, аппроксимация, Fortran, трехмерное обтекание, зона рассеивания.

Movement of launch vehicle stages for open alienation zones

Povalyayev P.P., Moscow Aviation Institute (National Research University), engineer

Kalenova N.V., Moscow Aviation Institute (National Research University), Associate Professor

Romanenkov A.M., Moscow Aviation Institute (National Research University), Associate Professor

Sharikov D.V., Moscow Aviation Institute (National Research University), engineer

Abstract: In this paper, we consider an approach that will allow us to investigate the movement of the stages of launch vehicles. For example, the launch of the Soyuz launch vehicle was considered, or rather the critical fall of the second stage. To develop recommendations for reducing the “alienation” zone, a method of numerical modeling of free oscillations of the model has been implemented. The basic program based on the Godunov – Rodionov method of integrating the non-stationary Euler equations was improved by adding the ordinary differential equation of plane body motion (including plane angular motion and center-of-mass motion) to the integration unit, taking into account changes in the parameters of the incident flow due to height changes. To solve the gas-dynamic problem, a program for constructing a three-dimensional structured computational grid, in a rectangular coordinate system, has been developed. Processed data and visualized in specialized software packages such as Tecplot and Grapher 2.0.

Keywords: Mach numbers, Euler method, dynamic equation, ideal gas motion, Riemann problem, Godunov scheme, approximation, Fortran, three-dimensional flow, scattering zone.

Источник