Characterization of 3D Displacement and Stress Fields in Coal Based on CT Scans

: Computed tomography (CT) scans were performed on samples of an outburst-prone coal seam at different loading stages. The area and roundness of the CT images were used to quantify the degree of the coal macroscopic deformation under different loads. A spatial matching algorithm was used to calculate the three-dimensional (3D) displacement ﬁelds of different regions of interest (ROIs, containing primary fractures, minerals, and only coal) under different loads. The presence of fractures and minerals were found to promote and inhibit displacement, respectively, and the 3D displacement ﬁeld data followed a normal distribution. A meshfree numerical simulation was used to determine the 3D maximum principal stress, shear stress and displacement ﬁelds under different loads. The following results were obtained: fractures and minerals signiﬁcantly affect the stress state and displacement ﬁeld distribution features, the maximum principal stresses and shear stresses in different matrices differ signiﬁcantly, and the presence of minerals and fractures induce a prevalent shear stress in coal and make coal prone to stress concentration.


Introduction
Coal is typically anisotropic, and high sample dispersion makes it challenging to determine the deformation features and stress field patterns of raw coal samples. Considering the importance of the stress and deformation fields in analysing coal stability, a computed tomography (CT) scan with a loading function was used to determine the deformation features of coal under different loads in this study. A novel numerical simulation method was applied to a three-dimensional (3D) CT model to analyse the coal stress state.
A variety of experimental and theoretical methods have been used to investigate the internal structure and deformation features of coal over the past few years [1][2][3][4][5][6][7], which has advanced the prevention and control of related hazards. CT scans enable the noncontact nondestructive testing of samples and can be combined with 3D reconstruction techniques to quantitatively analyse the coal mesostructure. Kawakata et al. [8] first used CT to scan the damage state of rock after triaxial loading but the samples were not loaded during the scanning process. Cao et al. [9] analysed the fracture evolution and seismic response features of hydraulic fracturing under triaxial loading conditions and performed CT scans to quantitatively characterise pre-existing and hydraulic fractures. Fan et al. [10] performed threshold segmentation on CT scan data and statistically analysed the pore size distribution, throat radius, throat length, and coordination number of the pore network model. Li et al. [11] experimentally analysed the pore features of coal samples with different bursting proneness by low-temperature nitrogen adsorption and desorption, scanning electron microscopy (SEM) and CT. The results showed that stress exacerbates the complexity of the pore structure. Li et al. [12] used CT to visualize the evolution of the fracture features of coal samples under uniaxial and triaxial loading conditions. Quantitative characterisation of the fractures demonstrated that the grey value, fracture volume and CT porosity can be used to quantitatively evaluate the degree of damage to coal. Lu et al. [13] found using CT scans on coal samples under different load conditions that the expansion of primary fractures and the development of new fractures increase the coal connectivity. Liu et al. [14] analysed the porosity of coal through CT scans and NMR experiments and quantitatively characterised the fractal dimension. Tao et al. [15] analysed the porosity and fracture heterogeneity of coal using a variety of experimental methods (SEM, carbon dioxide and nitrogen adsorption, and CT scans) and proposed an index to quantitatively predict coal macrolithotypes. Wang et al. [16][17][18] used CT scans to analyse the pore structure and seepage features of coal, which were then used to conduct a dynamic water injection simulation analysis. Wu et al. analysed the fractal features of a coal fracture network to find that the porosity is exponentially related to fractal dimension. Zhang et al. [19] used various experimental methods to perform a multiscale characterisation of coal samples of different ranks. The results of the aforementioned studies have revealed the variation features of fractures under stress and the corresponding influence on gas migration.
Wang et al. [20] analysed two-dimensional (2D) CT scan images of coal samples as well as 3D reconstruction models to find that the confining pressure and presence of minerals have significant effects on the features of fracture propagation. Mao et al. [21] compared digital volumetric speckle photographs with the deformation results obtained from CT images to determine the deformation features of coal at different loading stages. Sampath et al. [22] generated a 2D fracture network from CT imaging to study the gas flow and coal deformation process. Roslin et al. [23] used SEM and CT to analyse pore structure and carried out numerical simulation based on the reconstruction model.
The 3D deformation features of coal during different loading processes are investigated in this study. A feature point spatial matching algorithm is used to calculate and analyse the displacement field. A meshfree numerical simulation of the 3D reconstruction model is performed to analyse the 3D stress field and the stress states of different matrices (coal, minerals, and regions near fractures).

