Журнал вычислительной математики и математической физики, 2019, T. 59, № 9, стр. 1591-1604

Численное моделирование волновых процессов при горении неоднородно распределенного заряда

И. С. Меньшов 12**, М. Ю. Немцев 13***, И. В. Семенов 13*

1 ФГУ ФНЦ НИИСИ РАН
117218 Москва, Нахимовский пр-т, 36, к.1, Россия

2 ИПМ РАН
125047 Москва, Миусская пл., 4, Россия

3 ИАП РАН
123056 Москва, ул. 2-ая Брестская, 19/18, Россия

** E-mail: menshov@kiam.ru
*** E-mail: maks-ivant@mail.ru
* E-mail: semenov@icad.org.ru

Поступила в редакцию 14.05.2019
После доработки 14.05.2019
Принята к публикации 15.05.2019

Полный текст (PDF)

Аннотация

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

Ключевые слова: двухфазная неравновесная гетерогенная смесь, внутренняя баллистика, горение гранулированных зарядов.

ВВЕДЕНИЕ

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

Экспериментальные исследования внутрибаллистических процессов, как правило, сильно ограничены и позволяют получить лишь фрагментарные данные, такие как скорость снаряда на выходе из ствола, временные зависимости давления в некоторых точках наблюдения, а также ограниченный набор данных о температуре ствола. Математическое моделирование, напротив, дает значительно больше возможностей для того, чтобы более глубоко вникнуть в суть внутрибаллистического процесса и понять основные его закономерности. Уже использование простых квазиодномерных моделей позволяет качественно описать внутрикамерный волновой процесс и получить оценки основных его параметров [1]–[4]. Более детальный анализ должен основываться на более сложных моделях, учитывающих реальное пространственное распределение заряда, горение воспламенителей, неравновесные тепловые и силовые процессы между несгоревшей твердой компонентой пороховых элементов и газообразными продуктами горения.

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

Полученная система уравнений, описывающая неравновесные динамические процессы в двухфазной дисперсной смеси, относится к так называемым в литературе неконсервативным уравнениям Эйлера [5], [6]. Неконсервативные члены в уравнениях возникают из-за соплового эффекта, связанного с неравномерным распределением по пространству объемной доли твердой компоненты.

Неконсервативная часть в системе определяющих уравнений привносит определенные трудности в математическое исследование решений и в построение адекватных дискретных моделей. Так, например, разрывные решения требуют разработки специальных математических подходов для их исследований, так как стандартная теория на основе соотношений Ренкина–Гюгонио не работает [5], [7].

Дискретизация неконсервативных членов выполняется в настоящей работе с помощью общего подхода, предложенного в работах [7], [8]. Он основан на модифицикации схемы Русанова, с помощью которой обеспечивается точное выполнение свойства однородности (well-ballancing) для неконсервативных уравнений Эйлера общего вида.

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

1. МАТЕМАТИЧЕСКАЯ МОДЕЛЬ ГОРЕНИЯ ГРАНУЛИРОВАННОГО ЗАРЯДА

Рассматривается гетерогенная смесь дисперсной фазы пороховых элементов и газовой фазы продуктов горения, которая моделируется двумя взаимопроникающими континуумами [9]–[12]. Газовая фаза состоит из гетерогенной смеси $N$ компонентов, которые соответствуют продуктам горения $N$различных пороховых элементов, и одной инертной компоненте – воздуху, заполняющим в начальный момент времени межпороховое пространство. Пороховые газы нумеруются индексом от 1 до $N$. Индекс $N + 1$ соответствует инертному воздуху.

Дисперсная фаза состоит из $N$ сортов пороха, отличающихся друг от друга законом горения, характером механического и теплового взаимодействия с газовой фазой и другими параметрами. Предполагается, что газовая фаза представляет собой многокомпонентную гетерогенную равновесную смесь с одним давлением, температурой и вектором скорости. Движение дисперсной фазы описывается одним межгранулярным давлением и одной скоростью.

Обозначим объемную долю $j$-го компонента дисперсной фазы ${{\beta }_{j}}$, суммарную объемную долю дисперсной фазы $\beta $. Тогда объемная доля газовой фазы или пористость среды $\varphi = 1 - \beta $. Материал пороха предполагается жестким и недеформируемым и характеризуется постоянной плотностью ${{\delta }_{j}}$.

Массовый состав газовой фазы характеризуется средними по пористости плотностями компонент ${{\rho }_{j}} = \frac{{{{\varphi }_{j}}{{\rho }_{{0j}}}}}{\varphi },$ где ${{\rho }_{{0j}}}$ – истинная плотность, а ${{\varphi }_{j}}$ – объемная доля так, что $\sum\nolimits_{i = 1}^{N + 1} {{{\varphi }_{i}}} = \varphi .$

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

$p\left( {1 - {{b}_{j}}{{\rho }_{{0j}}}} \right) = \frac{{{{\rho }_{{0j}}}RT}}{{{{M}_{j}}}},\quad {{\varepsilon }_{j}} = \frac{{p\left( {1 - {{b}_{j}}{{\rho }_{{0j}}}} \right)}}{{\left( {{{\gamma }_{j}} - 1} \right){{\rho }_{{0j}}}}},$
где $\gamma $ – показатель адиабаты, $b$ – коволюм, $M$ – молекулярный вес, $R$ – универсальная газовая постоянная, $\varepsilon $ – удельная внутренняя энергия. Можно показать, что в этом случае рассматриваемую смесь можно считать некой эффективной газовой средой, описываемой также уравнением состояния Дюпре, но с переменными параметрами, зависящими от массового состава смеси:
$\frac{1}{M} = \frac{1}{\rho }\sum\limits_{j = 1}^{N + 1} {\frac{{{{\rho }_{j}}}}{{{{M}_{j}}}}} ,\quad b = \frac{1}{\rho }\sum\limits_{j = 1}^{N + 1} {{{\rho }_{j}}{{b}_{j}}} ,\quad \gamma = 1 + \frac{{\sum\limits_{j = 1}^{N + 1} {{{{{\rho }_{j}}} \mathord{\left/ {\vphantom {{{{\rho }_{j}}} {{{M}_{j}}}}} \right. \kern-0em} {{{M}_{j}}}}} }}{{\sum\limits_{j = 1}^{N + 1} {{{{{\rho }_{j}}} \mathord{\left/ {\vphantom {{{{\rho }_{j}}} {\left[ {{{M}_{j}}\left( {{{\gamma }_{j}} - 1} \right)} \right]}}} \right. \kern-0em} {\left[ {{{M}_{j}}\left( {{{\gamma }_{j}} - 1} \right)} \right]}}} }},\quad \rho = \sum\limits_{j = 1}^{N + 1} {{{\rho }_{j}}} ,{\text{ }}$
где $R$ – универсальная газовая постоянная, $\varepsilon $ – удельная внутренняя энергия газовой фазы.

