Численное исследование динамики внутрикамерных процессов ракетного двигателя на твёрдом топливе особой компоновочной схемы. Часть 1. Постановка вычислительного эксперимента Текст научной статьи по специальности «Механика и машиностроение»

DOI: 10.15593/2224-9982/2018.53.06 УДК: 519.9, 629.7

М.Ю. Егоров1, Д.М. Егоров2, С.М. Егоров2

Пермский национальный исследовательский политехнический университет, Пермь, Россия

АО «Научно-исследовательский институт полимерных материалов», Пермь, Россия

ЧИСЛЕННОЕ ИССЛЕДОВАНИЕ ДИНАМИКИ ВНУТРИКАМЕРНЫХ ПРОЦЕССОВ РАКЕТНОГО ДВИГАТЕЛЯ НА ТВЁРДОМ ТОПЛИВЕ ОСОБОЙ КОМПОНОВОЧНОЙ СХЕМЫ. ЧАСТЬ 1. ПОСТАНОВКА ВЫЧИСЛИТЕЛЬНОГО

ЭКСПЕРИМЕНТА

Численное исследование динамики внутрикамерных процессов является одной из главных задач при разработке и проектировании ракетного двигателя на твёрдом топливе (РДТТ). Современный РДТТ - сложная техническая система, в которой одновременно протекает ряд взаимосвязанных нестационарных и нелинейных физико-химических процессов. Рассматриваемый РДТТ имеет свои функциональные и конструктивные особенности. Для оптимизации параметров ракетного двигателя предпринята попытка прямого численного исследования динамики его внутрикамерных процессов. Формулируется сопряжённая постановка задачи, включающая в себя: срабатывание воспламенителя; прогрев, воспламенение и последующее нестационарное и турбулентное горение заряда твёрдого топлива; нестационарное ударно-волновое и вихревое гомогенно-гетерогенное течение воздуха и продуктов сгорания в камере сгорания, газоходах и сопловых блоках; разгерметизацию камеры сгорания. Каждая из подзадач рассматривается во взаимосвязи и разрешается одновременно - на одном шаге по времени. Для решения поставленной задачи разработаны физико-математическая модель, комплекс прикладных программ на ЭВМ и произведено их тестирование.

Ключевые слова: ракетный двигатель на твёрдом топливе, внутрикамерные процессы, горение, газовая динамика, постановка вычислительного эксперимента.

M.Yu. Egorov1, D.M. Egorov2, S.M. Egorov2

1

Perm National Research Polytechnic University, Perm, Russian Federation AO "Research Institute of Polymeric Materials", Perm, Russian Federation

NUMERICAL STUDY OF THE DYNAMICS OF INTRA-CHAMBER PROCESSES ROCKET ENGINE SOLID FUEL THE SPECIAL LAYOUT SCHEME. PART 1. STATEMENT OF COMPUTER EXPERIMENT

Numerical study of dynamics of intracameral processes is one of the main tasks in the development and design of a rocket engine on solid fuel (RESF). Modern RESF is a complex technical system in which a number of interrelated non-stationary and nonlinear physico-chemical processes take place simultaneously. Rocket engine solid fuel the special layout scheme has its own functional and design features. To optimize the parameters of the rocket engine, an attempt is made to directly numerical study the dynamics of its intracameral processes. The conjugate formulation of the problem, which includes: -operation of the igniter; - heating, ignition and subsequent non-stationary and turbulent combustion of a solid fuel charge; - non-stationary shock wave and vortex homogeneous-heterogeneous flow of air and combustion products in the combustion chamber, flues and nozzle blocks; - depressurization of the combustion chamber. Each of the subtasks is considered in the relationship and resolved at the same time-one step at a time. To solve this problem, a physical and mathematical model, a set of computer applications and their testing were developed.

Keywords: rocket engine solid fuel, intrachamber processes, combustion, gas dynamics, formulation of the numerical experiment.

Введение

Численное исследование динамики внутрикамерных процессов является сегодня одной из главных задач при разработке и проектировании ракетного двигателя на твёрдом топливе (РДТТ) любого типа и назначения [1-3]. В рамках этой задачи определяются рабочие парамет-

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

Современный РДТТ - сложная техническая система, в которой одновременно протекает ряд взаимосвязанных (сопряжённых) нестационарных и существенно нелинейных физико-химических процессов [4-7]. Рассматриваемый РДТТ, с учётом его функций, имеет свои конструктивные особенности. Общий вид РДТТ показан на рис. 1.

Рис. 1. РДТТ особой компоновочной схемы

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

Для оптимизации параметров РДТТ особой компоновочной схемы в предлагаемой работе предпринята попытка прямого численного исследования (на уровне постановки вычислительного эксперимента) динамики его внутрикамерных процессов.

Рассматривается сопряжённая постановка задачи, включающая в себя:

- срабатывание воспламенительного устройства;

- нестационарный прогрев, воспламенение и последующее нестационарное и турбулентное горение заряда твёрдого топлива;

- нестационарное ударно-волновое и вихревое гомогенно-гетерогенное течение воздуха, продуктов сгорания воспламенительного состава и заряда твёрдого топлива в камере сгорания, боковых газоходах и сопловых блоках;

