Abstract
The erosion of soluble rock and transformation of groundwater result in the high irregularity of cavities in tunnel. At present, however, karst cavities are mainly simplified as circular, rectangular, or elliptical shape in the numerical simulation. The purpose of this paper is to propose a new method to analyze the stability of irregular cavities. First of all, we used the drilling laser scanning method to reconstruct the three-dimensional point clouds model of irregular cavities. Furthermore, we proposed the method of determining the point density to reduce the computational error under the premise of ensuring the accuracy in engineering scale. We also developed the Geomagic-COMSOL coupling numerical method to conduct the stability analysis of irregular cavities. The results demonstrated that the geometrical shape of karst cavities has great effects on the deformation and mechanical properties of the surrounding rock. The displacement and equivalent plastic strain of simplified cavities exhibited symmetric characteristics, while those of irregular cavities are highly randomly distributed. We noted that a greater displacement value and more intensive plastic strain can be triggered by irregular cavities shape, compared with the simplified cavity shape. The results also showed that the larger displacement occurred at the long-edge part, while the plastic zone was concentrated at the sharp turning angle of the irregular cavities.
1. Introduction
Karst cavities will deteriorate the mechanical properties and bearing capacity of the surrounding rock during tunnel construction, resulting in rock instability and failure and then causing hazards such as collapse and water inrush [1–3]. The development of cavities is affected by many factors such as stratum orientation, geological structure, fracture degree, and groundwater action. Due to the complex mechanism of karst, the geometrical shape of the karst cavity can be anisotropic and highly irregular. Because of the difficulties in describing the geometrical shape of a natural irregular cavity, circular, rectangular, or elliptical geometries are assumed in most related studies. However, these simplifications do not reflect the real shape of cavities in nature. Therefore, it can be of great significance to propose a practical numerical simulation method to study the effect of the irregular cavity geometric characteristics on the surrounding rock.
Numerical analyses have been extensively employed for analyzing the effect of karst cavities on the surrounding rock stability and several factors including size, location, water pressure, and rock mechanical parameters are considered [4, 5]. Nevertheless, a few scholars have considered the cavity shape effect on the mechanical properties of surrounding rock. Abbo et al. [6] pointed out that the tunnel shape had a great impact on the rock mechanical properties. Du et al. [7] studied the influence of cavity geometry on the failure process of sandstone material. The results showed that the sandstone samples had the lowest bearing capacity when they contain a square cavity, and the highest bearing capacity when containing a rhombus cavity. Zhao et al. [8] emphasized the fact that the cavity shape needs to be considered in the simulation and discussed the difficulty of quantitative description of irregular cavities shape. Liu et al. [9] adopted the 3D laser technology to obtain the distribution of the goaf and proposed the Surpac-Flac coupling numerical interface through node transformation to analyze the goaf stability. Li and Qin [10] have constructed a real 3D mining patio model and analyzed the variation law of displacement and stress of surrounding rock. Zhang et al. [11, 12] pointed out that the bearing capacity of pile foundation was related to the shapes of karst cavities, thus suggesting that the karst cavities shape should be considered into the calculation. Zhang et al. [13] analyzed the effect of karst cavities with different shapes on the bearing capacity of rock-socketed piles and presented a suggestion table of cavity shape factors. However, only the flatness parameter was considered, making the table lack of comprehensive quantitative description.
At present, a few studies on the numerical simulation of irregular cavities have been conducted, and the reasons can be classified into the following three aspects. Firstly, it is difficult to obtain the real cavities geometrical shape in practical engineering. The mainly used methods including the comprehensive geophysical prospecting and drilling are both qualitatively and semiquantitatively. Therefore, detailed cavities parameters including the specific location, shape, size, and volume are not accurate. Moreover, the results are mostly 2D data with low visualization. The rapid development of noncontact measurement technology provides a solution for the refined exploration of karst cavities. Currently, CMS (Cavity Monitoring System) and C-ALS (Cavity Auto Laser Scanner) are the most extensively used laser scanning method for karst cavities [14–17]. Secondly, the cavity models obtained were only visualized and cannot be effectively converted into computable numerical models. Thirdly, due to the complex mechanism of karst, the geometrical shape of karst cavity can be anisotropic and highly irregular. Because of the difficulties in describing the geometrical shape of a natural irregular cavity, circular, rectangular or elliptical geometries are assumed in most related studies. It is obvious that these simplifications cannot reflect the irregular characteristics of natural cavities [18, 19]. Thus, the errors during the mesh generation process occurred, making it infeasible to conduct the following simulation procedure.
Based on the above analysis, we proposed a construction method of 3D irregular cavities points model using laser scanning. Additionally, the Geomagic-COMSOL coupling numerical interface to conduct the irregular cavities stability was further developed. Based on the approach proposed, the shape effects of irregular and simplified cavities on the deformation and mechanical properties of surrounding rock were analyzed. As will be shown, the proposed method is robust and efficient in simulation analysis that represents a realistic cavity. It hence offers a possible way to further quantitatively study the cavity shape effect on the mechanical properties of surrounding rock.
2. 3D Model Construction of Irregular Cavity
2.1. Irregular Cavities Modelling Using Laser Scanning
The laser scanning method is an efficient way to obtain the point coordinates of the cavity with high precision [20]. In this paper, we used the Cavity Auto Scanning Laser System (C-ALS) to perform cavity detection, which is the first commercial product using drilling laser scanning technology [15]. The diameter of the miniaturized laser range finder is only 5 cm, which can enter the internal cavity through the borehole to obtain points coordinates and reconstruct a 3D model, as shown in Figure 1. The probe has a 3D navigation system, which can record the depth in real time to determine the precise position of the cavity. The scanning probe rotated using a 3D mechanical servo system, measuring the 3D points coordinates of the cavity. Finally, the data processing software was further adopted to reconstruct the visualized cavity points cloud model with irregular shape including the location, shape, and volume.

