Microscale Precursors of Failure in Fine Grained Porous Sandstone Under Cyclic Thermo-Mechanical Coupling: A Grain Based Model (GBM) Study

The cyclical injection and withdrawal operations for geoenergy storage and/or extraction in subsurface rock induce localised thermo-mechanical (TM) stress fluctuations, which may alter the long-term integrity of the targeted lithologies. Accurate prediction of the behaviour of subsurface rock formations subjected to these stress fluctuations is therefore essential for sustainable energy applications but remains a scientific challenge. This study aims to address this challenge by employing Discrete Element Method (DEM) to simulate the behaviour of rock materials under cyclic TM coupling conditions. A two-dimensional, particle-based numerical model has been constructed to explore microscale precursors, such as mineral deformation, pre-existing microcracks and microstructure damage, that could contribute to rock failure. Particle motion and force-displacement relationships follow linear elastic behaviour. Each particle is connected by cement, with each cement entity represented as a bond (see Fig. 1). These bonds are assigned tensile and shear strengths to resist external stress applied at the particle contact points. If the stress on a bond exceeds its corresponding strength, the bond breaks, triggering a microcrack event. Voronoi tessellation, an effective method for simulating microstructure of rocks, is introduced to mimic the rock mineral shape, with particles filling each tessellated region while maintaining uniform properties (see Fig. 2). First, we calibrate and validate our model using previously published laboratory data (Woodman et al., 2021), achieving a good agreement between the simulated and experimental global response (see Fig. 3). Next, we define microstructural damage as the breaking of bonds, with microcrack events occurring as a direct result during the heating and mechanical cyclic loading stages. Furthermore, to correlate damage with the microstructural evolution, we conduct a spatial analysis of particle rotation angles, porosity distribution, and contact force orientation. Our simulation indicates that the linear contact model seems incapable of reproducing the crack closure stage during compression, the stress-strain curve typically remains linear up to the peak strength. However, introducing particle separation or breaking some bonds before running the model helps produce a gradual, nonlinear compaction (Ji, Zhang and Zhang, 2018). Spatial damage distribution is influenced by Voronoi tessellation geometry because, in GBM, grain boundaries typically have lower strength than grain themselves. As a result, bond breaking occurs preferentially along the grain boundaries. Therefore, the grain size effect needs to be considered to determine a reasonable GBM. Particle rotation evolution before failure can be divided into three stages: initially random directions, followed by local concentrations, and finally aligning along the macroscopic fracture plane.