ЧИСЛЕННОЕ ИССЛЕДОВАНИЕ ДИНАМИКИ ВНУТРИКАМЕРНЫХ ПРОЦЕССОВ В МНОГОСОПЛОВОМ РАКЕТНОМ ДВИГАТЕЛЕ НА ТВЕРДОМ ТОПЛИВЕ ЧАСТЬ 1. МЕТОДИКА РАСЧЁТА Текст научной статьи по специальности «Математика»

УДК 519.9, 629.7

DOI: 10.15593/2224-9982/2023.72.05

М.Ю. Егоров1, Д.М. Егоров1'2, С.М. Егоров1'2, В.И. Белов1,2

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

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

Исследуется динамика переходных внутрикамерных процессов бустерного (разгонного) многосоплового ракетного двигателя на твердом топливе. Метод исследования - постановка вычислительного эксперимента. Рассматриваться сопряжённая постановка задачи, включающая в себя:

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

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

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

- разгерметизацию камеры сгорания и вылет заглушек многосоплового блока ракетного двигателя на твердом топливе (уравнение движения заглушки соплового блока - второй закон Ньютона; метод решения задачи - метод Эйлера).

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

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

M.Yu. Egorov1, D.M. Egorov1,2, S.M. Egorov1,2, V.I. Belov1,2

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

NUMERICAL STUDY OF DYNAMICS IN-CHAMBER PROCESSES IN MULTI-NOZZLE SOLID ROCKET MOTORS. PART 1. CALCULATION METHOD

The dynamics of transient in-chamber processes of booster (overclocking) multi-hop of solid rocket motor is investigated. The research method is the formulation of a computational experiment. The conjugate formulation of the problem is considered, which includes:

- the ignition device is triggered (the combustion rate of the igniting composition is described on the basis of an experimental and theoretical approach that takes into account the burning of combustion products behind the body of the igniting device);

- heating, ignition and subsequent unsteady and turbulent combustion of a charge of a mixed solid fuel (a quasi-homogeneous combustion model is used, formulated on the basis of the equations of thermal conductivity and chemical kinetics

1 Работа была представлена в виде секционного доклада на XXIII Всероссийской научно-техни-

ческой конференции «Аэрокосмическая техника, высокие технологии и инновации - 2022» (Пермь, Пермский национальный исследовательский политехнический университет, 16-18 ноября 2022 г.).

recorded for the condensed phase (solid fuel), taking into account the influence of the gas phase (flare) on the combustion process in the condensed phase; the method of solving the problem is the finite difference method);

- unsteady three-dimensional homogeneous-heterogeneous three-phase flow of nitrogen (pre-pressurization gas of the combustion chamber) and combustion products of the igniting composition and mixed solid fuel in the combustion chamber, multi-core block and behind the multi-core block of the rocket engine (approaches of mechanics of continuous multiphase media are used; the basic system of equations is a system of vortex differential equations of gas dynamics; the solution method is multiparametric a class of difference splitting schemes by physical processes of the Davydov method);

- depressurization of the combustion chamber and the departure of the plugs of the multi-nozzle block of solid rocket motor (the equation of motion of the nozzle block plug is Newton's second law; the method of solving the problem is the Euler method).

Each of the subtasks is considered in relation and resolved simultaneously - at one-time step. To solve the formulated problem, a set of application programs has been developed using (for the main calculation module) the OpenCL multithreading information processing standard. The operability of the software product was checked.

Keywords: numerical study, multi-nozzle solid rocket motors, burning of solid propellant, gas dynamics, calculation method.

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

Рис. 1. Компоновочная схема многосоплового РДТТ

Метод исследования - постановка прямого вычислительного эксперимента. Рассматриваться сопряжённая постановка задачи, включающая в себя:

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

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

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

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

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

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

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

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

V = А0 + А1 ■ X1 + А2 ■ X2 + А ■ X1-X2 + А4 ■ X12 + А5 ■ X22 + А6 ■ X1 • X22;

X1 =

Г F Л — в.

V

\ в

/ В2; X 2 =

Г* >

^ - С1

\ Vв 1 )

