# Thermal Cracking Analysis of Microbial Cemented Sand under Various Strains Based on the DEM

Authors: Yang Tang; Guobin Xu; Yue Yan; Dengfeng Fu; Chunlai Qu; Zilong Li; Evance Chaima
Journal: Advances in Materials Science and Engineering
 (2018)

---

## Abstract

Microbial-induced calcite precipitation (MICP) is a novel ground improvement method to increase the strength and stiffness of sand. However, the influences of temperature load on the internal microstructure of microbial cemented sand (MCS) material under the experimented strain have always been a key concern for the extensive application. Three kinds of experiments, X-ray diffraction (XRD), X-ray computed tomography (XCT), and scanning electron microscopy (SEM), were conducted to explore the composition, shape, and bonding characteristics of physical assemblies in this paper. A precision DEM modelling of MCS, mainly composed of irregular particle modelling and a mesoparameter calibration algorithm, has been proposed for the thermal cracking analysis under various strains (i.e., 1.0‰–3.0‰). Research results indicate that three kinds of bonding (that is sand-calcite, calcite-calcite, and sand-sand) are present in the MCS material. The application of temperature has a superposition effect on the damage of MCS material with increasing strain. Moreover, as the heating duration gradually increases, the effect of thermal rupture produces a distinct quiet period. The length of thermal cracks in the transverse direction increases throughout the heating process.


---

## Body

## 1. Introduction

Microbial-induced calcite precipitation (MICP) is a novel ground improvement method that is used to increase the strength and stiffness of sand using natural biogeochemical processes. Currently, MICP technology is being gradually applied in the field of engineering [1–3]. As a reinforcement technology, microbial cemented sand (MCS) structure based on MICP technology will inevitably encounter an extremely harsh geological environment, such as high stress and high temperature, which will result in the failure surface at the internal yield point of the internal structure as the strain increases [4]. If the deformed microbial cement structure is exposed to a temperature load, the internal structure of the MCS material might be affected. At the mesolevel, various internal microdefects may occur and develop; finally, it will change the basic mechanical properties of the MCS structure. It is difficult to observe the development of the microcrack network of the structure under the influence of temperature from the macroscopic physical experiments.The discrete element method (DEM) models are well-suited for providing the necessary data to analyse the interactions in granular mixtures and are beneficial for improving our understanding of thermal cracking behaviour [5]. Recently, many studies have used the DEM method to perform thermal cracking analysis of rock-like materials [6–12]. However, the analysis of the fracture mechanism of deformed cemented structures under temperature load and, in particular, MCS material is still lacking, which is mainly due to multiple interactions along the interface among different particles (i.e., sand and calcite) in the MCS material [13]. This behaviour leads to the complexity of establishing both mesoscale bonding model and coupled thermomechanical model.This study is aimed at providing an improved understanding of the thermal cracking behaviour of MCS material on a microscopic scale, which is not easily elucidated by macroscopically experiments. The remainder of the following sections has been organized as follows: First, the physical experimental process of MCS has been introduced briefly. The results obtained using X-ray diffraction (XRD), X-ray computed tomography (XCT), and scanning electron microscopy (SEM) were analysed to explore the material physical characterizations at different scales. Second, the observed evidence was used to set up the mesoparameters of the MCS material model with the aid of the adaptive differential evolution (ADE) algorithm; the interactions of sand-sand, sand-calcite, and calcite-calcite were carefully considered. Based on this work, the extensive coupled thermomechanical model was finally constructed to model the thermal cracking behaviour of the MCS material, with five different strains (i.e., 1.0‰, 1.5‰, 2.0‰, 2.5‰, and 3.0‰) and various heating durations (i.e., 0.02 s, 0.05 s, 0.10 s, 0.25 s, and 0.50 s) being considered.


## 2. Experimental Setup and Experimental Characterization Analysis of Physical Assemblies

### 2.1. Experimental Setup

The experimental MICP treatment process and uniaxial compressive strength and mineralogical and textural analysis tests of the MCS material are summarized below.


#### 2.1.1. Sample Preparation

The basic physical properties of standard sand used in this study are shown in Table1. The paper briefly introduces the sample preparation method; however, the detailed method has been described by Xu et al. [14].Table 1
Physical properties of the sand used in this study.
Grain shape
Apparent density (g/cm3)
Water content (%)
Porosity (%)
Chemical composition
Grain-size distribution (mm)
Round
2.65
0
40.5
>98% SiO2
0.075∼0.60A visible light spectrophotometer (WFJ2000, UNICO) was used to detect the number of microorganisms. The wavelength was 600 nm, and the measured value was OD600. A cementation solution of urea-CaCl2 was used where the urea was the nitrogen and energy source of microbial growth and CaCl2 was the calcium source in the process of MICP. The test conditions selected in this paper were OD600 = 1.50 and [Ca2+] = 0.75 mol/l.Five additions of cementing solutions and one of the bacterial liquids were considered one cycle for the entire cementation process, and two cycles were performed for the test sample.


#### 2.1.2. Uniaxial Compressive Strength (UCS) Test

The UCS testing of the MCS material was conducted using a microcomputer-controlled electrohydraulic servo loading testing machine (TAW-2000) (shown in Figure1). The maximum loading capacity of TAW-2000 was 2000 kN. The uniaxial compressive strength test was set to be displacement controlled at a loading rate of 0.05 mm/min until failure occurred. A Micro-II digital acoustic emission (AE) system was used simultaneously to identify fracture and deformation events caused by damage to the sample.Figure 1
The UCS test of MCS material.


#### 2.1.3. Mineralogical and Textural Analysis Test

The mineralogical and textural analysis test consisted of three parts: XRD, XCT, and SEM. The cemented sand body on the surface of the column was examined using a scanning electron microscopy (SEM) and X-ray diffraction (XRD). The paper gives a brief introduction to the methods; however, detailed methods and relative parameters have been described by Xu et al. [14].A micro-XCT instrument (diondo d5) with the resolution of 27µm, a voltage of 180 kV, and tube current of 0.18 mA was used, and the exposure time of the detector was 1000 ms in this study. The scanning process (shown in Figure 2) consisted of four parts: (a) test sample selection, (b) placing and adjusting the sample position, (c) setting the parameters and scanning the sample, and (d) feature extraction.Figure 2
The micro-XCT experimental process. (a) Test sample selection, (b) placing and adjusting the sample position, (c) setting the parameters and scanning the sample, and (d) feature extraction.
(a)
(b)
(c)
(d)


### 2.2. Experimental Characterization Analysis of Physical Assemblies

#### 2.2.1. Composition of Physical Assemblies

The XRD patterns of the pure sand and cemented sand are shown in Figure3. According to the XRD patterns, the pure sand was composed of quartz, and the mineralized cemented sand was composed of quartz and calcite phases.Figure 3
XRD patterns of the pure sand and the cemented sand.
(a)
(b)


#### 2.2.2. Shapes of Physical Assemblies

The XCT experimental result is shown in Figure4. Obviously, the calcite crystals in the MCS samples were irregularly shaped and randomly distributed, and most of the particles were elongated and polygonal.Figure 4
The micro-CT experimental result. (a) Cross Section. (b) Longitudinal section.


#### 2.2.3. Bonding Characteristics of Physical Assemblies

According to the results of the microscopic images shown in Figure5, it was found that three kinds of bonding were inside the cemented sand: sand-sand bonding, sand-calcite bonding, and calcite-calcite bonding. The sand-sand bonding exhibited an inhomogeneous bonding between particles (Figure 5(b2)), with some visible holes at the interfaces of the bonding, while a tight bonding was exhibited for the sand-calcite particles and calcite-calcite particles (Figure 5(b1)). Considering the aforementioned bonding characteristics of the MCS material, it is necessary to take into account the different bonding characteristics and distribution of the MCS material when setting the mesoparameters of the numerical model.Figure 5
SEM micrographs and cementation formation in the physical specimen.


## 2.1. Experimental Setup

The experimental MICP treatment process and uniaxial compressive strength and mineralogical and textural analysis tests of the MCS material are summarized below.


### 2.1.1. Sample Preparation

The basic physical properties of standard sand used in this study are shown in Table1. The paper briefly introduces the sample preparation method; however, the detailed method has been described by Xu et al. [14].Table 1
Physical properties of the sand used in this study.
Grain shape
Apparent density (g/cm3)
Water content (%)
Porosity (%)
Chemical composition
Grain-size distribution (mm)
Round
2.65
0
40.5
>98% SiO2
0.075∼0.60A visible light spectrophotometer (WFJ2000, UNICO) was used to detect the number of microorganisms. The wavelength was 600 nm, and the measured value was OD600. A cementation solution of urea-CaCl2 was used where the urea was the nitrogen and energy source of microbial growth and CaCl2 was the calcium source in the process of MICP. The test conditions selected in this paper were OD600 = 1.50 and [Ca2+] = 0.75 mol/l.Five additions of cementing solutions and one of the bacterial liquids were considered one cycle for the entire cementation process, and two cycles were performed for the test sample.


