An investigation on prevalent strategies for XFEM-based numerical modeling of crack growth in porous media

Abstract

Crack growth modeling has always been one of the major challenges in fracture mechanics. Among all numerical methods, the extended finite element method (XFEM) has recently attracted much attention due to its ability to estimate the discontinuous deformation field. However, XFEM modeling does not directly lead to reliable results, and choosing a strategy of implementation is inevitable, especially in porous media. In this study, two prevalent XFEM strategies are evaluated: a) applying reduced Young’s modulus to pores and b) using different partitions to the model and enriching each part individually. We mention the advantages and limitations of each strategy via both analytical and experimental validations. Finally, the crack growth is modeled in a natural porous media (Fontainebleau sandstone). Our investigations proved that although both strategies can identically predict the stress distribution in the sample, the first strategy simulates only the initial crack propagation, while the second strategy could model multiple cracks growths. Both strategies are reliable and highly accurate in calculating the stress intensity factor, but the second strategy can compute a more reliable reaction force. Experimental tests showed that the second strategy is a more accurate strategy in predicting the preferred crack growth path and determining the maximum strength of the sample.

Graphical abstract

Keywords

numerical modeling / extended finite element method / porous media / crack growth / stress intensity factor

1 Introduction

Crack, and porosity are ubiquitous features in geomaterials, mainly controlling the physical and mechanical behaviors of rocky environments and corresponding engineering works. Any release of in situ stress is concentrated on the tip of these discontinuities, and micro-cracks are initiated when they overcome the rock resistance. Micro-cracks grow progressively, and after coalescence to each other, they emerge as a fracture plane and cause a failure in the rock, which is a critical issue from an engineering point of view. Fracture mechanism has been investigated both experimentally and numerically in rock mechanics for various applications such as tunneling, blasting, hydraulic fracturing, earthquake mechanics, and slope stability analysis [ 1- 3].

In fracture mechanics, numerical methods are divided into discrete, continuum, and hybrid methods. Since fracture characteristics are highly variable, selecting an efficient method is not a straightforward task and may vary in each problem [ 4]. In discrete methods such as Distinct Element Method (DEM), the rock medium is considered an assemblage of individual particles connected by combining elastic law and Coulomb friction [ 5]. On the contrary, rock is supposed to be a continuous body in continuum methods [ 6]. Among all continuum methods, the extended finite element method (XFEM), a derivative of the standard Finite Element Method (FEM), has gained many attractions in fracture mechanics, because it has resolved the FEM’s fundamental weakness, the inability to model discontinuities and singularity of the crack tip. The presence of porosity in the rock mass greatly affects its strength and failure [ 7]. Pre-existing cracks and pores in natural rocks reduce their mechanical properties and strength [ 8]. Porosity as a defect is an influential factor in the crack growth mechanism and, according to the location and angle of porosity to the crack, it changes the stress distribution, which can either increase or decrease the maximum strength of the rocks [ 9].

Numerical modeling of crack growth is not direct implementation of the XFEM and rarely leads to reliable results on a porous sample. Therefore, two different approaches, or in this study, let say strategies, could be theoretically used to model crack propagation [ 10]. The first strategy is assigning a loose material property to the pores with Young’s modulus: E _void << E _background. However, in the second strategy, multiple partitions are applied and, then, different enrichment regions are assigned with independent crack propagation.

Many researchers have used these two strategies to investigate the crack growth procedure, especially in rocks as porous media. For example, Duarte et al. [ 11, 12] investigated the mechanical behavior of concrete in both cubic and cylindrical specimens under uniaxial pressure and indirect tensile stress using the finite element technique. In the first step, they designed the overall geometry of the model using image processing algorithms, and, in the second step, they simulated the material heterogeneity distribution by applying the first strategy for XFEM modeling. They finally compared the numerical results with the experimental data and suggested that the XFEM could be a suitable numerical tool to precisely model the crack initiation and propagation. Supar et al. [ 13] modeled the effect of pores arrangement on the stress distribution in the composite plates using the XFEM. They used the second strategy in their research to model the pores and found that a combination of thicker and non-staggered pore arrangement gives a relatively better prediction than other multi-pore arrangements. Rezanezhad et al. [ 14] used the second strategy to explore pore shape effect on simulation. They found equeivalent ellips are reliable shapes to be used instead of real comples pores.

Hedjazi et al. [ 15, 16] simulated the crack growth in dense, brittle materials. They used the second strategy to investigate the crack propagation and deflection of the crack toward the porosity. Chen et al. [ 17] modeled the effect of porosity on the crack behavior in a porous ceramic. They used the second strategy to represent the pores and examined the stress intensity factor changes with increasing the pore-crack distance. In another study, Rodriguez-Florez et al. [ 18] numerically simulated the effect of porosity on the crack growth in bone using the XFEM. They used the both strategies to simulate porosity and found that the crack growth path and reaction force differed. As a recent work, the same authors (Rezanezhad et al. [ 19]) investigated the effects of pore-crack relative location on crack propagation based on the second strategy. Rabczuk et al. [ 20] introduced a simplified mesh-free method for randomly evolving cracks. In this method, the crack can be arbitrarily oriented, but its growth is shown by the activation of crack surfaces in individual particles and, therefore, there is no need to display the crack topology. In another study, Rabczuk et al. [ 21] examined the applicability of the previously proposed method for brittle fracture in three dimensions and compared the results with experimental studies. Also, they [ 22] presented a new robust and efficient approach for modeling discrete cracks in meshfree methods where a set of cracked segments models the crack. The main advantage of the proposed method compared to the method in Ref. [ 20] was the elimination of additional unknowns in the equations to obtain displacement discontinuities. Ren et al. [ 23] proposed a stable solution to varying horizons based on which the balance of momentum and angular momentum in peridynamics (PD) are naturally satisfied. They also investigated the crack pattern of random point distribution and the issue of multi-materials on PD. It was shown that dual-horizon-PD is less sensitive to the spatial PD formulation than the original one for selected benchmark problems. Also, Ren et al. [ 24] introduced an explicit phase-field method for brittle dynamic fracture, and with numerical examples, they indicated that the explicit algorithm could improve the efficiency of the traditional phase-field model. In another study, they [ 25] presented a higher-order nonlocal operator method for solving partial differential equations and demonstrated the accuracy and capabilities of this method using several numerical examples. Areias et al. [ 26] studied damage and failure pattern using the screened Poisson equation and local re-meshing. They proposed a model for the propagation algorithm independent of any particular constitution and specific element technology. To evaluate the accuracy of the proposed algorithm, they used both quasi-brittle benchmarks and ductile tests.

Although the above methods are capable of handling crack branching and share the simplicity of interelement separation models with acceptable accuracies, they still do not have the power and accuracy of the XFEM method [ 20]. For example, although the peridynamics method can model the crack growth path well, it can only calculate the stress tensor [ 27], but not the stress intensity factor. Therefore, one cannot calculate important parameters in crack propagation, such as the stress intensity factor and the effect of porosity on that, with those methods. However, in this paper, these parameters are calculated by the two different strategies using the XFEM, and their cons and pros are discussed.

Previous studies have extensively analyzed the strength, deformation, and interconnection process of cracks in rock samples containing cracks and porosity, and the XFEM use in all of these studies demonstrates this method’s ability to simulate the crack propagation. However, considering the natural rock medium with the exact porous structure is still a challenge in the crack propagation problem, and the use of an appropriate strategy has not yet been comprehensively evaluated in the XFEM. The stress distribution field of pores can affect the crack growth path causing the new cracks creation, a deflection on the existing cracks, or even ceasing the crack propagation in some cases. In addition, recent developments in three-dimensional micro tomography computed (µ-CT) imaging of rock samples prepared more accurate view into microscale porous structure of rocks [ 28]. It is believed that the physical and mechanical behaviors of rock, such as fluid flow and wave propagation, are controlled by microstructures of rock [ 29]. Similarly, microcrack initiation and coalescence occur in microscale from crack tips, known as fracture process zone (FPZ) [ 30]. Numerical modeling of crack growth in such a medium and scale could lead to a deeper understating of crack propagation and its dependency on the porous structure of the rock. However, the first step of gaining such a goal is to find a valid and reliable numerical method and strategy.

