Континуальная модель перидинамики для задач хрупкого разрушения

1.      Введение

Статья является продолжением исследования метода перидинамики [1$–$2] $–$ альтернативного подхода для решения задач механики твердых тел. Формулировка перидинамики построена на использовании интегральных уравнений, избегая пространственных производных [3]. Область континуума задается в виде системы взаимодействующих между собой частиц с заданным объемом. Взаимодействие частиц осуществляется посредством связи в пределах замкнутого горизонта $–$ масштаба длины $\delta$.

Non-ordinary state-based (NOSB) модель перидинамики [4] является обобщением подходов, рассмотренных в [1]. Она использует формулировку классической теории и перидинамики, которая позволяет избавиться от ограничений, свойственных моделям bond-based (BB) и ordinary state-based (OSB) [1$–$2]. Зависимость вектора силы частицы от тензоров деформаций и напряжений позволяет задавать любое направление действия силы, что гораздо более реалистично для моделирования континуума.

В нелокальной теории точность вычисления тензоров зависит от расстояния между частицами $\Delta x$ и горизонта взаимодействия $\delta$. Согласно теории сходимости [5], уравнения перидинамики превращаются в уравнения механики сплошной среды (МСС) при уменьшении $\delta$. Предельные численные решения этих уравнений также будут совпадать при увеличении числа частиц дискретной модели и измельчении сетки в МСС. В задачах разрушения для получения корректного результата моделирования при минимально возможных вычислительных затратах рекомендуется использовать $\delta \approx 3\Delta x$ [6].

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

2. Метод перидинамики

Теория перидинамики была введена Силингом [3] как способ решения задач механики твердого тела без использования пространственных производных, основываясь на интегральных уравнениях. Твердое тело представляется в виде системы взаимодействующих между собой частиц с заданными массой $—$ $m_i$, плотностью $—$ ${\rho }_i$, начальными координатами $—$ ${\overrightarrow{r}}_i^0\left(x_i^0\mathrm{,}{y}_i^0\mathrm{,}{z}_i^0\right)$ ). Каждая дискретная частица ${\overrightarrow{r}}_i$ взаимодействует с другими частицами ${\overrightarrow{r}}_k$ на конечном расстоянии ${\overrightarrow{\xi }}_i\mathrm{<mo>=</mo>}{\overrightarrow{r}}_k-{\overrightarrow{r}}_i$ посредством связи в пределах замкнутого горизонта $—$ $\delta$. Движение деформируемой среды описывается векторами относительных смещений ${\overrightarrow{\eta }}_i\mathrm{<mo>=</mo>}{\overrightarrow{u}}_k-{\overrightarrow{u}}_i$ и векторами относительных положений ${\overrightarrow{Y}}_i\mathrm{<mo>=</mo>}{\overrightarrow{\xi }}_i+{\overrightarrow{\eta }}_i$.

Уравнение движения частицы ${\overrightarrow{r}}_i\left(x_i\mathrm{,}{y}_i\mathrm{,}{z}_i\right)$ задается в виде

                                           ${\rho }_i\ddot{\overrightarrow{u}}\left({\overrightarrow{r}}_i\mathrm{,}t\right)\mathrm{<mo>=</mo>}{\displaystyle {\int }_{{\mathcal{H}}_{{\overrightarrow{r}}_i}}}{\overrightarrow{f}}_i\left({\overrightarrow{u}}_k-{\overrightarrow{u}}_i\mathrm{,}{\overrightarrow{r}}_k-{\overrightarrow{r}}_i\right)d{V}_{{\overrightarrow{r}}_k}+\overrightarrow{b}\left({\overrightarrow{r}}_i\mathrm{,}t\right)\mathrm{,}$ (2.1)

 где $\overrightarrow{b}$ $—$ граничные условия.

Выделяют две основные формулировки перидинамики $—$ на основе связи (bond-based) [3, 8] и на основе состояний (state-based) [4, 9] . Модель на основе состояний в свою очередь подразделяется на ordinary state-based (OSB) и NOSB. NOSB модель обобщает рассмотренные ранее PMB и LPS модели [1, 2]. Более подробно теория перидинамики представлена в [10, 11] .

3. NOSB модель

