Impact interaction between granular media and a rod with conical head

Full Text

Введение

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

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

При моделировании пористой среды необходимо учитывать площадь контакта и трение между контактирующими элементами (гранулами). Сыпучие материалы проявляют сложную реакцию при динамической нагрузке: в начальный момент взаимодействия происходит переукладка и уплотнение пористого слоя. Насыпные пористые материалы, свойства которых изменяются в зависимости от степени сжатия, привлекают внимание исследователей из различных областей [7–11].

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

1. Метод дискретных элементов

Метод дискретных элементов (DEM, Discrete Element Method) представляет собой численный инструмент для моделирования поведения пористых, сыпучих и дробимых материалов. В отличие от традиционных континуальных методов, таких как метод конечных элементов (МКЭ) или метод конечных разностей (МКР), DEM рассматривает материал как совокупность дискретных элементов, взаимодействующих посредством контактных сил. Такой подход позволяет учитывать межчастичные контакты, трение, упругие и пластические деформации, а также разрушение материала на уровне отдельных элементов.

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

Основная идея метода заключается в численном решении уравнений движения для каждой частицы с учетом внешних и контактных сил. Динамика системы определяется на основе законов Ньютона, а контактные взаимодействия описываются различными моделями, включая модель Герца–Миндлина и модели липкого контакта. В данной работе используется модель Герца для определения контактных сил и моментов. По своей природе DEM относится к бессеточным лагранжевым методам [12–15] и сводится к явному интегрированию уравнений движения для каждой частицы как твердого тела:
\[\begin{equation}
\tag{1}
m_i\ddot{\boldsymbol{x}}_i = \boldsymbol{F}_{i,n} + \boldsymbol{F}_{i,t} + \boldsymbol{F}_{i,f} + \boldsymbol{F}_{i,b},
\end{equation}\]
\[\begin{equation}
\tag{2}
I_i\frac{d\boldsymbol{\omega}_i}{dt} = \boldsymbol{r}_{i,c} \times \boldsymbol{F}_{i,t} + \boldsymbol{T}_{i,r},
\end{equation}\]
где $\boldsymbol{r}_{i,c}$ — вектор, соединяющий центр масс частицы с точкой контакта; $\boldsymbol{F}_{i,n}$ — нормальная составляющая контактной силы между частицами; $\boldsymbol{F}_{i,t}$ — тангенциальная составляющая контактной силы между частицами; $\boldsymbol{F}_{i,f}$ — сила со стороны окружающей среды (не рассматривается в данной работе); $\boldsymbol{F}_{i,b}$ — суммарное воздействие других центральных сил (гравитация, сила Кулона и др.); $\boldsymbol{T}_{i,r}$ — дополнительный момент сил, учитывающий неидеальность формы частицы.

В простейшей постановке каждая частица характеризуется радиус-вектором центра $\boldsymbol{x}_i$ и угловой скоростью $\boldsymbol{\omega}_i$. Линейные и угловые ускорения вычисляются из уравнений (1) и (2). В DEM частицы рассматриваются как упругие слабодеформируемые тела, где деформация математически описывается как пространственное пересечение частиц. Линейная модель пружинно-демпферного взаимодействия (рис. 1) выражает нормальную составляющую контактной силы через характеристический размер области пересечения $\delta$ и относительную скорость контактирующих частиц $\Delta u$:
\[\begin{equation}
\tag{3}
\boldsymbol{F}_{n} = -k_n\boldsymbol{\delta}_{n} + c_n\Delta\boldsymbol{u}_{n}.
\end{equation}\]

 

Рис. 1. Схема контактного взаимодействия двух частиц в модели дискретных элементов; показаны нормальная ($\boldsymbol{F}_n$) и тангенциальная ($\boldsymbol{F}_t$) составляющие контактной силы
[Figure 1. Schematic of contact interaction between two particles in the discrete element method; shown are the normal ($\boldsymbol{F}_n$) and tangential ($\boldsymbol{F}_t$) components of the contact force]

 

