Multiscale analysis-based peridynamic simulation of fracture in porous media

Abstract

The simulation of fracture in large-scale structures made of porous media remains a challenging task. Current techniques either assume a homogeneous model, disregarding the microstructure characteristics, or adopt a micro-mechanical model, which incurs an intractable computational cost due to its complex stochastic geometry and physical properties, as well as its nonlinear and multiscale features. In this study, we propose a multiscale analysis-based dual-variable-horizon peridynamics (PD) model to efficiently simulate macroscopic structural fracture. The influence of microstructures in porous media on macroscopic structural failure is represented by two PD parameters: the equivalent critical stretch and micro-modulus. The equivalent critical stretch is calculated using the microscale PD model, while the equivalent micro-modulus is obtained through the homogenization method and energy density equivalence between classical continuum mechanics and PD models. Numerical examples of porous media with various microstructures demonstrate the validity, accuracy, and efficiency of the proposed method.

Graphical abstract

Keywords

porous media / multiscale / variable-horizon peridynamic / equivalent critical stretch / equivalent micro-modulus

1 Introduction

The porous media composes of solid structures with the presence of pores [1–3]. These voids can be micro-scale holes, air pockets, or interconnected pore networks, giving the material a low density and highly porous characteristics. Porous media possess a range of unique properties and applications, including lightweight, sound absorption and thermal insulation. Their applications span across areas such as building materials, biomedical materials, and energy storage. Understanding the fracture behavior of porous medias is crucial to ensure their safety in practical engineering and applications. Furthermore, studying the fracture behavior of porous materials allows for the prediction of their performance, such as in terms of toughness and crack resistance, in real-world applications. This helps guide material design and optimization, improving the performance and reliability of the materials. The numerical simulation of fracture in porous materials holds significant importance in engineering and scientific fields. It provides a deep understanding of the fracture behavior of porous materials and offers valuable information for material design and engineering applications. The complexity and multiscale nature of porous materials make fracture simulation computationally demanding. In particular, when simulating large-scale porous material structures such as practical engineering components or material specimens, the computational cost further increases. Therefore, numerical simulation of fracture in porous media remains a challenge.

Numerical simulation techniques have been developed to analyze fracture and failure in porous media, utilizing two main types of models: local models based on partial differential equations and nonlocal models based on integral equations, such as peridynamics (PD). Local models include methods like the finite element method [4–7], extended finite element method [8–11], phase-field model [12–20] and cohesive zone model [21,22]. For example, Kang and Bolton [4] developed a two-dimensional elastic-absorption finite element model for isotropic elastic porous noise control materials. He [11] presented a hydromechanical model for fluid flow in fractured porous media. He et al. [12] conducted a numerical investigation of the effect of microporosity and micropore shapes on macrofracture behavior in porous amorphous silica. Zhou et al. [17] proposed a phase field model for fluid-driven dynamic crack propagation in poroelastic media. Cajuhi et al. [19] proposed a mechanical and computational model to describe the coupled problem of poromechanics and cracking in variably saturated porous media. Additionally, a finite element model was developed in [21] to analyze hydraulic fracturing in partially saturated porous media, where fracture propagation was governed by a cohesive fracture model.

PD theory, introduced by Silling [23], is a nonlocal model that employs spatial integration instead of differentiation to formulate the equation of motion. This approach resolves the singularity problem caused by spatial discontinuity in fracture problems. PD models have been extensively used to study crack nucleation and propagation in porous media. For instance, Sun and Fish [24] proposed a framework that coupled non-ordinary state-based PD with the finite element method to simulate fracture propagation in saturated porous media. Sun et al. [25] presented a fully coupled hydraulic fracture propagation simulation method using a hybrid finite element and PD approach for porous media. Shen et al. [26] developed a PD model with an improved energy-based surface correction to simulate fractures in porous materials. Other studies, such as those by Ni et al. [27–28], Mehrmashhadi et al. [29], Wu et al. [30], Wu et al. [31–33], and Chen et al. [34] introduced PD models for simulating brittle damage and fracture in elastic porous materials or concretes based on an intermediate homogenization approach. These models and methods contribute to our understanding of fracture behavior in porous media. However, the fracture simulation of large-scale structures composed of porous media poses a significant challenge. Unfortunately, an effective PD model that considers the multiscale structure of porous materials and efficiently simulates their fracture behavior is still lacking.