Sample Preparation
The studied coal samples were taken from the no. 3 coal seam of the Sihe coal mine in Shanxi Province, China. The coal seam is located at the southern end of the Qinshui coalfield and presents an outburst risk. The coal samples were cut into ϕ25 mm × h50 mm cylinders, of which the end faces were polished to ensure parallelism.

Experimental Device and Methods
The experimental device consists of a microfocus CT scan system and a triaxial stressdisplacement test system. A Vtomex L300 CT scan system with an observation resolution of up to 500 nm was used. The triaxial stress-displacement test system can perform uniaxial and triaxial loading tests on ϕ25 mm × h50 mm samples, with a maximum axial compression of 100 kN, a maximum confining pressure of 30 MPa, and a loading rate range of 0.1-3 mm/min. The attenuation coefficient of the X-ray intensity and the density of the material satisfy the following equation [24]:

of 16
where I o is the initial intensity of the X-ray; I is the intensity of the X-ray after passing through the material; µ is the attenuation coefficient of the material; and x is the penetration distance of the ray. Three sets of phased CT scan tests were performed on the Sihe coal samples under the loading process, with an axial loading rate of 0.2 mm/min and a CT scan resolution of 30 µm. The axial stress-strain data were recorded, and the axial compression was maintained constant during the CT scan. Figure 1 presents CT images showing highly discernible minerals inside the coal. Using the minerals as feature points, the Z-axis position in the CT images was recorded to delineate the upper and lower boundaries of the 3D model. Threshold segmentation was performed on the CT data to determine the boundaries of the 3D model, thus accurately delineating the range of the 3D CT model at each stage.

Boundary Delineation of 3D CT Model
where Io is the initial intensity of the X-ray; I′ is the intensity of the through the material; μ is the attenuation coefficient of the mat penetration distance of the ray.
Three sets of phased CT scan tests were performed on the Sihe the loading process, with an axial loading rate of 0.2 mm/min and a C 30 μm. The axial stress-strain data were recorded, and the axia maintained constant during the CT scan. Figure 1 presents CT images showing highly discernible mine Using the minerals as feature points, the Z-axis position in the CT im to delineate the upper and lower boundaries of the 3D model. Thre was performed on the CT data to determine the boundaries of t accurately delineating the range of the 3D CT model at each stage.

Displacement Field Algorithm and ROI Delineation
The feature points were used to select a region of interest (ROI) each stage, and the 3D displacement field was calculated. A spat algorithm was used to measure the displacement field of the 3D R loading conditions with high precision. The ROI was determined to mm × 4 mm × 2 mm, considering the computational burden, the size minerals and the excessive number of voxels (more than 1 billion) in t 2 shows the selected ROIs, where ROI1, ROI2, and ROI3 correspond t primary fractures, coal only, and minerals, respectively.

Displacement Field Algorithm and ROI Delineation
The feature points were used to select a region of interest (ROI) for the 3D model at each stage, and the 3D displacement field was calculated. A spatial image matching algorithm was used to measure the displacement field of the 3D ROI under different loading conditions with high precision. The ROI was determined to be a cuboid of size 4 mm × 4 mm × 2 mm, considering the computational burden, the size of the fractures and minerals and the excessive number of voxels (more than 1 billion) in the 3D model. Figure 2 shows the selected ROIs, where ROI1, ROI2, and ROI3 correspond to regions containing primary fractures, coal only, and minerals, respectively.

Deformation and Failure Features
Considering the large volume of data, Sihe sample 1 was analysed in this study. Figure 3 shows clear anisotropic features for the failure of the sample along the Zaxis at the post-peak stage. Figure 4 shows that the distribution of primary fractures and minerals in the sample had a determining effect on the failure features. Considering the post-peak 3D model in conjunction with the CT images in Figures 3 and 4 shows that a slip instability occurred in the sample along the mineral region. The yellow circles in Figure 4 show the damage in the range containing minerals

