# Research on Simulation Analysis Method of Microbial Cemented Sand Based on Discrete Element Method

Authors: Yang Tang; Guobin Xu; Jijian Lian; Yue Yan; Dengfeng Fu; Weichen Sun
Journal: Advances in Materials Science and Engineering
 (2019)

---

## Abstract

A simulation method for microbial cemented sand (MCS) based on the Two-Dimensional Particle Flow Code (PFC 2D) has been developed in this study. It consists of an identification of mesoscopic contact model for structural mesoscopic particles, development of morphological algorithm for irregular crystal particles, and classification and setting of mesoscopic parameters for compositional materials and particle size distribution simulation. Additionally, an acoustic emission algorithm based on moment tensor theory was developed for discrete element simulation analysis on fracture characteristic of cemented sand under the action of the force. The simulation method proposed in this article/paper may reflect the physical and mechanical characteristics of real microbial cemented sand based on physical experiments. Quantification of the fracture process of microbial cemented sand is possible by introducing a moment magnitude (MW) in discrete element simulation analysis. The material may have greater probability of fracture between the MW ranges of −6.8 and −6.6. The relationship between probability of fracture at different MWs and the MW follows the Gaussian curve. The research results are a new trial for fracture analysis on microbial cemented sand.


---

## Body

## 1. Introduction

Microbial-induced calcite precipitation (MICP) is an emerging technology that utilises environmental benign microbes to improve the microstructure and material composition of rocks and soils so as to enhance their engineering mechanical properties. This technology has been used in many different applications, including the reinforcement of rock and soil, prevention of sand liquefaction, and healing of concrete cracks. However, a detailed understanding of the mechanical properties and fracturing behaviour of MCS under load is still unclear due to the lack of effective testing methods, which severely hinders the further application and promotion of such technology.The discrete element method (DEM) is a numerical method that models the mesoscopic mechanical behaviour of granular materials. Many studies have made promising advances in the analysis of mechanical properties of nonbonded particle systems [1, 2]. However, for the numerical simulation of cemented materials under the action of microorganisms, there is relatively little literature to draw on. Evans and Khoubani [3, 4] developed novel methods to reduce the complexity of bonded particle systems, which provide inspiration to analyse such systems through numerical simulation. Wang [5] analysed the mechanical characteristics of artificial cemented sand from which the strong mechanical strength of cemented sand was attributed to the formation of stable force chains in the material structure. Feng [6, 7] developed a discrete element model to characterize cemented sands and analysed their mechanical behaviour considering the effect of particle bonding and particle shape. All these studies provide valuable guidelines for developing a numerical simulation model for microbial-induced cemented materials.Nevertheless, most of the developed DEM models so far are deficient from the in-depth details of the mesoscale bonding model in the cemented structure with regard to the following: (1) the coarse composition of the cemented sand body and during identification aspect of the mesoscopic contact model; only the mesoscopic contact between the sand particles and the sand particles is considered, i.e., it is not refined enough; (2) the simulation of the irregular particle shape inside the material being characterized by circular particles, thereby not truly reflecting on the irregular shape of the particle. Particles with irregular shapes exhibit more physical contacts and different contact models (mainly caused by the spinning of the particles) with other particles in comparison to those with spherical shapes [8]. Such contacts have a substantial impact on the mechanical behaviour of the materials [9, 10]; (3) the characterization of key bonding parameters and the distribution of the bonding strength at different positions inside the structure are unclear in the DEM models. Considering all the above discussed aspects, there is still considerable scope to further improve DEM modelling methods.Structural fractures in MCS materials obtained by the MICP technique are generally initiated by mesoscale crack nucleation and eventually expand to macroscopic fracture features during the loading process. Although it is possible to characterize the fracturing behaviour of an MCS material under the laboratory tests by obtaining acoustic emission events, such as uniaxial compression tests, triaxial compression tests, and tensile tests, etc., it is challenging to record those features through the whole fracturing process due to the limited storage capacity and recording speed of the experimental instruments [11]. The moment tensor theory, which was developed in the context of seismology, is a powerful tool that can be used to analyse the fracturing behaviour of MCS materials. So far, the moment tensor theory has been widely used in the geotechnical field in studying structural fractures [12–14]. However, there has been no study on the fracturing process of MCS materials based on the moment tensor theory.This study has developed a new method to construct a numerical model of MCS materials using the Two-Dimensional Particle Flow Code (PFC 2D), based on a cementation sample obtained under typical experimental conditions. The method includes an identification of mesoscopic contact model, a novel irregular aggregate-generation algorithm that characterizes the irregular shapes of crystalline particles, characterization of bonding parameters and setting of distribution of the bonding strength, etc. Finally, an acoustic emission algorithm was embedded in the constitutive model to simulate and analyse the fracturing behaviour of MCS materials.


## 2. Construction of Mesoscopic Model for MCS Materials Using PFC Software

Currently, there are two different versions of PFC 5.0 software, namely; PFC 2D and PFC 3D. Although PFC 3D is more helpful in modelling real 3D structures, it has the limitation of considerably more computational time and other controlling aspects associated with loading boundary conditions in 3D simulation [15]. On the contrary, 2D simulation can provide insights into the mechanisms of key phenomena with much less computational time while capturing the complex material mechanical properties with consideration of the material shape effect [16–19]. Therefore, the PFC 2D 5.0 version has been adopted to perform the analysis on the MCS material.In this section, various mesoscopic contact models and determination of the composition of the MCS material will first be discussed, followed by introduction of the simulation method for capturing the particle-shape effect. Furthermore, the model construction process will be demonstrated based on an MCS sample obtained under a typical experimental condition (OD600 = 1.50, [Ca2+] = 0.75 mol/l). Finally, the authors will propose a classification and setting method of mesoscopic parameters selection in the numerical model.


### 2.1. Methods for Determining the Mesoscopic Contact Model and the Composition of the MCS Material

#### 2.1.1. Determination of the Mesoscopic Contact Model

The microstructure of the microbial cemented sand structure observed by scanning electron microscopy is shown in Figure1 where three kinds of bonds exist inside the microbial cemented sand materials, namely, sand particle-sand particle (B1), sand particle-calcite (C), and calcite-calcite (C). These bonds share similar characteristics to the linear parallel bonding model (LPBM) described in PFC [20]. This is because the linear parallel bonding model includes two components: the linear component and the bonding component. The linear component describes the elastic relationship between the contact force and the relative displacement of the adjacent entities, which can be used to express the rigidity characteristics of MCS material, whereas the bonding component defines the physical connection at the interface of adjacent entities and can be used to simulate the cementation characteristics of the MCS material [21]. Moreover, the LPBM has been successfully used to explore the mechanical behaviour of MCS material [6]. As such, the linear parallel bonding model was adopted to simulate the bonding behaviour between sand particle-sand particle, sand particle-calcite, and calcite-calcite contacts.Figure 1
The microstructure of the microbial cemented sand structure observed by scanning electron microscopy.


#### 2.1.2. Methods for Determining the Composition Content of MCS Material

The main components of the MCS material include quartz and calcite for a given set of sample material obtained under typical testing conditions [22]. An accurate percentage composition of sand particles (quartz) and calcite particles in the MCS material was obtained through the following calculations:(1)
The pore filling percentP was calculated after the completion of microbial mineralisation: the pore filling percent was calculated as(1)P=VCaCO3V,where VCaCO3 is the total volume of calcium carbonate accumulated in the material (cm3) and V is the pore volume in the uncemented sand (cm3).The total volume of calcium carbonate,VCaCO3, was calculated as(2)VCaCO3=∑i=1nmiρCaCO3,where mi is the mass of calcium carbonate produced in the sand column during each cementation process, n is the total number of cementations, and ρCaCO3 is the density of calcium carbonate.(2)
Determination of the porosityn2d in the numerical simulation: the porosity obtained from the physical experiments under laboratory conditions was converted to a numerical value compatible with 2D simulation. The conversion was made with regard to the porosity conversion method proposed by Wang [23] since no study has previously focused on the porosity of MCS materials. The conversion equation is given by(3)n2d=0.42×nlab2+nlab×0.25,where n2d is the 2D porosity and nlab is the porosity measured under laboratory conditions.(3)
Determination of the percentage content of calcium carbonate (calcite): based on the 2D porosity calculated from equation (3), the pore area in the numerical model was calculated to obtain the filling content of calcium carbonate by multiplying the pore area by the pore filling percent derived from equation (1). Therefore, the percentage content of calcium carbonate in the numerical model CaCO3% was calculated as(4)CaCO3%=n2d×W×H×PW×H=P×n2d,where W and H are the width and height of the model, respectively. In this study, the dimensions of the numerical model are the same as the actual dimensions of the physical model.(4)
Determination of the percentage content of sand particles in the numerical model: the percentage content of sand particles was given by(5)Sand%=1−CaCO3%.It should be noted that the porosity of the sand particles nSand is equal to the difference between the porosity before pore filling and the percentage content of the filling calcium carbonate, namely,(6)nSand=n2d−n2d×P=n2d×1−P.


