Microstructural controls of anticrack nucleation in highly porous brittle solids

Porous brittle solids have the ability to collapse and fail even under compressive stresses. In fracture mechanics, this singular behavior, often referred to as anticrack, demands for appropriate continuum models to predict the catastrophic failure. To identify universal controls of anticracks, we link the microstructure of a porous solid with its yield surface at the onset of plastic flow. We utilize an assembly method for porous structures, which allows to independently vary microstructural properties (density and coordination number) and perform discrete element simulations under mixed-mode (shear-compression) loading. In rescaled stress coordinates, the concurrent influence of the microstructural properties can be cast into a universal, ellipsoidal form of the yield surface that reveals an associative plastic flow rule, as a common feature of these materials. Our results constitute a constructive approach for continuum modeling of anticrack nucleation and propagation in highly porous brittle, engineering and geo-materials.

The compression of porous brittle materials can lead to localization of compaction or compacting shear bands [1][2][3] . This singular behavior, referred to as anticrack, originates from microstructural failure processes and has been observed in the compression of porous sandstone 4,5 , submarine landslides 6 , firnquakes 7 or snow avalanches [8][9][10] . The latter is a particularly spectacular consequence of the nucleation of an anticrack under mixed mode loading: Snow as a highly porous material continuously changes its microstructure under thermodynamic forcing 11 . Under a high temperature gradient, a layer of small rounded crystals can turn into a weak layer of large faceted crystals 12 . Buried in a snowpack on a mountain slope, such a layer is always subjected to mixed-mode (shearcompressive) loading and the nucleation of an anticrack can lead to widespread anticrack propagation and dangerous slab avalanches 3,8 . To understand the generic mechanisms underlying these processes, irrespective of a particular material, it is necessary to decipher how the macroscopic behavior of mixed-mode anticrack nucleation originates from the microstructure.
From a continuum point of view, the nucleation of mixed-mode anticracks is controlled by the complex interplay between the yield surface, plastic flow rule and strain softening. It has been suggested that non-associated plasticity is necessary to reproduce anticracks in some porous sandstones 13 . In contrast, recent work on snow has demonstrated that anticrack nucleation and propagation in snow can be simulated using a continuum damage model with associated plasticity, if complemented by a modified strain softening law 3 . However, these question cannot be satisfactorily answered at the continuum level where involved assumptions cannot be traced back to and justified from the microstructure.
In the past decades, the extensive use of X-ray micro-computed tomography (XRCT) facilitated unprecedented insight into porous microstructures in all fields of material science 14 . Numerical simulations based on XRCT images allow to faithfully characterize and constrain the properties of the microstructural network (such as connectivity) on the mechanical behavior of different materials such as soils 15,16 , lunar soils 17 , snow [18][19][20][21][22] , rocks 23,24 or concrete [25][26][27] . Simulations based on XRCT microstructure images are nowadays seen as complementary (numerical) experiments which can be repeatedly performed with different material properties, loading states and boundary conditions. This opens excellent opportunities in addition to (destructive) laboratory experiments which naturally lag behind in view of versatility and parameter variability. However, high resolution microstructure-based simulations still come with considerable computational requirements. Thus systematic A common approach to investigate microstructural controls on the mechanical behavior of porous solids are particle based methods such as DEM 28 . To access generic mechanisms relevant for a wide range of microstructural conditions, particle assemblies with volume fractions close to a dense packing are, however, not suitable. Rather a significant variation of the (mechanically) relevant parameters, i.e. volume fraction and the coordination number, are required to understand the controls of porosity and matrix connectivity on anticracks. As demonstrated by Gaume et al. 29 , Baxter's sticky hard sphere (SHS) model 30 can be conveniently utilized as a generic assembly method for DEM to independently prescribe volume fraction and coordination number, for very low porosity aggregates. In addition, a mapping of particle properties to continuous two-phase microstructures acquired by XRCT is principally feasible in the sense of a stochastic reconstruction, by matching two-point correlation functions 31 .
Here, we employ this methodology to systematically investigate a large ensemble of diverse microstructures with a wide range of coordination numbers and volume fractions within DEM under mixed-mode loading conditions ( Fig. 1). From all simulations, we derive a single form of the continuum yield surface and the plastic flow rule relevant for anticrack nucleation. The results provide fundamental insight how the continuum mechanical behaviour of very loose and brittle solids is controlled by volume fraction and microstructural connectivity.