Касательная составляющая силы отталкивания определяется как
\[\begin{equation}
\tag{4}
\boldsymbol{F}_{t} = \min\biggl\{
\biggl|
k_t\int_{t_{c,0}}^t \Delta\boldsymbol{u}_{t}dt + c_t\Delta\boldsymbol{u}_{t}
\biggr|, \; \mu\boldsymbol{F}_{n}
\biggr\},
\end{equation}\]
где $\boldsymbol{\delta}_n = \delta\boldsymbol{n} = \delta(\boldsymbol{x}_i-\boldsymbol{x}_j)/\|\boldsymbol{x}_i-\boldsymbol{x}_j\|$; $\boldsymbol{n}$ — единичный вектор нормали, направленный вдоль линии, соединяющий центры масс частиц; $\Delta\boldsymbol{u}_n$ и $\Delta\boldsymbol{u}_t$ — нормальная и тангенциальная составляющие относительной скорости двух частиц, находящихся в контакте;  $k_n$, $k_t$ и $c_n$, $c_t$ — нормальные и касательные коэффициенты жесткости и демпфирования. Формулы (3) и (4) справедливы и для контакта частицы со стенкой с учетом ограничения кулоновской силой трения с коэффициентом $\mu$.

Интегральный член в (4) соответствует закону Гука для упругой тангенциальной деформации поверхностей двух частиц в контакте, который произошел в момент времени $t=t_{c,0}$, а демпфирующая составляющая учитывает диссипацию упругой энергии. Вариация коэффициентов $k_n$, $k_t$, $c_n$, $c_t$ как функций пространственного пересечения частиц позволяет получать различные линейные и нелинейные модели контактного взаимодействия [16–18].

Здесь и далее индексы $i$ и $j$ относятся к двум соседним взаимодействующим частицам. Для контактной модели Герца, которая используется в численных расчетах, коэффициенты $k_n$, $k_t$ и $c_n$, $c_t$ в случае взаимодействия двух произвольных частиц рассчитываются следующим образом  [13, 17, 20]:
\[\begin{align*}
k_n &= \frac{4}{3}Y^*\sqrt{R^*\delta_n}, \quad
c_n = -2\sqrt{\frac{5}{6}}\beta\sqrt{S_nm^*} \geqslant 0, \\
k_t &= 8G^*\sqrt{R^*\delta_n}, \quad
c_t = -2\sqrt{\frac{5}{6}}\beta\sqrt{S_tm^*} \geqslant 0.
\end{align*}\]

При этом в случае контакта двух частиц с массами $m_1$ и $m_2$, размерами $R_1$ и $R_2$, модулями Юнга $Y_1$ и $Y_2$, коэффициентами Пуассона $\nu_1$ и $\nu_2$:
\[\begin{gather*}
S_n = 2Y^*\sqrt{R^*\delta_n}, \quad
S_t = 8G^*\sqrt{R^*\delta_n}, \quad
\beta = \frac{\ln e}{\sqrt{\ln^2 e + \pi^2}}, 
\\
\frac{1}{Y^*} = \frac{1-\nu^2_1}{Y_1} + \frac{1-\nu^2_2}{Y_2}, \quad
\frac{1}{G^*} = \frac{2(2-\nu_1)(1+\nu_1)}{Y_1} + \frac{2(2-\nu_2)(1+\nu_2)}{Y_2}, \\
\frac{1}{R^*} = \frac{1}{R_1} + \frac{1}{R_2}, \quad
\frac{1}{m^*} = \frac{1}{m_1} + \frac{1}{m_2},
\end{gather*}\]
где $e$ — коэффициент упругого восстановления, который равен отношению скоростей частицы после и до контакта. Как правило, это пользовательская константа, значения которой задаются в диапазоне от 0 до 1.

