Microscale mechanical modeling of deformable geomaterials with dynamic contacts based on the numerical manifold method

Micromechanical modeling of geomaterials is challenging because of the complex geometry of discontinuities and potentially large number of deformable material bodies that contact each other dynamically. In this study, we have developed a numerical approach for micromechanical analysis of deformable geomaterials with dynamic contacts. In our approach, we detect contacts among multiple blocks with arbitrary shapes, enforce different contact constraints for three different contact states of separated, bonded, and sliding, and iterate within each time step to ensure convergence of contact states. With these features, we are able to simulate the dynamic contact evolution at the microscale for realistic geomaterials having arbitrary shapes of grains and interfaces. We demonstrate the capability with several examples, including a rough fracture with different geometric surface asperity characteristics, settling of clay aggregates, compaction of a loosely packed sand, and failure of an intact marble sample. With our model, we are able to accurately analyze (1) large displacements and/or deformation, (2) the process of high stress accumulated at contact areas, (3) the failure of a mineral cemented rock samples under high stress, and (4) post-failure fragmentation. The analysis highlights the importance of accurately capturing (1) the sequential evolution of geomaterials responding to stress as motion, deformation, and high stress; (2) large geometric features outside the norms (such as large asperities and sharp corners) as such features can dominate the micromechanical behavior; and (3) different mechanical behavior between loosely packed and tightly packed granular systems.


Introduction
Numerical modeling of microscale mechanical behavior of geomaterials (soils and rocks) is of great importance for understanding and predicting material constitutive and geomechanical behavior at larger scales in subsurface engineering activities such as unconventional hydrocarbon production [22], nuclear waste disposal [23], and CO 2 sequestration [21]. Some unique features of geomaterials are that they are naturally stressed and heated to variable degrees, and often fluid filled. The microscale structure of geomaterials contains minerals, pores, and fractures of complex shapes that evolve as a result of coupled fluid, heat, mechanics, and chemical reactions. In order to understand such dynamic multiphysics problems, in the past decade, new technologies have been developed for visualization and characterization of the multiphysics processes of geomaterials at the pore scale, microscale, and even smaller scales [16,20,31,33,38].
However, numerical modeling of mechanical processes in geomaterials at the microscale is challenging because of the computational geometry associated with (1) complex evolving geometric features that are discontinuous, and (2) multiple deformable material bodies with dynamic contacts. At the microscale, geomaterials exhibit more complex geometric features that cannot be simplified easily and therefore lead to discontinuities in physical fields (i.e., displacements and stress for mechanics). Moreover, in addition to the force balance, microscale mechanical processes may involve large deformation, large translational or rotational displacements, and dynamic contacts. The grand challenge of computing contacts is to identify when and where contacts occur among many blocks [27]. This is complicated by the fact that these blocks are moving, deforming, and in some cases breaking apart. In turn, motion, deformation, and breakage of blocks are impacted by contact forces, thus constituting a serial process. Inaccurate calculation of contacts, therefore, can lead to completely erroneous overall system behavior.
For systems involving both continua and discontinua such as single fractures at the asperity scale, various models have been developed to simulate fluid flow [1,39], heat transfer [18], reactive transport [28], or their couplings [3]. But numerical models that fully address contacts and deformation in combination are much rarer due to the aforementioned challenges. For computation of granular materials, a number of numerical approaches based on the discrete element method (DEM) have been developed, including those of Cundall [5], Houlsby [8] and Andrade et al. [2]. Though DEM is designed for computation of discontinuous bodies with dynamic contacts, the assumption of rigid bodies, the use of explicit time iteration, and the limitation of interpolation fields for continuum mechanics limit its accuracy for dealing with realistic geometric features or dynamic contacts. These limitations are in addition to its potentially high computational cost. On the other hand, a number of approaches have been developed using the finite element method (FEM) to enable contact calculations (Puso and Laursen, [40]; [24]). Though FEM is powerful for continuum mechanics possibly involving large deformation, the common contacting algorithms developed in those approaches assume that contacts are along prefixed contact pairs, or generally do not involve multi-body system or a large number of contacting pairs. Some other approaches have been developed for contacts in granular systems, but these do not overcome the fundamental limitations associated with prefixed contact pairs, a limited number of contacts, or the limited accuracy for either continuum or discontinuum calculations. Therefore, there is a need to develop a powerful toolset that can accurate capture both continuous and discontinuous behavior involving dynamic contacts that may involve a large number of interacting material bodies.
The numerical manifold method (NMM, [25,26]) is a promising method for analyzing both continuous and discontinuous media involving dynamic contacts. NMM is based on the theory of mathematical manifolds. The numerical meshes of NMM consist of mathematical covers and physical covers. The mathematical covers overlay the entire material domain and the physical covers are divided from the mathematical covers by boundaries and discontinuities. Based on this dual-mesh concept, both continuous and discontinuous problems can be rigorously solved by flexibly defining physical covers. In the past two decades, NMM has been successfully applied to mechanics analysis of both continuous and discontinuous geologic media [17], involving higher-order interpolation [4], fracture propagation [36], wave propagation through fractured media [6], analysis of slope stability [7,19], and microscale-macroscale modeling of fracturing of sandstone [34]. Most recently, the authors developed a number of models for analyzing flow and fully coupled hydromechanical processes of fractured porous media at different scales ( [13][14][15], [9]).
In this paper, we present development of an approach with rigorous treatment of dynamic contacts in deformable geomaterials for microscale mechanical analysis based on the NMM. We first present a general mathematical description of the problem of deformable geomaterial bodies with dynamic contacts in Section 2. Then, we introduce the approach of modeling continuum and discontinuum mechanics based on the NMM, including the fundamentals for global continuous and discontinuous interpolation, and a rigorous multi-step approach for contact calculations in Section 3. In Section 4, we apply this approach to a number of examples at the microscale, including a single fracture at the asperity scale, granular systems of aggregated clay settling, compaction of a loosely packed sand, and failure of a cemented marble sample. In Section 5, we conclude and provide perspectives for future directions.

