Stability Study of the Roof Plate of the Yuanjue Cave Based on the Equivalent Support Stiffness Method

Abstract

As precious cultural heritage sites, the state of preservation of cave temples is closely related to the geological and climatic conditions in which they are located. This paper constructed an analytical method of sized slate stability based on the equivalent support stiffness method. The stability analysis of the roof slab of Yuanjue Cave was carried out by establishing a three-dimensional numerical calculation model. Through comparative analysis of the results of stress and displacement fields under different conditions, the stress and deformation characteristics of the roof slab of Yuanjue Cave were revealed, as well as the study of the main factors affecting the stability of the roof slab of Yuanjue Cave and the key slate to be monitored. The main research results are as follows. The stress deformation of the roof plate of Yuanjue cave is mainly divided into the initial uniform change stage, the medium-term stable change stage or the medium-term accelerated change stage, and the later rapid change stage. With the increase in the number of overhanging and broken slates and the increase in the damage factor of cracked slates, the vertical stress extremum of the stones increases continuously, and the equivalent support stiffness decreases, which aggravates the uneven stress deformation of the roof of the Yuanjue Cave. This study provides a reliable reference basis for the stability analysis and evaluation of the roof slab of a large number of cave temples existing in the Sichuan and Chongqing areas in China.

1. Introduction

Cave temples are a special form of architecture that accompanied the introduction of Buddhism from India to China, with Buddhism as the main line, spreading from west to east and south along the Silk Road with a wide distribution [1,2]. The grotto temples with have a long history and not only constitute the eastward history of Buddhism but also record the historical and cultural characteristics of the Chinese nation in different periods. They are non-renewable resources and precious cultural relics [3]. The Dazu rock carvings are among the rarest surviving artefacts of the cave temples, an art form that continued into the Song dynasty in Sichuan, Chongqing, and other areas of China [4,5]. Cave temples are carved directly into the natural mountain and cliff face [6]. However, as the cliff face of the mountain is also usually the site of the development of the unloading zone, the cave body is usually cut by the fissures of the unloading zone, which, combined with thousands of years of internal and external camp forces, can easily produce unstable characteristics such as toppling, devolution, and landslides, threatening the safety of cave temples’ cultural relics [7,8]. Stone heritage and ancient sites account for over 50% of China’s cultural heritage [5]. In the field of rock stability assessment of cave temples, their heritage properties result in conventional assessment methods often being ineffective, and conditions such as non-destructive means of detection, the high spatial accuracy of assessment models, and precise assessment parameters need to be met as much as possible [9]. Taking Dunhuang Mogao Grottoes as the research objects, Fu et al. [10] analyzed the influencing factors of dynamic damage of the surrounding rock under seismic load and explored the damage of the surrounding rock and its subsidiary structure caused by earthquakes based on the finite element method. Wang et al. [11] explored the damage process of the roof slab of Longyou Grotto and numerical simulation and back analysis of rheologic fracture were carried out by using in situ damage pattern and deformation data. Chen et al. [12] analyzed the stress distribution and deformation characteristics of the surrounding rock mass under different tourist static loads on the long-term stability of the Mogao Grottoes using a finite difference method. The analysis showed that these will produce a stress concentration and large vertical displacement on the cave roof and densely distributed cave area. Over the centuries, with the constant action of natural forces and the erosion of various degradation factors, the side walls and roof of the cave have spalled, the lower support body has broken up, and the boundary conditions of the rock have deteriorated, posing a serious threat to the safety of the cave [13,14,15].

Much research has been done on the conservation of ancient cave temples in terms of the evolution of the surrounding rocks and the multi-source monitoring of the roof slab [16,17,18,19]. Genis et al. [20] used the slope quality rating, rock classification system, and top slab thickness-to-width ratio methods to evaluate the stability of Naga Mas cave. Guo et al. [21] used the finite difference method to analyze the stress distribution and deformation characteristics of the surrounding rock of the Mogao Caves under different static loads. Yang et al. [22] calculated the safe thickness of the cavern roof using the mechanical analysis method, cavern perimeter rock elasticity theory, and Griffith’s criterion. Meng et al. [23] obtained the deformation pattern of the surrounding rock by coupled analysis of the multi-source monitoring data of Yuanjue Cave. Gao et al. [24] adopted the methods of investigation factors, mechanical calculation, and analysis to carry out a study on the mechanism and law of the development of cracks in the rock roof of No. 3 cavern in Longyou Grotto, and found that the main reason for the cracks was the low fracture toughness of the surrounding rock. Brandi [25] predicted the deformation characteristics of a cavern and its roof slab based on finite element analysis and proposed that caverns with a roof slab thickness greater than 8.0 m can ensure the stability of the roof slab under the self-weight of the roof slab and the pile foundation load. Meng et al. [26] established a three-dimensional computational model of Yuanjue Cave to reveal the deformation law of the top plate of Yuanjue Cave under static loading conditions. Tao et al. [27] used a multi-source monitoring system to monitor the top slab of Yuanjue Cave to reveal the deformation pattern. Peng et al. [28] took the reinforcement project of Wan’an Temple Grottoes as an example, and quantitatively analyzed the contributions from the weathering degree of rock, vibration speed of anchor rod drilling, and distance between the anchor hole and grotto cave to the stability of the grotto cave. Meng et al. [29] studied a separate block with a large free space displacement in the Yuanjue Cave and the stress concentration of its surrounding rock, which provides a useful reference for the stability assessment of the surrounding rock in Yuanjue Cave. Huang et al. [30] used single full-scale loading tests, numerical simulations, and gypsum physical modeling tests to reveal the damage mechanisms of man-made caves. Abswlbaset et al. [31] employed both ground magnetic and GPR surveys for the investigation of the Abydos site to detect and determine the depth and the geometry of any unearthed archaeological features. Abswlbaset et al. [32] detected groundwater pathways to monitor their level of rise in Osirion at Abydos Archaeological Site, Sohag, Egypt using the electrical resistivity tomography technique to reduce the risk of deterioration of the stone foundations of the temples and any archaeological features. This technique has proven its dependability in the identification of groundwater aquifers and their distributions. Yang et al. [33] reported the identification of inorganic and organic materials applied in gildings of mural paintings in Kizil Grottoes located in the western part of China. These interesting findings provide scientific evidence for us to understand the materials and techniques of gold gildings on Kizil Grottoes’ murals.