Мы рассматриваем течение газопороховой смеси в осесимметричном канале переменного сечения. Координаты $\left( {z,r} \right)$ соответствуют продольному и радиальному направлению. Система определяющих уравнений выражает в дифференциальной форме фундаментальные законы сохранения массы, импульса и энергии для газовой и дисперсной фазы соответственно. Обозначив через $E = 0.5(u_{1}^{2} + u_{2}^{2}) + e$ удельную полную энергию газа, ${{u}_{1}}$ и ${{u}_{2}}$ – компоненты скорости газа, ${{v}_{1}}$ и ${{v}_{2}}$ – компоненты скорости дисперсной фазы, $J = E + {p \mathord{\left/ {\vphantom {p \rho }} \right. \kern-0em} \rho }$ – удельную полную энтальпию газовой фазы, ${{\dot {m}}_{j}}$, ${{{\mathbf{\dot {\Pi }}}}_{j}}$ и ${{\dot {Q}}_{j}},{\text{ }}j = 1,\;...,\;N$, соответственно, скорости изменения массы, импульса и энергии в единице объема газовой фазы в результате горения пороховых элементов и межфазного взаимодействия, получим следующую систему уравнений для газовой фазы:

$\begin{gathered} \frac{{\partial \varphi {{\rho }_{j}}}}{{\partial t}} + \frac{1}{r}\frac{{\partial r\varphi {{\rho }_{j}}{{u}_{2}}}}{{\partial r}} + \frac{{\partial \varphi {{\rho }_{j}}{{u}_{1}}}}{{\partial z}} = {{{\dot {m}}}_{j}},\quad j = 1,\;...,\;N, \\ \frac{{\partial \varphi {{\rho }_{{N + 1}}}}}{{\partial t}} + \frac{1}{r}\frac{{\partial r\varphi {{\rho }_{{N + 1}}}{{u}_{2}}}}{{\partial r}} + \frac{{\partial \varphi {{\rho }_{{N + 1}}}{{u}_{1}}}}{{\partial z}} = 0, \\ \end{gathered} $
(1)
$\frac{{\partial \varphi \rho {{u}_{1}}}}{{\partial t}} + \frac{1}{r}\frac{{\partial \left( {r\varphi \rho {{u}_{2}}{{u}_{1}}} \right)}}{{\partial r}} + \frac{{\partial \left( {\varphi \left( {\rho {{u}_{1}}{{u}_{1}} + P} \right)} \right)}}{{\partial z}} = P\frac{{\partial \varphi }}{{\partial z}} + \sum\limits_{j = 1}^{N + 1} {{{{\left( {{{{{\mathbf{\dot {\Pi }}}}}_{j}}} \right)}}_{z}}} ,$
$\begin{gathered} \frac{{\partial \varphi \rho {{u}_{2}}}}{{\partial t}} + \frac{1}{r}\frac{{\partial \left( {r\varphi \left( {\rho {{u}_{2}}{{u}_{2}} + P} \right)} \right)}}{{\partial r}} + \frac{{\partial \left( {\varphi \rho {{u}_{2}}{{u}_{1}}} \right)}}{{\partial z}} = \frac{P}{r}\frac{{\partial r\varphi }}{{\partial r}} + \sum\limits_{j = 1}^{N + 1} {{{{\left( {{{{{\mathbf{\dot {\Pi }}}}}_{j}}} \right)}}_{r}}} , \\ \frac{\partial }{{\partial t}}\left( {\varphi \rho E} \right) + \frac{1}{r}\frac{\partial }{{\partial r}}\left( {r\varphi \rho J{{u}_{2}}} \right) + \frac{\partial }{{\partial z}}\left( {\varphi \rho J{{u}_{1}}} \right) = - P\left( {\frac{1}{r}\frac{\partial }{{\partial r}}\left( {r\beta {{v}_{2}}} \right) + \frac{\partial }{{\partial z}}\left( {\beta {{v}_{1}}} \right)} \right) + \sum\limits_{j = 1}^{N + 1} {{{{\dot {Q}}}_{j}}.} \\ \end{gathered} $

Уравнения сохранения массы и импульса для дисперсной фазы имеют следующий вид:

(2)
$\begin{gathered} \frac{{\partial {{\beta }_{j}}\delta _{j}^{0}}}{{\partial t}} + \frac{1}{r}\frac{{\partial \left( {r{{\beta }_{j}}\delta _{j}^{0}{{v}_{2}}} \right)}}{{\partial r}} + \frac{{\partial \left( {{{\beta }_{j}}\delta _{j}^{0}{{v}_{1}}} \right)}}{{\partial z}} = - {{{\dot {m}}}_{j}},\quad j = 1,\;...,\;N, \\ \frac{{\partial \beta {{\delta }^{0}}{{v}_{1}}}}{{\partial t}} + \frac{1}{r}\frac{{\partial \left( {r\beta {{\delta }^{0}}{{v}_{2}}{{v}_{1}}} \right)}}{{\partial r}} + \frac{{\partial \left( {\beta \left( {{{\delta }^{0}}{{v}_{1}}{{v}_{1}} + \sigma } \right)} \right)}}{{\partial z}} = - \beta \frac{{\partial P}}{{\partial z}} - \sum\limits_{j = 1}^{N + 1} {{{{\left( {{{{{\mathbf{\dot {\Pi }}}}}_{j}}} \right)}}_{z}}} , \\ \frac{{\partial \beta {{\delta }^{0}}{{v}_{2}}}}{{\partial t}} + \frac{1}{r}\frac{{\partial \left( {r\beta \left( {{{\delta }^{0}}{{v}_{2}}{{v}_{2}} + \sigma } \right)} \right)}}{{\partial r}} + \frac{{\partial \left( {\beta {{\delta }^{0}}{{v}_{1}}{{v}_{2}}} \right)}}{{\partial z}} = - \beta \frac{{\partial P}}{{\partial r}} - \sum\limits_{j = 1}^{N + 1} {{{{\left( {{{{{\mathbf{\dot {\Pi }}}}}_{j}}} \right)}}_{r}}} + \frac{{\beta \sigma }}{r}, \\ \end{gathered} $
где ${{\delta }^{0}} = \sum {\delta _{j}^{0}{{\beta }_{j}}} {\text{/}}\beta $ – средняя истинная плотность дисперсной фазы.