Рассматриваемая модель основана на вычислении вектора силы каждой частицы в связи [7, 12, 13] . Описание модели начинается с задания энергии деформирования, как и в МСС. Для упругих деформаций:

                                                           $W_i\mathrm{<mo>=</mo>}\frac{{\lambda }_i}{2}{\displaystyle \sum _k}{E}_{kk}^i{E}_{kk}^i+G_i{\displaystyle \sum _{m\mathrm{,}l}}{E}_{ml}{E}_{ml}\mathrm{,}$ (3.1)

 где ${\lambda }_i\mathrm{<mo>=</mo>}{K}_i-\frac{2}{3}{G}_i$ $—$ коэффициенты Ламэ; $K_i$ $—$ объемный модуль упругости; $G_i$ $—$ модуль сдвига; $E_i$ $—$ упругий тензор деформаций Грина:

                                                                         $E_i\mathrm{<mo>=</mo>}\frac{1}{2}\left(F_i^T{F}_i-I\right)\mathrm{,}$ (3.2)

 где $F_i^T{F}_i$ $—$ упругий тензор деформаций Коши-Грина; $I$ $—$ единичная матрица; $F_i$ $—$ тензор градиента деформаций:

                                                  $F_i\mathrm{<mo>=</mo>}\left({\displaystyle \underset{{\mathcal{H}}_{{\overrightarrow{r}}_i}}{\int }}\omega \left(\mathrm{<mo>|</mo>}{\overrightarrow{\xi }}_i\mathrm{<mo>|</mo>}\right)\left({\overrightarrow{\xi }}_i+{\overrightarrow{\eta }}_i\right)\otimes {\overrightarrow{\xi }}_{i}d{V}_{{\overrightarrow{r}}_k}\right)\cdot K_i^{-1}\mathrm{,}$ (3.3)

 где $\omega \left(\mathrm{<mo>|</mo>}{\overrightarrow{\xi }}_i\mathrm{<mo>|</mo>}\right)$ $—$ функция влияния в виде кубического сплайна 6; $K_i$ $—$ шаровой тензор деформаций:

                                                             $K_i\mathrm{<mo>=</mo>}{\displaystyle \underset{{\mathcal{H}}_{{\overrightarrow{r}}_i}}{\int }}\omega \left(\mathrm{<mo>|</mo>}{\overrightarrow{\xi }}_i\mathrm{<mo>|</mo>}\right)\left({\overrightarrow{\xi }}_i\otimes {\overrightarrow{\xi }}_i\right)d{V}_{{\overrightarrow{r}}_k}\mathrm{,}$ (3.4)

                                  $\omega \left(\mathrm{<mo>|</mo>}{\overrightarrow{\xi }}_i\mathrm{<mo>|</mo>}\right)\mathrm{<mo>=</mo>}\left\{\begin{array}{rrr}\hfill \frac{2}{3}-4d^2+4d^3\mathrm{,}& \hfill d\le \frac{1}{2}\mathrm{,}& \hfill \\ \hfill \frac{4}{3}-4d+4d^2-\frac{4}{3}{d}^3\mathrm{,}& \hfill \frac{1}{2}\mathrm{<mo><</mo>}d\le \mathrm{1,}& \hfill d\mathrm{<mo>=</mo>}\frac{\mathrm{<mo>|</mo>}{\overrightarrow{\xi }}_i\mathrm{<mo>|</mo>}}{\delta }\mathrm{,}\\ \hfill \mathrm{0,}& \hfill d\mathrm{<mo>></mo>1,}& \hfill \end{array}\right.$ (3.5)

 Тензор напряжений для системы координат начального состояния (второй тензор напряжений Пиолы-Кирхгофа) определяется по закону Гука 3.6. Симметричный тензор напряжений в деформированном состоянии (тензор напряжений Коши) записывается в виде 3.7.

                                                          $S_i\mathrm{<mo>=</mo>}\frac{\partial W_i}{\partial E_i}\mathrm{<mo>=</mo>}{\lambda }_{i}Tr\mathrm{<mo\; stretchy="false">(</mo>}{E}_i\mathrm{<mo\; stretchy="false">)</mo>}{g}_{kl}+2G_i{E}_i\mathrm{,}$ (3.6)

 где $g_{kl}\mathrm{<mo>=</mo>}\left\{\begin{array}{rrr}\hfill \mathrm{1,}& \hfill k\mathrm{<mo>=</mo>}l& \hfill \\ \hfill \mathrm{0,}& \hfill k\ne l& \hfill \end{array}\right.$ $—$ дельта-функция,

                                                                  ${\sigma }_i\mathrm{<mo>=</mo>}{F}_i\left(\frac{1}{det\left(F_i\right)}{S}_i\right)F_i^T\mathrm{.}$ (3.7)

 Вычисление ${\sigma }_i$ относительно начальной конфигурации (3.2), (3.6), (3.7) используется для описания малых упругих деформаций [13]. Преобразование напряжения из начальной конфигурации в деформируемую конфигурацию некорректно описывает поведение материала при хрупком разрушении. В данной работе для описания кинематики частицы используется тензор градиента скорости деформаций [14]. Тензор напряжений Коши записывается в виде (3.11).

Тензор градиента скорости деформаций задается в виде:

                                                                                $L_i\mathrm{<mo>=</mo>}{\dot{F}}_i{F}^{-1}\mathrm{,}$ (3.8)

 где ${\dot{F}}_i$ $—$ производная тензора градиента деформаций по времени:

                                                         ${\dot{F}}_i\mathrm{<mo>=</mo>}\left({\displaystyle \underset{{\mathcal{H}}_{{\overrightarrow{r}}_i}}{\int }}\omega \left(\mathrm{<mo>|</mo>}{\overrightarrow{\xi }}_i\mathrm{<mo>|</mo>}\right){\dot{\overrightarrow{Y}}}_i\otimes {\overrightarrow{\xi }}_{i}d{V}_{{\overrightarrow{r}}_k}\right)\cdot K_i^{-1}\mathrm{.}$ (3.9)

 Симметричная часть тензора градиента скорости описывает скорость тензора деформаций:

                                                                           $d_i\mathrm{<mo>=</mo>}\frac{1}{2}\left(L_i+L_i^T\right)\mathrm{.}$ (3.10)

 Тензор напряжений Коши обновляется на каждом временном шаге:

                                                                          ${\sigma }_i^{n+1}\mathrm{<mo>=</mo>}{\sigma }_i^n+\Delta {\sigma }_i\mathrm{,}$ (3.11)

 где $\Delta {\sigma }_i$ $—$ тензор приращения напряжений 13.

Тензор приращения напряжений $\Delta {\sigma }_i$ (3.12) получается в зависимости от тензора приращения деформаций $\Delta E_i$ (3.13).

                                                             $\Delta {\sigma }_i\mathrm{<mo>=</mo>}{\lambda }_{i}Tr\left(\Delta E_i\right)g_{kl}+2G_i\Delta E_i$ (3.12)

                                                                               $\Delta E_i\mathrm{<mo>=</mo>}{d}_i\cdot \Delta t$ (3.13)

 Подставляя (3.8) и (3.10) в (3.13), получаем (3.14) $–$ зависимость $\Delta E_i$ от тензора приращения градиента деформаций $\Delta F_i$ (3.15).

                                                             $\Delta E_i\mathrm{<mo>=</mo>}\frac{1}{2}\left(\Delta F_i{F}_i^{-1}+F_i^{-T}\Delta F_i^T\right)\mathrm{,}$ (3.14)

                                                     $\Delta F_i\mathrm{<mo>=</mo>}\left({\displaystyle \underset{{\mathcal{H}}_{{\overrightarrow{r}}_i}}{\int }}\omega \left(\mathrm{<mo>|</mo>}{\overrightarrow{\xi }}_i\mathrm{<mo>|</mo>}\right)\Delta {\overrightarrow{Y}}_i\otimes {\overrightarrow{\xi }}_{i}d{V}_{{\overrightarrow{r}}_k}\right)\cdot K_i^{-1}\mathrm{,}$ (3.15)

 где $\Delta {\overrightarrow{Y}}_i\mathrm{<mo>=</mo>}\Delta {\overrightarrow{u}}_k-\Delta {\overrightarrow{u}}_i$ $—$ приращение относительного положения частицы $i$, $\Delta {\overrightarrow{u}}_i\mathrm{<mo>=</mo>}$ $\mathrm{<mo>=</mo>}\Delta {\overrightarrow{u}}_i^{n+1}-\Delta {\overrightarrow{u}}_i^n$ $—$ приращение смещения частицы $i$.

Вектор силы ${\overrightarrow{T}}_i$ для каждой частицы в связи задается в виде:

                                                                       ${\overrightarrow{T}}_i\mathrm{<mo>=</mo>}\omega \left(\mathrm{<mo>|</mo>}{\overrightarrow{\xi }}_i\mathrm{<mo>|</mo>}\right)P_i{K}_i^{-1}{\overrightarrow{\xi }}_i\mathrm{,}$ (3.16)

 где $P_i$ $—$ первый тензор напряжений Пиолы-Кирхгофа:

                                                               ${\overrightarrow{P}}_i\mathrm{<mo>=</mo>}J{\sigma }_i{F}_i^{-T}\mathrm{,}J\mathrm{<mo>=</mo>}det\left(F_i\right)\mathrm{.}$ (3.17)

 Тогда сила парного взаимодействия в связи ${\overrightarrow{f}}_i$ для (2.1) вычисляется как

                                                   ${\overrightarrow{f}}_i\mathrm{<mo>=</mo>}\omega \left(\mathrm{<mo>|</mo>}{\overrightarrow{\xi }}_i\mathrm{<mo>|</mo>}\right)P_i{K}_i^{-1}{\overrightarrow{\xi }}_i-\omega \left(\mathrm{<mo>|</mo>}{\overrightarrow{\xi }}_k\mathrm{<mo>|</mo>}\right)P_k{K}_k^{-1}{\overrightarrow{\xi }}_k\mathrm{,}$ (3.18)

 где $i$ и $k$ $—$ индексы данной частицы ${\overrightarrow{r}}_i$ и частицы ${\overrightarrow{r}}_k$ из множества соседей ${\mathcal{H}}_{{\overrightarrow{r}}_i}$ соответственно.

4. Объемы дискретных частиц

При численном интегрировании, уравнения движения (2.1) аппроксимируются дискретными объемами частиц в области ${\mathcal{H}}_{{\overrightarrow{r}}_i}$:

                                                  ${\rho }_i{\overrightarrow{\ddot{u}}}_i^n\mathrm{<mo>=</mo>}{\displaystyle \sum _{k\in {\mathcal{H}}_{{\overrightarrow{r}}_i}}}{\overrightarrow{f}}_i\left({\overrightarrow{u}}_k^n-{\overrightarrow{u}}_i^n\mathrm{,}{\overrightarrow{r}}_k^0-{\overrightarrow{r}}_i^0\right)\Delta V_k+{\overrightarrow{b}}_i^n$  (4.1)

 В (4.1) суммирование идет по частицам, представленным на рисунке 4.1a. В МСС область интегрирования сплошная, как показано на рисунке 4.1b. Для улучшения аппроксимации в данной работе используется корректировка объемов частиц, описанная в справочниках по перидинамике [10, 11]:

                                                                         $\Delta V_k\mathrm{<mo>=</mo>}k\left({\overrightarrow{\xi }}_i\right)\cdot \Delta x^3\mathrm{,}$ (4.2)

 где $k\left({\overrightarrow{\xi }}_i\right)$ $—$ линейная корректирующая функция:

                                            $k\mathrm{<mo\; stretchy="false">(</mo>}{\overrightarrow{\xi }}_i\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}\left\{\begin{array}{rrr}\hfill \frac{\delta +\frac{\Delta x}{2}-\mathrm{<mo>|</mo>}{\overrightarrow{\xi }}_i\mathrm{<mo>|</mo>}}{\Delta x}\mathrm{,}& \hfill \delta -\frac{\Delta x}{2}\le \mathrm{<mo>|</mo>}{\overrightarrow{\xi }}_i\mathrm{<mo>|</mo>}\le \delta \mathrm{,}& \hfill \\ \hfill \mathrm{1,}& \hfill \mathrm{<mo>|</mo>}{\overrightarrow{\xi }}_i\mathrm{<mo>|</mo>}\le \delta -\frac{\Delta x}{2}\mathrm{,}& \hfill \\ \hfill \mathrm{0,}& \hfill \mathrm{<mo>|</mo>}{\overrightarrow{\xi }}_i\mathrm{<mo>|</mo><mo>></mo>}\delta \mathrm{,}& \hfill \\ \hfill & \hfill & \hfill \end{array}\right.$ (4.3)

Рис. 4.1. Объемы дискретных частиц в ${\mathcal{H}}_{{\overrightarrow{r}}_i}$: $a$) $—$ перидинамика; $b$) $—$ МСС

Fig 4.1. Discret particles volumes in ${\mathcal{H}}_{{\overrightarrow{r}}_i}$: $a$) $—$ peridynamics; $b$) $—$ continuum mechanics

