{'fraction_non_alphanumeric': 0.0336105455674799, 'fraction_numerical': 0.015390455157182787, 'mean_word_length': 4.548970799425563, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Geometry evolution of mesoscopic mechanical structures during the rock fragmentation process induced by tunnel boring machine (TBM) cutters

We investigated the geometric evolution of mesoscopic mechanical structures in rocks during the rock fragmentation process induced by tunnel boring machine cutters. Numerical models were built using a grain-based discrete element method to accurately represent the mesoscopic structures and macroscopic mechanical behaviours of rocks, and the relationship between the mesoscopic evolution and the macroscopic response of rock was determined. The major results are as follows. First, the crushing and re-compaction of the grains were found to mainly occur in the thin crushed zone immediately beneath the cutter tip. Second, the reduction in the bearing ability of the dense core during cutter indentation was due to the order increment of the contact topological structure at the mesoscopic scale. Third, the area percentages of low- and high-order meso-loops decreased and increased, respectively, during the indentation process, and the volume expansion of the dense core was mainly caused by an increase in the internal pore area of high-order meso-loops that have low internal solid fractions. Fourth, the low-order meso-loops primarily bore and transferred the indentation force. Finally, the distribution contour of the meso-loops was found to be an appropriate and intuitive approach for representing the evolution of cracks on a macroscopic scale.

We investigated the geometric evolution of mesoscopic mechanical structures in rocks during the rock fragmentation process induced by tunnel boring machine cutters. Numerical models were built using a grain-based discrete element method to accurately represent the mesoscopic structures and macroscopic mechanical behaviours of rocks, and the relationship between the mesoscopic evolution and the macroscopic response of rock was determined. The major results are as follows. First, the crushing and re-compaction of the grains were found to mainly occur in the thin crushed zone immediately beneath the cutter tip. Second, the reduction in the bearing ability of the dense core during cutter indentation was due to the order increment of the contact topological structure at the mesoscopic scale. Third, the area percentages of low-and high-order meso-loops decreased and increased, respectively, during the indentation process, and the volume expansion of the dense core was mainly caused by an increase in the internal pore area of high-order meso-loops that have low internal solid fractions. Fourth, the low-order meso-loops primarily bore and transferred the indentation force. Finally, the distribution contour of the meso-loops was found to be an appropriate and intuitive approach for representing the evolution of cracks on a macroscopic scale.

Introduction
Full-face tunnel boring machines (TBMs) are major high-end equipment for rock tunnel excavation and are especially applied to large mountain tunnel projects and inter-basin water diversion works. The movements of the cutterhead system are described as follows: the cutterhead rotates and advances forward under the joint action of the driving and propulsion systems, resulting in the disc cutters installed on the cutterhead grinding on the tunnel face in concentric circles. Accordingly, the rock fragmentation mechanism can be demonstrated as follows: a high-pressure dense core forms beneath the cutter ring tip as the disc cutter penetrates the rock surface; this drives the propagation of the lateral and median cracks that initiate from the side of the dense core and their connection with cracks produced by neighbouring cutters, generating flat or elongated rock chips.
Investigation of the rock fragmentation mechanism is of vital importance as it directly affects the design of cutter-cutterhead systems and the selection of TBM operation parameters. During the past two decades, owing to the rapid advancement of numerical simulation technologies and computing power, several numerical studies have been conducted to investigate the rock fragmentation process induced by TBM cutters. The representative studies are summarized in table 1. It was found that the most applied numerical methods in these papers are the finite-element method (FEM) and discrete element method (DEM), which account for approximately 40% and 60% of the applied numerical methods, respectively. Focusing on the critical problem of accurate representation of rock material, numerical models have been continuously improved. First, many FEM-based rock-cutting simulations regard the rock as an isotropic material and represent the rock failure process using an element deletion approach controlled by a damage factor assigned to the elastic modulus. Thus, the critical rock re-compaction behaviour in the crushed zone immediately beneath the disc cutter cannot be well represented; as a result, the rock chipping performance and rock cutting forces cannot be reasonably obtained [24]. To overcome these shortcomings, some researchers have explored DEM/FEM coupling [21] or smoothed particle hydrodynamics/FEM coupling [26] approaches by representing a crushed zone with discrete particles and other parts with meshes. Second, compared with the FEMbased numerical approaches, a significant advantage of the DEM-based approaches, for example, the particle flow code (PFC), is that they provide visual representations of crack initiation and propagation. However, their simulation accuracy is usually restrained by the low rolling resistance between neighbouring spherical particles if the linear parallel bond contact model is directly assigned [31]. To solve this problem, some grain-based methods (GBM) [22,27,30,32] have been developed by grouping the particles into breakable clusters or unbreakable clumps according to the size and shape of the mineral grains; hence, the rock chipping behaviour and rock cutting forces can be accurately simulated.
It was also found that these numerical simulations mainly focus on the macroscopic scale, considering rock cutting performance indicators such as rock indentation or cutting forces, rock fragmentation or fracturing performance, macro crack evolution and rock cutting efficiency. This suggests that the rock cutting performance under different conditions considering rock materials, cutter geometries and cutting parameters can be well evaluated, which can help design the TBM cutterhead and select the operation parameters. However, investigations on the mesoscopic mechanical behaviours of the rock fragmentation process induced by TBM cutters are rare. In recent years, research on the multi-scale mechanical behaviour of geotechnical granular materials shows that the macroscopic response of materials under an external load is closely related to the evolution of the particle contact structure at the mesoscopic scale [33]. Herein, the representative studies are listed as follows: Qin & Chi [34] simulated the biaxial shear process of granular material using the software PFC 2D and discussed the relationship between the volume change of void cells and the shape evolutions of void cells; Zhu et al. [35,36] simulated the drained biaxial test of granular material using the DEM software YADE and analysed the critical state fabric in different failure modes on a mesoscopic scale considering force chains and meso-loops; Liu et al. [37][38][39] simulated the biaxial compression and biaxial shear tests of granular material using the software PFC 2D , and analysed the evolution characteristics of meso-mechanical structures such as porosity, coordination number and loop element; Ling et al. [33] simulated the biaxial compression process of soil-rock mixtures under flexible boundaries using the software PFC 2D and investigated the evolution of the micro-contact structure in the internal shear bands and its influence on the macroscopic deformation characteristics. These very interesting studies took advantage of the DEM method in investigating the mechanical behaviours of granular materials at the microscopic and mesoscopic scales, and inspired us to explore the evolution of mesoscopic indicators during the rock fragmentation process induced by TBM cutters, which may provide new insights into the evolution of the crushed zone and dense core.
royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 9: 211630 The major objective of this study is to investigate the geometric evolution of the mesoscopic mechanical structure in the rock fragmentation process induced by TBM cutters. In §2, twodimensional rock indentation models were built considering the constant cross-section (CCS)-and Table 1. Summary of the representative numerical studies on the rock fragmentation process induced by TBM cutters during the past two decades (FEM, finite-element method; DEM, discrete element method; SPH, smoothed particle hydrodynamics; VE-NMM, Voronoi element-based numerical manifold method GBM method, cracks, stress, indentation force; free face conditions Xu et al. [30] royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 9: 211630 V-type disc cutters. As many studies (e.g. [3,4,22]) have verified the appropriateness of the twodimensional equivalent model, disc cutter indentation was considered a plane problem, and tangential rolling was ignored. The rock sample was synthesized using a GB-DEM, where the assembled polycrystalline clusters were breakable from both grain bodies and grain boundaries. In §3, the rock indentation processes were analysed considering the macroscopic rock fragmentation performance and cutter indentation force. In §4, the evolution of the coordination number was first analysed during the entire process of cutter indentation. Then, the organizations among the particle contacts were characterized by the topology of the meso-structures called loops, and the evolution of the number, area and shape of the loops in the crushed zone, macro crack zone and other parts of the rock sample were analysed. Finally, the influence of particle packing parameters on mesoscopic structures was studied.
2. Set-up of the numerical models 2.1. Modelling of the rock sample using the GB-DEM approach The linear parallel bond model written as 'linearpbond' in PFC is a commonly used contact model for rocks, where neighbouring particles are connected and locked with parallel bonds that can be envisioned as elastic springs with limited shear and tensile strengths as well as rotation constraints [40]. However, it has been indicated that, under uniaxial loading, the ratio of compressive strength to tensile strength of synthesized rocks is obviously lower than that of real rocks when the spherical particles are directly connected via parallel bonds [31]. Such a discrepancy may be caused by the lack of interlocking friction along the irregular rock grain boundaries, which are approximated as smooth surfaces, in the spherical particle aggregate [22]. In this study, a modified Voronoi polygon algorithm was used to group the particles into grains to represent polycrystalline rocks and to avoid the uncontrolled randomness of the shape and size of grains generated using primary Voronoi polygons. As depicted in figure 1, the GB-DEM approach for modelling a synthesized rock comprises three steps. During the first step, the rock domain is filled with Voronoi polygons in three sub-steps. In sub-step 1, the rock domain is divided into Delaunay triangles. In sub-step 2, the centres of the triangles' circumferential circles are connected around the triangle vertices. In sub-step 3, the Delaunay triangles are deleted, leaving behind the Voronoi polygons. During the second step, the grains are generated based on the Voronoi polygons created in the first step. In sub-step 4, the rock domain is filled with spherical particles; in sub-step 5, the particles are divided into groups based on the Voronoi polygons, meaning that the particles that fall into the same Voronoi polygon share the same group name. During the third step, the parameters of the intragranular and intergranular contacts are calibrated and assigned. In sub-step 6, the linear parallel bond model is assigned to the contacts between the particles, and the linear model is assigned to the contacts between the particles and walls; in sub-step 7, the particle contact parameters are calibrated and assigned. The particle contact parameters were calibrated according to the following steps. First, a series of uniaxial compression strength (UCS) and Brazilian tensile strength (BTS) tests were conducted on granite collected from the Qinling Mountains, Shaanxi, China, according to the testing method suggested by the International Society of Rock Mechanics (ISRM), to obtain the macroscopic mechanical properties listed in table 2. Second, groups of UCS and BTS tests were simulated for the rock synthesized using the approach described in the preceding paragraph. To ensure that the sample remained in a quasi-static equilibrium throughout the UCS and BTS tests, the loading rate was set to 0.1 m s −1 and the time step was approximately 5.5 × 10 −8 s according to Cho et al. [31]. The calibration procedure was basically a trial-and-error process and the termination condition was that the errors between the Young's modulus, Poisson's ratio, UCS and BTS values of the synthesized rock sample royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 9: 211630 and those of the collected rock specimen were smaller than 5%. If so, the selected microscopic parameters were considered suitable. To reduce the number of independent parameters, several assumptions were made for the calibration procedure according to Li et al. [32]: (i) the density of the particles was determined by the realistic density of the rock material, (ii) the Young's modulus was determined by the effective modulus and the Poisson's ratio was determined by the stiffness ratio, (iii) the effective modulus of the particles and bonds was set the same and so was the normal-to-shear stiffness ratio, and (iv) the packing parameters including the particle radius, particle porosity and grain size were selected based on experience and previous studies [22,31,40]. As noted in the sub-step 5 of the preceding paragraph, the particles falling into the same Voronoi polygons share the same group name. It indicates that the two particles connected by the contact within a grain have the same group name, while the two particles connected by the contact in a grain boundary have different group names. This is the basic logic to distinguish the contacts in the grain boundaries and within the grains when assigning the contact parameters. Using the microscopic properties shown in table 3, the synthesized granite sample presents very similar failure patterns and macroscopic mechanical behaviours to the collected granite specimen (table 2 and figure 2), which indicates that the numerical method is reliable. The granite was crushed so that it could fail along an inclined fracture surface during the UCS experiment and simulation, and it was split into two intact pieces along the centre surface during the BTS experiment and simulation. Furthermore, the curves of the axial stress versus strain during the UCS experiment and simulation were close to each other. Moreover, the errors between the Young's modulus, Poisson's ratio, UCS and BTS of the synthesized granite sample and the collected granite specimen were all lower than 4%. Four types of cracks, that is, intergranular tensile cracks, intragranular tensile cracks, intergranular shear cracks and intragranular shear cracks, are illustrated in cyan, red, green and blue colours, respectively. During the UCS test simulation, the tensile cracks and intergranular tensile cracks accounted for approximately 84% and 61% of all the cracks, respectively. During the BTS test simulation, the tensile cracks and intergranular tensile cracks accounted for approximately 86% and 56% of all the cracks, respectively. This means that the major rock failure pattern is tensile failure of the contacts between neighbouring particles in neighbouring grains.

Rock indentation numerical models
The rock indentation numerical model is shown in figure 3. The tip width of the CCS-type cutter was 20 mm, the wedge angle of the V-type cutter was 100°, and the cutter spacing was 70 mm for both models. The length and width of the rock samples were 500 mm and 200 mm, respectively. The rock sample was composed of 113 055 particles that were grouped into 31 352 grains. The microscopic properties have been calibrated in §2.1 and are shown in table 3. The cutter indentation speed was set as 0.2 m s −1 (time step is 6.0 × 10 −8 s) to ensure the quasi-static equilibrium of the rock sample. The initial distance between the ring tip and the top surface of the rock sample was set as 0.3 mm. The indentation process was terminated when the indentation force dropped to approximately 0.

Macroscopic rock indentation performance
The macroscopic rock indentation performances of the CCS-and V-type cutters are shown in figure 3. According to the curves of the cutter indentation force (F N ) versus indentation depth (D), the Table 2. Mechanical properties of the synthesized rock samples for simulation and the collected rock specimen for experiment (note: sample-0 is for §3 and § §4.1-4.6, and sample-1 and -2 are for §4.7 of sensitivity analysis indentation process includes two stages: loading and unloading. Note that D equals (D N + 0.3) and D N is the net indentation depth. As shown in figure 3a,b, during the loading stage, the F N of CCS-type cutter 1 and 2 increased linearly to their peak values of 1.95 × 10 7 and 2.32 × 10 7 N m −1 , respectively, at D of 3.0 and 3.7 mm, respectively. During the unloading stage, the F N of the two cutters dropped rapidly to approximately 0 when D is 4.3 mm, which is called the critical indentation depth (D C ). There were many small fluctuations in both the loading and unloading stages of the curves. The change in F N with D was consistent with the jump crashing of rocks [12,14,41]. The rock beneath the cutters was compressed into many small fragments in the loading stage, and these fragments acted as a dense core to further bear and transfer the indentation force into the rock and drive the macro cracks to propagate. These results are similar to those of previous simulations [1,27,28]. As shown in figure 3a,b, the rock fragmentation results at D of 3.8 and 4.3 mm, respectively, are similar, indicating that the macro cracks could no longer propagate once the indentation entered the unloading stage. This is because the dense core collapsed during the unloading stage, resulting in stress release and rapid reduction of the indentation force. This is further evidenced by the distribution of the force chains, in which two inverted triangles appear beneath the cutter tips inside of which the force chains are almost zero at the D C . Macro cracks generally developed along the indentation direction, and the  royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 9: 211630 propagation of side cracks was inadequate. This is because the studied rock was a hard granite (UCS ≈ 180 MPa), and it was not confined. This agrees with the results of previous studies [1,15,27] that the rock failure process changes from the formation of rock chips to a vertical axial splitting failure with decreasing confining pressure. The cutter indentation specific energy (SE) calculated by dividing the work of the cutter indentation force by the total fragment area was used to evaluate the rock fragmentation efficiency [14]. The average SE of the two CCS-type cutters was calculated as 4.8 MJ m −3 . As shown in figure 3c,d, the F N of V-type cutter 1 and 2 increased linearly to their peak values of 1.00 × 10 7 and 1.14 × 10 7 N m −1 , respectively, at D of 2.0 and 2.3 mm, respectively, and D C is 2.7 mm. The dense cores were small. The main cracks were median cracks along the indentation direction, and the propagation of the side cracks was inadequate. Thus, no large chips were produced between the two V-type cutters. The average indentation SE of the two V-type cutters was calculated as 22.5 MJ m −3 . Compared with that of the CCS model, the rock fragmentation of the V model was less efficient, even though the peak F N and D C of the V-type cutters were much smaller. This is because the rock fragments produced by the V-type cutters were fewer than those generated with the CCS-type cutters. In addition, the median crack length of the V model was shorter than that of the CCS model. These differences may be because the dense cores collapse more easily under the wedge indenter than under the flat intender.

Evolution of the mesoscopic indicators during the cutter indentation process 4.1. Coordination number
The coordination number (Z ) represents the average number of particles in contact with each particle. This is an important factor used to evaluate the contact status, stress state and compactness of the PFC2D 6.00 PFC2D 6.00 PFC2D 6.00 royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 9: 211630 particle system. Z can be calculated as follows [42]: where N C and N P are the number of actual 'ball-ball' contacts and the total number of particles, respectively.
To illustrate the evolution of the coordination number of the rock sample during the cutter indentation process, four states were selected for both the CCS model and the V model. According to figure 3, in state 'A', the initial state, the cutter has not penetrated the rock. State 'B' represents the early loading stage where the net indentation depth (D N ) is small, state 'C' represents the later loading stage where the indentation force (F N ) almost reaches its peak value, and state 'D' is the final state, where D N equals to D C , and F N drops to 0. The distribution of the contours of Z is shown in figure 4.
The distribution of Z in the initial state 'A' changed in different areas of the rock sample. This was because the rock was synthesized by randomly distributed particles whose radii varied between 0.4 and 0.6 mm, which represents the anisotropy and heterogeneity of the rock constituents at the microscopic scale. As the cutter penetrated the rock, the Z of the rock area away from the cutters increased gradually, and the distribution became uniform. For the CCS model in state 'B', the Z in the two elongated areas beneath the cutter was lower than that in the other areas. In addition, the Z of the two thin areas immediately beneath the cutter tips was especially low. This indicates that the thin rock areas were crushed, and slight dislocations occurred in the particles in the elongated areas. Then, when the cutter D N reached 2.7 mm in state 'C', the Z of the thin areas became a little higher than that in state 'B', indicating the re-compaction of the crushed grains. In addition, the lengths of the elongated areas with low Z values were longer than those in state 'B', which is in accordance with the propagation of median cracks shown in figure 3. Finally, the model was in state 'D', where D N is 4.0 mm. In this state, the Z of the rock areas beneath the cutter tip was the lowest compared with other states. The areas with low Z values coincided with the two inverted triangles with small contact force chains, as shown in figure 3b. Moreover, the Z value of the median crack areas became larger in    state 'D' than that in state 'C'. This illustrates the closure of the median cracks due to the rebound deformation of the bottom rock area after the collapse of the rock areas beneath the cutter tip. For the V model, a similar progressive evolution of Z could be observed during the cutter indentation. However, the fragmented areas were smaller, and the lengths of the median crack areas were shorter than those in the CCS model. Accordingly, the rock beneath the cutter tip is divided into three areas, namely, the crushed zone, dense core and median crack area, as illustrated in figure 4a. The crushed zone is the thin area immediately beneath the cutter tip, where the grains are crushed and then re-compacted; the dense core represents the triangular area that bears the indentation force in the loading stage and then collapses in the unloading stage, and the median crack area represents the elongated area along the median crack that initiates from the bottom of the dense core. For further detailed analyses, the Z values of the crushed zones and dense cores of the CCS and V models were calculated, as shown in figure 5. The crushed zones of the CCS model are defined as two rectangles with widths and heights of 20 and 10 mm, respectively. For the V model, the width and height of the crushed zone rectangles were 12 and 4 mm, respectively. The dense cores of the CCS model are defined as two equilateral triangles with a bottom and height of 70 and 30 mm, respectively. For the V model, the bottom and height of the dense core triangles were 40 and 14 mm, respectively. As shown in figure 5, the variations in the Z values for the two crushed zones or dense cores beneath the same type of cutters are similar to each other, indicating a similar rock indentation performance of the neighbouring cutters. For both the crushed zones and the dense cores, the Z of the CCS model was slightly larger than that of the V model initially when D N was 0. This is caused by the difference in the sizes of the crushed zones and dense cores for the CCS and V models. As shown in figure 5a, the Z values of the crushed zones of the CCS model first decrease linearly, then tend to increase when D N is larger than 1.7 mm, and finally decreases when D N is larger than 3.2 mm. This variation in Z intuitively confirms the crushing and re-compaction of the grains in the crushed zone, which have been mentioned based on the distribution contours of Z in figure 4. Similar crushing and re-compaction phenomena can be found for the V model, but come earlier than the CCS model. As shown in figure 5b, the Z of the dense cores of the CCS model first decreases slowly and then decreases rapidly when the D N is larger than 3.2 mm. This variation in Z intuitively represents the gradual crushing and sudden collapse of the dense core during the loading and unloading stages, respectively. A similar variation can be found for the V model, but the sudden drop occurs earlier when the D N is approximately 1.8 mm, indicating that the dense core collapses more easily when indented by a V-type cutter than a CCS-type cutter. Comparing figure 5a,b, it can be summarized that (i) the re-compaction phenomenon only occurs in the thin crushed zone and (ii) the Z of the dense core is larger than that of the crushed zone.

Evolution of the meso-loop distribution
Meso-loops are basic structures that conduct the mechanical behaviour of granular materials at the mesoscale. As shown in the figure 6c, the loops are defined as enclosed polygons by contact branches; hence, The initial meso-loops were mainly of a low order (L 3 and L 4 ), and only a few were of the order of L 5 . As the cutter penetrated the rocks, higher-order meso-loops formed gradually in the 'dense core' domains and the 'median crack' domains. During the loading stage, as illustrated in figure 7b-d and figure 8b-d, the positions of the high-order loops in the dense cores were basically in coincidence with the macro cracks that distribute radially, and the shapes of the high-order loops were more elongated than those of the low-order loops. This indicates the separation of the dense core into several large fragments along the radial cracks, and results in a gradual reduction in the bearing ability of the dense core. During the unloading stage (figure 7e,f and figure 8e,f ), the number of high-order loops increased significantly in the dense cores, especially in the area directly beneath the cutter tip. This indicates that crushing and collapse of the dense core and rearrangement of the contact structure occurred. Therefore, as the internal solid fraction of the meso-loop tends to decrease, and more freedom is introduced into the loop structure to make it more deformable [35], as the loop order increases, the evolution of the bearing ability of the dense core is attributed to the order increment of the contact topological structure at the mesoscopic scale.

Evolution of the number proportions of meso-loops of different orders
The evolution of the number proportions of meso-loops of different orders is the most intuitive quantitative reflection of the evolution of the meso-structure [39]. The changing curves of the number proportions of the different-order meso-loops of the CCS model are shown in figure 9. The meso-  royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 9: 211630 loops were mostly of a low order initially when D N was zero. The number proportions of the L 3 -and L 4order meso-loops were approximately 66-67% and 27-29%, respectively, for all the studied domains. As shown in figure 9c, for the 'none crack' domains, the number proportions of meso-loops of different orders were almost constant during the indentation process, indicating that the rock meso-structure did not change. As shown in figure 9b, for the 'median crack' domains, the number proportions of low-order meso-loops slightly decreased, while the number proportions of high-order meso-loops slightly increased. For example, the number proportions of the L 3 -and L 4 -order meso-loops decreased by 1.4% and 0.2%, respectively, and the number proportions of the L 6 -and L 7 -order meso-loops increased from 0% to 0.37% and 0.24%, respectively. This is because the evolution of the rock mesostructure only occurred in the small domains along the median cracks. As shown in figure 9a, for the 'dense core' domains, the number proportions of meso-loops of different orders change significantly. The number proportions of different-order meso-loops in the initial stage and at the ends of the loading and unloading stages are listed in table 4, as are the calculated increments. During the loading stage, the proportion of the L 3 -order meso-loop decreased from 66.8% to 59.9%, the number proportion of the L 4 -order meso-loop increased from 28.3% to 31.1%, the number proportion of the L 5 -order meso-loop increased from 4.8% to 5.4%, and the number proportions of the L 6 -, L 7 -, L 8 -, L 9and L 10+ -order meso-loops increased from 0 or approximately 0 to 2.0%, 0.4%, 0.6%, 0.2% and 0.4%, respectively. This illustrates that only the most stable L 3 -order meso-loop decreased, while meso-loops of other orders increased during the loading stage. Meanwhile, the evolution of the number proportions of meso-loops in different orders was gradual and slow during the loading stage. By contrast, the evolution accelerated once the indentation entered the unloading stage, and obvious changes in the number proportions occurred in the short unloading stage, especially for the highorder meso-loops. For the L 4 -order meso-loop, the number proportion slightly decreased back to approximately the initial value during the unloading stage. This phenomenon reflects the destruction and recombination of mesostructures during the indentation process. The changing curves of the proportions of the different-order meso-loops of the V model are shown in figure 10. The general change rules were similar to those of the CCS model, but the detailed change percentages were different. During the entire indentation process for the 'dense core' domains, the number proportions of the L 3 -and L 4 -order meso-loop decreased by 19.6% and 9.7%, respectively, while the proportion of the L 5 -order meso-loop increased by 49.2%. royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 9: 211630 The results of the above analyses reveal that the evolution of the meso-structures of the rock sample during indentation mainly occur in the 'dense core' domains; hence, further investigation of the mesoloop evolution will be conducted by mainly considering these domains. For the CCS and V models, the proportion of low-order meso-loops (L 3 and L 4 ) decreased from 100% to 76.9% and 78.2%, respectively, during the entire indentation process, indicating that the general evolution of the mesostructures was the same when the rock was indented by CCS-or V-type disc cutters. The major royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 9: 211630 evolutions of the number proportions were a decrease in the L 3 -and L 4 -order meso-loops and an increase in the L 5 -and L 6 -order meso-loops, which mainly occurred in the unloading stage.

Evolution of the area percentages and average areas of meso-loops of different orders
In addition to the number proportions of meso-loops of different orders, the evolution of the area percentages and average areas of meso-loops of different orders are also needed in order to represent the evolution of the meso-structure. The evolutions of the area percentage of meso-loops in the 'dense core' domains for the CCS and V models are illustrated in figures 11a and 12a, the evolution of the average area of the meso-loops are illustrated in figures 11b and 12b, and the area percentages and increments of different meso-loops in the 'dense core' domains in different indentation stages are shown in table 5. The L 10+ -order meso-loops were not analysed because the variability of the area of the L 10+ -order meso-loops was too high. For the CCS model, the general change rules of the area percentages were similar to those of the number proportions of different-order meso-loops. The initial area percentages of the L 3 -and L 4 -order meso-loops were approximately 46% and 42%, respectively. During the loading stage, the area percentages of the low-order meso-loops (L 3 and L 4 ) decreased steadily and slowly to a sum of 77.8%, while the area percentages of the high-order meso-loops (L 5 , L 6 , L 7 , L 8 and L 9 ) increased steadily and slowly to a sum of 22.2%. By contrast, during the unloading stage, the area percentages of the low-order meso-loops (L 3 and L 4 ) decreased significantly to a sum of 51.4%, while the area percentages of the high-order meso-loops (L 5 , L 6 , L 7 , L 8 and L 9 ) increased significantly to a sum of 48.6%. As shown in figure 11b, the average areas of the low-order meso-loops are generally constant during the indentation process, while the average areas of the high-order meso-loops (especially L 9 ) improve as the cutters penetrate the rock, and the higher the meso-loop order, the more sensitive the average area to the indentation depth. In addition, as discussed in §4.3, the number proportions of the low-and high-order meso-loops decreased and increased, respectively, during the indentation process. These are the reasons for the decrease and increase in the area percentages of the low-and high-order meso-loops, respectively. As the average meso-loop area is a rough representation of the material porosity, the decrease in the number proportions and area percentages of low-order meso-loops and the increment of the number proportions and area percentages of high-order meso-loops result in the volume expansion and collapse of the dense core. Hence, the volume expansion of the dense core was not significant during the loading stage. In addition, it has been discussed in §3 that the propagation of macroscopic cracks mainly occurs in the loading stage. As a result, it may be inferred that the expansion of the dense core is not the principal cause of macroscopic crack propagation. For the V model, the breakover of the area percentage curves at the end of the loading stage was not as remarkable as that for the CCS model. During the entire indentation process, the area percentages of the low-order (L 3 and L 4 ) meso-loops continued to decrease from a sum of 100% to 55.0%. Table 4. Number proportions and increment of different-order loops in the 'dense core' domains in different indentation stages ('L_end' and 'U_end' mean the end of loading and unloading stage respectively, and 'INCR (%)' means the increment of the number proportion compared with the initial stage). 4.5. Evolution of the internal pore area and porosity of meso-loops of different orders As shown in figure 13, the area of a meso-loop (S L ) is composed of two parts: the covered particle area (S pi ) and the internal pore area (S V ). As the particles are usually assumed to be undeformable, it may be   royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 9: 211630 more rational and accurate to attribute the volume change of the 'dense core' domains to the evolution of the internal pore area of the meso-loops. According to Ling et al. [33], the internal pore area of a mesoloop was calculated using equation (4.2).
where S V denotes the area of the internal pore; S L denotes the area of a meso-loop that has been calculated in §4.4, S pi denotes the area of a particle region that is covered by the meso-loop and S oi denotes the overlapping area of neighbouring particles. Hence, the porosity of a meso-loop (PL) was calculated using equation (4.3).
As noted in §4.4, the internal pores of the L 3 -to L 9 -order meso-loops were studied, while the L 10order meso-loops were not considered. For the CCS model, the evolution of the average S V and P L of meso-loops of different orders in the 'dense core' domains during indentation are shown in figure 14.
The results revealed that the change rules of the average S V were similar to those of the average loop area illustrated in figure 11b, as the particles were undeformable. First, the average S V increased with increasing loop order, which agrees with the fact that the internal solid fraction of the meso-loop tends to decrease as the loop order increases. Furthermore, the average S V of the L 3 -, L 4 -and L 5 -order low: high: low: high: Figure 12. Evolution of the area percentage (a) and average area (b) of different meso-loops in the 'dense core' domains during cutter indentation process of the V model. meso-loops were constant during cutter indentation, while the average S V of the high-order meso-loops (especially L 9 ) improved as the cutters penetrated, and the higher the meso-loop order, the larger the increment of the average S V during the entire indentation process. This indicates that the volume expansion of the dense core is mainly caused by the increment in the internal pore area of the highorder meso-loops. The average P L increased with an increase in the loop order when the order was lower than L 6 , and the average P L of the L 7 -, L 8 -and L 9 -order meso-loops were not significantly different. During the loading stage, the average P L of different-order meso-loops generally decreased with an increase in D N ; subsequently, the average P L tended to increase with increasing D N during the unloading stage. This change rule was especially evident for the L 6 -and higher-order loops. This indicates that the meso-loops were compressed during the loading stage and then expanded during the unloading stage.
In the case of the V model, as shown in figure 15, the change rules of the average S V with D N were similar to those of the average loop area illustrated in figure 12b. First, the average S V increased with the loop order. Then, the average S V of the L 3 -to L 6 -order meso-loops did not change with D N , while the average S V of the L 6+ -order meso-loops increased significantly with an increase in D N . This means that the volume expansion of the dense core was mainly caused by the increment in the internal pore area of high-order meso-loops that had low internal solid fractions. The average P L of the L 3 -to L 6 -order meso-loops increased with the loop order, and the average P L of the L 7 -, L 8 -and L 9 -order meso-loops did not show a clear change rule with the loop order. The average P L of the L 3 -, L 4 -and L 5 -order meso-loops was generally constant during the indentation process, while the average P L of the L 6+ -order meso-loops generally increased with an increase in D N . This indicates that no obvious volume compression of the meso-loops occurred during the loading stage when the rock was indented by a V cutter.

Evolution of the shape of meso-loops of different orders
The shape of a meso-loop was evaluated using the elongation degree (β m ) calculated from the tensor (L m ) of the meso-loop [43]. The loop tensor L m is represented by formulae (4.4) and (4.5), and its trace measures the meso-loop perimeter, and the elongation degree (β m ) is calculated by formula (4.6). The components, l k and n b k , are shown in figure 16a. Figure 13. Sketch of the calculation of the internal pore area of a meso-loop.
royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 9: 211630 where l k is the branch vector at contact k; l k ¼ l k n k b , with length l k and unit vector n k b directed along branch k; r m is the number of sides of meso-loop m; and D m is the deviatoric part of loop tensor L m which characterizes the shape of the meso-loop under consideration. L m 1 and L m 2 are the major and minor eigenvalues of L m , respectively.
As shown in figure 16b, the elongation degree, β m , was zero for a meso-loop of regular shape, and the higher the elongation degree, β m , the greater the elongation of the meso-loop [43]. The β m of the mesoloops in the 'dense core' domains of the CCS model when D N was 4.0 mm are shown in figure 17. It can be seen that the number of meso-loops decreased with increasing meso-loop order, which has been discussed in detail in §4.3. Furthermore, β m generally increased with increasing loop order, indicating that the higher the loop order, the more elongated the meso-loop. This is because the high-order meso-loops mainly occur along macroscopic cracks. Moreover, it was found that the fluctuation amplitude of the L 3 -order meso-loop was much smaller than that of the other meso-loops, as the L 3 -order meso-loop was the most stable and usually had a regular shape.
To study the evolution of the shapes of different-order meso-loops during cutter indentation, the average values of β m of different-order meso-loops in the 'dense core' domains for the CCS and V models were calculated and are shown in figure 18. The results confirm the above finding that β m generally increases with increasing loop order. The average β m of the L 3 -order meso-loop was much low: high: low: high: Figure 15. The average internal pore area (S V ) (a) and porosity (P L ) (b) of meso-loops in the 'dense core' domains during cutter indentation process of the V model.
royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 9: 211630 smaller than those of the other meso-loops, and it was constant during the indentation process. The average β m of the L 4 -and L 5 -order meso-loops were also basically constant during cutter indentation, but they were evidently larger than that of the L 3 -order meso-loop. As meso-loops are the basic structure responsible for the mechanical behaviour of granular material, the constant shape and dominant proportion of low-order meso-loops during indentation indicate that the L 3 -, L 4 -and L 5 -order meso-loops are the principal meso-loops that bear and transfer the indentation force. For the CCS model, the average β m of the L 6 -order meso-loops increased from the initial value of 0.25 to 0.31 when D N was 0.6 mm. Subsequently, it continued to decrease to 0.25 as the cutters continued to penetrate the rock and as D N increased to 1.2 mm; then, it tended to increase to approximately 0.32 during the following loading stages. As the L 6 -order meso-loop is a link between the low-and highorder meso-loops, the evolution of its shape is discussed to explain the initiation and propagation of cracks in the 'dense core' domains at the mesoscopic scale. In general, the entire crack evolution process in the 'dense core' domains can be divided into three stages based on the indentation depths of 0-0.6, 0.6-1.2 and 1.2-3.0 mm, respectively, according to the change rule of the average β m of the L 6 -order meso-loops.
As shown in figure 19, during the first stage of local crack initiation, two dense and strong contact force regions with an outline similar to the Hertzian stress ellipse are formed directly beneath the cutter tips. Because of the displacement inconformity of the particles in the boundary of the strong contact force regions, the bonds between neighbouring particles are broken and illustrated in 'fractures'. The connections and constraints are thus reduced, and hence dislocation and shear dilation occur as the particles are driven to move by the contact force chains, causing the development of L 3 -, L 4 -and L 5 -order meso-loops to L 6 -and higher-order meso-loops (see figure 7b). As the newborn L 6 -order meso-loop is usually narrow and elongated in shape, the average β m increases accordingly. It has been noted previously that L 6 -and higher-order meso-loops are more deformable than the L 3 -, L 4and L 5 -order meso-loops as more particle freedom is introduced. Therefore, the scattered distributed high-order meso-loops can be regarded as the initiation of local cracks in the 'dense core' domains.
The second stage mainly includes the processes of local crack mergence and macro crack formation. As the cutters penetrates the rock, the elongated L 6 -order meso-loops are compressed to expand laterally and hence become more regular. Meanwhile, the low-order meso-loops around the local cracks are forced to develop into high-order meso-loops, causing the merging of neighbouring local cracks and the royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 9: 211630 formation of macro cracks, as illustrated in figure 7c. The third stage mainly includes the processes of macro crack propagation and local crack generation. As the macro cracks propagate, the 'dense core' domains are broken into several fragments (see figure 7d) and the force bearing ability reduces gradually, causing the generation of more elongated L 6 -and higher-order meso-loops from low-order meso-loops, which can explain the increase in the average β m . For the L 7+ -order meso-loops, the change rule of the average β m with D N was not evident, as there were only a few high-order mesoloops, and their shapes varied significantly during cutter indentation. For the V model, the average β m of the L 3 -, L 4 -and L 5 -order meso-loops were generally constant during the indentation process, and the average β m of the L 6 and higher-order meso-loops were higher than those of the low-order meso-loops; these results are similar to those of the CCS model. The variation in the average β m of the L 6 -order meso-loops with the D N of the V model was similar to that of the CCS model. The main difference existed at the beginning, as there were no initial L 6 -order meso-loops in the studied 'dense core' domains for the V model. The subsequent decrement and then the increment of the average β m of the L 6 -and higher-order meso-loops can be explained by the following assumptions: (i) local crack mergence and macro crack formation, and (ii) macro crack propagation and more local crack generation. The contact force regions and the initiation of the 'fractures' in the V model are illustrated in figure 19 when D N is 0.6 mm. As the contact area between the sharp V-type cutter and rock was very small, the contact force regions were small, and particle displacement mainly occurred. Even so, it confirms the conclusion drawn from the CCS model analysis that the initial 'fractures' come into being due to the displacement inconformity of the particles in the contact force region boundary.
For both the CCS and V models, the average β m of the L 7+ -order meso-loops decreased to low values once the indentation entered the unloading stage. This indicates that the shape of the L 7+ -order mesoloops becomes regular during the collapse of the 'dense core' domains, as the meso-loops evolve freely without the strong contact force chains, as shown in figure 3b,d. This result is in accordance with the results discussed in § §4.4 and 4.5 on the evolution of the average area of the meso-loop and the internal pore. For a more detailed analysis, the meso-loop can be simplified as an ellipse that is fitted by the vertices of the loop, according to Liu et al. [37][38][39]. The ratio i of the major axis length to the minor axis length can be regarded as an indicator of the meso-loop shape, and i is always not less than 1. The larger the value of i, the more elongated the meso-loop, which is similar to the trend for indicator β m . The area of an ellipse SE can be calculated using equation (4.7) once the ellipse perimeter L and ratio i are known. For a certain meso-loop simplified by an ellipse, the perimeter L is essentially constant during the shape evolution. Thus, it is easy to determine that the SE decreases with increasing i. Therefore, the results that (i) the average area of high-order meso-loops increases with D N , as illustrated in figures 11b and 12b, and that (ii) the average β m of the high-order mesoloops decreases with D N , as illustrated in figure 18, can confirm each other. royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 9: 211630 Based on the above analyses, some additional summaries and assumptions are discussed. For simulations using the PFC software, the macroscopic crack is usually represented by the built-in 'fracture', which is inserted in the contact position once the contact bond is broken. However, it might be inappropriate for representing macro cracks in the 'dense core' domains during cutter indentation. This is because the breakage of a contact bond is not equivalent to the disappearance of the contact itself. For example, as shown in figure 7b and figure 17a, many 'fractures' come into being during the early loading stage; however, most particles are compacted, and few macro cracks appear owing to the constraint of the cutter tip. By contrast, the distribution contour of the meso-loops that are coloured according to the loop order may be more appropriate and intuitive for representing the evolution of cracks on a macroscopic scale. The initiation and propagation of cracks in the 'dense core' domains are attributed to the combination of external and internal causes. The main external causes are the contact force chains and the boundary conditions, while the main internal cause is the evolution of mesoscopic contact structures, which are also called meso-loops. Therefore, the reason for crack growth on the mesoscopic scale is assumed to be the evolution of low-order meso-loops to highorder meso-loops driven by contact force chains.

Influence of particle packing parameters on mesoscopic structures
In the DEM, the mesoscopic mechanical structures, especially the coordination number and meso-loops, are strongly influenced by the factors of particle packing, including particle size or grain size, particle size distribution, porosity, etc. According to previous publications [22,31,32,40], there are some commonly used ranges or assumptions for the packing parameters of rock material aggregated by particles: (i) the particle radius is usually 0.1-0.9 mm considering the model accuracy and simulation cost, (ii) the ratio of maximum particle radius to minimum particle radius is usually 1.5-2.0, (iii) the porosity of particle-synthesized rock is usually 0.1-0.15 for two-dimensional models, and (ivd) the grain size is usually adjusted by the quantity of grains according to laboratory scan. Based on these assumptions, another two sets of particle packing parameters were selected different from those of the existing model (sample-0) to synthesize two new rock samples (sample-1 and -2), and the particle contact parameters were re-calibrated and were different from those of sample-0, as shown in tables 2 and 3. It indicates that the macroscopic mechanical properties of the rock sample are sensitive to particle packing and contact parameters.
To study the influence of particle packing parameters on mesoscopic structures, another two simulations of the rock indentation process by CCS cutters were conducted. For the 'dense core' domains of the three simulations, the coordination number (Z ) and number proportions of the mesoloops with different orders were selected as representative mesoscopic structure indicators and were analysed. As shown in figure 20, the initial Z values of the three models are different, indicating that the coordination number is sensitive to the particle packing parameters. However, the general trends of the evolution of the Z value for the three models are like each other. This variation in Z intuitively represents the gradual crushing and sudden collapse of the dense core during the loading and unloading stages, respectively. As shown in figure 21 and figure 9a, the initial number proportions of the meso-loops with different orders for the three simulations are different, indicating that the meso- 1.1000 × 10 +6 1.0000 × 10 +6 9.0000 × 10 +5 8.0000 × 10 +5 7.0000 × 10 +5 6.0000 × 10 +5 5.0000 × 10 +5 4.0000 × 10 +5 3.0000 × 10 +5 2.0000 × 10 +5 1.0000 × 10 +5 0.0000 × 10 PFC2D 6.00 (b) (a) Figure 19. Distribution of the ball displacement, contact force and fractures beneath the cutter tip when D N is 0.5 mm for the CCS model (a), and D N is 0.6 mm for the V model (b). ©2020 Itasca consulting group, Inc. royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 9: 211630 loop distribution is sensitive to the particle packing parameters. However, the L 3 -, L 4 -and L 5 -order mesoloops account for the most proportion for the three simulations. Besides, the general trends of the evolution of the number proportions of different-order meso-loops are similar for the three simulations. It reflects the destruction and recombination of mesostructures during the indentation process in the 'dense core' domains.
According to the above analysis, it is found that the mesoscopic mechanical structures are sensitive to the particle packing parameters. However, the general evolution trends of the mesoscopic structures are like each other for the simulations conducted on the rock samples with different particle packing parameters. It indicates that the conclusions drawn in this study are universal to some extent.

Conclusion
This study mainly investigated the geometric evolution of mesoscopic mechanical structures in rock fragmentation processes induced by TBM cutters. The rock indentation models were built using a GB-DEM considering CCS-and V-type cutters. The evolution of the mesoscopic indicators, including the coordination number, meso-loop distribution, number, area, porosity and shape, were investigated,  royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 9: 211630 and the relationship between the mesoscopic evolution and the macroscopic response of rock is discussed. The following conclusions can be drawn: (1) Macro cracks mainly grow in the loading stage as the collapse of the dense core in the unloading stage results in stress release and thus a rapid reduction in the indentation force. The rockbreaking efficiency evaluated by the indentation specific energy is much higher for the CCS model than for the V model. (2) According to the evolution of the coordination number, crushing and re-compaction of the mineral grains mainly occur in the thin crushed zone immediately beneath the cutter tip, while the 'dense core' is gradually crushed and then suddenly collapses during the loading and unloading stages respectively. (3) According to the evolution of the meso-loop distribution, it can be inferred that the reduction in the bearing ability of the dense core during cutter indentation is attributed to the order increment of the contact topological structure at the mesoscopic scale. (4) The evolution of the meso-structures during indentation mainly occurs in the 'dense core' domains.
The major evolutions of the number proportions are a decrease in L 3 -and L 4 -order meso-loops and an increase in L 5 -and L 6 -order meso-loops, which mainly occurs in the unloading stage. (5) The area percentages of the low-and high-order meso-loops decrease and increase, respectively, during the indentation process. The expansion of the dense core is inferred not to be the principal cause of macroscopic crack propagation. The volume expansion of the dense core is mainly caused by the increase in the internal pore area of high-order meso-loops with low internal solid fractions. (6) According to the shape analysis of the meso-loops, the higher the loop order, the greater the elongation of the meso-loop during cutter indentation. The constant shape and dominant proportion of the loworder meso-loops during cutter indentation indicate that the L 3 -, L 4 -and L 5 -order meso-loops are the principal meso-loops that bear and transfer the indentation force. The shape of the high-order mesoloops becomes regular as the cutter penetrates during the unloading stage. (7) The distribution contour of the meso-loops that are coloured according to the loop order is suggested to be an appropriate and intuitive approach for representing the evolution of cracks on a macroscopic scale. The reason for crack growth on the mesoscopic scale is assumed to be the evolution of loworder meso-loops to high-order meso-loops driven by contact force chains.
The main focus of this study was to provide new insights into the development of a dense core during the cutter indentation process from the perspective of the geometric evolution of the mesoscopic mechanical structure of rock. Limited by the length of this article, we were not able to thoroughly investigate the influence of parameters such as the cutter spacing, joint structure and confining stress.