A peridynamic model for damage and fracture in porous materials

: We introduce a peridynamic (PD) model for simulating damage and fracture in porous materials based on an Intermediate Homogenization (IH) approach. In this approach, instead of explicitly representing the detailed pore geometry, we use homogenization but maintain some information about the microstructure (porosity) in the model. Porosity is introduced in the model as initial peridynamic damage, implemented by stochastically pre-breaking peridynamic bonds to match the desired porosity value. We validate the model for elastodynamics with wave propagation in porous glasses, where we match observed wave propagation speeds and the apparent elastic moduli for various porosities. The model is then used to study the fracture behavior of Berea sandstone under three-point bending loading conditions. We validate the model for fracture problems using the case of failure in a sandstone sample with an off-center pre-notch under three-point bending. The IH-PD results agree very well with experiments: we obtain different failure modes depending on the length of the off-center pre-notch. When the pre-notch is short, most damage (and the subsequent


Introduction
Damage and fracture in natural or man-made porous materials (bone, wood, rock, sandstone, or metal foams, filters, concrete, composites, ceramics) present particular challenges: damage can be widespread or finely localized, as micro-cracks "jump" over pores and link together into a number of macro-cracks leading to final failure or they stay separated leading to microcrushing and eventual large scale fragmentation. The underlying physical principles governing their mechanical behavior, but not their failure behavior, is known as "poromechanics". [1][2][3] Poroelastic media are elastic materials with pores/voids of arbitrary shapes and sizes, usually randomly distributed and oriented through the material. Many porous materials can be considered as elastic quasi-brittle materials, as least locally. Their global failure response can, sometimes, appear as ductile failure, as locally brittle damage spreads over pores, flaws and other kinds of discontinuities. The geometry and distribution of pores have a close relationship to their initiation and crack propagation. Although many continuum-based and discontinuous (e.g. distinct element method) methods based on classically homogenized poromechanical approaches have been developed and applied to model porous materials fracture, it remains an open question whether they can correctly capture their fracture behavior. 12 A compromise between having a low computational cost (on par with that of homogenized models) and being able to reproduce the effects the porous microstructure has on the initiation and propagation of fracture and damage is sought in the present contribution. For this purpose, we introduce an intermediately-homogenized peridynamic (IH-PD) model and validate it against experimental results for wave propagation in porous media of different porosities. We then test the new model on Berea sandstone fracture from a three-point bending setup with an off-center prenotch. Experiments show that the sample fails from the pre-notch if the notch is sufficiently long, but if the notch is short, the porous stone fails from its center (the point of maximum tensile stress of an un-notched specimen). We compare the results from the IH-PD model with those from experiments, with other methods from the literature, and with a fully-homogenized peridynamic model. We comment on the results obtained and discuss the importance of preserving the effect of local heterogeneities on the damage initiation, distribution, and crack propagation. This can be accomplished with the new IH-PD model, and at a similar cost with a fully homogenized peridynamic model. Fully homogenized models (classical or peridynamics-based), however, lose important heterogeneity information that controls crack initiation and propagation and because of that will miss, in a significant way, the correct failure behavior of porous materials. Therein lie the benefits for the proposed modeling approach we present in this paper.

Literature review
Peridynamics, a new nonlocal continuum theory proposed by Stewart Silling, 15,16 eliminates spatial derivatives from the classical formulation of continuum mechanics, which makes it a consistent mathematical model for problems with discontinuities in the displacement field. The model is particularly well suited for dealing with cracks and damage in solid mechanics especially in situations where the crack path is not known in advance. These features have enabled peridynamics to be successfully applied to diverse engineering problems that involve fracture and damage evolution. [17][18][19][20][21][22][23] Several recent studies have applied peridynamics in porous materials-related problems: Katiyar et al. 24 and later Jabakhanji and Mohtar 25 developed a peridynamic approach to model flow in porous media. Peridynamics was also used to address hydraulic fracturing and reproduce the deformation induced by the fluid flow in a porous medium. [26][27][28] In the present study, we focus on the initiation and propagation of cracks in dry porous rock. Recently, Zhou et al. investigated how complex fracturing patterns initiate and propagate from a single flaw embedded in rock-like materials under compression using bond-based peridynamics. 29 With this homogeneous PD model, Zhou and co-authors simulated crack propagation in rocks with pre-existing flaws (e.g. [30][31][32] ). These PD models, when applied to failure of materials with relatively high porosity, would have the same issues as the fully-homogenized PD model discussed in section 5. Lee et al. applied bond-based peridynamics to understand crack coalescing processes in rock materials. 33 In these studies of rock fracture, the actual rock materials studied are dense materials with small porosity, such as marble 29 (porosity: 2% 34 ), and Longtan sandstone 14 (porosity: 1.38% 35 ). When the porosity is not negligible, such as in Berea sandstone (porosity range: 10.2~22.2% [36][37][38], damage initiation and crack propagation are strongly influenced by it. As we shall see in section 5, this aspect limits the applicability of such homogenized models.

Peridynamic model for porous materials
In this section, we first briefly review the peridynamic theory for elastic brittle materials.
Then, we describe the fully homogenized PD model (FH-PD) and introduce the new Intermediately Homogenized PeriDynamic (IH-PD) model for porous materials.

Brief review of peridynamic theory for elastic brittle materials
The peridynamic model is a framework for continuum mechanics based on the idea that pairs of material points exert forces on each other across a finite distance. This concept allows for the natural evolution (initiation, propagation, and interaction) of damage and cracks and can be viewed as an effective treatment of material length-scale induced by, for example, the material microstructure. The peridynamic equations of motion for the bond-based model are given as: 15 where f is the pairwise force function in the peridynamic bond that connects point � to x, u is the displacement vector field, ρ is the density and b(x,t) is the body force. The integral is defined over a region Hx called the "horizon region", or simply the "horizon". The horizon is the compactsupported domain of the pairwise force function around a point x (see Fig. 1). The horizon region is taken here to be a circular disk (sphere) of radius δ . We refer to δ also as the "horizon", and from the context there should be no confusion whether we refer to the region or its radius. A micro-elastic material is defined as one for which the pairwise force derives from a potential: where = � − is the relative position in the reference configuration and = (�, ) − ( , ) is the relative displacement between and �. A micropotential that leads to a linear microelastic material is given by: where = ‖ ‖, and = is the relative elongation of a bond, or bond strain. The function ) (ξ c is called the micro-modulus and has the meaning of bond's elastic stiffness. The pairwise force corresponding to the micropotential given above has the following form: Following the same procedure performed to calculate the micro-modulus functions in 1D (see 40 ), one obtains the conical micro-modulus function in 2D, plane stress conditions: 41 or, assuming a constant micro-modulus function over the horizon region, in 2D, plane stress conditions: where E is Young's modulus of the material. In the IH-PD model of porous material, E is the elastic modulus of the constituent material of the porous medium and this is discussed in Section 3.3. The material model in Eq. (4) is equivalent to the kernel function using = 1 for peridynamic kernel (�, )/|� − | by Chen et al. 42 They constructed a peridynamic kernel ( = 2) based on physical principles for dynamic elasticity and showed that, when the one-point Gauss quadrature is used for discretization, the = 2 model is the only one whose convergence to the classical solution does not depend on the fineness of the discretization grid. Models with = 0 or = 1 also converge, in the limit of horizon going to zero and ratio of horizon to grid spacing going to infinity, to the classical solution for problem with sufficient smoothness. No significant differences on crack patterns were observed between models with = 1 and = 2. 43 In this work we use the model with = 1.
Failure is introduced in peridynamics by considering that the peridynamic bonds break when they are deformed beyond a critical value, called the critical relative elongation or critical bond strain, s0, computed based on the material's fracture energy. In 2D, the energy per unit fracture length for complete separation of the two halves of the body is the fracture energy G0. Equating it to the work done in a PD material to accomplish the separation of the body into two halves gives: Substituting the micro-modulus functions from Eqs. (5) and (6) into Eq. (7), respectively, s0 is obtained for the conical micromodulus function (under plane stress conditions) as: δ π E G s 9 and for the constant micromodulus function as: For the IH-PD model of poroelastic materials, 0 G is the fracture energy of the constituent material of the porous medium (see details in Section 5).

Numerical discretization
In this section we describe the numerical implementation details, including both dynamic and static solutions. In principle, the peridynamic equations can be discretized using the finite element method, or any other method appropriate to compute solutions to integro-differential equations (or integral equations for the static case). These approaches, however, soon hit wellknown obstacles and difficulties for problems with evolving topologies, like those in which cracks initiate, grow, and interact with each other leading to damage, full fracture and/or fragmentation.
Instead, meshfree-types discretizations are preferred for peridynamic simulations of material failure and damage. The discretization proposed by Silling and Askari 44 uses the mid-point integration scheme (equivalent to one-point Gauss quadrature) for the domain integral. Numerical simulations are performed using the following discretized equation: where Fam(i) is the family of nodes j with their area (volume in 3D) covered, fully or partially, by the horizon region of node i, is the bond length between nodes i and j, is the relative elongation for the bond connecting nodes i and j, and is the area of node j estimated to be covered by the horizon of node i.
Note that node j may not be fully contained within the horizon of node i, so a "partial volume" integration, first introduced by Bobaru et al 45 and also shown in their following work, 46 is used here to improve the accuracy of the mid-point quadrature scheme. The main advantage of this algorithm compared with one that simply checks whether a node is inside or outside the horizon region is that, as the grid density increases (for a fixed horizon value), the numerical convergence (in terms of strain energy density, for example) is monotonic. 45 For a fixed horizon, the ratio = /Δ describes how accurate the numerical quadrature for the integral in Eq.
(1) will be. We call this ratio "the horizon factor". In the convergence study shown in Section 4, we study both m-convergence and δ-convergence 40 for an elastic wave propagation problem. We recall that in m-convergence we consider the horizon δ fixed and take → ∞. The numerical PD approximation will converge in this case to the exact nonlocal PD solution for the given δ. In the case of δ-convergence, the horizon → 0 while m is fixed or increases with decreasing δ. For δ-convergence and in problems with no singularities, the numerical PD solutions are expected to converge to the classical local solution (as m increases).
Both dynamic (see Section 4) and static (see Section 5) simulations are performed in this work. In the dynamic simulations for elastic wave propagation in a porous glass (see Section 4), we apply the Velocity-Verlet method with a time interval of 0.1 µs. For the quasi-static fracture tests in Section 5, the energy minimization method 47,48 is used, and the nonlinear conjugate gradient (CG) method with secant line search is adopted to minimize the strain energy of the system. For all static simulations in this paper, the nonlinear CG method is used with a convergence tolerance defined by: < 10 −6 , in which , and −1 are the total strain energy at the current (k-th) and previous (k-th-1) CG iterations.
(1), we obtain the nonlinear system of the discretized equations for quasi-static conditions: where bi is the body force at node i.

An intermediately-homogenized peridynamic model for poroelastic materials
In the fully-homogenized peridynamic (FH-PD) models, heterogeneous materials are understood as material models locally homogenized in terms of their elastic properties, density, and fracture energy. 19 The properties used in the FH-PD models can be calculated from the constitutive models of porous materials, such as Voigt or Reuss models with known porosity and properties of the constituent material, or the apparent parameters from direct measurements.
Porous materials, when the pores sizes are much smaller than the sample sizes, can also be viewed as locally homogeneous. Excellent simulation results are produced by the FH-PD models for rock with low porosity. When the rock porosity is not negligible, such as in Berea sandstone, local heterogeneities can highly influence damage initiation and crack propagation, which limits the applicability of fully-homogenized models for these cases. In this section, we introduce the intermediately-homogenized peridynamic (IH-PD) model for porous materials. The reader is referred to 49 for more details on using the IH-PD model for multi-phase materials.
Here we introduce a peridynamic model for porous materials based on an Intermediate Homogenization (IH) approach to simulate fracture and damage in such materials. The porosity is represented by peridynamic pre-damage, with mechanical bonds connected to a peridynamic node being pre-broken, stochastically, to achieve the desired porosity. The "intermediate homogenization" approach refers to the fact that we do not represent the explicit geometry of the actual pores in the material, and neither do we fully homogenize the porous medium. With the IH method we aim to maintain sufficient information about the porosity to allows us to more accurately compute the failure behavior of porous materials compared with the FH-PD. At the same time, the IH-PD model will be significantly more efficient computationally compared with a model that uses an explicit geometric representation of the microstructural individual pores, while hopefully, still being able to compute the macro-scale failure behavior with accuracy.
In the IH-PD model of porous materials, pores are treated as pre-existing material damage. In peridynamics, the damage index is computed as the ratio between the number of broken (or failed) bonds (Nf) and the total number of bonds (N) originally associated with a node ( ( , ) = ). To mimic the presence of pores, we randomly break a number of bonds at each node (see Fig. 2). We calibrate the number of broken bonds to the material's porosity. This procedure creates a "predamage index", related to material's porosity and computed like the damage index above but with Np, the number of pre-broken bonds at a node, replacing Nf. When the porosity reaches the critical porosity (the porosity beyond which the rock can exist only as a suspension 6 ), all bonds associated with that point should be broken ( = ), meaning that the pre-damage index is unity.
For zero porosity, no pre-broken bonds are introduced (pre-damage index is zero at all points).
To perform the calibration for the number of pre-broken bonds, we implement the initial damage representative of the pores by adopting the ideas used in the "concentration-dependent damage" (CDD) model originally proposed by Chen and Bobaru 18 and used for modeling of damage induced by corrosion processes [50][51][52][53] . Here, we assume that the pre-damage index at a point depends linearly on its porosity when zero-porosity material surrounds the point: iii. Go to next bond in the family of bonds of node x.
This procedure is applied for all nodes in the material. Here we assume that the pore shapes and distribution are "isotropic", therefore the peridynamic pre-damage representation of this case uses the same uniform distribution for bond-breaking independent of the node location or the bond direction. Special location-and/or orientation-dependent pre-damage can easily be introduced to mimic anisotropy in porous materials.
With the above algorithm, each PD bond (connecting material points and �) goes through the procedure twice. If porosities at and � are ( ) and (�), respectively, the chance for the . For materials with uniform porosity P, the change for any bond to stay intact is: (1 − / ) 2 . The pre-damage index, the ratio between the number of pre-broken bonds and the total number of bonds at a node, in this case, converges to: in the limit of the horizon factor (the ratio between the horizon size and grid spacing) going to infinity (also known as the m-convergence, see Section 3.2). Notice the nonlinear dependency, for a material with uniform porosity, between the pre-damage index and the given porosity.
In In the next section, we use the IH-PD model to simulate elastic wave propagation in porous glasses. By computationally measuring the wave speeds corresponding to different porosities, we calculate the apparent elastic moduli, and compare them with those measured in experiments.

Model verification and validation for elasto-dynamic problems
We verify (by performing convergence studies) and validate the IH-PD model for an elastodynamic problem: elastic wave propagation in porous glasses. We generate elastic waves by suddenly applying a force pulse over a small region in the sample and monitor their propagation through the sample. Experimental data regarding the apparent elastic modulus of these porous glasses in terms of varying degrees of porosity exists. 54 We test if the model can predict the experimentally measured relationship between elastic wave speed and material porosity. We arbitrarily choose a sample size of 1 m × 1 m. A suddenly applied load of 1 MPa in the vertical-up direction over the region (100 mm in length) at the bottom side is kept constant for 5 µs (see Fig. 3). After that, the load is suddenly removed. The loading generates elastic waves propagating through the sample.
Although the horizon-size-dependence of the loading region leads to wave patterns that depend on the horizon size, the wave speed does not change much with the horizon size and should converge to the measured value when the horizon size goes to zero. 40 Since our goal is to calculate the wave speeds for different porosities, the specific initial loading conditions (loading range, loading magnitude) are not important.
To compute the wave speed from the PD model results, we track the location of the pressure wave's peak over a time period that starts shortly after the pulse is applied and stops before waves return from the sample's boundaries. More specifically, we compute an average velocity by tracking the displacement of the crossing point between line x = 0 (shown in Fig. 4) and the crest of the front wave, over the time period 100 μs -200 μs from the application of the pulse. For this verification phase we are interested only in the elastic wave propagation, therefore we apply a "nofail" condition for all bonds in the model. Obviously, the pre-damage corresponding to different porosities is present in these calculations, but no "new damage" is allowed.   Fig. 4, the horizon size is 40 mm, and the horizon factor m=8.
One notable difference between the nonporous (P = 0) and the porous case (e.g., P = 0.5) that can be observed from these results, is that in the porous sample stress waves become less coherent, being locally dispersed ("noisier" velocity maps) by the more detailed representation (with the IH-   We now study m-convergence and -convergence to see the influence of the horizon factor and nonlocal horizon size ( ) on the wave speed at different porosities.  Fig. 6 shows the crest locations versus time, for different porosities. We use a linear fit for the vertical velocity data points within the time range mentioned above to extract the wave speed velocity shown in Fig. 7, for each porosity. As mentioned before, at high porosities, (0.7, 0.9), there is a relatively larger probability for peridynamic nodes to lose all of their connections with the surrounding nodes if m is relatively small, since the number of bonds is relatively small. This is the reason for which only the higher values of m are used in these cases (see Fig. 7b).
The results in Fig. 7 show that m = 4 is a safe choice for simulating elasto-dynamic problems in porous materials with the IH-PD model when the porosity is smaller than 0.5 (relative to a critical porosity of 1).
The apparent elastic modulus in glass samples with different porosities was experimentally measured. 54 To compare the PD results with experimental measurements, we compute the apparent modulus (EP) from the wave speeds numerically obtained in the IH-PD simulations above. Based on the theory of 2D (plane stress) wave propagation, 55 the relationship between longitudinal or pressure wave speed and the apparent elastic modulus (EP) is: where L C is the longitudinal waves speed. In the computation, we use 1/3 as the Poisson ratio, because the Poisson's ratio value for bond-based peridynamic model in 2D under plane stress conditions is 1/3. 56 The apparent modulus calculated from Eq. (13) by using the wave speed computed with the PD model matches well the experimental measurements (see Fig. 8).  With this model validation for elasto-dynamic problems, we next validate the IH-PD model for fracture using quasi-static fracture of a porous rock. We consider fracture induced by threepoint bending of Berea sandstone.

Peridynamic modeling of damage evolution in porous materials
We apply both the FH-PD and IH-PD models for quasi-static fracture of a porous rock sample. We find that only the IH-PD model delivers damage patterns and crack profiles similar to those seen in experiments in. 9

Description of experimental setup and results, and of the numerical model
The experimental setup shown in Fig. 10 Fig. 8 in Lin et al. 9 ) in mixed-mode. Fig. 7 in Lin et al. 9 also shows that for the long notch, some damage, measured via acoustic emission (AE) and electronic speckle pattern interferometry (ESPI), is recorded near the beam's center. In what follows, we investigate whether the observed damage sensitivity to notch length are reproduced by any of the two PD models described in this paper. We will compare our PD simulation results with the experimental results in terms of fracture paths, locations of microcracks, and peak loads. These comparisons are not possible through full sample failure simply because the experimental results shown 9 stop before that. In our simulations, the imposed displacement is 2μm between increments. We run 2,000 increments to ensure full splitting of the samples.
Berea sandstone is chosen for our study because of its fairly uniform grain size, 9 which makes it a good candidate for the isotropic porosity model, as well as for its brittle failure behavior.
Although Berea sandstone from different locations may have different porosities and thus different corresponding apparent moduli, we assume they all obey the relationship (see Section 4) = (1 − / ) 2 , with critical porosity PC and zero-porosity modulus E. The apparent fracture energy ( 0, ) has a similar relationship with porosity as the apparent modulus, thus Poisson ratio does not happen to match as well, can be performed using the state-based peridynamic formulation. 47 The state-based PD formulation eliminates the Poisson ratio restriction.
Because the Poisson ratio for the material considered here is close to 1/3 and the fact that our primary interest here is to observe the capabilities of a PD model in capturing the evolution of fracture and failure in poroelastic materials and to determine other factors that control crack growth in poroelastic materials, we use the bond-based peridynamic formulation. In our simulations, we use the conical micromodulus (see Eqs. (5) and (9)), δ = 2 mm and m = 4.

Results and discussion
We use the FH-PD and IH-PD models (with the quasi-static loading conditions specified in the previous section) to simulate the cases shown in Fig. 10, and to examine the effect the prenotch length has on the failure mode and crack patterns. In the IH-PD model, the effect of porosity is represented by inserting "pre-damage" in the material, while in the FH-PD model, porosity is indirectly represented by the apparent elastic and fracture properties. the long pre-notch leads to crack initiating from the off-center pre-notch tip, while in the short prenotch case, the crack initiates from the bottom-center of the beam.
Movies 6 and 7 show the damage evolution with and without pre-damage due to porosity, respectively, for the short pre-notched sample using the IH-PD model. Movies 8 and 9 show the damage evolution with and without pre-damage due to porosity, respectively, for the long prenotched sample using IH-PD model. In movies 6-9, the imposed displacement of the top-center of beam progresses from 0 to 2.6 mm. In these movies, the damage quantity shown is scaled to a 64 range between 0 and 1, and the same is done for the rest of such movies.
As displayed in movies 6-9, our PD simulations track the crack growth through full failure,   Because of the simplifying assumptions implicit in a full homogenization of a porous material (absence of microscale heterogeneities and defects), the high stress concentration at the tip of preexisting notch, whether long or short, causes fracture to initiate (Fig. 12). In reality, small scale heterogeneities (presence of pores) result in local stress concentrations at places other than the offcenter notch tip, and most likely in the regions near the bottom center of the beam where tensile bending stress is highest. This is what leads cracks to initiate from the beam's middle bottom rather than the pre-notch tip, in the short-notch case ( Fig. 11c and 11e).
c d e f Movies 10 and 11 show the damage evolution for the short and long pre-notched samples, respectively, using the FH-PD model with conical micromodulus. In these simulations the imposed displacement at the top-center of beam progresses from 0 to 1.9 mm. To make it easier to analyze the computed damage distribution in the material, we blank the material nodes with damage index of zero in Fig. 13. The snapshots shown in Fig. 13  13a and 13b) match very well the experiments (see Fig. 3 in Lin et al. 9 ). Note that the initiation of fracture in the specimen with the short notch is the same as that in the specimen without a notch (similar distribution of damage around the final crack path). Figs. 7 and 8 in Lin et al. 9 show that, experimentally, while for the long notch case the specimen fails from the notch tip, some damage is also recorded around the bottom center of the beam. As seen in Fig. 13b, the IH-PD model captures this micro-damage at the beam center for the long notch case, as observed in experiments, while the FH-PD model (Fig. 13d) does not. al. 9 , with the data produced by the IH-PD model for the long notch case. The available experimental results (red diamonds in Fig. 14) only show the early stages of the fracture process.

a b
The peridynamic results (from the IH-PD model) match well with the acoustic emission data 9 , showing damage accumulating at the right corner of the notch tip, while the numerical results from Lin et al. 9 , show crack initiation from the left corner of the notch. The damage evolution for the long notch case solved with the IH-PD model (seen from movie 13) shows distributed damage around the beam bottom-center occurring earlier than the damage that eventually leads to full failure, starting at the pre-notch tip.
We note that the location of long notch shown in Figs. 7 and 8 in Lin et al. 9 is incorrect. Based on Table 2 in Lin et al., 9 the notch is located at two-thirds of the half-span length (about 49 mm from the center, see Fig. 10); however, the notch is shown to be located close to 40 mm in Figs. 7 and 8 in Lin et al.. 9 We use the parameter in Table 2   color squares, and using the legend in Fig. 13), and the numerical model in Lin et al. 9 (black dots). The off-center black rectangle shows the notch geometry and location.
In Table 1  Compared with the numerical model presented in Lin et al., 9 the IH-PD model introduced here has following advantages: i.
The IH-PD model presented here uses only an elastic-and-brittle-failure model, matching the sandstone micro-scale behavior. The actual sandstone heterogeneities/porosity (represented in our model via pre-damage) lead to an effective, macroscale behavior that displays softening. In contrast, Lin et al.'s 9 use a softening contact bond model, in which the normal bond strength reduces linearly after its peak, to mimic the observed macroscale behavior.
ii. No trial-and-error calibration for material constants is necessary for the IH-PD model. In the IH-PD model, macroscopic measurable properties (e.g. porosity, etc.) are input parameters. On the other hand, the constitutive model in the contact bond model 9 is defined by many microscopic material constants. Through the interaction of many discrete element particles (with microscopic material constants), macroscopic properties emerge. Hence, trial and error steps (or calibration processes) are required to achieve appropriate macroscopic properties similar to those of a specific rock.