In this paper, the XFEM reliability as a numerical method for crack growth modeling in porous media is evaluated both analytically and experimentally. For this purpose, the strategies mentioned above are used, and the advantages and limitations of each one are demonstrated. Then, by comparing the stress intensity factor and the reaction force of each strategy with the analytical and experimental results, their accuracy and reliability are discussed. Finally, the initiation and preferred growth path of cracks in a porous medium are implemented on digital images of Fontainebleau sandstone as a natural porous rock medium.

2 Theoretical basis

2.1 The extended finite element method

Generally, the XFEM is based on the partition of unity method in enriching the approximation of standard finite element components, which involves the effects of singularity or discontinuities of the domain around a crack. Enrichment in the XFEM is done by a combination of two techniques for cracks. Generally, the nodes affected by the crack body are enriched by the Heaviside function, and the points that are located at the crack tip are enriched by the enrichment function. In fact, by adding the enrichment functions, the degree of freedom of each node in the finite element domain increases [ 31, 32].

A key advantage of the XFEM is that the finite element mesh does not need to be updated to track the discontinuities like crack path. The displacement approximation in the XFEM is expressed as [ 33]:

(1) u=\sum _{i=1}^n{N}_i\left(x\right)\left[u_i+H\left(x\right)a_i+\sum _{\alpha }^4{F}_{\alpha }\left(x\right)b_i^{\alpha }\right] ,

where n is the number of nodes in the mesh, N _i( x) is the shape function of node i, u _i is the classical degrees of freedom of node i, and a _i and b_i^a are the degrees of freedom associated with the Heaviside step function H( x) and the crack-tip functions F _α, respectively [ 34- 36]. In Eq. (1), the first term is considered for all nodes in the model, and the common finite element approximation that already existed is used. What is very important in this equation and play key roles in the extended finite elements is the second and third terms of the equation where the discontinuities can be modeled on the body and the crack’s tip, respectively.

The Heaviside function H ( X) is defined as follow:

(2) H\left(x\right)=\{\begin{array}{l}+1,\left(x-x^{\ast }\right)e_{\mathrm{n}}⩾0,\\ -1,\mathrm{o}\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{w}\mathrm{i}\mathrm{s}\mathrm{e},\end{array}

where x is a sample point, x^{\ast } is the point on the crack closest to x and e _n is the unit outward normal to the crack at x* [ 31, 33, 34, 37].

The crack-tip functions F _α( x) provide improved accuracy and are required if the crack-tip terminates inside an element.

(3) \begin{array}{rl}& \left[F_{\alpha }\left(x\right),\alpha =1-4\right]=[\sqrt{r}\sin \left(\frac{\theta }{2}\right),\sqrt{r}\cos \left(\frac{\theta }{2}\right),\\ & \sqrt{r}\sin \left(\theta \right)\sin \left(\frac{\theta }{2}\right),\sqrt{r}\sin \left(\theta \right)\cos \left(\frac{\theta }{2}\right)],\end{array}

where r and θ are the polar coordinates at each crack-tip.

The XFEM relies on the level-set method to describe the crack topology and its propagation. The functions \varphi \left(x\right) and \psi \left(x\right) are used to accurately describe the location of cracks as presented in Table 1.

The nodal value of \varphi \left(x\right) indicates the distance from the node to the crack surface. Hence, where the function \varphi \left(x\right) is zero indicates the crack face. The function \psi \left(x\right) represents the distance of nodes to the line perpendicular to the crack surface. Finally, the intersection of \varphi \left(x\right)=0 and \psi \left(x\right)= 0 shows the crack tip [ 10, 38].

Gauss quadrature law provides precise numerical solutions for polynomial integrals. Unfortunately, this is not the case for enriched finite element solutions because in the XFEM, cracks can be placed in an element by enrichment functions, and they contain very nonlinear and even discontinuous functions. In these cases, integration by conventional methods such as the use of Gaussian laws does not work, and there is no necessary accuracy in integration. Here, the element which contains cracks is divided into several parts in order to integrate. The method of division is to divide the element into sub-triangles or sub-quads. Both methods are suggested by [ 39]. Each part of the element on both sides of the crack is divided into a number of triangles or squares according to the selected method, and in each of them, Gauss’s law is applied for integration [ 32, 39]. From a different point of view, in the XFEM, the enrichment process is done locally because only a series of nodes are enriched to model an entire domain. So, elements may be classified into three categories: the region where all element nodes are enriched, a region where enrichment does not occur, and the third category called the blend region in which only some nodes of an element are enriched. The blend region is very important. Since the third type of element is incompletely enriched, they do not satisfy the conditions of the partition of unity method.

For this reason, it behaves inappropriately, and special attention should be paid to this region in order to obtain correct and accurate answers. The modified XFEM, introduced by Ref. [ 39], alters the global enrichment function in order to avoid these problems. Using shifted enrichment automatically removes the enrichment from the domain which is not required to be enriched. Thus, using the shifted enrichment within XFEM not only makes easy calculation of node displacement that simplifies the problem and formulation, but also helps in applying the Dirichlet boundary conditions as the enriched degree of freedom or, more precisely, the enrichment is zero at the node [ 32, 40, 41].

2.2 Modeling pores in extended finite element method

Due to defects such as porosity, cracks, joints, etc., precise analysis is required to investigate its strength and stability in a structure. Since porosity is one of the critical components of a structure and strongly affects its strength, it should be included for accurate modeling. Porosity modeling requires finite element meshes that correspond to the geometry of the porosity surface. This is both a time-consuming and computationally expensive approach. The first proposed strategy is based on this approach; however, in the second strategy, the voids (or pores) in the model play the role of porosities. Here, XFEM offers a great way to model arbitrary discontinuities, where the mesh does not require to be aligned with the boundaries of pores or material interfaces. This is done by combining the appropriate enrichment function in the finite element approximation space where these functions are matched using thepartition of unity method. As a result, the expanded approximate space or an enriched finite element can approximate the field variable at a lower computational cost. When a pore is created in an enriched domain, the Heaviside function is used because of its properties. In addition, a level-set function \varphi \left(x\right) is built on the domain to define the discontinuity geometry, evaluate the enrichment performance, and select the elements intersected by the discontinuity surface. In general, it can be said that the nodes that are outside the pore would have H \left(X\right)= 1 and the nodes inside the pore, H \left(X\right)= 0. In the next step, nodes inside the pore and nodes, which are not cut by the pore, are removed from the calculations. This is usually done by removing the degrees of freedom associated with those nodes from the system of equations, and the system is solved only for the remaining degrees of freedom.

2.3 XFEM weak formulation

A body C is considered with domain Ω ⊂ R ² and a surface is denoted by Γ and \overline{t} is the tension applied to the surface Γ _t [ 42]. Dirichlet boundary conditions are applied to the surface Γ _u such that Γ _t⊂ Γ, Γ _u⊂ Γ, and the body contains a pore as Γ _h. If U be the displacement field which should be approximate, we define the functional spaces in which we search for the solution. The acceptable displacement field is defined as follows:

(4) u\in U=\left\{uu\in H^1\left(\mathrm{\Omega }\right),u=\overline{u}\mathrm{o}\mathrm{n}{\mathrm{\Gamma }}_{\mathrm{u}}\right\},

where \overline{u} is the displacement prescribed and the test function is given as:

(5) w\in W=\left\{ww\in H^1\left(\mathrm{\Omega }\right),w=0\mathrm{o}\mathrm{n}{\mathrm{\Gamma }}_{\mathrm{u}}\right\}.

Then the weak form is obtained by multiplying the differential equation, that is, the strong form of the equilibrium equation by the weight function w, and then the integration in the domain Ω. Hence, the weak form can be described as u ∈ U such that ∀ w ∈ W [ 42]:

(6) {\int }_{\mathrm{\Omega }}\left[\mathrm{\nabla }w:\sigma \left(u\right)\right]\mathrm{d}\mathrm{\Omega }-\underset{{\mathrm{\Gamma }}_{\mathrm{t}}}{\int }w\cdot \overline{t}\mathrm{d}\mathrm{\Gamma }=0.

2.4 Analytical solution for interaction between pore(s) and crack

Consider a wide plate with an arbitrary stress distribution due to existing cracks, circular or elliptical pores. The analysis is somewhat complicated but similar to that of an infinite body. At infinity, the plate is assumed to be subjected to the following stresses:

(7) {\sigma }_x^{\mathrm{\infty }}=\sigma \left(\alpha +\frac{\delta Y}{d}\right),{\sigma }_y^{\mathrm{\infty }}=\sigma \left(\beta +\frac{\mu X}{d}\right),{\tau }^{\mathrm{\infty }}=\sigma \gamma ,

where \sigma is the applied stress, d is the pore distance from the center and \alpha, \beta, \gamma, \delta, and \mu are the constant values as described in Ref. [ 43].

The stress function is written as a summation:

(8) \chi ={\chi }_{{}_0}+ \sum _{k=1}^N{\chi }_k,

in which {\chi }_0 corresponds to the stress at infinity.

(9) {\chi }_{{}_0} =\sigma d^2\mathrm{R}\mathrm{e}\left\{\overline{z}{\varphi }_0\left(z\right)+{\psi }_0\left(z\right)\right\},

and

(10) {\varphi }_0 (z)=\frac{1}{4} (\beta +\alpha )z+\frac{1}{8}(\mu -i\delta )z^2,

(11) {\psi }_0 (z)=\frac{1}{4} (\beta -\alpha +2i\gamma ) z^2+\frac{1}{24}(\mu +i\delta )z^3.

The stress functions X_k\left(k=1,2,...,N\right) possess singularities within the pores and can be expanded in the following Laurent series:

(12) {\chi }_k =\sigma d^2\mathrm{R}\mathrm{e}\left\{\overline{z}{\varphi }_0\left(z\right)+{\psi }_0\left(z\right)\right\},

where

(13) {\varphi }_k \left(z_k\right)=\sum _{n=0}^{\mathrm{\infty }}\left(F_{n,k}^{\bullet }+i{F}_{n,k}^{\prime }\right)z_k^{-\left(n+1\right)},

(14) {\psi }_k \left(z_k\right)=-D_{0,k}^{\bullet }\log z_k+\sum _{n=1}^{\mathrm{\infty }}(D_{n,k}^{\bullet }+i D_{n,k}^{\prime })z_k^{-n}.

The stress function (Eq. (8)) is chosen to satisfy the boundary conditions at infinity, and thus only those cracks located along the pore edges are considered. To discuss about the jth crack, the total stress function is expressed using Eqs. (8) to (16) with the only z_j variable.

(15) z=z_j{\mathrm{e}}^{i\alpha j}+ r_j{\mathrm{e}}^{i\beta j},

(16) z_k=z_j{\mathrm{e}}^{i\left({\alpha }_j-{\alpha }_k\right)}-r_{jk}{\mathrm{e}}^{i\left({\beta }_{jk}-{\alpha }_k\right)}.

The results are written as follow:

(17) \chi =\sigma d^2\mathrm{R}\mathrm{e}\left\{{\overline{z}}_j{\varphi }_j\left(z_j\right)+{\psi }_j\left(z_j\right)\right\},

(18) {\varphi }_j \left(z_j\right)=\sum _{n=0}^{\mathrm{\infty }}\left\{\left(F_{n,j}^{\bullet }+i{F^{\prime }}_{n,j}\right)z_j^{-\left(n+1\right)}+\left(M_{n,j}^{\bullet }+i{M}_{n,j}^{\prime }\right)z_j^{n+1}\right\},

(19) \begin{array}{rl}{\psi }_j\left(z_j\right)=& -D_{0,j}^{\bullet }\log z_j+\sum _{n=1}^{\mathrm{\infty }}(D_{n,j}^{\bullet }+i{D}_{n,j}^{\prime })z_j^{-n}+\\ & \sum _{n=0}^{\mathrm{\infty }}\left(K_{n,j}^{\bullet }+i{K}_{n,j}^{\prime }\right)z_j^{n+2}.\end{array}

It is possible to find some relations between coefficients that reduce the stress function. The stress intensity factors of the crack tip X=a_j are given by the following series:

(20) \begin{array}{rl}{\left(k_1-i{k}_2\right)}_{A_j}& =2{\left(2d\right)}^{\frac{1}{2}}\underset{z_j\to {\lambda }_j}{lim}\left[{\left(z_j-{\lambda }_j\right)}^{\frac{1}{2}}{{\phi }^{\prime }}_j\left(z_j\right)\right]\\ & =\sigma {\left(a_j\right)}^{\frac{1}{2}}\left(F_{1,j}-i{F}_{2,j}\right),\end{array}

(21) F_{1,j}=\sum _{n=0}^M{C}_{n,j}^{\left(1\right)} {\lambda }_j^n,F_{2,j}=\sum _{n=0}^M{C}_{n,j}^{\left(2\right)} {\lambda }_j^n.

The corresponding values for the opposite crack tip x_j=- {\lambda }_j may be obtained upon replacing {\varphi }_j \left(z_j\right) by {\varphi }_j\left(-z_j\right) and {\lambda }_j by \left(-{\lambda }_j\right) in Eq. (17). The coefficients C_{n,j} are evaluated exactly by closed-form expressions. Typical stress intensity factor calculation results analytically using the above analysis proposed by Sih [ 43].

Using the above analytical results and graphs, the stress intensity factor at the crack tip under tensile loading can be obtained. One of the objectives of this paper is to evaluate and compare the accuracy of the proposed strategies to calculate the stress intensity factor. It is one of the critical parameters in fracture mechanics, and its value at the crack tip varies depending on the condition and the location of the porosity. In numerical models, the stress intensity factor is calculated from both strategies and compared with the analytical method, as explained in Section 3.2.

In the XFEM method, the calculation of the stress intensity factor is possible only in three-dimensional mode, so all the models are simulated as 3D models. However, due to the complexity of partitioning in XFEM, a large number of elements, and the long run time, the size of the third dimension is considered less than two other dimensions to reduce the computational costs. In these models, a plan stress condition is assumed, so the results are close to the 2D condition. To obtain accurate answers and correct comparison between the two strategies, in calculating the stress intensity factor, the initial crack length is critical, and the effect of the infinite plate must be observed [ 44]. In other words, the crack size should be small compared to the size of the sample so that the crack tip is not affected by external boundaries. Consequently, in all models, the initial crack size is approximately 1/10 of the main dimensions of the sample.

3 Comparison of two different strategies

There are some difficulties in simulating the growth of the cracks among several porosities using the XFEM method. In the enrichment zone where the criterion of damage has been applied, the initiation and growth of a new crack occurs when all pre-existing cracks are fully grown. Therefore, several cracks may grow suddenly after the porosity of the sample. This phenomenon is illustrated in Fig. 1, where a specimen containing a crack and porosity is subjected to tensile stress in the horizontal direction. The whole plate is enriched only as a region. The damage criterion is defined based on the maximum principal stress, and its value is 1 MPa. This means that after loading, new cracks must appear in elements where the stress is greater than 1 MPa and has reached the critical maximum principal stress. In frame 15 (Fig. 1(a)), the maximum principal stress at the initial crack tip is 1.108 MPa and reaches a critical value, and the crack starts to grow so that the upper crack tip hits the porosity. In frame 44 (Fig. 1(b)), the principal stress around the porosity exceeds the defined value and reaches 1.482 MPa, but no new crack initiates in this region.

On the other hand, the initial crack propagation continues until it reaches the lower edge of the sample, and the stress around the porosity is still increasing, and its value reaches 53.791 MPa (86 frames in Fig. 1(c)). In fact, in such conditions, porosity is considered the boundary of the enriched area. Thus, after the full growth of the initial crack, new cracks starts to grow. This causes a large region on the other side of the porosity where the elements have reached a critical maximum of the principal stress, resulting in the sudden appearance of a high number of cracks in frame 87 (Fig. 1(d)) with a stress value of 54.311 MPa.

To overcome this problem in this section, two different strategies are investigated for applying the XFEM-based crack propagation modeling in porous models. We aim to highlight each strategy’s advantages and limitations in calculating the stress intensity factor and reaction force or, more importantly, estimation of preferred cracks growth path. The point is that these strategies affect differently on each of these parameters. In what follows, we try to deeply dig such effects both analytically and experimentally, and obtain more applied criteria for a valid and purposeful crack growth modeling.

3.1 A quick look at advantages and limitations

As mentioned in Section 2.1, there are two different strategies to model crack propagation. The first strategy is assigning loose material properties to the pores, and the second one is applying multiple partitions to create different enrichment regions.

As the main limitation of the first strategy, multiple cracks growth could not be modeled in this strategy. The sample is more affected by the initial crack using the first strategy. It means, the whole model is considered as one enriched area, and no new crack begins to grow until the primary crack propagates through the whole model and reaches the model’s boundaries. Meanwhile, although some nodes may have a stress value more than the maximum principal stress, no failure initiates from those locations since the initial crack has not yet been fully developed. The effect of this issue is illustrated in Fig. 2, where a 2D sample with three pores and one crack is modeled, and the crack path is simulated using two strategies. According to this figure, samples are under tensional stress, and in the first strategy (Figs. 2(a) and 2(b)), only the initial crack propagates through the sample, and no other crack is generated even if there are nodes with critical stress values (as marked by arrows). However, in the second strategy (Figs. 2(c) and 2(d)), in addition to the initial crack, new cracks are initiated and propagated from the regions that reach the maximum principal stress.

This is one of the main limitations of the first strategy, which can strongly affect the accuracy of the final results. In this strategy, the possible error rate increases, with the increment of the number of the pores. Fortunately, this is not the same in the second strategy, where several different enrichment areas are generated. Contrary to the previous strategy, multiple cracks can grow simultaneously in the sample. However, this strategy has its own limitations. For example, where several cracks are propagated, two different cracks cannot coalesce in one element. In this case, the stress concentration in the element would be increased around that node, as shown in Fig. 3. On the other hand, the next limitation goes back to the user itself. The whole model should be split up and enriched into multiple sections, which is a very time-consuming and astringent task, especially when the model is complicated (see Section 4).

Although these results demonstrate that the output of each strategy can be completely different, it is noticeable that the main crack path is quite similar in both strategies.

3.2 Computation of the stress intensity factor and the reaction force

The stress intensity factor can be obtained in an arbitrary model of pores and cracks by solving the corresponding equations mentioned in Section 2.2. Sih [ 43] solved those equations for some predefined models and presented the results as plots. To examine the reliability of each strategy, we built similar samples as indicated in Fig. 4 and computed the stress intensity factor on the crack tip to compare them with the analytical results.

The physical and mechanical properties of these samples are given in Table 2, which is related to the Granite. They are used since our experimental samples are made of granite in the following sections (i.e., Section 3.3.2). We also explain the reasons for using this rock. Both samples are subjected to the tensile load horizontally.

In the first strategy, Young’s modulus assigned to the pore is 1000 times weaker than the other parts of the model. In both samples, the crack length, pore radius, and distance between pore(s) center and crack center are 12, 7.5, and 30 mm, respectively. So, their analytical modeling properties are a/ d = 0.2 and ρ/ d = 0.25, in which 2a is the length of the crack, d is the distance from the center of the crack to the center of the pore, and ρ is the pore radius (see Sih [ 43]).

The results, computed numerically from the XFEM with both strategies and obtained analytically from Sih [ 43], are compared in Table 3. As indicated in Table 2, the output results from both strategies are very close to the analytical solution in both samples. Therefore, both strategies are acceptable in the case of stress intensity factor estimation.

Although the results are promising, it should be noted that the choice of mesh size is generally a critical factor controlling the accuracy of the results in numerical continuum methods. Choosing the right element type and creating the optimal mesh size of elements would greatly impact the accuracy and convergence in the finite element analysis. An inappropriate element size could lead to considerable errors in the results of numerical modeling. On the other hand, the small size of the mesh increases the analysis time, accumulates errors, and finally induces stiffness in the model. Generally, convergence to the valid results has to be achieved for the validity of a finite element solution, which means that the answers must not change drastically if the element size becomes smaller. In our cases, two models have been solved using relatively coarse elements. It should be noted that the elements’ coarseness is relative and depends on the geometrical dimensions of the model constituents, details of the problem (such as crack and pore existence, etc.), and personal experiences. As a common approach, after completing the first solution, the size of elements is halved, and the problem is solved again. This process continues as long as the results from two consecutive solutions differ by less than 10%. Therefore, in our cases, the stress intensity factors were computed again by dividing the element size in half to compare the accuracy of the final solution results. As shown in Table 4, by half size of the previous elements, the obtained results changed less than 1% compared to the primary ones. According to the results in this section, it can be concluded that the XFEM is an accurate and reliable method to obtain the stress intensity factor comparable to analytical solutions.

The von Mises stress distribution of each sample simulated by two different strategies is illustrated in Fig. 5. As it can be seen in this figure, the stress distribution is quite identical in both strategies, and the major stresses concentrate between the crack tip(s) and the pore(s). In the first sample (Figs. 5(a) and 5(b)), the maximum stress is about 4.63 and 4.72 MPa in the first and second strategies. The second sample has one additional pore compared to the first sample, it is theoretically expected that the maximum stress increases. However, as recorded in Fig. 5(c), the maximum stress decreases to 4.16 MPa in the second sample modeled with the first strategy. On the other hand, the result of the second strategy is compatible with the theoretical expectation, and the maximum stress increases to 4.88 MPa, as shown in Fig. 5(d).

In the XFEM method, there are five criteria for creating cracks. Here we used the criterion of maximum principal stress for both strategies. In this method, when the maximum principal stress in the element is greater than its critical value, failure occurs in part.

(22) f=\left\{\frac{⟨{\sigma }_{max}⟩}{{\sigma }_{max}^o}\right\}.

Generally, the XFEM method simulates post-failure behavior in both displacement and energy. Different values for failure under normal or shear stresses cannot be defined in the displacement method, and only the displacement value after the initiation of failure must be defined. In the energy method, different behaviors and strengths in terms of normal and shear stress in different directions can be defined to simulate the failure behavior of the element. Therefore, this is a more accurate model to study the behavior of the model after failure. If different modes are independent, the amount of fracture energy is defined by the following formula: the ratio of the fracture rate of the element to the critical fracture rate of the element. When this ratio is greater than one, the crack tip is deboned [ 49, 50].

(23) f=\frac{G_{\mathrm{e}\mathrm{q}\mathrm{u}\mathrm{i}\mathrm{v}}}{G_{\mathrm{e}\mathrm{q}\mathrm{u}\mathrm{i}\mathrm{v}\mathrm{C}}}⩾ 1.

In the XFEM method, there are three criteria for calculating the fracture energy rate. The method used for both strategies in this paper is the power method, which is defined as the ratio between the final equivalent energy and the critical final equivalent energy.

(24) \frac{G_{\mathrm{e}\mathrm{q}\mathrm{u}\mathrm{i}\mathrm{v}}}{G_{\mathrm{e}\mathrm{q}\mathrm{u}\mathrm{i}\mathrm{v}\mathrm{C}}}={\left(\frac{G_{\mathrm{I}}}{G_{\mathrm{I}\mathrm{C}}}\right)}^{a_{\mathrm{m}}}+{\left(\frac{G_{\mathrm{I}\mathrm{I}}}{G_{\mathrm{I}\mathrm{I}\mathrm{C}}}\right)}^{a_{\mathrm{n}}}+{\left(\frac{G_{\mathrm{I}\mathrm{I}\mathrm{I}}}{G_{\mathrm{I}\mathrm{I}\mathrm{I}\mathrm{C}}}\right)}^{a_{\mathrm{o}}},

where G _IC, G _IIC, and G _IIIC are Mode I, II, and III fracture toughness of the delaminated material. Three coefficients a_{\mathrm{m}},a_{\mathrm{n}},a_{\mathrm{o}} are related to the parameters of the material which are considered the same in the definition of matter [ 51].

In Table 2, where the properties of materials are defined, Young’s modulus and Poisson’s ratio are the elastic properties of matter, and {\delta }_{\mathrm{t}} is tensile strength. However, since in fracture mechanics, it is assumed that each component has a series of defects, only the stress parameter is not sufficient to express the tolerance of a member, and stress intensity factor should be used. In general, the stress intensity factor in an infinite body containing crack in length 2 a, which is under tension \sigma is computed as:

(25) K_{\mathrm{I}}=\alpha \sigma \sqrt{\pi a}.

During loading, if stress reaches a critical value ( {\sigma }_{\mathrm{c}}), the critical stress intensity factor is obtained. The critical stress intensity factor is called the fracture toughness. This parameter represents a property of matter because it is the stress intensity factor at which a particular body exposed to tensile loading would break. That is, the stress intensity factor is equal to the critical value [ 52]:

(26) K_{\mathrm{e}\mathrm{q}}=K_{\mathrm{I}\mathrm{c}}.

Since we are going to simulate crack growth and propagation, in addition to density, Poisson’s ratio, and Young’s modulus, information about failure, damage, and material loss should also be defined. Using the maximum principal stress criterion, the value presented in the Table 2 is assumed to be the beginning of failure in the sample, and as shown in the stress-strain curve (Fig. 6(a)), this value is specified at point D = 0, where tension is equivalent to {\sigma }_{\mathrm{y}0}. Because the XFEM method is used and the crack tip can be placed inside the element, if no additional information is defined, the element’s behavior after reaching the final stress will follow the dashed line in Fig. 6(a). After the tension reaches the final limit, some energy is absorbed by the material before the complete failure. After reaching the yield strength, the area under the stress-strain curve is the fracture energy required to propagate the crack ( Gc). This range of material behavior is called the damage range. At the beginning of the damage range, D is 0 (Fig. 6(b)), and at the end, when the damage is complete, the value of D is 1. Therefore, the energy mentioned in Table 2 is the same area under the stress-strain curve, and its value is the critical amount of fracture energy needed to complete the failure process and crack growth and propagation. Therefore, as the amount of fracture energy increases, the sample’s energy absorption goes up and increases the maximum resistance of the sample.

To have a more detailed comparison of these strategies, force–displacement curves are also computed and investigated, as posted in Fig. 7. Generally, the maximum force, these 3D models could resist before failure, decreases by increasing the number of pores in front of a crack [ 9, 19]. Figure 7(a) shows the force-displacement curves in the first strategy. According to the results, the maximum force that the sample can tolerate with two pores is more (about 14%) than the sample with just one pore. This result could not be theoretically correct and is explained in the following text. However, both curves in Fig. 7(a) are similar overall, with force values as 170.94 and 199.06 N for the first and second samples, respectively (see Table 5). Force-displacement curves in the second strategy are also plotted in Fig. 7(b). In this figure, the reaction force value in the first sample with one pore equals 156.83 N, decreasing to 137.27 N in the second sample with two pores (see Table 5).

Table 5 demonstrates changes in the maximum stress and force-displacement curves for models in Figs. 11 and 13 using first and second strategies. Based on these results, in addition to the advantages and limitations of each strategy mentioned in Section 3.1, the second strategy seems to be the more accurate strategy to evaluate von Mises stress and reaction force.

For more detailed investigation, the crack growth steps in both samples using the first strategy are shown in Fig. 8. As shown in Fig. 8(a), in the first sample, the crack begins to grow until frame 35, where it reaches the pore. Although the regions after the pore are affected by a stress more than the maximum principal stress in frame 36 (indicated by yellow to red areas on top of the pore), no new crack has grown up in these excited areas as the growth of the previous crack has not yet been completed. This condition continues until frame 44, where the crack begins to grow upward from those previously excited areas.

Similarly, the same conditions happen for the second sample, as reported in Fig. 8(b). In this sample, the crack grows until frame 19, where it reaches the pores. In frame 21, the regions around the pores are affected by a stress more than the maximum principal stress, but no new cracks have grown up until the growth of the previous crack has been completed in frame 31. In the next step, the crack begins to grow from excited areas. This is why the excited areas that cannot grow before reaching the primary crack have higher stress in the second sample and so, the reaction force wrongly increased compared to the first sample. This problem is another weakness of the first strategy, which can significantly affect the accuracy of the final results. Unfortunately, by increasing the number of pores in the sample, this problem becomes even more severe, and so, the first strategy cannot correctly simulate the behavior of the real porous media (see Section 4).

In the continuous environment, the equation governing the behavior of matter is a nonlinear differential equation with no general solution, leading to be solved by discretization methods. The first step in discretizing and solving these equations is to create proper networking in the model. This networking is called meshing. The default element selection process is determined and selected in finite element methods by determining the problem conditions during modeling. Finite element analysis can be considered as including four basic steps [ 53]: dividing the analysis area into a limited number of elements, extracting the equations governing the behavior of each element, assembling the elements, and solving the resulting equation system. In this method, discretization of the physical domain of the problem, which is a continuous system with infinite degrees of freedom, is one of the primary sources of errors. Discretization error or model order error depends on the number and order of elements selected for modeling. The contribution of the error due to the number of elements can be minimized by increasing the mesh density (increasing the number of elements). To reduce the error due to the selection of the shape function, the order of the elements must be increased. This also increases the time required for analysis. Therefore, selecting the appropriate elements (quadrilateral, hexagonal, etc.) and their size is vital in three-dimensional finite element analysis. To calculate the stress intensity factor in the XFEM method, only regular discretization and the use of hexagonal elements would lead to accurate answers. Because the selected element must present the plane stress behavior well, and since the use of high-order elements is preferable to the low-order type (for a certain number of nodes), the hexagonal element can model the stresses very well due to loading and boundary conditions [ 54, 55]. Therefore, in case of irregular discretization and increasing the computational costs, the results are not consistent with the analytical solution, and the amount of stress is estimated to be lower than the actual value. Comparing the maximum strength obtained from regular and irregular discretization in the sample contains a crack, and the porosity posted in Figs. 9 and 10. Figure 9 shows the irregular discretization in both strategies.

According to Fig. 10, a more resistant model is estimated in irregular discretization, which is evident in both strategies. So that in the first strategy, the maximum resistance is 189.19 MPa (Fig. 10(a)) and suggests an increase of 10.6% of regular discretization (Fig. 7(a)). In the second strategy the maximum resistance is 164.95 MPa (Fig. 10(b)) and posts an increase of 5.1% (Fig. 7(b)). It should be noted that discretization in the XFEM method does not affect predicting the crack propagation path.

3.3 Estimation of crack growth path

Besides those limitations mentioned in previous sections, estimation of crack growth path is a challenging issue in crack growth modeling. We aim to evaluate these strategies to see how accurately they can predict the crack path via analytical and experimental validations.

3.3.1 Analytical validation of strategies

To investigate and compare the accuracy of the first and the second strategies in crack growth modeling of a porous medium, the analytical method proposed by Hedjazi et al. [ 15, 16] has been used. These results are also supported by experimental analysis. They investigated how the crack propagates and deviates toward a pore in dense, brittle materials using laboratory samples as well as numerical simulations. Figure 11 depicts the specifications of the laboratory sample, which we similarly used in XFEM modeling by two different strategies.

Figure 12 provides the comparison of the crack growth path in the laboratory sample reported by Hedjazi et al. [ 15] and our numerical modeling by each strategy. As can be seen in this figure, both modeling strategies have been able to predict the crack growth path well, and they are in good agreement with the experimental results. It should be noted that the crack branching is not considered numerically since the crack branching tracking requires an additional criterion such as a critical crack velocity [ 15, 33, 56, 57], which is not the case in this paper. Numerical simulation showed that the maximum von Mises stress equals to 3.38 and 3.41 MPa in the first and second strategies.

To investigate and compare the deflection of the crack toward the pore more precisely, an analytical method introduced by Ref. [ 16] is used. Here, y( x) is the theoretical crack deflection function for a circular pore with a radius of r, which is located at coordinate ( a, b):

(27) y\left(x\right)=\frac{r^2}{2b} \left[2-t\left(2+t-t^2\right)\right],

where t=b-x/\sqrt{b^2+\left(a-x^2\right)} and x is coordinate of the crack tip along y-axis.

Figure 13 illustrates the analytical and numerical (with two strategies) results of the crack deflection toward the pore with two different pore sizes. Here, the dimensions of the samples are equal to previous samples (i.e., 10 mm × 10 mm × 1 mm), the initial crack length is 1 mm and pore located at \left(a,b\right)=\left(2.57,2.02\right) with two different sizes of 2 and 1.7 mm.

As presented in Fig. 13, the results of the crack deviation to pore in both XFEM strategies are in good agreement with the analytical method. However, it seems the first strategy can estimate the crack deflection more precisely, and so this strategy could be more reliable in predicting preferred crack growth path (see Section 4). Also, these plots (Figs. 13(a) and 13(b)) show that the difference in analytical and numerical results in the calculation of the crack deflection to the pore has been reduced by decreasing the pore size in both XFEM strategies. This result becomes even more important in real rock media where the pore sizes are much smaller than the existing cracks, and therefore, the modeling results would be closer to the real samples in Section 4.

However, it should be noted that in the XFEM method, regular and irregular discretization does not affect the prediction of the crack propagation path. In regular discretization, the maximum value of critical stress is calculated in the middle of the element where the crack tip is located. On the other hand, in irregular discretization, the critical stress values are calculated at the crack tip, so the elements should be increased as much as possible to calculate the displacement in the nodes more accurately, increasing the processing time. Figure 14 suggests the crack growth path in irregular discretization, which agrees with the laboratory results in both the first (Figs. 14(a) and 14(b)) and the second (Figs. 14(c) and 14(d)) strategies. Numerical simulation showed that the maximum von Mises stress equals to 2.85 and 3.10 MPa in the first and second strategies. As mentioned before, these values are lower than the stress values in regular discretization.

3.3.2 Experimental validation of strategies

In this section, we compare two strategies of the XFEM with experimental tests for the sake of the crack growth path. However, due to the lack of high-tech facilities, conducting such an experiment in the pore scale is almost impossible in our case. Another choice is to resemble such a test in a macro scale. Therefore, we select a solid non-porous homogeneous rock as the background and, then make some holes and cracks on it with predefined pattern. Each hole plays the role of porosity in the rock with a known shape. Granite is one of the best choices to be used as the background rock in this study. It is a uniform assemblage of minerals with a mostly granular texture and very low porosity. Since the size of minerals is much smaller than the pores, it can be considered a homogenous medium. The only grain-to-grain contact may affect the crack path; however, it is neglected in this study for simplification. Granite is a relatively hard rock, which we used a water jet device to make the holes and cracks on it. Water jet can cut the rigid Granite with any pattern and size, such as a piece of cheese, without inducing micro-cracks in the sample.

Similar to analytical evaluation, two different Granite samples with length, width, and thickness of 80, 80, and 15 mm were prepared. The thickness in these samples has been chosen so that the jaws of the tensile machine can grip the sample without any damages due to the gripping pressure. The first sample contains one crack and pore in the middle of the sample, and the second sample consists of one crack in the middle of the sample and two pores on both sides of the crack, as described in Table 6 and shown in Fig. 15.

The STM-250 tensile device (Fig. 16) has been used to perform the tensile test and observe the crack propagation path in each sample. The sample is pulled away from both sides via gripping jaws, which are wider than the sample to generate a homogenous force. The device can implement tensile force either perpendicular or tilted to the sample. To induce more challenge in this study, we used an oblique tensile force with an angle of 30° relative to the sample edge. The samples were displaced by 3 mm from both sides. Crack growth is very unstable, and the failure of the sample happens instantaneously once the crack initiation begins. Therefore, we could barely check the crack path in the broken samples.

As shown in Fig. 15, two real samples were numerically modeled, as displayed in Fig. 17, with two different strategies of XFEM. Granite properties, presented in Table 2, were assigned to both models. The first sample is similar to Fig. 15(a), and the pore is located at its center. The distance between the middle of the crack to the pore center is 20 mm. This model is displaced by 3 mm from both sides by an angle of 30° (Fig. 17(a)). The second sample is also similar to Fig. 15(b). The crack located in the center of the sample and the distance between this point and the center of pores is 15 mm. Also, the distance of the center of the pores to the edges is 20 mm. This sample is displaced by 3 mm from both sides by an angle of 30° (Fig. 17(b)). In these samples, the maximum principal stress criterion is used as a failure criterion which expressed that the fracture is created when the primary stress in the element is more than the specific critical value.

Figure 18 shows the results of the crack growth in the experimental (Fig. 15) and the numerical models (Fig. 17). As shown in both samples, tensile stress is concentrated on tips of the primary crack, where new cracks begin to grow. Thus, cracks propagate through the samples almost perpendicularly to the tilted tensile axis [ 58]. As shown in Fig. 18(b), the first strategy in the sample with one pore is inefficient and failes to precisely predict the crack initiation point at the upper half of the pore. By increasing the number of pores in Fig. 18(e), the efficiency of this strategy decreases in predicting the crack growth path, and the result is not very similar to the experimental ones. On the other hand, Figs. 18(c) and 18(f) demonstrate that the second strategy can simulate the crack path similar to what happens in reality, and so, it is reliable in the predication of the crack growth path.

Since the simulated crack growth path in the second strategy is similar to experimental results, the stress distribution of both samples is posted in Fig. 19 to dig more on crack growth details. As it is observed from the first sample in Fig. 19(a)-frame 5, the stress concentrates on both sides of pore and crack tips. However, it reached the maximum critical stress level on the upper side of the pore and the crack tip near it. Thus, the initial crack starts to grow from these regions. As expected, the growth direction of these cracks is perpendicular to the tensile applied load. In frame 28, the other tip of the initial crack begins to grow while the previous cracks are moving in their earlier paths. When the crack is close to the pore, due to the influence of the pore on the stress field around the crack tip, it deflects toward the pore, as posted in frame 56. Finally, the initial crack reaches the pore in frame 96.

The crack growth in the second sample is also shown in Fig. 19(b). As seen in frame 4 of this figure, the tips of the initial crack begin to grow while the stress on the pores has not yet reached the maximum critical stress. In frame 13, the stress on the pores has reached the critical stress, and new cracks are created. Due to the effect of the pores on the initial crack stress field, this crack changes its way toward the pores, as described in frame 40. Finally, the initial crack reaches the pores in frame 79, and the final failure occurs. The overall path of the crack in this sample is also perpendicular to the tensile axis subjected to the sample. The crack initiation and propagation process from frame 4 to 79 is symmetric, due to a symmetric pattern of pores and crack and homogenous background.

On the other hand, in the first strategy, due to the influence of the initial crack on the preferred crack growth path, the excited areas around the pore(s) cannot create new cracks, and so, the final crack growth pattern does not fit well with the failure pattern of the experimental samples. As mentioned before (i.e., Fig. 8), the difference between the first strategy results and the actual answer increases by increasing the number of pores. Figure 20 shows the force-displacement curves of the first sample in both strategies and the experimental result. According to this figure, the maximum force that a specimen can withstand before failure in the laboratory is 1453.1 N. This value is less than the result of simulated models, which are 1593.6 and 1466.9 N in the first and second strategies, respectively. These differences are primarily due to the grain-to-grain contact, micro-cracks, and granular nature of the rock sample. However, the result of the second strategy is much closer to the experimental result.

4 A discussion on crack growth modeling in a real porous media

As a standard reservoir rock in digital rock physics, Fontainebleau sandstone is used in this study. Andrä et al. [ 59] and Madonna et al. [ 60] introduced it thoroughly; however, it dominantly consists of monodisperse Quartz grains. Laboratory analysis revealed that the Fontainebleau sandstone sample has an average porosity of about 15%. 3D µ-CT image of this sample contains 1024 × 1024 × 1024 elements with a voxel edge length of 0.74 μm.

In this study, a 2D image, extracted from the 3D cube, was used to crack growth investigation in a natural porous medium (Fig. 21(a)). This grayscale image was segmented using thresholding to produce binary images consisting of the pore (black areas) and mineral phases (white areas), as shown in Fig. 21(b).

One of the most essential parameters in image analysis methods is image resolution. Softening of the microstructure, edges, and sharp corners are some manifests of decreasing the image resolution. To investigate the effect of image resolution on crack growth in this sample (i.e., Fig. 21(b)), four different median filters that artificially reduced the original image resolution were used as shown in Figs. 21(c) to 21(f). Binary images were used as the base images to mae geometry of the corresponding XFEM models in both strategies. Each model was subjected to a displacement of 0.001 mm in the horizontal direction, and the Quartz properties, shown in Table 7, are assigned to them.

Although grain-to-grain contact of Quartz grains in sandstone may affect the crack path, we assumed them as continuum media for simplicity. However, it is an acceptable assumption for some applications in reservoir conditions, i.e., high pressure. Therefore, the maximum principal stress criterion has been used as the failure criterion, and the behavior of the samples after the failure is defined as energy adsorption.

The stress distribution and crack growth paths of the Font. #1 sample (Fig. 21(c)) with both strategies are presented in Fig. 22. As shown in Figs. 22(a) and 22(b) by arrows, crack initiation and propagation occur in the pores with the highest stress. It is observed that the crack initiation occurs more often near the large porosities, which have less cross-sectional area. Although the main crack’s growth is dominated in the first strategy (i.e., Fig. 22(c)), multiple cracks propagation are well illustrated in Fig. 22(c), which can show the modeling of the actual behavior of the sample in the second strategy.

The force-displacement curves of Font. #1 to 4 samples (Figs. 21(c) to 21(f)) in each strategy are shown in Fig. 23. As it can be seen, in the first strategy (Fig. 23(a)), the maximum force that the sandstone models can withstand before sudden failure is 536.6 kN for the Font. #1 sample, which is decreased to 245.1 kN for the Font. #4 sample. On the other hand, the maximum reaction force simulated by the second strategy (Fig. 23(b)) is much lower than the first strategy’s values and equals 105.3 and 146.5 kN for the first and the fourth sample, respectively.

Theoretically, reducing image resolution leads to more smoothed pores [ 62] and, in our case, reducing the porosity, making the sample much stronger against crack propagation. However, these results state that in the first strategy, the sample strength decreases by decreasing the resolution, whereas in the second strategy, this is an increasing trend and the sample strength increases by decreasing porosity. It is evident that results by the first strategy do not make sense in reality, because of the limitations we mentioned in Section 3.2. Also, according to Section 3.3.2, in the first strategy, the excited areas around the pores (i.e., the arrows in Fig. 22(a)) cannot create new cracks due to the influence of the primary crack on the crack growth path, and so the sample strength estimated is higher than the actual value. By decreasing the resolution (Figs. 21(d)-21(f)), the error rate is reduced, and the results are closer to the actual results.

As a result, we suggest implementing both strategies simultaneously and call it a hybrid strategy to improve the simulation results and use the advantages of both XFEM strategies. We demonstrate the hybrid strategy on another 2D image of the Fontainebleau sandstone (Fig. 24(a)). Similar to the previous image, a threshold value is used to segment the grayscale image (Fig. 24(b)). Then, both strategies are used simultaneously and in addition to sample partitioning, reduced Young’s modulus are assigned to the pores. As mentioned in Section 3.3.1, the first strategy can be a better strategy to estimate the crack deflection. On the other hand, according to Section 3.3.2, the reaction force and the multiple cracks growth paths can be simulated confidently by the second strategy. Therefore, the improved result is expected using the hybrid strategy. The results of this strategy are illustrated in Figs. 24(c) and 24(d). According to the stress distribution presented in Fig. 24(c), the sharp edges of the pores are susceptible to crack propagation. Also, the growth of multiple cracks in this sample (shown in Fig. 24(d)) represents the accuracy of this hybrid strategy.

Also, the force-displacement curve of this sample is shown in Fig. 25. According to this figure, the hybrid strategy’s maximum reaction force equals 132.6 kN, which is almost similar to the second strategy in Fig. 23(b).

5 Conclusions

In this paper, the modeling of the crack propagation in porous media is done using the XFEM. Two different strategies are used to model crack propagation in porous media. The first strategy is assigning loose material properties to the pores, while in the second strategy, multiple partitions are applied and enriched individually. The reliability of each XFEM strategy is evaluated both analytically and experimentally to find out the advantages and limitations of each strategy. Based on our investigations, the following results can be concluded.

1) Both strategies can identically predict the stress distribution in the sample. However, in the first strategy, only the initial crack propagation is simulated, and multiple cracks growth could be modeled in the second strategy only.