The objective of this paper is to introduce a multiscale model capable of directly capturing the influence of porous materials’ microstructure on macroscopic fracture evolution. To achieve this, we propose a dual-variable-horizon PD model based on multiscale analysis. Our approach focuses on extracting the influence of microstructures in porous media on macroscopic structural failure in the form of two equivalent parameters: the equivalent critical stretch and micro-modulus. By sequentially coupling the microscale PD and macroscale PD using these equivalent parameters, we establish a framework that enables efficient simulation of macroscopic structure fracture. Additionally, the dual-variable-horizon PD model we employ has the potential to account for surface effects in porous materials [35–40], enhancing the accuracy of our simulations. With the proposed multiscale analysis-based PD framework, we aim to simulate fracture in porous media structures with improved accuracy and reduced computational costs, particularly for large-scale structures.

The remainder of this paper is organized as follows. In Section 2, we provide a detailed description of the proposed multiscale analysis-based PD method. Section 3 presents numerical examples that assess the feasibility and effectiveness of our approach. Finally, in Section 4, we provide concluding remarks summarizing the contributions of this study and discuss potential avenues for future research.

2 Multiscale analysis-based peridynamic model

In our study, we focus on a porous media structure that can be viewed as a collection of representative volume elements (RVEs) at the microscale. These RVEs have a uniform size, denoted by \epsilon. Additionally, all the RVEs share the same probability distribution of pores, ensuring consistency across the structure. We denote the overall structure under investigation as Ω which is composed of multiple RVEs, each adhering to the same probability distribution model [41]. Importantly, the size of the RVEs is considerably larger than the characteristic size of the pores. However, the RVE size remains significantly smaller than the characteristic size of the macroscale structure (Fig.1).

At the microscale, we construct a dual-variable-horizon PD model to simulate the fracture behavior of the individual RVEs. Within this model, we define the equivalent critical stretch, which represents the strain level at which fracture occurs in the RVE. Furthermore, we employ the homogenization method to determine the equivalent elastic tensor defined on the RVE. This tensor characterizes the overall mechanical response of the RVE, taking into account its microstructure. Additionally, we define the equivalent micro-modulus based on the principle of energy density equivalence between classical continuum mechanics (CCM) and PD models. Moving to the macroscale, we construct a dual-variable-horizon PD model that incorporates the two equivalent parameters obtained from the microscale analysis. These parameters capture the essential microstructural features and their influence on the macroscopic fracture behavior of the porous media structure.

2.1 Microscale peridynamic model for equivalent critical stretch

To eliminate the influence of surface effects in internal pores of the RVE, a dual-variable-horizon PD model is introduced [35,36], and it is further applied to calculate the equivalent critical stretch. For convenience, the local coordinate system y=x/\epsilon \in Y is introduced on the unit RVE Y (Fig.1). For a material point y, its neighborhood is , where the neighborhood radius \delta \left(y\right) can be expressed as

(1)\delta \left(y\right)=min\left\{{\delta }_0,\underset{p\in \mathrm{\partial }Y}{min}\Vert y-p\Vert \right\},

where \mathrm{\partial }Y is the inner and outer boundaries of the RVE, {\delta }_0 is the pre-set size of the interaction radius. Similarly, we denote as the dual neighborhood of y. Considering a constitutive equation based on small deformation hypothesis, the force \hat{f}(y^{\prime },y) of material point y on y^{\prime }, and the reaction force \hat{f}(y^{\prime },y) of material point y^{\prime } on y satisfy the following equations

(2)\hat{f}\left(y^{\prime },y\right)=\frac{1}{2}c(y,\xi )\left(e_{\xi }\otimes e_{\xi }\right)\eta \left(y^{\prime },y\right),y^{\prime }\in H_y,

