Введение
Вычислительные алгоритмы, разработанные авторами, были реализованы в виде набора программ, включая подсистемы, предназначенные для установки начальных данных и определения термодинамических свойств веществ и механизмов химических превращений, а также систему визуализации результатов расчетов, которая не имеет аналогов. Оригинальная система визуализации позволяет отслеживать распределение параметров потока не только в фиксированные моменты времени, но и вдоль произвольных сеток линий, число которых может изменяться в процессе расчета.
Предложено использовать самоподобную проблему Г. М. Бам-Зеликовича в качестве тестовой задачи при моделировании одномерных и многомерных потоков многокомпонентного детонирующего газа за инцидентными и отраженными ударными волнами. Решение системы уравнений с частными производными, описывающей поток детонирующего газа, в среднем асимптотически стремится к решению проблемы Г. М. Бам-Зеликовича, что позволяет мгновенно определить параметры, характеризующие плоские детонационные волны до проведения многомерных расчетов.
Математическая модель
Здесь, щявляется временем, является пространственной координатой, , , , и являются плотностью, скоростью, давлением и удельной энтальпией смеси соответственно, = () – это площадь канала, — молярные концентрации массы -го химического компонента, – это количество учитываемых компонентов, – скорость изменения концентраций -го компонента на объем единицы смеси и – это скорость звука. Индексы ,P и обозначают частную дифференциацию соответствующего параметра.
Здесь, – это скорость SW (для CD, = u1 = u2), индекс 1 обозначает параметры до волны, а индекс 2 обозначает параметры после волны. При прохождении через ударную волну в неэквилибриальном потоке соотношения дополняются условиями постоянства концентрации:
Здесь, – универсальная газовая постоянная, 0 = 101,325 Па – стандартное давление, – стехиометрические коэффициенты, а – символ молекулы -го вещества. Число элементарных реакций в единицу времени в единице объема зависит от температуры, давления и объемных концентраций реагентов:
Для аппроксимации зависимости скоростных констант прямых реакций от температуры и давления используется обобщенная формула Аррениуса:
В случае равновесного процесса химических превращений концентрации являются неявно заданными через функции температуры и давления и могут быть рассчитаны путем решения уравнений химического равновесия и сохранения элементного состава:
Здесь, – количество элементов в смеси; – матрица состава; и – неизвестные значения.
Давайте кратко опишем особенности предложенной техники моделирования.
Определение шага интегрирования из условия Куранта, обеспечивающего стабильность и контроль ошибки при аппроксимации уравнений кинетики химических реакций;
Определение координат подвижных узлов на новом слое и анализ наличия пересечений линий разрыва между собой и траекторий газа. Если пересечение присутствует, то время самого раннего пересечения становится временем нового слоя, а шаг интеграции пересчитывается;
Проводится одна итерация для определения параметров газового поля и концентраций химических компонентов на всех узлах сетки;
Проверяется точность выполнения разностных уравнений: если заданная точность достигнута во всех точках разностной сетки, то переходим к шагу (е), а если хотя бы в одной точке она не достигнута, то выполняется другая итерация, что означает возвращение к шагу (б);
Если было установлено пересечение разрывов (см. шаг b), то в зависимости от физической природы взаимодействующих разрывов определяется структура материализованного потока в точке пересечения (например, если пересекаются два сильных разрыва, то задача Римана должна быть решена);
Переход к расчету для следующего временного слоя.
Передача параметра
Алгоритм передачи параметра (Т) является наиболее распространенным алгоритмом для взаимодействия сеточных линий, который используется в двух версиях. Первый случай – когда сеточные линии не являются линиями разрыва для параметров потока (рисунок 2а). Следует отметить, что в этом описании и далее координаты (время и пространственная координата) всех точек до и после взаимодействия равны, в то время как скорости движения сеточных линий разные. В нижней части находятся рассчитанные точки до обработки взаимодействия, а в верхней части – точки после обработки. Сеточные линии (1, 2) и (2, 1), указанные жирными стрелками, меняются местами, в то время как в некоторых случаях, чтобы избежать округлительных ошибок, ряд параметров в точках 1 и 2 сохраняются. Например, если линия (1-2) является траекторией газа, а линия (2-1) характеристикой (C−), то концентрации химических компонентов в точках 1 и 2 после обработки пересечения точно равны концентрациям химических компонентов в точке 1 до обработки пересечения.
Второй случай
Второй случай (рисунок 2b) – когда одна из сеточных линий является линией разрыва для параметров потока. Например, линия (1-2, 1-2) принадлежит семейству SW C+, а линия (2-1) – траектория газа; в этом случае параметры потока в точках 1 и 2 после обработки пересечения совпадают, в то время как скорость, давление, температура и другие параметры вдоль траектории резко изменяются от параметров в точке 2 к параметрам в точке 1 (включая комплексы, входящие в правые части характеристических соотношений), а концентрации химических компонентов остаются неизменными.
Взаимодействие ударной волны
При взаимодействии ударной волны с характеристической линией обработка пересечения зависит от семей, к которым они относятся. Например, конфигурация, показанная на рисунке 2b, может материализоваться, когда ударная волна SW (1-2, 1-2) и характеристика (1-2) принадлежат семейству C+. В этом случае используются два вида алгоритмов обработки. Первый – передача характеристики в точку 1 слева от ударной волны, в этом случае характеристика меняет свой тип с C+ на C− (нет характеристики C+ слева от SW, направленной вперед по времени), и, соответственно, изменяется скорость ее распространения. Второй – удаление характеристики C+. Линия 1-2 на рисунке 2b может составлять характеристику C−. В этом случае выполняется стандартная передача параметра, описанная выше.
Пересечение с свободной границей
На рисунке 2c показан случай пересечения сеточной линии с краем свободной границы и сеточной линии (GT, PC−). При пересечении сеточной линии с границей вычислительной области сеточная линия удаляется, и некоторые параметры удаленной сеточной линии передаются в точку границы.
Markdown table example
Код таблицы:
| Заголовок 1 | Заголовок 2 |
|-------------|-------------|
| Строка 1 | Строка 2 |
| Строка 3 | Строка 4 |
Результат:
| Заголовок 1 | Заголовок 2 |
|---|---|
| Строка 1 | Строка 2 |
| Строка 3 | Строка 4 |
Когда ударная волна (или контактное разрывное значение) взаимодействует с свободной границей (рисунок 2d, SW C+ или CR (1, 1′) и RHB (2-2) линия 1 после обработки—PC−), SW (CD) удаляются (R), и на новом временном слое появляется новая точка расчета: C− характеристика в случае правой свободной границы или C+ характеристика в случае левой свободной границы.
Рисунок 2е демонстрирует случай отражения SW+ от RHW. После отражения SW становится волной семейства C− (SW−), в то время как параметры потока в точке 1 остаются неизменными (за исключением скорости сетки), и параметры в точках 1′ и 2, вместе со скоростью SW, рассчитываются путем решения задачи отражения SW от стены. В случае отражения характеристики характеристика меняет свой тип (из C+ в C− или из C− в C+), а скорость движения, в то время как остальные параметры остаются неизменными.
Следует отметить, что метод моделирования предусматривает случаи дегенерации ударных волн и контактных разрывностей (рисунок 2h), если разница в параметрах, наблюдаемых на них, менее 0.001%. В этом случае ударные волны дегенерируют в характеристики, а контактные разрывности дегенерируют в траектории газа.
Случай пересечения характеристик, принадлежащих одной семье, который недопустим с точки зрения применяемой системы уравнений в частных производных, используется в численной реализации для описания возникновения ударных волн в изначально бесударном потоке (рисунок 2i). Здесь 1 и 2 являются характеристиками семьи C+. Параметры потока, полученные в точках 1 и 2 в результате решения системы разностных уравнений, различаются. Для разрешения кинематически несовместимой ситуации, возникшей, необходимо решить задачу Римана. В показанном случае (рисунок 2i) образуется ударная волна, распространяющаяся вправо SW+ (4, 4′), контактный разрыв (3, 3′), и слабая волна разряжения веера, описываемая двумя характеристиками 1 и 2. Следует отметить, что в реальных проблемах, моделированных через подробные разностные сетки, разница в параметрах на веере слабой волны и на контактной разрывности часто находится ниже порога дегенерации, и они не учитываются при расчете.
В принципе, задача Римана может быть использована для разрешения любого из взаимодействий между внутренними линиями сетки из таблицы, которые отмечены знаком -. Однако, для избежания случаев потери газовой массы, расположенной между изначально определенными линиями сетки (траекториями газа), если обнаружено пересечение между траекториями и контактными разрывностями, алгоритм прерывается с соответствующей диагностикой.
In SW + PW (point on the wall) interactions, the parameters of the reflected shock wave are determined via Rankine–Hugoniot type relations for a multicomponent gas using the condition that the velocity behind it is equal to zero, supplemented, if necessary, by the conditions of thermodynamic equilibrium.
Results
As a result of numerical simulation, the flow structure presented in Figure 1 was obtained. Here, the red lines indicate the trajectories of shock waves, and the green lines indicate contact discontinuities. During the observation time, which was ~0.2 s, the initial shock wave managed to be reflected four times from the ends of the installation (Figure 3a). Figure 3b shows an enlarged fragment of a complex shock wave structure that appears when the incident shock wave is first reflected from the wall and interacts with the initial contact discontinuity. The flow structure includes a set of secondary shock waves and contact discontinuities of varying intensity. It should be noted that, in contrast to the shock-capturing method, all generated shock waves and contact discontinuities are clearly distinguished; they are double computational grid nodes, which are interconnected by the exact fulfillment of the Rankine–Hugoniot relations. When these computational nodes interact with each other, the Riemann problem is solved with high accuracy. The shock waves and contact discontinuities formed as a result of the solution become new grid nodes.
The numerical modeling resulted in the flow structure presented in Figure 3. Here, “thick solid” lines denote SWs, “dotted” lines denote CDs, and “thin solid” lines denote characteristics. During the observation time, which was ~0.2 s, the initial SW managed to reflect four times from the ends of the setup.
Figure 4 provides a comparison of pressure sensor readings typical for this series of experiments with the results of the numerical modeling. Figure 4a demonstrates the readings of the pressure sensors located at the HPC end (the upper curve) and in the middle of the LPC (the lower curve). Figure 4b,c demonstrate the results of numerical modeling of the flow in the installation: the graph of the pressure at the LHW, which corresponds to the readings of the HPC end sensor, and the graph of the pressure in the middle of the LPC. The modeling technique under consideration makes it possible to “take” the readings of the sensors. To perform this, fixed nodes with coordinates corresponding to the location of the pressure sensors along the setup length are introduced into the computational grid.
The temperature behind the reflected SW exceeds the self-ignition temperature of the combustible mixture, and ignition occurs at the left wall (Figure 7, curve 1). An increase in temperature and pressure leads to the appearance of a compression wave, which, at t ~ 0.00245, forms a “hanging” shock wave (Figure 5 and Figure 6). It should be noted that the temperature drop at the hanging SW at the moment of its origination depends on the computational grid; in the calculation, it is less than 100 K (Figure 7, curves 2, 3). At subsequent time stamps, the hanging SW intensity increases, and a detonation wave is formed (Figure 7), which propagates behind the reflected shock wave. At t ~ 0.0031, the waves interact; an overcompressed detonation wave, a contact discontinuity, and a fan of expansion waves are generated (Figure 5). Over time, the intensity of the detonation wave decreases, and it enters the propagation regime with parameter fluctuations around a certain stationary state.
Conclusions
A grid-characteristic method has been developed that makes it possible to calculate quasi-one-dimensional flows of a reacting gas; the grid lines in this method are the gas trajectories, shock waves, contact discontinuities, and sound characteristics.
Examples of calculation of flows in experimental setups have been provided and demonstrate the effectiveness of the developed method.
Author Contributions
This research was funded by the Ministry of Science and Higher Education of Russia, state assignment number FSFF-2023-0008.
Institutional Review Board Statement
The authors declare no conflict of interest.
References
Figure 1. Design models.
Figure 2. Change of calculation nodes at the intersection of grid lines.
Figure 3. Time sweep of the flow (a) and its enlarged fragment (b).
Figure 5. Time sweep of detonation initiation behind the reflected SW ( = 3). (1—Reflected SW, 2—hanging SW, 3—CD, 4—expansion wave fan).
Figure 6. The origination of a hanging SW (2) as a result of the intersection of characteristics (1).
Figure 7. Dependences of temperature on time along the main grid lines (1—wall, 2—head SW, 3—to the left of the hanging SW, 4—to the right of the hanging SW, 5—to the left of the CD, 6—to the right of the CD).
Figure 8. Dependence of the speed of the reflected SW (1) on time: (a)— = 4; (b)— = 3;
, (2)—the speed of the equilibrium reflected wave; *, (3)—the speed of the Chapman–Jouguet DW.
Figure 9. Temperature distribution (1) behind the reflected SW: (a)— = 4, t = 0.00287 s; (b)— = 3, t = 0.0152 s;
, (2)—temperature behind the equilibrium reflected wave; *, (3)—temperature behind the Chapman–Jouguet wave.
Table 1. Interaction between grid lines. RP—“the Riemann problem”; T—“transfer of parameters”; TCC—“transfer of parameters” with a change of the characteristic family; RF—shock wave reflection from the wall; RFC—reflection of the characteristic with a change of its family; R—removal of the grid line; SWO—SW origination.
1SW++RP+RP+RP+R, TCC+TT−+RF−+R, RFC
4PC++RFC+T+T, RFC−SWO+TT−+RFC−+R, RFC
9LHB+R+R, RFC+R, RFC+R+R, RFC+R−+R−−
ReactionsA; m3, mole, cmE, K
1CH4 + O2 = 2/3CO + 4/3H2O + 1/3CH46 × 108−0.226422,660
4CO2 + H2 = CO + H2O1 × 109−120,998
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.