В системе уравнений (2) параметр $\sigma $ представляет собой так называемое межгранулярное давление, которое определяет сопротивление в твердой фазе, возникающее при плотном расположении частиц. Это давление возникает при достижении объемной доли дисперсной фазы критического значения, соответствующего плотной упаковке частиц. Межгранулярное давление зависит от объемной доли $\beta $ и аппроксимируется следующей эмпирической зависимостью:

$\sigma \left( \beta \right) = H\left[ {\beta - {{\beta }_{0}}} \right]B\left[ {{{{\left( {\frac{{1 - {{\beta }_{0}}}}{{1 - \beta }}} \right)}}^{k}} - 1} \right],$
где ${{\beta }_{0}}$ – объемная доля, соответствующая плотной упаковке, а параметры $B$ и $k$ – эмпирические константы, определяемые на основе экспериментальных данных.

Для описания процесса прогрева и горения пороха вводятся лагранжевые характеристики (переносимые вместе с движением твердой фазы). Это – относительная толщина сгоревшего свода порохового элемента [9] и параметры температурного профиля внутри порохового элемента. Более подробно речь о них пойдет ниже.

Величины ${{{\mathbf{\dot {\Pi }}}}_{j}}$ и ${{\dot {Q}}_{j}}$ в (1) и (2), которые определяют скорость приращения импульса и энергии в единице объема газовой фазы за счет горения пороховых элементов и межфазного взаимодействия, имеют следующий вид:

$\begin{gathered} {{{{\mathbf{\dot {\Pi }}}}}_{j}} = {{{\dot {m}}}_{j}}{\mathbf{v}} - {{{\mathbf{\tau }}}_{j}},{\text{ }}j = 1,\;...,\;N, \\ {{{{\mathbf{\dot {\Pi }}}}}_{{N + 1}}} = - {{{\mathbf{\tau }}}_{{N + 1}}}, \\ {{{\dot {Q}}}_{j}} = {{{\dot {m}}}_{j}}\left[ {\frac{{{{f}_{j}}}}{{{{\gamma }_{j}} - 1}} + 0.5\left( {v_{1}^{2} + v_{2}^{2}} \right)} \right] - \left( {{{{\mathbf{\tau }}}_{j}},{\mathbf{v}}} \right) - {{q}_{j}},\quad j = 1,\;...,\;N, \\ {{{\dot {Q}}}_{{N + 1}}} = - \left( {{{{\mathbf{\tau }}}_{{N + 1}}},{\mathbf{v}}} \right) - {{q}_{{N + 1}}}, \\ \end{gathered} $
где ${{f}_{j}}$ – сила j-го пороха, ${{q}_{j}}$ – тепловой поток между газовой и дисперсной фазами, ${{{\mathbf{\tau }}}_{j}}$ – сила межфазного трения, которые будут определены ниже.

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

В первом подходе для расчета прогрева порохового элемента решается следующая тепловая задача:

(3)
$\frac{{\partial {{T}_{j}}}}{{\partial t}} = \frac{{{{\lambda }_{j}}}}{{{{c}_{j}}{{\delta }_{j}}}}\frac{{{{\partial }^{2}}{{T}_{j}}}}{{\partial {{y}^{2}}}},\quad {{\left. {\frac{{\partial {{T}_{j}}}}{{\partial y}}} \right|}_{{y = {{e}_{j}}}}} = {{\alpha }_{j}}\left( {T - {{T}_{j}}} \right),\quad {{\left. {{{T}_{j}}} \right|}_{{y = 0}}} = {{T}_{{0j}}},\quad {{T}_{j}}\left( {y,0} \right) = {{T}_{{0j}}},$
где ${{T}_{j}}$, ${{\lambda }_{j}}$, ${{c}_{j}}$ – температура, коэффициент теплопроводности и теплоемкость $j$-го сорта пороха, ${{T}_{{0j}}}$ – его начальная температура, ${{e}_{j}}$ – половина толщины свода горения порохового элемента, $y$ – координата, которая отсчитывается от центра порохового свода к его поверхности. Воспламенение порохового элемента происходит в результате прогрева через время ${{t}_{j}}$, когда температура на его поверхности достигает значения, соответствующего заданной температуре воспламенения. После того, как температура поверхности порохового элемента достигает температуры воспламенения, температурный профиль внутри свода горения “замораживается” и может учитываться в дальнейшем для коррекции скорости горения прогретого свода [10].

Коэффициент ${{\alpha }_{j}}$ в (3) определяется через число Нуссельта:

(4)
${{\alpha }_{j}} = \frac{\lambda }{{{{\lambda }_{j}}}}{\text{N}}{{{\text{u}}}_{j}}\frac{{{{s}_{j}}}}{6},$
где ${{s}_{j}} = {{{{\chi }_{j}}} \mathord{\left/ {\vphantom {{{{\chi }_{j}}} {{{e}_{j}}}}} \right. \kern-0em} {{{e}_{j}}}}$ – площадь межфазной поверхности единицы объема дисперсной фазы для зерненых порохов, ${{\chi }_{j}}$ – безразмерный коэффициент формы порохового элемента [9]. В литературе [10]–[15] приводятся различные полуэмпирические аппроксимации числа Нуссельта ${\text{N}}{{{\text{u}}}_{j}}$ для учета конвективного теплообмена в дисперсных средах. В настоящей работе используется следующее соотношение:
${\text{Nu}}_{j}^{{}} = \left\{ \begin{gathered} 2.0 + 0.106\left( {\varphi {{{\operatorname{Re} }}_{j}}} \right){{\Pr }^{{{1 \mathord{\left/ {\vphantom {1 3}} \right. \kern-0em} 3}}}}\quad {\text{при}}\quad {{\operatorname{Re} }_{j}} \leqslant 200, \hfill \\ 2.27 + 0.6{{\left( {\varphi {{{\operatorname{Re} }}_{j}}} \right)}^{{{2 \mathord{\left/ {\vphantom {2 3}} \right. \kern-0em} 3}}}}{{\Pr }^{{{1 \mathord{\left/ {\vphantom {1 3}} \right. \kern-0em} 3}}}}\quad {\text{при}}\quad {{\operatorname{Re} }_{j}} > 200, \hfill \\ \end{gathered} \right.$
где числа Рейнольдса и Прандтля определяются соответственно как
${{\operatorname{Re} }_{j}} = \frac{{6{{e}_{j}}\rho \left| {{\mathbf{u}} - {\mathbf{v}}} \right|}}{{\mu {{\chi }_{j}}}},\quad \Pr = \frac{{\mu {{c}_{p}}}}{\lambda }.$
Здесь $\mu \left( T \right)$ – динамическая вязкость пороховых газов, ${{c}_{p}}$ – теплоемкость пороховых газов при постоянном давлении.

Приближенный аналитический подход для нахождения температуры на поверхности порохового элемента основан на аппроксимации профиля температуры внутри температурного слоя толщиной ${{h}_{j}}\left( t \right)$, на который прогревается пороховой элемент, полиномиальной зависимостью третьей степени [11]. Интегрирование (3) по пространственной переменной $\eta $ от 0 до глубины прогретого слоя ${{h}_{j}}\left( t \right)$ и по времени от 0 до текущего времени $t$ приводит к соотношению:

(5)
${{H}_{j}} \triangleq \int\limits_0^t {\left[ { - \frac{{{{\lambda }_{j}}}}{{{{c}_{j}}{{\delta }_{j}}}}\frac{{\partial {{T}_{j}}{\kern 1pt} \left( {t,0} \right)}}{{\partial x}}} \right]dt} = \int\limits_0^{h\left( t \right)} {{{T}_{j}}{\kern 1pt} \left( {t,\eta } \right)d\eta } - {{T}_{0}}{{h}_{j}}{\kern 1pt} \left( t \right).$

Неизвестные коэффициенты полинома ${{T}_{j}}{\kern 1pt} \left( {t,\eta } \right)$, аппроксимирующего температурную зависимость, находятся из условий того, что на левой границе задана первая производная температуры по $\eta $, а на правой известна температура, и предполагается, что первая и вторая производные по $\eta $ равны 0. Таким образом, коэффициенты кубического полинома оказываются выраженными через неизвестную пока толщину прогретого слоя пороха ${{h}_{j}}{\kern 1pt} \left( t \right)$.

Для ее нахождения зависимость ${{T}_{j}}{\kern 1pt} \left( {t,\eta } \right)$ с найденными коэффициентами подставляется в (5), и температура поверхности порохового элемента в рамках приближенного аналитического подхода определяется в виде

${{T}_{j}}{\kern 1pt} \left( {t,0} \right) = {{T}_{{0j}}} - \frac{{2{{H}_{j}}{{\alpha }_{j}}}}{3} + \sqrt {\frac{4}{9}H_{j}^{2}\alpha _{j}^{2} + \frac{4}{3}{{H}_{j}}{{\alpha }_{j}}{\kern 1pt} \left( {T - {{T}_{{0j}}}} \right)} {\kern 1pt} {\kern 1pt} .$
Здесь ${{\alpha }_{j}}$ – коэффициент из (4), а лагранжева переменная ${{H}_{j}}$, которая переносится вместе с пороховым элементом подобно температурному профилю, рассчитывается по (5) с использованием формул численного интегрирования.

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

${{q}_{j}} = \frac{1}{6}{{\beta }_{j}}s_{j}^{2}\lambda {\text{N}}{{{\text{u}}}_{j}}{\kern 1pt} \left( {T - {{T}_{{sj}}}} \right)H\left[ {{{t}_{j}} - t} \right],$
где ${{T}_{{sj}}}$ – температура поверхности пороха, $T$ – окружающего газа, $H\left[ . \right]$ – функция Хевисайда, ${{t}_{j}}$ – время до воспламенения j-го сорта пороха.

Сила сопротивления ${{{\mathbf{\tau }}}_{j}}$, действующая на газ со стороны j-го компонента дисперсной фазы, определяется формулой:

${{{\mathbf{\tau }}}_{j}} = \frac{1}{8}{{\beta }_{j}}{{s}_{j}}{{C}_{j}}\rho \left| {{\mathbf{u}} - {\mathbf{v}}} \right|\left( {{\mathbf{u}} - {\mathbf{v}}} \right),\quad j = 1,\;...,\;N,$
где ${{C}_{j}}$ – коэффициент сопротивления, для которого мы использовали формулу Эргана [16] в диапазоне пористости от 0.4 до 0.75. Вне этого диапазона использовались полуэмпирические зависимости [10].

Газоприход ${{\dot {m}}_{j}}$начинается по окончании воспламенительного периода и определяется степенным законом горения:

${{\dot {m}}_{j}} = {{A}_{j}}{{p}^{{{{\nu }_{j}}}}}H\left[ {t - {{t}_{j}}} \right],\quad j = 1,\;...,\;N,$
где ${{A}_{j}}$ – множитель, зависящий от геометрической формы порохового зерна, его начальной температуры, температуры внутри прогретого свода, параметров обдувающего потока (эрозионные эффекты) и химического состава пороха [10]–[17], ${{\nu }_{j}}$ – постоянный показатель степени в законе горения.

2. ЧИСЛЕННЫЙ МЕТОД

Численный метод для решения описанной выше математической модели основывается на интегральном представлении определяющей системы уравнений и методе конечных объемов [18]. Запишем систему уравнений (1), (2) в консервативном векторном виде:

(6)
$\frac{{\partial {\mathbf{q}}}}{{\partial t}} + \frac{{\partial {\mathbf{f}}}}{{\partial r}} + \frac{{\partial {\mathbf{g}}}}{{\partial z}} = {{{\mathbf{S}}}_{p}} + {{{\mathbf{S}}}_{m}} + {{{\mathbf{S}}}_{r}},$
где ${\mathbf{q}}$ – вектор консервативных переменных, ${\mathbf{f}}$ и ${\mathbf{g}}$ – потоковые векторы, ${{{\mathbf{S}}}_{p}}$, ${{{\mathbf{S}}}_{m}}$ и ${{{\mathbf{S}}}_{r}}$ – векторы правых частей, отвечающие членам с давлением (сопловые члены), межфазному обмену массой, импульсом и энергией и источниковым членам из-за осевой симметрии соответственно. Предполагая, что $\Omega $ – некоторая двумерная замкнутая область в плоскости $\left( {r,z} \right)$, ограниченная границей $\Gamma $, уравнения (6) сводятся к следующим интегральным соотношениям:
(7)
$\frac{\partial }{{\partial t}}\left[ {\int\limits_\Omega {{\mathbf{q}}d\Omega } } \right] + \int\limits_\Gamma {\left( {{\mathbf{f}}{{n}_{r}} + {\mathbf{g}}{{n}_{z}}} \right)ds} = \int\limits_\Omega {{{{\mathbf{S}}}_{p}}d\Omega } + \int\limits_\Omega {{{{\mathbf{S}}}_{m}}d\Omega } + \int\limits_\Omega {{{{\mathbf{S}}}_{r}}d\Omega } ,$
где ${\mathbf{n}} = \left( {{{n}_{r}},{{n}_{z}}} \right)$ – локальный вектор внешней нормали к границе $\Gamma $.

Применяя это интегральное соотношение к каждой ячейке расчетной сетки ${{\Omega }_{i}}$, приходим к системе дискретных алгебраических уравнений для определения средних по ячейке значений консервативного вектора ${\mathbf{q}}_{i}^{{n + 1}}$ на временном слое ${{t}^{{n + 1}}}$ по соответствующим значениям ${\mathbf{q}}_{i}^{n}$ с предыдущего слоя ${{t}^{n}}$:

(8)
${\mathbf{q}}_{i}^{{n + 1}} = {\mathbf{q}}_{i}^{n} - \frac{1}{{{{\omega }_{i}}}}\int\limits_{{{t}^{n}}}^{{{t}^{{n + 1}}}} {\left[ {\sum\limits_i {{{{\left( {{\mathbf{f}}{{n}_{r}} + {\mathbf{g}}{{n}_{z}}} \right)}}_{i}}{{s}_{i}}} } \right]dt} + \int\limits_{{{t}^{n}}}^{{{t}^{{n + 1}}}} {\left[ {{{{\left( {{{{\mathbf{S}}}_{p}}} \right)}}_{i}} + {{{\left( {{{{\mathbf{S}}}_{m}}} \right)}}_{i}} + {{{\left( {{{{\mathbf{S}}}_{r}}} \right)}}_{i}}} \right]dt} ,$
где ${{\omega }_{i}}$ – площадь контрольного объема ${{\Omega }_{i}}$, ${{s}_{i}}$ – длина бокового ребра, суммирование в правой части ведется по всем ребрам ячейки. Основным элементом расчетной схемы является аппроксимация потоковых членов через ребра ячеек (второе слагаемое в правой части уравнения (8)).

Дальнейшая дискретизация уравнений (8) строится на принципе разделения по физическим процессам. Процесс перехода с временного слоя ${{t}^{n}}$ на временной слой ${{t}^{{n + 1}}}$ разделяется на три составляющие:

1) динамика газовой фазы, описываемая системой уравнений:

(9)
${\mathbf{q}}_{{gi}}^{{n + 1}} = {\mathbf{q}}_{{gi}}^{n} - \frac{1}{{{{\omega }_{i}}}}\int\limits_{{{t}^{n}}}^{{{t}^{{n + 1}}}} {\left[ {\sum\limits_i {{{{\left( {{{{\mathbf{f}}}_{g}}{{n}_{r}} + {{{\mathbf{g}}}_{g}}{{n}_{z}}} \right)}}_{i}}{{s}_{i}}} } \right]dt} + \int\limits_{{{t}^{n}}}^{{{t}^{{n + 1}}}} {\left[ {{{{\left( {{{{\mathbf{S}}}_{{pg}}}} \right)}}_{i}} + {{{\left( {{{{\mathbf{S}}}_{{rg}}}} \right)}}_{i}}} \right]dt} ;$

2) динамика дисперсной фазы, описываемая системой уравнений:

(10)
${\mathbf{q}}_{{si}}^{{n + 1}} = {\mathbf{q}}_{{si}}^{n} - \frac{1}{{{{\omega }_{i}}}}\int\limits_{{{t}^{n}}}^{{{t}^{{n + 1}}}} {\left[ {\sum\limits_i {{{{\left( {{{{\mathbf{f}}}_{s}}{{n}_{r}} + {{{\mathbf{g}}}_{g}}{{n}_{z}}} \right)}}_{i}}{{s}_{i}}} } \right]dt} + \int\limits_{{{t}^{n}}}^{{{t}^{{n + 1}}}} {\left[ {{{{\left( {{{{\mathbf{S}}}_{{ps}}}} \right)}}_{i}} + {{{\left( {{{{\mathbf{S}}}_{{rs}}}} \right)}}_{i}}} \right]dt} ;$

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

(11)
${\mathbf{q}}_{i}^{{n + 1}} = {\mathbf{q}}_{i}^{n} + \int\limits_{{{t}^{n}}}^{{{t}^{{n + 1}}}} {{{{\left( {{{{\mathbf{S}}}_{m}}} \right)}}_{i}}dt} .$