Deformation and Failure Features
Considering the large volume of data, Sihe sample 1 was analysed in this study. Figure 3 shows clear anisotropic features for the failure of the sample along the Z-axis at the post-peak stage. Figure 4 shows that the distribution of primary fractures and minerals in the sample had a determining effect on the failure features. Considering the post-peak 3D model in conjunction with the CT images in Figures 3 and 4 shows that a slip instability occurred in the sample along the mineral region. The yellow circles in Figure 4 show the damage in the range containing minerals.

Deformation and Failure Features
Considering the large volume of data, Sihe sample 1 was analysed in this study. Figure 3 shows clear anisotropic features for the failure of the sample along the Zaxis at the post-peak stage. Figure 4 shows that the distribution of primary fractures and minerals in the sample had a determining effect on the failure features. Considering the post-peak 3D model in conjunction with the CT images in Figures 3 and 4 shows that a slip instability occurred in the sample along the mineral region. The yellow circles in Figure 4 show the damage in the range containing minerals (a) (b)    Figure 5 shows the stress-strain curve of Sihe sample 1 at the loading stages corresponding to CT scans and the red numbers correspond to the stages of CT scans. The CT images of the upper and lower end faces and the middle position at each stage were analysed to determine the changes in area and roundness.
The roundness is calculated as follows: where S is the image area and Lmajor is the length of the major axis of the smallest circumscribed ellipse fitted to the image.
The roundness quantifies the change in the shape of a CT image, indicating the degree of closeness to a perfect circle. A roundness of unity indicates that the image is a perfect circle. In Figure 6, the points of each curve correspond to the loading stage (red dots) in Figure 5.    Figure 5 shows the stress-strain curve of Sihe sample 1 at the loading stages corresponding to CT scans and the red numbers correspond to the stages of CT scans. The CT images of the upper and lower end faces and the middle position at each stage were analysed to determine the changes in area and roundness.
The roundness is calculated as follows: where S is the image area and Lmajor is the length of the major axis of the smallest circumscribed ellipse fitted to the image.
The roundness quantifies the change in the shape of a CT image, indicating the degree of closeness to a perfect circle. A roundness of unity indicates that the image is a perfect circle. In Figure 6, the points of each curve correspond to the loading stage (red dots) in Figure 5. The CT images of the upper and lower end faces and the middle position at each stage were analysed to determine the changes in area and roundness.
The roundness is calculated as follows: where S is the image area and L major is the length of the major axis of the smallest circumscribed ellipse fitted to the image. The roundness quantifies the change in the shape of a CT image, indicating the degree of closeness to a perfect circle. A roundness of unity indicates that the image is a perfect circle. In Figure 6, the points of each curve correspond to the loading stage (red dots) in Figure 5. The calculated roundness was a monotonically decreasing function, and the change in the roundness was generally consistent with the change in the area. The shape and area of each part exhibited a large degree of variation during the instability. The coefficient of variation of the image roundness was calculated to be 0.0473, 0.0097, and 0.0061 for the three slices, that is, the roundness of the lower slice was significantly larger than those of the middle and upper slices. The anisotropy of the deformation and failure features of the sample were very pronounced at the macroscopic scale.

Results for 3D Displacement Field
The ROIs for each scanning stage were compared with that of the initial state by using the digital volume correlation module of the Avizo software (i.e., scanning stage 1). Seven sets of calculation results were obtained for each ROI, as shown in Figure 7. The calculated roundness was a monotonically decreasing function, and the change in the roundness was generally consistent with the change in the area. The shape and area of each part exhibited a large degree of variation during the instability. The coefficient of variation of the image roundness was calculated to be 0.0473, 0.0097, and 0.0061 for the three slices, that is, the roundness of the lower slice was significantly larger than those of the middle and upper slices. The anisotropy of the deformation and failure features of the sample were very pronounced at the macroscopic scale.

Results for 3D Displacement Field
The ROIs for each scanning stage were compared with that of the initial state by using the digital volume correlation module of the Avizo software (i.e., scanning stage 1). Seven sets of calculation results were obtained for each ROI, as shown in Figure 7. Figure 7a corresponds to the displacement field of the second loading stage, and so on. The calculated roundness was a monotonically decreasing function, and the change in the roundness was generally consistent with the change in the area. The shape and area of each part exhibited a large degree of variation during the instability. The coefficient of variation of the image roundness was calculated to be 0.0473, 0.0097, and 0.0061 for the three slices, that is, the roundness of the lower slice was significantly larger than those of the middle and upper slices. The anisotropy of the deformation and failure features of the sample were very pronounced at the macroscopic scale.