(3)-\hat{f}\left(y^{\prime },y\right)=\{\begin{array}{l}-{\displaystyle \frac{1}{2}}c\left(y^{\prime }-\xi \right)\left(e_{\xi }\otimes e_{\xi }\right)\left(-\eta \left(y^{\prime },y\right)\right),y^{\prime }\in H_{y^{\prime }},\\ 0,\end{array}

where c(y,\xi ) is a scalar function representing the constitutive parameter on the RVE, \xi =y^{\prime }-y represents the bond vector between material points y and y^{\prime } inside the RVE, e_{\xi }=\xi /\parallel \xi \parallel denotes the unit directional vector of \xi, \eta \left(y^{\prime },y\right)=u\left(y^{\prime }\right)-u\left(y\right) represents the displacement of the bond vector, and u\left(y\right) denotes the displacement of material point y.

The static equilibrium equation for the dual-variable peridynamics (DVH-PD) model on the RVE is given as follows

(4){{\displaystyle \int }}_{H_{\text{y}}}\hat{f}\left(y^{\prime },y\right)\text{d}{V}_{y^{\prime }}-{{\displaystyle \int }}_{H_{y^{\prime }}}\hat{f}\left(y^{\prime },y\right)\text{d}{V}_{y^{\prime }}+b(y)=0,y,y^{\prime }\in Y

The constitutive scalar function c(y,\xi ) of the DVH-PD model is obtained through energy equivalence of CCM and PD models [42] as follows

(5)\text{A}(y)=\frac{1}{4}{{\displaystyle \int }}_{H_{\text{y}}}c(y,\xi )\parallel \xi {\parallel }^{-2}\xi \otimes \xi \otimes \xi \otimes \xi \text{d}{V}_{\text{y}},

where A\left(y\right) is the fourth-order elastic tensor of the CCM model in the RVE. The bond stretch, s, is defined by

(6)s=\frac{\Vert \xi +\eta \left(y^{\prime },y\right)\Vert -\Vert \xi \Vert }{\Vert \xi \Vert }.

To incorporate the failure behavior of the porous media, we implement a failure law by defining different history-dependent scalar-valued functions \mu (\xi ,t) for the bonds as follows

(7)

where t and t^{\prime } represent the loading steps during the calculation process, and s_0 denotes the critical stretch values. Once a bond fractures, \mu =0 and it is not irreversible.

According to the multiscale representation of porous media, as depicted in Fig.1, an RVE corresponds to a material point labeled x within the macroscopic structure \mathrm{\Omega }. Consequently, it becomes possible to define the critical stretch of the macroscopic structure by simulating fracture within the RVE. To illustrate this, let’s consider the example of uniaxial tension in the y_1 direction of an RVE, with a given random sample denoted as {\omega }^m. The uniaxial tension boundary conditions are applied to the left and right sides of the RVE, as shown in Fig.2(a). By employing the DVH-PD model, as the loading increases to {\tilde{u}}_1\left({\omega }^m\right), a penetrating crack perpendicular to the loading direction emerges within the RVE, as demonstrated in Fig.2(b). And the effective critical stretch {\hat{s}}_{01}\left({\omega }^m\right) corresponding to random sample {\omega }^m along the x_1 direction of macroscopic structure can be defined as

(8){\hat{s}}_{01}\left({\omega }^m\right)=\frac{\epsilon +{\tilde{u}}_1\left({\omega }^m\right)-\epsilon }{\epsilon }=\frac{{\tilde{u}}_1\left({\omega }^m\right)}{\epsilon }.

Similarly, applying uniaxial tensile boundary conditions in the y_2 and y_3 directions, the effective critical stretch {\hat{s}}_{02}\left({\omega }^m\right) and {\hat{s}}_{03}\left({\omega }^m\right) can also be obtained as Eq. (8).

Since the pores are randomly distributed within an RVE and their centroid positions follow a certain probability distribution function, we need to calculate the statistical average of effective critical stretches by taking different random samples. Based on Kolmogorov’s strong law of large numbers in probability theory, the statistical critical stretch can be obtained as follows

(9){\tilde{s}}_{0i}=\underset{M\to \mathrm{\infty }}{\mathrm{l}\mathrm{i}\mathrm{m}}\frac{{\displaystyle \sum _{m=1}^M{\hat{s}}_{0i}\left({\omega }^m\right)}}{M},

where M represents the number of random samples. Further, the critical stretch in any arbitrary direction of macroscopic structure, i.e., the equivalent critical stretch, can be defined through the ellipsoidal averaging [43] of the statistical critical stretch along the three orthogonal directions (Fig.3) as follows

(10){\overline{s}}_0\left(\xi \right)=\frac{1}{\sqrt{{\left({\displaystyle \frac{n_1}{{\overline{s}}_{01}}}\right)}^2+{\left({\displaystyle \frac{n_2}{{\overline{s}}_{02}}}\right)}^2+{\left(d{\displaystyle \frac{n_3}{{\overline{s}}_{03}}}\right)}^2}},

where (n_1,n_2,n_3{)}^{\mathrm{T}}=\xi /\Vert \xi \Vert.

2.2 Asymptotic homogenization model for equivalent micro-modulus

Based on the multiscale characterization of porous structure Ω, the continuum mechanics model (CMM) model can be described as

(11)-\frac{\mathrm{\partial }}{\mathrm{\partial }{x}_j}\left[a_{ijkl}^{\epsilon }\left(x\right)\frac{1}{2}\left(\frac{\mathrm{\partial }{u}_k^{\epsilon }}{\mathrm{\partial }{x}_l}+\frac{\mathrm{\partial }{u}_l^{\epsilon }}{\mathrm{\partial }{x}_k}\right)\right]=q_i\left(x\right),x\in \mathrm{\Omega },

where a_{ijkl}^{\epsilon }\left(x\right) denotes the components of the elasticity tensor, u^{\epsilon } is the displacement field, and q_i\left(x\right) denotes the external body force. Additionally, as depicted in Fig.1, x and y correspond to the global coordinates of the macroscopic structure \mathrm{\Omega }and the local coordinates of the microscopic RVE Y, respectively, and they are mathematically related as follows

(12)y=x/\epsilon \Rightarrow \frac{\mathrm{\partial }}{\mathrm{\partial }{x}_i}=\frac{\mathrm{\partial }}{\mathrm{\partial }{x}_i}+\frac{1}{\epsilon }\frac{\mathrm{\partial }}{\mathrm{\partial }{y}_i},i=1,2,3.

Then, we have a_{ijkl}^{\epsilon }\left(x\right)=a_{ijkl}\left(y\right).

According to the asymptotic homogenization method [41,44], the displacement field u^{\epsilon }\left(x\right) can be expanded as follows

(13)u^{\epsilon }\left(x\right)=u^0(x,y)+\epsilon u^1(x,y)+{\epsilon }^2{u}^2(x,y)+\cdots .

By substituting Eqs. (12) and (13) into Eq. (11) and equating the coefficients of {\epsilon }^i\left(i=1,2,3,\dots \right) on both sides, we obtain a series of equations. Specifically, the terms containing {\epsilon }^{-2} on both sides are equated to each other to obtain u^0\left(x,y\right)=u^0\left(x\right). And for all {\epsilon }^{-1} terms, we obtain

(14)u^1\left(x,y\right)=N_{\alpha }\left(y\right)\frac{\mathrm{\partial }{u}^0}{\mathrm{\partial }{x}_{\alpha }}\left(x\right)+{\tilde{u}}^1\left(x\right),

where {\tilde{u}}^1\left(x\right) represents a function that is independent of the microscopic coordinates y. In three-dimensional cases, N_{\alpha }\left(y\right) is a third-order matrix-valued cell function. For any given random sample {\omega }^m, the cell function N_{\alpha }(y,{\omega }^m) can be defined as follows:

(15)-\frac{\mathrm{\partial }}{\mathrm{\partial }{y}_j}\left[a_{ijkl}\left(y,{\omega }^m\right)\frac{1}{2}\left(\frac{\mathrm{\partial }{N}_{\alpha km}}{\mathrm{\partial }{y}_l}+\frac{\mathrm{\partial }{N}_{\alpha lm}}{\mathrm{\partial }{y}_k}\right)\right]=\frac{\mathrm{\partial }{a}_{ijm\alpha }\left(y,{\omega }^m\right)}{\mathrm{\partial }{y}_j},y\in Y^m

(16)N_{\alpha m}\left(y,{\omega }^m\right)=0,y\in \mathrm{\partial }{Y}^m.

After calculating N_{\alpha }(y,{\omega }^m), the homogenized elastic tensor for sample {\omega }^m can be defined as

(17)\begin{array}{rl}{\hat{a}}_{ijkl}\left({\omega }^m\right)=& \frac{1}{\left|Y^m\right|}{{\displaystyle \int }}_{Y^m}(a_{ijkl}\left(y,{\omega }^m\right)+a_{ijpq}\left(y,{\omega }^m\right)\\ & \cdot \frac{1}{2}(\frac{\mathrm{\partial }{N}_{kpl}\left(y,{\omega }^m\right)}{\mathrm{\partial }{y}_q}+\frac{\mathrm{\partial }{N}_{kql}\left(y,{\omega }^m\right)}{\mathrm{\partial }{y}_p}))\text{d}y.\end{array}

The homogenized micro-modulus {\hat{c}}_0(x,\xi ,{\omega }^m) can be obtained by equating the elastic energy density between the CCM and DVH-PD models

(18){\hat{a}}_{ijkl}\left({\omega }^m\right)=\frac{1}{4}{{\displaystyle \int }}_{H_x}{\hat{c}}_0\left(x,\xi ,{\omega }^m\right)\parallel \xi {\parallel }^{-2}{\xi }_i\otimes {\xi }_j\otimes {\xi }_k\otimes {\xi }_l\text{d}{V}_{x^{\mathrm{\prime }}}.

Based on Kolmogorov’s strong law of large numbers, the equivalent micro-modulus can be obtained by taking M random samples

(19){\overline{c}}_0\left(x,\xi \right)=\underset{M\to \mathrm{\infty }}{\mathrm{l}\mathrm{i}\mathrm{m}}\frac{{\displaystyle \sum _{m=1}^M{\hat{c}}_0\left(x,\xi ,{\omega }^m\right)}}{M}.

2.3 Macroscopic peridynamic model

The above analysis and derivation lead to the extraction of two equivalent parameters that capture the influence of microstructures: the equivalent critical stretch and micro-modulus. Then, we can construct the PD model on the macroscale homogenized structure as follows

(20){{\displaystyle \int }}_{H_x}\hat{f}\left(x^{\prime },x\right)\text{d}{V}_{x^{\prime }}-{{\displaystyle \int }}_{H_{x^{\prime }}}\hat{f}\left(x,x^{\prime }\right)\text{d}{V}_{x^{\prime }}+b(x)=0,x,x^{\prime }\in \overline{\text{Ω}},

where H_x and {H^{\prime }}_x denotes the neighborhood and dual neighborhood of x defined as

(21)

The size of the neighborhood and dual neighborhood is controlled by \delta \left(x\right), which is defined as follows:

(22)\delta \left(x\right)=min\left\{{\delta }_0,\underset{p\in \mathrm{\partial }\overline{\Omega }}{min}\parallel x-p\parallel \right\},

where \mathrm{\partial }\overline{\mathrm{\Omega }} represents the boundary of the homogenized structure, and {\delta }_0 is the default radius of the neighborhood size. The bond force functions are defined by

(23)\hat{f}\left(x^{\prime },x\right)=\frac{1}{2}{\overline{c}}_0(x,\xi )\left(e_{\xi }\otimes e_{\xi }\right)\eta \left(x^{\prime },x\right),x^{\prime }\in H_x,

(24)-\hat{f}\left(x,x^{\prime }\right)=\{\begin{array}{l}-{\displaystyle \frac{1}{2}}{\overline{c}}_0\left(x^{\prime },-\xi \right)\left(e_{\xi }\otimes e_{\xi }\right)\left(-\eta \left(x^{\prime },x\right)\right),x^{\prime }\in H_{x^{\prime }},\\ 0,\end{array}

where \eta \left(x^{\prime },x\right)=u\left(x^{\prime }\right)-u\left(x\right) represents the displacement of the bond vector. For the macroscale homogenized structure, the stretch and fracture criterion of the bond are defined as follows

(25)s=\frac{\Vert \xi +\eta \Vert -\Vert \xi \Vert }{\Vert \xi \Vert },

(26)

where t and t^{\prime } represent the loading steps, and {\overline{s}}_0\left(\xi \right) is determined by Eq. (10).

3 Numerical examples

This section presents examples to verify the effectiveness and efficiency of the proposed method and analyzes the influence of microstructures with different random distributions of pores on the macroscopic fracture of porous media. All examples are implemented using a self-developed finite element numerical software package that incorporates both continuous elements (CE) and discrete elements (DE) [45,46]. DEs are applied in the crack area, while CEs are used in the remaining structure. To counteract the mesh effects, the neighborhood radius {\delta }_0 was set to three times the average mesh size. The micro-modulus was chosen as c(x,\xi )=c_0\Vert \xi {\Vert }^2\exp (-\Vert \xi \Vert /L_c), with the characteristic length L_{\mathrm{c}} set to one-third of the average mesh size. All numerical experiments were conducted on a workstation equipped with the AMD EPYC 7453 32-core CPU and 256 GB of RAM.

3.1 Verification of proposed method

This section primarily focuses on validating the effectiveness and efficiency of the proposed method. Three kinds of periodic porous material consisting of 5\times 5 RVEs with the pore radii of 2.0, 2.5, and 3.0 mm are considered as shown in Fig.4. The matrix material is brittle with Young’s modulus of E=3\mathrm{G}\mathrm{P}\mathrm{a}, Poisson’s ratio of v=1/3, and tensile failure strength of {\sigma }_{\mathrm{f}}= 26 MPa. For the RVE, the equivalent critical stretch is calculated by considering the uniaxial tensile boundary condition of 0.15 mm along the y_1 direction, and the fracture process is numerically simulated using 200 incremental steps. For the macroscopic structure, the uniaxial tensile boundary condition of 0.6 mm is applied along x_1 direction. The numerical results are obtained using both the multiscale method and single-scale direct PD fracture analyses. Both simulations are carried out with 100 incremental steps.

Fig.5 shows the local damage and fracture of three RVEs with different pore sizes, and it demonstrates that pore sizes significantly affect the fracture of microscopic RVE. Fig.6 shows the curves of support reaction force versus imposed displacement, which are obtained by the single-scale direct PD (DPD) model and by the present model. Tab.1 presents the number of elements, nodes, and freedom degrees in the finite element method, along with the running times of the DPD simulation and the multiscale simulation. The multiscale simulation involves the PD simulation of the RVE at the microscale (Micro-PD) and the PD simulation of the homogenized structures at the macroscale (Macro-PD). From Fig.6 and Tab.1, it is evident that the results of both models exhibit good agreement, and the present method demonstrates significantly higher computational efficiency compared to that of the single-scale direct PD. This indicates that the proposed method not only obtains accurate results of fracture failure of macroscopic porous media structures, but also greatly improves the computational efficiency.

3.2 The effect of pore distribution

In this section, we conduct an extensive analysis of how microstructures within porous media influence macroscopic fracture behavior. Two RVEs with a side length of 1mm and with different random distributions of pores, as shown in Fig.7, are considered, including normal distribution and uniform distribution. Notably, each RVE contains 25 pores, each with a radius of 30 μm. The matrix material of the RVEs is characterized by a Young’s modulus of 192 GPa, a Poisson’s ratio of 0.3333, and a yield strength of {\sigma }_{\mathrm{f}}= 520 MPa. A total of 100 sets of random samples for each kind of RVE are generated and the fracture is simulated through 100 incremental steps. Fig.8 provides insight into the expected values of micro-modulus and critical stretch, considering different numbers of random samples. It is essential to note that different samples yield distinct outcomes. However, as the number of samples increases, the expected values gradually converge. Evidently, as illustrated in Fig.8, the scatter of data diminishes with an increasing number of samples. Furthermore, we explore the influence of microstructures on macroscopic fracture by examining two distinct porous media structures with varying boundary conditions. The assessment includes the three-point bending test and the tensile-shear test of double-notched specimens. These investigations are pivotal for comprehending the interplay between microstructural features and macroscopic fracture behavior.

3.2.1 Three-point bending test

The three-point bending test is performed on specimens using two types of RVEs as depicted in Fig.9. A displacement boundary condition of 0.25 mm is applied vertically at the midpoint of the specimen's upper surface. The specimen has two support bases at the bottom, with the left support being simply supported and the right support fixed. Additionally, a crack with a length of 23.30 mm is introduced. The entire simulation process consists of 100 incremental loading steps.

Fig.10 shows the fracture results of the specimens with different distribution of RVEs. The specimen with a normal distribution of fractures at the 27th step stops fracturing after the 67th step. The specimen with a random uniform distribution of fractures at the 33rd step completely cracks open after the 88th step. This is because different distributions result in different equivalent critical stretch, and the pores with a normal distribution are more prone to fracture. It is also worth noting that the equivalent critical stretches of these two specimens are very small, and the fracture patterns are similar to an opening-type crack (I-type crack). Specifically, because the pre-existing crack and the loading direction are not on the same straight line, the crack tip deviates to a certain degree. Then, as the stretch of the bonds near the deviated crack tip exceed {\overline{s}}_0, the crack rapidly develops along a straight line, almost parallel to the vertical direction. In the final stage, as the crack gradually propagates through the specimen, the shear effect becomes more pronounced, causing the crack to deviate toward the centerline. Tab.2 presents the number of elements, nodes, freedom degrees in the finite element method, and the respective running times for the single-scale direct PD simulation and the multiscale simulation. The computational efficiency of the present method is significantly higher than that of the single-scale direct PD. The latter results in an excessive number of elements, leading to an unmanageable computational cost.

3.2.2 Tensile-shear problem of biaxial notched specimens

Biaxial notched specimens are examined using two types of RVEs as depicted in Fig.11. Uniaxial tensile boundary conditions of 0.05 mm are applied along the vertical direction, while shear boundary conditions of 0.05 mm are applied along the horizontal direction. The entire simulation process consists of 100 incremental loading steps.

Fig.12 and Fig.13 show the fracture process of specimens with two kinds of RVEs. It can be observed that fracture is more likely to occur in the specimen with normal distribution of pores, and the cracks for the case with normal distribution are significantly longer than that with uniform distribution under the same loading. Moreover, the fracture mode of this example belongs to a mixed mode of opening and sliding (I + II type crack). It can be found from the fracture simulation process that the two cracks develop independently, almost parallel to each other, and remain basically symmetric. During the final stage of crack development, the crack direction deviates inward. Tab.3 presents the number of elements, nodes, and freedom degrees in the finite element method, along with the respective running times for both the single-scale direct PD simulation and the multiscale simulation. It is evident that using the single-scale PD model directly to simulate the fracture of macroscopic porous structures results in an excessively large number of degrees of freedom in the finite element equation, rendering it unsolvable. In contrast, the proposed method takes into account the influence of microstructures, allowing for efficient and accurate simulation of the macroscopic fracture behavior of porous media structures.

4 Conclusions

This paper proposes a sequential multiscale framework that couples the microscale PD and the macroscale PD through two equivalent parameters. The fracture of porous media is simulated using the macroscale dual-variable-horizon PD model with equivalent critical stretch and micro-modulus, which account for the influence of microstructure. Notably, the dual-variable-horizon PD model effectively eliminates surface effects on both the microscale and macroscale structures. Moreover, the homogenization of porous media significantly enhances the efficiency of macroscopic PD fracture simulation. The comparison between results from the multiscale method and direct fracture analysis using the PD model at a single scale verifies the validity and efficiency of this approach. The proposed multiscale analysis-based peridynamic framework has the potential to accurately simulate fracture in porous structures with reduced computational cost, particularly for large-scale structures.

References

[1] AdlerP. Porous Media: Geometry and Transports. Stoneham, MA: Elsevier, 2013

[2] EhlersWBluhmJ. Porous Media: Theory, Experiments and Numerical Applications. Berlin: Springer Science & Business Media, 2002

[3] LiuP SChenG F. Porous Materials: Processing and Applications. Waltham, MA: Elsevier, 2014

[4] Kang Y J, Bolton J S. Finite element modeling of isotropic elastic porous materials coupled with acoustical finite elements. Journal of the Acoustical Society of America, 1995, 98(1): 635–643

[5] Han N, Guo R. Two new Voronoi cell finite element models for fracture simulation in porous material under inner pressure. Engineering Fracture Mechanics, 2019, 211: 478–494

[6] Zhang R, Guo R. Voronoi cell finite element model to simulate crack propagation in porous materials. Theoretical and Applied Fracture Mechanics, 2021, 115: 103045

[7] Nilsen H M, Larsen I, Raynaud X. Combining the modified discrete element method with the virtual element method for fracturing of porous media. Computational Geosciences, 2017, 21(5–6): 1059–1073

[8] De BorstR. Computational Methods for Fracture in Porous Media: Isogeometric and Extended Finite Element Methods. Cambridge, MA: Elsevier, 2017

[9] Mohtarami E, Baghbanan A, Hashemolhosseini H, Bordas S P. Fracture mechanism simulation of inhomogeneous anisotropic rocks by extended finite element method. Theoretical and Applied Fracture Mechanics, 2019, 104: 102359

[10] Rezanezhad M, Lajevardi S A, Karimpouli S. An investigation on prevalent strategies for XFEM-based numerical modeling of crack growth in porous media. Frontiers of Structural and Civil Engineering, 2021, 15(4): 914–936

[11] He B. Hydromechanical model for hydraulic fractures using XFEM. Frontiers of Structural and Civil Engineering, 2019, 13(1): 240–249

[12] He B, Vo T, Newell P. Investigation of fracture in porous materials: A phase-field fracture study informed by ReaxFF. Engineering with Computers, 2022, 38(6): 5617–5633

[13] Lee S, Wheeler M F, Wick T. Pressure and fluid-driven fracture propagation in porous media using an adaptive finite element phase field model. Computer Methods in Applied Mechanics and Engineering, 2016, 305: 111–132

[14] He B, Schuler L, Newell P. A numerical-homogenization based phase-field fracture modeling of linear elastic heterogeneous porous media. Computational Materials Science, 2020, 176: 109519

[15] Heider Y, Reiche S, Siebert P, Markert B. Modeling of hydraulic fracturing using a porous-media phase-field approach with reference to experimental data. Engineering Fracture Mechanics, 2018, 202: 116–134

[16] Zeng Q, Liu W, Yao J, Liu J. A phase field based discrete fracture model (PFDFM) for fluid flow in fractured porous media. Journal of Petroleum Science Engineering, 2020, 191: 107191

[17] Zhou S, Zhuang X, Rabczuk T. Phase-field modeling of fluid-driven dynamic cracking in porous media. Computer Methods in Applied Mechanics and Engineering, 2019, 350: 169–198

[18] Heider Y, Markert B. A phase-field modeling approach of hydraulic fracture in saturated porous media. Mechanics Research Communications, 2017, 80: 38–46

[19] Cajuhi T, Sanavia L, de Lorenzis L. Phase-field modeling of fracture in variably saturated porous media. Computational Mechanics, 2018, 61(3): 299–318

[20] Zhou S, Zhuang X, Rabczuk T. Phase field method for quasi-static hydro-fracture in porous media under stress boundary condition considering the effect of initial stress field. Theoretical and Applied Fracture Mechanics, 2020, 107: 102523

[21] Dastjerdy F, Barani O, Kalantary F. Modeling of hydraulic fracture problem in partially saturated porous media using cohesive zone model. International Journal of Civil Engineering, 2015, 10: 86

[22] Komijani M, Gracie R, Yuan Y. Simulation of fracture propagation induced acoustic emission in porous media. Engineering Fracture Mechanics, 2020, 229: 106950

[23] Silling S A. Reformulation of elasticity theory for discontinuities and long-range forces. Journal of the Mechanics and Physics of Solids, 2000, 48(1): 175–209

[24] Sun W, Fish J. Coupling of non-ordinary state-based peridynamics and finite element method for fracture propagation in saturated porous media. International Journal for Numerical and Analytical Methods in Geomechanics, 2021, 45(9): 1260–1281

[25] Sun Y, Chen B, Edwards M G, Li C. Investigation of hydraulic fracture branching in porous media with a hybrid finite element and peridynamic approach. Theoretical and Applied Fracture Mechanics, 2021, 116: 103133

[26] Shen S, Yang Z, Han F, Cui J, Zhang J. Peridynamic modeling with energy-based surface correction for fracture simulation of random porous materials. Theoretical and Applied Fracture Mechanics, 2021, 114: 102987

[27] Ni T, Pesavento F, Zaccariotto M, Galvanetto U, Zhu Q Z, Schrefler B A. Hybrid FEM and peridynamic simulation of hydraulic fracture propagation in saturated porous media. Computer Methods in Applied Mechanics and Engineering, 2020, 366: 113101

[28] Ni T, Pesavento F, Zaccariotto M, Galvanetto U, Schrefler B A. Numerical simulation of forerunning fracture in saturated porous solids with hybrid FEM/Peridynamic model. Computers and Geotechnics, 2021, 133: 104024

[29] Mehrmashhadi J, Chen Z, Zhao J, Bobaru F. A stochastically homogenized peridynamic model for intraply fracture in fiber-reinforced composites. Composites Science and Technology, 2019, 182: 107770

[30] Wu L, Huang D, Wang H, Ma Q, Cai X, Guo J. A comparison study on numerical analysis for concrete dynamic failure using intermediately homogenized peridynamic model and meso-scale peridynamic model. International Journal of Impact Engineering, 2023, 179: 104657

[31] Wu P, Zhao J, Chen Z, Bobaru F. Validation of a stochastically homogenized peridynamic model for quasi-static fracture in concrete. Engineering Fracture Mechanics, 2020, 237: 107293

[32] Wu P, Yang F, Chen Z, Bobaru F. Stochastically homogenized peridynamic model for dynamic fracture analysis of concrete. Engineering Fracture Mechanics, 2021, 253: 107863

[33] Wu P, Chen Z. Peridynamic electromechanical modeling of damaging and cracking in conductive composites: A stochastically homogenized approach. Composite Structures, 2023, 305: 116528

[34] Chen Z, Niazi S, Bobaru F. A peridynamic model for brittle damage and fracture in porous materials. International Journal of Rock Mechanics and Mining Sciences, 2019, 122: 104059

[35] Ren H, Zhuang X, Cai Y, Rabczuk T. Dual-horizon peridynamics. International Journal for Numerical Methods in Engineering, 2016, 108(12): 1451–1476

[36] Ren H, Zhuang X, Rabczuk T. Dual-horizon peridynamics: A stable solution to varying horizons. Computer Methods in Applied Mechanics and Engineering, 2017, 318: 762–782

[37] Rabczuk T, Ren H, Zhuang X. A nonlocal operator method for partial differential equations with application to electromagnetic waveguide problem. Computers, Materials & Continua, 2019, 59(1): 31–55

[38] Ren H, Zhuang X, Rabczuk T. A nonlocal operator method for solving partial differential equations. Computer Methods in Applied Mechanics and Engineering, 2020, 358: 112621

[39] Bie Y, Liu Z, Yang H, Cui X Y. Abaqus implementation of dual peridynamics for brittle fracture. Computer Methods in Applied Mechanics and Engineering, 2020, 372: 113398

[40] Zhang Y, Yang X, Wang X, Zhuang X. A micropolar peridynamic model with non-uniform horizon for static damage of solids considering different nonlocal enhancements. Theoretical and Applied Fracture Mechanics, 2021, 113: 102930

[41] Yang Z, Guan X, Cui J, Dong H, Wu Y, Zhang J. Stochastic multiscale heat transfer analysis of heterogeneous materials with multiple random configurations. Communications in Computational Physics, 2020, 27(2): 431–459

[42] ShenS KYangZ HCuiJ ZZhangJ Q. Dual-variable-horizon peridynamics and continuum mechanics coupling modeling and adaptive fracture simulation in porous materials. Engineering with Computers, 2022, 1–21

[43] Le Q V, Bobaru F. Surface corrections for peridynamic models in elasticity and fracture. Computational Mechanics, 2018, 61(4): 499–518

[44] Yang Z, Zhang Y, Dong H, Cui J, Guan X, Yang Z. High-order three-scale method for mechanical behavior analysis of composite structures with multiple periodic configurations. Composites Science and Technology, 2017, 152: 198–210

[45] Yang Z, Zheng S, Han F, Shen S, Guan X. An improved peridynamic model with energy-based micromodulus correction method for fracture in particle reinforced composites. Communications in Computational Physics, 2022, 32(2): 424–449

[46] Li Z, Han F. The peridynamics-based finite element method (PeriFEM) with adaptive continuous/discrete element implementation for fracture simulation. Engineering Analysis with Boundary Elements, 2023, 146: 56–65

RIGHTS & PERMISSIONS

Higher Education Press