Mathematical statement of deformable geomaterial bodies with dynamic contacts
A single continuum or discontinuum material body in dynamic states satisfies force balance: where σ is the stress tensor, f is the body force vector, ρ is the density of the material body, u is the vector of displacements, and t is time. At steady state, the right-hand side representing the inertia terms = 0. The force term is not only a result of loading (F l ), but is also present when material bodies are discontinuous and are in contact (i.e., the contact force, F contact ). Thus, the complete expression of the force includes both terms as follows: In order to include continuum and discontinuum mechanics in a unified way, we express the displacement as: where the displacement includes not only the term that is due to deformation ∫εds, but also the translational u tr and rotational u r displacements. When the deformation is assumed small, strain is known as the spatial derivative of displacements.
With introduction of translational and rotational displacements, the motion of a material body can be described. When a material body does not interact with other material bodies, in response to loading, it may move or deform before the stress reaches the strength of the material. Various types of continuum constitutive laws with different stress-strain relationships have been developed to describe the continuum behavior of a material body, which can be expressed generally as: where g is a general function, which can be linear, nonlinear, or rate-dependent. For example, Hooke's law is most widely used law to describe the linear elastic relationship between the stress and strain tensors. When a material body interacts with other material bodies, conventional continuum constitutive laws are not sufficient to describe the relationships between contact forces and displacements and/or deformations associated with contacts. Between different material bodies, it is more natural to use the distance d and the relative displacement in the direction along the contacting face ⟦u s ⟧ to describe their state of contacts. Here the distance d is the distance between a contact pair (i.e., the exact contacting vertices and/or surfaces between two material bodies). The distance d includes time-dependent influences of relative displacement and relative deformation of the two material bodies between the potential contacting faces. Thus, a general expression of contact force as a function of d and ⟦u s ⟧ can be expressed as: where ⟦u s ⟧ denotes discontinuity of displacements (i.e., relative displacements) between a contact pair in the direction along the contacting face, and h is a function of d and ⟦u s ⟧, both of which may be dynamically changing before reaching equilibrium.
With or without considering motion or deformation of material bodies, there are always three possible relative positions between two material bodies A and B: (1) A and B are separated (no contact anywhere); (2) A and B are bonded to each other on the vertices or surfaces of the body; and (3) A is moving along B while A is in contact with B, and such a motion is sliding. For these three different contact states, their respective functions h are different.
When A and B are separated, in most cases there are no contact forces between them. Thus, we have: In rare cases, for example when there is an electrical double layer, there are forces between A and B, but this effect is not currently included in this study.
When A and B are bonded, the distance between A and B on the contacting pair (vertices and/or edges in 2D) should be zero, and the relative displacement along the contacting face ⟦u s ⟧ should be zero as well, satisfying: Note that in Eq. (7), we still consider the situations when A and B are moving and deforming, but remaining bonded to each other.
When A is sliding along B while A is in contact with B, the contact occurs at the sliding face. On this face, there is no distance between A and B, while the sliding satisfies certain friction laws. We use the Coulomb's law of friction to describe sliding along the contact face for the sliding state: where F s denotes the contacting force in the direction of the sliding face, F 0 n denotes the effective contact force in the direction normal to the contact face (here effective means eliminating force associated with water pressure), φ is the friction angle, and sgn(⟦u s ⟧) means the direction of F s that depends on the direction of relative shear displacement.
Although in Eqs. (7)- (8) there is no explicit description of the contact forces, the distance and/or relative displacement constraints function like implicit boundary conditions. These constraints on the contact faces can be imposed with different approaches, such as with a penalty method [25], Lagrange multiplier method [11], and the most advanced approach based on the variational inequality method [37], leading to force terms as contact forces. Numerical implementation of these boundary constraints will be introduced in more detail in the next section.
So far, we have described two individual material bodies in separated, bonded, and sliding contact states. In dynamic conditions, these contact states may be changed as follows: If material bodies A and B were not in contact (separated), but become bonded later, constraints in Eq. (7) should be added. If A and B were not in contact, but then they are in sliding state, constraints in Eq. (8) should be added. If material bodies A and B were bonded but become separated later, we need to consider one condition: In order to separate A from B, the effective force in the direction normal to the contacting face F 0 n needs to be larger than the force associated with the tensile strength T. That is because in a rather long geological period as a result of thermal-hydro-mechanical-chemical (THMC) processes, adhesion and/or tensile bond strength can be gained between A and B along their contacting faces. This criterion for separating two bonded bodies in the direction normal to their contact face can be expressed as: If Eq. (9) is satisfied, the two material bodies are no longer in contact; therefore, the constraints in Eq. (8) should be removed. If material bodies A and B were bonded but transfer to a sliding state, we need to consider one condition: In order to initiate sliding of A and B against each other, the force in the direction of contacting face F s needs to be larger than the shear strength S. The shear strength could be gained as a result of THMC processes, consisting of frictional force (satisfying Coulomb's law of friction) and cohesive force F cohe . This criterion for shearing two bonded bodies in the direction along their contact face can be expressed as: where φ′ is the internal friction angle. If Eq. (10) is satisfied, the two material bodies are transferred from bonded to sliding state. Comparing Eq. (7) and (8), we find that the constraint in the direction normal to contacting face should be retained, while in the sliding direction, the constraint needs to be changed. If A and B slide on each other but become separated, constraints in Eq. (9) should be removed. If they become bonded, the constraints should be modified in the direction of the contacting face so that no relative shear displacement will occur. With Eqs. (6)-(10), we are able to describe the three different contact states of two material bodies, dynamic changes of these contact states, and criteria that need to be satisfied for changes of the contact states.

Fundamentals of NMM: global approximation
The numeral manifold method (NMM) [25] is based on the concept of a "manifold" in topology. In NMM, independent meshes for interpolation and integration are defined separately. Based on this unique definition, a non-conforming mesh (not necessarily conforming with the physical boundaries) can be used as a mathematical mesh. Different local approximation functions can be constructed and averaged to establish global approximations for both continuous and discontinuous analysis. In NMM, independent mathematical and physical covers are defined. A mathematical cover is a set of connected patches that cover the entire material domain Ω. For example, we can use a quadrilateral patch, a circular patch, or a rectangular patch as a mathematical cover (such as the A, B, C in Fig. 1). Features such as density, shape, and sizes of these mathematical patches define the precision of the interpolation. The physical patches are mathematical patches divided by boundaries and discontinuities, determining the integration fields. The union of all the physical patches forms a physical cover. For example, mathematical patches B and C become smaller areas as physical patches B 1 and C 1 divided by the exterior boundaries of the material domain Ω, while mathematical patch A is divided into physical patches A 1 and A 2 by the inner discontinuity and the boundaries of Ω. The overlapping areas by multiple physical patches are defined as elements. As a result, the model domain Ω is discretized into five elements: A 1 B 1 C 1 (the overlap of physical patches A 1 , B 1 , and C 1 ), A 1 C 1 , A 2 C 1 , B 1 C 1 , and C 1 . From Fig. 1, we can see that the shape of the mathematical patches can be flexible; the relative location of the mathematical patches within the model domain can also be arbitrary (only if satisfying Ω ⊂ A ⋃ B ⋃ C), and the number of physical patches on each element can be flexible.
On each physical patch, a local function is assigned, such as one that is constant, one that is linear, or any function that is able to capture the behavior of the solution on the patch. The weighted average of the local patch functions constitutes the global approximation. For example, if using linear local functions, we can construct a global second-order approximation (Fig. 2a, [32]). If using a local function with a jump of the first derivative of the physical variable, we can simulate the weak discontinuity crossing patches and elements (Fig. 2b, [12]). Or most commonly, if using discontinuous local functions, we can simulate fractures (Fig. 2c, [13,15]). With this dualmesh concept, the NMM is capable of simulating both continua and discontinua with explicit geometric representation and flexible numerical approximation.
With the concept of global approximation, NMM can be related to other numerical methods as shown in Fig. 3. If using a bilinear weight function on a rectangular patch with a  [35]). If using a piecewise linear weight function in all directions (such as on a rectangular patch) with a constant local function, NMM simplifies to the finite volume method [9,10]. If using a constant weight function on an arbitrarily shaped patch with a constant local function, NMM simplifies to the discrete element method. If using a constant weight function on an arbitrarily shaped patch with a linear patch function (resulting in a notable difference from DEM with constant patch function), NMM simplifies to the DDA.
Here we do not attempt to compare all aspects of numerical methods including (1) interpolation/approximation, (2) construction of global equilibrium (transforming differential to integral equations), (3) approaches of integration, and (4) solving of linear or nonlinear global equations. We only compare aspect (1), i.e., approximation/ interpolation as it defines fundamentals of a numerical method. In this respect, NMM provides a flexible and general approach to include continuous and discontinuous methods in a unified form.
In this study, constant patch functions and linear weight functions composed of shape functions of mathematical triangular meshes are used to approximate the physical fields, which are generally expressed as follows: where φ, w, and φ pc are field variables (such as displacements), weight function, and physical patch functions.
The definition of independent mathematical and physical covers lays the basis for NMM to be able to simulate continuum mechanics (deformation) with sufficient accuracy and to simulate discontinuum behavior of blocks that can be flexibly discretized. This dual-mesh concept in combination with the contact algorithms (Section 3) makes it possible for the microscale NMM approach to model both continuum and discontinuum mechanical processes with a broad range of processes. This will be shown in Section 4.

