*Article* **Study on Asymmetric Failure and Control Measures of Lining in Deep Large Section Chamber**

**Yubing Huang 1,\*, Bei Jiang 2,\*, Yukun Ma 1, Huayong Wei 2, Jincheng Zang <sup>3</sup> and Xiang Gao <sup>3</sup>**


**Abstract:** Lining is often used as the last line of defense in deep large section chamber. Under the asymmetric load, it is easy to damage, resulting in the overall repair of the chamber. Aiming at this problem, taking the pump house in Wanfu Coal Mine under construction in China as an engineering example, we analyzed the asymmetric failure of pump house lining caused by construction disturbance, established the lining mechanical model and quantitative evaluation indexes, such as bending moment change rate, bending moment balance rate, displacement change rate and displacement balance rate, studied the influence mechanism of asymmetrical coefficient, section size and lining thickness on the mechanical behavior of lining, and proposed the control measures of deep large section chamber with the core of "strengthening asymmetric support, reducing section size and improving lining strength". The field monitoring shows that the asymmetric deformation of the pump house is effectively controlled, and the maximum displacement is only 7.3 mm, which ensures the long-term stability of the chamber.

**Keywords:** deep chamber; asymmetric failure; mechanical analysis; control measures

### **1. Introduction**

With the depletion of shallow coal resources and the improvement of coal mining modernization levels, the mine is developing onto a large scale, and there are more and more large section chambers [1,2]. At present, the combined support method of bolt mesh spray and high strength lining is mostly used in the deep chamber in the large section [3–8]. Especially in the core chamber with long service life, high-strength lining is regarded as the last line of defense for its advantages of high strength, high stiffness, strong sealing and favorable appearance [9,10]. When confronted with the high stress, inclined rock, secondary disturbance and other adverse factors, it is difficult to effectively control the stability of surrounding rock with a bolt mesh spray. As the last line of defense, the lining bears the load of surrounding rock; once damaged, the chamber will face a comprehensive repair. Under complex geological conditions, the loads of surrounding rock are often uneven and asymmetric, and the bearing capacity of the lining is far lower than the theoretical value. Consequently, the phenomenon of local cracking, spalling and even overall instability is very common, which poses a great threat to the safety of mine production [11–13].

Aiming at the deformation and failure of lining under asymmetric load, scholars have carried out a large number of studies [14–24]. Krzysztof Skrzypkowski [25] studied the influence of geological and mining factors on the stability of chamber excavation through roadway anchor support, and determined the key position of chamber failure. Xu Ying, Jing Shenguo et al. [26] established a mechanical model of chamber lining-anchor cable

**Citation:** Huang, Y.; Jiang, B.; Ma, Y.; Wei, H.; Zang, J.; Gao, X. Study on Asymmetric Failure and Control Measures of Lining in Deep Large Section Chamber. *Energies* **2021**, *14*, 4075. https://doi.org/10.3390/ en14144075

Academic Editor: Frede Blaabjerg

Received: 9 June 2021 Accepted: 2 July 2021 Published: 6 July 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

coupling support for a deep large-section pump house, studied the influence law of support parameters on the lining ultimate bearing capacity under non-uniform load, and proposed an effective reinforcement scheme. Wang Shuhong et al. [27] carried out an experimental study of the influence of the position of the cavity behind the wall on the lining failure mode, and summarized the failure modes and failure sequences of different cavity tunnels. Sun Xiaoming, Zhang guofeng, et al. [28] studied the asymmetric deformation and failure of roadway in inclined strata, obtained the key parts of the asymmetric deformation and failure of roadway in deep inclined strata, and proposed the asymmetric coupling control measures. Wang Jiong, Guo Zhibiao, et al. [29] studied the asymmetric deformation and failure characteristics of deep roadway crossings by means of numerical analysis and engineering tests, obtained the key parts of the first failure and put forward the asymmetric coupling support measures of anchor net cables. The above studies mainly focus on the asymmetric failure of lining caused by the cavity behind the lining and the attitude of the surrounding rock. The excavation of adjacent chambers will cause the construction disturbance, bringing the asymmetric load to the stable lining, which leads to deformation and failure, especially for the deep chamber with large section. At present, there are few studies on the mechanical behavior and control method of deep large section chamber lining caused by construction disturbance.