### 2.1.2. Uniaxial Compressive Strength (UCS) Test

The UCS testing of the MCS material was conducted using a microcomputer-controlled electrohydraulic servo loading testing machine (TAW-2000) (shown in Figure1). The maximum loading capacity of TAW-2000 was 2000 kN. The uniaxial compressive strength test was set to be displacement controlled at a loading rate of 0.05 mm/min until failure occurred. A Micro-II digital acoustic emission (AE) system was used simultaneously to identify fracture and deformation events caused by damage to the sample.Figure 1
The UCS test of MCS material.


### 2.1.3. Mineralogical and Textural Analysis Test

The mineralogical and textural analysis test consisted of three parts: XRD, XCT, and SEM. The cemented sand body on the surface of the column was examined using a scanning electron microscopy (SEM) and X-ray diffraction (XRD). The paper gives a brief introduction to the methods; however, detailed methods and relative parameters have been described by Xu et al. [14].A micro-XCT instrument (diondo d5) with the resolution of 27µm, a voltage of 180 kV, and tube current of 0.18 mA was used, and the exposure time of the detector was 1000 ms in this study. The scanning process (shown in Figure 2) consisted of four parts: (a) test sample selection, (b) placing and adjusting the sample position, (c) setting the parameters and scanning the sample, and (d) feature extraction.Figure 2
The micro-XCT experimental process. (a) Test sample selection, (b) placing and adjusting the sample position, (c) setting the parameters and scanning the sample, and (d) feature extraction.
(a)
(b)
(c)
(d)


## 2.1.1. Sample Preparation

The basic physical properties of standard sand used in this study are shown in Table1. The paper briefly introduces the sample preparation method; however, the detailed method has been described by Xu et al. [14].Table 1
Physical properties of the sand used in this study.
Grain shape
Apparent density (g/cm3)
Water content (%)
Porosity (%)
Chemical composition
Grain-size distribution (mm)
Round
2.65
0
40.5
>98% SiO2
0.075∼0.60A visible light spectrophotometer (WFJ2000, UNICO) was used to detect the number of microorganisms. The wavelength was 600 nm, and the measured value was OD600. A cementation solution of urea-CaCl2 was used where the urea was the nitrogen and energy source of microbial growth and CaCl2 was the calcium source in the process of MICP. The test conditions selected in this paper were OD600 = 1.50 and [Ca2+] = 0.75 mol/l.Five additions of cementing solutions and one of the bacterial liquids were considered one cycle for the entire cementation process, and two cycles were performed for the test sample.


## 2.1.2. Uniaxial Compressive Strength (UCS) Test

The UCS testing of the MCS material was conducted using a microcomputer-controlled electrohydraulic servo loading testing machine (TAW-2000) (shown in Figure1). The maximum loading capacity of TAW-2000 was 2000 kN. The uniaxial compressive strength test was set to be displacement controlled at a loading rate of 0.05 mm/min until failure occurred. A Micro-II digital acoustic emission (AE) system was used simultaneously to identify fracture and deformation events caused by damage to the sample.Figure 1
The UCS test of MCS material.


## 2.1.3. Mineralogical and Textural Analysis Test

The mineralogical and textural analysis test consisted of three parts: XRD, XCT, and SEM. The cemented sand body on the surface of the column was examined using a scanning electron microscopy (SEM) and X-ray diffraction (XRD). The paper gives a brief introduction to the methods; however, detailed methods and relative parameters have been described by Xu et al. [14].A micro-XCT instrument (diondo d5) with the resolution of 27µm, a voltage of 180 kV, and tube current of 0.18 mA was used, and the exposure time of the detector was 1000 ms in this study. The scanning process (shown in Figure 2) consisted of four parts: (a) test sample selection, (b) placing and adjusting the sample position, (c) setting the parameters and scanning the sample, and (d) feature extraction.Figure 2
The micro-XCT experimental process. (a) Test sample selection, (b) placing and adjusting the sample position, (c) setting the parameters and scanning the sample, and (d) feature extraction.
(a)
(b)
(c)
(d)


## 2.2. Experimental Characterization Analysis of Physical Assemblies

### 2.2.1. Composition of Physical Assemblies

The XRD patterns of the pure sand and cemented sand are shown in Figure3. According to the XRD patterns, the pure sand was composed of quartz, and the mineralized cemented sand was composed of quartz and calcite phases.Figure 3
XRD patterns of the pure sand and the cemented sand.
(a)
(b)


### 2.2.2. Shapes of Physical Assemblies

The XCT experimental result is shown in Figure4. Obviously, the calcite crystals in the MCS samples were irregularly shaped and randomly distributed, and most of the particles were elongated and polygonal.Figure 4
The micro-CT experimental result. (a) Cross Section. (b) Longitudinal section.


### 2.2.3. Bonding Characteristics of Physical Assemblies

According to the results of the microscopic images shown in Figure5, it was found that three kinds of bonding were inside the cemented sand: sand-sand bonding, sand-calcite bonding, and calcite-calcite bonding. The sand-sand bonding exhibited an inhomogeneous bonding between particles (Figure 5(b2)), with some visible holes at the interfaces of the bonding, while a tight bonding was exhibited for the sand-calcite particles and calcite-calcite particles (Figure 5(b1)). Considering the aforementioned bonding characteristics of the MCS material, it is necessary to take into account the different bonding characteristics and distribution of the MCS material when setting the mesoparameters of the numerical model.Figure 5
SEM micrographs and cementation formation in the physical specimen.


## 2.2.1. Composition of Physical Assemblies

The XRD patterns of the pure sand and cemented sand are shown in Figure3. According to the XRD patterns, the pure sand was composed of quartz, and the mineralized cemented sand was composed of quartz and calcite phases.Figure 3
XRD patterns of the pure sand and the cemented sand.
(a)
(b)


## 2.2.2. Shapes of Physical Assemblies

The XCT experimental result is shown in Figure4. Obviously, the calcite crystals in the MCS samples were irregularly shaped and randomly distributed, and most of the particles were elongated and polygonal.Figure 4
The micro-CT experimental result. (a) Cross Section. (b) Longitudinal section.


## 2.2.3. Bonding Characteristics of Physical Assemblies

According to the results of the microscopic images shown in Figure5, it was found that three kinds of bonding were inside the cemented sand: sand-sand bonding, sand-calcite bonding, and calcite-calcite bonding. The sand-sand bonding exhibited an inhomogeneous bonding between particles (Figure 5(b2)), with some visible holes at the interfaces of the bonding, while a tight bonding was exhibited for the sand-calcite particles and calcite-calcite particles (Figure 5(b1)). Considering the aforementioned bonding characteristics of the MCS material, it is necessary to take into account the different bonding characteristics and distribution of the MCS material when setting the mesoparameters of the numerical model.Figure 5
SEM micrographs and cementation formation in the physical specimen.


## 3. Numerical Specimen Generation and Parameter Calibration

In the current study, although the PFC 3D version could be employed for solving 3D problems, the PFC 3D simulations required very long computational times, and the boundary conditions were difficult to control [15]. On the contrary, the 2D simulations have the advantage of explaining key phenomena and mechanisms in less time than 3D modelling. In addition, the visual deformation patterns and force chains are easily captured by 2D simulations. Many previous studies [16–19] illustrated that the 2D simulations can adequately capture various complex mechanical features of the materials when considering the shape’s effects. Therefore, the numerical simulation for biocemented sand specimens was carried out using a PFC 2D model in this research.


### 3.1. Numerical Specimen Generation

The numerical specimen generation process was conducted through the following steps.


#### 3.1.1. Contact Constitutive Model Selection

According to the arrangement of MICP in the microscopic images in Figure5, the bond pattern was close to the parallel-bond model (Figure 6) of DEM [20]. The parallel-bond model (PBM) was depicted as a rectangle of cement-like material, which could transmit both force and the moment between the particles [21]. Moreover, the PBM has been successfully used to explore the mechanical behaviour of MCS material [22]. As such, the PBM was selected for the numerical study in this research.Figure 6
Parallel-bond model.


#### 3.1.2. Calcite and Sand Particles Modelling