5. Повреждение и критерий разрушения

Локальное повреждение (5.1) задается в зависимости от числа разорванных связей каждой частицы в пределах ее горизонта взаимодействия и принимает значения в диапазоне $0\le {\varphi }_i\le 1$, где 0 означает, что материал целый, а 1 означает завершенный разрыв связей частицы со всеми частицами, с которыми она изначально взаимодействовала [7]:

                                                                 ${\varphi }_i\mathrm{<mo>=</mo>1}-\frac{{\displaystyle \underset{{\mathcal{H}}_{{\overrightarrow{r}}_i}}{\int }}\mu \left({\overrightarrow{r}}_i\mathrm{,}t\mathrm{,}{\overrightarrow{\xi }}_i\right)d{V}_{{\overrightarrow{r}}_k}}{{\displaystyle \underset{{\mathcal{H}}_{{\overrightarrow{r}}_i}}{\int }}d{V}_{{\overrightarrow{r}}_k}}\mathrm{,}$ (5.1)

 где $\mu$ $—$ скалярная функция положения материальной точки, которая принимает значения $1$ или $0$, определяется формулой

                                                  $\mu \left({\overrightarrow{r}}_i\mathrm{,}t\mathrm{,}{\overrightarrow{\xi }}_i\right)\mathrm{<mo>=</mo>}\left\{\begin{array}{l}\mathrm{1,}{\sigma }_{ij}^I\mathrm{<mo><</mo>}{f}_t\forall 0\le t^{\mathrm{<mo>\text{'}</mo>}}\le t\\ \mathrm{0,}{\sigma }_{ij}^I\ge f_t\forall 0\le t^{\mathrm{<mo>\text{'}</mo>}}\le t\mathrm{,}\end{array}\right.$ (5.2)

 где $f_t$ $—$ критическое значение напряжения связи, пока предполагается константой; ${\sigma }_{ij}^I$ $—$ максимальное главное напряжение (5.4), полученное из тензора:

                                                                          ${\sigma }_{ij}\mathrm{<mo>=</mo>}\frac{1}{2}\left({\sigma }_i+{\sigma }_j\right)\mathrm{,}$ (5.3)

 где $j$ $—$ индекс частицы из ${\mathcal{H}}_{{\overrightarrow{r}}_i}$.

                                             ${\sigma }_{ij}^I\mathrm{<mo>=</mo>}\frac{1}{2}\left({\sigma }_{ij}^{11}+{\sigma }_{ij}^{22}\right)+\sqrt{{\left(\frac{{\sigma }_{ij}^{11}-{\sigma }_{ij}^{22}}{2}\right)}^2+{\left({\sigma }_{ij}^{12}\right)}^2}\mathrm{,}$ (5.4)

 где ${\sigma }_{ij}^{11}\mathrm{,}{\sigma }_{ij}^{22}\mathrm{,}{\sigma }_{ij}^{12}$ $—$ компоненты тензора ${\sigma }_{ij}$.

Частица ${\overrightarrow{r}}_i$ внутри ${\mathcal{H}}_{{\overrightarrow{r}}_i}$ взаимодействует со всеми своими соседями ${\overrightarrow{r}}_k$ до момента разрыва связи, как показано на рисунке 5.1 $a$. Вектор сил двух взаимодействующих частиц в связи вычисляется в зависимости от градиента деформаций, который является функцией относительных смещений. При достижении максимальным главным напряжением ${\sigma }_{ij}^I$ критического значения напряжения $f_t$, взаимодействие прекращается, как показано на рисунке 5.1 $b$. После этого тензор градиента деформаций вычисляется без учета разорванной связи.

Рис. 5.1. Схема взаимодействующих частиц: $a$) $—$ нет разорванных связей; $b$) $—$ есть разорванная связь