Based on this, this paper takes the permanent pump house of Wanfu Coal Mine in China as the engineering background, establishes the lining mechanical model and quantitative evaluation index based on an asymmetric coefficient, studies the influence mechanism of asymmetric coefficient, section size and lining thickness on lining mechanical behavior, proposes the control method of a deep large section chamber and guides the reinforcement of a permanent pump house in the Wanfu Coal Mine. It has an important guiding significance for engineering under similar conditions.

#### **2. Engineering Background**

#### *2.1. Engineering and Geology Background*

Wanfu Coal Mine is located at the southernmost point of Juye Coalfield in eastern China. It is a large-scale modern mine under construction, with a designed production capacity of 1.8 Mt/year and a design service life of 67 years. The shaft station is located at the level of −820 m, and the pump house is the largest chamber at this depth. The maximum height of the excavated section of the pump house is 10.4 m, and the maximum width is 7.2 m. The water distribution drift is located on the east side of the pump house, and is the nearest exploitation field to the pump house. The nearest horizontal distance is only 2.5 m, and the nearest vertical distance is only 1.8 m. The mine location and shaft station plan are shown in Figure 1.

⌥ᡯ᭟ᣔо⍻⛩⽪മ

**Figure 1.** Mine location and −820 m shaft station plan.

The pump house is located in the inclined rock strata with the interaction of mud and sand, and the dip angle of the rock strata is 13◦. From the bottom to the top, the lithology is mudstone, siltstone and fine sandstone. The maximum horizontal stress in the area is 37.1 MPa, and the vertical stress is 23.8 MPa. The pump house adopts the support form of anchor net spray + reinforced concrete lining. The parameter of anchor cable is Φ22 × 6200 mm, the inter-row spacing is 1600 × 1600 mm. Bolt parameters are Φ22 × 2500 mm, the inter-row spacing is 800 × 800 mm. Reinforced concrete lining is C40 concrete, 500 mm thick, double reinforced. The support and rock distribution of the pumping house are shown in Figure 2.


**Figure 2.** The support and rock distribution of the pumping house.

#### *2.2. Field Failure Analysis*

After lining work, long-term asymmetric deformation monitoring was carried out, including convergence of two walls (OA and OB), settlement of two arch shoulders (MC and ND) and displacement of roof and floor (OE and OF). The sensor used for convergence measurement is a mechanical convergent ruler, the accuracy is up to 0.1 mm and the monitoring frequency is once every two days by manual. The layout of measuring points is shown in Figure 2, the displacement curve and on-site damage are shown in Figure 3.

The deformation process of lining can be divided into three stages: <sup>1</sup> the initial stable stage, as of 328 days after lining working, the pump house is stable as a whole. <sup>2</sup> The disturbance stage, during the construction of the inside sump and water distribution drift, the deformation of the pumping house increased rapidly, and the deformation increment is as follows: west straight wall > west arch shoulder > arch top > east straight wall > arch bottom > east arch shoulder. The overall law is that the deformation on the west side is greater than that of the east side. <sup>3</sup> In the slower growth stage, the deformation development ratio slows down, but it is not completely stable.

In conclusion, due to the influence of water distribution drift and inside sump excavation, the original stable lining of the pump house shows asymmetric failure characteristics, and the deformation of the lining near the disturbed side is far greater than that of the other side. After the lining is damaged, the radial restraint effect on the surrounding rock is lost, and the deformation continues to develop slowly. If the reasonable reinforcement measures are not taken, it is difficult to meet the requirements.

**Figure 3.** The displacement curve and on-site damage.

#### **3. Asymmetric Mechanical Model and Evaluation Index**

In this section, the asymmetric mechanical model and quantitative evaluation index of lining are established to provide a basis for analyzing the deformation and stress rule of the lining under an asymmetric load.

#### *3.1. Asymmetric Mechanical Model*

Along the axial direction of the chamber, the lining per unit length is abstracted as a mechanical model, as shown in Figure 4a. In this mechanical model, the arch crown radius is r1, the arch bottom radius is r2, and the straight wall height is h. The asymmetrical coefficient of lining stress is λ. The uniform load on the left arch crown is q1, the uniform load on the left wall is q2, and the uniform load on the left bottom arch is q3. The uniform load on the right arch crown, straight wall and bottom arch is λq1, λq2 and λq3, respectively.

Based on the principle of structural mechanics [30], mechanic model of lining under asymmetric conditions is divided into positive symmetry and anti-symmetry, and the "force method" is used to solve them, respectively.