Continuum-discontinuum geomechanics calculation
In order to analyze the microscale mechanics associated with complex geometric features involving both continuous and discontinuous features, it is important to have a rigorous treatment for both continuum mechanics and discontinuum mechanics. Based on the dual-mesh system in NMM and flexible choices of local patch functions, NMM is fundamentally designed for both continuum and discontinuum analyses, including large deformation, large displacements, and multi-body movement and contacts. In a previous paper, the authors described the details of NMM modeling of continuum mechanics with nonlinear features and its full coupling with fluid flow [14]. Here in this paper, we focus on calculations of dynamic contacts among multiple material bodies.
In order to carry out contact calculations, simplifications are often made, for example by using simple shapes such as spheres or rectangles to approximate. The mechanisms involved are shown in Fig. 4, as well as the errors caused by   [32]. b A jump junction for a weak discontinuity [12]. c A discontinuous function for a fracture [13,15] these geometric approximations. For example, it is easy to detect if two spheres are in contact by comparing the distance of their centers with the summation of the radii of the two spheres. But when using a sphere to approximate an angle of a polygon, the sphere cannot capture the case when the angle already enters the boundary of another material body. Such an error leads to an unrealistic series of dynamic contacts that are far from the real world. Similar simplification based on rectangles for contacts is also used, leading to similar types of errors. Another common simplification for contact calculations is to use predefined and fixed contact pairs. With this assumption, the problem is significantly simplified because contacts only occur at those predefined contact pairs. With such simplifications, however, inaccurate calculation of contacts can lead to completely erroneous overall system behavior.
Depending on the states of contact, different constraints may apply, such as Eqs. (6)- (8). The grand challenge of computing contacts is to identify when and where contacts occur among many blocks, which is further complicated by the fact that these blocks are moving, deforming, and in some cases breaking apart. In NMM, a multi-step approach is developed, including detection of dynamic contact pairs, enforcement of contact constraints, and iteration for contact state convergence. Figure 5 shows the schematic of such a contact calculation.