Fig 5.1. Scheme of interacting particles: $a$) $—$ there are no broken bonds; $b$) $—$ there is a broken bond

6. Тестовые расчеты

Перидинамика относится к числу нелокальных методов с масштабом длины $\delta$. Точность вычисления тензоров зависит от способа дискретизации области континуума и от выбора горизонта взаимодействия. Согласно теории сходимости [5], уравнения перидинамики превращаются в уравнения МСС при уменьшении $\delta$, а также при увеличении числа частиц дискретной модели и измельчении сетки в МСС. В задачах разрушения для получения корректного результата моделирования при минимально возможных вычислительных затратах рекомендуется использовать $\delta \approx 3\Delta x$ [6].

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

Одноосное растяжение тонкого стержня в 2D-постановке Тонкий стержень подвергается растяжению с левого конца при фиксированном правом. Начальная постановка задачи и аналитическое решение взяты из [13]. Используются безразмерные величины.

Геометрия задачи: $\mathrm{0<mo><</mo>}x\mathrm{<mo><</mo>}L\mathrm{,}-\frac{H}{2}\mathrm{<mo><</mo>}y\mathrm{<mo><</mo>}\frac{H}{2}\mathrm{,}L\mathrm{<mo>=</mo>100,}H\mathrm{<mo>=</mo>10}$.