При решении уравнений (9) на временном интервале $[{{t}^{n}},{{t}^{{n + 1}}}]$ параметры дисперсной фазы ${{{\mathbf{q}}}_{s}}$ предполагают неизменными (“замороженными” во времени), а для уравнений (10), напротив, неизменными являются параметры газовой фазы ${{{\mathbf{q}}}_{g}}$.

Дискретизация уравнений (9) и (10) проводится единообразно с использованием интерполяционной предиктор-корректор схемы второго порядка точности по времени и пространственным переменным. Это – явная двухстадийная схема, на первом этапе которой считаются промежуточные величины на половинном шаге ${{t}^{{n + {1 \mathord{\left/ {\vphantom {1 2}} \right. \kern-0em} 2}}}}$:

(12)
${\mathbf{q}}_{i}^{ * } = {\mathbf{q}}_{i}^{n} - 0.5\frac{{\Delta t}}{{{{\omega }_{i}}}}\sum\limits_\sigma {{{{\left( {{\mathbf{f}}{\kern 1pt} \left( {{{{\mathbf{q}}}_{\sigma }}} \right){{n}_{r}} + {\mathbf{g}}{\kern 1pt} \left( {{{{\mathbf{q}}}_{\sigma }}} \right){{n}_{z}}} \right)}}_{\sigma }}{{s}_{\sigma }}} + 0.5\left[ {{{{\left( {{{{\mathbf{S}}}_{p}}} \right)}}_{i}} + {{{\left( {{{{\mathbf{S}}}_{r}}} \right)}}_{i}}} \right]\Delta t,$
где ${{{\mathbf{q}}}_{\sigma }}$ – значения компонентов вектора ${\mathbf{q}}$, полученные интерполяцией значений из центра ячейки с использованием градиента $\nabla {\mathbf{q}}_{i}^{n}$. Вектор градиента здесь рассчитывается по локальному шаблону методом наименьших квадратов [19]. Интерполяция основывается на MUSCL-подходе с использованием противопоточной схемы третьего порядка точности (на равномерной сетке) [20].

На втором этапе производится расчет значений компонент определяющего вектора консервативных переменных на новом временном слое ${{t}^{{n + 1}}}$:

${\mathbf{q}}_{i}^{{n + 1}} = {\mathbf{q}}_{i}^{n} - \frac{{\Delta t}}{{{{\omega }_{i}}}}\sum\limits_\sigma {T_{\sigma }^{{ - 1}}{\mathbf{F}}_{\sigma }^{ * }{{s}_{\sigma }}} + [{{({\mathbf{S}}_{p}^{*})}_{i}} + {{({\mathbf{S}}_{r}^{ * })}_{i}}]\Delta t,$
где ${{{\mathbf{F}}}_{\sigma }}$ – вектор численного потока в направлении внешней нормали к ребру счетной ячейки. Верхний индекс “*” означает, что значения вычисляются по значениям соответствующих величин с промежуточного временного слоя (уравнение (12)).

Аппроксимация численного потока ${{{\mathbf{F}}}_{\sigma }}$ в случае непрерывной пористости выполняется по методу С.К. Годунова [21]. На ребрах с разрывом пористости используется модифицированный поток Русанова. Соответственно в двухволновом приближении решаются задачи Римана о распаде произвольного разрыва для газовой и дисперсной фазы с начальными данными, полученными интерполяцией величин промежуточного слоя из центров ячеек на их ребра:

$\begin{gathered} {\mathbf{Q}}_{{g\sigma }}^{{ * R}} = {\mathbf{Q}}_{{g\sigma }}^{{ * R}}(\lambda ,{\mathbf{Q}}_{{g\sigma }}^{{ * - }},{\mathbf{Q}}_{{g\sigma }}^{{ * + }}), \hfill \\ {\mathbf{Q}}_{{s\sigma }}^{{ * R}} = {\mathbf{Q}}_{{s\sigma }}^{{ * R}}(\lambda ,{\mathbf{Q}}_{{s\sigma }}^{{ * - }},{\mathbf{Q}}_{{s\sigma }}^{{ * + }}), \hfill \\ \end{gathered} $
где $\lambda $ – автомодельная переменная в локальном базисе ребра ячейки, верхние индексы “–” и “+” обозначают экстраполированные на ребро значения из текущей и соседней ячейки соответственно. Детали аппроксимации численного потока с использованием модифицированной схемы Русанова можно найти в работах [8], [22].

Решение уравнений (11) третьего этапа выполняется с помощью квадратурных формул, что завершает вычислительный алгоритм перехода с временного слоя ${{t}^{n}}$ на слой ${{t}^{{n + 1}}}$.

При моделировании внутрибаллистического процесса в двумерной осесимметричной постановке необходимо рассчитывать взаимодействие пороховых газов с подвижным снарядом. Для этого предлагается использовать метод свободной границы [23]. Кратко опишем модификацию этого метода для уравнений двухфазной гидродинамики.

Рассмотрим некоторую область $\Omega \subset {{R}^{2}}$, ограниченную замкнутой гладкой кривой S. Пусть ${\mathbf{n}}$ – внешняя нормаль к S. Пусть также $\Gamma \in {{R}^{2}}$ представляет собой поверхность твердого тела, находящегося в контакте со средой. Учесть влияние твердого объекта на течение среды можно с помощью введения в правую часть уравнений эффективного вектора ${{{\mathbf{F}}}_{w}}$. Система уравнений (7) тогда модифицируется следующим образом:

(13)
$\frac{\partial }{{\partial t}}\left[ {\int\limits_\Omega {{\mathbf{q}}d\Omega } } \right] + \int\limits_S {\left( {{\mathbf{f}}{{n}_{r}} + {\mathbf{g}}{{n}_{z}}} \right)ds} = \int\limits_\Omega {{{{\mathbf{S}}}_{p}}d\Omega } + \int\limits_\Omega {{{{\mathbf{S}}}_{m}}d\Omega } + \int\limits_\Omega {{{{\mathbf{S}}}_{r}}d\Omega } - \int\limits_{\gamma (\Omega )} {{{{\mathbf{F}}}_{w}}dS} ,$
где $\gamma \left( \Omega \right) = \Omega \cap \Gamma $. Эффективный поток ${{{\mathbf{F}}}_{w}}$ должен препятствовать потерям массы, импульса и энергии на границе $\gamma \left( \Omega \right)$ из-за протекания газа и генерировать импульс давления, эквивалентный импульсу, порождаемому в результате взаимодействия среды с поверхностью тела:
(14)
$\begin{gathered} {{{\mathbf{F}}}_{{gw}}} = \left[ {\begin{array}{*{20}{c}} {\alpha {{\rho }_{i}}\left( {{\mathbf{u}} - {{{\mathbf{U}}}_{s}},{\mathbf{n}}} \right)} \\ {\alpha \rho \left( {{\mathbf{u}} - {{{\mathbf{U}}}_{s}},{\mathbf{n}}} \right){\mathbf{u}} + \alpha \left( {p - {{p}_{w}}} \right){\mathbf{n}}} \\ {\alpha \rho \left( {{\mathbf{u}} - {{{\mathbf{U}}}_{s}},{\mathbf{n}}} \right)E + \alpha \left( {p{\mathbf{u}} - {{p}_{w}}{{{\mathbf{U}}}_{s}},{\mathbf{n}}} \right)} \end{array}} \right], \\ {{{\mathbf{F}}}_{{sw}}} = \left[ {\begin{array}{*{20}{c}} {{{\delta }_{i}}{{\beta }_{i}}{{z}_{i}}\left( {{\mathbf{v}} - {{{\mathbf{U}}}_{s}},{\mathbf{n}}} \right)} \\ {{{\delta }_{i}}{{\beta }_{i}}\left( {{\mathbf{v}} - {{{\mathbf{U}}}_{s}},{\mathbf{n}}} \right)} \\ {{{\delta }^{0}}\beta \left( {{\mathbf{v}} - {{{\mathbf{U}}}_{s}},{\mathbf{n}}} \right){\mathbf{v}} + \beta \left( {\sigma - {{\sigma }_{w}}} \right){\mathbf{n}}} \end{array}} \right]. \\ \end{gathered} $
Здесь ${{p}_{w}}$ – давление, возникающее на поверхности $\Gamma $ в результате взаимодействия газовой фазы с поверхностью тела, ${{\sigma }_{w}}$ – межгранулярное псевдодавление, возникающее на поверхности $\Gamma $ в результате взаимодействия дисперсной фазы с поверхностью тела, ${{{\mathbf{U}}}_{s}}$ – вектор локальной скорости поверхности тела.

Расчетная схема при наличии в потоке среды твердых объектов строится следующим образом. Вначале проводится сквозной счет по базовой схеме, описанной выше. При этом не принимается во внимание присутствие поверхности твердого тела. Полученные значения вектора решения ${{{\mathbf{\tilde {q}}}}^{{n + 1}}}$ модифицируются затем с помощью уравнения:

(15)
${\mathbf{q}}_{i}^{{n + 1}} = {\mathbf{\tilde {q}}}_{i}^{{n + 1}} - \frac{{\Delta t}}{{{{\vartheta }_{i}}{{\omega }_{i}}}}\left( {{{{\mathbf{F}}}_{w}}} \right)_{\gamma }^{{n + 1}}{{s}_{\gamma }},$
где ${{s}_{\gamma }}$ – длина кривой $\gamma (\Omega )$, ${{\vartheta }_{i}}$ – часть объема ${{\omega }_{i}}$, занимаемого средой в пересекаемой ячейке. Уравнение (15) решается итерационным методом Ньютона.

3. ЧИСЛЕННЫЕ РЕЗУЛЬТАТЫ: ВНУТРЕННЯЯ БАЛЛИСТИКА СОСТАВНОГО ЗАРЯДА

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

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

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

На фиг. 1 представлены зависимости давления от времени для 6 различных точек наблюдения, расположенных по длине трубы. Хорошо видно, что для подвижного заряда амплитуда волнового процесса существенно ниже, чем для неподвижного. Локализация процесса горения в узкой области у затвора для неподвижного заряда приводит к более быстрому выгоранию заряда и большему росту давления по сравнению с подвижным зарядом. В результате торможения пороховых газов о неподвижный снаряд возникает отошедшая волна, которая для подвижного заряда имеет меньшую амплитуду и более пологий фронт, чем для неподвижного.

Фиг. 1.

Динамика давления в 6 различных положениях внутри трубы для заряда из 2-х модулей: (а) заряд подвижный, (б) заряд закреплен у затвора.

В случае подвижного заряда отошедшая волна раньше приходит к затвору и отражается от него. Для неподвижного заряда отошедшая от снаряда волна приходит к затвору в момент, когда весь заряд уже выгорел, а отток газов в сторону снаряда привел к падению давления у затвора до 200 атмосфер. Приход волны сопровождается быстрым ростом давления до 1400 атмосфер.

На фиг. 2 и 3 представлены распределения давления и объемной доли заряда вдоль продольной координаты на расстоянии 80 мм от оси симметрии трубы для подвижного и неподвижного зарядов. Видно, что перемещение заряда в трубе приводит к существенному снижению амплитуды волнового процесса за счет более равномерного газоприхода и снижения скорости пороховых газов во время горения заряда. Это приводит также к образованию более слабой отраженной волны, чем в случае с неподвижным зарядом.

Фиг. 2.

Распределения вдоль линии R = 80 мм, для нескольких моментов времени, подвижный заряд из 2-х модулей: (а) давление, (б) объемная доля основного заряда.

Фиг. 3.

Распределения вдоль линии R = 80 мм, для нескольких моментов времени, неподвижный заряд из 2-х модулей: (а) давление, (б) объeмная доля заряда.

Для исследования влияния начального расположения заряда и источника воспламенения на интенсивность волнового процесса провели численные эксперименты для пятимодульного заряда [24]. Размер и состав модулей аналогичны случаю двухмодульного заряда.

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

Фиг. 4.

Динамика давления в 5 различных положениях внутри трубы для заряда из 5 модулей: (а) заряд у затвора, (б) заряд у снаряда. Воспламенение начинается с первого модуля со стороны затвора.

