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

Б01: 10.15593/2224-9982/2016.47.03 УДК: 519.9, 629.7

М.Ю. Егоров1, С.М. Егоров2, Д.М. Егоров2, Р.В. Мормуль3

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

Пермь, Россия

3 ПАО «Научно-производственное объединение "Искра"», Пермь, Россия

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

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

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

M.Yu. Egorov1, S.M. Egorov2, D.M. Egorov2, R.V. Mormul3

1 Perm National Research Polytechnic University, Perm, Russian Federation

2 JSC "Research Institute of Polymeric Materials", Perm, Russian Federation

3 PJSC Scientific-Production Association "Iskra", Perm, Russian Federation

NUMERICAL SIMULATION OF UNSTEADY AND NONLINEAR INTRACHAMBER PROCESSES WHEN TRIGGERING ROCKET ENGINE ON SOLID FUEL FOR SPECIAL PURPOSE. PART 1. STATEMENT OF COMPUTATIONAL EXPERIMENT

Numerical modeling of intrachamber processes (so-called holding internal ballistics calculation) is one of the main challenges in the development and design of the rocket engine on solid fuel (SRM). Modern SRM is a complex technical system, which runs at the same time a number of interconnected non-stationary and non-linear physical-chemical processes. SRM for special purpose, considered in this paper, has peculiar functionality and design features. The attempt to optimize the rocket engine parameters was performed by direct numerical simulation of intra-chamber processes. Considered conjugate problem includes: 1) operation of the igniter; 2) heating, ignition and subsequent unsteady and turbulent combustion of solid fuel charge; 3) non-stationary shock wave and vortex homogeneous flow of air and combustion products in the combustion chamber (including igniter chamber), ducts, nozzle block and behind; 4) depressurization of the combustion chamber and the movement of nozzle closure diaphragms. Each of the subproblems is considered with reference to each other and at the same time step. To solve this problem a set of applications on a computer that uses a multi-threaded data processing the OpenCL standard was developed and tested.

Keywords: rocket engine on solid fuel for special purposes, intrachamber processes, combustion, gas dynamics, movement of nozzle closure diaphragm, statement of computational experiment.

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

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

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

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

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

;

Рис. 1. Конфигурация газодинамического тракта РДТТ: 1 - свободный объем камеры сгорания; 2 - центральный газоход; 3 - камера воспламенителя; 4 - боковой газоход; 5 - сопловой блок

- срабатывание воспламенителя;

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

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

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

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

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

Срабатывание воспламенителя

Скорость горения состава заряда воспламенительного устройства определяется по эмпирическим зависимостям вида

vb = Vl • pv1 •(0,9555 + k1-|W|), если p < p1;

vb = v2 • pv2-(0,9555 + kl-|W|), если p > pl. (1)

Приход продуктов сгорания с поверхности горения заряда в камеру воспламенительного устройства определяется по соотношению типа

Св = vB • -ри. (2)

В выражениях (1)-(2) и далее по тексту приняты следующие обозначения: V - скорость горения; p - давление; W - вектор скорости потока; vi, vi, k1, p1 - const (где i = 1, 2); G - приход продуктов сгорания; s - поверхность горения; р - плотность. Индексы: в - воспламенитель; и - истинное значение.

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

Описание процесса нестационарного прогрева, воспламенения и последующего нестационарного и турбулентного горения состава заряда РДТТ базируется на модели Мержанова-Дубовицкого с учетом, в рамках подхода Горохова-Липанова-Русяка, влияния газовой фазы

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

^ = к + у + Я ,ф (Т )

л 2 + Ук л + фк (тк ),

д Эук дУк ск

I = (Тк) ■

дх дУк

(3)

где фк (тк ) = 2к -ехр

Е,,

К0 Тк у

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

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

X = 0, Ук > 0, Тк = То, в = 0;

х > 0, Ук = 0, -Як- ^ = ат-(т, - Т,), р<р,; (4)

Эук *

х > 0, Ук =■

Тк = Т0,

Р = 0.

После воспламенения граничные условия запишутся в виде X > X*, Ук = 0, Тк = Т, ^ = Р*;

дТ

к-^ = ^-(Ср, - ск )-Рк-Ук'Т -

ЭУ,

а--

х > х*, Ук

Тк = Т0,

Рк р = 0.

•рк-^к-р;

(5)

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

оо

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

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

T +1 = A Tj1 +Bi, 2 < i < M-1. (6)

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

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