Свойства материала: модуль Юнга $\text{Е}\mathrm{<mo>=</mo>100}$, коэффициент Пуассона $\nu \mathrm{<mo>=</mo>0}$, плотность $\rho \mathrm{<mo>=</mo>1}$.

Граничные условия

                                                  $w_x\left(x\mathrm{<mo>=</mo>0,}t\right)\mathrm{<mo>=</mo>0.0001,}{w}_x\left(x\mathrm{<mo>=</mo>}L\mathrm{,}t\right)\mathrm{<mo>=</mo>0.}$ (7.1)

 применяются для одного слоя частиц с каждого конца стержня как показано на рисунке 7.1. К левому концу прикладывается постоянное давление $P\mathrm{<mo>=</mo>}-1$.

Рис. 7.1. Геометрия задачи и граничные условия

Fig 7.1. Problem geometry and boundary conditions

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

В расчетах исследуется сходимость реализованной модели для различных значений числа частиц и горизонта взаимодействия. Используются графики смещения и напряжения вдоль оси $Ox$ на начальный момент времени (рис. 7.2, 7.5), смещения и напряжения частицы в сечении $x\mathrm{<mo>=</mo>}L\mathrm{<mo>/</mo>2}$ на всем временном интервале (рис. 7.3, 7.6), а также диаграммы среднеквадратичного отклонения смещения от аналитического решения (рис. 7.4, 7.7).