### 2.2. Simulation Methods for Capturing the Shapes of Particles

According to our prior study [22], most sand particles selected under typical testing conditions have a spherical shape. Considering that previous researchers have assumed a spherical shape in characterization [5], similar spherical sand particle shapes were used in this study. For calcite crystals, studies have shown that they exhibit irregular geometries [22]; the shape of the particles directly affects the interaction between different particles and eventually the mechanical behaviour of the materials [9]. Therefore, particle clumps with irregular geometries were implemented in the simulation. Tang et al. [24] developed a method to generate particles with irregular geometries. However, postsimulation revealed that such a method can only produce irregular particles with high efficiency when the size of the model is relatively small. When the model size exceeds 100 mm × 50 mm, the method developed by Tang et al. becomes highly inefficient and yields coarse crystal grains. Thus, it is essential to optimise such a method for generating irregular particles. The general idea is to first pregenerate seeding particles with a uniform size in the model and subsequently regenerate irregular particles using irregular polygons generated randomly inside the model. To ensure both high efficiency and accuracy, the size of the seeding particles was set at half the minimal radius of the calcite particles in the physical sample. Figure 2 shows the evolution process for generating irregular particles based on spherical seeding particles. The detailed process is described below.Figure 2
Schematic of the process for generating irregular polygons. (a) Generating initial body. (b) Generating characteristic circle. (c) Determining vertices of the irregular polygons. (d) Generating irregular polygons.
(a)
(b)
(c)
(d)


#### 2.2.1. Generating Initial Particle Body

The initial circular body was generated with a characteristic length ofR and centre coordinate of A (xc, yc), as shown in Figure 2(a). The radius R is given by(7)R=R¯+Rstd2K−1,where R¯ is the average radius of the crystal (mean value of the maximum and minimum radius of crystal particles), Rstd is the variation factor of the crystal radius (half of the difference between the maximum and minimum radius of the crystal particles), and K is a random number generated between 0 and 1, and the distribution is uniformly distribution.


#### 2.2.2. Generating Characteristic Circle

An inner circleC2 and an outer circle C3 were generated by modifying the radius of the initial circle C1 located at A (xc, yc) with an increment of ΔR (shown in Figure 2(b)). The characteristic radii of the inner circle C2 and outer circle C3 are given by(8)Rin=R−ΔR,Rout=R+ΔR.The increment parameterΔR was calculated by(9)ΔR=γ×R×1+K2,where γ is a shape factor which represents different particle geometries with difference values and K is a random number generated between 0 and 1, and the distribution is uniform distribution.


#### 2.2.3. Determining Vertices of the Irregular Polygon

The internal angle of the inner circle C2 was divided into multiple identical central angles according to equation (10). With the introduction of the random number, the minimum and maximum numbers of central angles are N−Nstd and N+Nstd, respectively. For the minimum case, the number of central angles, k, was selected from the 1,N−Nstd interval; while for the maximum case, the number of central angles, k, was selected from the 1,N+Nstd interval.Subsequently, a local coordination system was constructed for defining the vertices of the irregular polygon. The characteristic radiusrk at these vertices was calculated by equation (11). It is worth noting that while the values of rk at different points may be different due to the introduction of the random number when calculating the characteristic radius, the overall characteristic radius must be greater than or equal to Rin (the radius of the inner circle) as shown in Figure 2(c). Tk and Tk+1 shown in Figure 2(c) represent the kth and (k+1)th vertices of the irregular polygon, while θk is the central angle formed between these two vertices with respect to the centre:(10)θk=2πN+2K−1×Nstd,k∈1,N−Nstd,ork∈1,N+Nstd,(11)rk=Rin+ΔRk×K,ΔRk=Rout−Rin,where N is the number of the sides of the irregular polygon, Nstd is the variation factor of the irregular polygon, and K is a random number generated between 0 and 1, and the distribution is uniform distribution.The local coordinates of the vertices of the irregular polygon are given by(12)xk=rk×cosθk×k,yk=rk×sinθk×k,where xk and yk are the x-coordinate and y-coordinate of the kth vertex of the irregular polygon, θk is the kth central angle of the polygon, and k is the number of the central angles.


#### 2.2.4. Generating Irregular Polygon

Sharing of the same coordinate system during the generation of irregular polygons will lead to a uniform pattern of irregular geometries. Therefore, the following approaches were used to enhance the randomness and disorderly direction of the polygons. First, a global coordinate system was constructed after generating all the vertices of the irregular geometries. Then, these polygon vertices were rotated randomly in the model using the rotating equation (13) given by(13)x′=xc+xk−xccosθ−yk−ycsinθ,y′=yc+xk−xcsinθ+yk−yccosθ,θ=2π×K,where xc and yc are the centre coordinates of the randomly generated polygon, xk and yk are the coordinates of the vertices in the local coordinate system before the rotation, x′ and y′ are the global coordinates of these vertices after rotation, θ is the rotating angle, and K is a random number generated between 0 and 1, and the distribution is uniform distribution.After all the vertices were rotated in the global coordinates, they were then connected in sequence to generate the irregular polygons shown in Figure2(d). Once these irregular polygons were formed, reconfiguring the seeding particles was started. If the centre of the seeding particles fell inside a polygon, it would then belong to that specific polygon and vice versa. Following this strategy, groups of irregular polygons were eventually generated in our simulation model.Another approach to increase the divergence of the irregular geometries in the model was to construct the initial body with an oval shape instead of a circular shape (the length ratio of the major and minor axes of a circle was fixed at one, while that of an oval can be any value). By introducing the length ratio of the major and minor axes as an additional parameter varying from 1 to 2, the initial bodies generated in our model will always contain two different shapes. The subsequent generation of irregular polygons based on oval initial geometries follows the same steps outlined above.


### 2.3. Methods for Constructing the Numerical Model for Typical Testing Conditions

Combining the methods discussed in Section2.1.2 with our prior testing result [22], for typical testing conditions, the pore filling percent of 18.2 was found for calcite crystals obtained from the sand column with two cycles’ cementations having the porosity of 0.405 which was measured in indoor experiments before mineralisation. Therefore, the 2D porosity and the percentage content of calcium carbonate were calculated as 0.17 and 3% according to equations (3) and (4), respectively. Furthermore, the percentage content and porosity of sand particles were calculated as 97% and 0.14 from equations (5) and (6), respectively. And so, our numerical model consists of 3% calcium carbonate (calcite) and 97% sand particles with a sand particle porosity of 0.14.During the modelling process, sand particles were simulated with spherical geometries in which the sizes of the particles were determined by the actual particle size distribution [22]. The calcite crystal particles were simulated using the methods described in section 2.2, where the maximum and minimum particle radii were calculated following the equivalent area principle [25] based on the XCT scan result of real testing samples. The maximum and minimum particle radii were found as 0.5 mm and 0.075 mm, respectively.The dimensions of the numerical model were 100 × 50 mm, equivalent to the actual dimensions of the physical model. The total number of sand particles and calcite particles in the numerical model were 165,244 and 33,955, respectively. The granular particle size distribution between the numerical model and actual samples is shown in Figure3, which demonstrates high consistency between the two. A closer examination of the calcite crystals obtained using the optimised irregular particle generation algorithm (γ = 0.25) suggests that such a method provides a better description of the crystal boundary features and therefore a more accurate characterization of the sample in comparison to the conventional algorithm [24] (Table 1). Furthermore, the morphology of the crystals generated using the new algorithm shows a strong similarity to the XCT scanned morphology of the particles obtained under typical testing conditions (Figure 4). Therefore, the optimised irregular particle generation method can accurately capture the morphological features of crystal particles.Figure 3
Comparison of granular particle size distribution of numerical model and actual sample.Table 1
Comparison of morphological features of particles generated with different algorithms.
The new algorithm
Tang et al. [24]Figure 4
Comparison of cross-sectional characteristic features of calcite particles. (a) Cross Section. (b) Longitudinal section.


### 2.4. Calibration Method and Initial Mesoscopic Contact Parameters Selection in the Numerical Model