Численное решение уравнений (1) и (2) получается их последовательным интегрированием с помощью методов трапеций или прямоугольников, в последнем в качестве опорной точки может быть выбрана точка посередине отрезка интегрирования. Для практического внедрения DEM важным аспектом является алгоритм поиска контактов для системы из многих частиц [19, 20]. В 1967 году Верлетом была представлена модель создания списков ближайших соседей. Идея состоит в том, что список потенциальных контактов формируется периодически. Далее на каждом временном шаге этот список проверяется и выполняются вычисления контактных сил. Такой подход позволяет заранее исключать пары частиц, которые находятся далеко друг от друга [21–23]. Две частицы добавляются в список контактов, если
\[
\|\boldsymbol{x}_i - \boldsymbol{x}_j\| < r_i + r_j + s,
\]
где $\boldsymbol{x}_i$, $\boldsymbol{x}_j$ — радиус-векторы центров взаимодействующих частиц; $r_i$, $r_j$ — радиусы взаимодействующих сфер; $s$ — параметр Верлета, или «скин» параметр, который определяет, как долго текущий список контактов может оставаться релевантным. Полагая постоянным временной шаг $\Delta t$, а также, обозначая максимальную скорость частиц как $v_{\max}$, текущий список контактов можно не обновлять следующее количество шагов:
\[
N_{\mathrm{verlet}} = \frac{s}{2v_{\max}\Delta t}.
\]

2. Экспериментальные результаты

С помощью методики обращенного эксперимента проведено исследование ударного взаимодействия тел с пористой средой [24]. Методика предназначена для определения силы сопротивления внедрению ударника путем анализа интегральных нагрузок на начальном участке проникания. Она сводится к следующему: контейнер заданных размеров, заполненный пористым материалом, который разгоняется до требуемой скорости и соударяется с неподвижным коническим индентором (полный угол раствора $60^\circ$), закрепленным на мерном стержне. Материал стержня и конуса подбирается таким образом, чтобы исключить пластические деформации при выбранных скоростях удара.

В ходе ударного взаимодействия в стержне регистрируется упругий импульс сжатия с деформацией $\varepsilon(t)$. Сила сопротивления прониканию $F(t)$, действующая на индентор при взаимодействии с пористой средой, определяется по зарегистрированному импульсу сжатия с использованием соотношения
\[
F(t) = E\varepsilon(t)S,
\]
где $E$ — модуль упругости материала мерного стержня; $S$ — площадь его поперечного сечения. Таким образом, задача измерения сил сопротивления сводится к регистрации продольных упругих деформаций в мерном стержне.

Регистрация продольного импульса осуществлялась датчиками, расположенными на боковой поверхности опорного стержня на расстоянии трех диаметров (60 мм) от места крепления конической части. Конический индентор имел диаметр основания 20 мм с полным углом раствора $60^\circ$. В эксперименте использовался пористый слой, размещенный в полипропиленовом контейнере со следующими параметрами: длина 70 мм, внешний диаметр 56.8 мм, внутренний диаметр 54.5 мм, толщина дна 2 мм. Глубина заполнения контейнера пористым материалом составляла 38 мм при массе 400 г. Длина конической части опорного стержня равнялась 17.3 мм.

Для анализа погрешностей измерений с помощью программного комплекса UPSGOD [40] было выполнено численное моделирование процесса распространения волн деформации от зоны контакта до тензодатчиков, расположенных на расстоянии трех диаметров стержня. Результаты моделирования (рис. 2, кривые 2 и 4) показали отсутствие значительных искажений сигнала на этом расстоянии. На рис. 2 кривые 1 и 3 соответствуют экспериментальным данным о силах сопротивления, зарегистрированным при скоростях соударения 179 м/с и 159 м/с соответственно.

 

Рис. 2. Фиксируемая сила в эксперименте (кривые 1, 3) и при численном расчете (кривые 2, 4) на боковой поверхности мерного стержня на расстоянии трех диаметров
[Figure 2. Measured force in the experiment (curves 1, 3) and in numerical simulation (curves 2, 4) on the lateral surface of the measuring rod at a distance of three diameters]

 

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

На рис. 3 представлена расчетная схема взаимного расположения упругого стержня с коническим контактным элементом и пористой среды. Схема соответствует обращенному эксперименту [24], проводившемуся для оценки силы сопротивления при внедрении опорного стержня путем анализа интегральных нагрузок на начальном участке проникновения. На рис. 3 показана конфигурация пористого слоя и опорного стержня с коническим элементом.

 