В первой части расчетов варьируется горизонт взаимодействия $\delta \mathrm{<mo>=</mo>3.015}\Delta x$, $\delta \mathrm{<mo>=</mo>}$ $\mathrm{<mo>=</mo>2.015}\Delta x$, $\delta \mathrm{<mo>=</mo>1.015}\Delta x$ при фиксированном числе частиц $N\mathrm{<mo>=</mo>25551}\mathrm{<mo\; stretchy="false">(</mo>}\Delta x\mathrm{<mo>=</mo>0.2<mo\; stretchy="false">)</mo>}$. Сходимость достигается при уменьшении горизонта взаимодействия (Рис. 7.5 $–$ 7.7).

Рис. 7.2. Смещение ( $a\mathrm{,}b\mathrm{,}c$ ) и напряжение ( $d$ ) вдоль оси $Ox$ в момент времени $t\mathrm{<mo>=</mo>4}$ для различных $\delta$.

Fig 7.2 Displacement ( $a\mathrm{,}b\mathrm{,}c$ ) and stress ( $d$ ) along the $Ox$ axis at time $t\mathrm{<mo>=</mo>4}$ for different $\delta$

Рис. 7.3. Зависимость смещения ( $a\mathrm{,}b\mathrm{,}c$ ) и напряжения ( $d$ ) частицы в сечении $x\mathrm{<mo>=</mo>}L\mathrm{<mo>/</mo>2}$ от времени для различных $\delta$ 

Fig 7.3. The dependence of the displacement ( $a\mathrm{,}b\mathrm{,}c$ ) and stress ( $d$ ) of the particle in the cross section $x\mathrm{<mo>=</mo>}L\mathrm{<mo>/</mo>2}$ on time for different $\delta$

Рис. 7.4. График среднеквадратичного отклонения смещения частиц от эталона на времена для различных $\delta$: 1 $—$ $\delta \mathrm{<mo>=</mo>1.015}\Delta x$, 2 $—$ $\delta \mathrm{<mo>=</mo>2.015}\Delta x$, 3 $—$ $\delta \mathrm{<mo>=</mo>3.015}\Delta x$

Fig 7.4. Mean-square deviation graph of particle diaplacement from analytical solution for times for different $\delta$: 1 $—$ $\delta \mathrm{<mo>=</mo>1.015}\Delta x$, 2 $—$ $\delta \mathrm{<mo>=</mo>2.015}\Delta x$, 3 $—$ $\delta \mathrm{<mo>=</mo>3.015}\Delta x$