2) By increasing the number of pores in front of the crack tip, the stress intensity factor is theoretically increased, similar to both strategies. However, only the second strategy can predict a valid reaction force.

3) The results of the crack deviation to pore in both strategies are in good agreement with the analytical method. However, it seems the first strategy can estimate the crack deflection better, and the accuracy increases with decreasing the pores’ sizes.

4) By comparing the laboratory sample results with the simulation results of each strategy, it was found that the second strategy is a more accurate strategy in predicting the preferred crack growth path and determining the maximum strength of the sample.

5) In the porous media, applying several different partitions to the sample (second strategy) is a better and more accurate strategy than assigning reduced Young’s modulus to the pores (first strategy). However, we suggest the hybrid strategy, which combines both strategies to benefit from each strategy and achieve the best results.

References

[1] ChangS H, LeeC I, JeonS. Measurement of rock fracture toughness under modes I and II and mixed-mode conditions by using disc-type specimens. Engineering Geology, 2002, 66( 1−2): 79– 97

[2] HoekE, MartinC D. Fracture initiation and propagation in intact rock—A review. Journal of Rock Mechanics and Geotechnical Engineering, 2014, 6( 4): 287– 300

[3] LisjakA, KaifoshP, HeL, TatoneB S A, MahabadiO K, GrasselliG. A 2D, fully-coupled, hydro-mechanical, FDEM formulation for modelling fracturing processes in discontinuous, porous rock masses. Computers and Geotechnics, 2017, 81 : 1– 18