Динамика движения заряда и его локализация в обоих случаях различны, как это следует из фиг. 5 и 6, где приведены распределения давления и объемной доли на разные моменты времени для разных начальных расположений заряда соответственно. Когда заряд расположен у снаряда, в каморе присутствуют две свободные от заряда области: между затвором и первым модулем, а также между 5 модулем и снарядом. Объем этих областей различен; со стороны затвора он больше, что приводит к более слабому темпу роста давления вблизи затвора. Присутствие свободной области вблизи затвора приводит также к затягиванию процесса воспламенения остальных модулей, большей неоднородности и снижению скорости горения заряда в целом.

Фиг. 5.

Распределения вдоль линии R = 80 мм, для нескольких моментов времени, заряд из 5 модулей у снаряда: (а) давление, (б) объемная доля заряда.

Фиг. 6.

Распределения вдоль линии R = 80 мм, для нескольких моментов времени, заряд из 5 модулей у затвора: (а) давление, (б) объемная доля заряда.

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

ЗАКЛЮЧЕНИЕ

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

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

2. Показано, что характер волнового процесса сильно зависит от расположения источника воспламенения, а также расположения самого заряда в трубе. Так, для случая пятимодульного заряда с воспламенением первого модуля расположение его со стороны затвора приводит к более интенсивному волновому процессу по сравнению с расположением аналогичного заряда у снаряда. Описаны механизм формирования волнового процесса и причины его усиления.

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

Список литературы

  1. Ищенко А.Н., Касимов В.З., Ушакова О.В. // Матер. конференции “Современная баллистика и смежные вопросы механики”. Томск, 2009. С. 85.

  2. Wang Y.-W., Guo Y.-H., Wang S.-C. // Proc. 23d Intern. Symp. on Ballistics. Tarragona, Spain, 2007. P. 401.

  3. Семенов И.В., Меньшов И.С., Уткин П.С., Ахмедьянов И.Ф. БАРС-1МП – программный комплекс для численного исследования внутрибаллистических процессов на многопроцессорных ЭВМ // Изв. высших учебных заведений. Физика. 2013. Т. 56. № 6–3. С. 61–63.

  4. Закаменных Г.И., Чернов В.В., Абдуллин А.К., Семенов И.В., Уткин П.С., Лебедева А.Ю., Ахмедьянов И.Ф. Экспериментальное и численное исследование влияния положения модульного заряда на характеристики выстрела // Сборник материалов Всероссийской научно-технической конференции “Фундаментальные основы баллистического проектирования”. Санкт-Петербург, 28 июня–2 июля 2010 г. Т. 1 / Под ред. д.т.н. проф. Кэрта Б.Э. СПб.: Балт. гос. техн. ун-т, 2010. С. 119–122.

  5. Dal Maso G., Le Floch P.G., Murat F. Definition and weak stability of a non-conservative product // J. Math. Pures. Appl. 1995. V. 74 (6). P. 483–548.

  6. Clain S., Rochette D. First- and second-order finite volume methods for the one-dimensional nonconservative Euler system // J. Comp. Phys. 2009. V. 228. P. 8214–8248.

  7. Gosse L. A well-balanced flux-vector splitting scheme designed for hyperbolic system of conservation laws with source terms // Comput. Math. Appl. 2000. V. 39. P. 135–159.

  8. Menshov I.S. Exact and approximate Riemann solvers for compressible two-phase flows // Math. Models and Comput. Simulat. 2017. V. 9. № 4. P. 405–422.

  9. Русяк И.Г., Ушаков В.М. Внутрикамерные гетерогенные процессы в ствольных системах. Екатеринбург: УрО РАН, 2001. 259 с.

  10. Хоменко Ю.П., Ищенко А.Н., Касимов В.З. Математическое моделирование внутрибаллистических процессов в ствольных системах. Новосибирск: Изд-во СО РАН, 1999. 255 с.

  11. Nusca M.J., Conroy P.J. Multiphase CFD Simulations of Solid Propellant Combustion in Gun Systems // Proc. of DoD High Performance Computing Modernization Program 2001 Users Group Conference, Biloxi, MS, 18–21 June 2001.

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

  13. Захаренков В.Ф. Полуэмпирический метод расчета теплообмена в гладких и шероховатых трубах при течении горячих газовых потоков // Труды Междунар. научно-практической конференции “Третьи Окуневские чтения”. Санкт-Петербург, 24–29 июня 2002 г. Т. 2. С. 176–185.

  14. Кутателадзе С.С., Боришанский В.М. Справочник по теплопередаче. М.: Госэнергоиздат, 1958.

  15. Чудновский А.Ф. Теплообмен в дисперсных средах. М.: Гостех., 1954.

  16. Ergun S. Fluid flow through packed columns // Chemical Engng Progress. 1952. V. 48. № 2. P. 89–94.

  17. Вилюнов В.Н. К теории эрозионного горения порохов // Докл. АН СССР. 1961. Т. 136. № 2. С. 381–383.

  18. Barth T., Ohlberger M. Finite volume methods: foundation and analysis // Encyclopedia of Comput. Mechanics. 2004. V. 1. P. 439–470.

  19. Gossler A. Moving Least Squares: A Numerical Differentiation Method for Irregularly Spaced Calculation Points // Sandia Report SAND2001-1669. 2001.

  20. Van Leer B. Towards the ultimate conservative difference scheme. V – A second-order sequel to Godunov’s method // J. of Comput. Phys. 1979. V. 32. P. 101–136.

  21. Годунов С.К., Забродин А.В., Иванов М.Я. и др. Численное решение многомерных задач газовой динамики. М.: Наука, 1976.

  22. Menshov I., Serezhkin A. A generalized Rusanov method for the Baer–Nunziato equations with application to DDT processes in condensed porous explosives // Int. J. Numer. Meth. Fluids. 2017. P. 1–19.

  23. Меньшов И.С., Павлухин П.В. Эффективный параллельный метод сквозного счета задач аэродинамики на несвязных декартовых сетках // Ж. вычисл. матем. и матем. физ. 2016. Т. 56. № 9. С. 1677–1691.

  24. Семенов И.В., Меньшов И.С., Немцев М.Ю. Математическое моделирование осесимметричных внутрибаллистических процессов // Препр. ИПМ им. М.В. Келдыша, 2017. № 143. 20 с.

Дополнительные материалы отсутствуют.