Conclusions
In this paper we presented an Intermediately-Homogenized Peridynamic (IH-PD) model for simulating the elastic and failure behavior of porous materials. In this model, porosity was represented by an initial peridynamic damage ("pre-damage"), introduced by breaking bonds, stochastically, to achieve a given porosity. The peridynamic micromodulus was computed using the elastic modulus of the porous medium constituent material. We validated the model for elastic wave propagation using experimental data for wave propagation speeds and the apparent elastic moduli in porous glass with different porosities.
While simple homogenization methods work well for linear problems (e.g., elastic response), for nonlinear problems (e.g., in damage and failure), where dissipation and weakest links play a role in determining the material behavior, such methods may fail to reproduce experimental observations. To answer the question "how much homogenization is too much?" when modeling fracture processes in porous materials, we compared two peridynamic models: one that homogenized the material to a greater extent (the Fully-Homogenized Peridynamic model, FH-PD) and the newly introduced IH-PD model. We studied a quasi-static crack growth problem in a brittle porous rock to understand the difference between the models' responses. We found that the fully-homogenized model (FH-PD) failed to capture the experimentally observed fracture behavior in Berea sandstones sample under three-point bending loading. For this problem, experiments show that fracture patterns are controlled by the size of the off-center pre-notch. The IH-PD model, however, reproduced perfectly the observed crack growth behavior and its dependency on the length of the pre-notch.
We conclude that for problems in which the microstructure has a large effect on the failure behavior, a fully homogenized strategy will not work to correctly capture the fracture and damage evolution. For a predictive model, some of the details of the microstructure and its role in controlling crack growth are needed. The new IH-PD model can predict the correct failure evolution in a porous material without requiring the explicit description of the microstructure geometry. This happens because some essential features (porosity represented via peridynamic pre-damage) of the material microstructure are incorporated in the model. It is this extra microstructure information that allows the computed damage to initiate and evolve in a way similar to that in the real porous medium.