As in the aforementioned XCT scanning results, the calcite crystals in the MCS samples were irregularly shaped and randomly distributed, and most of the particles were elongated and polygonal. There were a larger number of physical contacts for irregular particles than that for circular particles [23]. Moreover, the particle shape played a key role in the behaviour of the material [24]. To overcome these problems, the appropriate geometrical shape of irregular particles should be modelled. This study adopted an irregular generation algorithm to determine the shape effect of the MICP material. The irregular polygon generation process of calcite particle is shown in Figure 7, which was described in detail by Tang et al. [13].Figure 7
Irregular polygon generation algorithm process of calcite crystals in the MCS samples [13]. (a) Generate the initial cell; (b) generate the characteristic circle; (c) determine the vertex of irregular polygon; (d) irregular polygon generation.
(a)
(b)
(c)
(d)For the sand particles, the previous analysis results [14] have shown that the particles were mostly circular, and the particles for this particle size have been characterized by circular particles [25]. For this reason, the sand particles were simulated using circles.


#### 3.1.3. Numerical Specimen Generation

Using the parallel-bond model, a 2D numerical model with the height of 100 mm and the width of 50 mm was first built to simulate the laboratory behaviour of a cylindrical specimen of biocemented sand with the same size as that of the numerical model. According to the research results [13], the calcite content of cementation sand was 3% in the 2D model, and the numerical model consisted of 3% calcite and 97% quartz by area, combined with the generation method of calcite.Subsequently, the numerical particles were produced in the numerical specimen model based on the grain-size distribution (GSD) of the real particles. The numbers of quartz and calcite particles were 165244 and 33955, respectively. The GSD for real particles and numerical particles are shown in Figure8; a strongly consistent GSD was exhibited for real particles and numerical particles.Figure 8
Grain-size distribution of experimental and numerical particles.


### 3.2. Parameter Calibration

Realistic mesoparameters were required when using DEM to calculate the numerical value. However, the mesoparameters of the numerical model were generally adjusted by trial and error until the numerically simulated results corresponded to the macroparameters of the real specimens. This method usually requires a large number of trial calculations, and the calculation period is long; thus, the human error can be significant. Compared to the trial-and-error method, adaptive differential evolution algorithms (ADEs) have been widely used in many fields and have the advantages of simple control parameters, global search ability, strong robustness, fast convergence speed, and precision characteristics [26–28]. They would also be suitable for solving the problem of parameter inversion, especially for the inversion of discontinuous, nonlinear, nondifferentiable, and analytic functions.


#### 3.2.1. Calibration Method

In this paper, the mesoparameters of MCS material were calibrated using the ADE method. The calibrated mesoparameters includeE∗, E¯∗, k∗, k¯∗, fric, pb_ten, pb_coh, and R¯.The SEM characterisation of the material (shown in Figure5) makes it necessary to consider the different bonding characteristics and discreteness of MCS material when setting mesoparameters. According to the results of Feng et al. [22], different bonding characteristics can be characterized by different bonding radii. When the parallel bonding radius, R¯ (m), is determined, the particle flow theory suggested that the mesoscale bonding radius, R¯, was related to the radial coefficient, λ. For the sand particle-calcite and calcite-calcite bonds, λ was assigned a fixed value of 0.5 because these bonds were relatively strong according to the SEM analysis results. For the sand particle-sand particle bond, a random number to vary the parallel-bond radial coefficient λ between 0.25 and 0.5 was introduced based on the previous PFC study on microbial cemented sand column by Feng et al. [22]. Thus, the bonding radius was determined and is not included in the parameter inversion procedure.The mesoparametersM of the MCS material used in the simulation can be expressed as(1)M=E∗,E¯∗,k∗,k¯∗,fric,pb_ten,pb_cohT,where E∗ is the effective modulus (force/area), E¯∗ is the effective modulus of the bond (stress), k∗ is the normal-to-shear stiffness ratio, k¯∗ is the normal-to-shear stiffness ratio of the bond, fric is the friction coefficient, pb_ten is the tensile strength (stress), and pb_coh is the cohesion (stress).The adaptive differential evolution algorithm used to invert the mesoparameters is shown in Figure9, and the overall process is described in the following three steps.Figure 9
The realisation process of the inversion algorithm of mesoparameters.(1) Construction of the Characteristic Parameters. Using the PFC simulating uniaxial compression test, Nii≥20 groups of mesoparameters Mjj∈1,Ni were constructed, and the stress-strain and crack data corresponding to these parameters were obtained. It was noted from the laboratory test results of elastic modulus size of the MCS material that the elastic modulus of sand particles is relatively large and that of calcite is relatively small [29]. Thus, the initial elastic modulus parameters were set according to this size. Subsequently, based on the stress-strain and crack data, the characteristic parameters were calculated using MATLAB (R2015b) software. The obtained characteristic parameters were used as the initial population samples NPjj∈1,Ni for the adaptive differential evolution algorithm. The characteristic parameters KN−jj∈1,Ni included the peak stress (σc−j∗), peak strain (εc−j∗), elastic modulus (Ej), and the crack initiation stress (σn−j). Note that the crack initiation stress (σn−j) was identified as the axial stress where the number of cracks reached 1% of the total number of cracks at peak stress [20]:(2)Mj=E∗,E¯∗,k∗,k¯∗,fric,pb_ten,pb_cohjT⟶KN−j=σc−j∗,εc−j∗,Ej,σn−j.(2) Determination of the Target Inversion Value from the Uniaxial Compression Data. The experimental uniaxial compression tests provided the real stress-strain behaviour and AE monitoring data. According to the stress-strain relationship and the definition of the cracking stress, the peak stress σc, peak strain εc, elastic modulus (E), and crack initiation stress (σ1−p) of the MCS material were obtained and used as the target values for inversion. It should be noted that AE monitors low-energy high-frequency seismic events resulting from crack initiation, so it can be used for the identification of σ1−p [30]:(3)F=σc,εc,E,σ1−p.(3) Construction of the Fitness Function of the Inversion Algorithm. According to the combination of different mesoparameters Mjj∈1,Ni, the study focused on the relatively small errors between the simulated and experimental characteristic parameters KN−jj∈1,Ni, where F is the target value. The objective function of the inversion of the mesoparameters can be defined as(4)Gk=KN−j−F.Furthermore, the fitness function was constructed as follows:(5)fk=1Gk.It was found that the larger individual fitness values resulted in smaller objective function values; the mesoparameters reflected the real MCS material properties.


#### 3.2.2. Calibration Results and Comparison of Calibration Methods

According to the above calibration method, the mesoparameters with small errors were observed, as shown in Table2. The uniaxial stress-strain curves and failure characteristics calculated from the numerical simulation are compared to that of the experimental data in Figure 10.Table 2
Mesoparameters of material.
Contact properties
E
∗ (N/m2)
E
¯
∗ (N/m2)
k
∗
k
¯
∗
fric
pb_ten (MPa)
pb_coh (MPa)
λ
Sand-sand
1.50 × 109
1.50 × 109
1.2
1.2
0.7
13.05 ± 1.0
15.05 ± 1.5
0.25∼0.50
Sand-calcite
0.85 × 109
0.85 × 109
2.5
2.5
0.5
17.05 ± 1.0
24.05 ± 1.5
0.50
Calcite-calcite
1.25 × 109
1.25 × 109
2.5
2.5
0.5
14.05 ± 1.0
17.05 ± 1.5
0.50Figure 10
Comparison of numerical simulation and experimental stress-strain curves and failure modes.The results shown in Figure10 were found by comparing the numerical results to the results of physical experiments, with 1.95% and 4.0% as the relative errors in the peak stress and peak strain, respectively. Additionally, the numerical tests revealed a similar trend between the macroscale cracking and the overall fracturing feature with minor local deviations. The deviations are primarily due to the discrete connection of the microscale cracks in the macroscale cracks, which are inherently induced by the discrete contacts between particles in the model developed in this study. Meanwhile, some simulation results differ from the testing results; specifically, the elastic module obtained in each case differs from the others. This difference is justified by the model reconstruction of microscale particle structure and the related physical properties. In a numerical simulation, it is almost impossible to reconstruct a model sharing the same mesoscale structure feature and physical properties as those of a physical testing sample [31]. In view of these factors, the mesoscale parameters used in this study are appropriate for characterising the MCS material properties.Currently, the trial-and-error method has been the most commonly used calibration method for selecting appropriate mesoscale parameters [32, 33]. The stress-strain curve corresponding to the mesoparameters obtained using the trial-and-error method was plotted on the same graph (Figure 10) in order to compare the results of ADE and the trial-and-error method. Figure 10 clearly illustrates that the numerical curve obtained by the ADE method is similar to the curve obtained by Tang et al. [13] using the trial-and-error method. By comparison with the mesoparameters obtained by Tang et al. [13], most mesoparameters are almost the same, except for pb_ten and pb_coh.The mesoparameters used for the trial-and-error method usually have to be assigned according to the mesoscopic properties of the material. For relatively simple contact types, the trial-and-error time may be relatively short. However, when the material involves multiple contact types, the trial-and-error process is extremely complicated and takes more time. Compared to the trial-and-error method, ADE is able to create an initial database and obtains target parameter values through parameter inversion. The entire process only takes a certain amount of time during the initial database setup process. Once the initial database is established, corresponding mesoscopic parameter values can be quickly obtained. Furthermore, the established database can provide initial data for the subsequent acquisition of mesoparameters of related materials with similar mechanical behaviour.


