Численное моделирование взаимодействия деформируемого газопроницаемого фрагмента гранулированного слоя с ударной волной в трехмерной постановке

Полный текст

Введение

Гранулированные слои являются перспективным элементом для защиты ответственных конструкций от действия ударных волн. Защитные свойства слоев из различных материалов достаточно хорошо изучены для действия акустических и слабых ударных волн [1–5]. Исследования проводились экспериментальными, аналитическими и численными методами. Взаимодействия гранулированных слоев с интенсивными ударными волнами в ударных трубах рассматривались в основном экспериментальными методами в циклах работ [6–12] и численно в [13]. В этих работах исследовались песчаные гранулированные слои, в результате воздействия ударных волн на которые происходило разрушение и разлет частиц. В силу ограниченности возможностей экспериментальных подходов в части измерений и наблюдений многие особенности протекающих процессов остались невыясненными. Современные численные методы имеют значительный исследовательский потенциал и позволяют решать сложные нелинейные задачи, включающие распространение ударных волн в твердых телах, в газах и в гранулированных средах. В данной работе для моделирования трехмерных быстропротекающих процессов и в газе, и в деформируемом теле используется единый численный метод — модифицированный метод Годунова [14–23], основанный на интегрировании законов сохранения с использованием решения задачи «распада разрыва» в газе, деформируемом теле и на границе «газ -- деформируемое тело», имеющий второй порядок аппроксимации на гладких решениях и монотонный на разрывных. Повышение точности достигается только за счет модификации решения задачи распада разрыва на шаге «предиктор» схемы путем сближения областей влияния разностной и дифференциальной задач [14]. 

1. Постановка задачи

Постановка задачи соответствует условиям проведения экспериментальных исследований [12]. В квадратной стальной трубе (рис. 1) в левой ударной секции находится воздух под высоким давлением. В правой секции находится воздух низкого давления и песчаный гранулированный слой, расположенный поперек трубы. Его свободные поверхности с двух сторон закрыты бумажными экранами. Формирование ударной волны производится на границе двух камер путем разрыва диафрагмы. В расчетах использовались параметры ударной волны, близкие по давлению, плотности и скорости к экспериментальным. Требуется описать в связанной постановке процессы формирования, распространения ударной волны по длине трубы и ее взаимодействия с гранулированным деформируемым слоем.

Рис. 1. Экспериментальная установка
[Figure 1. Schematic front view of the shock tube]