(a)

(b)
The laser module is equipped with a laser emitter and a laser receiver. The distance can be calculated by recording the time and angle that each laser pulse emitted from the probe to the cavity surface; then, the 3D cavity coordinates can be obtained. The distance between the laser emitter and the karst cavity wall is expressed aswhere is the distance between the laser emitter and the karst cavity wall; is the laser propagation velocity; and is the time taken from laser emission to reception.
The relative coordinate of the cavity point can be solved based on the horizontal laser rotation angle and vertical laser rotation angle , as shown in Figure 2. It is presumed that points and are the laser emitting point and the cavity point, respectively. Thus, the relative coordinate of point can be expressed as

(a)

(b)
Then, the absolute coordinate of point can be further calculated using
The key of splicing points in the multiposition is the transformation relationship between different coordinate systems. Taking the two-station detection as an example, the two stations are and and the transformation relationship can be expressed aswhere the translation matrix can be denoted by the translation amount , , and on the 3D coordinate axis:
The above formulas contain 7 parameters including rotation parameters , translation parameters , and scale factor. The scale factor was equal to 1 without scale change. The rotation matrix along the coordinate axis is
The coordinates of the laser emission points in two stations are known, which means that the corresponding translation matrix is a known variable. Hence, the point coordinates in coordinate system can be calculated through equation (4). Based on the above analysis, the points cloud data of multiple stations can be converted and spliced to obtain points cloud model of large-scale cavities with an irregular shape.
2.2. Determining Method of Sampling Points Spacing
Theoretically, the points cloud model of the karst cavity should be a continuous and smooth surface. However, the actual one obtained from practical engineering often has error points. The main reasons are as follows: firstly, the cavity contains irregular fillings, corners, and other obstacles. Under such conditions, the probability of detection error during the operation increases. Secondly, the dim and dusty environment in the cavity reduces the detection accuracy of the points coordinate. Thirdly, during the mechanical process of laser scanning, a slight vibration of the laser emitter may cause the laser to undergo strong diffraction, generating some local points groups away from the based points model, that is, the “error points.”
Figure 3(a) shows the typical points model of irregular karst cavities and there are obvious isolated point groups and double-layer points clouds marked by red circle. Isolated point groups are prone to geometry spikes when performing points cloud encapsulation, and local irregular geometries probably appear during double-layer points cloud encapsulation. The encapsulation of the above two error points can generate easily distorted geometries, making it difficult to produce a closed cavity geometric surface, as shown in Figure 3(b).