[4] JingL, HudsonJ A. Numerical methods in rock mechanics. International Journal of Rock Mechanics and Mining Sciences, 2002, 39( 4): 409– 427

[5] Cundall P A. A computer model for simulating progressive large scale movements in blocky rock systems. In: Proceedings of the International Symposium on Rock Mechanics. Nancy: International Society for Rock Mechanics, 1971

[6] LisjakA, GrasselliG. A review of discrete modeling techniques for fracturing processes in discontinuous rock masses. Journal of Rock Mechanics and Geotechnical Engineering, 2014, 6( 4): 301– 314

[7] WangS Y, SloanS W, ShengD C, YangS Q, TangC A. Numerical study of failure behaviour of pre-cracked rock specimens under conventional triaxial compression. International Journal of Solids and Structures, 2014, 51( 5): 1132– 1148

[8] KatoT, NishiokaT. Analysis of micro–macro material properties and mechanical effects of damaged material containing periodically distributed elliptical microcracks. International Journal of Fracture, 2005, 131( 3): 247– 266

[9] RezanezhadM, LajevardiS A, KarimpouliS. Effects of pore(s)-crack locations and arrangements on crack growth modeling in porous media. Theoretical and Applied Fracture Mechanics, 2020, 107 : 102529–

[10] Rodriguez-Florez N. Mechanics of cortical bone: Exploring the micro- and nano-scale. Dissertation for the Doctoral Degree. London: Imperial College London, 2015