Для численного моделирования используются уравнения динамики сплошных сред в векторном виде [14]: 
\[ \begin{equation}
\frac\partial{\partial t } {\bf u}+
\frac\partial{\partial x } {\bf f}+
\frac\partial{\partial y } {\bf g}+
\frac\partial{\partial z } {\bf h} ={ \bf k},
\end{equation} \tag{1} \]
где
\[ \begin{equation*}
{\bf u}=\begin{vmatrix} 
\rho \\
\rho{u} \\
\rho{v} \\
\rho{w} \\
e \\
s_{xx} \\
s_{yy} \\
s_{zz} \\
s_{xy} \\
s_{xz} \\
s_{yz} 
\end{vmatrix}, 
\qquad\qquad
{\bf f}=\begin{vmatrix} 
\rho{u} \\
\rho{u^2}+p-s_{xx} \\
\rho{uv}-s_{xy} \\
\rho{uw}-s_{xz} \\
(e+p-s_{xx}){u}-s_{xy}v-s_{xz}w \\
u(s_{xx}-\frac4{3}\mu) \\
u(s_{yy}+\frac2{3}\mu) \\
u(s_{zz}+\frac2{3}\mu) \\
{u}s_{xy}-\mu{v} \\
{u}s_{xz}-\mu{w} \\
us_{yz} 
\end{vmatrix},
\end{equation*} \]
\[ \begin{equation*}
{\bf g}=\begin{vmatrix} 
\rho{v} \\
\rho{uv}+p-s_{xy} \\
\rho{v^2}+p-s_{yy} \\
\rho{vw}-s_{yz} \\
(e+p-s_{yy}){v}-s_{xy}u-s_{yz}w \\
v(s_{xx}+\frac2{3}\mu) \\
v(s_{yy}-\frac4{3}\mu) \\
v(s_{zz}+\frac2{3}\mu) \\
{v}s_{xy}-\mu{u} \\
{v}s_{xz} \\
{u}s_{yz}-\mu{w} 
\end{vmatrix},
\qquad
{\bf h}=\begin{vmatrix} 
\rho{w} \\
\rho{uw}-s_{xz} \\ 
\rho{vv}-s_{yz} \\ 
\rho{w^2}+p-s_{zz} \\ 
(e+p-s_{zz}){w}-s_{yz}{v}-s_{xz}u \\ 
w(s_{xx}+\frac2{3}\mu) \\ 
w(s_{yy}+\frac2{3}\mu) \\ 
w(s_{zz}-\frac4{3}\mu) \\ 
{w}s_{xy} \\ 
{u}s_{xz}-\mu{u} \\ 
us_{yz}-\mu{v} 
\end{vmatrix},
\end{equation*} \]
\[ \begin{equation*}
{\bf k}=
\begin{vmatrix} 
0 \\ 
0 \\ 
0 \\ 
0 \\ 
0 \\
s_{xx}(\frac{\partial{u}}{\partial{x}}+\frac{\partial{v}}{\partial{x}}+\frac{\partial{w}}{\partial{z}})+s_{xy}(\frac{\partial{u}}{\partial{y}}-\frac{\partial{v}}{\partial{x}})+s_{xz}(\frac{\partial{u}}{\partial{z}}-\frac{\partial{w}}{\partial{x}})-\lambda{s_{xx}} \\ 
s_{xy}(\frac{\partial{v}}{\partial{x}}-\frac{\partial{u}}{\partial{y}})+s_{yy}(\frac{\partial{u}}{\partial{x}}+\frac{\partial{v}}{\partial{x}}+\frac{\partial{w}}{\partial{z}})+s_{xz}(\frac{\partial{v}}{\partial{z}}-\frac{\partial{w}}{\partial{y}})-\lambda{s_{yy}} \\ 
s_{xz}(\frac{\partial{w}}{\partial{x}}-\frac{\partial{u}}{\partial{z}})+s_{yz}(\frac{\partial{w}}{\partial{y}}-\frac{\partial{v}}{\partial{z}})+s_{xz}(\frac{\partial{u}}{\partial{x}}+\frac{\partial{v}}{\partial{x}}+\frac{\partial{w}}{\partial{z}})-\lambda{s_{zz}} \\ 
s_{xy}(\frac{\partial{u}}{\partial{x}}+\frac{\partial{v}}{\partial{x}}+\frac{\partial{w}}{\partial{z}})-\frac1{2}(s_{xx}-s_{yy})(\frac{\partial{u}}{\partial{y}}-\frac{\partial{w}}{\partial{x}})+\frac1{2}s_{xz}(\frac{\partial{v}}{\partial{z}}-\frac{\partial{w}}{\partial{y}})+\frac1{2}s_{yz}(\frac{\partial{v}}{\partial{z}}-\frac{\partial{w}}{\partial{y}})-\lambda{s_{xy}} \\ 
s_{xy}(\frac{\partial{u}}{\partial{x}}+\frac{\partial{v}}{\partial{x}}+\frac{\partial{w}}{\partial{z}})-\frac1{2}(s_{xx}-s_{zz})(\frac{\partial{u}}{\partial{y}}-\frac{\partial{w}}{\partial{x}})+\frac1{2}s_{xy}(\frac{\partial{w}}{\partial{y}}-\frac{\partial{v}}{\partial{z}})+\frac1{2}s_{yz}(\frac{\partial{u}}{\partial{z}}-\frac{\partial{w}}{\partial{x}})-\lambda{s_{xz}} \\ 
s_{xy}(\frac{\partial{u}}{\partial{x}}+\frac{\partial{v}}{\partial{x}}+\frac{\partial{w}}{\partial{z}})-\frac1{2}(s_{yy}-s_{zz})(\frac{\partial{u}}{\partial{y}}-\frac{\partial{w}}{\partial{x}})+\frac1{2}s_{xz}(\frac{\partial{v}}{\partial{x}}-\frac{\partial{u}}{\partial{y}})+\frac1{2}s_{xy}(\frac{\partial{w}}{\partial{x}}-\frac{\partial{u}}{\partial{z}})-\lambda{s_{yz}} 
\end{vmatrix}.
\end{equation*} \]
Здесь $e=\rho\bigl(\varepsilon+\frac12(u^2+v^2+w^2)\bigr)$. Система (1) замыкается уравнением состояния (УРС) для шаровых компонент в форме 
\[ \begin{equation}
\varepsilon=\varepsilon(p,\rho).
\end{equation} \tag{2} \]