Half of the positive symmetric part is taken for analysis, as shown in Figure 4b, and it is simplified into a quadratic overstatically indetermined structure. Releasing the support constraint at point A of the vault, and replacing it with the redundant unknown forces *X*<sup>1</sup> and *X*<sup>2</sup> and *X*<sup>1</sup> represents the horizontal support reaction (Fax) at point A, and *X*<sup>2</sup> represents the constraint bending moment (*Ma*) at point A, as shown in Figure 4b.

(**a**) mechanical model (**b**) Positive symmetry calculation model (**c**) Antisymmetric calculation model

Half of the antisymmetric part is taken for analysis, as shown in Figure 4c, and it is simplified to a statically indeterminate structure. Releasing the support constraint at point A of the vault, and replacing it with the redundant unknown force *X*3, which represents the vertical support reaction (Fay) at point A, as shown in Figure 4c.

$$\begin{cases} \delta\_{11}X\_1 + \delta\_{12}X\_2 + \Delta\_{1P} = 0\\ \delta\_{21}X\_1 + \delta\_{22}X\_2 + \Delta\_{2P} = 0\\ \delta X\_3 + \Delta = 0 \end{cases} \tag{1}$$

According to Equations (2) and (3), the basic parameters in (1) can be obtained.

$$\begin{array}{rcl} \mathcal{S}\_{11} &=& \int \frac{\overline{M\_1 M\_2}}{EI} d\_s \\ &=& \frac{9\pi R\_1^3 - 24R\_1^2 + 12R\_1 h + 4h^2}{12EI} + \\ &\frac{4R\_2 a \left(R\_1 h\right)^2 + R\_2^3 \left(2\sin 2a \left(1 - \cos a\right) - 2\sin a \left(1 - \cos 2a\right) + 6a\right)}{4EI} \end{array} \tag{2}$$

Similarly, *δ*12, *δ*21, *δ*<sup>22</sup> and *δ* can be obtained.

$$\begin{array}{rcl} \Delta\_{2P} &=& \int \frac{\overline{M\_2} M\_P}{EI} d\_s \\ &=& \frac{q\_1 \mathcal{R}\_1 [R\_1^2 (3\pi - 6) + 6h + 3h^2 + 6\mathcal{R}\_2 \mathcal{R}\_1 a + 6h a] + q\_2 h^2 (3a + h)}{6EI} + \\ & \frac{\mathcal{R}\_2^2 [q\_3 \mathcal{R}\_2 (a - \sin a) + (q\_1 \mathcal{R}\_1 + q\_2 h)(\sin a - \cos a) - q\_1 \mathcal{R}\_1 (a \sin a + \cos a - 1)]}{EI} \end{array} \tag{3}$$

Similarly, Δ1*<sup>P</sup>* and Δ can be obtained.

For lining structure, the displacement is mainly caused by bending moment, and the displacement caused by axial force and shear force is small. Therefore, when analyzing the mechanical behavior of lining influenced by asymmetric load, only bending stress is considered. According to Equation (4), the bending moment value of any section can be obtained.

$$\begin{array}{l} M\_{\mathfrak{a}} = M\_{\mathfrak{p}\mathfrak{a}} + X\_1 \overline{M\_1} + X\_2 \overline{M\_2} \\ M\_{\mathfrak{b}} = M\_{\mathfrak{p}\mathfrak{b}} + X\_3 \overline{M} \\ M = M\_{\mathfrak{a}} + M\_{\mathfrak{b}} \end{array} \tag{4}$$

where (1)–(4), *δ* is the angle from the top of the structure to any section. *δij* represents the displacement of the basic structure along the *Xi* direction under the single action of the unit force, Δ represents the displacement of the basic structure along the *X* direction under the action of the load alone, EI is bending stiffness, which refers to the ability of an object to resist its bending deformation. *M*<sup>1</sup> and *M*<sup>2</sup> represent the bending moment of the unit force in any section of the basic structure, and *Mi* represents the bending moment of the load in any section of the basic structure [31].

#### *3.2. Analysis Scheme*

Mechanical behavior of lining under asymmetric load is most significantly affected by asymmetric coefficient, section size and thickness [32]. Based on the actual parameters of the pump house in Wanfu Coal Mine, different asymmetric coefficients, section size and lining thickness are selected to analyze the influence rule on lining bending moment and displacement, corresponding to the schemes A*i*, B*<sup>i</sup>* and C*i*, respectively, where *i* = 1~7. The scheme design is shown in Tables 1–3.