Results
Mixed-mode anticrack nucleation. Heterogeneous compaction in brittle, highly porous solids is highlighted in uniaxial compression DEM simulations (Fig. 1C). The deformation occurs in a localized manner within a compacted region (illustrated by large vertical strain rates and distributed broken bonds in Fig. 1B,C) while the rest of the sample remains undisturbed. The generic behavior under uniaxial compression is shown in Fig. 2A,C for different volume fractions and coordination numbers. Stress-strain curves follow a quasi elastic brittle behavior i.e. (1) an almost linear elastic phase including non-cascading bond breaking events for low deformations up to (2) catastrophic failure followed by (3) strain softening. Increasing values of volume fractions and/or increasing values of coordination number lead to increasing elastic modulus and compressive strength. In addition, the increase in the broken bonds percentage is strongly correlated with the amount of softening after failure. After the strain-softening phase, samples collapse under the load leading to the localization of a compaction band 2 (Fig. 1B) which is generally referred to as anticrack 4 . The residual stress after softening is the result of a competition between bond breaking and formation of new frictional contacts.
The yield surface of the samples (Fig. 2B,D) was evaluated from mixed-mode loading simulations. All porous solid samples can fail under tension, shear, compression, and mixed-mode loading states leading to a closed yield surface. Tensile and shear strength are similar in magnitude, and both lower than the compressive strength. While the shape of the yield surface does not appear to be significantly affected by the initial microstructures, its size drastically increases with increasing volume fraction and coordination number, in line with the observation made for the compressive strength. Figure 3 shows the volumetric response of the samples for different loading angles. A transition between volume increase (expansion), and volume reduction (compaction) after failure is found for a loading angle of ∼ 67.5 • . For the latter angle, the vertical displacement of the clump u z remains almost constant during shearing. For larger compressive stresses, i.e. lower loading angles, fracture is followed  On the other hand, we did not find any trend for the parameter M which varies between 0.6 and 1 with an average of ∼ 0.8. Given the above scaling, it appears natural to search for a universal form of the yield surface parameterized based on contact density only. To this end, we define a transformed stress coordinates (σ * , τ * ) in which the yield surface is a unit circle: To satisfy Eq. (1), σ * and τ * are defined as Using this normalization based on contact density only, all failure points collapse on a universal yield surface defined by Eq. (4) (Fig. 5A). Interestingly, the apex of the ellipse corresponds to a loading angle ψ between 45 • and 90 • with a mean value ∼ 67.5 • . This angle corresponds to the transition between samples' compaction and expansion. To evaluate the flow rule of the samples, we evaluate the plastic flow angle θ = arctan u z u x based on the simulations. To verify the assumption that the samples have an associative plastic flow rule, we relate θ to the derivative of the yield surface ∂τ ∂σ , which depends on the loading angle ψ . For the sake of universality, calculations are performed in the normalized stress space. Details of the calculations are presented in the Methods section. Figure 5B shows the normalized plastic flow angle θ * as a function of the normalized loading angle ψ * . All data points appear to collapse on a single universal linear curve given by Eq. (16) which indicates that highly porous brittle solids follow an associative flow rule. Note that the loading angle ψ +/− corresponding to a transition between expansion and compaction directly depends on the slope of the critical state line M, as well as the parameter β , which characterizes the cohesion of the material according to For a cohesionless material ( β = 0 ), ψ +/− = arctan M (critical state line) while for a porous solid with the same compressive and tensile strengths ( β = 1 ), ψ +/− = π/2 which is typical for metallic foams 32 . (1)