[11] DuarteA P C, SilvaB A, SilvestreN, deBrito J, JúlioE. Mechanical characterization of rubberized concrete using an Image-Processing/XFEM coupled procedure. Composites Part B: Engineering, 2015, 78 : 214– 226

[12] DuarteA P C, SilvestreN, deBrito J, JúlioE. Numerical study of the compressive mechanical behaviour of rubberized concrete using the extended finite element method (XFEM). Composite Structures, 2017, 179 : 132– 145

[13] Supar K, Ahmad H. XFEM Modelling of Multi-holes Plate with Single-row and Staggered Holes Configurations. In: International Symposium on Civil and Environmental Engineering 2016 (ISCEE 2016). Wuhan: MATEC Web of Conferences, 2017

[14] RezanezhadM, LajevardiS A, KarimpouliS. Application of equivalent circle and ellipse for pore shape modeling in crack growth problem: A numerical investigation in microscale. Engineering Fracture Mechanics, 2021, 253 : 107882–

[15] HedjaziL, GuessasmaS, Della ValleG, BenseddiqN. How cracks propagate in a vitreous dense biopolymer material. Engineering Fracture Mechanics, 2011, 78( 6): 1328– 1340

[16] HedjaziL, MartinC L, GuessasmaS, Della ValleG, DendievelR. Application of the Discrete Element Method to crack propagation and crack branching in a vitreous dense biopolymer material. International Journal of Solids and Structures, 2012, 49( 13): 1893– 1899

