Dual-variable-horizon peridynamics and continuum mechanics coupling modeling and adaptive fracture simulation in porous materials

Abstract

In this paper, we present a hybrid dual-variable-horizon peridynamics/continuum mechanics modeling approach and a strength-induced adaptive coupling algorithm to simulate brittle fractures in porous materials. Peridynamics theory is promising for fracture simulation since it allows discontinuities in the displacement field. However, they remain computationally expensive. Besides, there exists the surface effect in peridynamics due to the incomplete neighborhoods near the boundaries, including the outer boundaries and the boundaries of inner pores in porous materials. The proposed approach couples continuum mechanics and peridynamics into a closed equation system and an adaptive algorithm is developed to activate the peridynamics according to a strength criterion. In addition, the surface effect is corrected by introducing an improved peridynamic model with dual and variable horizons. We conduct the simulations using the relevant discretization scheme in each domain, i.e., the discontinuous Galerkin finite-element in the peridynamic domain and the continuous finite-element in the continuum mechanics domain. Two-dimensional numerical examples illustrate that successful fracture simulations of random porous materials can be achieved by this approach. In addition, the impact of distribution and size of pores on the fractures of porous materials is also investigated.

Time complexity analysis of the progressive fracture simulation

This section mainly introduces the detail analysis of the time complexity for the progressive fracture simulation algorithm. The main time costs of the algorithm includes assembling the element stiffness matrix and global stiffness matrix, solving linear equations ${\mathbf {K}}{\varvec{u}}= {\mathbf {F}}$, and processing broken bonds and reducing the stiffness matrix.

For the first part of the simulation, the element stiffness matrix is constructed by computing the bonds. So that if the assembling procedure need to process M bonds, the it will be a O(M) time complexity procedure. Assuming the assembling procedure will be performed S turns, and there shall be $M_i$ bonds that need to be process in the i-th turn. That is $M = \sum _i^S M_i$, and the total time complexity of this part will be O(M).

The second and third parts of the algorithm is hard to optimize but easy to estimate the time complexity. We use the LU decomposition method for sparse matrix in the SciPy package, which should have a ceiling of time complexity of $O(N^3)$ if and only if it is dealing a dense matrix, where N is the number of degree of the freedom (the dof). In fact, the sparsity of the matrix in our circumstance is very low, which means the time cost of our problem can be taken as a $O(N^{3/2})$ time complexity of the degree of the freedom (the dof) [63]. Say we process K times of solving linear equations, and the dof of our linear system is $N_{\text {P}}$, the total time complexity of this part will be $O(K) O(N_{\text {P}}^{3/2})$. In our numerical examples, one simulation need about 200-300 times to solving linear system, while the total steps is set to 100.

Since a rough estimate of time complexity of the whole simulation is obtained as $O(M) + O(K) O(N_{\text {P}}^{3/2})$, we need to specify the value of $M_i$ and S in running. To simplify the time complexity analysis, we take a single crack path simulation as an instance. In the PD model, damaged bonds in the microscale reveal the crack path in the macroscale. During the simulation, the propagation of the crack is progressive. To maintain the correctness of the simulation, the PD zone must cover all the damage zone. Therefore, the progressive algorithm will expand the PD zone step by step.

For example, we consider a rectangle plate with a final penetrating crack inside (see Fig. 4). The penetrating crack path has the length of L, and the width of $W_0$. The path will be covered by a set of horizon for applying the progressive fracture simulation algorithm. The minimum width of the damage zone is $W = 4\delta$, and the maximum length of each step during simulation is $\delta$. Say the crack will stop at step $S = L / \delta$, then for the specific step i, the current length of crack $L_i = i\delta$.

Let the grid size be $\Delta x$ and the number of finite elements is N, that is $N\Delta x^2 = L W_0$. In addition, the horizon radius of the PD model will be $\delta = m \Delta x$. Each horizon has $n_{\text {bond}} = \pi \delta ^2 / \Delta x^2$ elements. While for step i it will has $n_{\text {FE}} = L_i W / \Delta x^2 = 4i \delta ^2 / \Delta x^2 = 4 i m^2$. Therefore, the total number of bonds $M_i$ of step i will be

$$\begin{aligned} M_i = n_{\text {FE}} \cdot n_{\text {bond}} = 4 \pi m^4 \cdot i. \end{aligned}$$

(A1)

In addition, the total number of bonds M of all S step will be

$$\begin{aligned} M(S) = \sum _{i = 1}^S M_i = 4 \pi m^4 \Bigl( \sum _{i = 1}^S i \Bigr) = 2 \pi m^4 S(S+1). \end{aligned}$$

(A2)

In practice, the number of finite elements $N_{\text {E}}$ and the horizon size parameter m are determined manually. In our numerical examples, the value of m is set to 3, and $m \ll N$. Finally, we substitute the $S = L / \delta$ into M,

$$\begin{aligned} \begin{aligned} M&= 2 \pi m^4 \Bigl( \frac{L}{\delta } \Bigr) \Bigl( \frac{L}{\delta }+1 \Bigr) \le 2 \pi m^4 \Bigl( \frac{L + \delta }{\delta } \Bigr) ^2 \\&= 2 \pi m^4 \Bigl( \frac{L^2 + 2L\delta + \delta ^2}{m^2 \Delta x^2} \Bigr) \\&= 2 \pi m^2 \Bigl( \frac{L^2}{\Delta x^2} + 2m\frac{L}{\Delta x} + m^2 \Bigr) \\&= 2 \pi m^2 N_{\text {E}}+ 4 \pi m^3 \sqrt{N_{\text {E}}} + 2 \pi m^4 \\ \implies O(M)&= O\left( m^2\right) O(N_{\text {E}}). \end{aligned} \end{aligned}$$

(A3)

Thus, the algorithm has a linear time complexity in assembling stiffness.

In summary, the time complexity $T(m, K, N_{\text {P}}, N_{\text {E}})$ of the progressive fracture simulation algorithm is

$$\begin{aligned} T(m, K, N_{\text {P}}, N_{\text {E}}) = O(m^2) O(N_{\text {E}}) + O(K) O\left( N_{\text {P}}^{3/2}\right) , \end{aligned}$$

(A4)

where m is horizon size parameter in the PD model, K is the total number for solving linear system, $N_{\text {P}}$ and $N_{\text {E}}$ are the number of nodes and elements, respectively.