Рассматривается численное решение задачи определения параметров течения продуктов сгорания с учетом изменения поверхности горения с течением времени на основном участке работы РДТТ с зарядом осесимметричной формы, имеющем особенность в виде «зонтика». Для данной конфигурации двигателя произведены расчеты параметров внутрикамерных течений продуктов сгорания на основном участке работы с использованием разработанного численного алгоритма. Исследована аппроксимационная сходимость решения задачи как для случая неподвижной поверхности горения топлива, так и для случая движущейся поверхности горения топлива в течение всего времени работы двигателя. Приводятся результаты численного моделирования для разных углов наклона «зонтика».
Numerical simulation of intra-chamber processes in a solid rocket motor with account for burning surface motion.pdf К настоящему времени применительно к отдельным исследуемым процессам (внутренняя баллистика, газодинамика, прочность и др.) разработан достаточно мощный математический аппарат, и на его основе созданы современные пакеты прикладных программ. Тем не менее на практике довольно часто встречаются ситуации, когда параметры РДТТ, полученные в результате натурных испытаний, заметно отличаются от прогнозируемых параметров. Попытки уточнить прогнозируемые характеристики за счет применения численных схем более высокого порядка, более мелких вычислительных сеток, уменьшения шагов интегрирования и т.д. не всегда дают желаемый результат. В этих случаях к положительному результату приводит, как правило, построение математических моделей, учитывающих взаимное влияние процессов различной физической природы, которые формулируются в виде сопряженных задач [1]. Например, для задачи определения параметров течения продуктов сгорания с учетом изменения поверхности горения с течением времени на основном участке работы РДТТ, решение которой рассматривается в данной работе, имеют место, по крайней мере, два взаимовлияющих процесса различной физической природы: горение поверхности топлива и газодинамические процессы движения образовавшихся продуктов сгорания. Как правило, современные РДТТ имеют заряды сложной формы, характеризующиеся трехмерной либо осесимметричной геометрией. Течение газов внутри камеры сгорания может носить сложный характер: возможно наличие пульсаций давления, образование вихревых структур и пр. [2] На данный момент в открытых 1 Исследование выполнено при финансовой поддержке РФФИ в рамках научного проекта № 20-3190033 Моделирование внутрикамерных процессов в ракетном двигателе на твердом топливе 91 источниках преобладают работы, учитывающие сложную геометрию заряда, но предполагающие либо нульмерность, либо одномерность течения продуктов сгорания [1, 3-6]. Для некоторых конфигураций РДТТ такой подход может давать приемлемый результат, однако в общем случае проблематично адекватно описать сложные двумерные или трехмерные процессы течения с помощью нульмерного или одномерного приближения. Учет многомерности течения продуктов сгорания осложняется отслеживанием изменяющегося со временем положения поверхности топлива, которое изначально нельзя предсказать, так как оно в общем случае зависит от параметров течения внутри камеры сгорания. Движение трехмерной и осесимметричной поверхностей может сопровождаться топологическими изменениями, что значительно затрудняет алгоритмы интегрирования. Одним из возможных способов решения данной задачи является неявное представление отслеживаемой поверхности в виде нулевого уровня некоторой функции [7]. Такой способ представления позволяет определять положение любой поверхности, изменяющейся со временем, на достаточно мелкой декартовой сетке. Если неявное задание и отслеживание поверхности расчетной области на декартовой сетке не должно вызывать проблем в реализации, то при дискретизации уравнений, описывающих течение внутри области, могут возникнуть трудности. А именно, граница области не совпадает в общем случае с вычислительной декартовой сеткой, и для точек, находящихся рядом с границей, часть шаблона численной схемы лежит вне области. В данной статье используется численный алгоритм, разработанный авторами [8-10], который позволяет преодолеть вышеописанные трудности с дискретизацией уравнений течения и обладает, в общем случае, произвольным порядком точности, как по пространству, так и по времени. Для неявного представления поверхности используется метод уровней. Для задания значений параметров в фиктивных точках используется обратный метод Лакса - Вендроффа [11-13]. Для интегрирования по времени используется метод Рунге - Кутты [14, 15]. Физическая постановка задачи На рис. 1 приведено схематичное изображение РДТТ, состоящего из цилиндрической камеры сгорания 1 (длиной L и радиусом Rk), прочно скрепленного заряда твердого топлива 2 с прямым центральным каналом радиуса R0 и радиальным «зонтиком» 4, расположенным на расстоянии l от передней стенки камеры сгорания, имеющим высоту H и ширину h. Камера сгорания замыкается соплом 3 с радиусом критического сечения гкр. Рис. 1. Геометрические параметры расчетной области Fig. 1. Geometric parameters of a computational domain А.Е. Кирюшкин, Л.Л. Миньков 92 В начальный момент времени происходит зажигание поверхности твердого топлива, а затем последующее его горение, которое сопровождается образованием газообразных продуктов сгорания, истекающих под действием перепада давления через сопло в окружающую среду. Для упрощения математического описания физических процессов в камере сгорания РДТТ использовались следующие допущения: 1) Зажигание всей поверхности топлива от воспламенителя происходит мгновенно и одновременно. 2) Массовая скорость горения зависит от давления по закону pun = mtpv. 3) Линейная скорость горения равна отношению массовой скорости горения к плотности топлива ub = mtpv/pT. 4) Топливо гомогенное и образующиеся продукты сгорания представляют собой однофазную среду. 5) Образовавшийся газ считается невязким и идеальным в термодинамическом смысле. 6) Поле течения продуктов сгорания обладает осевой симметрией. 7) Время релаксации камеры сгорания много меньше характерного времени горения заряда топлива. Из допущений 2 и 3 видно, что скорость горения в общем случае для каждой точки поверхности горящего топлива различна и зависит от давления, которое может существенно меняться вдоль оси канала. Предположение 7 обеспечивает значительное ускорение времени расчетов при численной реализации, о чем подробнее будет сказано далее. Математическая постановка задачи течения продуктов сгорания Для описания течения внутри камеры сгорания и соплового блока используются уравнения Эйлера, которые с учетом осевой симметрии запишутся как Ut + F(U)x +G(U)y = S(U), (1) тот о где U = (р, pu, pv, pE) ; F(U) = (pu, pu + p, puv, puH) ; G(U) = (pv, puv, pv + p, pvH)T; S(U) = -1/y(pv, puv, pv2, pvH)T; x и y - осевая и радиальная координаты; p, u, v, p - плотность, осевая и радиальная составляющие скорости и давление; E = p/[(k - 1)p] + 0.5(u2 + v2) - полная энергия; H = E + p/p - энтальпия, а k - показатель адиабаты. Начальные условия, записанные в векторной форме, имеют вид (p(0, x, y), p(0, x, y), u(0, x, y), v(0, x, y))T = p,, p0, 0, 0)T. (2) Для данной задачи имеется 4 типа границ, отмеченных на рис. 2: твердая стенка - 1, ось симметрии - 2, выходная граница - 3, поверхность горения - 4. Область, соответствующая топливу, заштрихована. Рис. 2. Схематичное изображение расчетной области Fig. 2. Schematic representation of the computational domain 2 Моделирование внутрикамерных процессов в ракетном двигателе на твердом топливе 93 Граничные условия для каждого типа границы записываются следующим образом: Un = 0; (3) ii О (4) p = p*; (5) pun = mpv; (6) uT = 0; (7) и £ (8) где под un и uT понимается нормальная и тангенциальная составляющие скорости соответственно. На твердой стенке задается условие непротекания (3); на оси симметрии предполагается равенство нулю нормальной составляющей вектора скорости газа (3) и равенство нулю производных по нормали для всех других газодинамических параметров течения (4); на выходной границе в случае дозвукового истечения задается постоянное давление (5), а в случае сверхзвукового истечения граничные условия не ставятся; на поверхности горящего топлива задается приход массы газа (6), предполагается равенство нулю тангенциальной составляющей скорости (7) и задается полная энтальпия продуктов сгорания (8). Математическая постановка задачи изменения формы и геометрии расчетной области Уравнения (1) решаются в области Q, которая ограничена неподвижными границами 1 - 3 и изменяющейся со временем поверхностью горения топлива 4, рис. 2. При этом скорость горения твердого топлива (скорость движения подвижной поверхности) зависит от параметров течения в камере сгорания, которые, в свою очередь, зависят от положения границы твердого топлива. Таким образом, имеется взаимозависимость положения границы области и распределения параметров течения продуктов сгорания внутри области, разделить решение геометрической и газодинамической задач не представляется возможным, и эти задачи следует решать совместно. Рассмотрим уравнение, описывающее положение подвижной поверхности в момент времени t: ф(х, у, t) = 0. (9) Пусть все точки, удовлетворяющие уравнению (9) в некоторый момент времени t, лежат на поверхности рассматриваемой области Q(t). Тогда из условия d ф( х, у, t) d -^1 = 0 , получаемого при условии, что х = x(t) и у = y(t), где под - пониdt dt мается субстанциальная производная, следует дифференциальное уравнение, описывающее эволюцию подвижной границы: дф + дф дх + дфду = 0 dt дх dt ду dt В уравнении (10) производные координат по времени есть не что иное, как скорость движения границы: дх - = u dt ду х dt u у' А.Е. Кирюшкин, Л.Л. Миньков 94 Учитывая, что для задач внутренней баллистики при рассмотрении границ горящего топлива тангенциальная составляющая скорости газа uT равна нулю, то уравнение (10) можно переписать в более удобном виде: дф dt + u„ = 0. (11) Уравнение (11) является разновидностью уравнений Гамильтона - Якоби, которые имеют сходные особенности с гиперболическими уравнениями, и должно быть дополнено начальным условием (12) ф(х, у, 0) = фо(х, у). Скорость перемещения поверхности горения задается в виде u « Іф(х,y,t)=0 \\ub,ф(х,у, t) е Е, [0, ф(х, у, t )гЕ, (13) где Е - поверхность горения топлива в момент времени t, а ub - скорость горения топлива в данной точке. На функцию ф(х, у, t) накладываются определенные ограничения: она должна быть непрерывна и хотя бы один раз дифференцируема. Из этих условий следует, что внутри области решения системы уравнений (1) функция ф(х, у, t) не меняет знак. Можно положить, что она в этой области строго отрицательна. Тогда полная математическая постановка задачи будет состоять из уравнений (1), (9) с начальными условиями для них (2), (12), причем система уравнения (1) решается в области, где ф(х, у, t) < 0. Для системы уравнений (1) граничные условия для различных типов границ имеют вид (3) - (8) при ф(х, у, t) = 0. Скорость перемещения горящей поверхности определяется выражением (13). Численная реализация Поскольку численная реализация полученных уравнений газовой динамики и метода уровней подробно описана в работах [8-10], то далее будут приведены лишь основные этапы алгоритма расчета. Вычислительная область покрывается однородной декартовой сеткой с постоянными шагами Ах = Нх и Ду = Ну, а граница области в общем случае не совпадает ни с одной из координатных линий. Тогда полудискретная аппроксимация уравнений газовой динамики (1) запишется в виде d (Uj )=-Y ( - "-M/2J)-) ((„ - с,,,_„г)+- (М) х у где (х,, у) - точка, лежащая внутри области; FM/2j, Gy+i/2 - численные значения газодинамических потоков. Схематичное изображение области и шаблонов для нахождения потоков изображены на рис. 3. На рис. 3 показано, что часть точек шаблона может выходить за расчетную область. Такие точки будем называть фиктивными. Для корректной работы схемы (14) следует задать определенные значения параметров течения в фиктивных точках. Алгоритм задания значений параметров в фиктивных точках, разработанный С.-В. Шу, называется обратной процедурой Лакса - Вендроффа и подробно описан в работах [11-13]. Он заключается в построении полинома требуемого порядка и модификации получившихся коэффициентов с использованием граничных Моделирование внутрикамерных процессов в ракетном двигателе на твердом топливе 95 условий. Данная процедура является основополагающей в разработанном численном алгоритме, и особенности ее применения к задачам внутренней баллистики детально описаны в [8-10]. Рис. 3. Схематичное изображение границ области, расчетной сетки и шаблонов для нахождения численных значений потоков Fig. 3. Schematic representation of the domain boundaries, computational grid, and templates for numerical calculations of flux values В данной работе потоки газодинамических параметров вычислялись с использованием ENO-схемы третьего порядка точности [16, 17], а в качестве дискретизации уравнений по времени была выбрана схема Рунге - Кутты второго порядка точности, обладающая свойством TVD [14, 15]. Уравнение (11), описывающее изменение положения горящей поверхности, аппроксимировалось с пятым порядком точности по пространству и третьим порядком точности по времени [7]. Для численных расчетов используется следующая последовательность действий: 1) Численно решается система уравнений газовой динамики (1) до установления среднего внутрикамерного давления; 2) Вычисляются значения параметров в «фиктивных» ячейках с помощью обратной процедуры Лакса - Вендроффа; 3) Интегрируется уравнение изменения положения горящей поверхности (11). Время интегрирования уравнений газовой динамики в пункте 1 задавалось пропорционально разности теоретических давлений, вычисленных по формуле Бори для конфигураций на предыдущем и текущем шагах. Такой подход дает значительное ускорение расчетов на основном участке работы РДТТ, так как с течением времени давление внутри камеры сгорания изменяется плавно, что ведет к необходимости использовать малые времена интегрирования уравнений газовой динамики. Стоит заметить, что описанный численный алгоритм позволяет решать задачу и без использования допущения 7 (см. с. 92), интегрируя уравнение (11) каждый раз совместно с уравнениями газовой динамики. Однако такой подход требует намного больше вычислительных ресурсов. В общем случае можно варьировать количество итераций для решения системы уравнений газовой динамики на одну итерацию решения уравнения (11). Иными словами, можно сначала интегрировать уравнения газовой динамики на некоторый шаг по времени At, а затем на тот же шаг интегрировать уравнение поверхности горящего топлива. Такой способ решения мы назовем динамическим. В данной работе будет проведено сравнение результатов динамического подхода с квазистационарным. А.Е. Кирюшкин, Л.Л. Миньков 96 Анализ численных результатов для случая неподвижной поверхности горения Расчеты проводились при следующих значениях параметров: • L = 2.4 м, R0 = 0.2 м, rk = 0.1 м, а Rk, H1, l, h, а - варьируемые параметры. • ѵ = 0.5, mt = 5.3410-3 кг/(м2с-Паѵ), Н0 = 10251150 Дж/кг, • рт = 1700 кг/м3, к = 1.18, ср = 2628 Дж/(кг-К), p0 = р* = 101325 Па, р0 = 1.3 кг/м3. Уравнения (1) обезразмеривались с использованием в качестве характерных следующих масштабных величин: • Um к k ^ - масштаб скорости; • Pm = -Г t m 1 v - масштаб давления, где ST - площадь горящей поверх- Іг(к) Ър ) ности в начальный момент времени, Ъкр - площадь критического сечения сопла, а 2 \\ к+1 2 ' к-1 - газодинамический комплекс; г(к )=Шд • Pm = Pp - масштаб плотности; Um • xm = 1 м - масштаб длины; • t = - m U m масштаб времени. Дальнейшие результаты приведены для безразмерных параметров течения. Вычислительная сетка для всех приведенных результатов расчетов имеет одинаковые безразмерные пространственные шаги по пространству вдоль двух координатных осей h = hx = hy. Для проверки сходимости численного решения получим стационарные решения системы уравнений (1) с разными шагами по пространству, соответствующие трем промежуточным формам поверхности заряда с параметрами Rk = 0.7 м, Hi = 0.4 м, l = 0.7 м, h = 0.1 м, а = п/3, возникающим в процессе его прогорания и показанным на рис. 4. Рис. 4. Границы расчетной области: 1 - начальный момент времени, 2 -промежуточный момент времени, 3 - заключительный момент времени Fig. 4. Computational domain boundaries: 1, initial time instant; 2, intermediate time instant; and 3, final time instant Моделирование внутрикамерных процессов в ракетном двигателе на твердом топливе 97 На рис. 5 представлены зависимости осредненного по объему внутрикамерно-го давления в зависимости от обратного значения шага вычислительной сетки для трех геометрических конфигураций, изображенных на рис. 4. Из полученных результатов видно, что значение осредненного внутрикамерного давления перестает зависеть от обратного значения шага вычислительной сетки с точностью до 1% для значения 1/h = 375, или для шага сетки h = 0.0027, для всех трех геометрических конфигураций. Полученное с помощью численной схемы среднее давление достаточно близко к давлению, рассчитанному по формуле Бори (горизонтальные линии на рис. 5). Расхождение в значениях давления можно объяснить тем, что в формуле Бори коэффициент расхода сопла задавался равным единице. Р Р a 1.08 1.68 1.64 1.04 * • • • 1.6 1 1.56 0 250 500 750 1/h* b Р 1.36 1.32 1.28 0 250 500 750 1/h* 0 250 500 750 1/h* Рис. 5. Влияние обратной величины шага сетки 1/h на среднее внутрикамерное давление для: a - положения 1 на рис. 4, b - положения 2 на рис. 4, c - положения 3 на рис. 4 Fig. 5. Impact of a grid step reciprocal value 1/h on the average intra-chamber pressure for (a) location 1 (fig. 4), (b) location 2 (fig. 4), and (c) location 3 (fig. 4) На рис. 6 показаны распределения давления и числа Маха внутри камеры сгорания РДТТ для различных значений шага вычислительной сетки, соответствующие начальной форме заряда. Из рис. 6 можно заключить, что, начиная с шага по пространству h = 0.002, происходит образование вихревых структур и зон сжатия - 1 1 04 1 08 1 12 1 16 О 0.1 0.2 0 3 04 Рис. 6. Влияние шага расчетной сетки h на распределение давления (a) и числа Маха (b) внутри камеры сгорания Fig. 6. Impact of the computational grid step h on the (a) pressure and (b) Mach number distribution inside a combustion chamber А.Е. Кирюшкин, Л.Л. Миньков 98 разрежения, о чем свидетельствуют замкнутые изолинии поля чисел Маха и давления. Образованные вихревые структуры для начальной конфигурации имеют нестационарный характер, не исчезают с течением времени, а периодически появляются у изголовья «зонтика» и исчезают при истечении продуктов сгорания через сопло. Формирование вихревых структур можно объяснить тем, что при истечении струи газа из «зонтика» в основной канал происходит ее поджатие основным потоком к поверхности топлива, что приводит к формированию условий обтекания угла с острой кромкой, вызывающих, вследствие схемной вязкости, хоть и незначительной, возникновение вихрей. Подобное явление в условиях турбулентного течения описано в работе [2]. Стоит, однако, отметить, что нестационарный характер течения не влияет на среднее внутрикамерное давление. Таким образом, из вышеприведенных результатов видно, что численное решение для рассмотренных конфигураций, во-первых, сходится при измельчении сетки и, во-вторых, полученные значения давления хорошо согласуются с давлением, полученным по формуле Бори. Анализ численных результатов моделирования внутрикамерных процессов РДТТ на основном участке работы с подвижной поверхностью горения Рассмотрим следующую задачу: по заданным геометрическим характеристикам заряда топлива определить зависимость среднего давления внутри камеры сгорания от времени. Для того чтобы убедиться в достоверности численных результатов, определим шаг по пространству, при котором решение перестает от него зависеть, и в дальнейшем все результаты численного моделирования будут приведены при соответствующем шаге по пространству. Рассмотрим область со следующими геометрическими параметрами: Rk = 0.7 м, H = 0.4 м, l = 0.7 м, h = 0.1 м, а = п/3. На рис. 7 показано изменение осредненного давления в камере сгорания от времени. Решение, полученное для шага h = 0.0027, отличается от решения, полученного для шага h = 0.001, в метрике L2 на 1.8%. Fig. 7. Pressure as a function of time: h* = (1) 0.004, (2) 0.0027, (3) 0.002, (4) 0.0013, and (5) 0.001 Моделирование внутрикамерных процессов в ракетном двигателе на твердом топливе 99 В силу незначительного различия в численных решениях, полученных на разных шагах по пространству, и в силу того, что уменьшение шага разностной сетки по пространству в 2.7 раза ведет к увеличению времени расчета почти в 20 раз, можно заключить, что такой шаг по пространству h = 0.0027 является оптимальным для расчетов зависимости среднего давления внутри камеры сгорания от времени, и дальнейшие расчеты выполнялись при таком шаге. Время расчета задачи внутренней баллистики РДТТ на всем участке его работы с использованием технологии CUDA занимает приблизительно 1 час при расчете на видеокарте NVIDIA GeForce GTX 1080. Сравним зависимости среднего давления в камере сгорания, полученные с помощью квазистационарного и динамического подходов, которые представлены на рис. 8. Время расчета для динамического подхода составляет приблизительно 10 часов, что в 10 раз медленнее квазистационарного подхода. Относительная ошибка в метрике L2 для квазистационарного подхода составляет 1.3% по сравнению с динамическим. Из полученных результатов можно сделать вывод, что промежутки времени, которые берутся для квазистационарного подхода, в случае падения давления иногда оказываются недостаточно большими, однако разница в зависимости среднего давления от времени незначительная, поэтому в дальнейшем для расчетов использовался квазистационарный подход. Рис. 8. Зависимость давления от времени: 1 - динамический подход, 2 - квазистационарный подход Fig. 8. Pressure as a function of time: 1, dynamic approach and 2, quasi-stationary approach Результаты численного моделирования для различных геометрических форм топливных зарядов РДТТ Рассмотрим первую конфигурацию заряда, которая определяется следующими геометрическими параметрами: Rk = 0.9 м, H = 0.55 м, l = 0.5 м, h = 0.08 м, а = п/2. Для данной конфигурации «зонтик» расположен под прямым углом к оси симметрии. Как видно из представленных результатов на рис. 9, давление в начальный момент времени имеет достаточно сложное распределение внутри камеры сгорания с большими градиентами и вихревыми структурами. С течением времени А.Е. Кирюшкин, Л.Л. Миньков 100 вихревые структуры исчезают и градиент давления уменьшается. Когда все топливо на левой стенке выгорело, можно считать, что давление всюду в камере сгорания имеет постоянное значение. 0.91 0.935 0.96 0.985 1.01 1.035 1.55 1.56 1.57 1.58 1.59 1.6 1.84 1.86 1.88 1.9 1.05 1.065 1.08 Рис. 9. Распределение давления в камере сгорания: а - в начальный момент времени, b - в момент времени, когда топливо еще «не дошло» до верхней стенки, c - в момент времени, когда часть топлива на верхней стенке выгорела, d - момент времени, когда все топливо на левой стенке выгорело Fig. 9. Pressure distribution inside a combustion chamber: (a) at the initial time instant, (b) at a time instant when the propellant has not reached the top wall yet, (c) at a time instant when the propellant on the top wall has partly burned out, and (d) at a time instant when all the propellant on the left wall has burned out На рис. 10, а изображены линии поверхности твердого топлива через равные промежутки времени At = 8336. Из рис. 10, а можно увидеть, что топливо горит практически параллельными слоями с изменяющейся со временем скоростью горения, так как расстояния между линиями поверхности различно для одного и того же At. Горение почти параллельными слоями объясняется тем, что давление вдоль канала меняется не более чем на 20 %, и зависимость скорости от давления носит степенной характер со степенью ѵ < 1. На рис. 10, b представлены зависимости давления и площади горящей поверхности от времени. На данном графике Рис. 10. Эволюция горящей поверхности (а) и зависимость давления (1) и площади горящей поверхности (2) от времени (b) для заряда РДТТ с зонтиком, расположенным под прямым углом Fig. 10. (a) Evolution of a burning surface and (b) time dependences of (1) pressure and (2) burning surface area for SRM at a right inclination angle Моделирование внутрикамерных процессов в ракетном двигателе на твердом топливе 101 значительным уменьшениям площади горящей поверхности, а следовательно, и резким спадам давления, соответствует выгорание топлива у верхней стенки для t ~ 22000 и на левой стенке для t ~ 67000. В общем случае для всех зарядов твердого топлива можно наблюдать следующие особенности: по мере выгорания топлива и увеличения объема камеры распределение давление носит все более равномерный характер; несмотря на достаточно сильные градиенты давления в камере сгорания в начальный период работы РДТТ прогар топлива происходит параллельными слоями; резкие спады среднего по камере давления с течением времени соответствуют выгоранию топлива либо на верхней, либо на левой стенке. В остальные моменты времени изменение давления носит плавный характер. На рис. 11 представлены результаты расчетов эволюция поверхности (рис. 11, а) через равные промежутки времени At = 10000 и зависимости среднего внутрикамерного давления и площади горящей поверхности от времени (рис. 11, b). Для данной конфигурации топливо на левой границе выгорает одновременно с топливом на правой границе, поэтому в графике зависимости давления от времени имеется только один резкий спад давления при выгорании топлива на верхней стенке. Изолинии поверхности горения имеют схожие черты с формой канала заряда РДТТ, изображенной на рис. 10, а: горение происходит параллельными слоями, но расстояние между уровнями также различается для разных промежутков времени. Рассмотрим вторую конфигурацию заряда со следующими параметрами: Rk = 0.7 м, H = 0.4 м, l = 0.5 м, h = 0.1 м, а = 2п/3, в которой угол наклона «зонтика» острый. Рис. 11. Эволюция горящей поверхности (а) и зависимость давления (1) и площади горящей поверхности (2) от времени (b) для заряда РДТТ с «зонтиком», ориентированным под острым углом к оси абсцисс Fig. 11. (a) Evolution of a burning surface and (b) time dependences of (1) pressure and (2) burning surface area for SRM at an acute inclination angle to the axis of abscissas Наконец, рассмотрим третью конфигурацию, когда «зонтик» составляет острый угол с отрицательным направлением оси х. Форма шашки определяется следующими параметрами: Rk = 0.7 м, H = 0.35 м, l = 0.7 м, h = 0.2 м, а = п/6. На рис. 12 представлены результаты расчетов эволюции поверхности (рис. 12, a) через равные промежутки времени At = 10000 и зависимости среднего внутрикамерного давления и площади горящей поверхности от времени (рис. 12, b). Особенность данной конфигурации заключается в том, что топливо на верхней и левой стенках корпуса сгорает практически одновременно, что способствует А.Е. Кирюшкин, Л.Л. Миньков 102 сильному падению давления для t я 22000. Горение поверхности так же, как и для двух предыдущих конфигураций, происходит параллельными слоями с изменяющейся со временем скоростью горения. Рис. 12. Эволюция горящей поверхности (а) и зависимость давления (1) и площади горящей поверхности (2) от времени (b) для заряда РДТТ с «зонтиком», ориентированным под тупым углом к оси абсцисс Fig. 12. (a) Evolution of a burning surface and (b) time dependences of (1) pressure and (2) burning surface area for SRM at an obtuse inclination angle to the axis of abscissas Заключение В данной работе решена сопряженная задача определения параметров течения продуктов сгорания в проточном тракте РДТТ совместно с определением положения изменяющейся во времени горящей поверхности твердого топлива в осесимметричной постановке с помощью предложенного численного алгоритма. Численное решение было исследовано на аппроксимационную сходимость для случая неподвижной поверхности горящего топлива и для случая эволюционирующей поверхности горения во время работы двигателя РДТТ на основном участке его работы. Проведенные сравнения осредненного давления в камере сгорания РДТТ с оценками по формуле Бори показали, что результаты численного моделирования достоверно описывают изменение давления газа в камере сгорания РДТТ во времени совместно с эволюционирующей поверхностью горения заряда твердого топлива. В качестве примера разработанный численный алгоритм был применен для моделирования работы РДТТ на его основном участке для зарядов с «зонтиком», ориентированным под различными углами наклона к оси симметрии.
Милехин Ю.М., Ключников А.Н., Попов В.С. Сопряженная задача моделирования внутрибаллистических характеристик бессопловых РДТТ // Физика горения и взрыва. 2013. Т. 49. № 5. С. 77-85.
Глазунов А.А., Еремин И.В., Жильцов К.Н., Костюшин К.В., Тырышкин И.М., Шувариков В.А. Численное исследование определения величин пульсаций давления и собственных акустических частот в камерах сгорания с наполнителем сложной формы // Вестник Томского государственного университета. Математика и механика. 2018. № 53. С. 59-72. DOI: 10.17223/19988621/53/7.
Cavallini E. Modeling and Numerical Simulation of Solid Rocket Motors Internal Ballistics: PhD thesis. 2010. URL: https://core.ac.uk/download/pdf/74323997.pdf
Sultwald W. Grain regression analysis: Master’s thesis. 2014. URL: https://scholar.sun.ac.za/bitstream/handle/10019.1/86526/sullwald_grain_2014.pdf
Tshokotsha M.H. Internal Ballistic Modelling of Solid Rocket Motors Using Level Set Methods for Simulating Grain Burnback: Master’s thesis. 2016. URL: https://pdfs.semanticscholar.org/d0c7/5902ebacf32fc3c60e57158a9e040b9154f8.pdf
Lorente A.P. Study of Grain Burnback and Performance of Solid Rocket Motors: PhD Thesis. 2013. URL: https://upcommons.upc.edu/bitstream/handle/2099.1/17700/Memoria_Arnau_Pons_Lorente.pdf
Osher S., Fedkiw R. Level Set Methods and Dynamic Implicit Surfaces. New York: Springer, 2003. DOI: 10.1007/b98879.
Кирюшкин А.Е., Миньков Л.Л. Численное решение двумерных уравнений газовой динамики с подвижными границами на неподвижной вычислительной сетке на примере задач внутренней баллистики РДТТ // Всероссийская молодежная научная конференция «Все грани математики и механики»: сб. статей. Томск: Издательский Дом ТГУ, 2017. С. 168-176.
Кирюшкин А.Е., Миньков Л.Л. Численное решение задачи внутренней баллистики РДТТ на всем участке работы для зарядов сложной формы // V Всероссийская научнотехническая конференция с международным участием «Актуальные проблемы ракетно-космической техники» («V Козловские чтения»): сб. материалов, 2017. Т. 2. С. 55-64.
Kiryushkin A.E., Minkov L.L. Solution of internal ballistic problem for SRM with grain of complex shape during main firing phase // Journal of Physics: Conference series. 2017. V. 894(1). P. 012041-1-012041-7. DOI: 10.1088/1742-6596/894/1/012041.
Tan S., Shu C.-W. Inverse Lax-Wendroff procedure for numerical boundary conditions of conservation laws // Journal of Computational Physics. 2010. V. 229(26). P. 8144-8166. DOI: 10.1016/j.jcp.2010.07.014.
Tan S., Shu C.-W., Ning J. Efficient Implementation of high order inverse Lax-Wendroff boundary treatment for conservation laws // Journal of Computational Physics. 2012. V. 231(72). P. 2510-2527. DOI: 10.1016/j.jcp.2011.11.037.
Tan S., Shu C.-W. Inverse Lax-Wendroff Procedure for numerical boundary conditions of hyperbolic equations: Survey and new developments // Advances in Applied Mathematics, Modeling, and Computational Science. 2013. V. 66. P. 41-63. DOI: 10.1007/978-1-4614-5389-5_3.
Gottlieb S., Shu C.-W. Total variation diminishing Runge-Kutta schemes // Mathematics of Computation. 1998. V. 67. P. 73-85. DOI:10.1090/S0025-5718-98-00913-2.
Shu C.-W. Total-variation-diminishing time discretizations // SIAM Journal on Scientific and Statistical Computing. 1988. V. 9. P. 1073-1084. DOI: 10.1137/0909073.
Harten A., Engquist B., Osher S., Chakravarthy S.R. Uniformly high-order accurate nonoscillatory schemes // Journal of Computational Physics. 1987. V. 71(2). P. 231-303. DOI: 10.1137/0724022.
Shu C.-W., Osher S. Efficient implementation of essentially non-oscillatory shock capturing schemes // Journal of Computational Physics. 1989. V. 83(1). P. 32-78. dOi: 10.1016/0021-9991(89)90222-2.