**Table 1.** Analysis scheme of different asymmetric coefficients.

#### **Table 2.** Analysis scheme of different chamber section size.


**Table 3.** Analysis scheme of different lining thickness.


### *3.3. Evaluation Indexes*

In order to quantitatively analyze the influence rules of different factors on lining bending moment and displacement, evaluation indexes such as bending moment change ratio, bending moment balance ratio, displacement change ratio and displacement balance ratio are established.

• bending moment change ratio

$$\gamma\_{\rm xy-z} = \frac{|M\_{\rm x\_{\rm y}} - M\_{\rm x\_{\rm 1y}}|}{|M\_{\rm x\_{\rm 1y}}|} \times 100\% \tag{5}$$

• bending moment balance ratio

$$\delta\_{xy} = \frac{2\left|M\_{\mathbf{x}\_ny - I} - M\_{\mathbf{x}\_ny - r}\right|}{\left|M\_{\mathbf{x}\_ny - I}\right| + \left|M\_{\mathbf{x}\_ny - r}\right|} \times 100\% \tag{6}$$

• displacement change ratio

$$
\varepsilon\_{\rm xy-z} = \frac{|D\_{\rm x\_{\rm y}\rm y} - D\_{\rm x\_{\rm I}\rm y}|}{|D\_{\rm x\_{\rm I}\rm y}|} \times 100\% \tag{7}
$$

• displacement balance ratio

$$\varphi\_{xy} = \frac{2\left|D\_{x\_ny-l} - D\_{x\_ny-r}\right|}{\left(D\_{x\_ny-l} + D\_{x\_ny-r}\right)} \times 100\% \tag{8}$$

where, *Mxny*—In the *xn* scheme, the sum of absolute values of crown bending moments on both left and right sides at *y* distance from the arch bottom. *Mxny*−*l*—In the *xn* scheme, the bending moment value at *y* position of the left. In the *xn* scheme, the sum displacement of lining on both left and right sides is at *y* distance from the arch bottom. In the *xn* scheme, the displacement value of the right arch body is from the *y* position of the arch bottom. Where *x* is A~C, representing three types of schemes respectively; *n* is 2~7, representing 2~7 schemes respectively; *y* is the distance to the arch foot; *l* and *r* represent the left and right respectively.

The bending moment change ratio indicates the change of lining bending moment at a certain position. The greater the *γx*y−z, the greater the bending moment is affected. The bending moment balance ratio indicates the bending moment difference between the left and the right, and the larger the *δxy*, the more unbalanced the lining. The change ratio of displacement indicates the change of lining deformation at a certain position. The larger the *εx*y−z, the greater the deformation is influenced. Displacement balance ratio indicates the difference between left and right deformation. The greater the *δxy*, the worse the displacement balance.

#### **4. Analysis of Asymmetric Mechanical Behavior of Lining**

The section size of the −820 m pump house in Wanfu Coal Mine is 46.5 m2 (R1 = 3.5 m, R2 = 4.8 m, h = 5.3 m), the surrounding rock pressure is q1 = 100 kN/m, q2 = 150 kN/m, q3 = 50 kN/m, lining thickness is 0.5 m, and C40 reinforced concrete is adopted. The elastic modulus is 32.5 GPa, and the moment of inertia per unit length lining is 5.3 × <sup>10</sup>−<sup>3</sup> m4. By substituting the specific parameters into the corresponding Equations (1)–(4), the bending moment value M of each position can be calculated. To reduce the computational workload, 42 feature points are taken along the entire lining. On this basis, the influences of asymmetric coefficient, section size of chamber and lining thickness on lining bending moment and displacement are analyzed. In order to quantitatively evaluate the influence of various factors on the lining, the 45◦ position on both sides of the top arch is taken.

### *4.1. Asymmetric Coefficient*

#### 4.1.1. Bending Moment Analysis

The following result analysis can be obtained from Figure 5.

Under the action of symmetrical force, the bending moment of lining distributes symmetrically. With the increase of the asymmetry coefficient, the bending moment of the left arch crown is larger than the right arch crown, and the distribution of the bending moment of the bottom arch is opposite. The bending moment of the right straight wall is larger than the left.