There are two sets of mesoscopic contact parameters employed in the linear parallel bonding model in the PFC 2D 5.0, namely, (1) the particle elastic moduleE∗ (N/m2), stiffness ratio k∗, and friction coefficient μ associated with the linear contact component of the model; (2) the bonding elastic module E¯∗ (N/m2), bonding stiffness ratio k¯∗, normal bonding strength pb_ten (Pa), cohesive strength pb_coh (Pa), and bonding radius R¯ (m) associated with the parallel bonding component of the model.Since there is no direct quantitative relationship between the input mesoscopic contact parameters and the macroscopic properties of MCS material such as macroscopic elastic module, uniaxial compression strength, and peak strain, “trial and error” has been the most commonly used calibration method to obtain appropriate mesoscopic parameters [15, 24, 26, 27].As suggested by Potyondy [20], the same elastic module and stiffness ratio should be assigned to both the particle and bonding parameters (i.e., E∗=E¯∗, k∗=k¯∗) to reduce the number of unknowns. In this study, the same strategy was followed with the same values of elastic module and stiffness ratio parameters being assigned to the particle model and bonding model, respectively. When determining the elastic modules (i.e., E∗ and E¯∗) of the MCS material, reference was made to both the previous PFC study on cemented sand column [6] and the indoor elastic module testing results obtained by Porter [28] (sand particle > calcite > calcite-sand particle). Finally, the following values were assigned to the elastic module parameters for different material components: for sand particle, E∗=E¯∗ = 1.50 GPa; for calcite particle, E∗=E¯∗ = 1.25 GPa; and for calcite-sand particle, E∗=E¯∗ = 0.85 GPa. When determining the stiffness ratio, recent studies [29] have suggested that the macroscopic Poisson’s ratio can be derived from the stiffness ratio using equation (14). Therefore, the stiffness ratio can be back-calculated with known Poisson’s ratio. Poisson’s ratios for quartz and calcite were found to be 0.17 and 0.31 from a literature review [30, 31], respectively. Thus, the associated stiffness ratios were calculated as k∗=k¯∗ = 1.2 for sand particles and k∗=k¯∗ = 2.5 for calcite, respectively. These values were assigned to the stiffness ratio parameters in our model. Since the main component at the interface in the calcite-sand particle system is calcite, the stiffness ratio was approximated as the inherent property of calcite, which gives k∗=k¯∗ = 2.5:(14)v=0.14+0.20lnknks,where ν is Poisson’s ratio, kn is the normal stiffness, and ks is the tangential stiffness.When determining the parallel bonding radiusR¯, the particle flow theory suggests that the mesoscale bonding radius R¯ is related to the radial coefficient λ. The SEM characterization of the material (shown in Figure 1) indicates that there exist three different types of bonding inside the microbial cemented sand material, namely, sand particle-sand particle, sand particle-calcite, and calcite-calcite contact bonds. Among all three bonds, the sand particle-sand particle bond has the lowest strength but the highest discrete level. Therefore, for the sand particle-calcite and calcite-calcite bonds, λ is assigned with a fixed value of 0.5 because these bonds are relatively strong according to SEM analysis based on the previous PFC study on microbial cemented sand column by Feng [6]. For the sand particle-sand particle bond, a random number was introduced to vary the parallel bond radial coefficient λ between 0.25 and 0.5.Thepb_ten and pb_coh parameters of the MCS material were determined by setting up the initial values according to previous studies by Ali et al. [32] and Xue et al. [33]. Subsequently, these values were adjusted through “trial and error”. It should be noted that different MCS materials were found at different bonding locations and a porous structure was also observed at the bonding sites (as shown in Figure 1). This inevitably led to a nonuniform distribution of the mechanical strength. The mesoscale mechanical strength and cohesive strength were determined by introducing the Gaussian distribution function in the PFC model aimed at achieving a discrete strength distribution. Such a method has also been used in previous studies by Potyondy [20] and Tang [24]. For the friction coefficient μ, most scholars [5, 34] use a two-dimensional disc to simulate sand particles; the friction coefficient is generally chosen to be around 0.5. In this paper, considering the influence of residual cement between sand particles on the sliding between particles, it is considered that the friction coefficient between sand particles should be greater than 0.5. In addition, due to the irregular shape of the calcite particles, the bite interaction between the particles has been considered, so it is considered to be around 0.5. After a large number of trial and error tests, the friction coefficients between sand particles and calcite were determined to be 0.7 and 0.5, respectively.Subsequent to determination of all the initial mesoscopic contact parameters of the simulation of the microbial cemented sand structure, the mesoscopic contact parameters were calibrated by trial and error method until the true characteristics of microbial cemented sand structure were met. Thus, this process was achieved in accordance with the simulation results corresponding to the initial parameters and the physical test results based on typical working conditions.


## 2.1. Methods for Determining the Mesoscopic Contact Model and the Composition of the MCS Material

### 2.1.1. Determination of the Mesoscopic Contact Model

The microstructure of the microbial cemented sand structure observed by scanning electron microscopy is shown in Figure1 where three kinds of bonds exist inside the microbial cemented sand materials, namely, sand particle-sand particle (B1), sand particle-calcite (C), and calcite-calcite (C). These bonds share similar characteristics to the linear parallel bonding model (LPBM) described in PFC [20]. This is because the linear parallel bonding model includes two components: the linear component and the bonding component. The linear component describes the elastic relationship between the contact force and the relative displacement of the adjacent entities, which can be used to express the rigidity characteristics of MCS material, whereas the bonding component defines the physical connection at the interface of adjacent entities and can be used to simulate the cementation characteristics of the MCS material [21]. Moreover, the LPBM has been successfully used to explore the mechanical behaviour of MCS material [6]. As such, the linear parallel bonding model was adopted to simulate the bonding behaviour between sand particle-sand particle, sand particle-calcite, and calcite-calcite contacts.Figure 1
The microstructure of the microbial cemented sand structure observed by scanning electron microscopy.


### 2.1.2. Methods for Determining the Composition Content of MCS Material

The main components of the MCS material include quartz and calcite for a given set of sample material obtained under typical testing conditions [22]. An accurate percentage composition of sand particles (quartz) and calcite particles in the MCS material was obtained through the following calculations:(1)
The pore filling percentP was calculated after the completion of microbial mineralisation: the pore filling percent was calculated as(1)P=VCaCO3V,where VCaCO3 is the total volume of calcium carbonate accumulated in the material (cm3) and V is the pore volume in the uncemented sand (cm3).The total volume of calcium carbonate,VCaCO3, was calculated as(2)VCaCO3=∑i=1nmiρCaCO3,where mi is the mass of calcium carbonate produced in the sand column during each cementation process, n is the total number of cementations, and ρCaCO3 is the density of calcium carbonate.(2)
Determination of the porosityn2d in the numerical simulation: the porosity obtained from the physical experiments under laboratory conditions was converted to a numerical value compatible with 2D simulation. The conversion was made with regard to the porosity conversion method proposed by Wang [23] since no study has previously focused on the porosity of MCS materials. The conversion equation is given by(3)n2d=0.42×nlab2+nlab×0.25,where n2d is the 2D porosity and nlab is the porosity measured under laboratory conditions.(3)
Determination of the percentage content of calcium carbonate (calcite): based on the 2D porosity calculated from equation (3), the pore area in the numerical model was calculated to obtain the filling content of calcium carbonate by multiplying the pore area by the pore filling percent derived from equation (1). Therefore, the percentage content of calcium carbonate in the numerical model CaCO3% was calculated as(4)CaCO3%=n2d×W×H×PW×H=P×n2d,where W and H are the width and height of the model, respectively. In this study, the dimensions of the numerical model are the same as the actual dimensions of the physical model.(4)
Determination of the percentage content of sand particles in the numerical model: the percentage content of sand particles was given by(5)Sand%=1−CaCO3%.It should be noted that the porosity of the sand particles nSand is equal to the difference between the porosity before pore filling and the percentage content of the filling calcium carbonate, namely,(6)nSand=n2d−n2d×P=n2d×1−P.


## 2.1.1. Determination of the Mesoscopic Contact Model

The microstructure of the microbial cemented sand structure observed by scanning electron microscopy is shown in Figure1 where three kinds of bonds exist inside the microbial cemented sand materials, namely, sand particle-sand particle (B1), sand particle-calcite (C), and calcite-calcite (C). These bonds share similar characteristics to the linear parallel bonding model (LPBM) described in PFC [20]. This is because the linear parallel bonding model includes two components: the linear component and the bonding component. The linear component describes the elastic relationship between the contact force and the relative displacement of the adjacent entities, which can be used to express the rigidity characteristics of MCS material, whereas the bonding component defines the physical connection at the interface of adjacent entities and can be used to simulate the cementation characteristics of the MCS material [21]. Moreover, the LPBM has been successfully used to explore the mechanical behaviour of MCS material [6]. As such, the linear parallel bonding model was adopted to simulate the bonding behaviour between sand particle-sand particle, sand particle-calcite, and calcite-calcite contacts.Figure 1
The microstructure of the microbial cemented sand structure observed by scanning electron microscopy.


## 2.1.2. Methods for Determining the Composition Content of MCS Material

