Численное моделирование трехмерных газодинамических процессов в камерах сгорания РДТТ на основе метода геометрически-адаптивных сеток

ПРИКЛАДНАЯ МАТЕМАТИКА И МЕТОДЫ МАТЕМАТИЧЕСКОГО МОДЕЛИРОВАНИЯ

УДК 534.222.2

Ю. И. Димитриенко, С. Г. Изотова, С. Н. Ануфриев, А. А. Захаров

ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ ТРЕХМЕРНЫХ ГАЗОДИНАМИЧЕСКИХ ПРОЦЕССОВ В КАМЕРАХ СГОРАНИЯ РДТТ Ш. ОСНОВЕ МЕТОДА ГЕОМЕТРИЧЕСКИ-АДАПТИВНЫХ СЕТОК

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

Большинство существующих методов решения задач нестационарной газовой динамики основано на использовании так называемых регулярных прямоугольных сеток [1-4]. Однако для задач обтекания газовыми потоками поверхностей со сло ной геометрической формой эти методы становятся неэффективными, поскольку требуют чрезмерного сгущения сеток и приводят к появлению нефизических осцилля-ций, связанных с аппроксимацией границы регулярными участками. Более перспективными являются методы адаптивных сеток [5-10], в которых ячейки сетки становятся криволинейными, "подстраиваясь" под форму границы рассматриваемой области. Алгоритмы введения таких адаптивных сеток могут быть различными, но все они основаны на использовании специальных систем координат. В настоящей работе предложен такой алгоритм для решения трехмерных нестационарных задач течения продуктов сгорания в соплах.

Математическая формулировка задачи. Рассмотрим сопловой тракт ракетного двигателя на твердом топливе (РДТТ), состоящий из камеры сгорания и соплового блока, разделенных в начальный момент мембраной, расположенной в критическом сечении (рис. 1); мембрана разрушается после того, как давление в камере сгорания достигает определенного значения.

Задача течения продуктов сгорания в сопловом тракте Р ТТ (газ предполагается линейно-вязким, теплопроводным) имеет вид

+ V-рду9 = 0,

+ V •

-РдЕ + ^ (рд (Ед + Рд ).д) = ЛАвд + Ы

при следящих граничных условиях:

а) на поверхности горения топлива

V • п = -1(рз - р)Б, у • Т1 = 0, I = 1, 2, в = в*; (2) Р д

б) на остальных твердых поверхностях

V =0, п • ЛVвg = 0;

в) в выходном сечении рд = рде; и начальных условиях при £ = 0:

в = во, р = ро, V = 0.

Здесь рд — плотность газа; уд — вектор скорости; р — коэффициент

V2

вязкости; Ед — полная энергия, Ед = Су вд +--, Су — теплоемкость

2

при постоянном объеме, вд — температура газа, у2 — квадрат модуля скорости, определяемый следующим образом: у2 = у. + у. + у.;

Рис. 1. Расчетная схема задачи течения продуктов сгорания в РДТТ:

1 — твердотопливный заряд; 2 — корпус; 3 — сопло; 4 — мембрана; 5 — камера сгорания

p — давление, определяемое по формуле p = pgR Qg; pge — значение плотности в окружающем РДТТ газе; E — метрический тензор; n — вектор нормали; 77, I = 1, 2, — векторы в плоскости, касательной к поверхности горения.

Задача рассматривается в цилиндрической системе координат. Скорость горения топлива определяется по формуле

D = Ч РР У,

где D0, p0, v — константы, а функция ш определяется выражением