Рис. 3. Расчетная схема обращенного эксперимента: взаимодействие конического индентора с гранулированной средой
[Figure 3. Computational scheme of the inverse experiment: interaction between the conical indenter and granular medium]

 

Моделирование ударного взаимодействия осуществлялось с использованием собственного оригинального программного обеспечения, реализующего метод DEM [25]. Параметры пористой упаковки в расчетах:

• модуль Юнга $E = 200$ ГПа; 
• коэффициент Пуассона $\nu = 0.3$;
• объемная доля пор $\phi=0.36$.

Коэффициенты $k_n$, $k_t$ и $c_n$, $c_t$ рассчитывались на каждом временном шаге автоматически с применением контактной модели Герца [25–28] на основе заданных значений модуля Юнга $E$ и коэффициента Пуассона $\nu$. Численное моделирование отдельных частиц рассмотрено в работах [29–34], других пористых структур — в [35–40].

В DEM-моделировании мерный стержень и конический элемент рассматривались как упругие малодеформируемые тела. Расчеты выполнялись с поверхностной сеткой (размер ячейки 0.2 мм для конической части), коэффициент трения между пористым слоем и стержнем принят равным 0.3. Вариации коэффициента трения скольжения не оказывают существенного влияния на результаты. Полипропиленовый контейнер в модели не учитывался, начальная скорость пористой среды составляла 159 м/с и 179 м/с, при этом стержень с коническим элементом находился в состоянии покоя.

Трехмерная постановка задачи в декартовой системе координат $XYZ$ (начало в вершине конуса) предполагала неподвижность, но упругую деформируемость стержня и конуса (рис. 3). Движение гранулированной среды инициировалось вдоль оси $Z$. Контактное взаимодействие между соседними частицами с номерами $i$ и $j$ рассчитывалось в локальных системах координат для каждой сферы гранулированной среды на каждом временном шаге.

На рис. 4,  5 показаны измеренная сила воздействия DEM-частиц на коническую часть стержня вдоль оси $Z$ (кривая 1), усредненная сила (шаг усреднения — 100 временных интервалов, кривая 2) и экспериментальная зависимость (кривая 3).

 

Рис. 4. Расчетные силы нагрузки со стороны пористой среды в сравнении с экспериментом для скорости 159 м/с: 1 — сила воздействия DEM-частиц на конический индентор; 2 — усредненная сила с шагом 100 временных интервалов; 3 — экспериментальные данные
[Figure 4. Calculated loading forces from the porous medium compared with experimental data at 159 m/s: 1 — force exerted by DEM particles on the conical indenter; 2 — averaged force with 100-time-step interval; 3 — experimental measurements]

 

Рис. 5. Расчетные силы нагрузки со стороны пористой среды в сравнении с экспериментом для скорости 179 м/с: 1 — сила воздействия DEM-частиц на конический индентор; 2 — усредненная сила с шагом 100 временных интервалов; 3 — экспериментальные данные
[Figure 5.  Calculated loading forces from the porous medium compared with experimental data at 179 m/s: 1 — force exerted by DEM particles on the conical indenter; 2 — averaged force with 100-time-step interval; 3 — experimental measurements]

 

На рис. 6 и 7 представлены распределения скоростей в упаковке DEM-частиц и ее положения в моменты времени 40, 80 и 100 мкс соответственно для случая метания упаковки со скоростью 159 м/с и 20, 40 и 80 мкс соответственно для случая метания упаковки со скоростью 179 м/с (рис. 7). 

 

Рис. 6. Состояние упаковки частиц и их скорости (м/с), метаемой среды со скоростью 159 м/с, в моменты времени 40, 80 и 100 мкс соответственно
[Figure 6. State of particle packing and velocity distribution (m/s) for the impacted medium at 159 m/s, shown at time instants of 40, 80, and 100 $\mu s$ respectively]

 

Рис. 7. Состояние упаковки частиц и их скорости (м/с), метаемой среды со скоростью 179 м/с, в моменты времени 40, 80 и 100 мкс соответственно
[Figure 7. State of particle packing and velocity distribution (m/s) for the impacted medium at 179 m/s, shown at time instants of 40, 80, and 100 $\mu s$ respectively]

 

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

