Simulation of fracture coalescence in granite via the combined finite-discrete element method

Fracture coalescence is a critical phenomenon for creating large fractures from smaller flaws, affecting fracture network flow and seismic energy release potential. In this paper, simulations of fracture coalescence processes in granite specimens with pre-existing cracks are performed. These simulations utilize an in-house implementation of the Combined Finite-Discrete Element method (FDEM) known as the Hybrid Optimization Software Suite (HOSS). The pre-existing cracks within the specimens follow two geometric patterns: 1) a single crack oriented at different angles with respect to the loading direction, and 2) two cracks, where one crack is oriented perpendicular to the loading direction and the other crack is oriented at different angles. The intent of this study is to demonstrate the suitability of FDEM for modeling fracture coalescence processes including: crack initiation and propagation, tensile and shear fracture behavior, and patterns of fracture coalescence. The simulations provide insight into the evolution of fracture tensile and shear fracture behavior as a function of time. The single-crack simulations accurately reproduce experimentally measured peak stresses as a function of crack inclination angle. Both the single- and double-crack simulations exhibit a linear increase in strength with increasing crack angle; the double-crack specimens are systematically weaker than the single-crack specimens.

💡 Research Summary

This paper presents a comprehensive numerical investigation of fracture coalescence in granite specimens containing pre‑existing cracks, using an in‑house implementation of the Combined Finite‑Discrete Element Method (FDEM) called the Hybrid Optimization Software Suite (HOSS). The authors focus on two geometric configurations: (1) a single crack oriented at various angles relative to the loading direction, and (2) a pair of cracks where one is fixed perpendicular to the loading direction while the other varies in inclination. The study aims to demonstrate that FDEM can accurately capture crack initiation, propagation, tensile and shear fracture behavior, and the evolution of coalescence patterns over time.

The methodology section explains how FDEM merges continuum finite‑element analysis with discrete‑element contact detection and interaction. Solid domains are discretized into finite elements, and cohesive interfaces are introduced to model both normal (tensile) and tangential (shear) failure. Tensile failure occurs when the elastic threshold displacement is reached, while shear failure follows a Coulomb‑type criterion τ = c + σ_n tan φ, where c is cohesion and φ is the internal friction angle. Damage is quantified by a scalar D ranging from 0 (undamaged) to 1 (fully fractured), based on the relative displacement at the interface. This framework allows a clear separation between strain‑hardening (continuum response) and strain‑softening (damage localization) regimes.

Material properties are taken directly from the experimental work of Lee and Jeon (2013): Young’s modulus 55 GPa, Poisson’s ratio 0.15, density 2650 kg m⁻³, tensile strength 9.2 MPa, cohesion 55.4 MPa, and friction angle 35°. The authors first validate the model against uniaxial compression (UCS) and Brazilian disk (BD) tests, obtaining simulated UCS of 212 MPa (experimental 209 MPa) and BD tensile strength of 9.9 MPa (experimental 9.2 MPa), confirming the fidelity of the constitutive parameters and numerical implementation.

In the single‑crack simulations, the crack inclination angle α is varied from 0° to 90° in 15° increments. For α = 0°, tensile cracks nucleate at the crack tip and propagate vertically, aligned with the major principal stress. As α increases, the initial tensile cracks still follow the direction of maximum tensile stress, but the proportion of shear damage grows. At α = 45°, the simulation shows an initial tensile wing crack at 2α (≈90°) consistent with the Griffith elliptical model, followed by progressive shear‑tensile mixed damage as loading continues. For α ≥ 60°, shear dominates and the specimen fails primarily in shear. The transition from pure tension to mixed mode is driven by the development of high compressive principal stresses at the crack tips and the frictional interaction (μ = 0.5) between the specimen and loading platens, which promotes shear sliding.

The double‑crack simulations fix one crack perpendicular to the loading direction (α = 0°) and rotate the second crack. The interaction between the two stress concentration zones reduces the overall peak axial stress compared with the single‑crack cases. Although the peak stress still increases linearly with the inclination angle of the rotating crack, the double‑crack specimens are systematically weaker by about 5–10 % across all angles. This result illustrates the weakening effect of intersecting fracture networks, a phenomenon well documented in field observations of rock masses.

Overall, the study demonstrates that FDEM, as implemented in HOSS, can faithfully reproduce experimental fracture coalescence behavior, capture the nuanced transition between tensile and shear failure, and provide detailed temporal evolution of damage fields. The authors argue that this capability makes HOSS a powerful tool for a broad range of geomechanical applications, including block caving, rock blasting, seismic wave propagation, and large‑scale discrete fracture network modeling.