Discussion
An explanation as to why the yield surface is predominantly controlled by volume fraction and contact density is suggested by the fact that for arbitrary particle-based two-phase systems, the link between physical and structural properties is completely embodied in the hierarchy of n−point correlation functions 33 . The lowest point order ( n = 1 ) is only determined by the volume fraction as the widely accepted, most important parameter. In the next higher point order ( n = 2 ), the contact density ν c controls the short-range expansion of the corresponding correlation function 34 . This implies that the contact density must be considered as the the second-most important parameter (beyond the volume fraction) in the combined expansion of the microstructure, with respect to point-order and spatial range. This naturally raises the question to which extent these two parameters can explain the macroscopic yield surface of highly porous brittle solids.
To this end, we considered a zoo of highly porous microstructures generated from the sticky hard sphere ensemble 30 to independently control volume fraction and contact density. Within numerical uncertainties, both parameters are sufficient to quantify the yield surface and plastic flow rule of the porous solids. The yield surface of the samples collapse on a master curve (unit ellipse), after rescaling by the contact density. In addition, we have numerically shown that the plastic flow rule is associative. As a consequence, the volumetric response critically depends on the applied pressure. For large normal stresses ( σ > σ 0 (1−β) 2 , i.e. on the right side of the apex of the ellipse), the plastic flow angle is negative leading to the volumetric collapse of the samples (plastic compaction). This process is referred to as a mixed-mode (shear-compression) anticrack. Together with previous results 29 on the unique link between elastic modulus and contact density, our present results highlight the universal microstructural control of this quantity for the (pre-failure) mechanical behavior of highly porous brittle solids.
The proof of associativity presented here allows us to justify a strong assumption that was previously required to model anticrack nucleation leading to catastrophic slab avalanches in the recent work of Gaume et al. 3 . In this approach, a continuum model for dynamic anticrack propagation in snow has been developed based on critical state plasticity theory. These authors used a similar ellipsoidal (Modified Cam Clay) yield surface and made the assumption of an associative plastic flow rule to simulate the volume change during snow deformation. In addition to the proof of associativity, the present work allows us to evaluate the parameters of the model based on a single microstructural quantity: the contact density.
As aforementioned, the advantage of this approach is that parameters of the Sticky Hard Spheres model (volume fraction and coordination number) can be directly evaluated based on X-ray tomographic images by matching correlation functions 31 . While several recent numerical studies have analyzed the mechanical response of porous brittle solids like snow based on the real samples' microstructures 18,19,21 , it has been recently shown that simplified structures made of spherical particles can be used to reproduce accurately important features of snow mechanics in the brittle range for different processes: failure initiation 35,36 , crack propagation 9,10 , snowflake fragmentation 37 , wind blowing snow 38 , snow granulation 39 and avalanche impact pressures 40 . One the one hand, this simplification makes us loose important information about the microstructure. However, on the other hand, it allows us to significantly fasten the simulations and thus make detailed parameter sensitivity studies, which is not possible with highly detailed representations of the microstructure. In addition, the individual particle properties were chosen according to the ice mechanical behavior which constitute the solid matrix of snow. Yet, it was shown 29 that the bulk elasticity and strength of the samples linearly scales with the particle elasticity and bond strength, respectively. This is reflected by Equations (2) and (3) in which the yield surface parameters are presented in a normalized manner, which allows to apply our results to other highly porous materials with different solid matrix properties. The presented DEM model is able to reproduce the nucleation and propagation of anticracks observed in porous layers of snow in the brittle range 3 . The formation of new cohesive bonds during the simulation would allow us in the future to reproduce the ductile-to-brittle transition in snow and thus the repetitive formation and reflection of compaction bands, as observed by Barraclough et al. 2 . Additionally, www.nature.com/scientificreports/ implementing a particle breakage criterion 41,42 would enable us to simulate more complex types of localized deformation such as erratic and oscillatory compaction bands, observed in the compression of rice crispies 1 . Evaluating the conditions for the onset of localization of compacting shear bands or anticrack nucleation in porous rocks is a great challenge 5,43 and has implication for the understanding of deep earthquakes. Many associative plasticity models have been developed for porous rocks mechanics but localization is only possible with strain-softening, which occurs only under conditions of volume increase (dilation) according to classical critical state soil mechanics (CSM) models inspired by the behavior of granular materials. However, localized deformation was also reported under compressive stresses leading to compaction or compacting-shear bands [1][2][3][4][44][45][46] . Some researchers tried to overcome the inability of classical associative CSM models to reproduce anticracks or compaction bands through the development of non-associative flow rules. However, important discrepancies between experiments and model predictions have been found concerning localization features 13,43,47 . Wong et al. 47 attributed differences between experimental observations and localization analysis to the "inadequacy of the non-associative constitutive model to capture the partitioning of several damage mechanisms, including the growth and coalescence of stress-induced microcracks and pore collapse". The present analysis suggests that the plastic behavior of highly porous solids ( φ < 0.35 ) is associative. The fact is, highly porous brittle solids undergo significant softening, even under large compressive stresses 1-3 . Gaume et al. 3,29 suggested that the solid structure of porous solids under compression is actually under tension, which jeopardizes the continuum assumption. Hence, they proposed a modified hardening law based on the norm of the volumetric plastic strain, which leads to a shrinkage of the yield surface, even under compression until it corresponds to a point in the origin of the stress space, when cohesion is set to zero. The behavior under compression (on the cap of the yield surface) thus becomes similar to tensile mode I crack opening, which is refereed to as anticrack here (negative mode I). Here, simulations are performed for relatively low deformations ( < 5 %) corresponding to the nucleation of the mixedmode anticrack. For larger strain values, the yield surface changes significantly 3 . The value of β decreases due to cohesion loss and σ 0 increases because of the creation of new frictional contacts. The decrease of β induces a change in the value of the slope of the apparent critical state line (CSL). In the current formulation of the yield surface, M is the slope of the cohesionless CSL. Hence, the slope of the cohesive CSL which starts from the coordinate ( −βσ 0 , 0) is M ′ = M/ √ 1 + 2β . Most porous geomaterials, including rocks and snow have values of β between 0.1 and 0.5 typically 19,48,49 leading to an increase of the apparent critical state line after failure of 10 to 40%. Hence, the difference observed between the plastic flow direction and the normality condition to the initial yield surface of highly porous solids is related to a sudden change of the shape of the yield surface induced by the post-peak softening.
Finally, we note that the choice of volume fractions for the simulations is motivated by the applicability to highly porous systems below the close packing density where the used Monte Carlo method of mono-disperse spheres 50 may fail. At high volume fractions (close to the random close packing density, as obtained by Kun et al. 51 using a different packing algorithm) the observed compaction band regime is generally complemented by a regime of predominant shear bands. It would be interesting to study such crossover behavior in the future using a microstructure assembly method that is able to cover the whole range of volume fractions.
Methods initial states. The initial states of DEM simulations are generated using Baxter's model of Sticky Hard Spheres (SHS) 30 . SHS allows us to generate a random assembly of percolating spheres of given volume fraction φ and coordination number z c (average number of contacts per sphere). Note that the SHS ensemble contains a percolation transition where the fraction of particles in the load-bearing backbone of the microstructures approaches zero. Naturally, using SHS state in the non-percolating phase is meaningless for their interpretation as a porous brittle solid. The details of the SHS model can be found in Gaume et al. 29 .
We used the same microstructures as in Gaume et al. 29 . They are generated through Monte-Carlo simulations based on the SHS model. Through this statistical approach it is possible to generate different realizations of a sample with the same initial volume fraction φ and initial average coordination number z c but random microstructures. For every combination of φ and z c (Fig. 6), three realizations were simulated.
Discrete element simulations. We used PFC3D v5 by Itasca 52 to simulate the mixed-mode loading of porous cohesive granular samples inside a cubic box of 1 cm side-length. The simulations are performed without gravity which allowed us to use a large particle density ρ = 10,000 kg/m 3 to speed up the simulations.
Initial positions (x p0 , y p0 , z p0 ) and particle radii r are obtained from Monte Carlo simulations. Simulations are performed with N = 2,048 particles which was shown large enough to prevent size effects 29 .
All particles close to the bottom of the system ( z p0 < 3r ) were fixed in translational and rotational movement (green particles in Fig. 1A). On the side, periodic boundary conditions were applied (transparent gray domain in Fig. 1A). At the top, a rigid clump was created to apply the mixed-mode loading.
Initial contacts between particles are bonded using the parallel bond model described in detail in Gaume et al. 29 . During the simulation, bonds can only break and new frictional cohesionless contacts can occur. Bonds have an elastic modulus of E = 1 GPa and a normal-to-shear stiffness ratio κ = 3 . Bonds can break under shear and tension and have specific shear and tensile strength σ tb and τ sb , respectively. We chose σ tb = τ sb = 1 MPa. This leads to a bond failure strain of ǫ fb = 0.001 of the order of that of ice.
Loading and stress measure. Load-controlled simulations are performed by applying an increasing force F on the clump with different loading angles. The force is defined as � F = F(sin ψ � x − cos ψ� z) where ψ is the loading angle. A loading rate Ḟ = 0.05 N/s was chosen sufficiently small to ensure quasi-static conditions. In addition, we verified that the inertial number 53  First we identify a strain range in which failure occurred (transparent gray range in Fig. 7). We define the normalized kinetic energy as: where v c is the clump velocity vector. After failure, E * k significantly increases in all simulations. Hence, we define a strain range for failure by searching for E * k < ζ . In a second step, we search for the maximum equivalent von Mises stress q within the previously defined strain range. The equivalent von Mises stress is defined as with s as the deviatoric part of the stress tensor given by where I is the identity matrix, p the pressure and, σ is the Love homogenized stress tensor. The local maximum of q within the strain range defined above is defined as "failure point". We verified, based on the stress-strain curves, that a value ζ = 0.3 led a very robust detection of failure. We also checked that the presented results remain valid for other failure criteria such as the second order work 54 or a criterion based on the number of broken bonds.
Plastic flow in the normalized stress space. The normalized loading direction is characterized by with Combining Eqs. (4) and (12)    www.nature.com/scientificreports/ Hence, plugging µ * from Eq. (14) into Eq. (12) gives For an associative plastic flow, the plastic flow angle θ * is given by Given the following relationship between the original and the normalized yield surfaces: we obtain The last piece of the puzzle consists in relating µ * to µ in Eq. (16). This is done by combining Eqs. (5), (6) and (13) and leads to which can easily be solved for µ * .
Received: 13 January 2020; Accepted: 9 June 2020  Figure 7. Identification of failure: normalized kinetic energy E * k with the determined transparent gray search range and the equivalent von Mises stress q as well as normal stress σ . The identified failure is indicated by a black circle.