## 3.1. Numerical Specimen Generation

The numerical specimen generation process was conducted through the following steps.


### 3.1.1. Contact Constitutive Model Selection

According to the arrangement of MICP in the microscopic images in Figure5, the bond pattern was close to the parallel-bond model (Figure 6) of DEM [20]. The parallel-bond model (PBM) was depicted as a rectangle of cement-like material, which could transmit both force and the moment between the particles [21]. Moreover, the PBM has been successfully used to explore the mechanical behaviour of MCS material [22]. As such, the PBM was selected for the numerical study in this research.Figure 6
Parallel-bond model.


### 3.1.2. Calcite and Sand Particles Modelling

As in the aforementioned XCT scanning results, the calcite crystals in the MCS samples were irregularly shaped and randomly distributed, and most of the particles were elongated and polygonal. There were a larger number of physical contacts for irregular particles than that for circular particles [23]. Moreover, the particle shape played a key role in the behaviour of the material [24]. To overcome these problems, the appropriate geometrical shape of irregular particles should be modelled. This study adopted an irregular generation algorithm to determine the shape effect of the MICP material. The irregular polygon generation process of calcite particle is shown in Figure 7, which was described in detail by Tang et al. [13].Figure 7
Irregular polygon generation algorithm process of calcite crystals in the MCS samples [13]. (a) Generate the initial cell; (b) generate the characteristic circle; (c) determine the vertex of irregular polygon; (d) irregular polygon generation.
(a)
(b)
(c)
(d)For the sand particles, the previous analysis results [14] have shown that the particles were mostly circular, and the particles for this particle size have been characterized by circular particles [25]. For this reason, the sand particles were simulated using circles.


### 3.1.3. Numerical Specimen Generation

Using the parallel-bond model, a 2D numerical model with the height of 100 mm and the width of 50 mm was first built to simulate the laboratory behaviour of a cylindrical specimen of biocemented sand with the same size as that of the numerical model. According to the research results [13], the calcite content of cementation sand was 3% in the 2D model, and the numerical model consisted of 3% calcite and 97% quartz by area, combined with the generation method of calcite.Subsequently, the numerical particles were produced in the numerical specimen model based on the grain-size distribution (GSD) of the real particles. The numbers of quartz and calcite particles were 165244 and 33955, respectively. The GSD for real particles and numerical particles are shown in Figure8; a strongly consistent GSD was exhibited for real particles and numerical particles.Figure 8
Grain-size distribution of experimental and numerical particles.


## 3.1.1. Contact Constitutive Model Selection

According to the arrangement of MICP in the microscopic images in Figure5, the bond pattern was close to the parallel-bond model (Figure 6) of DEM [20]. The parallel-bond model (PBM) was depicted as a rectangle of cement-like material, which could transmit both force and the moment between the particles [21]. Moreover, the PBM has been successfully used to explore the mechanical behaviour of MCS material [22]. As such, the PBM was selected for the numerical study in this research.Figure 6
Parallel-bond model.


## 3.1.2. Calcite and Sand Particles Modelling

As in the aforementioned XCT scanning results, the calcite crystals in the MCS samples were irregularly shaped and randomly distributed, and most of the particles were elongated and polygonal. There were a larger number of physical contacts for irregular particles than that for circular particles [23]. Moreover, the particle shape played a key role in the behaviour of the material [24]. To overcome these problems, the appropriate geometrical shape of irregular particles should be modelled. This study adopted an irregular generation algorithm to determine the shape effect of the MICP material. The irregular polygon generation process of calcite particle is shown in Figure 7, which was described in detail by Tang et al. [13].Figure 7
Irregular polygon generation algorithm process of calcite crystals in the MCS samples [13]. (a) Generate the initial cell; (b) generate the characteristic circle; (c) determine the vertex of irregular polygon; (d) irregular polygon generation.
(a)
(b)
(c)
(d)For the sand particles, the previous analysis results [14] have shown that the particles were mostly circular, and the particles for this particle size have been characterized by circular particles [25]. For this reason, the sand particles were simulated using circles.


## 3.1.3. Numerical Specimen Generation

Using the parallel-bond model, a 2D numerical model with the height of 100 mm and the width of 50 mm was first built to simulate the laboratory behaviour of a cylindrical specimen of biocemented sand with the same size as that of the numerical model. According to the research results [13], the calcite content of cementation sand was 3% in the 2D model, and the numerical model consisted of 3% calcite and 97% quartz by area, combined with the generation method of calcite.Subsequently, the numerical particles were produced in the numerical specimen model based on the grain-size distribution (GSD) of the real particles. The numbers of quartz and calcite particles were 165244 and 33955, respectively. The GSD for real particles and numerical particles are shown in Figure8; a strongly consistent GSD was exhibited for real particles and numerical particles.Figure 8
Grain-size distribution of experimental and numerical particles.


## 3.2. Parameter Calibration

Realistic mesoparameters were required when using DEM to calculate the numerical value. However, the mesoparameters of the numerical model were generally adjusted by trial and error until the numerically simulated results corresponded to the macroparameters of the real specimens. This method usually requires a large number of trial calculations, and the calculation period is long; thus, the human error can be significant. Compared to the trial-and-error method, adaptive differential evolution algorithms (ADEs) have been widely used in many fields and have the advantages of simple control parameters, global search ability, strong robustness, fast convergence speed, and precision characteristics [26–28]. They would also be suitable for solving the problem of parameter inversion, especially for the inversion of discontinuous, nonlinear, nondifferentiable, and analytic functions.


### 3.2.1. Calibration Method

In this paper, the mesoparameters of MCS material were calibrated using the ADE method. The calibrated mesoparameters includeE∗, E¯∗, k∗, k¯∗, fric, pb_ten, pb_coh, and R¯.The SEM characterisation of the material (shown in Figure5) makes it necessary to consider the different bonding characteristics and discreteness of MCS material when setting mesoparameters. According to the results of Feng et al. [22], different bonding characteristics can be characterized by different bonding radii. When the parallel bonding radius, R¯ (m), is determined, the particle flow theory suggested that the mesoscale bonding radius, R¯, was related to the radial coefficient, λ. For the sand particle-calcite and calcite-calcite bonds, λ was assigned a fixed value of 0.5 because these bonds were relatively strong according to the SEM analysis results. For the sand particle-sand particle bond, a random number to vary the parallel-bond radial coefficient λ between 0.25 and 0.5 was introduced based on the previous PFC study on microbial cemented sand column by Feng et al. [22]. Thus, the bonding radius was determined and is not included in the parameter inversion procedure.The mesoparametersM of the MCS material used in the simulation can be expressed as(1)M=E∗,E¯∗,k∗,k¯∗,fric,pb_ten,pb_cohT,where E∗ is the effective modulus (force/area), E¯∗ is the effective modulus of the bond (stress), k∗ is the normal-to-shear stiffness ratio, k¯∗ is the normal-to-shear stiffness ratio of the bond, fric is the friction coefficient, pb_ten is the tensile strength (stress), and pb_coh is the cohesion (stress).The adaptive differential evolution algorithm used to invert the mesoparameters is shown in Figure9, and the overall process is described in the following three steps.Figure 9
The realisation process of the inversion algorithm of mesoparameters.(1) Construction of the Characteristic Parameters. Using the PFC simulating uniaxial compression test, Nii≥20 groups of mesoparameters Mjj∈1,Ni were constructed, and the stress-strain and crack data corresponding to these parameters were obtained. It was noted from the laboratory test results of elastic modulus size of the MCS material that the elastic modulus of sand particles is relatively large and that of calcite is relatively small [29]. Thus, the initial elastic modulus parameters were set according to this size. Subsequently, based on the stress-strain and crack data, the characteristic parameters were calculated using MATLAB (R2015b) software. The obtained characteristic parameters were used as the initial population samples NPjj∈1,Ni for the adaptive differential evolution algorithm. The characteristic parameters KN−jj∈1,Ni included the peak stress (σc−j∗), peak strain (εc−j∗), elastic modulus (Ej), and the crack initiation stress (σn−j). Note that the crack initiation stress (σn−j) was identified as the axial stress where the number of cracks reached 1% of the total number of cracks at peak stress [20]:(2)Mj=E∗,E¯∗,k∗,k¯∗,fric,pb_ten,pb_cohjT⟶KN−j=σc−j∗,εc−j∗,Ej,σn−j.(2) Determination of the Target Inversion Value from the Uniaxial Compression Data. The experimental uniaxial compression tests provided the real stress-strain behaviour and AE monitoring data. According to the stress-strain relationship and the definition of the cracking stress, the peak stress σc, peak strain εc, elastic modulus (E), and crack initiation stress (σ1−p) of the MCS material were obtained and used as the target values for inversion. It should be noted that AE monitors low-energy high-frequency seismic events resulting from crack initiation, so it can be used for the identification of σ1−p [30]:(3)F=σc,εc,E,σ1−p.(3) Construction of the Fitness Function of the Inversion Algorithm. According to the combination of different mesoparameters Mjj∈1,Ni, the study focused on the relatively small errors between the simulated and experimental characteristic parameters KN−jj∈1,Ni, where F is the target value. The objective function of the inversion of the mesoparameters can be defined as(4)Gk=KN−j−F.Furthermore, the fitness function was constructed as follows:(5)fk=1Gk.It was found that the larger individual fitness values resulted in smaller objective function values; the mesoparameters reflected the real MCS material properties.