(**a**) Bending moment distribution diagram (kN/m) (**b**) Bending moment of left and right shoulder

**Figure 5.** Bending moment analysis with different asymmetric coefficients. (**a**) Bending moment distribution diagram (**b**) Bending moment of left and right shoulder.

> Taking the left and right shoulder as an example, the bending moment of the left shoulder increases from 920 kN/m to 2104 kN/m and the bending moment of the right shoulder decreases from 920 kN/m to 625 kN/m, with the increase of the asymmetry coefficient from 1 to 1.6. The change ratio of bending moment of the shoulder is positively correlated with the asymmetric coefficient, which increases from 15% to 48%. The bending moment balance ratio is negatively correlated with the asymmetry coefficient, which decreases from 75% to 30%.

#### 4.1.2. Displacement Analysis

The following result analysis can be obtained from Figure 6.

(**a**) Displacement distribution diagram (mm) (**b**) Displacement of left and right shoulder

**Figure 6.** Displacement analysis of different asymmetric coefficients. (**a**) Displacement distribution diagram (**b**) Displacement of left and right shoulder.

> Under the action of the symmetric force, the displacement is distributed symmetrically. The maximum displacement of the straight wall is 35 mm, and the maximum displacement of the top arch is 26 mm. As the asymmetry coefficient increased, the overall lining

displacement increased, and the increase of the right lining displacement (the larger force side) was greater than the left.

Taking the left and right shoulder displacements as an example, the left shoulder displacements increased from 18 mm to 89 mm and the right shoulder displacements increased from 18 mm to 112 mm, with the increase of the asymmetry coefficient from 1 to 1.6. The displacement change ratio of the shoulder is positively correlated with the asymmetry coefficient, which increases from 70% to 437%. The displacement balance ratio is negatively correlated with the asymmetry coefficient and decreased from 60% to 42%.

*4.2. Section Size*

4.2.1. Bending Moment Analysis

The following result analysis can be obtained from Figure 7.

**Figure 7.** Bending moment analysis of different section sizes. (**a**) Bending moment distribution diagram (**b**) Bending moment of left and right shoulder.

> With the increase of the section size of the chamber, the lining bending moment increases, but the position of maximum bending moment does not change. The maximum bending moment of the arch crown is located at the left side of the arch top. The maximum bending moment of the straight wall is located at 2 h/3 to the wall foot. The maximum bending moment of the arch bottom is on the right of the arch bottom middle.

> Taking the left and right shoulder as an example, with the section size increasing from 4 m to 10 m, the bending moment of the left shoulder increases from 161 kN/m to 2006 kN/m, and the bending moment of the right shoulder increases from 12 kN/m to 298 kN/m. The bending moment change ratio increases from 118% to 1403% with the increase of section size. The bending moment balance ratio is negatively related with the section size, which decreases from 142% to 108%.

#### 4.2.2. Displacement Analysis

The following result analysis can be obtained from Figure 8.

With the increase in the section size, the displacement of the lining increases rapidly, and the displacement value of the right side increases more obviously than that left, but the maximum displacement position remains unchanged.

Taking the displacement of the left and right shoulder as an example, with the increase of the section width from 4 m to 10 m, the displacement of the left shoulder increases from 1.8 mm to 155.1 mm, and that of the right shoulder increases from 3.1 mm to 220.9 mm. The displacement change ratio increases from 12% to 810%, and the displacement balance ratio decreases from 51% to 42%. ս〫ᒣ㺑⦷

ᐖᤡ㛙 ਣᤡ㛙 ս〫ਈॆ⦷

(**a**) Displacement distribution diagram (mm) (**b**) Displacement of left and right shoulder

**Figure 8.** Displacement analysis of different section sizes. (**a**) Displacement distribution diagram (**b**) Displacement of left and right shoulder.

#### *4.3. Lining Thickness*

The following result analysis can be obtained from Figure 9.

(**a**) Displacement distribution diagram (mm) (**b**) Displacement of left and right shoulder

**Figure 9.** Displacement analysis of different lining thickness. (**a**) Displacement distribution diagram (**b**) Displacement of left and right shoulder.

> With the increase of the lining thickness, the overall displacement decreased, and the displacement of the right side (the side with larger force) decreases more significantly than that of the left, and the maximum displacement position does not change.