The main components of the MCS material include quartz and calcite for a given set of sample material obtained under typical testing conditions [22]. An accurate percentage composition of sand particles (quartz) and calcite particles in the MCS material was obtained through the following calculations:(1)
The pore filling percentP was calculated after the completion of microbial mineralisation: the pore filling percent was calculated as(1)P=VCaCO3V,where VCaCO3 is the total volume of calcium carbonate accumulated in the material (cm3) and V is the pore volume in the uncemented sand (cm3).The total volume of calcium carbonate,VCaCO3, was calculated as(2)VCaCO3=∑i=1nmiρCaCO3,where mi is the mass of calcium carbonate produced in the sand column during each cementation process, n is the total number of cementations, and ρCaCO3 is the density of calcium carbonate.(2)
Determination of the porosityn2d in the numerical simulation: the porosity obtained from the physical experiments under laboratory conditions was converted to a numerical value compatible with 2D simulation. The conversion was made with regard to the porosity conversion method proposed by Wang [23] since no study has previously focused on the porosity of MCS materials. The conversion equation is given by(3)n2d=0.42×nlab2+nlab×0.25,where n2d is the 2D porosity and nlab is the porosity measured under laboratory conditions.(3)
Determination of the percentage content of calcium carbonate (calcite): based on the 2D porosity calculated from equation (3), the pore area in the numerical model was calculated to obtain the filling content of calcium carbonate by multiplying the pore area by the pore filling percent derived from equation (1). Therefore, the percentage content of calcium carbonate in the numerical model CaCO3% was calculated as(4)CaCO3%=n2d×W×H×PW×H=P×n2d,where W and H are the width and height of the model, respectively. In this study, the dimensions of the numerical model are the same as the actual dimensions of the physical model.(4)
Determination of the percentage content of sand particles in the numerical model: the percentage content of sand particles was given by(5)Sand%=1−CaCO3%.It should be noted that the porosity of the sand particles nSand is equal to the difference between the porosity before pore filling and the percentage content of the filling calcium carbonate, namely,(6)nSand=n2d−n2d×P=n2d×1−P.


## 2.2. Simulation Methods for Capturing the Shapes of Particles

According to our prior study [22], most sand particles selected under typical testing conditions have a spherical shape. Considering that previous researchers have assumed a spherical shape in characterization [5], similar spherical sand particle shapes were used in this study. For calcite crystals, studies have shown that they exhibit irregular geometries [22]; the shape of the particles directly affects the interaction between different particles and eventually the mechanical behaviour of the materials [9]. Therefore, particle clumps with irregular geometries were implemented in the simulation. Tang et al. [24] developed a method to generate particles with irregular geometries. However, postsimulation revealed that such a method can only produce irregular particles with high efficiency when the size of the model is relatively small. When the model size exceeds 100 mm × 50 mm, the method developed by Tang et al. becomes highly inefficient and yields coarse crystal grains. Thus, it is essential to optimise such a method for generating irregular particles. The general idea is to first pregenerate seeding particles with a uniform size in the model and subsequently regenerate irregular particles using irregular polygons generated randomly inside the model. To ensure both high efficiency and accuracy, the size of the seeding particles was set at half the minimal radius of the calcite particles in the physical sample. Figure 2 shows the evolution process for generating irregular particles based on spherical seeding particles. The detailed process is described below.Figure 2
Schematic of the process for generating irregular polygons. (a) Generating initial body. (b) Generating characteristic circle. (c) Determining vertices of the irregular polygons. (d) Generating irregular polygons.
(a)
(b)
(c)
(d)


### 2.2.1. Generating Initial Particle Body

The initial circular body was generated with a characteristic length ofR and centre coordinate of A (xc, yc), as shown in Figure 2(a). The radius R is given by(7)R=R¯+Rstd2K−1,where R¯ is the average radius of the crystal (mean value of the maximum and minimum radius of crystal particles), Rstd is the variation factor of the crystal radius (half of the difference between the maximum and minimum radius of the crystal particles), and K is a random number generated between 0 and 1, and the distribution is uniformly distribution.


### 2.2.2. Generating Characteristic Circle

An inner circleC2 and an outer circle C3 were generated by modifying the radius of the initial circle C1 located at A (xc, yc) with an increment of ΔR (shown in Figure 2(b)). The characteristic radii of the inner circle C2 and outer circle C3 are given by(8)Rin=R−ΔR,Rout=R+ΔR.The increment parameterΔR was calculated by(9)ΔR=γ×R×1+K2,where γ is a shape factor which represents different particle geometries with difference values and K is a random number generated between 0 and 1, and the distribution is uniform distribution.


### 2.2.3. Determining Vertices of the Irregular Polygon

The internal angle of the inner circle C2 was divided into multiple identical central angles according to equation (10). With the introduction of the random number, the minimum and maximum numbers of central angles are N−Nstd and N+Nstd, respectively. For the minimum case, the number of central angles, k, was selected from the 1,N−Nstd interval; while for the maximum case, the number of central angles, k, was selected from the 1,N+Nstd interval.Subsequently, a local coordination system was constructed for defining the vertices of the irregular polygon. The characteristic radiusrk at these vertices was calculated by equation (11). It is worth noting that while the values of rk at different points may be different due to the introduction of the random number when calculating the characteristic radius, the overall characteristic radius must be greater than or equal to Rin (the radius of the inner circle) as shown in Figure 2(c). Tk and Tk+1 shown in Figure 2(c) represent the kth and (k+1)th vertices of the irregular polygon, while θk is the central angle formed between these two vertices with respect to the centre:(10)θk=2πN+2K−1×Nstd,k∈1,N−Nstd,ork∈1,N+Nstd,(11)rk=Rin+ΔRk×K,ΔRk=Rout−Rin,where N is the number of the sides of the irregular polygon, Nstd is the variation factor of the irregular polygon, and K is a random number generated between 0 and 1, and the distribution is uniform distribution.The local coordinates of the vertices of the irregular polygon are given by(12)xk=rk×cosθk×k,yk=rk×sinθk×k,where xk and yk are the x-coordinate and y-coordinate of the kth vertex of the irregular polygon, θk is the kth central angle of the polygon, and k is the number of the central angles.


### 2.2.4. Generating Irregular Polygon

Sharing of the same coordinate system during the generation of irregular polygons will lead to a uniform pattern of irregular geometries. Therefore, the following approaches were used to enhance the randomness and disorderly direction of the polygons. First, a global coordinate system was constructed after generating all the vertices of the irregular geometries. Then, these polygon vertices were rotated randomly in the model using the rotating equation (13) given by(13)x′=xc+xk−xccosθ−yk−ycsinθ,y′=yc+xk−xcsinθ+yk−yccosθ,θ=2π×K,where xc and yc are the centre coordinates of the randomly generated polygon, xk and yk are the coordinates of the vertices in the local coordinate system before the rotation, x′ and y′ are the global coordinates of these vertices after rotation, θ is the rotating angle, and K is a random number generated between 0 and 1, and the distribution is uniform distribution.After all the vertices were rotated in the global coordinates, they were then connected in sequence to generate the irregular polygons shown in Figure2(d). Once these irregular polygons were formed, reconfiguring the seeding particles was started. If the centre of the seeding particles fell inside a polygon, it would then belong to that specific polygon and vice versa. Following this strategy, groups of irregular polygons were eventually generated in our simulation model.Another approach to increase the divergence of the irregular geometries in the model was to construct the initial body with an oval shape instead of a circular shape (the length ratio of the major and minor axes of a circle was fixed at one, while that of an oval can be any value). By introducing the length ratio of the major and minor axes as an additional parameter varying from 1 to 2, the initial bodies generated in our model will always contain two different shapes. The subsequent generation of irregular polygons based on oval initial geometries follows the same steps outlined above.


## 2.2.1. Generating Initial Particle Body

The initial circular body was generated with a characteristic length ofR and centre coordinate of A (xc, yc), as shown in Figure 2(a). The radius R is given by(7)R=R¯+Rstd2K−1,where R¯ is the average radius of the crystal (mean value of the maximum and minimum radius of crystal particles), Rstd is the variation factor of the crystal radius (half of the difference between the maximum and minimum radius of the crystal particles), and K is a random number generated between 0 and 1, and the distribution is uniformly distribution.


## 2.2.2. Generating Characteristic Circle

An inner circleC2 and an outer circle C3 were generated by modifying the radius of the initial circle C1 located at A (xc, yc) with an increment of ΔR (shown in Figure 2(b)). The characteristic radii of the inner circle C2 and outer circle C3 are given by(8)Rin=R−ΔR,Rout=R+ΔR.The increment parameterΔR was calculated by(9)ΔR=γ×R×1+K2,where γ is a shape factor which represents different particle geometries with difference values and K is a random number generated between 0 and 1, and the distribution is uniform distribution.


## 2.2.3. Determining Vertices of the Irregular Polygon