### 3.2.2. Calibration Results and Comparison of Calibration Methods

According to the above calibration method, the mesoparameters with small errors were observed, as shown in Table2. The uniaxial stress-strain curves and failure characteristics calculated from the numerical simulation are compared to that of the experimental data in Figure 10.Table 2
Mesoparameters of material.
Contact properties
E
∗ (N/m2)
E
¯
∗ (N/m2)
k
∗
k
¯
∗
fric
pb_ten (MPa)
pb_coh (MPa)
λ
Sand-sand
1.50 × 109
1.50 × 109
1.2
1.2
0.7
13.05 ± 1.0
15.05 ± 1.5
0.25∼0.50
Sand-calcite
0.85 × 109
0.85 × 109
2.5
2.5
0.5
17.05 ± 1.0
24.05 ± 1.5
0.50
Calcite-calcite
1.25 × 109
1.25 × 109
2.5
2.5
0.5
14.05 ± 1.0
17.05 ± 1.5
0.50Figure 10
Comparison of numerical simulation and experimental stress-strain curves and failure modes.The results shown in Figure10 were found by comparing the numerical results to the results of physical experiments, with 1.95% and 4.0% as the relative errors in the peak stress and peak strain, respectively. Additionally, the numerical tests revealed a similar trend between the macroscale cracking and the overall fracturing feature with minor local deviations. The deviations are primarily due to the discrete connection of the microscale cracks in the macroscale cracks, which are inherently induced by the discrete contacts between particles in the model developed in this study. Meanwhile, some simulation results differ from the testing results; specifically, the elastic module obtained in each case differs from the others. This difference is justified by the model reconstruction of microscale particle structure and the related physical properties. In a numerical simulation, it is almost impossible to reconstruct a model sharing the same mesoscale structure feature and physical properties as those of a physical testing sample [31]. In view of these factors, the mesoscale parameters used in this study are appropriate for characterising the MCS material properties.Currently, the trial-and-error method has been the most commonly used calibration method for selecting appropriate mesoscale parameters [32, 33]. The stress-strain curve corresponding to the mesoparameters obtained using the trial-and-error method was plotted on the same graph (Figure 10) in order to compare the results of ADE and the trial-and-error method. Figure 10 clearly illustrates that the numerical curve obtained by the ADE method is similar to the curve obtained by Tang et al. [13] using the trial-and-error method. By comparison with the mesoparameters obtained by Tang et al. [13], most mesoparameters are almost the same, except for pb_ten and pb_coh.The mesoparameters used for the trial-and-error method usually have to be assigned according to the mesoscopic properties of the material. For relatively simple contact types, the trial-and-error time may be relatively short. However, when the material involves multiple contact types, the trial-and-error process is extremely complicated and takes more time. Compared to the trial-and-error method, ADE is able to create an initial database and obtains target parameter values through parameter inversion. The entire process only takes a certain amount of time during the initial database setup process. Once the initial database is established, corresponding mesoscopic parameter values can be quickly obtained. Furthermore, the established database can provide initial data for the subsequent acquisition of mesoparameters of related materials with similar mechanical behaviour.


## 3.2.1. Calibration Method

In this paper, the mesoparameters of MCS material were calibrated using the ADE method. The calibrated mesoparameters includeE∗, E¯∗, k∗, k¯∗, fric, pb_ten, pb_coh, and R¯.The SEM characterisation of the material (shown in Figure5) makes it necessary to consider the different bonding characteristics and discreteness of MCS material when setting mesoparameters. According to the results of Feng et al. [22], different bonding characteristics can be characterized by different bonding radii. When the parallel bonding radius, R¯ (m), is determined, the particle flow theory suggested that the mesoscale bonding radius, R¯, was related to the radial coefficient, λ. For the sand particle-calcite and calcite-calcite bonds, λ was assigned a fixed value of 0.5 because these bonds were relatively strong according to the SEM analysis results. For the sand particle-sand particle bond, a random number to vary the parallel-bond radial coefficient λ between 0.25 and 0.5 was introduced based on the previous PFC study on microbial cemented sand column by Feng et al. [22]. Thus, the bonding radius was determined and is not included in the parameter inversion procedure.The mesoparametersM of the MCS material used in the simulation can be expressed as(1)M=E∗,E¯∗,k∗,k¯∗,fric,pb_ten,pb_cohT,where E∗ is the effective modulus (force/area), E¯∗ is the effective modulus of the bond (stress), k∗ is the normal-to-shear stiffness ratio, k¯∗ is the normal-to-shear stiffness ratio of the bond, fric is the friction coefficient, pb_ten is the tensile strength (stress), and pb_coh is the cohesion (stress).The adaptive differential evolution algorithm used to invert the mesoparameters is shown in Figure9, and the overall process is described in the following three steps.Figure 9
The realisation process of the inversion algorithm of mesoparameters.(1) Construction of the Characteristic Parameters. Using the PFC simulating uniaxial compression test, Nii≥20 groups of mesoparameters Mjj∈1,Ni were constructed, and the stress-strain and crack data corresponding to these parameters were obtained. It was noted from the laboratory test results of elastic modulus size of the MCS material that the elastic modulus of sand particles is relatively large and that of calcite is relatively small [29]. Thus, the initial elastic modulus parameters were set according to this size. Subsequently, based on the stress-strain and crack data, the characteristic parameters were calculated using MATLAB (R2015b) software. The obtained characteristic parameters were used as the initial population samples NPjj∈1,Ni for the adaptive differential evolution algorithm. The characteristic parameters KN−jj∈1,Ni included the peak stress (σc−j∗), peak strain (εc−j∗), elastic modulus (Ej), and the crack initiation stress (σn−j). Note that the crack initiation stress (σn−j) was identified as the axial stress where the number of cracks reached 1% of the total number of cracks at peak stress [20]:(2)Mj=E∗,E¯∗,k∗,k¯∗,fric,pb_ten,pb_cohjT⟶KN−j=σc−j∗,εc−j∗,Ej,σn−j.(2) Determination of the Target Inversion Value from the Uniaxial Compression Data. The experimental uniaxial compression tests provided the real stress-strain behaviour and AE monitoring data. According to the stress-strain relationship and the definition of the cracking stress, the peak stress σc, peak strain εc, elastic modulus (E), and crack initiation stress (σ1−p) of the MCS material were obtained and used as the target values for inversion. It should be noted that AE monitors low-energy high-frequency seismic events resulting from crack initiation, so it can be used for the identification of σ1−p [30]:(3)F=σc,εc,E,σ1−p.(3) Construction of the Fitness Function of the Inversion Algorithm. According to the combination of different mesoparameters Mjj∈1,Ni, the study focused on the relatively small errors between the simulated and experimental characteristic parameters KN−jj∈1,Ni, where F is the target value. The objective function of the inversion of the mesoparameters can be defined as(4)Gk=KN−j−F.Furthermore, the fitness function was constructed as follows:(5)fk=1Gk.It was found that the larger individual fitness values resulted in smaller objective function values; the mesoparameters reflected the real MCS material properties.


## 3.2.2. Calibration Results and Comparison of Calibration Methods

