Abstract
To accurately detect the development height of the water flowing fractured zone (WFFZ) in the overlying strata of the working face after mining under water and to ensure the safety and reliability of coal mining, the coal seam located under Weishanhu Lake in the Jisan coal mine was used as the experimental system. A similar laboratory simulation and water injection-based fracturing test system were used with the working face before and after mining activity to calculate, quantitatively detect, and qualitatively analyze the development height of the WFFZ in the overlying strata. Meanwhile, a flow-stress-damage model and its criterion of fracture expansion were established based on the Mohr-Coulomb criterion, and the FLAC 3D software was used to simulate the deformation and failure of the overlying strata and the evolution of WFFZ during the mining process. The results showed that the height ranges of the WFFZ beneath Weishanhu Lake of the Jisan coal mine as established by the above three methods are 30-45 m, 30-48 m, and 30-50 m. In the process of mining, the caving zone and fractured zone are, respectively, subjected to tensile failure and shear failure. The development height of the water flowing through the fractured zone in the overlying strata is basically consistent with the range of the “breaking arch.” The flow-stress-damage model and its criterion of fracture expansion can be applied to the fracture law of overlying strata under water under similar geological conditions.
1. Introduction
China’s coal resources are widely distributed and are rather common under buildings, railways, and water bodies. The mining of coal resources under water is particularly seriously threatened by the overlying aquifer; as such, the mining of coal under water has become a focus of research. The primary hazard is mainly caused by the movement and failure of overlying strata and the development of water flowing fractured zone (WFFZ) during the mining process. The continuous development of the WFFZ will cause it to be linked to the main aquifer, and even the surface water system, of the mining area. Because of this, the surface, underground water system near the mining area, and the surface ecological environment will change. The key to achieving safe, green, and efficient production of coal resources under water is to ensure that the WFFZ is not affected by the main aquifer [1]. Therefore, it is of great practical significance to carefully study the development process of the WFFZ and the movement of overlying strata while mining under water.
Both domestic and international scholars have studied the damage of overlying strata and the accuracy of improving predictions of the height of the WFFZ from different perspectives. These studies have provided technical support for safe mining under water. The height and evolution of the WFFZ are often studied using empirical formulas, physical experiments, numerical simulations, and field tests. The hydraulic fracturing process involves the coupling of fluid movement, rock deformation, and fracture propagation. Fractures may form whenever stress in any part of the material exceeds its strength, which could be either shear stress, normal stress, or a combination of both [2]. The deformation of flexible strata during complete failure is not brittle; rather, brittle failures occur only as mining is in progress [3]. The mining-induced fracture zone provides pressure-relief passageways for gas that flows from coal seams and the surrounding coal/rock strata as well as spaces for gas storage [4, 5]. The mining-induced fracture zone, together with the caving zone, is defined as the WFFZ. Bai et al. [6] acquired an empirical formula for the height of the WFFZ through field tests in multiple coal mines under the condition of a stable gob. Majdi et al. [7] presented five new simple, yet conclusive, mathematical approaches to estimate the height of the distressed zone. Liu et al. [8] used on-site measurements, mechanical theory calculations, and numerical simulations to analyze the regularity of the height of the WFFZ. As the face continuously advances, the overlying strata collapse, and destruction and dissociation occur on the weak side [9]. Zhang et al. [10] analyzed the regularity of the development of the WFFZ under the influence of different factors in the process of block mining and established a height prediction model of the WFFZ for block mining. Other scholars have calculated the height of the WFFZ in the mining process according to the prediction formula presented in the Regulations of Mine Hydrogeology issued by the Ministry of Coal Industry. The above research has provided a foundation for coal mining under rivers and lakes. However, mining-induced fracturing and the deformation of overlying strata require further exploration from the perspective of mining under surface water and aquifers.
Based on the mining that occurs under Weishanhu Lake in the Jisan coal mine, as well as on the theoretical analysis of mining-induced fractures, deformation laws, and structural mechanics models of overlying strata, this paper provides a simulation-based analysis of the development height of the WFFZ on a working face before and after mining occurs. The established fracture expansion criterion based on the Mohr-Coulomb criterion is imbedded in the FLAC 3D software to simulate and analyze mining-induced fracturing and deformation laws of overlying strata and the regularity of the development of the WFFZ and to verify the engineering according to field-measured data in the mining process.
2. The Breaking Arch Model and Fracture Expansion Criterion
2.1. Mechanical Structure Model of the Breaking Arch
A stope mechanical structure model is established [11–13] as shown in Figure 1.
According to an existing theory, the overlying strata can be divided into three zones in the vertical direction: the caving, the fractured, and the sagging [14]. The breaking arch consists of the caving zone and the fractured zone, and the height of the WFFZ is basically consistent with the range of the breaking arch. (1)The Caving Zone. The caving and final arrangement of blocks in the stope after breaking are irregular. The bulking coefficient is relatively high, generally up to 1.3-1.5. However, after they are recompacted by overlying strata, the bulking coefficient can reach approximately 1.03. This zone is directly above the roof of the mining area, and in many cases, it is caused by immediate roof caving.(2)The Fractured Zone. After the strata are fractured, the area where rock blocks remain arranged in an orderly manner is called the fractured zone, which is located above the caving zone. Since the rock blocks are arranged in an orderly fashion, the bulking coefficient is relatively low.(3)The Sagging Zone. All strata from the top of the fractured zone to the surface are called the sagging zone. The notable characteristic of stratal movement in this zone is that the movement is continuous; it is also integral that the movement of strata from the top of the fractured zone to the surface is stratified and integral. In the vertical section, subsidence values of all upper and lower parts are quite low. If there is a thick and hard key layer, a separation area may emerge in the sagging zone.
The WFFZ is distributed in the shape of a “saddle” [15], including the caving zone and the fractured zone. If the overlying aquifer is located within the WFFZ, it will cause the fracture water to flow from the fractured strata into the gob and working face.
2.2. Fracture Expansion Criterion
Due to the interaction of seepage pressure and coal seam mining under water, the strata will fracture and form a fracture zone that utilizes existing holes in the strata, which causes the strata to shear and slide along the fracture zone, and eventually causes the destruction of the entire structure.
2.2.1. The Range of the Fractured Zone Generated during Mining
The working face mining process is divided into two stages, including (1) the stage of incomplete mining, during which the mining distance of the working face is less than the width of the working face, and (2) the full mining stage, during which the mining distance of the working face is greater than the width of the working face.
The fracture height of the overlying strata in the stope is determined by the width of the working face. When this is fixed, the maximum height of the WFFZ of the overlying strata is definite. When the mining distance of the working face does not reach the width of the working face, the development height of the overlying strata space structure is related to the length of the working face; overall, the height of the overlying strata space structure in the stope develops linearly as the work face is mined. When the mining distance of the working face reaches the width of the working face, the development height of the overlying strata space structure is approximately one-half of the width of the working face [16].
2.2.2. Fracture Expansion Criterion
As a natural material, pores and cracks exist in rock masses and have significant influence on the mechanical properties of the host rock mass. When the rock mass is subjected to external forces, it often breaks down in such defective areas first, followed by instability and failure of the entire rock mass [17]. Fracturing can be classified as three basic types based on how the fractures form under external forces: opening mode (Type I), sliding mode (Type II), and tearing mode (Type III). The splitting that is caused by water pressure is different from other mechanical splitting because any newly generated fracture will be immediately filled with water and that water pressure is applied to the fracture surface. Fractures formed under tension are Type I fractures. According to the Griffith theory [18, 19], the appearance of a fracture will produce a new surface for the solid material, which has the same surface energy as the surface of the liquid, and a part of the energy (U) released by the system will be converted to surface energy. Assume that the surface energy on the surface of the material is , namely, surface tension, and the maximum length of the rock is ; then, the surface energy is . When the fracture is in an equilibrium state , the left side of the formula is the force that drives fracture expansion, while the right side of the formula is the force with which the material hinders fracture expansion. Griffith believes that if the applied stress exceeds the critical stress value or the input energy exceeds the critical energy, the material will fracture. If the fluid in the fracture is continuously supplemented through permeation, the fracture pressure is maintained so that the external input energy continues to be greater than the critical energy; thus, part of the fracture will continue to expand.
2.2.3. Numerical Simulation Method for Fracture Expansion
From the perspective of energy, the fracturing of rock material mainly occurs when the input energy exceeds the maximum value of the material itself; thus, the material will begin to break down; the more energy the material bears, the higher the degree of damage. In the initial fracture, since the crack tip is not very sharp and the degree of energy accumulation is low, the cracking is slow; however, the crack tip becomes very sharp after cracking, and a large amount of energy will accumulate at the crack tip, resulting in fracturing of the rock. Wang et al. [20] presented the acoustic emission hits can be used to analyze the damage in rock material, but the number of acoustic emission hits cannot be used alone to determine the degree of rock damage directly. Therefore, it is an appropriate method to judge the degree of cracking from the perspective of energy.
As the rock is damaged, the microstructure changes, the material strength and the elastic modulus decrease, and the permeability becomes higher. The damage variable is defined to estimate the degree of damage. According to the ratio of dissipated energy to the maximum strain energy of rock, the degree of damage to the rock is estimated. The strain energy of the rock is as follows:
For the convenience of simple calculations, we believe that the stress-strain curve for the rock can be simplified into a bilinear stress-strain model that consists of two straight lines (Figure 2). The area surrounded by the stress-strain curve is the energy that the unit body bears. The input energy for the element is partly used to increase the elastic energy of the material and is also partly consumed by other kinds of energy, such as sound waves and crack propagation. Theoretically, this is called the dissipated strain energy density of the element, which is shown as follows: where is Young’s modulus after reduction, which is related to the stress path. When the energy of the element is greater than the critical energy, cracks appear, which are generally oriented perpendicular to or the maximum tensile stress.
Based on the built-in Fish function of FLAC 3D, the damage model is built with the energy defined as the determining criterion of fracturing. The simulation progress is shown in Figure 3. At the initial stage, the strain energy density of the element is set to 0 (that is, ), and the strain energy density of the element calculation step I is as follows:
The strain energy accumulation of the element during the simulation process is recorded. , , , , , and are the element stresses in step I; , , , , , and are the element stresses prior to step I. Stress-strain is indicated in the same way.
If the element energy meets , the element will not be damaged.
If the element energy meets , the damage variable will be defined as the ratio of to .
If , the element is considered to be completely damaged, and the bearing capacity is lost. To ensure that the element does not produce singular points, the element is given a residual volume modulus and a larger permeability .
It is necessary to note that the loss of the bearing capacity of the element does not indicate the host rock sample is destroyed. Under three-directional pressure, the cohesion of the material is lost, but overall damage will not occur; the rock will be destroyed only when the confining pressure is withdrawn. Therefore, it is practical and appropriate to give the damage element a smaller volume variable and higher permeability.
Therefore, the FLAC 3D rock mass energy-determining criterion can be used to comprehensively estimate the rock mass damage caused by mining-induced stress and fracture water pressure. Meanwhile, the development of the WFFZ can be simulated well.
3. Methods and Results
3.1. Numerical Simulation of the Development Regularity of the WFFZ
3.1.1. Model Configuration
Based on the 183upper04 working face that is located under Weishanhu Lake of the Jisan coal mine, the FLAC 3D simulation software was used to build a numerical model to study the movement of overlying strata and the development of the WFFZ during mining. According to the test results of rock samples, the strata overlying the coal seam are mainly sandstone, with fine-grained sandstone and medium-grained sandstone in the upper part and thick siltstone and mudstone in the lower part, where fractures are well-developed. Table 1 lists the main lithologies of the 183upper04 working face. The fracture parameters are obtained according to the properties of the rock materials [21].
The length, width, and height of the design model were 800 m, 400 m, and 200 m, respectively, as shown in Figure 4. The boundary conditions of the model mechanics were as follows: the top boundary of the model was a free boundary; the bottom, front, back, left, and right boundaries of the model were set as stress boundaries.
In the FLAC 3D model, the S-B method was applied to initialize the ground stress by applying the stress boundary; the equivalent load, namely, the gravity stress, was applied at the top of the model. The load is obtained by the following formula: where indicates the unit weight of overlying strata and indicates the depth from the top boundary of the model to the surface. The lateral stress generated by gravity stress in the horizontal direction is determined by the following formula: where is the side pressure coefficient, determined by the formula , where is Poisson’s ratio.
After modeling, the deformation and failure conditions of overlying strata were simulated, and the height and shape of the caving zone and WFFZ during mining were analyzed. During the simulation, the mining height is 3 m, and the open-off cuts are made from elevations of -715 m, -710 m, and -705 m, with each step of the excavation being 50 m. Then, the development of the WFFZ was analyzed.
3.1.2. Analysis of Simulation Results
(1) The Distribution of the Failure Field at Different Mining Distances. By analyzing the failure field, the distribution of failure areas in overlying strata after mining can be directly observed. Based on Figure 5, which shows the simulation results of the distribution of the plastic zone at different mining distances of the working face, when the mining of the working face reaches 50 m, the plastic zone is concentrated directly above the upper gob, which gradually decreases and mainly belongs to shear failure. The siltstone strata below are dominated by tensile failure. As the mining of the working face reaches 100 m and the gob area increases, the coal seam roof falls and sinks further. Tensile failure occurs mainly in the immediate roof above the gob, and tensile shear failure above the coal pillar indicates that the collapse of the immediate roof occurs directly from this area. Tensile shear failure occurs in the main roof stratum above the immediate roof, and shear failure becomes more obvious upward. The range of the plastic zone increases in the horizontal and vertical directions as the working face is mined further. Finally, the heights of the caving zone and fracture zone are approximately 15 m and 50 m, respectively, which develops into the mudstone with 8.5 m height. As the mining of the working face reaches 150 m or more, the development height of the caving zone and the WFFZ basically remains stable. In summary, the distribution of the plastic zone above the gob and above the coal wall in the front and back of the gob basically corresponds to the development height of the caving zone and the WFFZ. When open-off cuts are made from elevations of -715 m, -710 m, and -705 m, there is little change in geological conditions, and the regularity of the WFFZ is similar.
(a)
(b)
(c)
(2) Distribution of the Failure Field in Different Mining Upper Limits. For open-off cuts made at different elevations, the simulation results of the distribution of the plastic zone in overlying strata after excavating 200 m are shown in Figure 6. From the simulation results of different mining upper limits, the displacement vector map and contour of displacement of each mining upper limit are similar in shape with consistent stratal movement. After mining, stress is concentrated at the edge of the gob and shear failure occurs in the rock mass above the coal pillar. According to the distribution of the maximum and minimum principal stresses, the coal seam roof area is a destressed zone, where stress is redistributed after the coal seam is mined. The rock mass is dominated by tensile failure. On the upper gob, tensile failure, tensile shear failure, and shear failure zones develop successively from the bottom to the top. As the mining area increases, the range of the maximum principal stress contour constantly expands, the shape and height constantly change, and the plastic failure zone of the strata continues to develop forward and upward, indicating that the caving zone and the WFFZ are increasing continuously.
(a)
(b)
(c)
Overall, if the mining height is 3 m, the height of the caving zone is approximately 15 m and the height of the WFFZ is approximately 50 m.
We investigate two other methods to calculate, quantitatively detect, and qualitatively analyze the development height of the WFFZ in overlying strata to verify the practicality of the numerical model.
3.2. Physical Experiment on the Development Regularity of the WFFZ
3.2.1. Similar Simulation Test Design
A similar simulation test model was built based on the 3upper coal seam in the Jisan coal mine. During the simulation, the movement, failure, and distribution of stress in the overlying strata were continuously monitored. The mechanical phenomena and rules in site mining were obtained by inversion. Eighty-meter boundaries were left on both sides of the 183upper04 working face at the tail entry and head entry. The physical similarity model is shown in Figure 7. The length, width, and height of the model were , the geometric similarity parameter of the model was , the unit weight similarity parameter was , the stress similarity parameter was , and the elastic modulus similarity parameter was . Table 2 lists the main model similar parameters.
3.2.2. Mining and Field Monitoring
The mining height in the model was 1.5 cm, equivalent to 3 m in an actual mining scenario. Forty centimeter-wide boundaries were left on the left and right sides of the model, equivalent to 80 m boundaries in an actual mining scenario. Mining was performed from left to right.
To monitor the settlement of overlying strata, the coordinate grid and the light lens displacement meter were used for calibration, and the plane displacement monitoring grid was arranged on the front side of the model. The layout of monitoring spots is shown with red dots in Figure 7.
The mining process of the model was monitored, the movement of overlying strata during mining and the parameters in the failure state of the overlying strata were recorded, and the images were processed.
3.2.3. Experimental Results and Analysis
Excavation is begun from 40 cm on the left boundary in 5 cm increments, which is equivalent to actual excavation increments of 10 m. The height variation and location of the boundary of the WFFZ during the process of model mining are measured and recorded, and the relationship between the development height and the mining distance is also recorded. The physical model is shown in Figure 8.
Each time stratal movement is stabilized after mining, all points of displacement and failure on the model are observed and calculated comprehensively. The analysis of experimental results for the 183upper04 working face is shown in Figure 9.
(a)
(b)
(c)
(d)
(e)
(f)
During mining, the development height of the WFFZ increases with the mining distance, but when the mining distance exceeds a certain value, it no longer develops upward. Figure 10 shows the relationship between the failure height of the overlying strata and the mined distance of the working face.
When the mined distance of the working face is short, the overlying strata can maintain a state of stability without extensive rock collapse; the development height of the WFFZ also remains low, as shown in Figure 9(b). As mining progresses, the development heights of the caving zone and the WFFZ increase rapidly with the first weighing of the working face and the rapid collapse of overlying strata. Afterward, the development height of the caving zone gradually falls and stabilizes, with a maximum height reaching 12 m, as shown in Figure 9(c). As the mining distance increases, the development height of the WFFZ continues to increase steadily, and the fractures rapidly spread upward. When the working face is mined to 150 m, its maximum development height reaches 44.14 m, as shown in Figure 9(e). The ratio of the height of the fractured zone to the mining height is 14.71. Meanwhile, an obvious curved subsidence zone is formed above the WFFZ. Subsequently, with a further increase in the mining length, the development height of the WFFZ gradually decreases and stabilizes between 30 and 45 m.
From Figure 9(e), the development of the WFFZ is distributed in the shape of a “saddle,” in which the shape is slightly higher on both sides and lower in the middle. The reduction area appears in the middle of the gob, which corresponds to the results of the theoretical study of the failure of overlying strata. The broken overlying strata in the gob are gradually compacted as the working face advances, the gob boundary is supported by the coal wall, and the degree of reconsolidation of the fractured and caved strata in the middle of the gob is greater than that of the gob boundary.
As Figure 10 shows, the caving zone and the WFFZ follow a pattern of increasing first and then decreasing and stabilizing as the distance of coal mining increases. The height of the caving zone of the model is stable at a range of 8-12 m, and the development height of the WFFZ is 44.14 m. Compared to the strike model, the tilt angle is the main cause of the higher height of the WFFZ.
3.3. Field Measurements of the WFFZ
3.3.1. Application of Water Injection Fracturing Test System
The “underground up-hole observation device” was arranged in a chamber in a certain position around the downhole mining face, and the upward-inclined drill holes were pitted from the chamber to the WFFZ in the strata overlying the gob of the working face. A double-sided water plugging device was used to observe the development height of the WFFZ in the overlying strata before and after mining and to analyze the deformation and failure of the overlying strata, ultimately providing scientific technical parameters for the safety mining under water [22].
The observation method used for borehole injection water leakage in the surface is a traditional and reliable method. However, when mining under large water bodies, the traditional method of observing the height of the WFFZ through ground drilling cannot be performed due to the presence of a water body. In this study, underground water leakage was measured by up-hole injection. That is, at the periphery of the working face, up-holes were drilled upward to the gob; leakage measuring by subsection, sealing, and water injection was performed in the drill holes to determine the height of the WFFZ according to seepage in the drill holes at each depth. According to the existing roadway layout, the height observation station was selected on the side of the stopping line of the working face. In the auxiliary haulage gate of the 183upper04 working face, appropriate observation profiles were selected to arrange premining observation holes and postmining observation holes. The observation results from premining observation hole 1 and postmining observation hole 2 are shown in Figures 11 and 12, respectively. The injection water leakage of observation hole 1 is generally between 0.2 and 0.5 L/min. The drilling is not affected by mining action, and the data mainly fluctuate because the overlying strata of the 3upper coal seam include interbedded siltstone, fine sandstone, medium sandstone, and mudstone, for which the fracture development and water conductivity of each stratum are different, leading to different injection water leakages in each stratum. The entire pore segment shows different injection water leakages, which indicates that the upper strata that were exposed by premining drilling are relatively intact and fractures are not developed. The injection water leakage in observation hole 2 is obviously higher than that in observation hole 1, which indicates that fractures are better developed within the drilling depth of 60 m. As the drilling depth increases, the injection water leakage in the pore segment slowly increases until it is approximately consistent with observation hole 1, indicating that the degree of fracture development is reduced. Therefore, from observation hole 2, the upper limit of water conductivity is the depth of the hole at 60 m, and the height of the WFFZ is calculated by the following formula:
Based on the test results of all observation holes, the height of WFFZ in the strata overlying the working face is approximately 40 m.
According to the observation results, the development height of the WFFZ on the first mining face of the 3upper coal seam located under Weishanhu Lake of the Jisan coal mine is saddle-shaped, as shown in Figure 12. The width of the side expansion (or extrusion) of the WFFZ is approximately 10 m. The upper surface of the WFFZ is in the mudstone, which is in the eighth layer overlying the 3upper coal seam with a thickness of 8.5 m. It is believed that its upper surface is the boundary between the 8.5 m thick mudstone and the 1.6 m thick fine and medium-grained sandstone.
The reason for the formation of the saddle-shape and side boundary convex-shape is that the curvature of the overlying strata is the largest at the mining boundary, so the height of the WFFZ at the mining boundary is also the highest. On the outside of the mining boundary, that is, above the coal pillar, the overlying strata are in a state of tensile stress, which easily produces opening fractures. Therefore, the outer convex type is caused on the side boundary of the WFFZ.
3.3.2. Visual Observation of Boreholes
The fracture development of the 183upper04 working face roof was observed by a YTJ20 borehole television (Figure 13); the observation depth was 13.43 m due to limited equipment conditions. According to the nearest C4-11 bore hole columnar section to the chamber, the borehole camera passes through the 1.75 m thick siltstone of the roof, the 2.95 m thick fine siltstone interbed, and 5.22 m thick claystone in sequence, and the ultimate depth is at 5.22 m in the claystone. Based on the analysis of the borehole video, the height of the caving zone on the 183upper04 working face is greater than or equal to 9.92 m.
The heights of the WFFZ obtained using the above three methods were approximately 50 m, 30-45 m, and 40 m. The respective caving zone heights were approximately 15 m, 8-12 m, and 9.92 m. The results of the physical experiment and field measurements showed that the numerical simulation results are correct. The results verified the practicality of the flow-stress-damage model to study the fracture of overlying strata under water during mining. However, there is an obvious bias in the data obtained from the three methods; the data obtained from the physical experiment and field measurements were slightly lower than those of numerical simulation. The possible reasons mainly include the selected mechanical parameters of the strata, which have a significant influence on the experimental results. Another possible reason is that the caving zone height results are inaccurate because of limited equipment conditions that were presented in the field measurement section.
4. Conclusions
Based on the Griffith theory, energy was used as the fracture expansion criterion. Through the use of the Fish language, the flow-stress-damage model and its criterion were embedded in FLAC 3D software, and a calculation program was applied to simulate the mining process under a water body. The failure of the surrounding rock in the mining process under the water body was analyzed. Combined with the field observation results and laboratory simulation results for the height of the WFFZ, the effectiveness of the determination method was verified and the evolution law of the height of the WFFZ was obtained.
The WFFZ consists of the caving zone and the fracture zone, and its development height is basically the same as that of the “breaking arch.” As the working face is mined, caving occurs in the overlying strata and the degree of fracture development increases, but the height of the WFFZ will not increase if the overlying strata are supported by waste rock, which is formed after caving and expansion of the working face. The height of the WFFZ obtained by numerical simulations was approximately 50 m, the development height obtained by laboratory simulations was in the range of 30 m to 45 m, and the development height obtained by the water injection fracturing test system was approximately 40 m.
The key to safe mining under water is to ensure that the WFFZ does not affect the overlying water. The height of the WFFZ measured in the field was within the safe range, and thus, safe mining under water can be realized. The research results provide theoretical support for coal mining under water in the Jisan coal mine and a reference for determining the height of the WFFZ in similar geological conditions.
Data Availability
The tables, some figures, and data used to support the findings of this study are currently under embargo while the research findings are commercialized. Requests for data six months after the publication of this article will be considered by the corresponding author.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
This research was financially supported by the State Key Research Development Program of China (No. 2016YFC0600708), the Open Fund for State Key Laboratory (SHGF-18-13-30), the Tai’shan Scholar Engineering Construction Fund of the Shandong Province of China, the Tai’shan Scholar Talent Team Support Plan for Advanced and Unique Discipline Areas, the Shandong Province Higher Educational Science and Technology Program (No. J15LH04), the Natural Science Foundation of Shandong Province (No. ZR2018 MEE001), the State Key Laboratory of Open Funds (No. MDPC201601), the Source Innovation Programme (No. 18-2-2-68-jch).