Рис. 7.5. Смещение ( $a\mathrm{,}b\mathrm{,}c$ ) и напряжение ( $d$ ) вдоль оси $Ox$ в момент времени $t\mathrm{<mo>=</mo>4}$ для различного числа частиц

Fig 7.5. Displacement ( $a\mathrm{,}b\mathrm{,}c$ ) and stress ( $d$ ) along the $Ox$ axis at time $t\mathrm{<mo>=</mo>4}$ for different particles number

Рис. 7.6. Зависимость смещения ( $a\mathrm{,}b\mathrm{,}c$ ) и напряжения ( $d$ ) частицы в сечении $x\mathrm{<mo>=</mo>}L\mathrm{<mo>/</mo>2}$ от времени для различного числа частиц.

Fig 7.6. The dependence of the displacement ( $a\mathrm{,}b\mathrm{,}c$ ) and stress ( $d$ ) of the particle in the cross section $x\mathrm{<mo>=</mo>}L\mathrm{<mo>/</mo>2}$ on time for different particles number

Рис. 7.7. График среднеквадратичного отклонения смещения частиц от эталона на времена для различного числа частиц: 1 - $N\mathrm{<mo>=</mo>1111}$, 2 - $N\mathrm{<mo>=</mo>4223}$, 3 - $N\mathrm{<mo>=</mo>16441}$, 4 - $N\mathrm{<mo>=</mo>25551}$

Fig 7.7. Mean-square deviation graph of particle diaplacement from analytical solution for times for different particles number: 1 - $N\mathrm{<mo>=</mo>1111}$, 2 - $N\mathrm{<mo>=</mo>4223}$, 3 - $N\mathrm{<mo>=</mo>16441}$, 4 - $N\mathrm{<mo>=</mo>25551}$

Во второй части расчетов варьируется число частиц $N\mathrm{<mo>=</mo>1111}$ $\mathrm{<mo\; stretchy="false">(</mo>}\Delta x\mathrm{<mo>=</mo>1<mo\; stretchy="false">)</mo>}$, N = $4221$ $\mathrm{<mo\; stretchy="false">(</mo>}\Delta x\mathrm{<mo>=</mo>0.5<mo\; stretchy="false">)</mo>,}N\mathrm{<mo>=</mo>16441<mo\; stretchy="false">(</mo>}\Delta x\mathrm{<mo>=</mo>0.25<mo\; stretchy="false">)</mo>,}N\mathrm{<mo>=</mo>25551<mo\; stretchy="false">(</mo>}\Delta x\mathrm{<mo>=</mo>0.2<mo\; stretchy="false">)</mo>}$ при фиксированном значении $\delta \mathrm{<mo>=</mo>3.015}\Delta x$. Сходимость наблюдается при увеличении числа частиц.

8. Разрушение пластины с горизонтальной трещиной при одноосной нагрузке

Прямоугольная пластина с горизонтальным дефектом подвергается постоянной симметричной нагрузке как показано на рисунке 10. В качестве хрупкого материала используется стекло. Постановка задачи взята из [7].

Рис. 8.1. Геометрия и условия нагрузки для прямоугольной пластины

Fig 8.1. Geometry and load conditions for a rectangular plate

Геометрия задачи: $-\frac{L}{2}\mathrm{<mo><</mo>}x\mathrm{<mo><</mo>}\frac{L}{2}$, $-\frac{h}{2}\mathrm{<mo><</mo>}y\mathrm{<mo><</mo>}\frac{h}{2}$, $L\mathrm{<mo>=</mo>0.1}\text{м}$, $h\mathrm{<mo>=</mo>0.04}\text{м}$, $d\mathrm{<mo>=</mo>0.05}\text{м}$.

Свойства материала: модуль Юнга $\text{Е}\mathrm{<mo>=</mo>65}\text{ГПа}$, коэффициент Пуассона $\nu \mathrm{<mo>=</mo>0.2}$, плотность $\rho \mathrm{<mo>=</mo>2}235\text{кг}\mathrm{<mo>/</mo>}{\text{м}}^3$.