> Taking the displacement of left and right arches as an example, with the increase of arch thickness from 0.2 m to 0.8 m, the left arch shoulder displacement decreases from 550 mm to 10 mm, and the right arch shoulder displacement decreases from 816 mm to 15 mm. The displacement change ratio increases from 70% to 98%, and the displacement balance ratio increases from 39% to 42%.

> When the lining thickness exceeds 0.5 m, both displacement change ratio and displacement balance ratio have little change with the continuous increase of lining thickness.

#### *4.4. Summary and Engineering Recommendations*

Under the action of asymmetric load, the bending moment and displacement of the lining show an obvious asymmetry. The larger the asymmetry coefficient is, the more significant the increase of the bending moment and displacement is, and the more serious the failure is. Therefore, for the lining chamber subjected to secondary disturbance, the asymmetric design of bolt mesh spray should be carried out, the support strength of the disturbed side should be strengthened, and the strength of the surrounding rock should be restored through grouting and other measures to reduce the asymmetric load of the lining.

The larger the section size is, the more severely the lining is affected by the asymmetric load, the change rate of bending moment and displacement of the lining increases rapidly, and the symmetry decreases continuously. Therefore, for the lining under asymmetric load, measures such as optimizing equipment selection should be taken to reduce the chamber section size to minimize the size effect.

Increasing the lining thickness can improve the ability of lining to resist asymmetric load, but when the lining thickness exceeds a certain value, the effect of increasing the lining thickness is no longer obvious. Therefore, the section size and the load size of the chamber should be considered comprehensively, and the lining thickness should be reasonably selected by using the evaluation index.

#### **5. Control Measures and Application Effect**

#### *5.1. Control Measures*

Under the condition of asymmetric load, the core of stability control of the deep large section chamber lining is "strengthening asymmetric support, reducing the section size and improving the lining strength". In view of the asymmetric failure of the permanent pump house in Wanfu Coal Mine, the control method design is shown in Figure 10.

**Figure 10.** Schematic diagram of reinforcement scheme.


end of grouting anchor cable is anchored, 150 kN preload is applied first, and then grouting is carried out one by one to reinforce the broken surrounding rock.

<sup>3</sup> On the outside of the existing section, double-layer reinforcement is bound on the whole section, with the diameter of transverse reinforcement of 28 mm and longitudinal reinforcement of 32 mm. C40 concrete is used for pouring, and the lining thickness is 500 mm.

#### *5.2. Application Effect*

After the reinforcement following the above scheme, a deformation monitoring section is set, and the position of measuring point is the same as before. Figure 11 shows the displacement change curve of the pumping house within 90 days of reinforcement. It can be obtained that, after the reinforcement of the pumping house chamber, no obvious asymmetrical deformation continues to occur, and the maximum deformation is 7.3 mm, indicating a significant control effect.

**Figure 11.** Displacement monitoring of pump house after reinforcement.

For the permanent pump house, continuous deformation is not allowed during the reinforcement phase, and the rigid lining support is adopted to effectively avoid the continuous deformation of the pump house, which meets the support requirements. Figure 12 is the effect drawing of the pump house after reinforcement.

**Figure 12.** Field effect of pump house after reinforcement.

#### **6. Conclusions**

(1) Aiming at the asymmetric failure of lining in Wanfu Coal Mine pump house due to the secondary disturbance, the asymmetric mechanical model of the lining and the quantitative evaluation indexes such as the change rate of bending moment, the balance rate of bending moment, the change rate of displacement and the balance rate of displacement are established, and the influence laws of the asymmetric coefficient, the section size and the thickness on the bending moment and displacement of the lining are studied.


**Author Contributions:** All the authors contributed to this paper. Conceptualization, Y.H. and B.J.; methodology, H.W.; investigation, Y.H.; resources, J.Z.; data curation, X.G., Y.M.; writing original draft preparation, Y.H., Y.M. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was supported by the National Natural Science Foundation of China (Grant Nos. 52074164 and 42077267); the Natural Science Foundation of Shandong Province, China (Grant Nos. 2019SDZY04 and ZR2020JQ23); the Project of Shandong Province Higher Educational Youth Innovation Science and Technology Program, China (Grant No. 2019KJG013).

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The data supporting the findings can be found in the figures, tables and manuscript.

**Acknowledgments:** This work is supported by Wanfu Energy Co. Ltd. and Yankuang Group Limited, which are gratefully acknowledged.

**Conflicts of Interest:** The authors declare no conflict of interest.