Detection of dynamic contact pairs
First, we detect possible contacting blocks among a number of material bodies, considering they might be moving or deforming, or some of them might be already contacting with other material bodies. This involves a bit of an estimation because the motion and deformation of blocks may occur continuously. Setting a range of possibilities enables precise and complete detection of all possible block contacts. By finding the smallest and largest coordinates in horizontal and vertical directions, all the arbitrarily shaped material bodies can be contained in rectangles (shown in Fig. 5). By comparing the shortest distance between each two rectangles and the distance thresholds of contacts, possible contacts between every two material bodies can be detected.
After detecting every two possible contacting material bodies, it is critical to identify contact pairs, i.e., to calculate where exactly contacts could possibly occur between these two material bodies. All of these possibilities of contact pairs are accounted for in the NMM code. In 2D, contact pairs are categorized into three possibilities: vertex-to-vertex, vertex-to-edge, and edgeto-edge contacts. An edge-to-edge contact is a special case of a vertex-to-edge contact where one edge of the vertex is parallel to the edge of another material body. On the other hand, a vertex-to-vertex contact can be transformed to a vertex-to-edge contact. In Fig. 5, we list all the possibilities of vertex-to-vertex contacts where the thicker lines represent the faces of contacts. In the top row, these include possible contacts of two acute angles, one acute and one obtuse angle, two acute angles with parallel edges. In the bottom row, these include possible contacts of two obtuse angles with parallel edges, two obtuse angels, and one obtuse and one reflex angles. With these possibilities, we are able to consider material bodies in any shapes with any types of corners. Because (1) a vertex-to-vertex contact can be transformed to a contact between a vertex and a contact face, and (2) an edge-to-edge contact is a special case of vertexto-edge contact, in the code we uniformly record each contact pair with the vertex of one material body and the contact face on the other material body.