/ С2,

(1)

а для нахождения коэффициентов полинома в (1) применить теорию планирования эксперимента. Здесь А0,... А6, В1, В2, С1, С2 — коэффициенты полинома, зависящие от конкретного вида пиротехнического воспламенительного состава и интервалов варьирования конструктивных параметров ВУ. Более подробно методика расчёта изложена в [6].

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

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

Описание процесса нестационарного прогрева, воспламенения и последующего нестационарного и турбулентного горения заряда смесевого твёрдого топлива базируется на модели Мержанова - Дубовицкого - Соркина с учётом, в рамках подхода Горохова - Липанова - Русяка, влияния газовой фазы на процесс горения в конденсированной фазе (к-фазе) [1, 2, 10]. Используется так называемая квазигомогенная модель горения твёрдого топлива. Основное допущение данного подхода - протяжённость прогретого теплового слоя смесевого твёрдого топлива должна быть сопоставима (или больше) с максимальным размером частиц наполнителя смесевого твёрдого топлива.

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

к), ^дв+ф'Ю (2)

Г Е Л

где фк (тк ) = гк ■ ехР —.

\ К0 ■ Тк )

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

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

t = 0, Л > 0, Tk = To, в = 0;

t > 0, Jk = 0, -T), P<P.; (3)

t > 0, Ук =~, Tk = T0, в = 0.

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

t > t., Jk = 0, Tk = T„ ps = P.,

dT f 1 ^

-Xk=qs-ck)pk• vk■ Ts- a — -Pk• vk■ p; (4)