В соотношениях (1), (2) используются следующие обозначения: $p$ — давление; $\rho$ — плотность; $u$, $v$, $w$ — компоненты скоростей по $x$, $y$ и $z$; $\varepsilon$ — внутренняя энергия единицы массы; $e$ — полная механическая энергия единицы объема сплошной среды; $s_{xx}$, $s_{yy}$, $s_{zz}$, $s_{xy}$, $s_{xz}$, $s_{yz}$ – компоненты девиатора тензора истинных напряжений Эйлера; $\mu$ — модуль сдвига. 

Первые пять уравнений системы (1) представляют собой законы сохранения массы, импульса и энергии. Следующие шесть уравнений — физические соотношения упругости и пластичности с учетом поворота тензора напряжений в эйлеровых переменных (производная Яумана), записанные в дифференциальной форме. К системе (1), (2) добавляются начальные и краевые условия. В случае отсутствия сдвиговых напряжений система (1) переходит в уравнения Эйлера для движения сжимаемого газа [14].

Для описания процессов в плотных сжимаемых средах применяется баротропное УРС вида $p=p(\rho)$, позволяющее избежать интегрирования уравнения сохранения энергии. Для песчаных частиц плотностью $\rho_0$ это будет УРС идеального упругопластического тела $p=K\varepsilon_v$, где $K$ — модуль объемного сжатия, $\varepsilon_v=1- {\rho_0}/{\rho}$ — объемная деформация.

Критерием перехода из упругого напряженно-деформированного состояния в пластическое является условие текучести Мизеса: 
\[ \begin{equation*}
J_2=\frac12 s_{ij}s_{ij}\leqslant \sigma_T^2,
\end{equation*} \]
где $J_2$ — второй инвариант девиатора тензора напряжений $s_{ij}$, $\sigma_T$ — предел текучести. В соответствии с [20] в этом случае происходит коррекция компонент девиатора напряжений умножением на $\lambda=\sigma_T/\sqrt{3J_2}$. Для воздуха принимается $s_{xx}=s_{yy}=s_{zz}=s_{xy}=s_{yz}=s_{xz}=0$, $\varepsilon=p/[\rho(\gamma-1)]$, где $\gamma$ — показатель адиабаты.

Интегральная форма (1), на базе которой строится разностная схема, имеет вид
\[ \begin{equation}
\int\!\!\!\int\!\!\!\int _{\omega} \bigl({\bf u}\, dxdydz+{\bf f}\, dydzdt+{\bf g}\, dxdzdt+
{\bf h}\, dxdydt \bigr)= \int\!\!\!\int\!\!\!\int\!\!\!\int _\Omega {\bf k}\, dxdydzdt,
\end{equation} \tag{3} \]
где $\Omega$ — любой замкнутый объем, поверхность $\omega$ которого гомеоморфна сфере в четырехмерном пространстве $(x, y, z, t)$.

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

2. Многосеточный метод численного решения

Решение уравнений (1)–(3) производится методом Годунова повышенной точности [14, 17], единым как для газодинамических, так и для упругопластических течений.

Данная модификация позволяет повысить порядок аппроксимации схемы до второго на гладких решениях, сохранив монотонность на разрывных, без изменения разностного шаблона явной двухшаговой схемы, изменив только шаг «предиктор» [14].

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