The internal angle of the inner circle C2 was divided into multiple identical central angles according to equation (10). With the introduction of the random number, the minimum and maximum numbers of central angles are N−Nstd and N+Nstd, respectively. For the minimum case, the number of central angles, k, was selected from the 1,N−Nstd interval; while for the maximum case, the number of central angles, k, was selected from the 1,N+Nstd interval.Subsequently, a local coordination system was constructed for defining the vertices of the irregular polygon. The characteristic radiusrk at these vertices was calculated by equation (11). It is worth noting that while the values of rk at different points may be different due to the introduction of the random number when calculating the characteristic radius, the overall characteristic radius must be greater than or equal to Rin (the radius of the inner circle) as shown in Figure 2(c). Tk and Tk+1 shown in Figure 2(c) represent the kth and (k+1)th vertices of the irregular polygon, while θk is the central angle formed between these two vertices with respect to the centre:(10)θk=2πN+2K−1×Nstd,k∈1,N−Nstd,ork∈1,N+Nstd,(11)rk=Rin+ΔRk×K,ΔRk=Rout−Rin,where N is the number of the sides of the irregular polygon, Nstd is the variation factor of the irregular polygon, and K is a random number generated between 0 and 1, and the distribution is uniform distribution.The local coordinates of the vertices of the irregular polygon are given by(12)xk=rk×cosθk×k,yk=rk×sinθk×k,where xk and yk are the x-coordinate and y-coordinate of the kth vertex of the irregular polygon, θk is the kth central angle of the polygon, and k is the number of the central angles.


## 2.2.4. Generating Irregular Polygon

Sharing of the same coordinate system during the generation of irregular polygons will lead to a uniform pattern of irregular geometries. Therefore, the following approaches were used to enhance the randomness and disorderly direction of the polygons. First, a global coordinate system was constructed after generating all the vertices of the irregular geometries. Then, these polygon vertices were rotated randomly in the model using the rotating equation (13) given by(13)x′=xc+xk−xccosθ−yk−ycsinθ,y′=yc+xk−xcsinθ+yk−yccosθ,θ=2π×K,where xc and yc are the centre coordinates of the randomly generated polygon, xk and yk are the coordinates of the vertices in the local coordinate system before the rotation, x′ and y′ are the global coordinates of these vertices after rotation, θ is the rotating angle, and K is a random number generated between 0 and 1, and the distribution is uniform distribution.After all the vertices were rotated in the global coordinates, they were then connected in sequence to generate the irregular polygons shown in Figure2(d). Once these irregular polygons were formed, reconfiguring the seeding particles was started. If the centre of the seeding particles fell inside a polygon, it would then belong to that specific polygon and vice versa. Following this strategy, groups of irregular polygons were eventually generated in our simulation model.Another approach to increase the divergence of the irregular geometries in the model was to construct the initial body with an oval shape instead of a circular shape (the length ratio of the major and minor axes of a circle was fixed at one, while that of an oval can be any value). By introducing the length ratio of the major and minor axes as an additional parameter varying from 1 to 2, the initial bodies generated in our model will always contain two different shapes. The subsequent generation of irregular polygons based on oval initial geometries follows the same steps outlined above.


## 2.3. Methods for Constructing the Numerical Model for Typical Testing Conditions

Combining the methods discussed in Section2.1.2 with our prior testing result [22], for typical testing conditions, the pore filling percent of 18.2 was found for calcite crystals obtained from the sand column with two cycles’ cementations having the porosity of 0.405 which was measured in indoor experiments before mineralisation. Therefore, the 2D porosity and the percentage content of calcium carbonate were calculated as 0.17 and 3% according to equations (3) and (4), respectively. Furthermore, the percentage content and porosity of sand particles were calculated as 97% and 0.14 from equations (5) and (6), respectively. And so, our numerical model consists of 3% calcium carbonate (calcite) and 97% sand particles with a sand particle porosity of 0.14.During the modelling process, sand particles were simulated with spherical geometries in which the sizes of the particles were determined by the actual particle size distribution [22]. The calcite crystal particles were simulated using the methods described in section 2.2, where the maximum and minimum particle radii were calculated following the equivalent area principle [25] based on the XCT scan result of real testing samples. The maximum and minimum particle radii were found as 0.5 mm and 0.075 mm, respectively.The dimensions of the numerical model were 100 × 50 mm, equivalent to the actual dimensions of the physical model. The total number of sand particles and calcite particles in the numerical model were 165,244 and 33,955, respectively. The granular particle size distribution between the numerical model and actual samples is shown in Figure3, which demonstrates high consistency between the two. A closer examination of the calcite crystals obtained using the optimised irregular particle generation algorithm (γ = 0.25) suggests that such a method provides a better description of the crystal boundary features and therefore a more accurate characterization of the sample in comparison to the conventional algorithm [24] (Table 1). Furthermore, the morphology of the crystals generated using the new algorithm shows a strong similarity to the XCT scanned morphology of the particles obtained under typical testing conditions (Figure 4). Therefore, the optimised irregular particle generation method can accurately capture the morphological features of crystal particles.Figure 3
Comparison of granular particle size distribution of numerical model and actual sample.Table 1
Comparison of morphological features of particles generated with different algorithms.
The new algorithm
Tang et al. [24]Figure 4
Comparison of cross-sectional characteristic features of calcite particles. (a) Cross Section. (b) Longitudinal section.


## 2.4. Calibration Method and Initial Mesoscopic Contact Parameters Selection in the Numerical Model

There are two sets of mesoscopic contact parameters employed in the linear parallel bonding model in the PFC 2D 5.0, namely, (1) the particle elastic moduleE∗ (N/m2), stiffness ratio k∗, and friction coefficient μ associated with the linear contact component of the model; (2) the bonding elastic module E¯∗ (N/m2), bonding stiffness ratio k¯∗, normal bonding strength pb_ten (Pa), cohesive strength pb_coh (Pa), and bonding radius R¯ (m) associated with the parallel bonding component of the model.Since there is no direct quantitative relationship between the input mesoscopic contact parameters and the macroscopic properties of MCS material such as macroscopic elastic module, uniaxial compression strength, and peak strain, “trial and error” has been the most commonly used calibration method to obtain appropriate mesoscopic parameters [15, 24, 26, 27].As suggested by Potyondy [20], the same elastic module and stiffness ratio should be assigned to both the particle and bonding parameters (i.e., E∗=E¯∗, k∗=k¯∗) to reduce the number of unknowns. In this study, the same strategy was followed with the same values of elastic module and stiffness ratio parameters being assigned to the particle model and bonding model, respectively. When determining the elastic modules (i.e., E∗ and E¯∗) of the MCS material, reference was made to both the previous PFC study on cemented sand column [6] and the indoor elastic module testing results obtained by Porter [28] (sand particle > calcite > calcite-sand particle). Finally, the following values were assigned to the elastic module parameters for different material components: for sand particle, E∗=E¯∗ = 1.50 GPa; for calcite particle, E∗=E¯∗ = 1.25 GPa; and for calcite-sand particle, E∗=E¯∗ = 0.85 GPa. When determining the stiffness ratio, recent studies [29] have suggested that the macroscopic Poisson’s ratio can be derived from the stiffness ratio using equation (14). Therefore, the stiffness ratio can be back-calculated with known Poisson’s ratio. Poisson’s ratios for quartz and calcite were found to be 0.17 and 0.31 from a literature review [30, 31], respectively. Thus, the associated stiffness ratios were calculated as k∗=k¯∗ = 1.2 for sand particles and k∗=k¯∗ = 2.5 for calcite, respectively. These values were assigned to the stiffness ratio parameters in our model. Since the main component at the interface in the calcite-sand particle system is calcite, the stiffness ratio was approximated as the inherent property of calcite, which gives k∗=k¯∗ = 2.5:(14)v=0.14+0.20lnknks,where ν is Poisson’s ratio, kn is the normal stiffness, and ks is the tangential stiffness.When determining the parallel bonding radiusR¯, the particle flow theory suggests that the mesoscale bonding radius R¯ is related to the radial coefficient λ. The SEM characterization of the material (shown in Figure 1) indicates that there exist three different types of bonding inside the microbial cemented sand material, namely, sand particle-sand particle, sand particle-calcite, and calcite-calcite contact bonds. Among all three bonds, the sand particle-sand particle bond has the lowest strength but the highest discrete level. Therefore, for the sand particle-calcite and calcite-calcite bonds, λ is assigned with a fixed value of 0.5 because these bonds are relatively strong according to SEM analysis based on the previous PFC study on microbial cemented sand column by Feng [6]. For the sand particle-sand particle bond, a random number was introduced to vary the parallel bond radial coefficient λ between 0.25 and 0.5.Thepb_ten and pb_coh parameters of the MCS material were determined by setting up the initial values according to previous studies by Ali et al. [32] and Xue et al. [33]. Subsequently, these values were adjusted through “trial and error”. It should be noted that different MCS materials were found at different bonding locations and a porous structure was also observed at the bonding sites (as shown in Figure 1). This inevitably led to a nonuniform distribution of the mechanical strength. The mesoscale mechanical strength and cohesive strength were determined by introducing the Gaussian distribution function in the PFC model aimed at achieving a discrete strength distribution. Such a method has also been used in previous studies by Potyondy [20] and Tang [24]. For the friction coefficient μ, most scholars [5, 34] use a two-dimensional disc to simulate sand particles; the friction coefficient is generally chosen to be around 0.5. In this paper, considering the influence of residual cement between sand particles on the sliding between particles, it is considered that the friction coefficient between sand particles should be greater than 0.5. In addition, due to the irregular shape of the calcite particles, the bite interaction between the particles has been considered, so it is considered to be around 0.5. After a large number of trial and error tests, the friction coefficients between sand particles and calcite were determined to be 0.7 and 0.5, respectively.Subsequent to determination of all the initial mesoscopic contact parameters of the simulation of the microbial cemented sand structure, the mesoscopic contact parameters were calibrated by trial and error method until the true characteristics of microbial cemented sand structure were met. Thus, this process was achieved in accordance with the simulation results corresponding to the initial parameters and the physical test results based on typical working conditions.