According to the above calibration method, the mesoparameters with small errors were observed, as shown in Table2. The uniaxial stress-strain curves and failure characteristics calculated from the numerical simulation are compared to that of the experimental data in Figure 10.Table 2
Mesoparameters of material.
Contact properties
E
∗ (N/m2)
E
¯
∗ (N/m2)
k
∗
k
¯
∗
fric
pb_ten (MPa)
pb_coh (MPa)
λ
Sand-sand
1.50 × 109
1.50 × 109
1.2
1.2
0.7
13.05 ± 1.0
15.05 ± 1.5
0.25∼0.50
Sand-calcite
0.85 × 109
0.85 × 109
2.5
2.5
0.5
17.05 ± 1.0
24.05 ± 1.5
0.50
Calcite-calcite
1.25 × 109
1.25 × 109
2.5
2.5
0.5
14.05 ± 1.0
17.05 ± 1.5
0.50Figure 10
Comparison of numerical simulation and experimental stress-strain curves and failure modes.The results shown in Figure10 were found by comparing the numerical results to the results of physical experiments, with 1.95% and 4.0% as the relative errors in the peak stress and peak strain, respectively. Additionally, the numerical tests revealed a similar trend between the macroscale cracking and the overall fracturing feature with minor local deviations. The deviations are primarily due to the discrete connection of the microscale cracks in the macroscale cracks, which are inherently induced by the discrete contacts between particles in the model developed in this study. Meanwhile, some simulation results differ from the testing results; specifically, the elastic module obtained in each case differs from the others. This difference is justified by the model reconstruction of microscale particle structure and the related physical properties. In a numerical simulation, it is almost impossible to reconstruct a model sharing the same mesoscale structure feature and physical properties as those of a physical testing sample [31]. In view of these factors, the mesoscale parameters used in this study are appropriate for characterising the MCS material properties.Currently, the trial-and-error method has been the most commonly used calibration method for selecting appropriate mesoscale parameters [32, 33]. The stress-strain curve corresponding to the mesoparameters obtained using the trial-and-error method was plotted on the same graph (Figure 10) in order to compare the results of ADE and the trial-and-error method. Figure 10 clearly illustrates that the numerical curve obtained by the ADE method is similar to the curve obtained by Tang et al. [13] using the trial-and-error method. By comparison with the mesoparameters obtained by Tang et al. [13], most mesoparameters are almost the same, except for pb_ten and pb_coh.The mesoparameters used for the trial-and-error method usually have to be assigned according to the mesoscopic properties of the material. For relatively simple contact types, the trial-and-error time may be relatively short. However, when the material involves multiple contact types, the trial-and-error process is extremely complicated and takes more time. Compared to the trial-and-error method, ADE is able to create an initial database and obtains target parameter values through parameter inversion. The entire process only takes a certain amount of time during the initial database setup process. Once the initial database is established, corresponding mesoscopic parameter values can be quickly obtained. Furthermore, the established database can provide initial data for the subsequent acquisition of mesoparameters of related materials with similar mechanical behaviour.


## 4. Coupled Thermomechanical Model

In the following section, the thermal-mechanical model is constructed.


### 4.1. Thermal Model

The heat conduction in a continuum is given by(6)−∂qi∂xi+qv=ρC∂T∂t,where qi is the heat flux vector (W/m2), qv is the volumetric heat intensity (W/m3), ρ is the density of the material (kg/m3), C is the heat capacity under constant pressure (J/kg·K), and T is the temperature (K).The heat outflow per unit volume is given by the divergence ofqi. The average value of the divergence of qi in the reservoir is given by(7)∂qi∂xi=1V∭V∂qi∂xidV.The volume integral can be replaced with a surface integral by the application of the Gauss divergence theorem to the reservoir. In addition, if the heat flow is restricted only to theN thermal pipes associated with the reservoir, the surface integral can be replaced by the following summation:(8)∭V∂qi∂xidV=∯SqieidS=∑P=1NqiPeiPΔSP=∑P=1NQP.Hence,(9)∂qi∂xi=1V∑P=1NQP.If each pipep is regarded as a one-dimensional object with a thermal resistance per unit length of η, then the power in a pipe would be given by(10)Qp=−ΔTηlp,where ΔT is the temperature difference between the two reservoirs on each end of the pipe and lp is the pipe length.A heat source,Qp, and thermal resistance, η, are associated with each pipe, which has a finite length equal to the distance between the centroids of the two particles that it connects. For a two-dimensional model, the thermal resistance, η, is given by(11)η=12K1−n∑NbVb∑Nlp,where n is the porosity of the numerical specimen, Vb is the volume of the ball (b), lp is the length of the thermal pipe p, Nb is the number of balls, K is the thermal conductivity, and N is the number of thermal pipes.For a two-dimensional system, the ball disk in the PFC is unit 1. Hence,Vb=Sb. From equations (6)–(11), the governing equation can be rewritten as follows:(12)∑P=1N2K∑NbSbΔT1−n∑Nlplp+Qv=mC∂T∂t.The amount of thermal energy deposited into a material due to microwave heating, which is equated to the volumetric heat intensity, is given by [7](13)Pd=2πfε0εr″E02.where Pd is the power density (W/m3), f is the frequency of the microwave irradiation (Hz), ε0 is the permittivity of free space (8.854 × 10−12 F/m), εr″ is the dielectric loss factor of the mineral, and E0 is the electric field intensity resulting from the microwave power (V/m).SubstitutingPd in place of Qv in equation (12) gives(14)∑P=1N2K∑NbSbΔT1−n∑NlPlp+Pd=mC∂T∂t.Thus, the temperature increase can be obtained by solving equation (14) for a given power density and heating duration.


### 4.2. Thermal-Mechanical Coupling

The temperature field can be simulated using the governing equations for microwave heating. Based on the definition of different material properties, the conversion of microwave energy was achieved by setting the power density and heating duration in PFC modelling.When the microwave energy was converted into heat energy, thermal strains were produced in the material, resulting in changes in the particle radii and the force transmitted by each parallel bond [34]. This should be taken into consideration in the heating of both the particles and the bonding material. For a temperature change ΔT, the change in the particle radius R is given by(15)ΔR=αRΔT,where α is the coefficient of thermal expansion of the particle.If a parallel bond is present at the contact associated with a pipe, then(16)ΔF¯n=−K¯nAα¯L¯ΔT,where K¯n is the normal bond stiffness, A is the cross-sectional area of the bond, α¯ is the expansion coefficient, L¯ is the bond length (assumed to be equal to the distance between the centroids of the two particles connected by the pipe), and ΔT is the temperature change.


### 4.3. Boundary Conditions and Thermal Parameters

The definition of the boundary conditions generally requires a description of the local boundary in terms of an outward normal unit vector. More realistic boundary conditions also consider the convection and radiative boundary properties between the specimen and the surrounding air [35]. However, a material modelled in the PFC 2D software has no explicit representation of its boundary (unlike in a continuum code, where zone surfaces comprise the boundary). This problem is overcome by specifying the power input to a boundary surface instead of the flux. The external boundaries of the thermal model were considered to be thermally insulated. Therefore, the convection and radiation were not incorporated into the model since the microwave heating durations were quite short.Considering that the thermal input properties of an MCS was determined by the microstructural components, two different sets of thermal parameters were obtained from the literature [36–38], as presented in Table 3. Based on a previous study [33], the microwave heating durations were between 0.02 s and 0.50 s, and the power density was equal to 1 × 109 W/m3 in this study.Table 3
Thermal input data for sand and calcite.
Property
Quartz
Calcite
Specific heat (J/kg·K)
700
817.7
Thermal conductivity (W/m·K)
7.69
3.16
Thermal expansion coefficient (1/K)
8.1 × 10−6
13.1 × 10−6
Density (kg/m3)
2648
2712


## 4.1. Thermal Model