Этап «корректор» численного интегрирования уравнений (3) остается неизменным и совпадает с классической схемой.

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

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

3. Исследование взаимодействия ударной волны с гранулированным слоем

Математическая постановка задачи основывается на следующих предположениях.

1. Гранулированный слой состоит из совокупности одинаковых шаровых деформируемых кварцевых частиц, представляющей собой кубическую упаковку (рис. 2).
2. Пространство между частицами заполнено сжимаемой газовой средой (воздухом).
3. Кубическая упаковка предполагает симметрию расположения частиц вдоль осей координат.
4. Воздействующая ударная волна распространятся вдоль оси $X$. 

Выделяется симметричный элемент упаковки в виде последовательности шаровых частиц диаметром 1 мм.

Для демонстрации работоспособности численной методики предполагается, что многослойная гранулированная среда в направлении оси $X$ состоит из трех слоев частиц, ограниченных в канале квадратного сечения с жесткими стенками. Фрагмент расчетной области представлен на рис. 3. Фактически рассматривается взаимодействие шаровых тел и газа в канале квадратного сечения с жесткими стенками. На рис. 3, b сетка внутри окружности — сетка по шаровой песчинке.

Рис. 2. Вид гранулированного слоя
[Figure 2. View of a granular layer]

Рис. 3. Расчетная область (a) и сетка в поперечном сечении (b)
[Figure 3. Computational domain (a) and mesh in cross section (b)]

Коэффициент газопроницаемости поперечного сечения, то есть отношение минимальной площади, занимаемой газом, к площади бокового сечения составляет 0.215, объемное содержание воздуха — 0.476. На боковых стенках канала общей длиной 21 мм по оси $X$ реализуются условия жесткой стенки, а на граничных плоскостях $X=-1.5$ мм и $X=19.5$ мм выполняются условия «свободного вытока» [16]. Для кварца, из которого состоят тела, принимаем, что модуль объемного сжатия $K=109$ ГПа, модуль сдвига $\mu=43$ ГПа, предел текучести $\sigma_T=0.17$ ГПа, плотность $\rho=2650$ кг/м$^3$. Для воздуха принимаем, что показатель адиабаты $\gamma=1.4$, плотность $\rho_0=1.23$ кг/м$^3$.

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

В начальный момент времени все среды предполагаются невозмущенными, кроме области газа. Перед первым шариком длиной 1.5 мм по оси $X$ в начальный момент задаются параметры, как за фронтом отраженной ударной волны: давление 1.6 МПа, плотность 13.2 кг/м$^3$ и нулевая скорость.

Размер расчетной области по оси $X$ составляет 21 мм, по оси $Y$ — 1 мм, по оси $Z$ — 1 мм, размер ячеек в основной сетке второго типа 0.03 мм, всего использовалось около 900 тыс. ячеек, в том числе порядка 100 тыс. ячеек по кварцевым телам.

Расчеты проводились на 4-ядерных компьютерах с процессором Intel I7. Время расчета составило порядка 100 ч.

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

Рис. 4. Поля давлений (Па) в газе в диагональном сечении в различные моменты времени: a — 0.2 мкс, b — 2.5 мкс, c — 10 мкс
[Figure 4. Pressure fields (Pa) in gas in a diagonal section at different times: a — 0.2 μs, b — 2.5 μs, c — 10 μs]

На рис. 5 показано распределение давлений в диагональном сечении вдоль двух линий, параллельных оси $X$, в момент времени 10.4 мкс: кривая 1 проходит через центры шариков, кривая 2 проходит через всю расчетную область газа в диагональном сечении. Разрывы в кривых связаны с тем, что линии проходят и через газ, и через твердые частицы. Амплитуда на фронте выходящей квазиодномерной ударной волны составляет около 0.4 МПа, что примерно на 10 % превышает аналогичные экспериментальные значения [12]. Несмотря на используемый в расчетах на порядок более тонкий гранулированный слой (3 мм), чем в экспериментах 
(20 мм), наблюдается соответствие численных и экспериментальных данных.