+1 = A-Tj1+b{, (7)

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

T1j+1 = A* • т/+1+B*. (8)

Глубина превращения твердого топлива:

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

в/+1 =P/ + At-a3 •¥/, 1 < i < M -1; (9)

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

P/+1 = A** •p/+ + B;**, 1 < i < M -1. (10)

Скорость горения твердого топлива vk определяется из формулы

(10), с учетом граничных условий р/+1 =р/ =Р*, PM+1 = PM = 0, итерационным способом (методом секущих).

В выражениях (3)-(10) и далее по тексту приняты следующие обозначения: T - температура; t - время; к - коэффициент температуропроводности; y - координата; Q - тепловой эффект реакции; c -удельная теплоемкость; в - глубина превращения топлива; Z - пред-экспонент; E - энергия активации; R0 - универсальная газовая постоянная; X - коэффициент теплопроводности; ат - коэффициент теплоотдачи; q - плотность теплового потока; a - коволюм; A, B - коэффициенты прогонки; At - шаг по времени; a3 - const; ¥ - безразмерная

экспоненциальная функция температуры. Индексы: k - конденсированная фаза (твердое топливо); g - газ; s - поверхность горения; p -давление; i - номер расчетной точки по координате; / - номер расчетной точки по времени; * - специальное значение.

Приход продуктов сгорания с поверхности горения заряда в камеру сгорания РДТТ также определяется по соотношению типа (2).

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

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

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

- уравнения неразрывности (сохранения массы)

дх+<Иу ^ ) = +Св,

д(рф)+^ ) = +Фв С, (11)

дХ

Ф = к, ср, а;

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

д(ри) + ё1у(риЩ) + др = WXw■Gw + Wxв Св,

+ <Иу (рvW) + Зр = Wyw■Gw + Wyв • Св, (12)

+ (рwW) + др = WZw■Gw + Wzв Св;

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

Щ1^ + ё1у(рЕЩ) + ШУ(рЩ) — Ек-Ок + Ев Ов, (13)

ш

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

^ (^ )=дМ+дМ+дМ,

дх ду дг ^ — [р, ри,..., рЕ, р].

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

р — (к-1)-р.5- 5 — (е(14)

1 - а-р ^ 2 )

Система уравнений (11)—(14), с учетом дополнительных соотношений, интегрируется численно с помощью метода Давыдова (метода крупных частиц) - метода постановки вычислительного эксперимента [6, 7, 10]. Область интегрирования покрывается фиксированной в трехмерном пространстве (эйлеровой) равномерной (однородной и полностью изотропной) расчетной сеткой с ячейками АххАу хАг. Значения целых чисел «г» (вдоль оси 0Х), «}» (вдоль оси 0У) и «к» (вдоль оси 07) обозначают геометрический центр ячейки.

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

рп _ и .

рг+0,5,} ,к рг-0,5, },к А?

пп

Щ,},к — Щ,},к ^п

(аналогично по другим направлениям ■ к, \%П1 к); (15)

Ах Р}

},к, уу(, 1 ,к);

п —п _ п —п А

п _ Т^п рг+0,5,},к'"г+0,5,},к рг-0,5,},к ' "г-0,5,},к А?

Еп — Еп

Ах Рп, },к

рп —п _ рп —п .

рг, 1+0,5,к ч 1+0,5,к рг, }-0,5,к Ч}-0,5,к А? -

АУ Р}

п _п п _п

рЦ,к+0,5 - ™1},к+0,5 - рг,},к-0,5 - ^г,},к-0,5 А?

Аг Р}

Моделирование нестационарных и нелинейных внутрикамерных процессов. Часть 1 nn + nn

n г i, jk -ri+1, j,k —n !л ix \ n . ix ~n

где Pi+0,5,j,k = j 2——; ui,jk =(1 - alfa)• ui,j,k+j,k и т.д.

В выражениях (15) alfa - сеточный параметр, alfa = 3,0.

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

№ )n+0,5, j,k =

(1 - beta )-pn, j л + beta• pn+1, м ]• j л-йГ+одм, если H^j > 0,

i (16) (1 - Ma) • pn+1,j,k + bern • pn, м , если ü^j,k < 0,

£ = (l, fc, cp, a, U, v, , £);

(аналогично для потоков по другим направлениям (p^v)!\,

(p^vv )п+05) и т.д. В выражениях (16) beta - сеточный параметр, beta = -0,15.

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

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