- разгерметизацию камеры сгорания (вскрытие заглушек сопловых блоков).

Каждая из подзадач рассматривается во взаимосвязи и разрешается одновременно - на одном шаге по времени.

Для решения поставленной задачи разработаны физико-математическая модель, комплекс прикладных программ на ЭВМ, и произведено их тестирование.

Срабатывание воспламенительного устройства

Процесс срабатывания воспламенительного устройства (ВУ) с перфорированным корпусом и вкладным зарядом воспламенительного состава описывается на основе экспериментально-теоретического подхода. Путём решения обратной задачи внутренней баллистики для системы «ВУ - имитатор свободного объёма камеры сгорания РДТТ» рассчитывается реальный газоприход от ВУ с учётом особенности горения заряда воспламенительного состава - догорания продуктов сгорания за корпусом ВУ. Продукты сгорания воспламенительного состава рассматриваются как газ с «эффективными» показателем адиабаты и газовой постоянной, учитывающими наличие в газе твёрдой фазы.

Основная система дифференциальных уравнений процесса срабатывания ВУ РДТТ имеет следующий вид:

- изменение давления в корпусе ВУ

К ■ Тв

л Гв

- изменение давления в камере сгорания РДТТ

Лк = Ккс ■ Ткс Л Г

•(рГ-рв)-»2 ; С1)

т. (2)

В (1)-(2) и далее по тексту приняты обозначения: т - массовый расход (приход); Р -давление; К - газовая постоянная; 5 - площадь поверхности горения; Т - температура; I -время; V - объём; V - скорость горения; р - плотность; в - воспламенитель, кс - камера сгорания; ш - шашка.

Выражая из (1) массовый секундный расход из корпуса ВУ в камеру сгорания РДТТ - т, подставляя его в (2) и разрешая полученное относительно скорости горения заряда ВУ, имеем

1 ( V ЛР ^ ЛР„„}

V3/