fyk { Pk

t > t., Jk = ~, тк = T0, p=0.

Задача (2)-(4), дополненная замыкающими соотношениями, обезразмеренная и записанная в неравномерной (экспоненциальной) системе координат, интегрируется численным сеточным методом по явным и неявным конечно-разностным схемам [10]. Неявные схемы разрешаются способом прогонки.

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

Т/+1 = At • T++ + Д.; 2 < i < M -1. (5)

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

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

T/+1 = A, - T/+1 + B{; (6)

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

j = Aj. - T/+1 + Br. (7)

Глубина превращения пороха:

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

P/+1 = Pj +Д? • a3 ; 1 <i <M-1; (8)

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

Pi+1 = A*. -Pi+ + B,1 < i <M-1. (9)

Скорость горения заряда смесевого твёрдого топлива vk определяется из (9) с учётом граничных условий Pi+1 = P1 = P,; PM+ = PM = 0 итерационным способом (методом секущих).

В (2)-(9) и далее по тексту приняты обозначения: T - температура, t - время, К - коэффициент температуропроводности, y - координата, Q - тепловой эффект реакции, c - удельная теплоёмкость, в - глубина превращения пороха, Z - предэкспонент, E - энергия активации, R - универсальная газовая постоянная, X - коэффициент теплопроводности, - коэффициент теплоотдачи, q - плотность теплового потока, р - плотность, a - коволюм, p - давление, A, B - коэффициенты прогонки, Дt - шаг по времени, a3 - const, ¥ - безразмерная экспоненциальная функция температуры. Символы: k - конденсированная фаза (твёрдое топливо), g -

газ, 5 - поверхность горения, р - давление, г - номер расчётной точки по координате, у - номер расчётной точки по времени, * - специальное значение.

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

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

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

смеси: а1, а2,а3. Здесь а,. = рг / р" (г = 1...3), где р - истинная плотность фазы. Движение

этих фаз рассматривается как движение взаимопроникающих и взаимодействующих сред.

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

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

^ + Лгу (р^1 ) = ^в + ;

дР дг

+ Лу С*т ) = срт; т = (2,3);

(10)

д(р1ф)

+ Лгу (Дф^) = фв ■ Сев + ф„ ■ С^, ф = ', Ср, цЛ а;

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

+ Лгу(р^Ж ) + адр = —т12 — т13 + и ■ С + и • С ;

д^ 1 и 1 дх х х ев ев е™ е™

^ + Лгу(р1^1 ) + аг| = —ту + у„ Сев + у^;

+ Лгу ( р1 W1) + аг | = — Т + ^ • Сев + ^ • Се™

д(ртит) + Лгу(р и W ) + а •др = Т1т + и С ;

д^ У *т т т / т д^ х рт рт'

д(рт^) + Лгу(р у W ) + а •др = т1т + у С ;

д^ V ^ т т т / т У рт рт'

+ Лгу (рт ™т Wm ) + ат ■ | = Ъ + ' С

рт рт'

т = (2,3);

(11)

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

Э(рт3т )

д( '- + Лу(рт^т ) = д\т + д1т + 3рт • Орт; т = (2,3); (12)

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

Э(р£) + Э(Д) + 3(^1 + ^ ) + Лу (р2£^2 ) + ЛУ (РЗВД )+ (13)

О» ( а1 pWl ) + ( а2рЩ ) + ( азpWз ) = Егв • О8в + Е^ • О^ + Ер2 • Ор2 + Ер3 • Ор

Р 3'

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

л» +ЭМ,

^ Эх Эу Эг (14)

4 = [р/, р1Ф, ри, Р»,р,^, Р, рД, ар]; ¿' = (1-3); / = (3).

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

р = ( - 1)ри •

Е - ^

1 2 ^ У

1

1 - « •рГ

(15)

В (10)-(15) и далее по тексту приняты обозначения: ср — удельная теплоёмкость при постоянном давлении, Е - полная удельная энергия, 3 — внутренняя удельная энергия, О -расходно-приходный комплекс, к — показатель адиабаты, д - функция теплового межфазного взаимодействия, и, v, w — проекции вектора скорости по осям координат 0Х, 0У, 02, W — вектор скорости, — — модуль вектора скорости, х, у, 7 — текущие координаты вдоль осей 0Х, 0У, 0Z, Д — динамическая вязкость, т - функция силового межфазного взаимодействия. Символы: g — газ, к — конвективный, I — лучистый, р — частицы, в - воспламенитель, w - заряд смесе-вого твёрдого топлива, х - вдоль оси 0Х, у - вдоль оси 0У, г - вдоль оси 0Х.

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

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

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

и, = и, — а,

Р¿+0,5, /, к Р 1—0,5, /, к Д

'•,/," ';/,к ';/,к Дх рП

/, к

I = (1-3);

Еп =J" +

mi, j, k щ,,, k

(u ;,j, k )2+(Lk )2.

m

= (2,3);

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

(16)

EE" = E" - V

1,, j ,k 1i, j ,k ¿-U

- V

_ П П _ " "

a ' p i+o 5 ¿k' и - a ' p ,'-0 5 j k • u д t

1 +0,5,j ,k J ^^J ,k 1 +0,5,j ,k ^ -0,5, j ,k J ' U,5,j,k h -0,5,j ,k

( - En )•

\ mi,j,k mi,j,k)

-П П

Pi'

Дх " p1,j,k

" а1 h, j+0,5,k Pi,j+0,5,k —" " V¡ - a¡ h,j +0,5 ,k h ,j -0,5,k D" •v" Pi, j -0,5,k V1,_ j-0 5kk Дt

Ду " p1-j ,k

" ai l¡, j ,k+0,5 Pi,j ,k+0,5 w" - a" h ,j ,k +0,5 h ,j ,k-0,5 • Pi,j,k-0,5 • W1¡jJc-0,5 Дt

Дг

Pi,

где pni+05 jk =—^—^ '+ ; Ü" t =(1 - alfa )u" t + alfa' Ü" k; и т.д. В (16) alfa - сеточный параметр, значение которого изменяется в пределах alfa = 1...3.

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

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

(Pi^i)

i+0,5, j ,k

(1 - beta)-pT + beta • p" • • й" , если и" > 0; (1 - beta ) + beta • p" • • U1 ,если и" < 0; ^ = (1, k, Cp, ц, a, U1, % #1, E1);

(аналогично по другим направлениям (pj^v?!))j+0,5,k, (pj^Wj ),k+0,5); - вторая и третья фаза вдоль оси 0Х:

(17)

(p Ей )

\гш i т ',

если й" > 0;

m,+05, j,k '

р" -Е/ • u

[Pm1+1,j,k ^ Е'+1,j,k

Е = (1, u , v , w , J , E ); m = (2,3);

~ \ ? m? m? m ' m ' m^' V'/'

_ ,k ,jmi +0,5,j ,k '

^j= V •£"+,.* • йи , если Г < 0;

^m1+1,j,k ~<+1, j,k m,+05,jk m,+05, j,k '

2

m=2

1=1

(аналогично по другим направлениям (рm)) 05к, (pm)+05) и т.д. В (17) beta — сеточный параметр, значение которого изменяется в пределах beta = —0,2...0.

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

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

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

(п+1 п (р1и1 ^+0,5,/,к —(р1и1 ^—0,5,/,к д (Р1»1 ))/+0,5,к -(р1»1 ))/ —0,5,к ^ —

р1/ р1/ Дх ' Ду '

(plWl Тим );

\ gв,,/ ,к эт./,к) '

Дг

(18)

(Р тит );+0,5,/,к (( т^ т );-0,5,/,к д, (Рт»т );,/+0,5,к (((тУт );,/-0,5,к д,

Рт/ = Рт/ дх ду —

(р w )п —(р w )п

\гт т ,/ к+0,5 т т ,/

,к—0,5

Дг

•Д, + ^О'Рт,, к ' Д,; т = (2,3);

фп+1 п (р1Фи1 ))+0,5,/,к ~(Р1Фи1 ))0,5,/,к Д, _

^ ф',/,к рп+1 Дх рп+1

Г 1, / ,к Г 1, / ,к

_ (Р1ФУ1 )п,/+0,5,к —(Р1Ф»1 )п,/—0,5,к Д, (Р1Ф^>1 )п,/,к+0,5 ~ ( Р1Ф^>1 )п,/,к —0,5 Д , +

ДУ рпп+1к Дг С

+ (( -о;,.,к +<*/ )); ф=(Ср,дя,а);

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

Рп (р1и1и1 )п — (р1и1и1 )п Д,

^1,/ ,к ^ 1 1> ¿+0,5,/,к м^1 1 1> ¿—0,5, /,к д(

п+1 ~п Г11,/ ,к *>г1 1 1 / ¿+0,5, /,к *>г1 1 1 /1—0,5,/,к

А = и1--/---—-—---:--

■ 1,'/к рп+1 Дх рГ1

Г 1, / ,к г 1, / ,к

(РМ )п,/ +0,5,к —(Р1и1У1 )п,/-0,5,к Д, — (РlU1W1 )п,/,к+0,5 ~(Р1и^1 )п,/,к—0.5 Д,

' ДУ Д2 С

Д? / „ ™ . „ ™ \ Д,

( +т13" )-Д- + ( • О" + и" • ОИ )•

\ х,/,к х,/,к. ! рп+1 \ gвl,/,к gвl,/,к gwl,/,к gwl,/,к )

1,,/,к ^1,,/,к

рп (р и и )п , —(р и и )п , д,

п+1 _ ~п Гт1,/,к \Гт т т )¿+0,5,/,к \Гт т т ) 1-0,5,/,к

"■■/к "-/к рт Дх рт

Г т,/ ,к г т,/ ,к

(р и V )п -(р и V )п д, (р и W )п -(р и W )п д,

УГт т т ,/+0,5,к т т ,/■ -0,5,к Уг»г т т Л,/,к+0,5 ^Гиг т т Л,/,к-0.5

(19)

Ду р:1 Дг Рп

+

,1т» , Д, . п ™ д,

п

+ т1 т--— + ип • ОИ--—; т = (2,3);