- уравнения неразрывности (сохранения массы)

рп+1 — рп (РЩ)п+0,5, },к -(ри)п-0,5,1 ,к

},к —рп },к АХ А?-

(Р—)п 1 +0,5,к -(Р—)

г,} +0,5,к Уг '¡,}-0,5,к а?

АУ "

(р#)" п-(р>~)п.,,

У1^ >г, 1,к +0,5 У^ >г,},к-I , -А? + ( О

Аг у ¡к '

г,},к +0,5 ,к-0,5-А? + (О" + Овп )-А?,

фп+1 —тп Р1 (р^)п+0,5,1,к -(рФи)п-0,5,1 ,к А? (17)

фг,} к — фг, ---7---п+1— (17)

} 1 рп+к Ах рЩл

(рФ— )п 1+0,5,к -(рФ— )пп 1 -0,5,к А?

АУ рп^^1к

(рФ")пп 1,к +0,5 -(рФ")п 1,к-0,5 А?

Аг р1

■ +

+(Ф"г, 1,-0;, +Ф;М-о;м ^ Ф—( к, , а);

рг, 1,к

- уравнение сохранения импульса по осям координат рп., (рйй)" . ,-(рий)" ., а?

п+1 _ ~п г г,],к У^ /г+0,5,1,к У^ /г-0,5,1,к А?

йй т, — и,-

г, 1 ,к ^ 1 ,к -п+1 Ах -п+1

рг,1,к Ах рг, 1,к

(рй— )п 1+0,5,к -(рй— )п 1 -0,5,к А?

г, 1+0,5,к У~ Н, 1 -0,5,к

АУ рп+1

г,1,к

(рй1~ )п 1,к+0,5 -(рй")пп ;,к-0,5 А?

Аг рп!+1к

+

+(^п -оп + -Оп )--А?-

\ ^ ],к "г, 1 к хвг, ],к в1,1 к ) пп+1

рг, 1 ,к

(аналогично по другим направлениям —п+1к, "гп+1к );

(18)

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

ЕП+1 ] ,к

■ Еп

7 РП,+1к

(рЕй)"+0,5,к ~(рЁй)

г-0,5, ¡,к

Ах

Аг

Рп+1

рг, ],к

(РЕ* 17+0,5,к )

г, ]-0,5,к

Ау

Аг Р«+1

рг, },к

(19)

(рЕ*){+0,5 -(рЁ*)

г, ¡,к-0,5

+(ЕI,,, ,к

Аг

Сп + Еп Сп I

^Ч,к + Ев;,,к ' к )

- +

Аг

Рп+1

рг,} ,к Аг

п+1

г, ],к

Р

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

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

Рис. 2. Постановка граничных условий (схема)

п

п

п

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

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

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

Фа = Е(V 'Ф;), = 1, Ф = Р,к,,а,р,Ё. (20)

I I

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

иа = '(-Мп- Мх + Шг - Мх2)],

I

= \-Ип-Ыу + Шг-Му 2)],

I

= \-Wn-N, + Ш2-2)],

I

= \-Wn-N, + Ш2-Мх 2)],

I

= Е[у; {-Мп-Ыу + МП2-Му 2 )] ,

I

-(-т-м<+шг- м12 )], = 1.

w

а

иа

(21)

w

а

В уравнениях (11)—(21) приняты следующие обозначения: к - показатель адиабаты; и - проекция вектора скорости вдоль оси 0Х; V -проекция вектора скорости вдоль оси 07; w - проекция вектора скорости вдоль оси 0Z (ось 0Z совпадает с осью симметрии РДТТ); И - мо-

дуль вектора скорости; Е - полная удельная энергия; ] - внутренняя удельная энергия; V - объем; ^п, Wn1, ^п2 - проекции вектора скорости W на направляющие оси локальной системы координат (№пг — 0); Ых, Ыу, Ыг - проекции нормали к плоскости отображения в базовой системе координат; Ых2, Ыу2, Ыг2 - проекции оси локальной системы

координат в базовой системе координат. Индексы: " - заряд твердого топлива; х - вдоль оси 0Х; у - вдоль оси 07; г - вдоль оси 02; п - номер шага по времени; а - фиктивная ячейка.

Движение заглушки соплового блока

Поступательное движение каждой заглушки соплового блока РДТТ специального назначения описывается своим уравнением (второй закон Ньютона) вида

г "г

тг'—ГГ — | Р^ - I (22)