+ Г~дв^Тгв Vr + Tze Vz + Твв v^ + dz (Trz Vr + Tzz Vz + Tzd v^;

где Trr, Tzz, Tee, Trz, Tze, Tre — компоненты тензора вязких напряжений:

^ dvr dvz (1 dvr ve \

Trr = 2^—, Tzz = 2^——, Tee = 2^ —— +--,

or oz \r дв r J

T = (d^z + dM T = (1 d^r + r—(vi

rz \dr dz / r \r дв dr\r

T = (d^e + 1 dM

z \dz r дв J

Метод геометрически-адаптивных сеток для решения задачи (1)-(4). Введем новую систему криволинейных координат, связанную с поверхность соплового тракта:

r' = k(t, z)r, Z = z, r" = k*(t, z, U)r, в' = в,

_ r" Ok* _ dr'' _ dr''

x = k*~дt, v = ^ =

где

R0 R0 k(t, z) = fR V k*(t, z, в) = ff. w flV

f

f (t, z), g(z, в) — функции, описывающие форму поверхности рассматриваемого соплового тракта, для которой в цилиндрических координатах задается уравнение r = f (t, z)g(z^), а R0 — некоторый приведенный радиус (выбирается совпадающим с радиусом критического сечения).

В координатах (t, r", z', 9') система (1) принимает вид

dU dU dW (U) dV (U) dV (U) ~öt + Vd?> + dr'' + dz' +X dr'' +

+

k* dS(U) k* dS(U) k

9'

+ С

r'

+ -W

где координатные столбцы имеют следующий вид:

V =

Рд

Рд vgr

Рд VgO Рд vgz , Рд Eg >

W (U) =

Р<д vgr Рд v2r + Р

Рд vgr Ще

Р д V д r V д Z

Р<1°д r

Eg + Р

Р

V (U ) =

Р V Z

Р V ZV r

Р V ZV е 2

Рдv2z + Р

Р V Z

Eg + Р

Р

S (U ) =

Р g v g е

i)gvg е vg r

Р gv2e + Р Рg v g еvgz

Рgvgе ( Eg + ^^ Р

V (U ) =

г 0

pAv(g r pAvg е

pAvg z

ХАв

, H (U ) =

В координатах (т'', г', в') рассматриваемая область интегрирования представляет собой прямоугольный параллелепипед, грани которого в точности соответству т границам рассматриваемой области в исходных координатах (т, г, в). Введем в координатах (т'', г', в') разностную сетку, а алгоритм решения системы (6) разделим на три этапа.

На этапе 1 рассматривается система уравнений (6) с нулевой правой частью и применяется разностная схема типа Мак-Кормака, состоящая из двух шагов.

Шаг 1.1 (предиктор). Вычисляем промежуточное значение следу-ей функции:

Um+l/6 = U-rn m А

Uijk = Uijk Vijk Ar„

; Ar'

*

0

At At j At

^^^ ^ / т -riyyi т -riyyi \ m ^^^ ^ / т -riyyi т -riyyi \ j k ^^^

- (~\fm _\im\_ш /т^ш _\fm\ J'" (Qm _Qm \_

Az! ( Vi(3+1)k - V ijk ) -Xijk Ar// ( V(i+1)jk - Vijk )- r„ AQ' - Sijk )-

At k*m k*m

At jk j аш пш \ jk Tirm

Ar'' r'' (S(i+1)3k Sijk) At rn Wijk.

i i

Шаг 1.2 (корректор). Вычисляем промежуточное значение следующей функции:

ит+1/3 _ 1 (ттт+1/6 + ит \ _ пт А (пт+1/6 _ ит+1/6 N

и13к _ 2 \ичк + иг]к) Чг3к ^АТ'' гзк и

_ к*т А (ЛЛТ т+1/6 _ ш т+1/6 ч А т+1/6 _ л/т+1/6\_

кзк 2Аг"

_ т А (ут+1/6 _ Vт+1/6 \ _ к,]к А / ет+1/6 _ ет+1/6

Хгзк 2Аг" ^Угзк у

_ ст АТ кзк /ет+1/6 _ ет+1/6 \ кзк Аг шт+1/6

З 2АГ" Г" ^^гзк е

Здесь итк, ит+1/6, Щ+1/3 — значения функций на временных шагах т, т + 1/6, т + 1/3 соответственно в узле (г,], к) разностной сетки (индекс г соответствует координате т", индекс ] — координате г', а индекс к — координате 9'). Для коэффициентов х, применялась следующая разностная аппроксимация:

Vijk

¿ш _

^ijk

r'' 1 г kjk kjk

jk 2At '

r'' 1 i k(j+1)k k(j- 1)k

jk 2 A '

r'' i kj

jk 2AQ'

хш _ ' г

Xijk _ Ь,*ш о Л ' ^ '

На этапе 2 для учета вязкости и теплопроводности в системе (6) применяется метод расщепления по координатам, согласно которому рассматривается следующая система уравнений:

рд_ ?

где ит = (удГ,ь9в,Удг, СуОд), а координатные столбцы Г(й),Н(и) определяются по формуле (7), но рассматриваются как функции от и.

Для решения системы уравнений применяется четырехшаговая схема расщепления.

Шаг 2.1. рт+1/3

д 2

Шаг 2.2. рт+1/3

д2

Шаг 2.3. рт+1/3

Шаг 2.4. рт+1/3(ит+1 - ит+1/3) =

= Аг

Здесь оператор й(и) представлен в виде суммы трех координатных (прогоночных) операторов:

й

где Л — линейный дифференциальный оператор второго порядка по соответствующей координате. Разностные системы (13) - (15) решаются методами прогонки относительно ит+1/2, ит+2/3 и ит+5/6соответственно, а система (16) — по явной схеме относительно ит+1. Для учета функции ш вводится оде один расчетный шаг. Шаг 2.5. итк+1 = итт++1 + АЬНт+\ где Щ^1 — искомое решение, а Щт+1 — решение, полученное на (т + 1)-м временном шаге схемы расщепления.

На этапе 3 для подавления нефизических осцилляций применяется оператор искусственной вязкости

итт?1 = итт?1+аьз)+ъьз)+азьз (з1),

где

Т (ТТт+1) _ ттт+1 _ 2ТТт+1 + ттт+1

Ь1(иг]к ) = и(г+1)3к 2иЦк + и(г-1) з к, т+1\ _ ттт+1 г\ттт+1 , ттт+1

Г (ТТгп+1\ _ ттт+1 _ Г) ттт+1 , т тт+1

L2(Uijk ) = Ui(j+1)k 2Uijk + Ui(j-l)k

T (Tjm+L\ _ 7-7-m+l _ 9 7-7-m+l 1 т тт+1

L3

— операторы сглаживания, а aL, a2, — параметры сглаживания, зависящие от Um+l.

Предложенная разностная схема с наличием оператора искусственной вязкости, как известно, приводит к искажению решения задачи. Для устранения этого эффекта параметры сглаживания выбирались отличными от нуля только в окрестности ударных волн. При численной реализации ирина этой окрестности составляла 4-5 разностных ячеек. Нахождение окрестностей ударных волн осуществлялось автоматически, с помощью отношения вторых и первых производных рассматриваемых функций итк+1. Для проверки точности решения было проведено тестирование разностной схемы на классической одномерной задаче о распаде произвольного разрыва, которое осуществлялось на продолжительном временном интервале (около 20 тысяч шагов по времени с шагом —Ь, соответствующим шагу основной трехмерной задачи

(8)-(11)). Было установлено, что для одномерной задачи искусственная вязкость вносит незначительные искажения решения в окрестности ударной волны, отличие точного решения от численного составляло не более 2-3 %.

Отметим, что схема типа Мак-Кормака (8), (9) является неконсервативной, однако проведенные тестовые решения для цилиндрической области показали вполне приемлемую точность решения, получаемую с помощью этой схемы, с незначительной погрешностью численных результатов (в пределах 3-4 %). В то же время предложенная схема (8),

(9) обладает определенными преимуществами по сравнению с консервативными схемами. Оказалось, что для многомерных адаптивных сеток схема (8), (9) значительно менее чувствительна к "качеству" разностной сетки: она требует лишь непрерывности и достаточно приблизительной гладкости сетки. В то е время консервативная схема обеспечивает возможность вычислений только лишь на сетках с "идеальной" гладкостью, генерация таких сеток для многомерных областей сложной формы, подобных рассмотренным в настоящей работе, является достаточно непростой дополнительной задачей. Кроме того, известно (см., например, [11]), что на консервативных схемах часто хуже, чем на неконсервативных, осуществляется расчет косых скачков уплотнения, характерных для данной задачи.

Аппроксимация граничных условий. Для граничных условий (2), (3) применялась следующая разностная аппроксимация.

Граничные условия на жестких стенках:

т+1/3 _ „ т+1/3 _ „ т+1/3 _ „

= 0, 'идгг0к = 0, 'двг0к = 0,

т+1/3 т+1/3 т+1/3 т+1/3

Одг0к = Одг1к , рдг1к = рдг0к .

Условия на 0си симметрии:

т+1/3 т+1/3 т+1/3 т+1/3 т+1/3

'дг 0зк = 0, 'дг 0зк = 'дг 1зк , 'дв 0зк = 'дв 1зк ,

©m+1/3 _ ©m+1/3 ^m+1/3 _ ^m+1/3

'}9 0зк _ ^9 Цк , Ря 0зк _ р91^к

Условия на поверхности горения: 1

Ря

Граничные условия по окружной координате:

1 ( \ т\ г\ гл* m+1/3 m+1/3

v • n _ - — (PS - pg)D, ©g , PgjVi/fc _ Pg(N1-1)Jfc.

т+1/3 _ ^ т+1/3 _ т+1/3 т+1/3 _ т+1/3

УдО^Ыз _ 0, У9'ГгуКз _ ^9гг?(Жз-1), ^Я^Ыз _ Уд^(Ыз-1)5

/-ч т+1/3 _ /-ч т+1/3 т+1/3 _ т+1/3

^Яг^з _ "яг^(Ыз-1)' р9г^з _ р9гЯЫ3-1)'

Граничные условия в выходном сечении сопла:

/ \т+1/3 _ / \т+1/3 / \т+1/3 _ / \т+1/3

(У9Г)гЫ2к _ (У9Г)г(Ы2-1)к, (У9О/г^к _ (У9О)г(^2-1)к,

(у )т + 1/3 — (V )т + 1/3 От + 1/3 _ От + 1/3 рт+1/3 _ р

Л^гк _ )г(Ы2-1)к, ^Яг^к _ ^Яг(Ы2-1)к, рЯгЫ2к _ рде-

Результаты численного моделирования. При численном решении

задачи использовались следующие данные: Ро—1 • 107 Па, До_0,01 м/с,

п

V _ 0,4. Угол 0 звездности заряда (см. рис. 1) был выбран равным —.

16

На рис. 2-7 представлены некоторые результата расчетов нестационарной трехмерной задачи газовой динамики с вязким газом в различные моменты времени. В табл. 1-3 представлены значения давления, температуры, осевой и окружной составляющих скорости в различные моменты времени, полученные в результате решения двумерной (без учета звездности области) и трехмерной задач газовой динамики в случаях идеального и вязкого газов.

Рис. 2. Вид конечно-разностной сетки, использованной в задаче:

КС — камера сгорания; ДО — донная область, КР — критическое сечение; СБ — сопловой блок

Рис. 3. Распределения температуры в камере сгорания РДТТ для моментов времени 0,3,0,5, 2 и 4 мс соответственно

Таблица 1

Сравнительные значения параметров газа при различных моделях расчета

и при I = 0,00016 с

Параметры 2D, идеальный газ 2D, вязкий газ 3D, идеальный газ 3D, вязкий газ

р, в камере сго- 30,26 29,84 30,79 29,13

рания, атм

р, в критическом 3,71 2,80 3,67 2,64

сечении, атм

©5, К 1624,7 1631,1 1741,5 1689,2

'Ог, м/с 308,2 305,6 319,3 313,7

О0, м/с 0 0 0,7 0,7

Рис. 4. Распределения поперечной скорости в камере сгорания РДТТ для моментов времени 0,3, 0,5,2 и 4 мс соответственно

Рис. 5. Распределения давления в камере сгорания для моментов времени 0,3, 0,5, 2 и 4 мс соответственно

Рис. 6. Распределения продольной скорости в камере сгорания РДТТ для моментов времени 0,3, 0,5,2 и 4 соответственно

Таблица 2

Сравнительные значения параметров газа при различных моделях расчета

и при г = 0,0007 с

Параметры 2D, идеальный газ 2D, вязкий газ 3D, идеальный газ 3D, вязкий газ

р, в камере сго- 37,92 37,63 37,30 35,35

рания, атм

р, в критическом 26 26,68 26,87 26,20

сечении, атм

©5, К 2171,5 2213,2 2210,8 2207,1

'0г, м/с 866,9 819 881,8 796,5

Ьв, м/с 0 0 0,3 0,3

Графики на рис. 7 показывают распределение давления в зависимости от времени решения нестационарной трехмерной задачи с учетом вязкости газа. Значения рассчитаны в точках 1-5, представленных на рис. 1.

В целом, ход развития газодинамических процессов подобен тому, который был установлен ранее в работах [6-8] в рамках осесимметрич-

Таблица 3

Сравнительные значения параметров газа при различных моделях расчета и при t = 0,00423 с

Параметры 2D, идеальный газ 2D, вязкий газ 3D, идеальный газ 3D, вязкий газ

р, в камере сго- 87,7 99,2 93,97 103,68

рания, атм

р, в критическом 73,66 90,62 78,83 92,96

сечении, атм

Од , К 3505,3 3607,9 3639,1 3691,5

, м/с 1859,7 1794,7 1925,9 1812,7

У$, м/с 0 0 0,2 0,2

0,00371 0,00372 0,00404 0,00420

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

Рис. 7. Зависимость давления в рабочем тракте РДТТ от времени в фиксированных точках 1—5, показанных на рис. 1

эта скорость оказывается в случае вязкого газа примерно на 5 % ниже (различие около 100 м/с) по сравнению с моделью идеального газа (см. табл. 3). Максимальное даление в камере сгорания, наоборот, в случае вязкого газа выше примерно на 7 % (абсолютное значение 8-9 атм).

В результате численного моделирования установлено, что наличие в рабочем тракте донной области (см. рис. 2), что характерно для рассматриваемого типа двигателей, существенно влияет на распределение параметров в РДТТ. Так, давление в донной области примерно на 78 % выше, чем в камере сгорания на установившемся режиме, а радиальная компонента скорости уг в донной области имеет значительные отрицательные значения (- 112 м/с) при этом режиме. Максимального значения компонент уг достигает в сопловом блоке и на установившемся режиме оно составляет 316 м/с.

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

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

СПИСОК ЛИТЕРАТУРЫ

1. П и р у м о в У. Г., Росляков Г. С. Численные методы газовой динамики. -М.: Высшая школа, 1987.

2. Шевелев Ю. Д. Пространственные задачи вычислительной аэрогидродинамики. - М.: Наука, 1986.

3. Самарский А. А., Попов Ю. П. Разностные методы решения задач газовой динамики. - М.: Наука, 1992.

4. Рихтмайер Р., Мортон К. Разностные методы решения краевых задач. - М.: Мир, 1972.

5. Гильманов А. Н. Методы адаптивных сеток в задачах газовой динамики. -М.: Наука, 2000.

6. Димитриенко Ю. И., Димитриенко И. Д. Численное моделирование процессов горения смесевых твердых топлив // Вестник МГТУ им. Н.Э. Баумана. Сер. "Естественные науки". - 2001. - № 2. - С. 9-23.

7. Димитриенко Ю. И., Димитриенко И. Д. Численное моделирование нестационарных процессов горения в модельных РДТТ // Вопросы оборонной техники. Сер. 15. - 2002. - Вып. 1-2. - С. 3-10.

8. Д имитриенко Ю .И., И з о т о в а С. Г. Численное исследование нестационарных газодинамических процессов горения твердых топлив в камере РДТТ // Аэрокосмические технологии: Труды Всероссийской научно-технической конференции / Под ред. Р.П. Симоньянца. -М.: Изд-во МГТУ им. Н.Э. Баумана, 2002. - С. 144-154.

9. Димитриенко Ю. И., Ануфриев С. Н., Изотова С. Г. Разработка метода решения трехмерной нестационарной внутренней задачи газовой динамики на многопроцессорных вычислительных системах // Математика в современном мире / Под ред. Ю.А. Дробышева. - Калуга: Изд-во КГПУ, 2004. - С. 139-146.

10. Димитриенко Ю. И., Ануфриев С. Н., Изотова С. Г. Разработка технологии параллельных вычислений для решения задач газовой динамики в соплах РДТТ // Аэрокосмические технологии. - М.: Изд-во МГТУ им. Н.Э. Баумана, 2004.- С. 109-111.

11. Роуч П. Вычислительная гидродинамика. -М.: Мир, 1980.

Статья поступила в редакцию 19.11.2004

Юрий Иванович Димитриенко родился в 1962 г., окончил в 1984 г. МГУ им. М.В. Ломоносова. Д-р физ.-мат. наук, профессор, зав. кафедрой "Вычислительная математика и математическая физика" МГТУ им. Н.Э. Баумана, академик РАИН. Автор более 100 научных работ в области вычислительной механики, нелинейного тензорного анализа, термомеханики композитов, математического моделирования в материаловедении.

Yu.I. Dimitrienko (b. 1962) graduated from the Lomonosov Moscow State University in 1984. D. Sc. (Phys.-Math.), professor, head of "Computing Mathematics and Mathematical Physics" department of the Bauman Moscow State Technical University, full member of the Russian Academy of Engineering Sciences. Author of more than 100 publications in the field of computing mechanics, nonlinear tensor analysis, thermomechanics of composites, mathematical simulation in science of materials.

Светлана Геннадиевна Изотова родилась в 1980 г., окончила в 2003 г. МГТУ им. Н.Э. Баумана. Аспирант кафедры "Вычислительная математика и математическая физика" МГТУ им. Н.Э. Баумана. Автор ряда работ по численным методам газовой динамики.

S.G. Izotova (b. 1980) graduated from the Bauman Moscow State Technical University in 2003. Post-graduate of "Computing Mathematics and Mathematical Physics" department of the Bauman Moscow State Technical University. Author of a number of publications in the field of numerical methods of gas dynamics.

Андрей Алексеевич Захаров родился в 1982 г., окончил в 2005 г. МГТУ им. Н.Э. Баумана. Аспирант кафедры "Вычислительная математика и математическая физика" МГТУ им. Н.Э. Баумана. Специализируется в области вычислительной газовой динамики.

A.A. Zacharov A.A. (b. 1982) graduated from the Bauman Moscow State Technical University in 2005. Post-graduate of "Computing Mathematics and Mathematical Physics" department of the Bauman Moscow State Technical University. Are specialised in the field of numerical methods of gas dynamics.

Сергей Николаевич Ануфриев родился в 1982 г., окончил в 2005 г. МГТУ им. Н.Э. Баумана. Аспирант кафедры "Вычислительная математика и математическая физика" МГТУ им. Н.Э. Баумана. Автор научной работы по численным методам газовой динамики.

S.N. Anufriev (b. 1982) graduated from the Bauman Moscow State Technical University in 2005. Post-graduate of "Computing Mathematics and Mathematical Physics" department of the Bauman Moscow State Technical University. Author of a publication in the field of numerical methods of gas dynamics.