[17] ChenM, WangH, JinH, PanX, JinZ. Effect of pores on crack propagation behavior for porous Si3N4 ceramics. Ceramics International, 2016, 42( 5): 5642– 5649

[18] Rodriguez-FlorezN, CarrieroA, ShefelbineS J. The use of XFEM to assess the influence of intra-cortical porosity on crack propagation. Computer Methods in Biomechanics and Biomedical Engineering, 2017, 20( 4): 385– 392

[19] RezanezhadM, LajevardiS A, KarimpouliS. Effects of pore-crack relative location on crack propagation in porous media using XFEM method. Theoretical and Applied Fracture Mechanics, 2019, 103 : 102241–

[20] RabczukT, BelytschkoT. Cracking particles: A simplified meshfree method for arbitrary evolving cracks. International Journal for Numerical Methods in Engineering, 2004, 61( 13): 2316– 2343

[21] RabczukT, BelytschkoT. A three-dimensional large deformation meshfree method for arbitrary evolving cracks. Computer Methods in Applied Mechanics and Engineering, 2007, 196( 29−30): 2777– 2799

[22] RabczukT, ZiG, BordasS, Nguyen-XuanH. A simple and robust three-dimensional cracking-particle method without enrichment. Computer Methods in Applied Mechanics and Engineering, 2010, 199( 37−40): 2437– 2455