(a)

(b)
In Geomagic-Studio processing platform, the probability of error points and geometric distortion can be reduced by setting the sensitivity or manually deleting the obvious error points, especially including the isolated points groups and the affiliated error point group. The points model after removing the significant error points and the cavity model after the noise-reducing encapsulation is shown in Figure 4. Comparing with Figure 3(b), the spikes and affiliated distorted geometries in Figure 4(b) reduced significantly, greatly improving the accuracy and reliability of the irregular cavity model. However, a certain number of geometric distortions marked by yellow area still appear around the surface of irregular cavities, making it difficult to form a single closed layer. Therefore, unified points reduction should be conducted before the numerical simulation of irregular cavities. Nevertheless, the points reduction will inevitably affect the cavity modelling accuracy. The determining of the sampling points interval is of significance under the premise of ensuring the engineering calculation accuracy.

(a)

(b)
Local distortion optimization of the cavity model can be performed by setting a certain ratio or spacing value for sampling, under the premise of ensuring that the geometrical model can be packaged into a continuous closed surface and then transformed into a solid model. There are four main methods for points sampling including curvature sampling, unified sampling, random sampling, and raster sampling. Unified sampling, which is the extensively used method, can reduce the number of points clouds on flat surfaces and optimize the local complex points clouds. In this paper, the unified sampling method was adopted to perform the error points reduction. We proposed the determining formula of the minimum distance between the horizontal and vertical directions of the karst cavity points, which can be regarded as a criterion for points cloud sampling. The cavity coordinates are shown in Figure 5(a), where point O is the vertical projection of the laser emission point on the horizontal plane; points P, M, and H are the continuous points in the vertical direction; and is the adjacent angle along the vertical direction. Theoretically, the connecting line of points P, M, and H was an irregular curve. However, the idea of a straight replacing curve was adopted. PM and MH lines were the spacing between adjacent points in the vertical direction. The vertical adjacent spacing can be expressed as

(a)

(b)

(c)
Similarly, the horizontal adjacent spacing can be calculated using
Due to the irregularity of the cavity shape, the distance between the cavity surface points and the center point B was not fixed. However, from the perspective of statistical sampling, the average value of the cavity in a 3D scale can be taken as the diameter and the orthogonal sections are shown in Figure 6. The single-dimensional projected area of the cavity is , equivalent to the circle with the same area. Then, the equivalent radius can be calculated asthat is, the average radius of the irregular cavity is

(a)

(b)

(c)

(d)
The average distance between points in the vertical direction can be expressed as
The average distance between points in the horizontal direction is
The average distance between the points clouds is proportional to the average cavity radius and the progressive rotation angle. Since the average distance between the vertical and horizontal points was obtained, the average spacing of the cavity points can be calculated as
2.3. Cavity Points Encapsulation and Reconstruction
The statistical results show that the error points reduction ratio can be low when the sampling interval was less than 1 time of the average spacing, resulting in the generation of self-crossing surface, isolated points, and distortion points. While the sampling interval was larger than 2 times the average spacing, the shape cannot exactly reflect the geometrical characteristic in practical engineering. Therefore, when the sampling interval ranged from 1 to 2 times the average spacing, the success rate of points encapsulation improved greatly.
For example, the average radius of a cavity model is 2.67 m; its vertical progressive rotation angle is 0.00174 rad; its vertical average points spacing is 0.00465 m; its horizontal progressive rotation angle is 0.0035 rad; and its horizontal average points spacing is 0.0093 m. The average distance between the points is 0.00698 m. Sampling tests with different spacing were conducted on the above cavity model. The sampling intervals were 0.5 times, 1 time, 2 times, 5 times, and 10 times the average points distance, and the cavity points and the packaged model after unified sampling are shown in Figure 7. From Figure 7(a), we can see that the cavity model was consistent with the actual cavity shape. However, there was still plenty of geometric distortion in some irregular regions with large tortuosity, resulting in an error in subsequent model meshing. The disparity between the model and the actual one increased as the sampling spacing increased. However, the corresponding triangle unit division was rougher and the geometry smoothness reduced.

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