Enforcement of contact constraints
After contact pairs are identified at each time step, we calculate the normal distance and relative distance along the contact face to pre-estimate the contact state for each contact pair. Then, we enforce the constraint depending on the contact state accordingly. A contact pair may be separated, bonded, or sliding, as described by Eqs. (6)- (8). Because there are no constraints for the separated state, it is naturally calculated by the code.
The distance and displacement constraints for bonded and sliding states can be imposed by using penalty methods [25] or Lagrange multiplier methods [11]. Lagrange multiplier approaches need to introduce additional topology or to relate the Lagrange multiplier to the contact stress. Introducing additional topology associated with Lagrange multipliers is rigorous, but could lead to an ill-conditioned matrix that dramatically slows down the computation [11]. Using contact stress to represent Lagrange multiplier requires the use of currentstep displacements to calculate the stress, which is rigorous but not straightforward. In the current model, we use the penalty method to enforce the constraints associated with distance and displacement. The penalty method is based on the concept of constructing a penalty function g to penalize the deviation from the displacement constraint c(u). The key to effectively using a penalty method is to choose a reasonable value of the penalty parameter. We assume a significantly stiff spring applied to the deviation of a constraint associated with displacements. Therefore, the stiffness of the spring becomes the penalty parameter p. Minimization of deviation of the displacement constraint can be achieved by minimizing the potential energy associated with the work done by the penalty spring. The potential energy Π c by the penalty spring to enforce the constraint can be generally expressed as: The second constraint, i.e., Coulomb's law of friction, is imposed directly by constructing the potential energy: Note that (1) because Eq. (13) includes a dot product, it requires calculation of projecting shearing force on the sliding face, and (2) on the two sides of the sliding face, the absolute shear displacements are different due to sliding (as a relative displacement).
With the above formulae, we are able to consider separated, bonded, and frictional sliding states of contacting material bodies. Dynamic transformation of contact states can be computed by adjusting contact constraints between time steps, or within a single time step before convergence of the iterations is achieved. In some cases, we also need to use criteria described by Eqs. (9)-(10) to judge whether such changes from one contact state to another would happen before switching to the other contact state.

Iteration for contact state convergence
Within each time step, iterations may be carried out several times and pre-estimated contact states may be adjusted until the enforced contacts reach convergence. In a dynamic process, contact pairs may change continuously. For the same contact pair, the three possible states may change dynamically. Thus, for every iteration within a time step, the global equilibrium equations are solved with enforcement of contact constraints from the pre-estimated contact states. At this point, the contact states are re-evaluated to check for consistency with the pre-estimated contact states as shown in Fig. 5. If they are consistent for each contact pair, the calculation proceeds to the next time step. If they are not consistent, the contact states will be adjusted according to the criteria in Eqs. (9)-(10) and the calculation is looped back to re-solve the equilibrium equations until such consistency (i.e., convergence) is achieved.