[23] RenH, ZhuangX, RabczukT. Dual-horizon peridynamics: A stable solution to varying horizons. Computer Methods in Applied Mechanics and Engineering, 2017, 318 : 762– 782

[24] RenH L, ZhuangX Y, AnitescuC, RabczukT. An explicit phase field method for brittle dynamic fracture. Computers & Structures, 2019, 217 : 45– 56

[25] RenH, ZhuangX, RabczukT. A higher order nonlocal operator method for solving partial differential equations. Computer Methods in Applied Mechanics and Engineering, 2020, 367 : 113132–

[26] AreiasP, MsekhM A, RabczukT. Damage and fracture algorithm using the screened Poisson equation and local remeshing. Engineering Fracture Mechanics, 2016, 158 : 116– 143

[27] LehoucqR B, SillingS A. Force flux and the peridynamic stress tensor. Journal of the Mechanics and Physics of Solids, 2008, 56( 4): 1566– 1577

[28] KarimpouliS, TahmasebiP. A hierarchical sampling for capturing permeability trend in rock physics. Transport in Porous Media, 2017, 116( 3): 1057– 1072

[29] KarimpouliS, TahmasebiP, SaengerE H. Estimating 3D elastic moduli of rock from 2D thin-section images using differential effective medium theory. Geophysics, 2018, 83( 4): MR211– MR219

[30] HillerborgA, ModéerM, PeterssonP E. Analysis of crack formation and crack growth in concrete by means of fracture mechanics and finite elements. Cement and Concrete Research, 1976, 6( 6): 773– 781

[31] AsadpoureA, MohammadiS, VafaiA. Crack analysis in orthotropic media using the extended finite element method. Thin-walled Structures, 2006, 44( 9): 1031– 1038

[32] Mohammadi S. Extended Finite Element Method: For Fracture Analysis of Structures, Oxford: Blackwell Publishing Ltd, 2008

[33] SharafisafaM, NazemM. Application of the distinct element method and the extended finite element method in modelling cracks and coalescence in brittle materials. Computational Materials Science, 2014, 91 : 102– 121

[34] MoësN, BelytschkoT. Extended finite element method for cohesive crack growth. Engineering fracture mechanics, 2002, 69( 7): 813– 833

[35] GinerE, SukumarN, TarancónJ E, FuenmayorF J. An Abaqus implementation of the extended finite element method. Engineering Fracture Mechanics, 2009, 76( 3): 347– 368

[36] BelytschkoT, BlackT. Elastic crack growth in finite elements with minimal remeshing. International Journal for Numerical Methods in Engineering, 1999, 45( 5): 601– 620

[37] DolbowJ, MoësN, BelytschkoT. An extended finite element method for modeling crack growth with frictional contact. Computer Methods in Applied Mechanics and Engineering, 2001, 190( 51−52): 6825– 6846

[38] LiL, WangM Y, WeiP. XFEM schemes for level set based structural optimization. Frontiers of Mechanical Engineering, 2012, 7( 4): 335– 356

[39] MoësN, DolbowJ, BelytschkoT. A finite element method for crack growth without remeshing. International Journal for Numerical Methods in Engineering, 1999, 46( 1): 131– 150

[40] MoësN, CloirecM, CartraudP, RemacleJ F. A computational approach to handle complex microstructure geometries. Computer Methods in Applied Mechanics and Engineering, 2003, 192( 28−30): 3163– 3177

[41] ZhuQ Z. On enrichment functions in the extended finite element method. International Journal for Numerical Methods in Engineering, 2012, 91( 2): 186– 217

[42] AgathosK, ChatziE, BordasS P A. Multiple crack detection in 3D using a stable XFEM and global optimization. Computational Mechanics, 2018, 62( 4): 835– 852

[43] Sih G C. Methods of Analysis and Solution of Crack Problems. Leyden: Noordhoff International Publishing, 1973

[44] Anderson T L. Fracture Mechanics: Fundamentals and Applications. 3rd ed. Boca Raton: Taylor and Francis, 2005

[45] ArshadnejadS. Analysis of the first cracks generating between two holes under incremental static loading with an innovation method by numerical modelling. Mathematics and Computer Science, 2017, 2( 6): 120– 129

[46] ZhangZ. An empirical relation between mode I fracture toughness and the tensile strength of rock. International Journal of Rock Mechanics and Mining Sciences, 2002, 39( 3): 401– 406

[47] BažantZ P, KazemiM T. Size effect in fracture of ceramics and its use to determine fracture energy and effective process zone length. Journal of the American Ceramic Society, 1990, 73( 7): 1841– 1853

[48] MarshallG P, WilliamsJ G, TurnerC E. Fracture toughness and absorbed energy measurements in impact tests on brittle materials G. Journal of Materials Science, 1973, 8( 7): 949– 956

[49] Nasaj MoghaddamH, KeyhaniA, AghayanI. Modelling of crack propagation in layered structures using extended finite element method. Civil Engineering Journal, 2016, 2( 5): 180– 188

[50] ZhangC, CaoP, CaoY, LiJ. Using finite element software to simulation fracture behavior of three-point bending beam with initial crack. Journal of Software, 2013, 8( 5): 1145– 1150

[51] AbdellahM Y. Delamination modeling of double cantilever beam of unidirectional composite laminates. Journal of Failure Analysis and Prevention, 2017, 17( 5): 1011– 1018

[52] GrigoriuM, SaifM T A, El BorgiS, IngraffeaA R. Mixed mode fracture initiation and trajectory prediction under random stresses. International Journal of Fracture, 1990, 45( 1): 19– 34

[53] S. Moaveni, Finite Element Analysis: Theory and Application with ANSYS. Hoboken: Prentice Hall, 1999

[54] TroyaniN, PérezA, BaízP. Effect of finite element mesh orientation on solution accuracy for torsional problems. Finite Elements in Analysis and Design, 2005, 41( 14): 1377– 1383

[55] Logan D L. A First Course in the Finite Element Method. 4th ed. Toronto: Nelson, 2007

[56] SongJ, BelytschkoT. Cracking node method for dynamic fracture with finite elements. International Journal for Numerical Methods in Engineering, 2009, 77( 3): 360– 385

[57] LinderC, ArmeroF. Finite elements with embedded branching. Finite Elements in Analysis and Design, 2009, 45( 4): 280– 293

[58] LiX, KonietzkyH. Simulation of time-dependent crack growth in brittle rocks under constant loading conditions. Engineering Fracture Mechanics, 2014, 119 : 53– 65

[59] AndräH, CombaretN, DvorkinJ, GlattE, HanJ, KabelM, KeehmY, KrzikallaF, LeeM, MadonnaC, MarshM, MukerjiT, SaengerE H, SainR, SaxenaN, RickerS, WiegmannA, ZhanX. Digital rock physics benchmarks—Part I: Imaging and segmentation. Computers & Geosciences, 2013, 50 : 25– 32

[60] MadonnaC, QuintalB, FrehnerM, AlmqvistB S G, TisatoN, PistoneM, MaroneF, SaengerE H. Synchrotron-based X-ray tomographic microscopy for rock physics investigations. Geophysics, 2013, 78( 1): D53– D64

[61] Huang J Q, Huang Q A, Qin M, Dong W J, Chen X W. Experimental study on the dielectrostriction of SiO₂ with a micro-fabricated cantilever. In: IEEE Sensors 2009 Conference. Christchurch: IEEE, 2009

[62] KarimpouliS, TahmasebiP, SaengerE H. Coal cleat/fracture segmentation using convolutional neural networks. Natural Resources Research, 2020, 29( 3): 1675– 1685

RIGHTS & PERMISSIONS

Higher Education Press 2021.