В качестве критерия разрушения используется критическое значение напряжения при растяжении $f_t\mathrm{<mo>=</mo>10}$ МПа. Дефект задается как локальное повреждение по формуле (5.1) с заданным числом разорванных связей. Начальное значение прикладываемой нагрузки ${\sigma }_0\mathrm{<mo>=</mo>1}$ МПа, горизонт взаимодействия $\delta \mathrm{<mo>=</mo>3.015}\Delta x$, шаг сетки $\Delta x\mathrm{<mo>=</mo>5}\cdot {10}^{-4}$ м, число частиц $N\mathrm{<mo>=</mo>16}000$, счетный шаг $\Delta t\mathrm{<mo>=</mo>0.01}$ мкс, конец счета при $t\mathrm{<mo>=</mo>60}$ мкс.

Результаты расчета представлены в виде растровых картин разрушения и напряжения пластины, полученные в ParaView 5.9.0. На Рис. 8.2 $a$ видно, что предельное значение напряжения достигается в центре пластины в области начального дефекта на момент времени 7.6 мкс, где происходит зарождение трещины, то есть первый разрыв связей под нагрузкой в области концентрации максимальных главных напряжений. При дальнейшей нагрузке повреждение развивается горизонтально вдоль прямой линии (Рис. 8.2 $b$ ) до момента начала ветвления (Рис. 8.2 $c$ ). Образовавшиеся ветви трещины распространяются симметрично относительно оси $Ox$ до конца счета задачи (Рис. 8.2 $d$ ). Качественное сравнение результата разрушения представлено на рисунке 8.3. Наблюдается адекватное поведение эволюции разрушения с сохранением симметрии ветвления. На рисунке 8.4 показана скорость образования трещины в сравнении с данными из открытых источников. Видно, что начало зарождения разрушения совпадает с результатом, полученным по методу конечных элементов [16], однако дальнейшее распространение трещины протекает медленнее. С этим связано и разное конечное время счета задачи (60 мкс по NOSB модели и 45 мкс по методу конечных элементов).

Рис. 8.2. Растровые картины максимального главного напряжения (первый столбец, Па) и повреждения (второй столбец) на времена: $a$ ) $t\mathrm{<mo>=</mo>7.6}$ мкс; $b$ ) $t\mathrm{<mo>=</mo>20}$ мкс; $c$ ) $t\mathrm{<mo>=</mo>29}$ мкс; $d$ ) $t\mathrm{<mo>=</mo>60}$ мкс

Fig 8.2. Raster pictures of the maximum main stress (first column, Pa) and damage (second column) at times: $a$ ) $t\mathrm{<mo>=</mo>7.6}\mu s$; $b$ ) $t\mathrm{<mo>=</mo>20}\mu s$; $c$ ) $t\mathrm{<mo>=</mo>29}\mu s$; $d$ ) $t\mathrm{<mo>=</mo>60}\mu s$

Рис. 8.3. Сравнение картины ветвления трещины на конечный момент времени: $a$ ) $—$ эксперимент [15]; $b$ ) $—$ численная модель из [7]; $c$ ) $—$ метод конечных элементов XFEM из [16]; $d$ ) $—$ bond-based модель из [17]; $e$ ) $—$ NOSB модель

Fig 8.3. Comparison of the crack branching pattern at a finite point in time: $a$ ) $—$ experiment [15]; $b$ ) $—$ numerical model in [7]; $c$ ) $—$ XFEM finite element method in [16]; $d$ ) $—$ bond-based model in [17]; $e$ ) $—$ NOSB model

 

Рис. 8.4. Скорость эволюции трещины: 1 $—$ эксперимент [15]; 2 $—$ численная модель из [7]; 3 $—$ метод конечных элементов XFEM [16]; 4 $—$ bond-based модель из [17]; 5 $—$ NOSB модель

Fig 8.4. The rate of crack evolution: 1 $—$ experiment [15]; 2 $—$ numerical model in [7]; 3 $—$ XFEM finite element method in [16]; 4 $—$ bond-based model из [17]; 5 $—$ NOSB model

9. Заключение

 В результате проделанной работы реализована и протестирована NOSB модель перидинамики. Показаны сходимость численного решения упругой задачи к аналитике и возможность реализованной модели описывать хрупкое разрушение с применением критерия максимального напряжения при растяжении. Использование большого горизонта взаимодействия $\delta \approx 3\Delta x$ для задачи с дефектом с одной стороны позволяет адекватно моделировать процесс разрушения, с другой стороны вносит погрешность численного интегрирования при вычислении тензора градиента деформаций.