х. ■ к п+1 рт; ,к рт; ,к Лп+1 ' 4 '

(аналогично по другим направлениям , w,'"4+1k);

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

Рт■ к )г-+0.5, у к {Рш^т^ш X—0.5,у к А?

ТП+1 _ ТП

т у к = т ■ к ' РП+1 г т. ■

рШС,

{рш^ш \, /+0.5,к \,/-0.5,к А? {рт^т^т \,/ к+0.5 {Рт^т^'т \

■ к-0.5

А?

Ау

А?

+ Г • вП

ь ■ к А?

Аг

+ 1 а1т + а1т'

■ кг у к Ч/г ■ к! лп+1 "Рт,- ,,■ ,к ^РЩ у к пП+1

■; т = (2,3);

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

{ршЕШ - {РшЕш ,■

Р" ш=3

ДП+1 = Д " •■> ,к — Ё

", / к 1, / ,к ЛП+1 ¿-I

Р1., у к т=2

Ё

{р/Е/й/ ) —(рДи )

\Г/ 1 '/г+0,5,■ к \ 1 1 1 ' г—0,5,/,к

А?

Ах РП+1

(р/Д/ С„0Лк -(р/Д<)",у— 0,5,к А?

АУ рг1

(рДл )П/к+0,5 -(РД*/ )Пу к—0,5 А?

+

■+ (20)

Аг рП+1

+ ({ к ,к+■ к „к+Ек ■ к • ^ ■ к+^ „к • }

А?

(21)

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

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

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

/=1

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

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

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

ф/ = Е(-ф,); ЕХ = 1; ф=р,-,к,••• а,а,р,-у,Е; у=1-3- (22)

г г

Для вычисления векторных газодинамических параметров фиктивной ячейки ( = иу,,; у = 1...3); на базе вектора скорости W{и,V,м>} потоковой ячейки и нормали к

граничной поверхности N = , Ny, } строится локальная система координат с направляющими ортами и вычисляются проекции вектора скорости на эти орты (см. рис. 2). Согласно условиям непротекания получаем:

N = {Nx,Ny,N¡}; Wn = и-Ых +у-Ыу + w■Nг;

N = WXN = {^у -, -w■Nx, -и^у};

Wn1 = 0;

N2 = N1 XN ={N-N1-Ny, N,1 N-N,1 N N,1 N -N1

(23)

Wn2 = и^х 2 + 2 + w■Nz 2;

иу = +^ у-К 2)];

/

У/у = +Щ гЪ 2)];

/

WfJ = Е^-Кп^+Wn2 ^)];

г

ЕК= 1; У = 1-з.

г

В (22)-(23) приняты обозначения: V - доля объёма отражённой в поток фиктивной ячейки в объёме расчётной ячейки, Wn, Wn1, Wn2 - проекции вектора скорости W на направляющие орты локальной системы координат. Символы: /- фиктивная ячейка, г - отображённая в поток фиктивная ячейка.

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

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

dw7 , 1 7 г г т,- — = ] р>-] рг

(24)

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

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

w„

< + ■

Аг

т.