(i)

(j)
3. Geomagic-COMSOL Coupling Numerical Interface of Irregular Cavities
COMSOL Multiphysics® is versatile multifield coupling simulation software, and the geotechnical module has been adopted by many scholars to investigate geotechnical engineering stability. For the stability analysis of irregular cavities, the key problem in numerical simulation is to transform the “visualized” cavity model into an entity model with “calculability” and further conduct the grid meshing effectively. The karst cavity points can be encapsulated by Geomagic Studio to form a closed surface. Additionally, COMSOL software reserves a points scanning model interface, which can support the divided grid and the STL format model. The specific flow chart is shown in Figure 8.

Taking a typical karst cavity as an example, the cavity points model simplified with 1.5 times the sampling spacing was imported into Geomagic software. A closed cavity surface model was reconstructed, instead of an entity karst cavity, as shown in Figure 9(a). Therefore, it was further converted into a STL format file and then imported into the geometry module in COMSOL software through the Geomagic-COMSOL interface. The cavity entity was shown in Figure 9(b). It needs to be noted that the hexahedral element for meshing was infeasible because of the extreme irregularity of the karst cavity. Therefore, the mesh generation of the cavity model used tetrahedral elements. The model was sparsely meshed in a smooth area, and the meshing accuracy was improved at the extremely irregular part, as seen in Figure 9(c).

(a)

(b)

(c)
4. Results between Irregular Cavities and Simplified Model
In order to investigate the difference between simplified and irregular karst cavity, the effect of karst cavities with a different shape on the mechanical properties of the surrounding rock was analyzed. Additionally, we have investigated the mechanics characteristics of karst cavities with extremely irregular shapes.
4.1. Geometric Model
Among karst cavities with irregular shapes, sphere-like and cube-like cavities were selected as examples to conduct numerical simulation and comparison with simplified cavity models. Additionally, the cavity model with an extremely irregular shape was selected for further analysis. The simplified model scale was determined based on the principle of volume equal. Assuming that the volume of the real cavity is , then, the formula for determining scales of sphere and cube can be expressed aswhere is the volume of the irregular cavity; is the radius of the simplified sphere; and is the side length of the standard square model.
Typical sphere-like, cube-like, and extremely irregular cavities were constructed by laser scanning in practical engineering. The geometrical models of sphere-like and cube-like cavities are shown in Figures 10(a) and 10(c), respectively. The volume of the cube-like cavity was 7.5184 m3 and the radius of the sphere cavity was 1.196 m using (14), as shown in Figure 10(b). The volume of the cube-like cavity was 0.207 m3 and the side length was 0.592 m in equation (15), as shown in Figure 10(d). Additionally, the two cavity models with extremely irregular shapes were also taken into consideration in this paper, as shown in Figures 10(e) and 10(f), respectively.

(a)

(b)

(c)

(d)

(e)