Results for 3D Displacement Field
The ROIs for each scanning stage were compared with that of the initial state by using the digital volume correlation module of the Avizo software (i.e., scanning stage 1). Seven sets of calculation results were obtained for each ROI, as shown in Figure 7. Figure 7a corresponds to the displacement field of the second loading stage, and so on.   Figure 7a shows the displacement field at the compaction stage: the overall displacement was low, and an analysis of the vector direction shows that the transverse displacement was larger than the longitudinal displacement. Comparing Figure 7b with Figure 7a show that the vector angle was closer to the longitudinal direction than the transverse direction, indicating a significant compression effect of the load. Compared to the displacement field in Figure 7a,b, the results in Figure 7c,d show little change, whereas those in Figure 7e show a pronounced change. The norm of the vector in and near the fracture region is significantly larger than that in fracture-free regions, indicating a significant increase in the displacement in the fracture region. Figure 7f corresponds to the peak stress state: a large displacement was distributed over a significantly wide range, including regions near the fractures and at the boundary. This result shows that under the peak strength loading, the internal displacement of the coal sample increased and was distributed over a wide range. The most significant change in the vector field is shown in Figure 7g: the norm of the vector field increases dramatically, especially over the fracture region and the nearby area, where the displacement generally exceeds 300 µm. The displacement in fracture-free regions generally remains larger than 150 µm.
The displacement fields of ROI2 and ROI3 were calculated, and the post-peak displacement fields were comparatively analysed. Figure 8 shows significant differences among the post-peak displacement field features of the different ROIs. The displacement field of ROI2 (basically coal) was most uniform, with a displacement of approximately 150 µm and large displacements in only small well-defined regions. The distribution features of the displacement field of ROI3 (containing minerals) were similar to those of ROI1 at the corresponding stages: the displacement was large over the region near the minerals, and smaller in the nearby coal region.
Mathematics 2022, 10,2512 8 of 17 Figure 7a shows the displacement field at the compaction stage: the overall displacement was low, and an analysis of the vector direction shows that the transverse displacement was larger than the longitudinal displacement. Comparing Figure 7b with Figure 7a show that the vector angle was closer to the longitudinal direction than the transverse direction, indicating a significant compression effect of the load. Compared to the displacement field in Figure 7a,b, the results in Figure 7c,d show little change, whereas those in Figure 7e show a pronounced change. The norm of the vector in and near the fracture region is significantly larger than that in fracture-free regions, indicating a significant increase in the displacement in the fracture region. Figure 7f corresponds to the peak stress state: a large displacement was distributed over a significantly wide range, including regions near the fractures and at the boundary. This result shows that under the peak strength loading, the internal displacement of the coal sample increased and was distributed over a wide range. The most significant change in the vector field is shown in Figure 7g: the norm of the vector field increases dramatically, especially over the fracture region and the nearby area, where the displacement generally exceeds 300 μm. The displacement in fracture-free regions generally remains larger than 150 μm.
The displacement fields of ROI2 and ROI3 were calculated, and the post-peak displacement fields were comparatively analysed. Figure 8 shows significant differences among the post-peak displacement field features of the different ROIs. The displacement field of ROI2 (basically coal) was most uniform, with a displacement of approximately 150 μm and large displacements in only small well-defined regions. The distribution features of the displacement field of ROI3 (containing minerals) were similar to those of ROI1 at the corresponding stages: the displacement was large over the region near the minerals, and smaller in the nearby coal region.   A statistical analysis showed that the displacement field of each ROI follows a normal distribution. The displacement fields for ROI1-3 were fitted to obtain the following equation: where x is the position parameter for the normal distribution and corresponds to the mean displacement, and w describes the dispersion of the data. The fitting parameters x and w of    Table 1 lists the mean displacements for the post-peak stage and the full stress-strain process of ROI1-3.  A statistical analysis showed that the displacement field of each ROI follows a normal distribution. The displacement fields for ROI1-3 were fitted to obtain the following equation: where x is the position parameter for the normal distribution and corresponds to the mean displacement, and w describes the dispersion of the data. The fitting parameters x and w of the displacement field of ROI1-3 were statistically calculated at each stage (stage 2 to 8), and the results are presented in Figure 12.  Table 1 lists the mean displacements for the post-peak stage and the full stress-strain process of ROI1-3.  A statistical analysis showed that the displacement field of each ROI follows a normal distribution. The displacement fields for ROI1-3 were fitted to obtain the following equation: where x is the position parameter for the normal distribution and corresponds to the mean displacement, and w describes the dispersion of the data. The fitting parameters x and w of the displacement field of ROI1-3 were statistically calculated at each stage (stage 2 to 8), and the results are presented in Figure 12.  Table 1 lists the mean displacements for the post-peak stage and the full stress-strain process of ROI1-3.  Figure 12 shows the mean and dispersion of the displacement field versus the axial strain for ROI1-3. Considering these results in conjunction with those presented in Table 1 shows that (1) there is a turning point in the statistics (i.e., the mean and dispersion) of the displacement field at the compaction stage; (2) the displacement increases monotonically with the load after the compaction stage; (3) the dispersion in the statistical results of the displacement field increases steadily with the loading, reflecting an increase in the number of extreme displacements in the sample; and (4) the displacement fields of different ROIs at the same stage differ significantly.