Cave temples are a form of Buddhist temple architecture and important elements of the world’s cultural heritage. Wang et al. [34] conducted a detailed investigations to analyze the mechanisms of the weathering deterioration of the stone relics based on different aspects of geo-environmental characteristics. The weathering sandstone samples were analyzed using X-ray diffraction and X-ray fluorescence, and it was found that the rock outside of the boundary zone had much stronger weathering than that inside. Goh et al. [35] evaluated the stability of Naga Mas natural cave temples by using slope quality rating (SMR), a Q rock mass classification system and the thickness–width ratio of the cave top, and the relationship between the rock quality, Q values, and cave width showed that all parts inside the cave need support except for the cave chamber at the northern part. Wang et al. [36] carried out a preliminary study to assess the condition of the façade of rock-hewn cave temples using in situ non-invasive portable methods commonly used on stone-built heritage to explore the stability of rock-hewn grotto temples. The Buddhist sacred site, Taya Caves, is a gem hidden underground in Yokohama, Japan. Germinario et al. [37] carried out the first scientific investigation of Taya Caves, dealing with the stone properties, deterioration, and the environmental setting. It is considered that rock–water interaction is particularly damaging even considering just the physical mechanisms. Wang et al. [38] proposed a method that is mainly based on the analytic hierarchy process (AHP) to evaluate the preservation state and risk of Mogao cliff, meaning that numerical simulation was conducted to quantitatively evaluate the stability and effectiveness of protective measures for Mogao cliff. Wu et al. [39] investigated the diversity, distributions, ecological functions, and interaction patterns of both the fungal and microalgal (including cyanobacteria and algae) communities on sandstone in Beishiku Temple, located on the ancient Silk Road, by using high-throughput sequencing analysis. Palka et al. [40] examined the archaeological patterning and Indigenous religious beliefs regarding fortifications, sanctuaries, and deity communication. Ethnohistoric and ethnographic information points to the importance of deities and ancestors for community well-being and success in war.

In this paper, the parameters of rock mechanics of Yuanjue Cave were obtained by a field geological survey, engineering geophysical prospecting, and field rock point load test. An equivalent stiffness method applied to the calculation of strip stability was proposed according to the definition of equivalent stiffness [41,42,43]. Combining the stratigraphic distribution characteristics of the cavern and the physical and mechanical parameters of the rock body, the finite difference method was used to simulate the effects of overhanging, broken, and cracked slates on the stress and deformation of the roof slab to refine the evaluation of the stability of the roof slab of Yuanjue Cave and to classify the sized slates that need reinforcement and attention, providing a new idea for the conservation research of similar cave relics.

2. Study Area

Yuanjue Cave is the core cultural heritage site of the Dazu Rock Carving World Heritage Site and has a status and role in the history of cave construction, sculptural art, and religious culture in China that cannot be ignored. Yuanjue Cave is the 29th cave in Dafo Bay, Baodingshan Mountain, facing west to east, with a slightly arched roof, and is 5.75 m in height, 10 m in width, and 12.13 m in depth. There is a canal cut into the entrance of the cave, which is 4.68 m high, 1.64 m wide at the front entrance, 2.49 m wide at the back entrance, 3.92 m deep at the left entrance, and 4.80 m deep at the right entrance (Figure 1). The stone of the cave is gray sandstone, which is soft and easy to weather off, making it difficult to preserve. The overall depth of the cave is long and wide, and the top is prone to collapse, but it is carefully designed with a slightly curved arch span and the use of four stone walls to disperse the forces, reducing the pressure on the rocks at the top of the cave and dispersing the overall forces on the top of the cave, resulting in the internal structure of the cave remaining intact despite hundreds of years of natural environmental changes and erosion.

2.1. Stratigraphic Lithology

The stratigraphy of Yuanjue Cave is typical of the red-bedded sandstone of the Sichuan Basin, which is a river–lake sedimentary environment, with horizontal fine-grained sandstone and mudstone interlayers as the main stratigraphic layers, and local sandstone layers which contain more schistose mica in thin layers. The lithological distribution of the stratigraphy of the south cliff where Yuanjue Cave is located and the relative relationship between Yuanjue Cave and the cliff wall are shown in Figure 2.

The cliff’s stratigraphy is mainly horizontal fine-grained sandstone and mudstone interstratification, which can be divided into 13 geological strata. The natures of hard rock and soft rock are extremely different, and the development of the cliff structure fissure derivation is closely related to the distribution of these soft and weak strata, and the weak layer mainly plays the role of blocking and offsetting the fissure development. The rock structure of the surrounding rock of Yuanjue Cave is dominated by tension-type fissures, and no fault or other large geological structure traces are found to exist. Through on-site geological investigation, there are 23 fissures in total on the south cliff wall where Yuanjue Cave is located. These fissures are mainly damaged by tension, with the outcrop showing tension and the width of the opening ranging from a few millimeters to several meters; among them, the giant fissures are mostly filled (powder and sand left with water migration, etc.), but the fissures with smaller openings are usually filled with less or no filling.