Examples
The above formulation based on NMM for contact calculation considering large displacements, deformation, and different contact states in dynamic situations was implemented into a new numerical simulator. Here we show a number of examples to demonstrate the capability of the simulator. We first show an example involving rough fractures with explicit geometric representation of their surfaces and asperities with different boundary conditions. Then, we show relocation, compaction, and fracturing of granular systems including clay, loosely packed sand, and cemented marble.

Calculation of rough fractures with explicit geometric representation
In this example, we calculate mechanical compression of rough fractures with explicit representation of their asperities for different degrees of lateral confinement. We first consider a rough fracture with no lateral confinement and then three cases with differing roughness characteristics: (1) The fracture is relatively smooth with small asperities (Fig. 6a); (2) the fracture is relatively rough with larger asperities (Fig. 6b); (3) the fracture has two dominant asperities and a few smaller asperities, leading to an uneven distribution with dominant outliers (Fig. 6c). The sizes of the domains are the same: 10 mm × 5 mm. The meshes that are used for the calculation are shown in Fig. 6a-c. In the first scenario, we apply a vertical loading of 0.42 MPa on the top of the domain. In this scenario, the lower portion of each fracture is fixed, while the upper portion can move freely. We assume a Young's modulus of 400 MPa and a Poisson ratio of 0.3. Figure 6d-f show calculated results of vertical stress. For the first case when the loading is applied, the upper portion of the fracture moves toward and makes contacts with the lower segment on the asperities. As the asperities are relatively small and evenly distributed across the sample, there is not much shearing or sliding on individual asperities. When static equilibrium has been reached, the fracture has closed evenly without much displacement laterally. We can see that high compressive stress concentration occurs at the three contacting areas (Fig. 6d). For the second case in which the fracture has larger asperities, the contact occurs first on the left, leaving an open space on the right side. Then, the entire upper portion slides toward the left as a response to the stress on the right until the entire upper portion contacts with the lower portion and reaches stability, leading to uneven closure and a large offset. Stress concentration is observed at the three contacting areas (Fig. 6e). For the last case, we see that the upper portion first contacts with the lower portion at the dominant asperity on the right, and then the right portion bends as a response to the stress. When the shear stress exceeds the threshold at this dominant asperity, the entire upper part starts to slide along this asperity to the right until it finds the next contacting point as a supporting point to be stabilized. Such a contact point is the top point of the asperity on the right of the lower point (Fig. 6f). As we can see, the concentration of stress is mostly at the dominant asperities. This example involved a very soft system without lateral confinement in order to demonstrate the capability of the simulator to deal with complex micromechanical fracture closure involving large displacement, deformation, and multiple evolving contacts. From this example, we conclude that (1) when single fractures do not have confinement, they tend to move and deform before developing high contacting stress in responding to loading; (2) large asperities dominate contacts and can lead to highly heterogeneous closure of fractures and stress distribution.
In the second scenario, we apply loading on the top of the domain and laterally confine the fracture with stiff columns on each side of a fracture. In this case, we assume a Young's modulus of 4 GPa and a Poisson ratio of 0.3. We assume a Young's modulus for the columns of 40 GPa. Figure 7 shows the vertical and shear stress for each case. We do not observe large displacements in this case because of the confinement, while rather high shear stress is observed at the contacting surface between the fracture and the columns. The average values of closure for these three cases are 0.8 mm, 0.5 mm, and 0.5 mm, respectively. Figures 7 shows that both the vertical stress and the shear stress concentrate at the contacting areas evenly through the fracture for case (1), whereas for cases (2) and (3), the dominant contacting asperities govern the closure as well as stress concentrations. Such a high stress for cases (2) and (3) could lead to a number of responses such as fracture initiation, plastic deformation, or pressure solution if the chemicalmechanical conditions are satisfied.
With this example, we have shown the capability of the NMM code to model the mechanical behavior of fractures with realistic surface geometry at the microscale. Our modeling illustrates that on the local scale, the dominant asperities, which may not be uniformly distributed, govern the closure and stress of a rough fracture. Such asperities may be outliers not captured by the statistical distribution or average roughness.

