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 represent ing 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 v alidate the model for elast odynamic s with wave propagation in porous glass es, where we match observed wave propagation speeds and the apparent elastic modul i 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 s ample 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 2 failure) happens in the center of the beam (largest tensile stresses in an un -notched sample ). In this case, material porosity make s the fracture behavior “insensitive” to the presence of the off -center pre-notch. When the pre -notch is long, damage starts from near its tip , and progress es in mixed mode towards the loading point. T he model capture s the location of crack initiation to be the right - corner of the notch, exactly as shown by the experimental acoust ic emission results. A fully - homogenized model for the porous sandstone fails to reproduce the failure mode sensitivity on the pre-notch length. The results show the importance of using the IH -PD model when trying to capture the important effect s local heterogeneities (like pores) hav e on damage initiation and growth in porous materials . Keywords: Porous materials; Porosity ; Peridynamics; Fracture ; Rock; Damage. 1. Introduction Damage and fracture in natural or man -made p orous 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 t o final failure or they stay separated leading to micro - crushing and eventual large scale fragmentation. The underlying physical principles governing their mechanical behavior , but not their failure behavior, is known as “poromechanics”.1-3 Poroelastic media are elastic materials with pores /voids of arbitrary shapes and sizes, usually randomly distributed and oriented through the m aterial . Many p orous 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 ha ve a close relationship to the ir 3 mechanical properties , such as mechanical strength, Young's modulus, permeability and sonic velocity.4-7 Many constitutive models have been proposed to investigate and predict elastic properties of porous materials , for ins tance: Voigt, Reuss, Wood, Hashin- Shtrikman (HS), Gassmann, Modified Voigt and Modified HS models.5 These poromechanical models capture the dependence of apparent, or effective, elastic properties on porosity; they do not, however, say anything about the dependence of “effective” fracture and damage properties on pores .8 Our goal is to introduce a model which can directly capture the influence of porosity on the fracture and damage evolution, without the need of the explicit ly modeling the porous geometry. Damage initiation , progression, and final failure in rock, a porous material, is an important topic in rock engineering , hydraulic fracturing, and mining sciences . Substantial efforts, experimental , theoretical , and computational, have been dedicated to stud ying damage evolution in rock under various loadings .9-11 Some recent computational methods have been helpful in predict ing the initiation and propagation of cracks in porous materials , especially rocks .12 Two main types of computational models have been used to investigate damage and fracture in porous materials: models that use a detailed description of the microstructure geometry, and homogenized approach es. The explicit microstructure models involve either digitized porous microstructures from experimental images , or generated by researchers to mimic actual ones . These models have sever e limitations in terms of the size of the sample that can be modeled, due to the high computational cost. When the sample size is much larger than the pore size, it is computationally unfeasible to consider the exact geometry and location for every single pore.2 Consequently, averaged, or “homogenized” models have been used as an alternative for simulating porous materials ’ behavior under mechanical loading conditions .13, 14 These homogenization approach es, however, may not be able to captur e the strong influence local heterogeneities have on damage 4 initiation and crack propagation. Although many continuum -based and discontinuous (e.g. distinct element method) methods based on classically homogenized poromechanical approaches ha ve been developed and applied to model porous materials fracture, it remains an open question whether the y can correctly capture the ir 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 medi a of different porosit ies. We then test the new model on Berea sandstone fracture from a three- point bending setup with an off -center pre - notch. Experiments show that the sample fails from the pre -notch if the notch is suf ficiently 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 literatu re, and with a fully -homogenize d 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 peridynamic s-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. 5 2. 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. The se features have enabled peridynamics to be successfully appli ed to diverse engineering problems that involve fracture and damage evolution. 17-23 Several recent studies have applied peridy namics in porous materials -related problems : Katiyar et al.24 and later Jabakhanji and Mohtar25 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-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 us ing 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-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 process es in rock materials.33 In these studies of rock fracture, the actual rock materials studied are dense materials with small porosity, such as marble29 (porosity : 2%34), and Lon gtan sandstone14 (porosity: 1.38%35). When the porosity is not negligible, such as in Berea sandstone 6 (porosit y range: 10.2~22.2%36-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 homogeniz ed models. 3. 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 H omogenized P eriDynamic (IH-PD) model for porous materials. 3.1. 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 micros tructure. The peridynamic equations of motion for the bond- based model are given as:15 𝜌𝜌𝒖𝒖̈(𝒙𝒙,𝑡𝑡)=∫𝒇𝒇�𝒙𝒙�−𝒙𝒙,𝒖𝒖(𝒙𝒙�,𝑡𝑡)−𝒖𝒖(𝒙𝒙,𝑡𝑡)�d𝑉𝑉𝒙𝒙� 𝐻𝐻𝒙𝒙+𝒃𝒃(𝒙𝒙,𝑡𝑡) , (1) 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 H x called the “ horizon region”, or simply the “horizon”. The horizon is the compact - supported 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. 7 Fig. 1. Each point x interacts directly with any point x� in the horizon region delimitated by the red circle . A micro -elastic material is defined as one for which the pairwise force derives from a potential: 𝒇𝒇(𝝃𝝃,𝜼𝜼)=𝜕𝜕𝜕𝜕(𝝃𝝃,𝜼𝜼) 𝜕𝜕𝜼𝜼 , (2) 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: ξξω2)(21sc= , (3) 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: δx Hx 8  >≤++ = δξδξξ )(),( 0scηξηξ ξηf (4) Following the same procedure performed to calculate the micro -modulus functions in 1D (see40), one obtains the conical micro- modulus function in 2D, plane stress conditions:41  −= −= δξ πδδξξ 1361 )(3 1Ec c , (5) or, assuming a constant micro -modulus function over the horizon region, in 2D, plane stress conditions: 39() ,Ecξ πδ= (6) 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 li mit of horizon going to zero and ratio of horizon to grid spacing going to infinity, 9 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 G 0. Equating it to the work done in a PD material to accomplish the separation of the body into two halves gives: ∫∫∫−  =δδξ ξθξξξ 0)( cos 02 0 01 2)(2 dzddscG zz . (7) Substituting the micro -modulus functions from Eqs. (5) and (6) into Eq. (7), respectively, s 0 is obtained for the conical micromodulus function (under plane stress conditions) as: δπ EGs 950 0= , (8) and for the constant micromodulus function as: 0 04 9GsEπ δ= . (9) 10 For the IH -PD model of poroelastic materials, 0G is the fracture energy of the constituent material of the porous medium (see detail s in Section 5) . 3.2. 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 equa tions (or integra l equations for the static case). These approaches, however, soon hit well - known obstacles and difficulties for problems with evolving topologies, like those in which cracks initiate , grow, and interact with each other lead ing to damage, f ull 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: �𝒇𝒇�𝒙𝒙�−𝒙𝒙𝑖𝑖,𝒖𝒖(𝒙𝒙�,𝑡𝑡)−𝒖𝒖(𝒙𝒙𝑖𝑖,𝑡𝑡)�d𝑉𝑉𝒙𝒙� 𝐻𝐻𝒙𝒙 ≃�𝑐𝑐(𝜉𝜉𝑖𝑖𝑖𝑖)𝑠𝑠𝑖𝑖𝑖𝑖𝑉𝑉𝑖𝑖𝑖𝑖 𝑖𝑖∈Fam (𝑖𝑖) (10) 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 al45 and also shown in their following work,46 is used 11 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 e xample) 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 S ection 4, w e 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 method47, 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 m ethod is used with a convergence tolerance defined by: |𝑊𝑊𝑘𝑘−𝑊𝑊𝑘𝑘−1| 𝑊𝑊𝑘𝑘−1<10−6, in which 𝑊𝑊𝑘𝑘, and 𝑊𝑊𝑘𝑘−1 are the total strain energy at the current ( k-th) and previous ( k-th-1) CG iterations. Substituting Eq. (10) into Eq. (1) , we obtain the nonlinear system of the discretized equation s for quasi -static conditions : 12 ∑ 𝑐𝑐(𝜉𝜉𝑖𝑖𝑖𝑖)𝑠𝑠𝑖𝑖𝑖𝑖𝑉𝑉𝑖𝑖𝑖𝑖 𝑖𝑖∈Fam (𝑖𝑖) +𝑏𝑏𝑖𝑖=0, (11) where b i is the body force at node i . 3.3. An intermediate ly-homogeniz ed peridynamic model for poro elastic 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 F H-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 porosit y. 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 full y-homogeniz ed models f or these cases. In this section, we introduce the intermediately -homogenized peridynamic (IH -PD) model for porous materials . The reader is referred to49 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 t o allows us to more accurately compute the failure behavior of porous material s compared with the F H-PD. At the 13 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) origi nally 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 “pre- damage 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 suspension6), all bonds associated with that point should be broken (𝑁𝑁𝑝𝑝=𝑁𝑁), mean ing 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 adopt ing the ideas used in the “ concentration -dependent damage ” (CDD) model originally proposed by Chen and Bobaru18 and used for modeling of damage induced by corrosion process es50-53. Here, we assume that the pre-damag e index at a point depends linearly on its porosity when zero- porosity material surrounds the point : 𝑑𝑑𝑝𝑝(𝒙𝒙)= 𝑃𝑃(𝒙𝒙)/𝑃𝑃𝑐𝑐, where 𝑃𝑃(𝒙𝒙) is the given local material porosity. The algorithm for implementing the pre-damage, at a node x, is: i. For each intact bond in the family of th e node, generate a random number r from a uniform distribution in (0, 1); 14 ii. If r is smaller than 𝑃𝑃(𝒙𝒙)/𝑃𝑃𝑐𝑐, then break this bond. Since each bond connects two nodes , updat e the pre -damage information at both nodes ; iii. Go to next bond in the family of bonds of node x. This procedure is applied for all nodes in the material. Here w e assume that the pore shape s 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 anisotrop y in porous materials . With the above algori thm, each PD bond (connecting material points 𝒙𝒙 and 𝒙𝒙�) goes through the procedure twice. If porosities at 𝒙𝒙 and 𝒙𝒙� are 𝑃𝑃(𝒙𝒙) and 𝑃𝑃 (𝒙𝒙�), respectively, the chance for the 𝒙𝒙−𝒙𝒙� bond to remain intact is (1−𝑃𝑃(𝒙𝒙)/𝑃𝑃𝑐𝑐)(1−𝑃𝑃(𝒙𝒙�)/𝑃𝑃𝑐𝑐). 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: 𝑑𝑑𝑝𝑝=1−(1−𝑃𝑃/𝑃𝑃𝑐𝑐)2 (12) in the limit of the horizon factor (the ratio between the horizon size and grid spacing) going to infinit y (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 Fig. 2 b we show the computed pre -damage index value for the input uniform porosities of 0.1, 0.3, 0.5 and 0.7, with a critical porosity of 1. When the horizon factor goes to infinit y, the pre-damage index value , according to Eq. ( 12), converge s to 0.19, 0.51, 0.75 and 0.91, respectively. Observe that t hese values match well the dominant colors in the pre -damage index maps shown in 15 Fig. 2b for the particular value of the horizon size equal to 4 cm (sample size is 1 m by 1 m) and a horizon factor m =4. In the next section, we use the IH -PD model to simula te elastic wave propagation in porous glass es. By computationally measuring the wave speeds corresponding to different porosities, we calculate the apparent elastic moduli , and compare them with those measured in experiments . Fig. 2. Peridynamic modeling of porosity via random pre -breaking of bonds; in (b): peridynamic pre -damage maps corresponding to different given porosities: 0.1, 0.3 , 0.5 and 0.7. The unit s for the axes is meter. 16 4. Model verification and validation for elast o-dynamic problems We verify (by performing convergence studies) and validate the IH-PD model for an elast o- dynamic problem: elastic wave propagation in porous glasses . We generate elastic waves by suddenly appl ying a force pulse over a small region in the sample and monitor their propagation through the sample . Experimen tal data regarding the apparent elastic modulus of the se porous glass es in terms of varying degrees of porosity exists .54 Fig. 3 show s the problem setting for the dynamic test with the IH-PD model . An elastic wave is induced in the porous glass (2D plane stress) sample by a sudden pulse applied (as a body force) over a region at the bottom. For zero porosity glass , elastic modulus is E = 46.1 GPa, and density 𝜌𝜌 = 2511 kg/m3.54 For porous glass or glass foam, the critical porosity can be higher than 0.9.6 Here we consider 1CP≈ for the porous glass simulations. 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 l ength) 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 patter ns 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 condition s (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 17 return from the sample’s boundaries. More specifically, we compute an average velocity by tracking the displacement of t he 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 “no- fail” condition for all bonds in the model. Obviously, the pre -damage corresponding to different porosities is present in the se calculations , but no “new damage” is allowed. Fig. 3. The problem setting used for verification and validation of the IH -PD model for elasto -dynamic conditions. Wave propagation is induced by a suddenly applied load over the “punch” zone shown in the right figure . Left: pulse profile; Right: sample dimension s and loading location. Fig. 4 shows the elastic waves propagating in the porous glass with the assumed porosities of P = 0 and 0.5, at different times: t = 4 µs, 40 µs, 100 µ s, and 160 µs, respectively. The colors indicat e the magnitude of the vertical component of the nodal velocity vector generated by the sudden pulse. For cases shown in Fig. 4, the horizon size is 40 m m, 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 -5 µstimepulse 1 MPa 0 pulse100 mm1000 mm1000 mm 18 PD model) of the heterogeneous micr ostructure. When porosity (pre -damage) is higher, the possibility for material nodes in the PD model to have no connections with surrounding nodes increases, especially when the horizon factor m is small. Because of this, t o compare the wave speeds in glass plate s with different porosities , especially for high porosities such as 0.9, we need to use a sufficiently large horizon factor m . Here we used m =16. Snapshot s at 100 μ s from the initial time step for the vertical component of the nodal velocities in porous glass sa mples are shown in Fig. 5 (Movies 1- 5 display the time -evolution for the vertical component of the nodal velocity vector over the first 200 μs , showing the progression of wave propagation for porosities of: 0, 0.1, 0.3, 0.5 and 0.7, respectively) . For all results shown in Figs. 4 and 5 the horizon size was 40 mm. Fig. 5 shows that higher porosity leads to lower wave propagation speed , as expected . In the IH -PD model, the pre -damage amount (corresponding to the given, input porosity) leads to this lower apparent modulus . In contrast, in FH-PD models , the lower input effective elastic modulus leads to this slowing in propagation of elastic waves in materials with higher and higher porosity . 19 Fig. 4. The vertical component of the nodal velocity vector, showing the progression of wave propagation in time for porous and non- porous materials : (a) P = 0, (b) P = 0.5. See also movies 1 and 3. 20 Fig. 5. The vertical component of the nodal velocity vector, showing the wave propagation at t = 100 μs for samples with different porosit y P. See also movies 1 -5. 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. Locations of the crest of the front wave versus time. Results obtained with a horizon size of 40 m m, and horizon factor m=8. time (s)0 50x10-6100x10-6150x10-6200x10-6Location of the first crest (m) -0.6-0.4-0.20.00.20.4 P = 0.0 P = 0.1 P = 0.3 P = 0.5 P = 0.7 21 Fig. 7 presents the convergence results from the IH -PD model for elastic wave speed in the porous glass versus porosity. To calculate the wave speed, we compute an averaged velocity by tracking the location of the crest ( along the line of x = 0) of the front wave during the time range 100~200 µs from the application of the pulse load. Fig. 6 shows the crest locations versus time, for different porosities . We use a linear fit for the vertical velocity data points with in the time range mentioned above to extract the wave speed velocit y shown in Fig. 7, for each porosity. Fig. 7. Convergence results in terms of the pressure wave propagation speed for the porous glass versus porosity: (a) δ-convergence for m=8, and (b) m -convergence for δ= 40 mm . As mentioned before, at high porosities, ( 0.7, 0.9), there is a relatively larger probability for peridynamic nodes to lose all of the ir 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 problem s 22 in porous material s with the IH -PD model w hen 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, w e compute the apparent modulus (EP) from the wave speeds numerically obtained in the IH -PD simulation s 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: 2(1 )P LECρν=− (13) where LC 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 match es well the experimental measurements (see Fig. 8 ). Fig. 8 shows the m-convergence for the normalized apparent modulus (EP/E, where E is the modulus at zero porosity) of the porous glass versus porosity. The limit of m -convergence approaches asymptotically the following functional behavior : 𝐸𝐸𝑃𝑃 𝐸𝐸=(1−𝑃𝑃/𝑃𝑃𝑐𝑐)2 . Note that the pre-damage index introduced by the given porosity was 𝑑𝑑𝑃𝑃=1−(1−𝑃𝑃/𝑃𝑃𝑐𝑐)2. Combining the se relationships lead s to 𝑑𝑑𝑃𝑃=1−𝐸𝐸𝑃𝑃 𝐸𝐸, which matches the classical formula for damage fro m continuum damage mechanics theory .57 This is another verification of the soundness of our proposed model. Note that, i n practice, one usually measures the apparent modulus, while the critical porosity is unknown ( or the modulus/porosity relation does not follow the equation 𝐸𝐸𝑃𝑃 𝐸𝐸 = (1−𝑃𝑃/𝑃𝑃𝑐𝑐)2). In 23 such cases , one can use the apparent modulus value and apply the IH-PD model to generate the pre-damage ind uced by the existing por osity, according to the formula: 𝑑𝑑 𝑃𝑃=1−𝐸𝐸𝑃𝑃 𝐸𝐸. Fig. 8. Normalized apparent modulus for porous glasses versus porosity , computed from the IH -PD model results shown in Fig. 7 . Critical porosity used here was taken as 1. Fig. 9 compares the normalized apparent modulus obtained from the IH-PD model simulation s with that from the experimental measurements .7 In the IH -PD model, we select m = 16, and δ = 40 m m. The normalized modulus/porosity relationship obtained by the IH-PD model matches well the experimental measurements. The computations are able to test porosity values outside of those that can be manufactured, and in the limit of porosity converging to the critical porosity we confirm the virtual stopping of wave propagation in the medium. With this model validation for elasto -dynamic problem s, we next validate the IH -PD model for fracture using quasi -static fracture of a porous rock. We consider fracture induced by three - point bending of Berea sandstone. 24 Fig. 9. Comparison between the IH -PD model results and experiments54 for the apparent modulus/modulu s ratio of porous glass (critical porosity taken as 1 in computations ). 5. Peridynamic modeling of damage evolution in porous materials We apply both the F H-PD and IH -PD model s for quasi -static fracture of a porous rock sample. We find that only the IH -PD model delivers damage pattern s and crack profiles similar to those seen in experiment s in.9 5.1. Description of experimental setup and results, and of the numerical model The experiment al setup shown in Fig. 10 was used9 for testing damage in Berea sandstone samples with an off -center notch and under quasi -static three -point bending. The experimental tests were conducted in specimens with notches of different lengths. The notch length determines the failure mode: when the notch i s short (5.65 mm), the main crack initiates from near the beam’s center and grows towa rds the loading point (see section 3.4 and Fig. 6 in Lin et al.9) in mostly mode I failure ; when the notch is longer ( 12.23 mm ), the beam fails by cracking from the tip of off-center ; this crack grows at an angle towards the loading region on the top boundary ( see section 3.4 and Fig. 8 in Lin et al.9) in mixed -mode. Fig. 7 in Lin et al.9 also shows that for the long notch, 25 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 th e observed damage sensitivity to notch length are reproduced by any of the two PD models described in this paper. Fig. 10. Geometry and loading configuration9 used in the peridynamic simulations. In Lin et al. ,9 a discrete element model (the “slightly overlapped circular particles interaction ” (SOCPI) model , see Fig. 3. i n Fakhimi, et al.58) was used to simulate the failure process. In the ir model,9 a synthetic rock is composed of rigid circular particles bonded together using normal and shear bonds at contact points (see Fig. 1. i n Fakhimi59). Each contact bond breaks if the applied normal or shear contact force exceeds the corresponding normal or shear bond strength s. With the breakage of a bond, a Coulomb failure criterion is applied to a contact provided the contact is under a compressive normal loading. This contact bond model is defined by five microscopic parameters: normal and shear spring constants, normal and shear bonds, and a friction coefficient. Through the interaction of many rigid circular particles with given micro mechanical properties, the macroscopic properties emerge. Hence, a process of trial and error is required to achieve appropriate macroscopic properties (e.g. Young’s modulus ) similar to those of a specific 26 rock, i.e., to calibrate the numerical model.9 In section 5.3 w e highlight advantages of our IH -PD model compared to the numerical model used by Lin et al.9 5.2. Problem setup We consider the 2D problem of quasi -static crack growth in a notched Berea sandstone specimen under three -point bending loading . The geometry and boundary configuration used in experiments9 and in our PD simulations are as shown in Fig. 10. Since the thickness of the specimen (26 mm) is small compared to its height and length, we model the problem as 2D plane stress. In experiments, the loading conditions were shown as point loads f or the three -point bending tests, and no details about the size/shape of the loaded areas were given. Personal communication with Fakhimi60 revealed the values shown in Fig. 10. Few details about the loading method were provided in Lin et al.9, regarding both their expe riments and simulations . Personal communication with Labuz61, clarified that in experiments the crack mouth opening displacement (CMOD) was contro lled with a constant rate of 1~2×10−4 mm/s . This rate allows us to consider the crack growth as quasi -static in our model . From Fakhimi ,62 we found that the discrete element method simulations9 used displacement -control led conditions applied at the top middle points of the beam to simulate the three -point bending loading. We use a similar approach to represent the loading . 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 shown9 stop before that. In our simulations, the imposed displacement is 2 μm between increments . We run 2,000 increment s to ensure full splitting of the samples. 27 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 porosit ies and thus different corresponding apparent modul i, 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 𝐺𝐺0,𝑃𝑃 𝐺𝐺0=(1−𝑃𝑃/𝑃𝑃𝑐𝑐)2. In Lin et al. ,9 only the apparent modulus ( 𝐸𝐸𝑃𝑃=10~12 GPa) was mentioned, while other parameters ( P, 𝑃𝑃𝑐𝑐, E, 𝐺𝐺0, 𝐺𝐺0,𝑃𝑃) were not given . From Jaeger et al.,34 we have the apparent modulus 𝐸𝐸𝑃𝑃=14.4 GPa and a corresponding porosity P = 19% for a Berea sandstone sample. With these values and assuming a critical porosity of 40%6, we calculate E = 52.24 GPa according to 𝐸𝐸𝑃𝑃 𝐸𝐸=(1−𝑃𝑃/𝑃𝑃𝑐𝑐)2. For the sandstone beam used in Lin et al. ,9 we apply this modulus/porosity relationship to calculate its porosity. By assuming 𝐸𝐸𝑃𝑃=11 GPa , we get the relative porosity 𝑃𝑃 𝑃𝑃𝑐𝑐= 54.1%. Next, we deter mine the zero -porosity fracture energy, 𝐺𝐺0,𝑃𝑃. Measured values of Berea sandstone’s apparent modulus ( 𝐸𝐸𝑃𝑃=14 GPa ) and fracture toughness ( 𝐾𝐾IC= 0.3 MPa m1/2) have been reported in the literature .63 The apparent fracture energy is found as 𝐺𝐺0,𝑃𝑃=𝐾𝐾IC2 𝐸𝐸𝑃𝑃=6.43 J/m2 (for plane stress problems) .64 Since 𝐺𝐺0,𝑃𝑃 𝐺𝐺0=𝐸𝐸𝑃𝑃 𝐸𝐸=(1−𝑃𝑃/𝑃𝑃𝑐𝑐)2, we obtain 𝐺𝐺0=24 J/m2. We use th e material properties of the constituent material of the porous medium ( E = 52.24 GPa and 𝐺𝐺0=24 J/m2) to calculate the bond information ( micro -modulus and critical bond strain) 28 for the IH -PD model , with relative porosity 𝑃𝑃 𝑃𝑃𝑐𝑐=54.1%. For the FH -PD model, we use the apparent mechanical properties: 𝐸𝐸 𝑃𝑃=11 GPa , 𝐺𝐺0,𝑃𝑃=5.05 J/m2. For the two- dimensional bond- based peridynamic model of an isotropic and homogeneous material in plane -stress conditions, the Poisson ratio is 1/3. T he Poisson ratio of the Berea sandstone varies between 0.3 and 0.35.9 While Poisson’s ratio can influence the path of a growing crack, these variation s are not expected to play a significant role in the f racture behavior of porous rocks.65 Peridynamic modeling of other poroelastic materials , for which bond- based model’s fixed 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 g rowth in poroelastic materials, we use the bond- based peridynamic formulation. In our sim ulations, we use the conical micromodulus (see Eqs. (5) and (9)), δ = 2 mm and m = 4. 5.3. Results and discussion We use the FH-PD and I H-PD model s (with the quasi -static loading condition s specified in the previous section) to simulate the cases shown in Fig. 10, and to examine the effect the pre- notch length has on the failure mo de and crack pattern s. 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 . Fig. 11 shows the simulation results generated by the IH-PD model in the short and long pre- notched Berea sandstone samples with 𝑃𝑃 𝑃𝑃𝑐𝑐=54.1%. Figs. 11c and 11d show the damage maps due to mechanical damage only, without considering the pre-damage. The fracture paths seen from 29 these damage maps are consistent with the experiment al findings on the effect of the notch length: the long pre-notch leads to crack initiating from the off -center pre-notch tip, while in th e short pre- notch 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 pre - notched 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 , when the sample splits in two . The numerical results in Lin et al. ,9 are limited to the intermediate stages of crack growth and do not include the crack paths through final failure (see Figs. 10 and 12 in Lin et al.9). It is not clear whether this is a limitation of the model they used or the authors only wanted to track the reported experimental results, which also are only reported well before full splitting . 30 Fig. 11. Damage maps from IH -PD model simulation: Pre -damage maps in the short (a) and long (b) pre -notched sample s; (c) and (d) : fracture paths shown by mechanical damage due only to loading (pre -damage not shown) at the total imposed displacement of 0.4 mm at the beam’s top- center ; (e) and (f) : damage maps from combin ing pre- damage with mechanical damage. Evolution of damage with increased loading is shown in movies 6-9. Fig. 12 presents the damage index maps for the short and long pre -notched samples using the FH-PD model. T he FH-PD model does not capture the influence of the pre-notch length on the crack initiation and propagation in the porous sample , observed experiment ally. With the F H-PD model, the crack i nitiates from the tip of pre -notch , even in the short notch case . Because of the simplif ying 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 off - center 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 ). a b c d e f 31 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 t o 1.9 mm. Fig. 12. Damage maps for the short ( a) and long ( b) pre-notch cases using the FH-PD model . See movies 1 0 and 1 1 for the evolution of damage . 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 are at the applied displacement of 0.4 mm at the top- center of the beam, using the IH-PD (Figs. 13a and 13b) and FH -PD (Figs. 13c and 13d) models , respectively . Results obtained with the IH -PD model ( Figs. 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 . Fig. 14 overlaps the data extracted from the numerical model and experiments used in Lin et al.9, with th e 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 32 The peridynamic results (from the IH -PD model) match well with the acoustic emission data9, showing damage accumulating at the right corner of the notch tip, while the numerical result s from Lin et al.9, show crack initiation from the left corner of the notch. Movies 12 and 13 show the mechanical damage evolution ( blanking nodes with zero mechanical damage index ) for the short and long pre -notched samples, respectively, using IH -PD model . In movies 14 and 15, the mechanical damage evolution for the short and long pre -notched samples, respectively, using FH -PD model are displayed. In movies 12- 15, the imposed displacement of the top- center of beam progresses from 0.6 μm to 0.4 mm. 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 occur ring earlie r than the damage that eventually leads to full failure, starting at the pre -notch tip. We n ote 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 in Lin et al. ,9 which is also the data shown in their other paper .58 a b c d 33 Fig. 13. Damage index maps resulting from mechanical damage (pre- damage and nodes with zero mechanical damage not shown) . (a) and (b): IH-PD model with short and long pre -notch, respectively (see movies 12 and 13) ; (c) and (d): FH -PD model with short and long pre -notch, respectively (see movies 14 and 15) . Red rectangles show notch geometr y. Fig. 14. Overlay of damage obtained by acoustic emission from Lin et al.9 (red diamonds), the IH -PD model ( multi - color squares , and using the legend in Fig. 1 3), and the numerical model in Lin et al.9 (black dots) . The o ff-center black rectangle shows the notch geometry and location . In Table 1, we compare the peak load s (the maximum load values) reported in the experiments with those obtained with our models. The results from the IH-PD model match well the experimental data. The peak load obtained by FH -PD model for the short notch case deviates significantly from the experiments because of the different mode of failure obtained in this case . Table 1. Comparison of peak load obtained from experiments with those from the peridynamic models Peak load/thickness by experiment9 Peak load/thickness by IH -PD model Peak load/thickness by FH -PD model Short notch 59.535kN m 49.952 kN m (Relative error=16%) 77.476 kN m (Relative error=30.1%) Long notch 46.48kN m 49.694 kN m (Relative error=6.9%) 49.376 kN m (Relative error=6.2%) 34 Compar ed 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 /porosit y (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, t he constitutive model in the contact bond model9 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 process es) are required to achieve appropriate macroscopic properties similar to those of a specific rock. 6. Conclusions In this paper we presented an Intermediately -Homogenized Peridynamic (IH -PD) model for simulating the elastic and failure behavior of porous material s. In this model, porosity wa s represented by an initial peridynamic damage (“pre -damage”) , introduced by breaking bonds, stochastically, to achieve a given porosity. The peridynamic micromodulus wa s computed using the elastic modulus of the porous medium constituent material. We validate d the model for elastic wave propagation using experimental data for wave propagation speeds and the apparent elastic 35 modul i in porous glass with different porosities . While simple homogenization methods work well for linea r 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 experiment al observations . To answer the question “ how much homogenization is too much?” when modeling fracture processes in porous materials , we compar ed two peridynamic models : one that homogenize d the material to a greater extent (the Fully -Homogenized Peridynamic model, F H- PD) and the newly introduced IH -PD model. W e stud ied a quasi -static crack growth problem in a brittle porous rock to understand the difference between the models’ responses . We f ound that the fully -homogenized model (FH -PD) failed to capture the experimentally obs erved fracture behavior in Berea sandstone s 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, reproduce d perfectly the observed crack growth behavior and its dependency on the length of the pre -notch. We conclud e that for problems in which the microstructure ha s a large effect on the failure behavior, a ful ly homogenized strategy will not work to correctly capture the fracture and damage evolution. For a predictive model, s ome 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 th is extra microstructure information that allow s the computed damage to initiate and evolve in a way similar to that in the real porous medium . 36 Acknowledgements: This work was supported by a grant from ONR (program manager: William Nickerson ) and by the AFOSR MURI Center for Material Failure Prediction through Peridynamics (program manager s: Drs. Jaimie Tiley, David Stargel , Ali Sayir , Fariba Fahroo, and James Fillerup ). This work was completed utilizing the Holland Computing Center of the University of Nebraska, which receives support from the Nebraska Research Initiative.