(рш-Рв К К ■ Тв Л Ккс ■ ткс л

(Рв -р

Выражение (3) является основным расчётным соотношением, в котором давление в корпусе ВУ - Р и давление в камере сгорания РДТТ - Р и их производные по времени определяются экспериментально на специально созданной лабораторной установке. Однако такой способ определения скорости горения неудобен, так как для каждого конкретного ВУ необходимо провести эксперимент по замеру рабочего давления в корпусе ВУ и в камере сгорания РДТТ и в последующем произвести расчёт. Удобнее обобщить и напрямую связать значение скорости горения заряда ВУ ув с параметрами, определяющими процесс горения конкретного воспламенительного состава в корпусе ВУ (площадью перфорации корпуса ВУ - ^кр, площадью поверхности горения заряда ВУ - 5В и свободным объёмом в корпусе ВУ - V), в виде

V = / (^р, V )• Такого рода связи, не вникая в сущность сложного физико-химического

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

V = А + А • XI + А • X2 + А • XI-X2 + А • XI2 + А • X22 + А • XIX22;

XI =

Г

- д

V V

V В

/ Д; X 2 =

^ - с

V в У

/ с

(4)

2

а для нахождения коэффициентов полинома применить теорию планирования эксперимента [8]. Здесь А,•••, А,Д,Д,С,С - коэффициенты полинома, зависящие от конкретного вида воспламенительного состава и интервалов варьирования конструктивных параметров ВУ.

По полученной скорости горения заряда воспламенительного состава (4) определяется газоприход от ВУ в камеру сгорания РДТТ.

Воспламенение и горение заряда твёрдого топлива

Описание процесса нестационарного прогрева, воспламенения и последующего нестационарного и турбулентного горения состава заряда твёрдого топлива (ТТ) базируется на модели Мержанова-Дубовицкого с учётом (в рамках подхода Горохова-Липанова-Русяка) влияния газовой фазы на процесс горения в конденсированной фазе (£-фазе) [9 и др.]. Будем рассматривать ТТ как твёрдое тело, к которому применимы известные уравнения теплопроводности и химической кинетики. Для удобства будем рассматривать эти уравнения в системе координат, связанной с поверхностью горения, направив ось от поверхности в ТТ. Считаем, что реакции в £-фазе удовлетворяют закону Аррениуса. Тогда в предположении «0» - мерности порядка химических реакций данная система уравнений, описывающая процесс, имеет вид

8х

8%+л; Т , а

8У1 8Ук

8Т -к 8Т+^ 8Тк+^ (тк); 8т ^ ^+Фк (Т),

_8Р

8Ук

(5)

где

Фк (Тк ) = ^ • ехР

Л

Я • Т

к у

До воспламенения необходимо положить ук = 0. Условие воспламенения и последующего горения = р* = 1.

До воспламенения начальные и граничные условия для системы (5) имеют вид

X = 0, Ук > 0, Тк = Т0, Р = 0; 8Т

X > 0, Ук = 0, -к- = ат(Тя - Т ), р<р„;

8У,

X > 0, Ук = <ю, Тк = Т), р = 0. После воспламенения граничные условия запишутся в виде

X > X*, Ук = 0, Тк = Т, р* =р.,

(6)

^к= Я -(Ср* -ск)рк • V ■Т -

8У1

г 1 ^

а -- р

V Рку

Т, Т 0, р= 0.

рк- ч •р;

(7)

Задача (5)-(7), дополненная замыкающими соотношениями, обезразмеренная и записанная в неравномерной (экспоненциальной) системе координат, интегрируется численным сеточ-

ным методом по явным и неявным конечно-разностным схемам [9, 10 и др.]. Неявные схемы разрешаются способом прогонки.

Температура по своду к-фазы, исключая поверхность горения твёрдого топлива,

T+1 = a • j + в; 2 < i < M-1. (8)

Температура на поверхности горения ТТ: до воспламенения

T+1 = A • Ti+1 + B; (9)

после воспламенения

TJ+1 = A* • t/+1 + в* (10)

Глубина превращения ТТ: до воспламенения

Р/+1 = p'+A • а3 Wj; 1 < i < M -1; (11)

после воспламенения

Pj+1 = A** + B**; 1 < i < M -1. (12)

Скорость горения ТТ vk определяется из (12), с учётом граничных условий

Pi+1 = Pi = Р*; Рм1 = Рм = 0, итерационным способом (методом секущих).

В (5)-(12) и далее по тексту приняты обозначения: к - коэффициент температуропроводности; y - координата; Q - тепловой эффект реакции; c - удельная теплоёмкость; в - глубина превращения топлива; Z - предэкспонент; E - энергия активации; R - универсальная газовая постоянная; X - коэффициент теплопроводности; ах - коэффициент теплоотдачи; q -плотность теплового потока; a - коволюм; A, B - коэффициенты прогонки; At - шаг по времени; а - const; W - безразмерная экспоненциальная функция температуры; к - конденсированная фаза (твёрдое топливо); g - газ; 5 - поверхность горения; p - давление; i - номер расчётной точки по координате; j - номер расчётной точки по времени; * - специальное значение.

По полученной скорости горения vt определяется газоприход с поверхности горения заряда ТТ в камеру сгорания РДТТ.

Газовая динамика в камере сгорания РДТТ

Для математического описания процесса течения в камере сгорания, боковых газоходах и сопловых блоках РДТТ используем подходы механики сплошных многофазных сред [11, 12]. Согласно этим подходам воздух в камере сгорания, газообразные продукты сгорания воспла-менительного состава и заряда ТТ назовём первой фазой, мелкодисперсные частицы (жидкие капли) в продуктах сгорания воспламенительного состава и заряда ТТ - второй фазой. Первую и вторую фазы будем считать гомогенно-гетерогенной смесью со своими температурами и скоростями движения. В такой системе каждая фаза занимает часть объёма смеси: а, а2. Здесь

а =рг- /p"c (i = 1,2), где рги° - истинная плотность фазы. Движение этих фаз рассматривается как движение взаимопроникающих и взаимодействующих сред.

Тогда полная (нестационарная и трёхмерная) система вихревых дифференциальных уравнений газовой динамики для гомогенно-гетерогенного потока в камере сгорания, боковых газоходах и сопловых блоках ракетного двигателя запишется в виде: уравнения неразрывности (сохранения массы)

^ + div (P1W1 ) = + Gge,

я (13)

^ + div (р2 W2 ) = Gp,

5(р,ф)

-+ ^ (р1фW1) = ф^ • О^ + ф6 • О^, ф = к, Ср, ц Я, а;

&

уравнения сохранения импульса по осям координат

а(рЛ) + а1у(рЛ ^) + «! ^ = Ч2 + ^ ■ О^ + Жхе ■

5х

12

^^ + а1у(рlVlWl) + «1 -Р = -х^2 + ^ • О^^ + Гув • ; (14)

5(Р1^1) + ё1у(р^) + «1 ^ = -х1/ + Жги, • О^ + • ОЯв;

5Л 1 1 " 1 -г г ™ **

+ d1v(р2U2W2)+«2 ■ -Р = х!2 + ■ Ор;

^ + d1v(р2V2W2)+ а2 ■ ^ = х- + ^ ■ Ор;

^ + d1v(р2^W2)+«2 ■ -Р = х- + ■ Ор;

уравнения сохранения внутренней удельной энергии второй фазы

5(р2 Л )

^ +div(р2J2W2) = д]? + д]2 +Зр ■ Ор; (15)

уравнения сохранения полной удельной энергии смеси

—+ + d1v (р^ ) + div (р2Е2 ^ ) +

+ d1v(«р^) + d1v(а2PW2 ) = Е&в ■ О?е + Е^ ■ О^ + Ер ■ Ор, (16)

где для декартовой системы координат

d1v ) = + + , Ир* , ри, PiVi, р^., р^2, PiEi, ауР]; i = (1,2). (17)

Для замыкания системы дифференциальных уравнений (13)-(17) используем уравнение состояния в виде

р = (к-1 )рис

ж

2 Л

Е

V 1 2 У

1

(18)

1 ... ис

1 - а Ф!