(f)
4.2. Input Parameters and Boundary Conditions
In the analysis of simplified and irregular shape cavities, the structural mechanics and Darcy modules were used to investigate the fluid-solid coupling steady stability. The karst model was imported into the abovementioned Geomagic-COMSOL coupling method as shown in Figure 10. The geotechnical material ranging in 3∼5 times of the cavity diameter was guaranteed. Due to the high irregularity of the karst cavity, the model was meshed using tetrahedral elements generation. In terms of the boundary conditions, a fixed constraint was applied at the bottom of the model body. A roller support was adopted to the surrounding boundary. A free boundary was formed at the upper portion. Moreover, a fixed water pressure loads on the boundary of the cavity. The yield criterion of the rock mass was Drucker–Prager matching the Mohr–Coulomb, which can be expressed as
The parameters input in the numerical simulation are shown in Table 1.
4.3. Numerical Simulation Results of Sphere and Sphere-Like Cavity Models
The geometric models of the sphere and sphere-like cavities are shown in Figures 11(a) and 11(b), respectively. Both models have the same element grid size, boundary conditions, rock parameters, yield criteria, and water pressure in the karst cavity. The only difference was the geometric shape of the cavity, in order to compare the shape effects on the mechanical properties. The displacement and plastic zone of the surrounding rock were the main research indexes in this paper.

(a)

(b)
Figure 12 plots the displacement field of the surrounding rock with spherical and the sphere-like cavity. The constant water pressure in the cavity loaded on the cavity surface wall as boundary loading, thus exerting force on the surrounding rock mass. From Figures 12(a) and 12(b), we can see that the displacement in the upper part presented an upward trend, due to the hydraulic pressure, while the displacement in other parts went downward. The obvious difference was the distribution of the displacement field around the irregular cavity. Figures 12(c) and 12(d) present clearly the 2D displacement field at the selected sections. The displacement field with the regular sphere model was symmetric, while the one of the sphere-like cavity was randomly distributed. The maximum displacement value with irregular cavities was approximately 2 times that for regular cavities.

(a)

(b)

(c)

(d)
Figure 13 plots the effective plastic strain in 3D and 2D scales at the selected sections. It can be seen from Figures 13(a) and 13(b) that the effective plastic strain region around the cavity was uniformly distributed due to the water pressure in the cavity. The effective plastic strain region of the regular sphere model was circularly symmetrically distributed, while the one of the sphere-like cavities was dissymmetric, especially for the extremely irregular part. In terms of numerical values, the maximum plastic strain of irregular cavities was approximately 3 times that for regular karst cavities.

(a)

(b)

(c)

(d)
4.4. Numerical Computation Results of Cube and Cube-Like Cavity Models
The geometrical models for the cube and cube-like cavities are shown in Figures 14(a) and 14(b), respectively. The above models had the same element grid size, boundary conditions, mechanics parameters, yield criteria, and water pressure. The only difference lied in the cavity shape.

(a)

(b)
Figure 15 plots the displacement field with the simplified cube and the cube-like cavity in 3D and 2D scale at the section. It can be seen from Figures 15(a) and 15(b) that the displacement in the upper part of the cavity moved upward due to the hydraulic pressure, which was similar to the results of simplified sphere cavities. In terms of displacement values, the variation of displacement around the cube-like cavity was slightly larger than that of the standard cube cavity. It can be observed that the displacement field of the regular cube model presented symmetric, while the one of the cube-like karst cavity showed little symmetry, as shown in Figures 15(c) and 15(d). From the perspective of numerical values, the displacement results were basically consistent, but the distribution areas were totally different.

(a)

(b)

(c)

(d)
Figure 16 plots the effective plastic strain of the surrounding rock triggered by the simplified cube and the cube-like cavities. The effective plastic strain region around the cube cavity was uniformly distributed. From Figures 16(c) and 16(d), we can also see that the effective plastic strain region of the regular cube cavity model was mainly concentrated at the sharp angle and exhibited a standard symmetry distribution. For the cube-like cavity, the effective plastic strain region was also concentrated at the turning angle but was dissymmetric due to the difference in shapes. In terms of numerical values, the maximum value induced by the cube-like cavity was five times that of the regular cavity.

(a)

(b)

(c)