The heat conduction in a continuum is given by(6)−∂qi∂xi+qv=ρC∂T∂t,where qi is the heat flux vector (W/m2), qv is the volumetric heat intensity (W/m3), ρ is the density of the material (kg/m3), C is the heat capacity under constant pressure (J/kg·K), and T is the temperature (K).The heat outflow per unit volume is given by the divergence ofqi. The average value of the divergence of qi in the reservoir is given by(7)∂qi∂xi=1V∭V∂qi∂xidV.The volume integral can be replaced with a surface integral by the application of the Gauss divergence theorem to the reservoir. In addition, if the heat flow is restricted only to theN thermal pipes associated with the reservoir, the surface integral can be replaced by the following summation:(8)∭V∂qi∂xidV=∯SqieidS=∑P=1NqiPeiPΔSP=∑P=1NQP.Hence,(9)∂qi∂xi=1V∑P=1NQP.If each pipep is regarded as a one-dimensional object with a thermal resistance per unit length of η, then the power in a pipe would be given by(10)Qp=−ΔTηlp,where ΔT is the temperature difference between the two reservoirs on each end of the pipe and lp is the pipe length.A heat source,Qp, and thermal resistance, η, are associated with each pipe, which has a finite length equal to the distance between the centroids of the two particles that it connects. For a two-dimensional model, the thermal resistance, η, is given by(11)η=12K1−n∑NbVb∑Nlp,where n is the porosity of the numerical specimen, Vb is the volume of the ball (b), lp is the length of the thermal pipe p, Nb is the number of balls, K is the thermal conductivity, and N is the number of thermal pipes.For a two-dimensional system, the ball disk in the PFC is unit 1. Hence,Vb=Sb. From equations (6)–(11), the governing equation can be rewritten as follows:(12)∑P=1N2K∑NbSbΔT1−n∑Nlplp+Qv=mC∂T∂t.The amount of thermal energy deposited into a material due to microwave heating, which is equated to the volumetric heat intensity, is given by [7](13)Pd=2πfε0εr″E02.where Pd is the power density (W/m3), f is the frequency of the microwave irradiation (Hz), ε0 is the permittivity of free space (8.854 × 10−12 F/m), εr″ is the dielectric loss factor of the mineral, and E0 is the electric field intensity resulting from the microwave power (V/m).SubstitutingPd in place of Qv in equation (12) gives(14)∑P=1N2K∑NbSbΔT1−n∑NlPlp+Pd=mC∂T∂t.Thus, the temperature increase can be obtained by solving equation (14) for a given power density and heating duration.


## 4.2. Thermal-Mechanical Coupling

The temperature field can be simulated using the governing equations for microwave heating. Based on the definition of different material properties, the conversion of microwave energy was achieved by setting the power density and heating duration in PFC modelling.When the microwave energy was converted into heat energy, thermal strains were produced in the material, resulting in changes in the particle radii and the force transmitted by each parallel bond [34]. This should be taken into consideration in the heating of both the particles and the bonding material. For a temperature change ΔT, the change in the particle radius R is given by(15)ΔR=αRΔT,where α is the coefficient of thermal expansion of the particle.If a parallel bond is present at the contact associated with a pipe, then(16)ΔF¯n=−K¯nAα¯L¯ΔT,where K¯n is the normal bond stiffness, A is the cross-sectional area of the bond, α¯ is the expansion coefficient, L¯ is the bond length (assumed to be equal to the distance between the centroids of the two particles connected by the pipe), and ΔT is the temperature change.


## 4.3. Boundary Conditions and Thermal Parameters

The definition of the boundary conditions generally requires a description of the local boundary in terms of an outward normal unit vector. More realistic boundary conditions also consider the convection and radiative boundary properties between the specimen and the surrounding air [35]. However, a material modelled in the PFC 2D software has no explicit representation of its boundary (unlike in a continuum code, where zone surfaces comprise the boundary). This problem is overcome by specifying the power input to a boundary surface instead of the flux. The external boundaries of the thermal model were considered to be thermally insulated. Therefore, the convection and radiation were not incorporated into the model since the microwave heating durations were quite short.Considering that the thermal input properties of an MCS was determined by the microstructural components, two different sets of thermal parameters were obtained from the literature [36–38], as presented in Table 3. Based on a previous study [33], the microwave heating durations were between 0.02 s and 0.50 s, and the power density was equal to 1 × 109 W/m3 in this study.Table 3
Thermal input data for sand and calcite.
Property
Quartz
Calcite
Specific heat (J/kg·K)
700
817.7
Thermal conductivity (W/m·K)
7.69
3.16
Thermal expansion coefficient (1/K)
8.1 × 10−6
13.1 × 10−6
Density (kg/m3)
2648
2712


## 5. Heating Effect Analysis Using Coupled Thermomechanical Model

### 5.1. Heating Effects

The temperature fields as a function of the MCS material under different heating durations are shown in Figure11, where a gradual increase in the sample temperature as a result of longer simulation time of heating can be observed. The distributions of different particles at the computational time of 0.50 s in the testing sample were illustrated with sand in crimson and calcite in light red. It should be noted that the heating-induced cracks in the sample have not been accounted for in this study. The focus is on the influence of short-term heating on the deformed samples under different magnitudes of strains.Figure 11
The temperature field under the different heating durations.


### 5.2. Influence of Strain Magnitudes on the Failure Characteristics

The results of the simulation using the thermomechanical model were interpolated using a python language compiler with the uncoupled model being included for a direct comparison (Figure12). This article defines two types of cracks: cracks that occur only under mechanical loading are called mechanical cracks (shown in blue), and cracks that occur under temperature loading are called thermal cracks. It should be noted that the same heating duration of 0.50 s (corresponding to T = 210°C) is adopted under five different strain magnitudes.Figure 12
The failure characteristics under the different strains. The failure distribution of MCS material (a) without considering temperature load and (b) considering temperature load 0.50 s.
(a)
(b)As shown from Figure12(a), the results captured from uncoupled modelling can be summarized as follows: (1) when the strain was 1.5‰, few discrete microcracks were present throughout the sample; (2) when the strain exceeded 1.5‰, there were mechanical cracks in the structure (inside there exist damage) and mainly in the quartz matrix; (3) under larger strain (3.0‰), more microcracks gradually commonly coalesced, with few cracks in the calcite.As shown from Figure12(b), thermal influence on the cracks can be observed; when the strain exceeded 1.5‰, the number of microcracks (shown in cyan) caused by thermal cracking continuously increased. Specifically, it was found that a large number of microthermal cracks (shown in the green frame) formed and coalesced in the mechanical crack tip (shown in blue) and its surroundings when the strain was 3.0‰. These indicated that the mechanical crack (shown in blue) acted as a “rock bridge” with the increasing strain, where the damage to the material played a guiding role in the thermal crack (shown in cyan).Besides, as shown from Figure12(b), the thermal crack in the material mainly extended to the quartz matrix. Although the hardness of quartz is clearly higher than that of calcite, the bonding between the quartz particles in the MCS material is determined by the cementation. Compared to the calcite particle (Figure 5), weak cementation occurred in the quartz sample (some holes were observed at the bond interface that decreased the cementation strength and degraded the cementation performance), resulting in uneven bond strengths. Therefore, the bond rupture was more likely to occur, and the crack initially propagated through the quartz matrix.According to the simulation results shown in Figure12(b), we can also see that the crack surface was rough during crack propagation and that the crack bifurcation was easily observed when the strain was 3.0‰. For the crack bifurcation phenomenon, it showed that, under the action of stress, the crack propagation path with the minimum impedance is taken [39]. This demonstrates the heterogeneity of the material and confirms that the numerical model is reasonable for considering the strength distribution of the actual sample.


### 5.3. Influence of Increase in Temperature on the Failure Characteristics

Figure13 shows that thermal cracks (shown in cyan) for a sample, with a strain of 3.0‰, were observed as the temperature increased. During the heating process, a tensile stress caused by microwave energy is formed in the material [40–42].Figure 13
The failure characteristics under the different heating durations (3.0‰).At a relatively lower temperature, the tensile stress generated is smaller. The fracture effect for this thermal crack developed a quiet period with an increase in time for the duration of 0.25 s (corresponding toT = 110°C) of the quiet period. With the increase in temperature, the tensile stress rapidly increased after the experience of higher temperature, and thermal cracks gradually converged and interacted, resulting in the appearance of obvious failure zones (Figure 13: 0.50 s (corresponding to T = 210°C)).These results were obtained because the increase in heating duration gradually raised the temperature of the material. The temperature in the MCS material was not very high for a shorter heating period, and the tensile stress generated potentially promoted the fracture at a location of the model with lower internal bond strength. A location with higher bond strength was not sufficient to induce bond fracture. Therefore, in the short term, the fracture of the model was not very obvious. Subsequently, the temperature rose rapidly as the heating duration increased, and the generated tensile stress exceeded the bond strength between the particles. The new thermal cracks developed rapidly, and many small thermal cracks gradually converged and interacted, resulting in the appearance of obvious failure zones (shown in the cyan lines in the final sample). Finally, the thermal cracks exhibited the behavioural characteristics across the structural plane in the transverse direction.In summary, the temperature gradually changed with an increase in heating duration during the heating process, and the material inevitably produced a rupture, which further led to a change in the macromechanical behaviour of the MCS material. The application of temperature has a superposition effect on the damage of the MCS material.


### 5.4. Comparison between Calculated and Experiment Results