Рис. 5. Давление в поровом газе в диагональном сечении в момент времени 10.4 мкс
[Figure 5. Pressure in the pore gas in the diagonal section at a time of 10.4 μs along lines parallel to the $X$ axis: curve 1 passes through the centers of the balls, curve 2 passes through the entire computational domain of the gas]

Рис. 6. Поля давлений (Па) в кварцевых шариках в различные моменты времени: a — 0.2 мкс, b — 2.5 мкс, c — 10 мкс
[Figure 6. Pressure fields (Pa) in quartz balls at different times: a — 0.2 μs, b — 2.5 μs, c — 10 μs]

На рис. 6 показано распределение давления в кварцевых шариках в моменты времени 0.2, 2.5, 10 мкс. Наблюдается сложная картина взаимодействия деформационных волн при контактном взаимодействии шариков между собой, окружающим газом и жесткими границами. Уже на ранних этапах взаимодействия возникают кратковременные отрывы шариков друг от друга. Первый отрыв формируется к моменту времени 10 мкс. В результате схождения упругих волн с оси симметрии шариков образуются локальные области высокого давления до 2 ГПа, но время их существования незначительно, порядка сотых долей микросекунды, не оказывающих заметного влияния на контактное взаимодействие шариков, их формоизменение и газодинамические процессы.

Заключение

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

Конкурирующие интересы. Заявляем, что в отношении авторства и публикации этой статьи конфликта интересов не имеем.
Авторский вклад и ответственность. И. А. Модин—проведение численных расчетов, обработка и анализ результатов, работа с черновиком и переработанным вариантом рукописи. А. В. Кочетков—идея исследования, формулировка целей и задач исследования, визуализация и верификация результатов. Е. Г. Глазова—проведение численных расчетов, работа с черновиком и переработанным вариантом рукописи. А. А. Лисицын—проведение численных расчетов, визуализация и верификация результатов, работа с черновиком и переработанным вариантом рукописи. Авторы несут полную ответственность за предоставление окончательной рукописи в печать. Окончательная версия рукописи была одобрена всеми авторами.
Финансирование. Исследование выполнено за счет гранта Российского научного фонда № 22–29–00672, https://rscf.ru/project/22-29-00672/.

Об авторах

Елена Геннадьевна Глазова

Научно-исследовательский институт механики Нижегородского государственного университета им. Н. И. Лобачевского

Автор, ответственный за переписку.
Email: glazova@mech.unn.ru
ORCID iD: 0000-0003-4351-889X
https://www.mathnet.ru/person163935

кандидат физико-математических наук; ученый секретарь; лаб. динамики многокомпонентных сред

Россия, Россия, 603022, Нижний Новгород, пр. Гагарина, 23, корпус 6

Анатолий Васильевич Кочетков

Научно-исследовательский институт механики Нижегородского государственного университета им. Н. И. Лобачевского

Email: kochetkov@mech.unn.ru
ORCID iD: 0000-0001-7939-8207
https://www.mathnet.ru/person32889

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

Россия, Россия, 603022, Нижний Новгород, пр. Гагарина, 23, корпус 6

Артем Александрович Лисицын

Научно-исследовательский институт механики Нижегородского государственного университета им. Н. И. Лобачевского

Email: artem.lisitsyn@unn.ru
ORCID iD: 0009-0006-9397-6257
https://www.mathnet.ru/person195058

аспирант; младший научный сотрудник; лаб. динамики многокомпонентных сред

Россия, Россия, 603022, Нижний Новгород, пр. Гагарина, 23, корпус 6

Иван Александрович Модин

Научно-исследовательский институт механики Нижегородского государственного университета им. Н. И. Лобачевского

Email: mianet@mail.ru
ORCID iD: 0000-0002-3561-4606
https://www.mathnet.ru/person138504

кандидат технических наук; старший научный сотрудник; лаб. физико-механических испытаний материалов

Россия, Россия, 603022, Нижний Новгород, пр. Гагарина, 23, корпус 6

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

Дополнительные файлы