(d)
4.5. Mechanical Properties of Surrounding Rock with Extremely Irregular Karst Cavities
In this paper, two typical cavity models with high irregularity were selected for numerical simulation, in order to reveal the deformation and stress characteristics of rock mass around extremely irregular cavities. The karst cavity model in case 1 was a narrow karst cavity, flat with a relatively smooth boundary shape. The larger displacement region was mainly distributed in the long-side region, while the larger plastic strain region was mainly concentrated on the short-edge with a larger turning angle. In case 2, the model was a triangular similar cavity. Due to the boundary irregularity, the larger displacement was concentrated on the long side. Apart from the plastic zone concentrated at a turning angle, there were plastic zones around the entire cavity boundary, as shown in Figure 17.

(a)

(b)

(c)

(d)
5. Discussion
At present, the shapes of the karst cavities have been simplified into spheres, cubes, or ellipsoids in most numerical simulation, considering the difficulties in modelling real irregular cavities, resulting in low accuracy and reliability of numerical results. In this paper, the noncontact measurement method of laser scanning was adopted to obtain the points model of irregular shape cavities and reconstruct a visual solid model. Additionally, the optimal interval of points sampling was determined and then the Geomagic-COMSOL coupling interface was proposed to conduct simulation analysis, by which, the stability analysis of irregular karst cavities shape has been realized. However, there are still difficulties in the 3D information investigation of irregular karst cavities. The environment inside the cavity is complicated, and problems such as gloom, dust, and water-filled laser attenuation can reduce the points density and reliability. For the cavity model with particularly high irregularity, errors in meshing generation may still remain.
Based on the above analysis, this paper analyses the effect of the cavities shape on the mechanical properties. However, further study should quantitatively determine the response relationship between the mechanical properties and the shape descriptors of the karst cavity. Parameters of 2D or even 3D scale need to be proposed to quantitatively describe the geometric characteristics, which have a certain physical meaning and the number should be as small as possible. Furthermore, the coupling interface proposed in this paper can be regarded as an effective method to study the relationship between critical safety distance and the cavity shape. The safe treatment distance will determine the safety of tunnel excavation and the economic cost. Conversely, if the safety distance is determined too small, high-risk cavities can not be effectively treated, resulting in water inrush, karst collapse, and other engineering hazards. At present, the method for determining the safe treatment distance of karst cavities is still lack of consideration for the irregular shape.
6. Conclusion
(1)In this paper, we proposed a construction method of 3D irregular cavities points model using laser scanning. The sampling interval was determined, ranging from 1 to 2 times the average spacing, and the success rate of points encapsulation improved greatly, under the premise of ensuring the accuracy in engineering scale.(2)The Geomagic-COMSOL coupling interface was proposed to conduct the stability analysis of irregular karst cavities. The method was proved to be robust and efficient in simulation analysis that represents a realistic cavity. It hence offers a possible way to further quantitatively study the cavity shape effect on the mechanical properties of surrounding rock.(3)The effects of simplified and irregular cavity shapes on the deformation and mechanical properties of surrounding rock were analyzed. The results demonstrated that the geometrical shape of karst cavities has great effects on the deformation and mechanical properties of the surrounding rock. The displacement and equivalent plastic strain of simplified cavities exhibited symmetric characteristics, while those of irregular cavities are highly randomly distributed. The results also showed that the larger displacement occurred at the long-edge part, while the plastic zone was concentrated at the sharp turning angle of the irregular cavities.(4)In the future study, the characterization method of the irregular cavity shape needs to be further considered in the stability analysis.Data Availability
Some or all data, models, or code that support the findings of this study are available from the corresponding author upon reasonable request.
Conflicts of Interest
The authors declare no conflicts of interest.
Acknowledgments
This research was supported by the National Natural Science Foundation of China (no. 52009076), Shandong Provincial Key Research and Development Program (no. 2019JZZY010428), The National Science Foundation for Distinguished Young Scholars of China (no. 52025091), and Shandong Province Science Foundation for Young Scholars (no. ZR2020QE290).