## 3. Acoustic Emission Simulation Algorithm Based on the Moment Tensor Theory

The acoustic emission simulation method based on the moment tensor theory is an important tool that has been widely used in the field of geotechnical engineering to analyse many structural fracturing problems [12–14, 35]. In this study, the acoustic emission simulation method was employed based on the moment tensor theory to analyse the fracturing characteristics of the MCS material.


### 3.1. Determining Acoustic Emission Simulation Parameters

There are two primary parameters for setting up the acoustic emission simulation in PFC: the quality indexQs and the shear wave speed Vs. These two parameters can be calculated by the following equations:(15)Qs=π2α,(16)Vs=E21+νρ,where α is the damping coefficient, E is the elastic module (N/m2), ν is Poisson’s ratio, and ρ is the density (kg/m3). The quality index Qs is set as 100 for this study by reference [36].


### 3.2. Solving Moment Tensor and Moment Magnitude

For a bonding model, the particles in contact with each other at a microcrack before fracturing are defined as source particles. Once the bond of source particles breaks, the source particles will shift in location and induce further variation in the contact force. Therefore, based on the change of contact force and contact position, the individual component of the moment tensor can be calculated by multiplying the magnitude of variation in each single source particle contact force by its associated moment arm and summing up the products (shown in equation (17)). The contact force and the induced motion of the particles can be obtained directly from the simulation model:(17)Mij=∑sΔFiRj,where ΔFi is the ith component of the contact force variation and Rj is the jth component of the distance between the contact point and the microcrack centre.The moment tensor was calculated instantaneously during a period since the acoustic emission lasts for a period, i.e., the moment tensor was calculated with respect to every single time step during the whole acoustic emission period, hence taking up huge storage resources. Therefore, the computation efficiency was improved by characterizing the acoustic emission using the scalar torque. The scalar torque was calculated based on the moment tensor with the maximum scale torque magnitude within each acoustic emission period, as shown in the following equation [37]:(18)M0=∑j=13mj221/2,where mj is the jth eigenvalue of the moment tensor matrix.The eigenvaluemj (with j equal to 1, 2, and 3) was calculated by the following equation [38]:(19)m1=−2Qcosθ3−a3,m2=−2Qcosθ+2π3−a3,m3=−2Qcosθ−2π3−a3,Q=a2−3b9,R=2a3−9ab+27c54,θ=arctanQ3−R2R.Parameters a,b, and c in equation (19) were calculated by the following equation:(20)a=−M11−M22−M33,b=−M23×M32−M21×M12−M31×M13+M11×M22+M33×M11×M33×M22,c=−M11×M22×M33−M21×M13×M32−M31×M12×M23+M11×M23×M32+M21×M12×M33+M31×M13×M22,where Mij is the individual component of the moment tensor (with i and j equalling 1, 2, and 3). It should be noted that since the 2D numerical model is used in this paper, therefore M13 = M31 = M23 = M32 = M33 = 0.The moment magnitudeMw of the acoustic emission event was calculated by equation (21) [39] in combination with equation (17) and the scale torque calculated by equation (18):(21)Mw=23logM0−6.The moment magnitude indicates the breaking intensity of the fracturing event. It characterizes the displacement magnitude of the fracture surface and the energy released during fracturing [40].


## 3.1. Determining Acoustic Emission Simulation Parameters

There are two primary parameters for setting up the acoustic emission simulation in PFC: the quality indexQs and the shear wave speed Vs. These two parameters can be calculated by the following equations:(15)Qs=π2α,(16)Vs=E21+νρ,where α is the damping coefficient, E is the elastic module (N/m2), ν is Poisson’s ratio, and ρ is the density (kg/m3). The quality index Qs is set as 100 for this study by reference [36].


## 3.2. Solving Moment Tensor and Moment Magnitude

For a bonding model, the particles in contact with each other at a microcrack before fracturing are defined as source particles. Once the bond of source particles breaks, the source particles will shift in location and induce further variation in the contact force. Therefore, based on the change of contact force and contact position, the individual component of the moment tensor can be calculated by multiplying the magnitude of variation in each single source particle contact force by its associated moment arm and summing up the products (shown in equation (17)). The contact force and the induced motion of the particles can be obtained directly from the simulation model:(17)Mij=∑sΔFiRj,where ΔFi is the ith component of the contact force variation and Rj is the jth component of the distance between the contact point and the microcrack centre.The moment tensor was calculated instantaneously during a period since the acoustic emission lasts for a period, i.e., the moment tensor was calculated with respect to every single time step during the whole acoustic emission period, hence taking up huge storage resources. Therefore, the computation efficiency was improved by characterizing the acoustic emission using the scalar torque. The scalar torque was calculated based on the moment tensor with the maximum scale torque magnitude within each acoustic emission period, as shown in the following equation [37]:(18)M0=∑j=13mj221/2,where mj is the jth eigenvalue of the moment tensor matrix.The eigenvaluemj (with j equal to 1, 2, and 3) was calculated by the following equation [38]:(19)m1=−2Qcosθ3−a3,m2=−2Qcosθ+2π3−a3,m3=−2Qcosθ−2π3−a3,Q=a2−3b9,R=2a3−9ab+27c54,θ=arctanQ3−R2R.Parameters a,b, and c in equation (19) were calculated by the following equation:(20)a=−M11−M22−M33,b=−M23×M32−M21×M12−M31×M13+M11×M22+M33×M11×M33×M22,c=−M11×M22×M33−M21×M13×M32−M31×M12×M23+M11×M23×M32+M21×M12×M33+M31×M13×M22,where Mij is the individual component of the moment tensor (with i and j equalling 1, 2, and 3). It should be noted that since the 2D numerical model is used in this paper, therefore M13 = M31 = M23 = M32 = M33 = 0.The moment magnitudeMw of the acoustic emission event was calculated by equation (21) [39] in combination with equation (17) and the scale torque calculated by equation (18):(21)Mw=23logM0−6.The moment magnitude indicates the breaking intensity of the fracturing event. It characterizes the displacement magnitude of the fracture surface and the energy released during fracturing [40].


## 4. The Analysis of the Fracturing Characteristics Based on Typical Working Conditions

The fracturing characteristics of the MCS material were analysed based on typical test working conditions after establishing the modelling method developed in the previous sections. The uniaxial compression process was simulated with integration of a moment tensor algorithm that calculates the momentum tensor every time a bond breaks for each single numerical sample. This allows calculation of the moment magnitude and total number of breaking bonds in one acoustic emission event.


### 4.1. Calibration Results of MCS Material