2.2. Top Plate Catastrophe Mode

The characteristics of the roof slab and supporting structure of Yuanjue Cave were obtained through detailed on-site investigation and research, which were greatly influenced by the process of cave excavation and the long-surviving process, and the history of the force evolution of the roof slab of Yuanjue Cave was influenced by the joint action of man and nature. There are several fissures in a direct cutting relationship with Yuanjue Cave, in which the L6 fissure has developed in the top slab of Yuanjue Cave, and the J10 fissure cuts the top slab rock into independent blocks (Figure 3); the volume of this independent block is 256.03 m³ with a total weight of 570.88 tons. There are a large number of overhanging, broken, and fractured slates in Yuanjue Cave, which have a great impact on the stability of the roof of Yuanjue Cave. It is necessary to carry out theoretical calculations and numerical simulations on the support and the effect of vertical stress on the slates.

3. Analytic Methods

To calculate the equivalent stiffness of each group of boulders in Yuanjue Cave, the distributions of boulders in Yuanjue Cave are different, and there are cases of overhanging crushing and fracturing, etc. The stress conditions in the cave are complicated. In this paper, we proposed an equivalent stiffness method for calculating the stability of boulders [44,45]. Based on the definition of equivalent bearing stiffness and the main distribution of supporting boulders in the roof of the cave, we equated the boulders in Yuanjue Cave into series and parallel arrangements and established analytical solutions for the effective bearing stiffness of boulders in series and parallel forms under different contact areas on the upper and lower surfaces of the boulders. The effective bearing stiffness of each stone grouping was obtained after the contact relations of the stones in the cave were fully counted. The effect of the number of overhanging, broken, and damaged slates on the force–deformation characteristics and stability of the top slab of Yuanjue Cave was revealed through numerical simulation, which is shown in Figure 4.

(1)

Calculation of equivalent stiffness of the sized slate

Knowing the strip area A, the strip height H, the strip modulus of elasticity E, and the stress S on the strip, the equivalence of the support stiffness is calculated as follows (1):

$$K=\frac{F}{\Delta H}=\frac{SA}{\epsilon H}=\frac{EA}{H}$$

(1)

The slates in direct contact with the roof slab were grouped according to their distribution (as shown in Figure 5), and the overall stiffness of the slate and all its lower supporting slates is taken as the equivalent supporting stiffness of the group.

The method of direct contact between the strips and the top plate can be mainly divided into two basic forms: series and parallel connection, as shown in Figure 6.

According to the definition of equivalent support stiffness, the calculation formulae are derived as follows:

① Tandem arrangement

$$\mathsf{\Delta }H=\frac{F}{K_1}+\frac{F}{K_2}+\frac{F}{K_j}$$

(2)

Substituting Equation (2) into Equation (1), we obtain

$$K=\frac{K_1{K}_2{K}_j}{K_1{K}_2+K_2{K}_j+K_1{K}_j}$$

(3)

If there is no joint, its equivalent elastic modulus is

$$E=\frac{E_1{E}_2}{{\alpha }_1{E}_2+\left(1-{\alpha }_1\right)E_1}$$

(4)

where, ${\alpha }_1=h_1/(h_1+h_2)$, and $E_i$ is the modulus of elasticity of the 𝑖th strip of stone.

② Parallel arrangement

It is assumed that when arranged in parallel, the deformation of the two stones below is equal, i.e., only vertical displacement of the upper slate occurs, and no deflection occurs. This leads to

$$F=K_2\Delta H_2+K_3\Delta H_2$$

(5)

$$\Delta H_2=\frac{F}{K_2+K_3}$$

(6)

Therefore, the total deformation is

$$\Delta H=\frac{F}{K_1}+\frac{F}{K_2+K_3}+\frac{F}{K_j}$$

(7)

Substituting Equation (7) into Equation (1), we obtain

$$K=\frac{K_1{K}_j\left(K_2+K_3\right)}{(K_1+K_j)+(K_2+K_3)+K_1{K}_j}$$

(8)

If the height of the two stones below, the elastic mode, and the contact area are known, we can obtain

$$\Delta H_2=\epsilon h_2$$

(9)

$$E_2\epsilon A_2+E_3\epsilon A_3=F$$

(10)

$$\epsilon =\frac{F}{E_2{A}_2+E_3{A}_3}$$

(11)

Substituting Equation (11) into Equation (9), we obtain

$$\Delta H_2=\frac{F{h}_2}{E_2{A}_2+E_3{A}_3}$$

(12)

Substituting Equation (12) into Equation (6), the equivalent stiffness of the two stones below can be obtained as follows:

$$K_{2-3}=K_2+K_3=\frac{E_2{A}_2+E_3{A}_3}{h_2}$$

(13)

where $E_2$, $E_3$ are the modulus of elasticity of the slates below, $A_2$, $A_3$ are the contact area of the slate below, and $h_2$ is the height of the slate below.

Both Equation (3) and Equation (8) are derived with the same contact area on the upper and lower surfaces of the sized slates. The contact conditions on the upper and lower surfaces of the strips are not the same and the contact areas vary. Therefore, the effective support stiffnesses of the series and parallel forms need to be analyzed separately for different contact area cases.

③ Series arrangement with different contact areas

In the tandem arrangement, if the total area of the upper slates is $A$, the ratio of the effective contact area of each contact surface to the total area of the upper slates are ${\alpha }_1$, ${\alpha }_2, {\alpha }_3$, and the modulus of elasticity of the slates and joints are $E_1,E_2,E_j$. For the convenience of calculation, it is assumed that the internal stress influence range varies linearly, and the stress is uniformly distributed within the influence range, and its calculation sketch is shown in Figure 7.