Е ((,, к)-Е

г, у ,к г, у ,к

(25)

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

w = wz; дф = 0; ф = р ,к,...X,а, ау,р,3.,Еу; у = 1...3. д,

(26)

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

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

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

Wi,у к+1 Wi,У к-2

п , о

wi■, у к-1 +2 К ;

Ф к+1 = а ■ w 2 + Ь ■ w к

тг, у ,к+1 г, у ,к-2 г ,у ,к-1'

(27)

а + Ь = 1; ф = ру,к,...X,а, ау,р,,Е.; у = 1...3.

В (24)-(27) приняты обозначения: т - масса, а, Ь - весовые коэффициенты. Символы: г -заглушка, 1 - внутренняя сторона, г - внешняя сторона.

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

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

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

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

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

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

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

Наиболее критичной с точки зрения устойчивости и точности численного решения является задача газодинамического течения. Здесь для анализа свойств разностных схем метода Давыдова (см. выражения (16)-(21)) использовался эвристический подход, основанный на рассмотрении параболической формы их дифференциальных приближений [15, 17, 18]. В этом подходе оценивается знак коэффициентов диффузии а^ у диссипативных членов дифференциальных приближений, содержащих частные производные второго порядка по пространственным переменным. Эти коэффициенты обычно группируются в виде матрицы - так называемой матрицы аппроксимационной вязкости. Положительность детерминанта, либо диагональных элементов, либо следа матрицы аппроксимационной вязкости является условием вычислительной устойчивости исследуемой конечно-разностной схемы. Приближая значения этих параметров к нулю - повышаем точность вычислений.

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

Итак, имеем:

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

1 2 Р 2 а + — • -3р-и -р-рр---pr-и • pT -At;

a33 = — • р • u -Ax---р • ux ■ Ax2---(1 - beta ) • u • p x • Ax

1

1

1

.2

2 ^ ^ ^ " p jj

ГДе Ux =3", PpJ

du

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

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

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

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

1. Соркин Р.Е. Газотермодинамика ракетных двигателей на твердом топливе. - М.: Наука,

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

3. Внутренняя баллистика РДТТ / РАРАН; А.В. Алиев, Г.Н. Амарантов и др.; под ред. А.М. Липа-нова и Ю.М. Милёхина. - М.: Машиностроение, 2007. - 504 с.

4. Газовые течения в соплах энергоустановок / под ред. проф. В.Н. Емельянова. - М.: ФИЗМАТЛИТ,

5. Энергетика и внутренняя баллистика ракетных двигателей на твердом топливе / Ю.М. Милёхин, Г.В. Бурский, Г.С. Лавров, В.С. Попов, Д.Н. Садовничий. - М.: Наука, 2018. - 359 с.

6. Численное моделирование внутрикамерных процессов при выходе на режим работы ракетного двигателя твёрдого топлива / Г.Н. Амарантов, М.Ю. Егоров, С.М. Егоров, Д.М. Егоров, В.И. Некрасов // Вычислительная механика сплошных сред. - 2010. - Т. 3, № 3. - С. 5-17.

7. Численное моделирование нестационарных и нелинейных внутрикамерных процессов при срабатывании ракетного двигателя на твёрдом топливе специального назначения. Часть 1. Постановка вычислительного эксперимента / М.Ю. Егоров, С.М. Егоров, Д.М. Егоров, Р.В. Мормуль // Вестник Пермского национального исследовательского политехнического университета. Аэрокосмическая техника. - 2016. - № 47. - С. 53-72.

8. Егоров М.Ю., Егоров Д.М., Егоров С.М. Численное исследование динамики внутрикамерных процессов ракетного двигателя на твёрдом топливе особой компоновочной схемы. Часть 1. Постановка вычислительного эксперимента // Вестник Пермского национального исследовательского политехнического университета. Аэрокосмическая техника. - 2018. - № 53. - С. 63-76.

9. Егоров М.Ю., Егоров Д.М., Егоров С.М. Численное исследование динамики внутрикамерных процессов в маршевом ракетном двигателе на твёрдом топливе с учётом полётных перегрузок. Часть 1. Методика расчёта // Вестник Пермского национального исследовательского политехнического университета. Аэрокосмическая техника. - 2021. - № 64. - С. 91-103.

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

1967. - 368 с.

2017. - 328 с.

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

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

12. Нигматулин Р.И. Динамика многофазных сред: в 2 ч. - М.: Наука, 1987.

13. Давыдов Ю.М. Нестационарный метод расчёта газодинамических задач. - Отчёт ВЦ АН СССР и МФТИ, № 173. - М.: ВЦ АН СССР, 1968. - 29 с.

14. Давыдов Ю.М. Многопараметрические схемы расщепления для решения пространственно-трёхмерных нестационарных задач // Доклады академии наук СССР. - 1979. - Т. 247, № 6. - С. 1346-1350.

15. Давыдов Ю.М. Численный эксперимент в газовой динамике: дис. ... д-ра физ.-мат. наук. - М.: ВЦ АН СССР и МФТИ. - 1981.

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

17. Давыдов Ю.М. Современная нелинейная теория разностных схем газовой динамики. - М.: НИИ Парашютостроения, 1991. - 104 с.

18. Численное исследование актуальных проблем машиностроения и механики сплошных и сыпучих сред методом крупных частиц: 5 т. / под ред. Ю.М. Давыдова. - М.: Национальная Академия прикладных наук, 1995. - 1595 с.

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

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

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

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

23. Programming Guide. AMD Accelerated Parallel Processing OpenCL. Advanced Micro Devices. -2011. - 294 p.

References

1. Sorkin R.E. Gazotermodinamika raketnykh dvigateley na tverdom toplive [Gas-Thermodynamics of solid-fuel rocket engines]. Moscow: Nauka, 1967, 368 p.

2. Lipanov A.M. et al. Chislennyy eksperiment v teorii RDTT [Numerical experiment in the theory of solid propellant motors]. Yekaterinburg: Nauka, 1994, 301 p.

3. A.V. Aliyev and others. Vnutrennyaya ballistika RDTT [Internal ballistics solid propellant motors]. Ed. by A.M. Lipanov and Y.M. Milekhina. Moscow: Mashinostroyeniye, 2007, 504 p.

4. V.N. Yemelyanov. Gazovyye techeniya v soplakh energoustanovok [Gas flows in the nozzles of power plants]. Moscow, FIZMATLIT, 2017, 328 p.

5. Milekhin Yu. M., Bursky G.V., Lavrov G.S., et. al. Energetika i vnutrennyaya ballistika raketnykh dvigateley na tverdom toplive [Power Engineering and internal ballistics of rocket engines on solid fuel]. Moscow, Nauka, 2018, 359 p.

6. Amarantov G.N., Egorov M.Yu., Egorov S.M., Egorov D.M., Nekrasov V.I. Chislennoye modeliro-vaniye vnutrikamernykh protsessov pri vykhode na rezhim raboty raketnogo dvigatelya ^rdogo topliva [Numerical modeling of intra-chamber processes when entering the operating mode of a solid fuel rocket engine]. Computational continuum mechanics, 2010, Vol. 3, no. 3, pp. 5-17.

7. Egorov M.Yu., Egorov S.M., Egorov D.M., et.al. Numerical simulation of non-stationary and nonlinear in-chamber processes when a rocket engine is fired on a special-purpose solid fuel. Part 1. Setting up a computational experiment. PNRPU Aerospace Engineering Bulletin, 2016, no. 47, pp. 53-72.

8. Egorov M.Yu., Egorov D.M., Egorov S.M. Numerical study of the dynamics of in-chamber processes of a rocket engine on solid fuel of a special layout scheme. Part 1. Statement of computational experiment. PNRPU Engineering Aerospace Bulletin, 2018, no. 53, pp. 63-76.

9. Egorov M.Y., Egorov D.M., Egorov S.M. Numerical study of the dynamics of in-chamber processes in solid propellant sustainer taking into account flight overloads. Part 1. Calculation method. PNRPU Engineering Aerospace Bulletin, 2021, no. 64, pp. 91-103.

10. Davydov Yu.M., Egorov M.Yu. Chislennoye modelirovaniye nestatsionarnykh perekhodnykh prot-sessov v aktivnykh i reaktivnykh dvigatelyakh [Numerical modeling of non-stationary transient processes in active and jet engines]. Moscow: NAPN, 1999, 272 p.

11. Rakhmatulin H.A. Osnovy gazovoy dinamiki vzaimopronikayushchikh dvizheniy sploshnykh sred [Fundamentals of gas dynamics of interpenetrating motions of continuous media]. Journal of Applied Mathematics and Mechanics, 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. Davydov Yu.M. Nestatsionarnyy metod ras^ta gazodinamicheskikh zadach [Nonstationary method of calculation of gas-dynamic problems]. - Report of the Central Committee of the USSR Academy of Sciences and MIPT, No. 173, Moscow: Central Committee of the USSR Academy of Sciences, 1968, 29 p.

14. Davydov Yu.M. Mnogoparametricheskiye skhemy rasshchepleniya dlya resheniya prostranstvenno-trekhmernykh nestatsionarnykh zadach [Multiparametric splitting schemes for solving spatial-three-dimensional nonstationary problems]. Reports of the Academy of Sciences of the USSR, 1979, vol. 247, no. 6, pp. 1346-1350.

15. Davydov Yu.M. Chislennyy eksperiment v gazovoy dinamike [Numerical experiment in gas dynamics]. Dissertation for the degree of doctors in physical and mathematical sciences. Moscow, VC of the USSR Academy of Sciences and MIPT, 1981.

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

17. Davydov Yu.M. Sovremennaya nelineynaya teoriya raznostnykh skhem gazovoy dinamiki [Modern nonlinear theory of difference schemes of gas dynamics]. Moscow, Research Institute of Parachute Engineering, 1991, 104 p.

18. Yu.M. Davydov. [Numerical study of actual problems of mechanical engineering and mechanics of continuous and bulk media by the method of large particles]. Moscow, National Academy of Applied Sciences, 1995, 1595 p.

19. Egorov M.Yu. Davydov's method - a modern method of setting up a computational experiment in rocket solid-fuel engine building. PNRPU Engineering Aerospace Bulletin, 2014, no. 37, pp. 6-70.

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

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

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

23. Programming Guide. AMD Accelerated Parallel Processing OpenCL. Advanced Micro Devices, 2011, 294 p.

Об авторах

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

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

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

Белов Василий Игоревич (Пермь, Россия) - старший научный сотрудник расчётного отдела, Научно-исследовательский институт полимерных материалов (Пермь, 614113, ул. Чистопольская, 16, e-mail: vasili.belov.1995@gmail.com).

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., 614990, Perm, e-mail: egorov-m-j@yandex.ru).