Settling of aggregates due to a body force
In this example, we calculate settling of suspended clay aggregates induced by a body force. The simulation might be viewed as analogous to dehydration-driven clay shrinkage [29] or to gravitational settling. Considering a hydrostatic situation in which water is drained from the bottom, the clay particles gradually settle due to loss of water. We set up such a problem with the initial structure and geometry shown in Fig. 8a. The domain is 10 μm × 5 μm. The Young's modulus of the clay minerals (suggested by [30]) is 20 GPa and Poisson ratio is 0.3. The mesh that is used for the calculation is shown in Fig. 8a.
Because the particles are not subjected to loading and are loosely packed, we do not see a significant stress developed on the contacts. Such a system forms an unconfined, large-displacement, multi-body system. From Fig. 8 b-c, we see that the system reaches static equilibrium when all the particles settle and the maximum vertical displacement reaches 2.5 μm. We estimate the loss of porosity in this simple case is 60%. In a realistic situation in which the clay is completely dehydrated, depletion of water would cause more complex mechanical-chemical processes (e.g., pressure solution) to occur, rather than resulting in the relatively simple process of mechanical settling. But this example demonstrates the capability for such an application with the key component of modeling multi-body reorganization.

Compaction of a loosely packed sand
In this example, we extract an image of an assembly of grains from a high-porosity sandstone reservoir image [33]. A natural sandstone is cemented along the grain boundaries. However, in order to demonstrate the capability to model grain reorganization, we neglect cementation between grains apparent in the image and calculate compaction of a loosely packed system with drained conditions during the compaction. Our computational domain is 1 mm × 1 mm with left, right, and bottom boundaries fixed and a 1.0 MPa loading on the top. The image and approximated geometric representation of the sand aggregates in NMM with the mesh are shown in Fig. 9. As we can see, such a system consists of grains that have various types of realistic shapes. In our NMM simulation, we consider that all the minerals are 100% quartz with a Young's modulus 50 GPa and a Poisson ratio of 0.3. The friction angle between grains is 30°.
As the system here is rather loosely packed, the particles first move from top to the bottom as a response to loading. The particles on the top contact with neighboring particles stabilize and then the contact force drives the lower particles to move downward until the entire system reaches static equilibrium. Figure 10 shows the horizontal and vertical displacements. Because the loading is vertical, horizontal displacements occur mostly due to rotation of particles induced by contact Fig. 9 a Image of high-porosity sandstone [33]. b Approximated representation of the image in the numerical modeling with the NMM mesh forces. The largest vertical compaction occurs on the right side of the domain because the system is looser on the right and there is more space to accommodate particle movements.
We also show the calculated vertical and shear stresses in Fig. 11. As is apparent, both vertical and shear stresses are the highest at grain contact areas, especially those with sharp corners. Because we use a relatively coarse mesh, we see discontinuity of stress (as a result of linear displacement interpolation) at individual elements, especially when the size of a grain is small relative to its boundary area that is contacting with other grains. This example demonstrates the capability of the model for analyzing a granular system in a dynamic contacting process, and where the shapes of the grains are realistic and can be arbitrary and may contain a number of sharp corners. In this dynamic process, the grains first move with possible large displacement, and then gradually form a tightly packed system while high stress is developed along the contacting boundaries.

Failure of cohesive geomaterial at grain scale
In this example, we extract a picture of an intact marble sample [20] with its pattern of grains (Fig. 12) and simulate the failure of this sample. The simulated domain is 1 mm × 1 mm. We approximate the geometry of all the grain boundaries and generate the NMM mesh shown in Fig. 12, right. In our simulation, we consider that the sample contains 100% calcite with a Young's modulus of 85 GPa, a Poisson ratio of 0.3, and a friction angle of 40°. As marble typically has relatively strong bonds between grains, we set the cohesion as 18 MPa  [20]; right: approximated representation of the image in the numerical modeling with the NMM mesh and tensile strength of grain boundaries as 28 MPa. We simulated vertical loading from the top while the left and right boundaries are free to move and the bottom boundary is fixed. We carried out a number of simulations with two different magnitudes of loading, 11 MPa and 110 MPa, with different scenarios to investigate impacts of cohesion and tensile strength of grain boundaries. Figure 13 shows the simulation results of vertical and shear stress when the loading is 11 MPa without considering cohesion or tensile strength of grain boundaries (Fig. 13a) and considering both effects (Fig. 13b). For the case of no cohesion or tensile strength of grain boundaries, the fractures open easily and slide when contacting forces in the direction of the sliding faces reach the threshold of the Coulomb's law of friction represented by Eq. (8). We can see that a fracture is developed on the left of the sample from the top to the bottom along the grain boundaries (Fig. 13a). On the right of the sample, because of a lack of confinement, the calcite grains break into individual grains along their boundaries, forming a granular system that is free to move. For the case when cohesion and adhesion are applied, we see that the sample remains intact with unbroken grain boundaries (Fig. 13b). Comparing stress for two cases, we find that in the case with no bond strength, a higher vertical stress occurs at the left boundary as a result of the nearby vertical fracture. Figure 14 shows the results of stress when the loading is 110 MPa. As is apparent, the fracture openings for the case with no strength for grain boundaries are larger when the loading is increased (Fig. 14a). For the case with high cohesion and tensile strength of grain boundaries, the sample still breaks as the stress exceeds the cohesive and tensile strengths on the right side, but remains intact on the left side. Except for these two effects, we do not see significant difference in terms of patterns of initiated fracture or stress distribution for this higher loading case.
Using this example, we have demonstrated the capability of the model for handling compressive failure of rocks consisting of strongly cemented grains, from intact to postfailure fragmentation at the microscale. As expected, the intact rock breaks after the stress reaches strengths at the grain boundaries, forming fractures and possibly becoming a granular system when no confinement is present.