аг 0 0

Изменение распределения давления во времени вокруг каждой заглушки берется из задачи газодинамического течения в камере сгорания и газоходах ракетного двигателя.

Уравнение движения (22) интегрируется численно методом Эйлера [8] по следующей конечно-разностной схеме:

— " + А

г1, ],к гг, 1 ,к т

V ( р'п ^ )-V ( ргп .5г )

г, 1,к г, 1 ,к

(23)

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

" — " , дф — 0, ф —(р, к,..., и, V,..., Е, й, V, Е ). (24)

г Эг у '

При конечно-разностной реализации подвижных граничных условий (24) также используются фиктивные ячейки [10] (рис. 3).

Рис. 3. Постановка подвижных граничных условий на внутренней поверхности заглушки соплового блока

Например, для внутренней поверхности заглушки соплового блока записываются следующие разностные соотношения:

< j,k +1 ="«•< j,k-2

jk-1 + 2

ф1 jkk+1 = «•< j k-2 + b^wl jk a + b = 1; ф = (р,..., E).

(25)

В уравнениях (22)-(25) приняты следующие обозначения: m -масса; a, b - весовые коэффициенты. Индексы: г - заглушка; l - внутренняя сторона; r - внешняя сторона.

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

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

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

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

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

Комплекс прикладных программ написан в среде программирования Kdevelop 4.1 для ЭВМ с операционной системой Linux Open SUSE 13.1 x86-64 на алгоритмическом языке C/C++ с использованием (для основного расчетного модуля (main-модуля)) стандарта многопотоковой обработки информации OpenCL [11-14]. Основная идея стандарта

ОреиСЬ состоит в реализации многопотокового выполнения кода по схеме «одна команда - много данных», т.е. одна операция одновременно применяется к большому массиву данных. Подобный подход позволяет в полной мере использовать вычислительный потенциал устройств с множеством относительно простых исполнительных модулей, например таких, как современные графические процессоры (видеокарты), и существенно (на порядок!) повысить производительность вычислений.

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

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

Наиболее критичной с точки зрения устойчивости и точности численного решения является задача газодинамического течения. Здесь для анализа свойств разностных схем метода Давыдова (см. выражения (15)-(21)) используется эвристический подход, основанный на рассмотрении параболической формы их дифференциальных приближений [7, 10]. В этом подходе оценивается знак коэффициентов диффузии а1 у диссипативных членов дифференциального приближения,

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

Согласно трактовке лагранжева этапа метода Давыдова [10] каждая проекция конвективных членов исходной системы дифференциальных уравнений (11)—(13) на координатную ось может аппроксимироваться независимо. Это позволяет, в рамках определенного приближения, при анализе свойств разностной схемы ограничиваться одной координатой. Выпишем для одномерного аналога используемой выше конечно-разностной схемы метода (15)-(21) диагональные элементы матрицы аппроксимационной вязкости.

Будем считать, что поток течет слева направо. Для противоположного направления потока достаточно поменять Ах на - Ах. Итак, имеем

аи = 1'и 'Ах -их ■Ах2 + -1- (рр- и2 )-Аг,

1 л 1 л 2 1 Г 3 , 1 . 2

а22 = — - р - и- Ах---р -их - Ах---1 — Ьега I - и - р х - Ах +

2

2

2 I 2

+ 2-(и-Рл^х + и-Р1 р-Р х + Р1'их )-Ах 2 +

+ — • 2

-3-р-и2-р-Рр-Р'РJ -и2-pJ

Р

- Аг,

(26)

а33 = 1- р-и - Ах-1- р-их •Ах2 -1-(1 -Ьега)-и -рх - Ах2

1 2 1 / \ 2 - 4-Р1'их •Ах - (и-Рл-]х + и-Рр/ -рх )-Ах +

1

+ — 2

-р- и + Р1 •

г Р ^

£-(1 - 2-а!/а )-р

■Аг,

где

ди Эр

и- = & • Рр' = ЗрЭ7 И т'д'

Условие положительности следа матрицы аппроксимационной вязкости (а11 + а22 + а33 )> 0 из уравнений (26) рассматривалось в качестве критерия устойчивости выбранной конечно-разностной схемы метода.

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

Дополнительно результаты расчетов сравнивались с имеющимися данными натурных стендовых испытаний РДТТ специального

1

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

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

Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований (проект № 01-96023/14-16_р_Урал_а).

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

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

2. Численный эксперимент в теории РДТТ / А.М. Липанов,

B.П. Бобрышев, А.В. Алиев [и др.]; под ред. А.М. Липанова. - Екатеринбург: Наука, 1994. - 301 с.

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

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

C. 41-45.

5. Егоров М.Ю., Егоров Д.М. Численное моделирование внутри-камерных процессов в бессопловом РДТТ // XXIII семинар по струйным, отрывным и нестационарным течениям (с междунар. участи-

ем): сб. тр. - Томск: Изд-во Нац. Томск. политехн. ун-та, 2012. -С.124-127.

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

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

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

9. Лойцянский Л.Г. Механика жидкости и газа. - 6-е изд., пере-раб. и доп. - М.: Наука, 1987. - 840 с.

10. Давыдов Ю.М. Крупных частиц метод // Математический энциклопедический словарь. - М.: Сов. энциклопедия, 1988. - С. 303-304

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

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

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

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

References

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

2. Lipanov A.M., Bobryshev V.P., Aliev A.V. [et al.]. Chislennyy ek-speriment v teorii RDTT [Numerical experiment in the theory of solid pro-pellant rocket motors]. Ekaterinburg: Nauka, 1994. 301 p.

3. Egorov M.Yu., Egorov Ya.V., Egorov S.M. Issledovanie neus-toychivosti rabochego protsessa v dvukhkamernom RDTT [Study of working process instability in the two chamber solid propellant rocket engine]. Izvestiya vysshikh uchebnykh zavedeniy. Aviatsionnaya tekhnika, 2007, no. 4, pp. 39-43.

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

5. Egorov M.Yu., Egorov D.M. Chislennoe modelirovanie vnu-trikamernykh protsessov v bessoplovom RDTT [Numerical modeling inter-chamber processes in nozzleless solid propellant rocket motors]. Sbornik trudov XXIII seminara po struynym, otryvnym i nestatsionarnym techeniyam (s mezhdunarodnym uchastiem). Natsionalnyy issledovatelskiy Tomskiy politekhnicheskiy universitet, 2012, pp. 124-127.

6. Egorov M.Yu. Metod Davydova - sovremennyy metod postanov-ki vychislitelnogo eksperimenta v raketnom tverdotoplivnom dvigate-lestroenii [Davydov method - modern method of statement of computing experiment in solid rocket motors industry]. PNRPU Aerospace Engineering Bulletin, 2014, no. 37, pp. 6-70.

7. Davydov Yu.M., Egorov M.Yu. Chislennoe modelirovanie ne-statsionarnykh perekhodnykh protsessov v aktivnykh i reaktivnykh dvigatelyakh [Numerical simulation of unsteady transition processes in active and jet engines]. Moscow: Natsionalnaya akademiya prikladnykh nauk Rossiyskoy Federatsii, 1999. 272 p.

8. Rikhtmayer R.D., Morton Kh. Raznostnye metody resheniya kraevykh zadach [Difference methods for solving boundary value problems]. Moscow: Mir, 1972. 420 p.

9. Loytsyanskiy L.G. Mekhanika zhidkosti i gaza [Fluid mechanics]. Moscow: Nauka, 1987. 840 p.

10. Davydov Yu.M. Kruphykh chastits metod [Large particle method]. Matematicheskiy entsiklopedicheskiy slovar. Moscow: Sovetskaya encyclopedia, 1988, pp. 303-304.

11. Stakhnov A.A. Linux. Saint Peterburg: BKhV-Peterburg, 2004.

912 p.

12. Derk L. Si i Si++. Spravochnik [C and C ++. Directory]. Moscow: VKK, 1997. 592 p.

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

14. Egorov M.Yu., Egorov S.M., Egorov D.M. Primenenie grafich-eskikh uskoriteley dlya povysheniya proizvoditelnosti vychisleniy pri chislennom modelirovanii funktsionirovaniya slozhnykh tekhnicheskikh sis-

tem [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: know_nothing@bk.ru).

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

Мормуль Роман Викторович (Пермь, Россия) - инженер-конструктор 2-й категории ПАО «Научно-производственное объединение "Искра"» (614038, г. Пермь, ул. Академика Веденеева, д. 28).

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).

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

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

Roman V. Mormul (Perm, Russian Federation) - Designer of the 2nd Category, PJSC Scientific-Production Association "Iskra" (28, Academica Vedeneyeva st., Perm, 614038, Russian Federation).

Получено 29.09.2016