According to the calibration method for mesoscopic contact parameters, uniaxial compression tests were performed using the initial mesoscopic parameter values given in section2.4. The numerical tests showed same loading rate in comparison to indoor experimental tests. The final calibrated parameters are summarized in Table 2. The mechanical parameters and testing results obtained from the numerical tests using calibrated parameters were compared with those obtained from experiments conducted under typical testing conditions (Table 3). Table 3 presents the comparison between the numerical and indoor experimental tests from which a close match of the mechanical parameters was found. Figure 5 shows the stress-strain curves and structural fracturing modes from both physical experiments and numerical simulations. The monitoring points A, B, C, D, and E shown on the numerical curve in Figure 5 correspond to strains of 1.0‰, 2.0‰, 3.0‰, 3.4‰, and 3.6‰, respectively. Among all these values, the peak strain is at 3.4‰. It can also be observed from Figure 5 that combined tension-shear fracturing is the primary mode of bond breaking in the numerical model. The numerical tests and physical tests have the following characteristics: (1) they show an identical developing trend although there is a small local deviation between the stress-stain curves. (2) The compressive strength and strain amplitude are basically the same. (3) The numerical tests revealed a similar trend between the macroscale cracking and the overall fracturing feature with minor local deviations.Table 2
Mesoscopic contact parameters.
Mesoscopic contact
E
∗ (N/m2)
E
¯
∗ (N/m2)
k
∗
k
¯
∗
µ
p
b
_
t
e
n (MPa)
p
b
_
c
o
h (MPa)
λ
Sand-sand
1.50 × 109
1.50 × 109
1.2
1.2
0.7
13.50 ± 1.0
15.50 ± 1.5
0.25–0.50
Sand-calcite
0.85 × 109
0.85 × 109
2.5
2.5
0.5
17.50 ± 1.0
24.50 ± 1.5
0.50
Calcite-calcite
1.25 × 109
1.25 × 109
2.5
2.5
0.5
14.50 ± 1.0
17.50 ± 1.5
0.50Table 3
Comparison of characteristic parameters.
Physical experiment
Numerical experiment
σ
p/MPa
4.61
4.80
E/GPa
1.32
1.42
ε
p/‰
3.5
3.4Figure 5
Comparison of stress-strain curves and structural fracturing modes from both physical experiments and numerical simulations.For the small local deviation between the stress-strain curves, the difference is derived from the mesostructural difference between the physical sample and the model sample. Although it is as close as possible to the real mesostructure, it is difficult to consider all the mesostructures of materials [20]. It is almost impossible to reconstruct a model in a numerical simulation sharing the same micro/mesoscale structure feature and physical properties as those of a physical testing sample [41]. For the difference in failure mode, it is more derived from the randomness of the model itself, which reflects the nonuniformity of the material. Due to the nonuniformity of the material itself, even the same material will not respond exactly the same under the same load. In addition, the macroscopic cracks in the numerical test are formed by microscopic cracks. Due to the discreteness of the particle contact in the numerical model, the generated microcracks are also discrete [42].Based on the above analysis results, the mesoscopic parameters used in general are reasonable and can be used to characterize the characteristic properties of MCS material under selected conditions.


### 4.2. The Fracturing Process and Crack Growth of MCS Material

Figure6 illustrates the fracturing process and crack growth of MCS material, where the coloured legend represents the contact force. Three stages have been illustrated in the fracturing process (Figure 6), namely, the stress-bearing stage, the stable crack growth stage, and the unstable crack growth stage.(1)
During the stress-bearing stage (A1), the MCS material was deformed internally and the internal contact force in the model was relatively small at around 0.5 kN. This level of contact force was insufficient to induce fracture of the structure. Consequently, no crack was observed in the material structure, and the strain increased up to 1.0‰.(2)
During the stable crack growth stage (A2) with further increase in the load, a stronger contact force is induced inside the material (e.g., a contact force yielding a strain of 2.0‰). Thus, once the contact force exceeded the bonding strength between the particles, a microcrack was formed inside the material. The number of cracks increased linearly during this stage.(3)
During the unstable crack growth stage, the contact force was distributed approximately between 0.5 and 2.5 kN. Owing to the significant change in the contact force, microcracks continued to force, expand, merge, and eventually turn into a macroscale fracturing zone inside the modelled material (e.g., the fracturing zone F1 with a strain of 3.6‰). Further, a sharp increase in the number and coverage of the cracks was observed, which expanded/grew/propagated nonlinearly with time with apparent bifurcations being observed at the fracturing zone. The bifurcations found in the fracturing process of the MCS material reflects the nonuniformity of the MCS material where the cracks started to expand in favourable paths with minimum resistive force under the applied stress. This implies that implementing a bonding strength distribution yields a more accurate result in numerical modelling.Figure 6
Distribution of contact force and crack growth during the fracturing process of MCS material.It is also worth mentioning that during the unstable crack growth stage when macroscale cracks start to form, a more concentrated stress was observed at some local cracking tip regions (e.g., the K region). However, the cracks did not expand further towards these areas because the dynamic stress intensification factor at the cracking tip was smaller than the minimum fracture toughness required to induce a crack. Previous studies on fracture mechanics [43] have shown that the expansion of a crack creates a stress field with the accumulation of elastic energy. Therefore, when the dynamic stress intensification factor at the cracking tip is smaller than the critical fracture toughness, the crack ceases to expand. Whereas when the dynamic stress intensification factor at the cracking tip exceeds the critical fracture toughness, the crack continues to expand. The simulation results obtained in this study are in accordance with the first case, i.e., cracks did not expand further towards the stress intensified regions.By performing regression analysis of the crack growth rate with respect to the time step, the following exponential growth expression was obtained which shows the trend of crack growth during the fracturing process:(22)y=−37.65+0.056ex+1721.5/16402.1,R2=0.992.where y is the increase in number of cracks and x is the time step.


### 4.3. The Distribution Characteristic of Moment Magnitude for MCS Material

Figure7 shows a cloud map of the moment magnitude calculated at each individual space through the interpolation method. The numerical sample with a strain of 1.0‰ is not included in the moment magnitude distribution map because no crack formed inside the material.Figure 7
Moment magnitude cloud map under uniaxial compression.By comparing the coloured legend associated with the difference moment magnitude to the cloud map shown in Figure7, it was found that the moment magnitude was mainly distributed between ‒7.6 and ‒6.1 during the fracturing process. Prior to reaching a peak strain of 3.4‰ (including the stable crack growth stage), the moment magnitude cloud map is covered with light colours, that is, mostly light yellow and yellow. This indicates a small moment magnitude, small displacement between fracturing surfaces, and a low level of energy released during the structural fracturing process. As the strain further increases from the peak value of 3.4‰ until complete fracture (unstable crack growth stage), the moment magnitude cloud map is covered with dark colours with interconnected dark red blocks. These blocks exhibit typical features of the fracturing zone (e.g., F1 and F2). Thus, they indicate a high moment magnitude, large displacement between fracturing surfaces, and a high level of energy released during the unstable crack growth stage of the structural fracturing process. Magnifying the moment magnitude cloud map near F1 reveals a gradual decreasing trend of the moment magnitude, which is represented by a colour gradient from dark red, light red, light yellow, to yellow in the localised cloud map (shown as M in Figure 7).In summary, with increasing mismatching areas between the structural surfaces, an increase in the moment magnitude was observed in the MCS material where the fracture strength and energy release became significantly higher when approaching the fracture zone. The released energy propagated outwards around the source of the damage and became weaker with increasing distance. In particular, at certain locations far away from the macroscale fracturing zone, the local moment magnitude was relatively small (e.g., the N zone) together with a low level of released energy and destruction. Therefore, it is an inherently challenging task to detect any fracturing behaviour through acoustic emission at these regions. However, the moment magnitude was successfully captured at these locations using the simulated acoustic emission method. This implies that the acoustic emission simulation is a powerful tool to detect the local structural change in unknown regions.Finally, the fracturing probability associated with different crack numbers was analysed based on the moment magnitude results obtained from the acoustic emission event. The probability associated with the crack number follows a normal distribution with respect to the moment magnitude as shown in Figure8. The fracturing probability follows a Gaussian distribution with respect to the moment magnitude under different levels. In the end, it was found that the material is most prone to fracture when the moment magnitude lies between ‒6.8 and ‒6.6.Figure 8
Distribution of the moment magnitude probability in the acoustic emission event under uniaxial compression.


## 4.1. Calibration Results of MCS Material