Computational efficiency
The simulation time of the code varies depending on the complexity of problems (e.g., the number of contact pairs, the number of contact iterations that is required to reach convergence, and dynamic changes of these factors). For example, for loosely packed granular systems, dynamic changes of contact pairs and contact states may be significant. This requires more loops in the code to detect contact pairs among a number of blocks, and more iterations to reach convergence of contact states within each time step. The simulation time for the examples provided in this paper ranges from minutes to hours on a standard laptop computer. Micromechanical modeling of geomaterials is challenging because of the complex geometry of discontinuities and multiple deformable bodies that contact each other dynamically. In this study, we developed a numerical approach based on the NMM to tackle these challenges. The NMM approach enables us to simulate continuous and discontinuous problems with flexible choices of physical cover functions. In order to calculate dynamic contacts, we consider contacts in separated, bonded, and sliding states that may dynamically change during multi-body movement and deformation. With a multi-step contact calculation involving detection of contacts among multiple blocks with arbitrary shapes, enforcement of different contact constraints for three different contact states, and iteration for contact state convergence, we are able to implicitly calculate dynamic contact behavior rigorously. With these features, our model is capable of simulating geomaterials with arbitrary shapes of grains and interfaces with dynamic contacts at the microscale. We demonstrate the capability with several examples, including a rough fracture with different geometric surface asperity characteristics, settling of clay aggregates, compaction of a loosely packed sand, and failure of an intact marble sample. With our model, we are able to accurately analyze (1) large displacements at a single fracture which is not confined, at loosely packaged granular systems (such as settling of the clay aggregates, the early stage of sand compaction), and at the post-failure stage of geomaterials without confinement (such as the failure of the marble); (2) processes of accumulating high stress at the geomaterial contacting areas, such as at major asperities of a single fracture, and in a confined and tightly packed granular system (the later stages of sand compaction), at grain boundaries of an intact marble before failure; (3) failure of an intact sample due to high stress such as in a mineral cemented marble sample; and (4) postfailure fragmentation such as in the marble sample.
With our simulations, we found the following common features in these different types of mechanical processes of geomaterials: & When geomaterials are not confined, they tend to move and deform before developing high contacting stress in responding to loading. Because of this feature, we can use different types of numerical interpolation for different stages of calculation. For example, when the system is undergoing large displacements of individual grains with dynamic contact alteration, we can use coarser meshes with fewer degrees of freedom to enable efficient calculation of contacts without sacrificing the accuracy, as the stress is relatively small in this case. Once the system reaches a stable contact state and develops high stress, we can use denser meshes to compute the stress accurately. & Large geometric deviations from norms (such as large asperities of a rough fracture and sharp corners in a granular system) dominate contacts and can lead to rather significant heterogeneity of displacement, deformation, and stress distributions in the systems. & For loosely packed granular systems, large displacements may be the dominant response to stress. After large displacements and deformation, such loosely packed systems become tightly packed and high stress can develop at contact areas.
Future work will involve combination of the NMM approach with laboratory experiments for analyzing deformation and contacts of geomaterials with realistic geometric representation, possibly with extension of being coupled with fluid flow, heat transfer, and chemical reaction for multiphysics analysis in complex geosystems at the microscale.