Dmitry M. Egorov (Perm, Russian Federation) - CSc in Technical Sciences, General Director, Scientific-Research Institute of Polymeric Materials (16, Chistopolskaya str., 614113, Perm, e-mail: egorovdimitriy@mail.ru).

Sergey M. Egorov (Perm, Russian Federation) - CSc in Physical and Mathematical Sciences, Head of the Calculation Department, Scientific-Research Institute of Polymeric Materials (16, Chistopolskaya str., 614113, Perm, e-mail: know_nothing@bk.ru).

Vasily I. Belov (Perm, Russian Federation) - Senior Researcher, Department of the Calculation, Research Institute of Polymer Materials (16, Chistopolskaya str., 614113, Perm, e-mail: vasili.belov.1995@gmail.com).

Финансирование. Исследование не имело спонсорской поддержки.

Конфликт интересов. Авторы заявляют об отсутствии конфликта интересов.

Вклад авторов. Все авторы сделали равный вклад в подготовку публикации.

Поступила: 15.01.2023

Одобрена: 20.01.2023

Принята к публикации: 28.02.2023

Просьба ссылаться на эту статью в русскоязычных источниках следующим образом: Численное исследование динамики внутрикамерных процессов в многосопловом ракетном двигателе на твердом топливе. Часть 1. Методика расчёта / М.Ю. Егоров, Д.М. Егоров, С.М. Егоров, В.И. Белов // Вестник Пермского национального исследовательского политехнического университета. Аэрокосмическая техника. - 2023. - № 72. - С. 59-74. DOI: 10.15593/22249982/2023.72.05

Please cite this article in English as: Egorov M.Yu., D.M. Egorov, S.M. Egorov, V.I. Belov Numerical study of dynamics in-chamber processes in multi-nozzle solid rocket motors. Part 1. Calculation method. PNRPU Aerospace Engineering Bulletin, 2023, no. 72, pp. 59-74. DOI: 10.15593/2224-9982/2023.72.05