According to the calibration method for mesoscopic contact parameters, uniaxial compression tests were performed using the initial mesoscopic parameter values given in section2.4. The numerical tests showed same loading rate in comparison to indoor experimental tests. The final calibrated parameters are summarized in Table 2. The mechanical parameters and testing results obtained from the numerical tests using calibrated parameters were compared with those obtained from experiments conducted under typical testing conditions (Table 3). Table 3 presents the comparison between the numerical and indoor experimental tests from which a close match of the mechanical parameters was found. Figure 5 shows the stress-strain curves and structural fracturing modes from both physical experiments and numerical simulations. The monitoring points A, B, C, D, and E shown on the numerical curve in Figure 5 correspond to strains of 1.0‰, 2.0‰, 3.0‰, 3.4‰, and 3.6‰, respectively. Among all these values, the peak strain is at 3.4‰. It can also be observed from Figure 5 that combined tension-shear fracturing is the primary mode of bond breaking in the numerical model. The numerical tests and physical tests have the following characteristics: (1) they show an identical developing trend although there is a small local deviation between the stress-stain curves. (2) The compressive strength and strain amplitude are basically the same. (3) The numerical tests revealed a similar trend between the macroscale cracking and the overall fracturing feature with minor local deviations.Table 2
Mesoscopic contact parameters.
Mesoscopic contact
E
∗ (N/m2)
E
¯
∗ (N/m2)
k
∗
k
¯
∗
µ
p
b
_
t
e
n (MPa)
p
b
_
c
o
h (MPa)
λ
Sand-sand
1.50 × 109
1.50 × 109
1.2
1.2
0.7
13.50 ± 1.0
15.50 ± 1.5
0.25–0.50
Sand-calcite
0.85 × 109
0.85 × 109
2.5
2.5
0.5
17.50 ± 1.0
24.50 ± 1.5
0.50
Calcite-calcite
1.25 × 109
1.25 × 109
2.5
2.5
0.5
14.50 ± 1.0
17.50 ± 1.5
0.50Table 3
Comparison of characteristic parameters.
Physical experiment
Numerical experiment
σ
p/MPa
4.61
4.80
E/GPa
1.32
1.42
ε
p/‰
3.5
3.4Figure 5
Comparison of stress-strain curves and structural fracturing modes from both physical experiments and numerical simulations.For the small local deviation between the stress-strain curves, the difference is derived from the mesostructural difference between the physical sample and the model sample. Although it is as close as possible to the real mesostructure, it is difficult to consider all the mesostructures of materials [20]. It is almost impossible to reconstruct a model in a numerical simulation sharing the same micro/mesoscale structure feature and physical properties as those of a physical testing sample [41]. For the difference in failure mode, it is more derived from the randomness of the model itself, which reflects the nonuniformity of the material. Due to the nonuniformity of the material itself, even the same material will not respond exactly the same under the same load. In addition, the macroscopic cracks in the numerical test are formed by microscopic cracks. Due to the discreteness of the particle contact in the numerical model, the generated microcracks are also discrete [42].Based on the above analysis results, the mesoscopic parameters used in general are reasonable and can be used to characterize the characteristic properties of MCS material under selected conditions.


## 4.2. The Fracturing Process and Crack Growth of MCS Material

Figure6 illustrates the fracturing process and crack growth of MCS material, where the coloured legend represents the contact force. Three stages have been illustrated in the fracturing process (Figure 6), namely, the stress-bearing stage, the stable crack growth stage, and the unstable crack growth stage.(1)
During the stress-bearing stage (A1), the MCS material was deformed internally and the internal contact force in the model was relatively small at around 0.5 kN. This level of contact force was insufficient to induce fracture of the structure. Consequently, no crack was observed in the material structure, and the strain increased up to 1.0‰.(2)
During the stable crack growth stage (A2) with further increase in the load, a stronger contact force is induced inside the material (e.g., a contact force yielding a strain of 2.0‰). Thus, once the contact force exceeded the bonding strength between the particles, a microcrack was formed inside the material. The number of cracks increased linearly during this stage.(3)
During the unstable crack growth stage, the contact force was distributed approximately between 0.5 and 2.5 kN. Owing to the significant change in the contact force, microcracks continued to force, expand, merge, and eventually turn into a macroscale fracturing zone inside the modelled material (e.g., the fracturing zone F1 with a strain of 3.6‰). Further, a sharp increase in the number and coverage of the cracks was observed, which expanded/grew/propagated nonlinearly with time with apparent bifurcations being observed at the fracturing zone. The bifurcations found in the fracturing process of the MCS material reflects the nonuniformity of the MCS material where the cracks started to expand in favourable paths with minimum resistive force under the applied stress. This implies that implementing a bonding strength distribution yields a more accurate result in numerical modelling.Figure 6
Distribution of contact force and crack growth during the fracturing process of MCS material.It is also worth mentioning that during the unstable crack growth stage when macroscale cracks start to form, a more concentrated stress was observed at some local cracking tip regions (e.g., the K region). However, the cracks did not expand further towards these areas because the dynamic stress intensification factor at the cracking tip was smaller than the minimum fracture toughness required to induce a crack. Previous studies on fracture mechanics [43] have shown that the expansion of a crack creates a stress field with the accumulation of elastic energy. Therefore, when the dynamic stress intensification factor at the cracking tip is smaller than the critical fracture toughness, the crack ceases to expand. Whereas when the dynamic stress intensification factor at the cracking tip exceeds the critical fracture toughness, the crack continues to expand. The simulation results obtained in this study are in accordance with the first case, i.e., cracks did not expand further towards the stress intensified regions.By performing regression analysis of the crack growth rate with respect to the time step, the following exponential growth expression was obtained which shows the trend of crack growth during the fracturing process:(22)y=−37.65+0.056ex+1721.5/16402.1,R2=0.992.where y is the increase in number of cracks and x is the time step.


## 4.3. The Distribution Characteristic of Moment Magnitude for MCS Material

Figure7 shows a cloud map of the moment magnitude calculated at each individual space through the interpolation method. The numerical sample with a strain of 1.0‰ is not included in the moment magnitude distribution map because no crack formed inside the material.Figure 7
Moment magnitude cloud map under uniaxial compression.By comparing the coloured legend associated with the difference moment magnitude to the cloud map shown in Figure7, it was found that the moment magnitude was mainly distributed between ‒7.6 and ‒6.1 during the fracturing process. Prior to reaching a peak strain of 3.4‰ (including the stable crack growth stage), the moment magnitude cloud map is covered with light colours, that is, mostly light yellow and yellow. This indicates a small moment magnitude, small displacement between fracturing surfaces, and a low level of energy released during the structural fracturing process. As the strain further increases from the peak value of 3.4‰ until complete fracture (unstable crack growth stage), the moment magnitude cloud map is covered with dark colours with interconnected dark red blocks. These blocks exhibit typical features of the fracturing zone (e.g., F1 and F2). Thus, they indicate a high moment magnitude, large displacement between fracturing surfaces, and a high level of energy released during the unstable crack growth stage of the structural fracturing process. Magnifying the moment magnitude cloud map near F1 reveals a gradual decreasing trend of the moment magnitude, which is represented by a colour gradient from dark red, light red, light yellow, to yellow in the localised cloud map (shown as M in Figure 7).In summary, with increasing mismatching areas between the structural surfaces, an increase in the moment magnitude was observed in the MCS material where the fracture strength and energy release became significantly higher when approaching the fracture zone. The released energy propagated outwards around the source of the damage and became weaker with increasing distance. In particular, at certain locations far away from the macroscale fracturing zone, the local moment magnitude was relatively small (e.g., the N zone) together with a low level of released energy and destruction. Therefore, it is an inherently challenging task to detect any fracturing behaviour through acoustic emission at these regions. However, the moment magnitude was successfully captured at these locations using the simulated acoustic emission method. This implies that the acoustic emission simulation is a powerful tool to detect the local structural change in unknown regions.Finally, the fracturing probability associated with different crack numbers was analysed based on the moment magnitude results obtained from the acoustic emission event. The probability associated with the crack number follows a normal distribution with respect to the moment magnitude as shown in Figure8. The fracturing probability follows a Gaussian distribution with respect to the moment magnitude under different levels. In the end, it was found that the material is most prone to fracture when the moment magnitude lies between ‒6.8 and ‒6.6.Figure 8
Distribution of the moment magnitude probability in the acoustic emission event under uniaxial compression.


## 5. Conclusion

In this study, a detailed simulation analysis method for MCS materials numerically using the PFC 2D software has been explored. Three major conclusions have been drawn from the study:(1)
In consideration of bonding characteristic of particles, concentration and proportion of compositional structure, particle size distribution of sand particles, and characterization of crystal structure and morphological characteristics, a simulation method for microbial cemented sand based on the Two-Dimensional Particle Flow Code (PFC 2D) has been developed, mainly consisting of fine identification of mesoscopic contact model for mesoscopic particles of microbial cemented sand, development of morphological algorithm for irregular crystal particles, and particle size distribution. XCT scanning test has shown that the calcite crystal generated with morphological algorithm for irregular crystal particles has better consistence with the morphology of real crystal particles, which means the method may really reflect physical characteristics of cemented sand.(2)
Based on the relation between macroscopic and mesoscopic theory of discrete element, the method for setting of mesoscopic contact parameters of different particles within microbial cemented sand has been proposed. As per uniaxial compression test results, mesoscopic contact parameters in selected operating conditions have been determined. In combination with physical tests in such operating conditions, the method has been demonstrated to effectively reflect mechanical characteristics of cemented sand.(3)
The acoustic emission method based on moment tensor theory has been developed for discrete element analysis on cemented sand fracture. The moment magnitude (MW) has been used to reflect fracture strength of microbial cemented sand materials, resulting in quantitative analysis of cemented sand fracture. Analysis shows that the probability of fracture of the material is greater when the MW ranges between −6.8 and −6.6. The relationship between fracture probability at each of different MWs and the MW follows the Gaussian curve.The research is of much significance for the simulation of damage process of microbial cemented sand under the action of real load and mesoscopic analysis of microbial cemented sand. Further research will consider conducting an in-depth study with respect to fine characterization of crystal structure outline and boundary and optimisation of calibration method for mesoscopic contact parameters.



---

*Source: 7173414-2019-03-11.xml*