На рис. 8 представлено трехмерное изображение, иллюстрирующее конфигурацию внедряемого тела и состояние пористого слоя в момент времени $t = 250$ мкс для случая начальной скорости упаковки 159 м/с.

 

Рис. 8. Состояние упаковки пористого слоя, соответствующее моменту времени 250 мкс для случая начальной скорости 159 м/с
[Figure 8. Configuration of the porous layer packing at 250 $\mu s$ for the initial impact velocity of 159 m/s]

 

Заключение

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

Основные результаты исследования могут быть сформулированы следующим образом:

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

Конкурирующие интересы. Авторы заявляют об отсутствии каких-либо конфликтов интересов, связанных с публикацией данной статьи.
Авторский вклад и ответственность. А.В. Кочетков  — концепция исследования, постановка целей и задач, визуализация данных, верификация результатов, работа с черновиком рукописи. И.А. Модин  — концепция исследования, выполнение численного моделирования, подготовка и редактирование рукописи. В.В. Баландин  — планирование и проведение экспериментальных исследований, обработка и интерпретация данных,  работа с черновиком рукописи. Вл.Вл. Баландин  — планирование и проведение экспериментальных исследований, обработка и интерпретация данных, работа с черновиком рукописи. К.Д. Бессмертный  — выполнение численных расчетов, обработка и анализ результатов, работа с черновиком рукописи. Все авторы принимали участие в утверждении окончательной версии рукописи и несут полную ответственность за содержание публикации.
Финансирование. Исследование выполнено при  поддержке Российского научного фонда (проект № 22–79–10076, https://rscf.ru/project/22-79-10076/).
Благодарности. Авторы выражают признательность анонимным рецензентам за внимательное рассмотрение работы и конструктивные замечания, способствовавшие улучшению качества публикации.

About the authors

Anatoliy V. Kochetkov

Research Institute of Mechanics, National Research Lobachevsky State University of Nizhny Novgorod

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

Doctor of Physical and Mathematical Sciences; Head of the Laboratory; Lab. of Dynamics of Multicomponent Media

Russian Federation, 603022, Nizhny Novgorod, pr. Gagarina, 23

Ivan A. Modin

Research Institute of Mechanics, National Research Lobachevsky State University of Nizhny Novgorod

Author for correspondence.
Email: mianet@mail.ru
ORCID iD: 0000-0002-3561-4606
Scopus Author ID: 57192279101
ResearcherId: E-9088-2019
https://www.mathnet.ru/rus/person138504

Candidate of Technical Sciences; Senior Researcher; Lab. of Physical and Mechanical Testing of Materials

Russian Federation, 603022, Nizhny Novgorod, pr. Gagarina, 23

Vladimir V. Balandin

Research Institute of Mechanics, National Research Lobachevsky State University of Nizhny Novgorod

Email: vbalandin99@gmail.com
ORCID iD: 0000-0002-3907-3480
https://www.mathnet.ru/rus/person189695

Candidate of Physical and Mathematical Sciences; Leading Researcher; Lab. of Multicomponent Media Dynamics

Russian Federation, 603022, Nizhny Novgorod, pr. Gagarina, 23

Vladimir Vl. Balandin

Research Institute of Mechanics, National Research Lobachevsky State University of Nizhny Novgorod

Email: rustydog2007@yandex.ru
ORCID iD: 0000-0003-3061-8875
https://www.mathnet.ru/rus/person189694

Candidate of Technical Sciences; Leading Researcher; Lab. of Dynamic Testing of Materials

Russian Federation, 603022, Nizhny Novgorod, pr. Gagarina, 23

Kirill D. Bessmertnyi

Research Institute of Mechanics, National Research Lobachevsky State University of Nizhny Novgorod

Email: besskirill@gmail.com
ORCID iD: 0009-0003-0515-4375
https://www.mathnet.ru/rus/person231443

Junior Researcher; Lab. of Multicomponent Media Dynamics

Russian Federation, 603022, Nizhny Novgorod, pr. Gagarina, 23

References

Supplementary files