Analysis of Simulation Results
Mechanical tests using the MTS816 rock test system were used to measure the Young's modulus and Poisson's ratio of each sample. The grey value ranges for different matrices (coal, minerals, and fractures) were determined, and the mechanical test results were used to set the variable physical parameters (Young's modulus and Poisson's ratio) for different matrices, as shown in Figure 13. A meshfree numerical simulation of the 3D CT model under mechanical loading (with stress values corresponding to stages 2-7 in Figure 5) was carried out by using VG software.
in the number of extreme displacements in the sample; and (4) the displacement fields of different ROIs at the same stage differ significantly.

Analysis of Simulation Results
Mechanical tests using the MTS816 rock test system were used to measure the Young's modulus and Poisson's ratio of each sample. The grey value ranges for different matrices (coal, minerals, and fractures) were determined, and the mechanical test results were used to set the variable physical parameters (Young's modulus and Poisson's ratio) for different matrices, as shown in Figure 13. A meshfree numerical simulation of the 3D CT model under mechanical loading (with stress values corresponding to stages 2-7 in Figure 5) was carried out by using VG software. The 3D CT model was directly simulated to obtain the 3D shear stress, principal stress, and displacement fields under different loads, as shown in Figures 14-16. Combining these results with those presented in Figure 17 clearly shows that the mineral distribution had a significant effect on the stress field. The 3D CT model was directly simulated to obtain the 3D shear stress, principal stress, and displacement fields under different loads, as shown in Figures 14-16. Combining these results with those presented in Figure 17 clearly shows that the mineral distribution had a significant effect on the stress field.        The shear stress distribution was bimodal (except in the fracture region): the minimum shear stress was 2.517 MPa, and 2.3% of the shear stress ranged between 2.517 and 3 MPa. More than 95% of the shear stresses were between 3 and 6.4 MPa, and only 1.6% of the shear stresses were greater than 6.4 MPa, with the maximum shear stress being 11.6 MPa.
The maximum principal stress followed a trimodal distribution. Only 0.89% of the maximum principal stresses were below 6 MPa, and 3.57% of the principal stresses were between 6 and 8.57 MPa. The percentage of maximum principal stresses between 8.5 and 8.9 MPa increased rapidly, and 79.5% of the maximum principal stresses were between 8.9 and 10.7 MPa. The maximum principal stresses in 14.34% of the regions were between 10.7 and 12.87 MPa, whereas the maximum principal stresses were higher than 12.87 MPa in only 1.67% of the regions.
The changes in the proportions of the maximum principal stress and shear stress between 10 and 20 MPa and above 20 MPa for each stress state were statistically analysed. Figure 19 shows that the proportions of the maximum principal stress and the shear stress above 20 MPa in the sample increased monotonically with the load. For loads above 15.46 MPa, the proportion of the maximum principal stresses higher than 20 MPa increased  The shear stress distribution was bimodal (except in the fracture region): the minimum shear stress was 2.517 MPa, and 2.3% of the shear stress ranged between 2.517 and 3 MPa. More than 95% of the shear stresses were between 3 and 6.4 MPa, and only 1.6% of the shear stresses were greater than 6.4 MPa, with the maximum shear stress being 11.6 MPa.
The maximum principal stress followed a trimodal distribution. Only 0.89% of the maximum principal stresses were below 6 MPa, and 3.57% of the principal stresses were between 6 and 8.57 MPa. The percentage of maximum principal stresses between 8.5 and 8.9 MPa increased rapidly, and 79.5% of the maximum principal stresses were between 8.9 and 10.7 MPa. The maximum principal stresses in 14.34% of the regions were between 10.7 and 12.87 MPa, whereas the maximum principal stresses were higher than 12.87 MPa in only 1.67% of the regions.
The changes in the proportions of the maximum principal stress and shear stress between 10 and 20 MPa and above 20 MPa for each stress state were statistically analysed. Figure 19 shows that the proportions of the maximum principal stress and the shear stress above 20 MPa in the sample increased monotonically with the load. For loads above 15.46 MPa, the proportion of the maximum principal stresses higher than 20 MPa increased rapidly from 1.628% to 51.683%. At the peak stress, the maximum principal stress exceeded 20 MPa in 96.35% of the regions. The proportion of regions with shear stresses above 20 MPa increased rapidly for loads above 15.46 MPa, but the proportion itself remained very low. The proportion of the shear stresses between 10 and 20 MPa increased significantly with the load from 1.5% to 43.457% to 82.739%. Two points in the 3D model were selected for each matrix (coal, minerals, and near fractures), and Figure 20 shows the corresponding maximum principal stress and shear stress. The different matrices within the sample exhibited significantly different stress states. The maximum principal stresses in the mineral matrix were higher than those in the coal matrix, indicating stress concentration in the mineral regions. By comparison, the stresses in the regions near the fractures were all low.