The internal stress distribution function of the rock is

$$\sigma \left(y\right)=\frac{F}{A}·\left\{\begin{array}{l}\left(\frac{1}{{\alpha }_2}-\frac{1}{{\alpha }_3}\right)\frac{1}{h_2}y+\frac{1}{{\alpha }_3}, 0<y<h_2\\ \frac{1}{{\alpha }_2}, h_2<y<h_2+h_j\\ \left(\frac{1}{{\alpha }_2}-\frac{1}{{\alpha }_2}\right)\frac{1}{h_1}\left(y-h_2-h_j\right)+\frac{1}{{\alpha }_2},h_2+h_j<y<h_1+h_2+h_j\end{array}\right.$$

(14)

The strain from Equation (14) can be obtained as follows:

$$\epsilon \left(y\right)=\frac{F}{A}·\left\{\begin{array}{l}\frac{1}{E_2}\left[\left(\frac{1}{{\alpha }_2}-\frac{1}{{\alpha }_3}\right)\frac{1}{h_2}y+\frac{1}{{\alpha }_3}\right],0<y<h_2\\ \frac{1}{E_i}\frac{1}{{\alpha }_2}, h_2<y<h_2+h_j\\ \frac{1}{E_1}\left[\left(\frac{1}{{\alpha }_1}-\frac{1}{{\alpha }_2}\right)\frac{1}{h_1}\left(y-h_2-h_j\right)+\frac{1}{{\alpha }_2}\right],h_2+h_j<y<h_1+h_2+h_j\end{array}\right.$$

(15)

From the relationship between the displacement and strain,

$$\Delta H={{\displaystyle \int }}_0^{h_1+h_2+h_j}\epsilon \left(y\right)dy$$

(16)

Substituting Equation (15) into Equation (16), we obtain

$$\begin{gathered} \Delta H={\left.\frac{F}{A{E}_2}\left[\frac{1}{2}\left(\frac{1}{{\alpha }_2}-\frac{1}{{\alpha }_3}\right)\frac{1}{h_2}{y}^2+\frac{1}{{\alpha }_3}y\right]\right|}_0^{h_2}+{\left.\frac{F}{A{E}_j}\frac{1}{{\alpha }_2}y\right|}_{h_2}^{h_2+h_j}+{\left.\frac{F}{A{E}_1}\left[\frac{1}{2}\left(\frac{1}{{\alpha }_1}-\frac{1}{{\alpha }_2}\right)\frac{1}{h_1}\left(y^2-2h_{2}y-2h_{j}y\right)+\frac{1}{{\alpha }_2}y\right]\right|}_{h_2+h_j}^{h_1+h_2+h_j}= \\ \frac{F{h}_2}{A{E}_2}\left(\frac{1}{2{\alpha }_2}+\frac{1}{2{\alpha }_3}\right)+\frac{F{h}_j}{A{E}_j}\frac{1}{{\alpha }_2}+\frac{F{h}_1}{A{E}_1}\left(\frac{1}{2{\alpha }_1}+\frac{1}{2{\alpha }_2}\right)=\frac{F}{A}\left[\frac{h_2\left({\alpha }_2+{\alpha }_3\right)}{2E{\alpha }_2{\alpha }_3}+\frac{h_j}{E_j{\alpha }_2}+\frac{h_1\left({\alpha }_1+{\alpha }_2\right)}{2E_1{\alpha }_1{\alpha }_1}\right] \end{gathered}$$

(17)

Substituting Equation (17) into Equation (1), we obtain

$$K=\frac{2A{E}_1{E}_2{E}_j{\alpha }_1{\alpha }_2{\alpha }_3}{h_2\left({\alpha }_2+{\alpha }_3\right)E_1{E}_j{\alpha }_1+2h_j{E}_1{E}_2{\alpha }_1{\alpha }_3+h_1\left({\alpha }_1+{\alpha }_2\right)E_2{E}_j{\alpha }_3}$$

(18)

④ Parallel arrangement with different contact areas

In parallel arrangement, if the total area of the upper slates is A, the ratio of the effective contact areas of each contact surface to the total area of the upper slates are ${\alpha }_1$, ${\alpha }_{21}$, ${\alpha }_{22}$, ${\alpha }_{31}$, ${\alpha }_{32}$, and the moduli of elasticity of the sized slates and joints are $E_1$, $E_2$, $E_3$, $E_j$, and the ratio of the height of slate 3 to that of slate 2 is ${\beta }_h$. For the convenience of calculation, it is still assumed that the internal stress influence range varies linearly and the stress is uniformly distributed within the influence range, whose calculation sketch is shown in Figure 8, where the contact area of the lower surface of the upper slate is equal to the sum of the contact area of the upper surface of the two stones below.

According to the method described in the previous section, the deformation of the upper slates and joints can be obtained as follows, respectively:

$$\Delta H_1={\left.\frac{F}{A{E}_1}\left[\frac{1}{2}\left(\frac{1}{{\alpha }_1}-\frac{1}{{\alpha }_{21}+{\alpha }_{31}}\right)\frac{1}{h_1}\left(y^2-2h_{2}y-2h_{j}y\right)+\frac{1}{{\alpha }_{21}+{\alpha }_{31}}y\right]\right|}_{h_2+h_j}^{h_1+h_2+h_j}=\frac{F{h}_1}{A{E}_1}\left(\frac{1}{2{\alpha }_1}+\frac{1}{2{\alpha }_{21}+2{\alpha }_{31}}\right)$$

(19)

$$\Delta H_j=\frac{F{h}_j}{A{E}_j}\frac{1}{{\alpha }_{21}+{\alpha }_{31}}$$

(20)

Since the deformations of the two stones below are equal when arranged in parallel, it follows that

$$F=K_2\Delta H_2+K_3\Delta H_3=\left(K_2+K_3\right)\Delta H_2$$

(21)

The displacement of lower slate 2 is

$$\Delta H_2={\left.\frac{F_2}{A{E}_2}\left[\frac{1}{2}\left(\frac{1}{{\alpha }_{21}}-\frac{1}{{\alpha }_{22}}\right)\frac{1}{h_2}{y}^2+\frac{1}{{\alpha }_{23}}y\right]\right|}_0^{h_2}=\frac{F_2{h}_2}{A{E}_2}\left(\frac{1}{2{\alpha }_{21}}+\frac{1}{2{\alpha }_{22}}\right)$$

(22)

From the definition of stiffness, it follows that

$$K_2=\frac{F_2}{\Delta H_2}=\frac{2A{E}_2{\alpha }_{21}{\alpha }_{22}}{h_2\left({\alpha }_{21}+{\alpha }_{22}\right)}$$

(23)

Similarly, the stiffness of lower slate 3 can be obtained as

$$K_3=\frac{F_2}{\Delta H_3}=\frac{2A{E}_3{\alpha }_{31}{\alpha }_{32}}{{\beta }_h{h}_2\left({\alpha }_{31}+{\alpha }_{32}\right)}$$

(24)

Substituting Equation (23) and Equation (24) into Equation (21), we obtain

$$\Delta H_2=\frac{F}{K_2+K_2}=\frac{F}{A}\frac{1}{\frac{2E_2{\alpha }_{21}{\alpha }_{22}}{h_2\left({\alpha }_{21}+{\alpha }_{22}\right)}+\frac{2E_3{\alpha }_{31}{\alpha }_{32}}{{\beta }_h{h}_2\left({\alpha }_{31}+{\alpha }_{32}\right)}}$$

(25)

Associating Equations (19), (20), and (25), from the definition of equivalent stiffness, it follows that

$$K=\frac{F}{\Delta H_1+\Delta H_2+\Delta H_j}$$

(26)

The arrangement of all the slates in Yuanjue Cave can be reduced to one or more series and parallel combinations of ①, ②, ③, and ④. According to the abovementioned method, the equivalent support stiffness of all the rocks in Yuanjue Cave can be obtained.

(2)

Determination of the contact area of the sized slates

The contact area of the upper and lower surfaces of the strip is calculated according to the contact between the strips (or the top plate), combined with parameters such as the length and depth of the strips. The calculation method is as follows:

① If the upper and lower surfaces of the slate are “tight” or “tightly closed”, the upper and lower surfaces of the slate are considered to be in close contact with the slate (or top plate), and the contact area of the upper and lower sides is calculated accordingly.

$$S_{ts}=l·d$$

(27)

where $l$ is the contact length and $d$ is the contact depth. Among them, when determining the contact length and contact depth, in addition to considering the spatial position relationship between the strips, the influence of the gap between the left and right strips and the overhanging depth of the strips should also be considered.

② If there is a gap between the upper and lower surfaces of the slate, and there is no filling, the contact area between the upper and lower parts of the slate is 0.

③ If there is a gap between the upper and lower surfaces of the strip, but there is a certain degree of filling, it will be regarded as a joint with a certain bearing capacity according to the specific filling, and its contact area is the same as that between the adjacent strip (or top plate).

(3)

Slate and joint parameter determination

According to practical experience in rock engineering, the equivalent stiffness calculation requires the following principles for the selection of parameters related to slates and joints:

① Selecting the stratigraphic parameters similar to the lithology of the slates as the parameters of the slates (Table 1);

② When the joint between the slates (or the top slab) is tightly closed, the joint height is taken as 0;

③ When the gap between the boulders (or the top slab) is filled, the elastic mode of the joint is 1/5 of the elastic mode of the lower boulder;

④ When there is a gap between the slate and the stone (or the top plate) and the filling is small, the jointed elastic mold is 1/50 of the lower slate’s elastic mold;

⑤ When there is a gap between the slate and the stone (or the top plate) and there is no filling, the two do not come into contact, and the jointed elastic mode is 0;

According to the physical and mechanical indexes of the strips, the strips are divided into five groups.

Table 1. Slate-related parameters.

Grouping	Elastic Modulus (Pa)	Poisson Ratio	Cohesion (MPa)	The Angle of Internal Friction (°)	Natural Density (kg/m3)	Tensile Strength (MPa)
Group TS04	2.76 × 1010	0.19	1.38	38.1	2230	0.2
Group TS05	1.38 × 1010	0.25	1.02	28.2	2150	0.2
Group TS09	2.50 × 1010	0.21	1.1	30.6	2200	0.1
Group TS56	1.38 × 1010	0.25	1.02	28.2	2150	0.1
Group TS01	2.50 × 1010	0.21	1.1	30.6	2200	0.2

Table 1. Slate-related parameters.

Grouping	Elastic Modulus (Pa)	Poisson Ratio	Cohesion (MPa)	The Angle of Internal Friction (°)	Natural Density (kg/m3)	Tensile Strength (MPa)
Group TS04	2.76 × 1010	0.19	1.38	38.1	2230	0.2
Group TS05	1.38 × 1010	0.25	1.02	28.2	2150	0.2
Group TS09	2.50 × 1010	0.21	1.1	30.6	2200	0.1
Group TS56	1.38 × 1010	0.25	1.02	28.2	2150	0.1
Group TS01	2.50 × 1010	0.21	1.1	30.6	2200	0.2

(4)

Equivalent support stiffness calculation

From the survey results, it can be seen that the strips TS04, TS11, TS15, TS20, TS55, TS56, and TS57 all have a certain gap with the top plate of Yuanjue Cave, with a minimum of 0.5–1 cm and a maximum of 3–5 cm, and are not filled. Therefore, the overhang of the top slab occurs at the abovementioned slate positions, and the equivalent stiffness of each corresponding group of overhanging slates is 0. Based on this, the equivalent stiffness of each group of slates can be obtained by combining the series and parallel laws in different modes (

Table 2. Equivalent support stiffness of cracked boulders.

Grouping	K(N/m)	Grouping	K(N/m)
TS01	3.42 × 1010	TS38	5.22 × 109
TS04	0	TS39	6.94 × 109
TS10	1.55 × 1010	TS40	5.16 × 109
TS11	0	TS41	5.69 × 109
TS12	3.56 × 1010	TS42	1.19 × 109
TS13	9.38 × 109	TS43	9.10 × 109
TS14	6.00 × 109	TS46	1.01 × 107
TS15	0	TS49	1.12 × 1010
TS16	9.99 × 109	TS51	1.81 × 1010
TS20	0	TS52	2.14 × 1010
TS25	8.51 × 109	TS55	0
TS29	2.68 × 1010	TS56	0
TS30	6.86 × 109	TS57	0
TS34	5.00 × 109	TS59	3.07 × 1010

).

4. Numerical Simulation and Results

The parameters of rock mechanics required for the numerical simulation were obtained after the aforementioned evolutionary law and geomechanical qualitative analysis (

Table 3. Stratigraphic lithological parameters.

Sequence Stratigraphy	Formation Lithology	Bulk Modulus (GPa)	Shear Modulus (GPa)	Cohesion (MPa)	Friction Angle (°)	Density (kg/m3)	Tensile Strength (MPa)
1	Fine sandstone	10	5	1.38	38.14	2230	0.2
2	Slate	10	7	1.5	30	2300	0.2

). The boundary conditions of the rock body of the roof slab of Yuanjue Cave are very clear, and this chapter analyzes them from the perspective of 3D modeling and numerical simulation for a precise quantitative perspective [46]. Combined with the aforementioned stratigraphic distribution characteristics and physical and mechanical parameters of the rock body [47,48], the final stability study was conducted by finite difference software for the effect of overhanging, crushing, and damage of the roof slab on the overall top slab of Yuanjue Cave, respectively. In the aforementioned study, it is given that the top plate of the Yuanjue Cave has formed a relatively independent block structure in the form of a “triangular thick plate” (Figure 3), in which the J10 fissure forms a natural cutting boundary on the east side, so the model boundary on the east side can be extended 5 m eastward as the boundary of the J10 fissure. The rock body of the top slab on the west side forms a fracture boundary along the L6 fissure, but this boundary intersects the axis of the cave chamber horizontally at a large angle to accommodate the space of the cave chamber, so the west side extends westward to the end of the cliff wall of the L6 fissure with the west boundary of the cave chamber. The northern boundary is the hollow surface of the cliff wall, and the entrance of the cave canal can be used as the boundary; the southern side is bounded by the convergence angle of the roof slab. The rock supporting the top slab is the image formation, which is relatively stable and intact, so it can be extended below the top slab to the bottom part of Yuanjue Cave, and the three-dimensional model is shown in Figure 9. The model of the top slab of the Yuanjue Cave was established with 109,524 units and 19,856 nodes.

During the numerical simulation, the vertical displacement and the maximum principal stress location distribution of the roof slab were monitored at the four edge measurement points above the roof slab of Yuanjue Cave. By comparing and analyzing the vertical displacement and the maximum principal stress development trend of the roof slab at the four measurement points, the effect of the number of overhanging, broken, and damaged boulders on the force deformation characteristics and stability of the roof slab of Yuanjue Cave was revealed. In addition, since the L6 crack developed in the roof of the Yuanjue Cave, it was formed by the gradual cracking of the filling slate after the excavation of Yuanjue Cave and the increase of the bending moment of the roof slab, which can reflect the force deformation and stable state of the roof slab. Therefore, in the process of numerical simulation, the L6 crack width was monitored as a focus, and the changing trend of the L6 crack width was mainly analyzed.

4.1. Analysis of the Effect of Slate Overhang on the Force Deformation of the Roof Slab

The study on the effect of overhanging slates showed an increase in the number of overhanging slates sequentially from the drainage of J10 behind the cavern of Yuanjue Cave to the direction of the light hole under the current condition. The empty model in FLAC 3D numerical analysis software was used to simulate the state of the overhanging slates without bearing pressure, and the change in the number of overhanging stones was simulated by changing the intrinsic structure of the slates (He et al. [15]).

Figure 10 shows, from left to right, the vertical displacement clouds of the top slab of Yuanjue Cave in the case of 1, 4, 10, and 15 overhanging slates. It can be seen that with the increase in the number of overhanging boulders, the top slab of Yuanjue Cave underwent uneven settlement, and the maximum settlement position was shifted from the central area to the left rear of Yuanjue Cave (boulder overhanging side).

According to the numerical simulation results, the variation curve of vertical displacement of each measurement point with the number of overhanging slate blocks was made with the number of overhanging slate blocks as the horizontal coordinate (Figure 11a). The horizontal displacement in the x and y directions of the two pairs of measurement points above and below the L6 crack was decomposed and differentiated along the vertical direction of the L6 crack to obtain the variation curve of the width above and below the L6 crack with the number of overhanging slate blocks (Figure 11b).

With the increase of the number of overhanging slates, the settlement of Y6 and Y10 measuring points on the roof of the Yuanjue Cave decreased, the settlement of Y8 and Y9 measuring points increased, the width of L6 cracks kept increasing, and the change rate increased with the increase of the number of overhanging boulder blocks. The settlement of the Y8 measuring point was the largest and the most significant change (Figure 11a), and the settlement of Y6, Y9, and Y10 were smaller because the Y8 measuring point was close to the development area of overhanging slates. The settlement of Y8 and Y9 began to increase sharply when the overhanging boulders exceeded 12 (TS30) and entered the stage of rapid change.

The L6 crack pattern was “wide at the top and narrow at the bottom”, which is a typical tensile crack, and the L6 crack width increased with the increase of overhanging slates, and when the overhanging slates exceeded 12 pieces, the L6 crack width started to increase sharply and entered the rapid change stage (Figure 11b). The development process of the settlement of the roof slab and the L6 crack width of Yuanjue Cave was similar, and both can be divided into three stages, i.e., the initial uniform change stage (the number of overhanging sized slates is less than 3), the middle stable change stage (the number of overhanging sized slates is between 3 and 12), and the later rapid change stage (the number of overhanging sized slates exceeds 12). When the number of overhanging slates exceeds 12, the vertical displacement and L6 crack width of each measurement point enter the late rapid change stage. Based on the principle of preventive protection, the number of overhanging slates reaching 12 can be taken as a sign that the roof of Yuanjue Cave enters the early warning state of instability. The critical number of overhanging stones was 12 (TS30 location in Figure 12).

4.2. Analysis of the Effect of Slate Crushing on the Force Deformation of the Roof Slab

The analysis scheme was used to increase the number of broken slates in sequence from the light hole to the direction of the J10 drainage place behind Yuanjue Cave and simulated the change of the number of broken- slates by changing the slate structure relationship. Figure 13 shows the vertical displacement of the roof slab of Yuanjue Cave under the action of different crushed rocks. With the increase in the number of crushed rocks, the roof slab of Yuanjue Cave underwent uneven settlement, and the maximum settlement position was shifted from the central area to the left front of Yuanjue Cave (the side of crushed rocks), and continuously developed to the left rear side. The above analysis shows that the number of broken boulders has a significant influence on the force deformation and stability of the roof slab of Yuanjue Cave.

With the increase in the number of broken slate blocks, the settlement of Y6, Y8, and Y9 on the top plate of Yuanjue Cave gradually increased, the settlement of Y10 decreased, and the width of L6 cracks kept increasing. The settlement of the four monitored ones had the characteristics of Y6 > Y8 > Y9 > Y10 (Figure 14a), which is mainly since the monitoring points of Y6, Y9, and Y8 were closer to the broken slate area; when the broken slate exceeded 18 pieces (TS46), the settlement of Y6, Y8, and Y9 started to increase sharply and entered the stage of rapid change.

The L6 crack still showed the characteristic of being “wide at the top and narrow at the bottom”, and the width of the L6 crack increased with the increase of broken boulders, and when the number of broken boulders exceeded 18 (TS46), the width of the L6 crack stared to increase sharply and entered the rapid change stage (Figure 14b). According to the evolution law of vertical displacement and the L6 width change rate of measurement points, the change process can be divided into three stages: the initial uniform change stage (less than 11 pieces of broken slates), the middle stable change stage (the number of broken slates is between 11 and 18 pieces), and the late rapid change stage (broken slates have more than 18 pieces).

When the broken slate exceeded 18 pieces, the vertical displacement and L6 crack width of each measurement point entered the late rapid change stage. Based on the principle of preventive protection, the broken slate reaching 18 pieces can be taken as a sign that the top plate of the Yuanjue Cave has entered the destabilization warning state, i.e., the critical number of the broken slate is 18 pieces (Figure 15).

4.3. Analysis of the Influence of the Degree of Damage to the Cracked Slates on the Force Deformation of the Top Slab

Through the preliminary cave investigation, it was found that the single blocks of cracked boulders included TS01, TS02, TS03, TS04, TS07, TS08, TS10, TS11, TS12, TS14, TS17, TS31, TS32, TS36, TS52, and TS53. According to the contact area and grouping of the boulders, the groups where the cracked boulders are located were TS01, TS10, TS12, TS14, TS16, TS29, TS30, TS34, TS49, TS51, TS52, as shown in Figure 16. The cracked and damaged boulders were mainly distributed in the middle and front parts of the boulder-filled area.

To analyze the influence of different degrees of damage on the stability of the roof plate of Yuanjue Cave, the mechanical parameters such as the modulus of elasticity of the cracked damaged slates were reduced to 1/2, 1/5, 1/10, 1/50, and 1 × 10⁻⁶ of the original ones (where 1 × 10⁻⁶ can be equated to complete damage). The overall equivalent stiffness of each group of slates was calculated, and the results are shown in

Table 4. Equivalent bearing stiffness of the cracked boulders.

Damage Factor	0.5	0.8	0.9	0.98	1
Grouping	K(N/m)	K(N/m)	K(N/m)	K(N/m)	K(N/m)
TS01	1.71 × 1010	6.84 × 109	3.42 × 109	6.84 × 108	≈0
TS10	7.77 × 109	3.11 × 109	1.55 × 109	3.11 × 108	≈0
TS12	1.78 × 1010	7.12 × 109	3.56 × 109	7.12 × 108	≈0
TS14	4.71 × 109	2.87 × 109	1.73 × 109	4.17 × 108	≈0
TS16	8.91 × 109	6.74 × 109	4.79 × 109	1.45 × 109	≈0
TS29	1.89 × 1010	1.09 × 1010	6.66 × 109	1.64 × 109	≈0
TS30	5.23 × 109	3.30 × 109	2.11 × 109	5.52 × 108	≈0
TS34	4.00 × 109	2.50 × 109	1.54 × 109	3.78 × 108	≈0
TS49	1.04 × 1010	9.55 × 109	9.13 × 109	8.73 × 109	8.61 × 109
TS51	1.21 × 1010	6.05 × 109	3.30 × 109	7.13 × 108	≈0
TS52	1.28 × 1010	6.43 × 109	3.66 × 109	8.48 × 108	≈0

(only the slates with changed equivalent stiffness are listed).

Figure 17 shows the vertical displacement of the roof slab of the cave from left to right when the damage factor was 0.5, 0.8, 0.9, and complete damage. With the increase of the damage of the cracked slate, the uneven settlement of the roof plate of the cave occurred, and the maximum settlement position was shifted from the central area to the left front of the cave, which had a significant effect on the force deformation and stability of the roof plate of the cave. According to the numerical simulation results, the vertical displacement curve of each measurement point with the damage factor was made in the horizontal coordinate of the cracked stone damage factor (Figure 18a). The horizontal displacement in the x and y directions of the two pairs of measurement points above and below the L6 crack was decomposed and differentiated along the vertical direction of the L6 crack. The change curve of the width of the L6 crack above and below the damage factor was obtained (Figure 18b).

The settlement of Y6, Y8, and Y9 gradually increased and the settlement of Y10 decreased with the increase of damage factor, and the change rate increased with the increase of the damage factor of L6. Y9 settlement started to increase sharply and entered the stage of rapid change. At the same time, the L6 cracks still showed the typical tensile fracture characteristics of “wide at the top and narrow at the bottom”, and the width of the cracks kept increasing as the damage of the cracked rocks intensified. When the damage factor of the cracked slate exceeded 0.9, the L6 crack width started to increase sharply and entered the stage of rapid change.

According to the evolution law of vertical displacement and the L6 width change rate of measurement points, the change process can be divided into three stages: the initial uniform change stage (damage factor less than 0.8), the middle accelerated change stage (damage factor between 0.8 and 0.9), and the late rapid change stage (damage factor greater than 0.9). When the damage factor exceeded 0.9, the vertical displacement of each measurement point and L6 crack width entered the late rapid change stage; based on the principle of preventive protection, the damage factor could reach 0.9 blocks, serving as a sign that the top plate of the Yuanjue Cave is entering the early warning state of instability, that is, the critical value of the damage factor is 0.9.

4.4. Reinforcement Measures

Based on the analysis of the stability of the roof of Yuanjue Cave and the calculation results of the support stress of the lower stone, combined with the requirements of the space in the cave for the supporting technology, the following two reasonable support schemes in the cave are proposed:

(1) Slate replacement and gap filling: firstly, we analyzed the influence of the reinforcement measures such as the size of slate replacement and gap filling on the stability of the roof slab and the vertical stress of the stone in Yuanjue Cave under the current conditions. A more reasonable reinforcement plan of stone replacement and gap filling to improve the stability of Yuanjue Cave was determined. Reinforcement measures such as replacing these developmentally damaged size slates with intact bars and filling gaps were used to improve the overall stability of the roof slab of Yuanjue Cave;

(2) Set the support column: firstly, the position where the roof support stress is the largest or close to the bearing capacity was analyzed, and the reasonable support position in the cave was determined according to the position requirement of the support member in the cave space (Figure 19). According to the previous survey results, the severely damaged stones are mainly distributed in the front and middle of the left side of the entrance of Yuanjue Cave, and the maximum displacement of the damaged slates in the numerical simulation results is also in the area. Therefore, it is proposed to adopt a support column with strong vertical stability as the support scheme. A circular support column should be set between the first and second Buddha statues on the right side of the cave entrance.

5. Conclusions

In this paper, the equivalent support stiffness was analytically calculated for three cases of overhanging, broken, and cracked slates, respectively, and the analytical solutions were obtained. Numerical simulations were used to analyze the stress field and displacement field of the roof slab of Yuanjue Cave. To reveal the stress and deformation characteristics of the roof slab of Yuanjue Cave, as well as the main factors affecting the stability of the roof slab and the slates that need to be monitored as a focus, the conclusions obtained are summarized as follows.

(1) The influence law of overhanging, broken stone block number change, and cracked stone damage change on the force deformation and stability of the top plate of Yuanjue Cave can be roughly divided into three stages, including the initial stage of uniform velocity change, the middle stage of stable change or the middle stage of accelerated change, and the late stage of rapid change. Therefore, the rapid change stage can be taken as a sign that the roof plate of Yuanjue Cave is entering the destabilization warning state; the critical number of overhanging slates is 12, the critical number of broken slates is 18, and the critical value of the damage factor is 0.9;

(2) An increase in the number of overhanging and broken slates and the increase in the damage factor of cracked slates would aggravate the degree of uneven force deformation of the roof slab. With an increase in the number of overhanging slates, the maximum settlement position of the top slab was shifted from the central area to the left rear (slate-overhanging side). With an increase in the number of crushed slates, the maximum settlement position of the top slab was shifted from the central area to the left front (bar-crushing side) and kept developing to the left rear side. With an increase in the damage factor of the currently cracked slates, the maximum settlement position of the top slab was shifted from the central area to the left front of the hole;

(3) With an increase in the number of overhanging and broken boulders, the vertical stresses of the boulders with extreme values of stress increased continuously and changed significantly, while the changes of other boulders were relatively small. With an increase of the damage factor of the cracked boulders, their equivalent support stiffness decreased, and the vertical stresses of TS12, TS16, and TS34 decreased continuously, while the vertical stresses of TS13, TS25, TS43, and TS49 increased continuously.