В (13)-(18) и далее по тексту приняты обозначения: с - удельная теплоёмкость при постоянном давлении; Е - полная удельная энергия; J - внутренняя удельная энергия; О - рас-ходно-приходный комплекс; к - показатель адиабаты; д - функция теплового межфазного взаимодействия; и, V, м> - проекции вектора скорости по осям координат 0X,0Y,0Z; Ж - модуль вектора скорости; W - вектор скорости; х, у, г - текущие координаты вдоль осей 0Х, 07, 0Z; ц - динамическая вязкость; х - функция силового межфазного взаимодействия; g - газ; к - конвективный; I - лучистый; р - частицы (жидкие капли); ^ - заряд ТТ; х -вдоль оси 0Х; у - вдоль оси 07; г - вдоль оси 0Z.

Выражения в (13)-(16) для расходно-приходных комплексов, функций силового и теплового межфазного взаимодействия, а также используемые в расчётах дополнительные соотношения подробно изложены в [9].

Система уравнений (13)-(18) с учётом дополнительных соотношений, интегрируется численно с помощью метода Давыдова (метода крупных частиц) - метода постановки вычислительного эксперимента [4-7, 9, 13-15]. Область интегрирования покрывается фиксированной в трёхмерном пространстве (эйлеровой) равномерной (однородной и полностью изотропной) кубической расчётной сеткой с ячейками Ах х Ау х Az. Для такой сетки геометрический центр ячейки совпадает с центром масс ячейки. Значения целых чисел <«» (вдоль оси 0Х), «у» (вдоль оси 07) и «к» (вдоль оси 0. обозначают геометрический центр ячейки, «п» - порядковый номер шага по времени.

Эйлеров этап метода. На этом этапе расчёта изменяются величины, относящиеся к ячейке в целом, а исследуемая среда предполагается заторможенной:

pw pw .

~п и п р+0.5,М ~ р-0.5,М At

Um — ot„ ----

Ax

mi,j,k

E" =Jn

2¡J,k 2iJ,k

К J2+teJ2.

m

= (i,2);

(аналогично по другим направлениям для v"n . w"n );

/ \ P?

Ёп =En -1Ё" -E" | iJ'k

Pli,. i,k

oc" P" ы" ос" P" ы"

1i+0.5, j,k í+°.5. j.k 1i+0.5,j ,k 1i-0.5,j,k 'j.k 1i—0.5,j,k At

Ax

P1

i,j,k

a" • P" • v" —a" • P" • v" л

1i, j+0.5,k i,J+0-5,k li, j+0.5,k 1i, j—0.5,k j-0.5.k 1i,j—0.5,k At

Ay

n" P" <T," n" P" <7,"

a • P í t+n ^ • w — a • p . t n с • w-,

1i, j ,k +0.5 .k+a5 1i, j,k+0.5 1i,j,k—0.5 г. М-0.5 1

P4,,k

i, j,k—0.5 At

Az

Pi

i,j ,k

где

Ax P"u ,k

a" 2i, j+0.5,k • p" p i,j+0.5,k • v" —a" 2i, j+0.5,k 2i,j—0.5,k p" ' P, j-0.5,k ' 4, —0.5,k At

Ay

a" 2i, j,k+0.5 • P" p í,j,k+0.5 • w" —a" 2i, j,k +0.5 2i,j ,k—0.5 p" • P, j,k - 0.5 • w". ,j,k—0.5 At

Az

Pj

p" i p" p ¡ ъ + p-

p"

p+0.5, j,k

2

м" ^ = (1 - alfa) -u" ,k+ alfa ■ ü[

n

i.j.k

(19)

и т.д. В (19) a//a - сеточный параметр, значение которого alfa = 3,0.

Лагранжев этап метода. На данном этапе метода вычисляются эффекты переноса, учитывающие обмен между ячейками при их перестройке на прежнюю эйлерову сетку:

2

(рМ)

первая фаза вдоль оси 0Х

(l-beta)pljk+beta.plijk\i (I-beta) plijk +beta-pljk1f+l]k-u-

i+0.5,j,k

" -zT , если й" >0;

, , если й," < 0;

(20)

(аналогично для потоков по другим направлениям (р^у,5 к,; ¿+0 5); вторая фаза вдоль оси 0Х

и к и ~г

i+0.5,j,k

если и, > 0;

И i-H —и ^ п

если «2i+05Jji<05

(21)

(аналогично для потоков по другим направлениям (p2i;v'2 . (| _ ^ ,(р2<;и,2 . ^ (| _) и т.д. В (20)

beta - сеточный параметр, значение которого —0,15.

На лагранжевом этапе метода вычисляются также приходные комплексы и функции силового и теплового межфазного взаимодействия, входящие в (13)-(16), с учётом изменения параметров потока на эйлеровом этапе.

Заключительный этап метода. Здесь происходит перераспределение массы, импульса и энергии по пространству, и определяются окончательные поля эйлеровых параметров потока на фиксированной сетке в новый момент времени: уравнения неразрывности (сохранения массы)

и+1 _ n

P4j,k = РЧj,k

i-0.5,j,k

Ax

A t —

i,j-0.5,k

Ay

A t-

Ay

n+1 w

= P2j

^^At + lc" +G" ) A/:

(P 2^2)i+0.5,j,k (P2^2 )i-Q.5J,k д 1 /ij+05k

Ax

Af

A t-

(рА)г"М+0.5 "(pA)

i,j,k-0.5

Az

At + G;uyAt;

(22)

r.n+1 ,„n

ФЧj,k = Ф,М --"n+r P4. a

Pita Ai (Pi^i)I"/+o.5,it-(Pi(Pi7i)"y-o.5,it A t

Ax

p1,jk

Ay

n+1

i,j,k-0.5

Az

Pi,

At

Ф = ( k, cp, ц, A,, a ,; уравнения сохранения импульса по осям координат

P"u,k (pl"l"l )Г+0.5. /.Jt - (Pl"l«l )i

и+1 _~„ " ij,k \r^h+0.5,J,k

Ъ1л — Z-ii Л

\j,k \j,k »+1 д^

i-0.5,j,k

Д/

n+l

Рчм

_ <

w

(р^)1]ш5,к~(р^)г,]-0.5,к &

" Ау р^

Az Р;+1

,1 ,4

-X12" А1 + (ж" ■ О" + Ж" • а" ) А/ • (23)

х,1,к пп+1 \ ,к Щ,1,к gвi,j,к ) п+1 ' К '