Conclusions
The deformation features of coal under stress at macro-and mesoscales was investigated, and a meshless simulation was performed using a 3D CT model to obtain the 3D principal stress, shear stress, and displacement fields. The main conclusions are given below. Two points in the 3D model were selected for each matrix (coal, minerals, and near fractures), and Figure 20 shows the corresponding maximum principal stress and shear stress. The different matrices within the sample exhibited significantly different stress states. The maximum principal stresses in the mineral matrix were higher than those in the coal matrix, indicating stress concentration in the mineral regions. By comparison, the stresses in the regions near the fractures were all low. Two points in the 3D model were selected for each matrix (coal, minerals, and near fractures), and Figure 20 shows the corresponding maximum principal stress and shear stress. The different matrices within the sample exhibited significantly different stress states. The maximum principal stresses in the mineral matrix were higher than those in the coal matrix, indicating stress concentration in the mineral regions. By comparison, the stresses in the regions near the fractures were all low.

Conclusions
The deformation features of coal under stress at macro-and mesoscales was investigated, and a meshless simulation was performed using a 3D CT model to obtain the 3D principal stress, shear stress, and displacement fields. The main conclusions are given below.

Conclusions
The deformation features of coal under stress at macro-and mesoscales was investigated, and a meshless simulation was performed using a 3D CT model to obtain the 3D principal stress, shear stress, and displacement fields. The main conclusions are given below.

1.
Primary fractures and minerals have a determining influence on the deformation and failure features of coal and are the key factors producing anisotropic coal deformation. The primary fractures propagate as compression increases and crisscross new fractures.

2.
The ROI displacement field data of different matrices in the coal sample follow a normal distribution, and the differences in the corresponding deformations under load are statistically significant. The deformation at the primary fractures increases substantially with stress, and the mineral region is also prone to large deformation from stress concentration. The increase in the displacement is promoted by the presence of fractures and inhibited in the presence of minerals.

3.
A highly complex coal stress state forms under an applied stress, and the distribution of minerals and fractures significantly affects the stress field distribution. The presence of minerals and fractures produces a prevalent shear stress in the coal that is mainly concentrated in the vicinity of where these entities are located and is highly unfavourable for sample stability.

4.
To ensure safe field operation, the mineral distribution features in coal seams and the development level and distribution features of primary fractures should be determined. Targeted measures can then be taken to improve the stress state and enhance the safety level.