As mentioned above, for the MCS material, the cracks were preferentially located in the internal damage area of the structure under the action of temperature. With the increase of temperature, many small thermal cracks gradually converged and interacted with the obvious failure zone. To validate the accuracy of numerical simulation results, a heating test was carried out. The heating test has been simplified to verify the resulting effect of temperature on the existing damage zone.The experimental specimens were selected after uniaxial compression, and the heating experiment was conducted using heating equipment (KW-202). The longest heating duration, 285 min, was chosen to allow the temperature to reach at approximately 210°C, which corresponded to the maximum value in the numerical analysis. Subsequently, the damage images of the MCS material under the different heating durations via 3D Super Depth Digital Microscope (VHX-5000) were observed. Figure14 shows the heating test process with the test results shown in Figure 15.Figure 14
The heating test process of MCS material.Figure 15
The damage images of the MCS material under the different heating durations. (a) 0 min, (b) 15 min, (c) 45 min, (d) 105 min, (e) 195 min, and (f) 285 min.
(a)
(b)
(c)
(d)
(e)
(f)From Figure15, it can be observed that when the heating duration was 15 min, the particles in the damaged area were closely bonded, with a few local cracks (Figure 15(b)). As the heating durations increased, more localized cracks were observed, and the interaction of the localized cracks occurred (Figure 15(e)). When the heating duration reached 285 min, the crack had spread from the damage boundary to the surrounding area (as shown in Figure 15(f)). According to the test results, the failure mode of the experimental results is consistent with the simulation results, and it confirms that the simulation results are reasonable.


## 5.1. Heating Effects

The temperature fields as a function of the MCS material under different heating durations are shown in Figure11, where a gradual increase in the sample temperature as a result of longer simulation time of heating can be observed. The distributions of different particles at the computational time of 0.50 s in the testing sample were illustrated with sand in crimson and calcite in light red. It should be noted that the heating-induced cracks in the sample have not been accounted for in this study. The focus is on the influence of short-term heating on the deformed samples under different magnitudes of strains.Figure 11
The temperature field under the different heating durations.


## 5.2. Influence of Strain Magnitudes on the Failure Characteristics

The results of the simulation using the thermomechanical model were interpolated using a python language compiler with the uncoupled model being included for a direct comparison (Figure12). This article defines two types of cracks: cracks that occur only under mechanical loading are called mechanical cracks (shown in blue), and cracks that occur under temperature loading are called thermal cracks. It should be noted that the same heating duration of 0.50 s (corresponding to T = 210°C) is adopted under five different strain magnitudes.Figure 12
The failure characteristics under the different strains. The failure distribution of MCS material (a) without considering temperature load and (b) considering temperature load 0.50 s.
(a)
(b)As shown from Figure12(a), the results captured from uncoupled modelling can be summarized as follows: (1) when the strain was 1.5‰, few discrete microcracks were present throughout the sample; (2) when the strain exceeded 1.5‰, there were mechanical cracks in the structure (inside there exist damage) and mainly in the quartz matrix; (3) under larger strain (3.0‰), more microcracks gradually commonly coalesced, with few cracks in the calcite.As shown from Figure12(b), thermal influence on the cracks can be observed; when the strain exceeded 1.5‰, the number of microcracks (shown in cyan) caused by thermal cracking continuously increased. Specifically, it was found that a large number of microthermal cracks (shown in the green frame) formed and coalesced in the mechanical crack tip (shown in blue) and its surroundings when the strain was 3.0‰. These indicated that the mechanical crack (shown in blue) acted as a “rock bridge” with the increasing strain, where the damage to the material played a guiding role in the thermal crack (shown in cyan).Besides, as shown from Figure12(b), the thermal crack in the material mainly extended to the quartz matrix. Although the hardness of quartz is clearly higher than that of calcite, the bonding between the quartz particles in the MCS material is determined by the cementation. Compared to the calcite particle (Figure 5), weak cementation occurred in the quartz sample (some holes were observed at the bond interface that decreased the cementation strength and degraded the cementation performance), resulting in uneven bond strengths. Therefore, the bond rupture was more likely to occur, and the crack initially propagated through the quartz matrix.According to the simulation results shown in Figure12(b), we can also see that the crack surface was rough during crack propagation and that the crack bifurcation was easily observed when the strain was 3.0‰. For the crack bifurcation phenomenon, it showed that, under the action of stress, the crack propagation path with the minimum impedance is taken [39]. This demonstrates the heterogeneity of the material and confirms that the numerical model is reasonable for considering the strength distribution of the actual sample.


## 5.3. Influence of Increase in Temperature on the Failure Characteristics

Figure13 shows that thermal cracks (shown in cyan) for a sample, with a strain of 3.0‰, were observed as the temperature increased. During the heating process, a tensile stress caused by microwave energy is formed in the material [40–42].Figure 13
The failure characteristics under the different heating durations (3.0‰).At a relatively lower temperature, the tensile stress generated is smaller. The fracture effect for this thermal crack developed a quiet period with an increase in time for the duration of 0.25 s (corresponding toT = 110°C) of the quiet period. With the increase in temperature, the tensile stress rapidly increased after the experience of higher temperature, and thermal cracks gradually converged and interacted, resulting in the appearance of obvious failure zones (Figure 13: 0.50 s (corresponding to T = 210°C)).These results were obtained because the increase in heating duration gradually raised the temperature of the material. The temperature in the MCS material was not very high for a shorter heating period, and the tensile stress generated potentially promoted the fracture at a location of the model with lower internal bond strength. A location with higher bond strength was not sufficient to induce bond fracture. Therefore, in the short term, the fracture of the model was not very obvious. Subsequently, the temperature rose rapidly as the heating duration increased, and the generated tensile stress exceeded the bond strength between the particles. The new thermal cracks developed rapidly, and many small thermal cracks gradually converged and interacted, resulting in the appearance of obvious failure zones (shown in the cyan lines in the final sample). Finally, the thermal cracks exhibited the behavioural characteristics across the structural plane in the transverse direction.In summary, the temperature gradually changed with an increase in heating duration during the heating process, and the material inevitably produced a rupture, which further led to a change in the macromechanical behaviour of the MCS material. The application of temperature has a superposition effect on the damage of the MCS material.


## 5.4. Comparison between Calculated and Experiment Results

As mentioned above, for the MCS material, the cracks were preferentially located in the internal damage area of the structure under the action of temperature. With the increase of temperature, many small thermal cracks gradually converged and interacted with the obvious failure zone. To validate the accuracy of numerical simulation results, a heating test was carried out. The heating test has been simplified to verify the resulting effect of temperature on the existing damage zone.The experimental specimens were selected after uniaxial compression, and the heating experiment was conducted using heating equipment (KW-202). The longest heating duration, 285 min, was chosen to allow the temperature to reach at approximately 210°C, which corresponded to the maximum value in the numerical analysis. Subsequently, the damage images of the MCS material under the different heating durations via 3D Super Depth Digital Microscope (VHX-5000) were observed. Figure14 shows the heating test process with the test results shown in Figure 15.Figure 14
The heating test process of MCS material.Figure 15
The damage images of the MCS material under the different heating durations. (a) 0 min, (b) 15 min, (c) 45 min, (d) 105 min, (e) 195 min, and (f) 285 min.
(a)
(b)
(c)
(d)
(e)
(f)From Figure15, it can be observed that when the heating duration was 15 min, the particles in the damaged area were closely bonded, with a few local cracks (Figure 15(b)). As the heating durations increased, more localized cracks were observed, and the interaction of the localized cracks occurred (Figure 15(e)). When the heating duration reached 285 min, the crack had spread from the damage boundary to the surrounding area (as shown in Figure 15(f)). According to the test results, the failure mode of the experimental results is consistent with the simulation results, and it confirms that the simulation results are reasonable.


## 6. Conclusions

The material physical characterizations at different scales were analysed in this paper. A method for calibrating the input MCS material properties was designed, and a coupled thermomechanical model was constructed. The behaviour of the MCS material subjected to temperature loading was investigated under strains varying from 1.0‰ to 3.0‰. The following conclusions were drawn from the study:(1)
From a comparative analysis between ADE and trial-and-error method, the adaptive differential evolution algorithm was successfully used to explore the MCS mesoparameters and can also be used to establish the comprehensive database of parameters for the MCS material in DEM analysis.(2)
For the wide range of the applied strain magnitudes (lower than that corresponds to the failure state), the elastic deformation can be observed under low-strain magnitudes (<1.5‰), while the progressive development of the localized cracks is more evident under high-strain magnitudes. These cracks occur more frequently at the sand-sand interface due to lower bonding strength among them.(3)
For the typical strain of 3.0‰ associated with the prefailure state under the different temperatures, thermal cracks behaviour experience two phases. (1) A quiet period: local crack progressively developed and the thermal effect on the crack might be ignored. (2) An explosion period: the local cracks gradually converged and interacted within the distinct failure zone. The length of the thermal cracks in the transverse direction increased throughout the heating process.



---

*Source: 7528746-2018-12-13.xml*