р 111,к V 1 1 1 р 1,к

„+1 (р2»2Щ)"+0^к-(р2«2«2).

1-0.5,],к

At

Р 2,

¡,1 ,к

Ах

и+1

Р 2+

¡,1 ,к

1,]-0.5,к

Д/

Ау

(р2«2^2)г"М+о.5-(р2«2^2);

п+1

Р2,

г,У,¿-0.5

¡,1 ,к

At

Az

п+1

Р2

¡,1,к

_ 1 2"

А/

А/

+ Х12--+ Ж" ■ О"

х',1 ,к "+1 хж',1,к Р',1,к +1

Р 2',.1,к

Р2<„ 1 ,к

(аналогично по другим направлениям , , , );

уравнение сохранения внутренней удельной энергии второй фазы

л+1 = л

/¿-0.5J.it

А г

2',1,к 2',1,к „"+1

Р 2

¡',1 ,к

Ах

"+1 Р 2+

•,1,к

1,]-0.5,к

А t

Ау

„"+1 Р2+.,к

(р2-/2^2),"М+о.5-(р2^2)"М-о.5 А/ , 12» 12» . ^ _

"+1 ( ,1,к ^ ,1,к Р, 1,к Р, 1 ,к ) ^ "+1

А

Az

Р 2

'¡и ,к

Р2,

',1 ,к

уравнение сохранения полной удельной энергии смеси

1 77^+1 И

Р, р9 • — О, •

1; 1 к ~ 1; г,А _и+1

Г11,1,4 Г1!,1 ,к

-I

Ш=1

{ртЁтйт )(+о 5 - (?тЁтйт )

1-0.5,],к

Д/

Ах

/ /г ~ V Р ~ ^

уРш Ш^Ш Д. у+0 5 £ ш^ш у

1,]-0.5,к

А(

Ау р"

(р^Л ). . к+о 5 - (ри£иУи ). м_о 5 д,

Аг

¡,1 ,к

Р1

+(Е" ■ О" + Е" ■ О" + Е" ■ О"

¡, 1,к gвi, 1,к gwi,j,к ,к Pi,j,к Pi,1,к} "+1

РЧ 1,к

)

А/

(24)

(25)

Для повышения точности вычислений (обеспечения условия полной консервативности) в схему метода при расчёте давления по уравнению состояния (18) вводится поправка, обеспечивающая (уточняющая) баланс по внутренней удельной энергии смеси [9].

Постановка граничных условий. Для описания граничных условий вокруг расчётной области (свободный объём камеры сгорания, боковых газоходов и сопловых блоков (см. рис. 1)

вводятся слои фиктивных ячеек [9, 13, 14] (рис. 2). На закрытых поверхностях (стенки камеры ВУ, поверхность горения заряда ТТ, стенки камеры сгорания, газоходов и сопловых блоков, стенки заглушек сопловых блоков) используются условия непротекания: нормальная к поверхности компонента вектора скорости из приграничных ячеек сносится в слой фиктивных ячеек с противоположным знаком, остальные параметры потока из приграничных ячеек сносятся в слой фиктивных ячеек без изменения. На открытых поверхностях (расчётная область за сопловыми блоками) используется экстраполяция параметров потока. Приход продуктов сгорания от ВУ осуществляется в расчётные ячейки, центры масс которых расположены в потоке у ВУ. Приход продуктов сгорания с поверхности горения заряда ТТ осуществляется в расчётные ячейки, центры масс которых расположены в потоке у поверхности горения заряда ТТ.

На нерегулярных (несовпадающих с координатной сеткой) криволинейных границах расчётной области применяется предложенный Ю.М. Давыдовым аппарат дробных ячеек. Используется процедура нормального отображения фиктивной ячейки относительно границы расчётной области в поток. Везде применяются расчётные формулы только для целых ячеек.

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

Ф IX =1; ч = Ррк,-\а,арР^рЁр /=1,2. (26)

г г

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

% '(-^ • N + ^ • К2)]; = Дг, \-Wnj • Ny + Ш2] • Ny2)];

I г

= 1|Х-(-^ N + Wn2, * Nz 2)]; IV, = 1;

(27)

% =1Х {-Щ-лд]; **=Т[Гь {-щ-лд];

г г

7=1,2.

г

В (26)-(27) приняты обозначения: Wn, , Wn2 - проекции вектора скорости W на направляющие оси локальной системы координат (Wn = 0), N, N' N - проекции нормали к плоскости отображения в базовой системе координат, N2, N2, N2 - проекции оси локальной системы координат в базовой системе координат; а - фиктивная ячейка; в - отражённая фиктивная ячейка.

Разгерметизация камеры сгорания

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

Комплекс прикладных программ

Для проведения численных расчётов на ЭВМ разработан комплекс прикладных программ, включающий в себя следующие основные модули:

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

- основной расчётный модуль (main-модуль) - расчёт газодинамического течения в камере сгорания, боковых газоходах и сопловых блоках РДТТ с учётом работы воспламенительного устройства, прогрева, воспламенения и последующего горения заряда ТТ и вскрытия заглушек сопловых блоков;

- модуль визуализации полученной расчётной информации.

Комплекс прикладных программ написан в среде программирования Kdevelop 4.6 для рабочей станции высокой производительности с операционной системой Linux Mint 17.3 x86-64 на алгоритмическом языке C/C++ с использованием (для основного расчётного модуля) стандарта многопотоковой обработки информации OpenCL [16-19]. Основная идея стандарта OpenCL состоит в реализации многопотокового выполнения кода по схеме «одна команда — много данных», то есть одна операция одновременно применяется к большому массиву данных. Подобный подход позволяет в полной мере использовать вычислительный потенциал устройств с множеством относительно простых исполнительных модулей, например, таких как современные графические процессоры, и существенно (на порядок и более) повысить производительность вычислений.

Проверка работоспособности программного продукта

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

Наиболее критичной с точки зрения устойчивости и точности численного решения является задача газодинамического течения. Здесь для анализа свойств разностных схем метода Давыдова (см. выражения (19)-(25)) использовался эвристический подход, основанный на рассмотрении параболической формы их дифференциальных приближений [9, 13, 14]. В этом подходе оценивается знак коэффициентов диффузии aу диссипативных членов дифференциальных приближений, содержащих частные производные второго порядка по пространственным переменным. Эти коэффициенты обычно группируются в виде матрицы - т.н. матрицы аппрок-симационной вязкости. Положительность детерминанта, либо диагональных элементов, либо следа матрицы аппроксимационной вязкости является условием вычислительной устойчивости исследуемой конечно-разностной схемы. Приближая значения этих параметров к нулю, повышаем точность вычислений.

Согласно трактовке лагранжева этапа метода Давыдова [13, 14] каждая проекция конвективных членов исходной системы дифференциальных уравнений, (13)-(16) на координатную ось может аппроксимироваться независимо. Это позволяет, в рамках определённого приближения, при анализе свойств разностной схемы ограничиваться одной координатой. Выпишем для одномерного однофазного аналога используемой выше конечно-разностной схемы метода (19)—(25) диагональные элементы матрицы аппроксимационной вязкости. Будем считать, что поток течёт слева направо. Для противоположного направления потока достаточно поменять Ах на —Ах.

Итак, имеем

а,, =1 u -Ax -1 и -Ax2 +—I п - u2 ii 2 4 x 2\

1 (pP- u 2 )At;

1 A 1 A 2 1 f 3 L 1 л 2

а22 = ^'P 'U 'Ax ^2'P u 'Ax "eta I и 'P*'

+ 1 (U'Pjj'Jx + U ' Pj P ' P x + Pj'ux )■Ax" +

+ -2

j P x

-3' P' ll2 -P' Pp- — ' Pj - u 2' Pj

At;

(28)

a33 = 1P 'u 'Ax -1P u ' Ax2 -1 (1 - beta) u фх ' Ax2 -

-1 Pj • ux • Д^2 -1 (u'Pjjjx + u ' Ppj • Px ^ +

1

+ — 2

-P 'u + Pj

E -(1 - 2' alfa ) P

At,

где

5u 5 p

u~ =—< —pj и т-Д-

5x 5p5j

Условие положительности следа матрицы аппроксимационной вязкости (an + a22 + a33 )>

> 0 из (28) рассматривалось в качестве критерия устойчивости выбранной конечно-разностной схемы метода.

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

Заключение

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

Библиографический список

1. Алемасов В.Е., Дрегалин А.Ф., Тишин А.П. Теория ракетных двигателей. - М.: Машиностроение, 1980. - 533 с.

2. Численный эксперимент в теории РДТТ / А.М. Липанов, В.П. Бобрышев, А.В. Алиев [и др.]; под ред. А.М. Липанова. - Екатеринбург: Наука, 1994. - 301 с.

3. Газовые течения в соплах энергоустановок / под ред. проф. В.Н. Емельянова. - М.: Физматлит, 2017. - 328 с.

4. Егоров М.Ю., Егоров Я.В., Егоров С.М. Исследование неустойчивости рабочего процесса в двухкамерном РДТТ // Изв. вузов. Авиационная техника. - 2007. - № 4. - С. 39-43.

5. Егоров М.Ю., Егоров С.М., Егоров Д.М. Численное исследование переходных внутрикамер-ных процессов при выходе на режим работы РДТТ // Изв. вузов. Авиационная техника. - 2010. - № 3. -С. 41-45.

6. Егоров М.Ю., Егоров Д.М. Численное моделирование внутрикамерных процессов в бессопловом РДТТ / под ред. Г.В. Кузнецова, В.Н. Ускова, С.К. Матвеева, В.И. Запрягаева, И.К. Жарова, Е.Е. Бульба, Б.В. Борисова, А.В. Захаревич, Е.А. Маслова // XXIII семинар по струйным, отрывным и нестационарным течениям (с международным участием): сб. тр. - Томск: Изд-во Нац. исслед. Том. поли-техн. ун-та, 2012. - С. 124-127.

7. Егоров М.Ю. Численное исследование динамики внутрикамерных процессов при срабатывании специализированного РДТТ // Изв. вузов. Авиационная техника. - 2017. - № 4. - С. 104-111.

8. Спиридонов А.А. Планирование эксперимента при исследовании технологических процессов. -М.: Машиностроение, 1981. - 184 с.

9. Давыдов Ю.М., Егоров М.Ю. Численное моделирование нестационарных переходных процессов в активных и реактивных двигателях / под ред. Ю.М. Давыдова. - М.: НАПН РФ, 1999. - 272 с.

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

11. Рахматулин Х.А. Основы газовой динамики взаимопроникающих движений сплошных сред // ПММ. - 1956. - Т. XX, № 2. - С. 184-195.

12. Нигматулин Р.И. Динамика многофазных сред. - М.: Наука, 1987. - Ч. I. - 464 с.; - Ч. II. - 360 с.

13. Белоцерковский О.М., Давыдов Ю.М. Метод крупных частиц в газовой динамике. Вычислительный эксперимент. - М.: Наука, 1982. - 392 с.

14. Давыдов Ю.М. Крупных частиц метод // Математика. Большой энциклопедический словарь. -Изд-е 3-е. - М.: Российская энциклопедия, 1998/2000. - С. 303-304.

15. Егоров М.Ю. Метод Давыдова - современный метод постановки вычислительного эксперимента в ракетном твёрдотопливном двигателестроении // Вестник Пермского национального исследовательского политехнического университета. Аэрокосмическая техника. - 2014. - № 37. - С. 6-70.

16. Стахнов А.А. Linux. - СПб.: БХВ-Петербург, 2004. - 912 с.

17. Дерк Л. С и С++. Справочник. - М.: ВКК, 1997. - 592 с.

18. Programming Guide. AMD Accelerated Parallel Processing OpenCL. Advanced Micro Devices. -

2011.

19. Егоров М.Ю., Егоров С.М., Егоров Д.М. Применение графических ускорителей для повышения производительности вычислений при численном моделировании функционирования сложных технических систем // Вестник Пермского национального исследовательского политехнического университета. Аэрокосмическая техника. - 2015. - № 40. - С. 81-91.

References

1. Alemasov V.E., Dregalin A.F., Tishin A.P. Teoriya raketnykh dvigateley [Theory of rocket engines]. Moscow: Mashinostroyeniye, 1980, 533 p.

2. Lipanov A.M., Bobryshev V.P., Aliyev A.V. Chislennyy eksperiment v teorii RDTT [Numerical experiment in the theory of solid propellant rocket motors]. Ekaterinburg: UIF «Nauka», 1994, 301 p.

3. Emelyanov V.N. Gazovyye techeniya v soplakh energoustanovok [Gas flows in nozzles of power plants]. Moscow: FIZMATLIT, 2017, 328 p.

4. Egorov M.Yu., Egorov Ya.V., Egorov S.M. Issledovanie neustoychivosti rabochego protsessa v dvu-khkamernom RDTT [Study of working process instability in the two chamber solid propellant rocket engine]. Izvestiya vysshikh uchebnykh zavedeniy. Russian Aeronautics, 2007, no. 4, pp. 39-43.

5. Egorov M.Yu., Egorov S.M., Egorov D.M. Chislennoe issledovanie perekhodnykh vnutrikamernykh protsessov pri vykhode na rezhim raboty RDTT [Numerical study of transient interchamber processes when reaching the SPRE]. Izvestiya vysshikh uchebnykh zavedeniy. Russian Aeronautics, 2010, no. 3, pp. 41-45.

6. Egorov M.Yu., Egorov D.M. Chislennoye modelirovaniye vnutrikamernykh protsessov v bessoplovom RDTT [Numerical Modeling Interchamber processes in at nozzleless solid propellant rocket motors]. Editors: Kuznetsov G.V., Uskov V.N., Matveev S.K., Zapryahaev V.I., Zharov I.K., Tuber E.E., Borisov B.V., Zakha-

revych A.V., Maslov E.A. In the collections of: XXIII seminar on the jet, tear and unsteady flow (with international participation). Collections of Labor. National Research Tomsk Polytechnic University, 2012, pp. 124-127.

7. Egorov M.Yu. Chislennoe issledovanie dinamiky vnutrikamernykh protsessov pri srabatyvaniy specializirovannogo RDTT [Numerical Research of Intra-Chamber Processes Dynamics during Startup of Special Solid Propellant Engine]. Izvestiya vysshikh uchebnykh zavedeniy. Russian Aeronautics, 2017, no. 4, pp. 104-111.

8. Spiridonov A.A. Planirovaniye eksperimenta pri issledovanii tekhnologicheskikh protsessov [Experiment Planning in the study of technological processes]. Moscow: Mashinostroyeniye, 1981, 184 p.

9. Davydov Yu.M., Egorov M.Yu. CHislennoye modelirovaniye nestatsionarnykh perekhodnykh protsessov v aktivnykh i reaktivnykh dvigatelyakh [Numerical simulation of unsteady transition processes in active and jet engines]. Moscow: NAPS Russia, 1999, 272 p.

10. Ryhtmayer R.D., Morton H. Raznostnyye metody resheniya krayevykh zadach [Difference methods for solving boundary value problems]. Moscow: Mir, 1972, 420 p.

11. Rakhmatulin H.A. Osnovy gazovoy dinamiki vzaimopronikayushchikh dvizheniy sploshnykh sred [Fundamentals of gas dynamics of interpenetrating motions of continuous media]. PMM, 1956, vol. XX, no. 2, pp. 184-195.

12. Nigmatulin R.I. Dinamika mnogofaznykh sred [Dynamics of multiphase media]. Moscow: Nauka, 1987, Part I, 464 p, Part II, 360 p.

13. Belotserkovsky O. M., Davydov Yu. M. Metod krupnykh chastits v gazovoy dinamike. Vychislitelnyy eksperiment [Method of large particles in gas dynamics. Computational experiment]. Moscow: Nauka, 1982, 392 p.

14. Davydov Y.M. Krupnykh chastits metod [Large particle method]. In: Mathematics. Large Encyclopedic Dictionary. Publishing is 3rd ed. Moscow: Rossiyskaya entsiklopediya, 1998/2000, pp. 303-304.

15. Egorov M.Y. Metod Davydova - sovremennyy metod postanovki vychislitelnogo eksperimenta v raketnom ^rdotoplivnom dvigatelestroyenii [Davydov method - modern method of statement of computing experiment in solid rocket motors industry/. PNRPUAerospace Engineering Bulletin, 2014, no. 37, pp. 6-70.

16. Stahnov A.A. Linux. St.Peterburg.: BHV-Petersburg, 2004, 912 p.

17. Derk L. Si and Si++. Spravochnik [C and C ++. Directory]. Moscow: WCC, 1997, 592 p.

18. Programming Guide. AMD Accelerated Parallel Processing OpenCL. Advanced Micro Devices, 2011

19. Egorov M.Y., Egorov S.M., Egorov D.M. Primeneniye graficheskikh uskoriteley dlya povysheniya proizvoditelnosti vychisleniy pri chislennom modelirovanii funktsionirovaniya slozhnykh tekhnicheskikh sistem [The use of graphics accelerators to enhance computing performance for numerical simulation of the functioning of complex technical systems]. PNRPU Aerospace Engineering Bulletin, 2015, no. 40, pp. 81-91.

Об авторах

Егоров Михаил Юрьевич (Пермь, Россия) - доктор физико-математических наук, профессор, профессор кафедры «Высшая математика» ФПММ, ФГБОУ ВО ПНИПУ (614990, г. Пермь, Комсомольский пр., д. 29, e-mail: egorov-m-j@yandex.ru).

Егоров Дмитрий Михайлович (Пермь, Россия) - кандидат технических наук, зам. главного конструктора, АО «Научно-исследовательский институт полимерных материалов» (614113, г. Пермь, ул. Чистопольская, д. 16, e-mail: egorovdimitriy@mail.ru).

Егоров Сергей Михайлович (Пермь, Россия) - кандидат физико-математических наук, ведущий научный сотрудник, АО «Научно-исследовательский институт полимерных материалов» (614113, г. Пермь, ул. Чистопольская, д. 16, e-mail: know_nothing@bk.ru).

About the authors

Mikhail Yu. Egorov (Perm, Russian Federation) - Doctor of Physical and Mathematical Sciences, Professor, Department of Higher Mathematics, Perm National Research Polytechnic University (29, Komsomolsky av., Perm, 614990, Russian Federation, e-mail: egorov-m-j@yandex.ru).

Dmitry M. Egorov (Perm, Russian Federation) - CSc in Technical Sciences. Deputy Leading Engineer, JSC "Scientific-Research Institute of Polymeric Materials" (16, Chistopolskaya st., Perm, 614113, Russian Federation, e-mail: egorovdimitriy@mail.ru).

Sergey M. Egorov (Perm, Russian Federation) - CSc in of Physical and Mathematical Sciences, Chief Researcher, JSC "Scientific-Research Institute of Polymeric Materials" (16, Chistopolskaya st., Perm, 614113, Russian Federation, e-mail: know_nothing@bk.ru).

Получено 08.05.2017