**1. Introduction**

Hydrogen bonding, both inter- and intramolecular, has profound effects on the properties of small molecules, reactivity, pKa values, polarity, solubility and with that, e.g., penetration of membranes and biological effects in general. Most drugs and bioactive molecules are small molecules and often depend on binding to receptors partly via hydrogen bonds [1,2]. Related to intramolecular hydrogen bonding is tautomerism. One of the effective tools in the study of hydrogen bonding is isotope effects on NMR chemical shifts. Both primary and secondary isotope effects may be used. The most common isotope is 2H, deuterium (D), but also 13C, 15N and 18O as well as more rare isotopes are used. For observation of secondary isotope effects, the most common nuclei are 1H, 13C, 15N, 19F and 18O. In case of primary isotope effects, the pair 1H, 2H is mostly used. The intrinsic isotope effects are defined as:

Primary: PΔH(D) = δ(H) − δ(D);

Secondary: nΔX(h) = δX(l) − δX(h).

l = light, h = heavy and X is the nucleus under investigation, n is the number of covalent bonds between the isotope and the investigated nucleus. For intramolecular hydrogen bonded cases, the isotope effects may be transmitted via the hydrogen bond. In the cases in which to opposite sign convention has been used, the sign is changed and the number is marked with an asterisk.

Deuterium isotope effects can be related to hydrogen bond strength and to hydrogen bond energies. Isotope effects are very useful in the equilibria of symmetric systems (lifting of degeneracy). They respond to nearby charges.

Most of the studies were performed in the liquid state, but isotope effects on chemical shifts can also in favorable cases be studied in the solid state.

The present review will primarily cover the last ten years. Other reviews that cover the subject are [3–6].

**Citation:** Hansen, P.E. Isotope Effects on Chemical Shifts in the Study of Hydrogen Bonds in Small Molecules. *Molecules* **2022**, *27*, 2405. https:// doi.org/10.3390/molecules27082405

Academic Editor: Miroslaw Jablonski

Received: 11 March 2022 Accepted: 5 April 2022 Published: 8 April 2022

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

**Copyright:** © 2022 by the author. 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/).


*Article* **Numerical Study on Characteristics of Bedrock and Surface Failure in Mining of Shallow-Buried MCS**

**Guangchun Liu 1,2, Wenzhi Zhang 1,\*, Youfeng Zou 1, Huabin Chai 1,\* and Yongan Xue 1**

**Abstract:** Coal is one of the important energy sources for industry. When it is mined, it will cause the destruction of bedrock and surface. However, it is more severe in mining shallow-buried multi coal seams (SBMCS). To better reveal the characteristics of the bedrock and surface damage, we have carried out a theoretical analysis, as well as used numerical simulations and field monitoring methods to study the surface and bedrock damage caused by the mining of SBMCS. The characteristics of bedrock and surface failure structure, settlement, and stress distribution were studied and analyzed. The findings show that the collapsed block, formed by the rupture of the overlying stratum, interacts with the surrounding rock to form large cavities and gaps, and the stress concentration occurs between them. The maximum downward vertical concentration stress is about 9.79 MPa. The mining of the lower coal seam can lead to repeated failure of the upper bedrock and goaf. The settlement of bedrock presents gradient change, and the settlement of upper bedrock is large, about 8.0 m, and the maximum settlement is 8.183 m, while that of lower bedrock is small and about 3.5–4.0 m. The weak rock stratum in the bedrock is crushed by the change stress of repeated mining, and formed a broken rock stratum. The cracks in the bedrock develop directly to the ground. On the ground, tensile cracks, compression uplift, stepped cracks, and even collapse pits are easy to cause in mining SBMCS. Affected by repeated mining, the variation of surface vertical stress is complex and disorderly in the middle of the basin, and the variation of horizontal stress is mainly concentrated on the edge of the basin. The maximum stress reaches 100 KPa, and the minimum stress is about 78 KPa. Through theoretical analysis and discussion, the size of the key blocks is directly related to the thickness and strength of the rock stratum.

**Keywords:** coal mining; multi coal seam; CDEM; numerical simulation; bedrock failure; surface damage

#### **1. Introduction**

Coal is a non-renewable fuel energy, which is gradually replaced by renewable energy. However, in today's world, it is still one of the most important industrial energy sources. Coal resource exploitation still plays a vital role in the energy planning and management of all countries in the world. The safe and efficient mining of coal is related to the national industrial and economic development. After coal mining, the bedrock layer and surface are damaged, affecting and threatening the safety of working face and surface. Coal fields are widely distributed in China, and the occurrence conditions of coal seams are different [1,2]. Qinshui coalfield, in Shanxi Province, is the basin coalfield with the most abundant coal in North China. It belongs to the coal-bearing strata of the Benxi Formation, Taiyuan Formation, and Shanxi formation developed in the late Paleozoic [3,4]. It is a typical geological condition for the occurrence of multiple coal seams. The mining thickness of coal seams in Qinshui coalfield varies greatly. According to [1], coal seam in Benxi

**Citation:** Liu, G.; Zhang, W.; Zou, Y.; Chai, H.; Xue, Y. Numerical Study on Characteristics of Bedrock and Surface Failure in Mining of Shallow-Buried MCS. *Energies* **2022**, *15*, 3446. https://doi.org/10.3390/ en15093446

Academic Editors: Longjun Dong, Yanlin Zhao, Wenxue Chen and Vamegh Rasouli

Received: 9 April 2022 Accepted: 6 May 2022 Published: 9 May 2022

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

**Copyright:** © 2022 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/).

Formation is thick in the Northeast and West, and thin in the central and South; coal seam in Taiyuan formation is thick in the East and thin in the West; coal seam in Shanxi formation is the thickest coal bearing group, which is thick in the North, thin in the South, thick in the East, and thin in the West. The geological conditions of shallow burial multi coal seams (SBMCS) refer to the geological conditions in which the coal seams are buried in a shallow depth, the composition and structure of the bedrock are relatively simple, and there are three or more coal seam occurrence geological conditions. The mining of SBMCS is affected by factors such as coal seam depth, coal seam thickness, bedrock and interlayer lithology, geological structure, and fault [5]. Many scholars have conducted in-depth research on multi coal seam (MCS) mining, under different geological mining conditions, from the perspectives of bedrock failure structure characteristics, roof structure characteristics, rock movement law and stability, and fracture development characteristics, and they achieved fruitful research results. Zhang [6], Lin [7] and others studied the characteristics of overburden movement and fracture development, affected by repeated mining in MCS, by using a similar simulation method, and they determined the relationship between fracture development height, overburden lithology, and thickness. The fracture development, rock fracture, and water loss of MCS mining are studied. The stress field distribution of bedrock and surface needs to be further discussed. Qin [8], Dong [9], Yang [10], and others used a combination of various methods to study the structural characteristics of mining MCS and analyzed roof instability and working face pressure, the calculation method of coal mine support resistance and the breaking impact effect of roof structure, and obtained the laws of MCS bedrock failure and roof structure characteristics. Based on the study of coal seam group, the stability and rock failure mechanical characteristics are analyzed. The bedrock and surface damage caused by SBMCS mining need to be further studied. Chen [11], Yan [12], Dong [13], and others studied the surface movement and deformation characteristics, surface dynamic prediction, bedrock damage, and migration law of SBMCS mining by means of numerical simulation, similarity test, and field monitoring. They found the relationship between the surface movement, bedrock damage, and bedrock fracture development of SBMCS and the bedrock buried depth and lithology. Their research on SBMCS is more comprehensive, through on–site monitoring, numerical simulation and theoretical analysis. However, they are aimed at the geological conditions with thick loose layers, and most of the research areas are located in the desert area of Northwest China. According to the geological conditions of Qinshui coalfield, further research is needed. Xie [14] et al. studied the basic roof structure characteristics, movement, and deformation law of MCS mining, providing the basis for the mechanical calculation of bedrock movement and deformation. The equivalent basic roof structure and activity law of overburden between MCS are studied. The mechanical analysis of surface failure and the fracture structure of bedrock in MCS are not studied in detail. Zhao [15] and others used the discrete element UDEC numerical simulation software and mining subsidence prediction software to simulate and study the surface residual settlement after mining, which provided a theoretical feasibility demonstration for coal mining under a high-speed railway. Wu [16] et al. Used FLAC3D to simulate the conditions of different bedrock heights and determined the impact of repeated mining on the height of "two zones". UDEC is a discrete element method, which can better simulate the characteristics of an overburden collapse zone caused by mining, and it has certain limitations on the simulation of initial stress before mining and the interaction between broken rock and soil. FLAC belongs to the finite element method and adopts the algorithm based on the Lagrange continuity equation. It also has certain limitations in the simulation of collapse and cracking. The above research shows that the mining of SBMCS is a complex problem related to geological conditions. The bedrock and surface failure characteristics and stress distribution law of SBMCS mining, under the geological conditions of Qinshui coalfield, need to be studied further. In numerical simulation, the discreteness of bedrock failure and the continuity of initial stress after mining cannot be well treated by using the finite element or discrete element method, respectively. Discrete element and finite element

methods are widely used in geotechnical engineering and underground engineering. Then, before mining, the rock stratum is a continuous medium, and it becomes a discrete block structure after mining. Therefore, a new method of coupling a finite element and discrete element, which is Continuous Discontinuous Element Method (CDEM), can solve this problem [17,18].

This paper attempts to use the CDEM and Generalized Digital Emulation Manufacturer(GDEM) software numerical simulation platform to study the bedrock failure characteristics, structural form, and bedrock stress distribution characteristics of SBMCS [19]. We take Majiazhuang coal mine in the southeast of Qinshui coalfield as the engineering background. It provides a numerical analysis basis for SBMCS mining in the southeast of Qinshui coalfield, and it provides a theoretical basis for safe mining and green mining of mines with MCS geological conditions. The remainder of this manuscript is organized as follows. In Section 2, Methodology, we mainly expound the mechanical analysis of key blocks and CDEM algorithm. In Section 3, Materials and Methods, we mainly describe the engineering background of the engineering research area, modeling, and simulation process. In Section 4, Results, we analyzed the failure characteristics of bedrock, surface stress distribution, and settlement, according to the results of numerical simulation and field monitoring. In Section 5, Discussion, we obtain bedrock failure characteristics and processes, and we discuss the impact of water on mining and the CDEM numerical simulation method. In Section 6, Conclusions, we describe the research results for the research content, as well as the conclusions obtained.

#### **2. Methodology**

#### *2.1. The Calculation of Key Stratum Failure*

After the coal seam is mined, the key layer in the rock stratum can support the overlying rock. However, the key layer will collapse when it reaches the ultimate bearing capacity [20]. The key block structure is formed after collapse, which interacts with the rock stratum and floor to form a stable structure. According to relevant literature research [21], it is found that, after the coal seam is mined, the key block plays a very important role. It forms a relatively stable structure to support the rock mass on it. We need to establish a mechanical model to deeply study and discuss the process.

Mechanical analysis shall be carried out from the initial stage of key block collapse. With the continuous advancement of the working face, the direct roof of the upper coal seam reaches the limit span, resulting in separation and collapse [22]. The key stratum bears the upper bedrock load. When the key stratum reaches the ultimate bearing capacity, cracking and disconnection will occur at both ends, as well as midspan of the rock beam, in the key stratum [23]. Under the action of friction and extrusion force between key blocks, the blocks produce rotary motion and finally reach a static equilibrium state, as shown in Figure 1.

**Figure 1.** Force analysis of the basic structure of the key stratum failure.

After bedrock collapse, the key blocks interact with the coal seam floor and surrounding rock to form a relatively stable structure. According to the "masonry beam" theory [24,25], there is a height difference between the horizontal thrusts *TA* and *TB* on both sides, and the compressive stress passes through the overhanging fault block to form the rear arch. The compression force *TC* in the fault block increases with the increase in the overhanging area of the block, where the two blocks are squeezed and broken, and the friction force of the block decreases. When the friction force between the two blocks is less than the shear force on the contact surface, the broken block slides down and sinks along the contact surface, resulting in rotational deformation and instability. The key stratum fractured for the first time, and the formed key block constituted an asymmetric three–hinge arch structure [26,27].

As for the geometry in Figure 1, assuming that the rotation angle of the key block is *β*, the subsidence amount is *w*, the length of the key block is equal to *L*, the upper part of the block is uniformly distributed with load *q*, and the generated horizontal thrusts are *TA* and *TB*, it functions with points *A* and *B*. The thickness of the block is *hm*, point *A* is subjected to frictional force *RA*, and point *B* is subjected to a supporting force *qB* of upward distance *d*. Point *C* is the hinge point of the fractured block, which transmits the internal force *TC* between the two block systems. It is a force per unit length (1 m).

According to the overall equilibrium state of the two key block structures, the equilibrium equation can be obtained:

$$\begin{aligned} \begin{cases} \sum F = 0\\ \sum M = 0 \end{cases} \end{aligned} \tag{1}$$

where *F* is the force in different directions, N;

*M* is the moment generated by each force, N·m;

According to the force analysis of the key blocks, as shown in Figure 1 and Formula (1), *TC* can be calculated that:

$$T\_{\mathbb{C}} = \frac{qL(L\cos\beta + h\_m\sin\beta)}{2(h\_m\cos\beta - w)}\tag{2}$$

where *TC* is the internal force between the two blocks, N/m;

*q* is the load above the block, MPa;

*β* is the tilting rotation angle of the block, ◦;

*hm* is the height of the block rock beam, m;

*W* is the subsidence caused by block rotation, m.

According to the geometric relationship, as shown in Figure 1, we can be known that

$$\begin{cases} \begin{array}{c} \sin \beta = \frac{w}{L} \\ \cos \beta = \frac{\sqrt{L^2 - w^2}}{L} \end{array} \tag{3} \end{cases} \tag{3}$$

Substituting Formula (3) into Formula (2), the relationship between thrust *TC* and settlement *w* can be calculated:

$$T\_C = \frac{qL}{2} \frac{\sqrt{L^2 - w^2} + \frac{h\_{\text{max}}}{L}w}{\frac{h\_{\text{m}}}{L}\sqrt{L^2 - w^2} - w} \tag{4}$$

In the initial stage of key stratum failure, when periodic failure occurs, the rotation angle *β* is usually small, and the ratio of rock layer thickness to failure distance (*hm/L*) is usually small and can be ignored in the calculation process. We can get:

$$\begin{array}{c} \sqrt{L^2 - w^2} \approx L \\ \frac{\eta\_{\text{fin}}}{L} \approx (\frac{1}{90} \sim \frac{1}{20}) \end{array} \tag{5}$$

Therefore, Formula (4) can be approximately simplified to

$$T\_{\mathbb{C}} = \frac{qL}{2} \frac{L}{h\_m - w} \tag{6}$$

According to the relevant knowledge of mechanics, rock beam collapse is usually caused by tensile failure in the middle, so the relationship between the initial collapse distance *Lini* and the tensile strength of rock stratum *σ*max can be calculated as follows [11,24].

$$L\_{\rm ini} = h\_{\rm m} \sqrt{\frac{\sigma\_{\rm max}}{q}} \tag{7}$$

where *Lini* is the initial collapse distance of bedrock, m;

*q* is the load applied on the collapse block, Mpa;

*σ*max is the maximum tensile strength of broken block, Mpa.

From the above analysis, it can be concluded that the size of the key blocks is directly related to the thickness and strength of the rock stratum, as well as the load on the bedrock.

The failure of key stratum has a direct relationship with the properties of the rock formation, and the failure distance to reach the initial failure can be calculated, as shown in Formula (7). Therefore, the horizontal thrust *TC* can be expressed by the tensile strength of the rock stratum *σ*max, and the relationship between *hm* and *L* can be obtained.

$$T\_{\mathbb{C}} = \frac{1}{2} \frac{h\_m^2}{h\_m - w} \sigma\_{\text{max}} \tag{8}$$

The large tensile strength and thickness of bedrock make the thrust *TC* between blocks larger. Therefore, the friction between collapses is large, and it is not easy to slip. When the horizontal thrust reaches the limit, the indirect contact surface of the blocks is broken, resulting in slip and collapse. This also shows that the length and thickness of the collapsed block, as well as the tensile strength of the block itself, affect the size of the space under it and even the duration of the cavity.

#### *2.2. CDEM Algorithm*

Before coal mining, each rock mass layer and coal seam are considered to be rock layers with continuous properties. When the coal seam is mined, it first causes the separation and collapse of the direct roof, and then, with the advance of the working face, the basic roof separation and collapse form a discrete block structure. After the coal seam is mined out, the bedrock fractures are developed, which shows the discreteness of cracking. The collapsed blocks fall into the goaf to form a discrete block structure. Therefore, the CDEM, based on the Lagrange equation, should be selected as the numerical simulation method of coal mining. This method can not only calculate the deformation state of coal and rock mass but also simulate the plasticity, damage, and the fracture process of coal and rock mass after applying the Mohr Coulomb criterion of brittle fracture on the virtual contact surface.

CDEM algorithm is an algorithm for solving the explicit solution of the dynamic equation of fractured elements, according to the theory of the generalized Lagrange system [17,18]. The algorithm is based on the strict generalized Lagrange equation, as shown in Formula (9), combined with the global dynamic relaxation algorithm, to iteratively solve the explicit solution of the equation.

$$\frac{d}{dt}(\frac{\partial L}{\partial \mathbf{x}\_i}) + (\frac{\partial L}{\partial y\_i}) = Q\_i \tag{9}$$

where *xi*, *yi* are generalized coordinates;

*L* is the energy of the Lagrangian system;

*Qi* is the work done by non-conservative force.

The basic structure of the CDEM algorithm is composed of two parts: blocks and interface. The block element is usually composed of quadrangles. It has only elastic deformation and no plastic deformation. The plastic deformation is concentrated on the virtual crack. For the direction of crack development and actual cracking, the theoretical direction is perpendicular to the direction of the maximum principal stress of the surrounding element, while the actual direction is the closest element boundary or along a diagonal direction [23]. Interface element is a constitutive model describing the relationship between adjacent blocks in the block discrete element, also known as contact surface. Three attributes, connecting two block elements and the normal vector of the contact surface, are used to represent the force between block elements. CDEM algorithm is the same grid shared by the Lagrangian element algorithm and the discrete element algorithm.

The core calculation formula of CDEM is shown in Formula (10) [28]:

$$M\{\ddot{u}(t)\} + \mathcal{C}\{\dot{u}(t)\} + K\{u(t)\} = \{Q(t)\} \tag{10}$$

where *M* is quality matrix;

.. *u*(*t*) is acceleration matrix;

*C* is damping matrix;

 . *u*(*t*) is velocity matrix;

*K is* stiffness matrix;

{*u*(*t*)} is displacement matrix;

When the cell is not embedded before the medium cracks, the node only contains unbalanced force (only elastic force and external force); when the medium is cracked or embedded, the node force changes, and the corresponding parameters of node velocity and acceleration motion state are generated. CDEM algorithm couples continuous and discontinuous media algorithms, and it realizes the perfect unity of the finite element, block discrete element, and particle discrete element [29]. GDEM numerical simulation platform, developed based on CDEM algorithm, is independently developed by the Institute of mechanics, Chinese Academy of Sciences. It is a constitutive model based on CDEM method. By means of CPU/GPU and behavior acceleration, it can not only simulate the elastic, plastic, damage, and fracture processes of geological bodies and artificial materials under multi field coupling but also develop unit constitutive and contact constitutive by using C++ or JavaScript, according to actual needs. The platform software contains 14 builtin block constitutive models and 12 built-in interface constitutive models [19]. In this study, the elastic damage fracture constitutive model is adopted, and the mechanical parameters of rock stratum include density, elastic modulus, Poisson's ratio, tensile strength, and internal friction angle. The brittle damage fracture constitutive model is applied on the virtual interface, and the mechanical parameters of the interface include normal stiffness, tangential stiffness, cohesion, internal friction angle, and tensile strength [30].

#### **3. Materials and Methods**

#### *3.1. Geological and Mining Conditions*

The study area is located in the southeast of Qinshui coalfield and the east of Jincheng mining area, belonging to the edge of Qinshui Basin. Qinshui coalfield has simple geological structure conditions, a flat coal seam dip angle, multi coal seam occurrence, and strong regularity, and most of the hydrogeological conditions are not complex. The geological conditions of the mine are the most superior mining conditions in China. Geographical location of the study area is shown in Figure 2.

**Figure 2.** Geographical location of the study area.

The coal bearing strata of coal mines in this area are the Permian Shanxi Formation and Carboniferous Taiyuan formation, with a total of 13–15 layers of coal, with a total thickness of nearly 13 m. According to the analysis of the geological data of Majiazhuang coal mine, the stratigraphic structure is mainly sandstone and mudstone, including limestone and siltstone. There are three stable minable coal seams, of which the buried depth of 3# coal seam is 31.4 m, and the average coal thickness is 2.9 m; the buried depth of 9# coal seam is 84.6 m, and the average coal thickness is 1.0 m; the burial depth of 15# coal seam is 120.4 m, and the average coal thickness is 4.3 m. Other coal seams are unstable coal seams or coal lines, which are not minable or are partially minable. The coal seam dip angle of the mine is less than 3◦, which is a nearly horizontal coal seam, belonging to the typical SBMCS occurrence condition on the southeast boundary of Qinshui coalfield, as shown in Figure 3.

**Figure 3.** Comprehensive geological histogram of Majiazhuang mine.

#### *3.2. Model and Mechanical Parameters*

This simulation takes the geological conditions in the southeast of Qinshui coalfield as the background, and it establishes a 1:1 numerical simulation model. A 2D model of 35 coal and rock strata is designed, including seven groups of rock strata with different properties and three coal seams. The design dimension of the model is 127.5 m high and 600 m wide. The mining length of the simulated working face is 300 m. The mining thickness is designed according to the actual average mining thickness of each layer of coal. According to the boundary conditions of the model, the left and right boundaries are fixed by means of zero velocity in the X-axis direction; in the Y-axis direction, the bottom boundary is completely fixed, and the upper boundary is free, as shown in Figure 4.

**Figure 4.** Schematic of model size and boundary conditions.

The initial in-situ stress of this simulation is calculated by using the linear elastic constitutive model, and the unbalance rate is set to 1 × 10−5. The mechanical parameters required for calculation are shown in Table 1.

**Table 1.** Calculation parameters of rock mechanics.

#### *3.3. Simulation Scheme and Steps*

This time simulates the mining process of three stable and minable coal layers in the southeast of Qinshui coalfield, studies the characteristics of bedrock failure and surface movement, and analyzes the characteristics of bedrock stress distribution. The numerical calculation is divided into two stages and three times of coal seam excavation. In the initial geostress calculation stage, the virtual mass calculation method is adopted. During the simulation, the coal and rock strata are discretized as a whole. The strength and stiffness of the contact surface are inherited from the cell strength, and the strength of the contact surface is set as 10% of the cell strength to obtain the initial geostress state under the action of gravity. The initial strain of coal is converted into the constitutive model of brittle rock after the initial strain of the coal-softening element is converted into the constitutive model of brittle rock. In the stage of the coal seam excavation solution, in order to better simulate the bedrock breaking process, the displacement and change gradient within the calculation range are zeroed. Then, calculate the coal seam excavation stage, excavate in three times, and set it as "*None*" model, according to the designed coordinate range of coal seam goaf, to realize coal seam excavation. After each coal seam excavation, the solution time step is set to 10,000. In order to master the variation law of surface settlement, monitoring points with an interval of 10 m are set on the surface. When simulating the excavation of different coal seams, the characteristics of bedrock fracture and surface failure are recorded.

The numerical model is corrected and adjusted by using surface monitoring data. During 3# coal mining, surface safety monitoring was conducted, and monitoring points were arranged along the advancing direction of the working face to monitor surface settlement.

#### **4. Results**

#### *4.1. Analysis of Bedrock Failure Characteristics*

This simulation adopts a 2D plane model. After the excavation calculation of three layers of coal, the failure characteristics of bedrock are obtained, as shown in Figure 5. Figure 5a shows that, when mining the upper coal seam, the coal seam floor near the coal pillar at the edge of the goaf is uplifted, and the floor is cracked and damaged. The key bedrock layer above the edge of the goaf breaks down to form a key block structure. The inclined broken block contacts the coal seam floor and the upper rock mass to form a triangular suspended space. The key blocks in the middle of the goaf are unstable and broken, and the "masonry beam" structure is formed by stacking in the goaf [19]. Affected by the failure of bedrock, there are the tensile zone, compression zone, and compressive uplift zone on the surface.

With the mining of the middle coal seam, the floor uplift of the upper coal seam gradually disappears. The mining thickness of the middle coal seam is small, and the change of bedrock failure characteristics is not obvious. Figure 5 shows the bedrock failure characteristics after mining of the lower coal seam. When the lower coal seam is mined, the coal seam floor is broken and lifted. The failure of the lower direct roof rock layer develops upward and affects the goaf and bedrock of the middle and upper coal seams, resulting in its instability and fracture again, and the floor rock layer of the middle coal seam breaks down, which makes the key rock layer above the middle coal seam unstable and collapse. A key block structure is formed at the edge of the goaf to support the overlying strata. The rock stratum breaks and collapses, and the cracks of the bedrock develop directly to the surface, resulting in great damage to the surface, forming huge staggered platforms, cracks, or collapse pits, with large widths and depths of damage. The development of fissures connects the goaf and the surface of the three layers of coal, resulting in the discharge of harmful gases. At the same time, it also forms a water diversion channel, and the water and sediment are fed into the working face and goaf. This not only pollutes the air and water near the mining area but also threatens the safety of the working face and surface water.

(**a**)

**Figure 5.** *Cont*.

(**d**)

**Figure 5.** Characteristics of bedrock and surface damage. (**a**) Bedrock failure in upper coal seam mining; (**b**) Bedrock failure of all coal mining; (**c**) Cloud map of stress distribution of bedrock; (**d**) Cloud map of bedrock settlement (left is vertical, right is horizontal).

According to the analysis of Figure 5c, the stress change of bedrock has obvious characteristics, and the failure of rock stratum starts from the surface of rock stratum. The bedrock above the edge of the goaf is mainly broken, cracked, and with tensile damage, while the middle part is mainly compressive and shear damage. The weak rock stratum in the bedrock is damaged by shear and tension, and the broken rock stratum belt is formed under the action of repeated mining. The damage boundary in bedrock decreases suddenly with the increase in mining depth. The damage boundary line presents a parabola shape and is symmetrically distributed on both sides of the goaf. The content shown in Figure 5d is the cloud map of bedrock settlement, and the bedrock settlement shows three gradients with the change of depth. The settlement of the lower bedrock is only affected by the mining of 15# coal seam, and the settlement is small, about 3.5–4.0 m; the subsidence of the middle bedrock is larger than that of the lower bedrock. Affected by the mining of the middle and lower coal seams, the subsidence caused by the bedrock has a superposition effect, and the subsidence is about 4.88–5.5 m; the upper bedrock has the largest settlement, which is affected by the mining of all coal seams, and the cumulative effect is more obvious. The maximum settlement is 8.183 m. From the gradient settlement caused by bedrock failure, the bedrock settlement caused by MCS mining is related to the mining thickness of the coal seam, and the influence is cumulative. The mining of the lower coal seam will lead to the repeated destruction of the rock mass of the upper bedrock and goaf, and the subsidence of the bedrock and surface bedrock, caused by the mining of the lower coal seam, is the most severe.

#### *4.2. Bedrock Failure Stress Analysis*

The stress balance of bedrock is broken during coal mining, resulting in bedrock instability, fracture, and fragmentation. After all the coal is mined, the stress distribution nephogram of bedrock failure shows left-right symmetry. The horizontal stress and vertical stress are selected for comparative analysis. It can be seen from Figure 6 that the left part is the cloud diagram of vertical stress distribution, and the right part is the cloud diagram of horizontal stress distribution. Before the instability and collapse of the key bedrock layer, it exists in the form of "cantilever beam" and suddenly breaks after reaching the ultimate bearing capacity. The key block is subjected to the pressure of the upper rock mass, and the vertical stress concentration occurs in the rock stratum above the coal pillar. The maximum downward vertical stress is about 9.79 MPa, and the maximum upward vertical stress is about 2.73 MPa. At the same time, the coal pillar is also affected by vertical stress, resulting in bending deformation and even coal explosion and rock explosion accidents. Before the bedrock is broken, it is stretched and bent, and there is a horizontal stress concentration area in the upper part of each rock stratum. After the rock stratum is broken, the stress on the upper and lower sides of the key block changes: the upper part is subjected to compressive stress, the lower part is subjected to tensile stress, and then, there is a horizontal stress concentration area. In the broken rock stratum, the interaction between rock blocks also produces stress concentration points. After the bedrock is broken, the maximum horizontal stress value is 4.9–6.1 MPa, which is symmetrically distributed left and right, in the opposite direction, and points to the middle of the basin.

**Figure 6.** Cloud map bedrock stress distribution (**left** is vertical, **right** is horizontal).

#### *4.3. Surface Deformation and Stress Characteristics*

There is a tensile zone, a compression zone, and a compressive fracture zone in the surface basin. The tensile zone is located at the edge of the surface subsidence basin, and the surface fractures are relatively developed, resulting in step fractures and tensile fractures [31]. The compression zone is located at the bottom of the edge of the subsidence basin, the cracks are closed, and the topsoil is compressed without crushing damage. The compression and fracture area is in the middle of the surface basin [32]. The topsoil is in the state of compression and fracture, resulting in compression and shear failure, and compression uplift occurs in some areas. These surface damage phenomena are consistent with the characteristics of surface damage caused by SBMCS [31,33].

In order to better analyze the characteristics of surface movement and deformation, in the process of numerical simulation, a monitoring point is set every 10 m on the surface, and a total of 61 surface monitoring points are arranged to monitor the settlement value and stress distribution of the surface during the simulation process, as shown in Figure 7. The characteristics of surface subsidence accord with the law of surface subsidence in the actual mining process. The surface presents a flat bottom moving basin, and the characteristics of surface subsidence caused by MCS mining are also obvious. According to Figure 7a, the three dense areas of the curve in the middle of the basin show the maximum subsidence curves of the three coal seam excavation stages, respectively. According to the 3# coal surface monitoring data collected, the maximum settlement is 2.622 m, and the maximum settlement obtained by numerical simulation is 3.0 m. The simulation state is ideal, and the difference of settlement is 0.378 m, which is a reasonable range. According to Figure 7b, the distribution law of surface horizontal stress is that the stress concentration fluctuates up and down at 10 KPa (downward), the stress in the middle of the basin (within 200~400 m) changes frequently, complexly, disorderly, and greatly, the stress value changes greatly, and the difference between the maximum value and the minimum value is greater than 26 KPa; however, the stress changes at the edge and outside of the basin are small [34].

According to Figure 7c, the overall variation characteristics of surface horizontal stress are obvious. Near the edge of the mobile basin (within 100 m), the horizontal stress changes significantly, and the horizontal stress changes slightly in the middle and outside of the mobile basin. The maximum stress reaches 100 KPa, and the minimum stress is about 77 KPa. Therefore, within a certain range of the surface, at the edge of the goaf, the horizontal stress changes actively, and the phenomena of tension and compression are obvious. The opening and closing movement of surface cracks in this area is obvious. Combined with the change of vertical stress, it is easy to produce staggered cracks and surface uplift. The stress change has obvious stratification phenomenon. The mining activity of the lower coal seam has an impact on the upper goaf and its bedrock. The stress increases due to the influence of lower coal seam mining.

**Figure 7.** *Cont*.

**Figure 7.** Characteristic analysis curve of surface movement and deformation. (**a**) Surface subsidence characteristic curve; (**b**) Characteristic of surface vertical stress distribution curve; (**c**) Characteristic of surface horizontal stress distribution curve.

#### **5. Discussion**

#### *5.1. The Bedrock Failure Process and Characteristics of MCS*

After each layer of coal is mined, a caving zone is formed. Then, the roof collapses directly, forming a collapse zone. The rock blocks falling into the goaf are relatively broken, and the particle size is small, which is mainly affected by the impact force of the mining. The volume of the broken rock and soil mass of the roof becomes larger, due to the influence of rock formation fragmentation. In the bedrock between different coal seams, when the key stratum reaches the failure limit, it presents a cantilever beam structure to support the broken rock layer above [20,35]. As the working face advances, a sudden brittle fracture occurs when the cantilever beam exceeds the ultimate bearing capacity. As the height increases, the filling space becomes narrow, and after the key layer above breaks down, an ordered block structure is formed as shown in Figure 8. The blocks are discrete and stacked on each other to form a masonry beam structure [12,36]. Therefore, the settlement of the broken bedrock layer, composed of the collapse zone and "masonry beam", needs many years to be gradually compacted, and the surface of the goaf can reach a relatively stable

state. The bedrock damage caused by MCS mining is more serious than that caused by single seam mining, and the damage recovery time may be longer.

**Figure 8.** Schematic of bedrock failure in mining of MCS.

The mining of SBMCS causes the damage state of bedrock and surface, which is directly related to the mining geological conditions. In Northwest China, the coal seam is thicker, and the bedrock is thinner. The surface of many coal mines is covered with thick loose aeolian sand layer, and the bedrock is thin [11,30]. The destruction of bedrock and surface caused by mining is more serious. However, the loose layer can repair the surface damage automatically and slowly. There are few buildings on the surface, which is suitable for mechanized high-intensity mining [37]. There are many differences between the multi coal seam mining in the Qinshui coalfield and the MCS mining in the Shendong mining area. Qinshui Coalfield does not have a thick loose layer on the surface, but it has thinner Loess and Sand Gravel. Therefore, the self-healing ability of the surface is poor. As a result, the failure characteristics and stress field distribution of the bedrock are different.

#### *5.2. Impact of Water on Mining Safety*

Underground mining engineering is a very dangerous thing if it is under the city or underwater. The water body has a great impact on the mining of coal mine. The mining causes the destruction of bedrock, and the crack is a good channel for the water body, which poses a great threat to the mining face. The mining of shallowly buried multi coal seams leads to bedrock collapse and relatively developed fractures, which leads to the direct connection between the surface and the working face and forms a water diversion channel. Surface water, groundwater, underground phreatic water, and other water bodies are affected by bedrock fissures and feed into the mining face, resulting in a flood in the working face.

The interaction between water body and coal rock mass is a complex physical and chemical process. Chemical actions such as dissolution, hydration, dissolution, and oxidation-reduction may occur between the water body and coal rock mass, which will not only change the composition and structure of coal rock mass but also affect the mechanical properties of coal rock mass [38] The infiltration of water into the rock and soil will also change its physical properties. The water in the fractured rock mass will produce lubrication at the fracture contact surface, reduce the friction resistance at the fracture surface, enhance the shear stress effect, and induce the shear movement along the fracture surface. Under the action of water softening and argillization, the filler in the fracture surface of rock mass begins to develop from solid state to plastic state to liquid state with the increase in water content, which reduces the mechanical properties of fractured rock mass, cohesion, and friction angle [39]. At the same time, the surface water is fed into the goaf, so that a large number of water bodies are stored in the goaf and overburden fractures, and the

harmful substances in the goaf will also pollute the water bodies. For the research on the influence of overburden on water, it is necessary to establish water monitoring (water level monitoring, water quality monitoring, pore water pressure monitoring) and relevant tests of various rock masses under different saturated water conditions. The research on the coupling effect of water body and rock mass needs further research.

#### *5.3. Numerical Simulation Method*

In terms of numerical simulation, the CDEM algorithm has many advantages in the process. The GDEM platform can organically combine the finite element and discrete element, and it has a good effect on the numerical simulation of bedrock and surface damage in coal mining. UEC software, based on the discrete element, and FLAC software, based on the finite element, can simulate and analyze underground mining. However, compared with the GDEM platform, based on the CDEM method, the GDEM platform is closer to engineering practice, has more constitutive models, and the mining process can be redeveloped. If there is a water-resisting layer in the goaf geological conditions, corresponding measures shall be taken to protect the water resisting layer. In mining MCS, attention should be paid to controlling the development height of bedrock fractures and keeping the water resisting layer from being damaged. We can regard the aquifer as an elastic template, build a balanced differential equation according to the theory of elastic mechanics and elastic foundation, and solve the differential equation, according to the relevant boundary conditions in the mining process. Furthermore, the mechanical analysis and calculation of multiple water barriers can be realized.

#### **6. Conclusions**

Through the research of this paper, it is found that the bedrock and surface damage caused by coal mining have strong regional characteristics. Different geological conditions have different characteristics. Bedrock damage and fracture development, caused by SBMCS mining, have regional characteristics. Different lithology, different water conditions, different burial depths, mining thickness conditions, and method of mining are all important factors affecting the failure characteristics of overlying rocks. The conclusions of this paper are as follows:

(4) The CDEM algorithm has many advantages in the process of numerical simulation. The GDEM platform can organically combine a finite element and a discrete element, and it has a good effect on the numerical simulation of bedrock damage in coal mining. In this study, it is necessary to further study and analyze the coupling effect between water body and bedrock. The coupling mechanical model of water body and broken rock mass is established to study the influence of the coupling model on the stress and settlement of bedrock and surface.

**Author Contributions:** Conceptualization, Y.Z. and H.C.; methodology, G.L.; software, G.L. and W.Z.; validation, Y.X., H.C. and G.L.; formal analysis, Y.Z.; investigation, W.Z.; resources, G.L.; data curation, G.L.; writing—original draft preparation, G.L.; writing—review and editing, Y.Z.; visualization, G.L.; supervision, W.Z.; project administration, W.Z.; funding acquisition, Y.Z. and W.Z. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the National Natural Science Foundation of China, grant number U1810203 and U21A20108, and This research was funded by Henan province science and technology tackling key project, China, grant number<PHONE_NUMBER>04. And This research was supported by the Fundamental Research Funds for the Universities of Henan Province, China, grant number NSFRF220432.

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

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The data used to support the findings of this study are included within the article.

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

#### **References**


*Communication* **Multiplexed Weak Waist-Enlarged Fiber Taper Curvature Sensor and Its Rapid Inline Fabrication**

**Duo Yi 1 , Lina Wang 1 , Youfu Geng 1,\* , Yu Du 1 , Xuejin Li 2 and Xueming Hong 1**

**Abstract:** This study proposes a multiplexed weak waist-enlarged fiber taper (WWFT) curvature sensor and its rapid fabrication method. Compared with other types of fiber taper, the proposed WWFT has no difference in appearance with the single mode fiber and has ultralow insertion loss. The fabrication of WWFT also does not need the repeated cleaving and splicing process, and thereby could be rapidly embedded into the inline sensing fiber without splicing point, which greatly enhances the sensor solidity. Owing to the ultralow insertion loss (as low as 0.15 dB), the WWFTbased interferometer is further used for multiplexed curvature sensing. The results show that the different curvatures can be individually detected by the multiplexed interferometers. Furthermore, it also shows that diverse responses for the curvature changes exist in two orthogonal directions, and the corresponding sensitivities are determined to be 79.1◦/m−1 and –48.0◦/m−1 respectively. This feature can be potentially applied for vector curvature sensing.

**Keywords:** curvature sensor; fiber interferometer; waist-enlarged fiber taper

**Citation:** Yi, D.; Wang, L.; Geng, Y.; Du, Y.; Li, X.; Hong, X. Multiplexed Weak Waist-Enlarged Fiber Taper Curvature Sensor and Its Rapid Inline Fabrication. *Sensors* **2021**, *21*, 6782. https://doi.org/10.3390/s21206782

Academic Editor: Jin Li

Received: 15 September 2021 Accepted: 11 October 2021 Published: 13 October 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/).

#### **1. Introduction**

Curvature sensing is essential in various fields such as smart robot [1], engineering structure [2] and medical image [3], etc. When compared with the traditional resistive curvature sensors, the optical fiber curvature sensor owns the advantages such as flexibility, small size, low cost, immunity to electromagnetic interference and easiness to form sensing network [4,5]. Nowadays, optical fiber curvature sensors based on intermodal interferometer have been widely reported [6–8]. Particularly in recent years, the use of the microstructured fiber as the transducing platform promotes the exploitation of various curvature sensors based on photonic crystal fiber [9], hollow-core fiber [7], few-mode fiber [10], four-core fiber [11] or Fabry-Perot [12], etc. However, for these intermodal interferometers, post-processing techniques such as fiber tapering [13], side-polishing [14] and cladding etching [15] are always used to modify the fiber structure and convert part of fundamental mode into the high-order mode in order to construct an optical arm. It inevitably results in a large insertion loss and structure fragileness, which makes the essentially multiplexed curvature sensing difficult to be achieved with these techniques.

In our previous work [16], a waist-enlarged fiber taper sensor has been proposed for high-temperature sensing. It has a remarkable and distinguishable fiber structure, a strongly waist-enlarged taper, to excite high-order fiber modes. This type of fiber taper can be named as a strong waist-enlarged fiber taper (SWFT) and has been comprehensively used to construct fiber sensors for the sensing applications such as temperature [16], strain [17], refractive index [18,19] and humidity [20]. It owns the advantages of robust structure, easy fabrication and flexible mode excitation, etc. However, it needs repeated cleaving and splicing during the fabrication process, and the structural health of the sensor still needs

to be improved due to the embedded splicing points in the fiber. Besides, similar to the intermodal interferometers mentioned in the previous paragraph, the fabricated SWFT also faces the problem of large transmission loss, which results in a poor signal-noise ratio, but in contrast, presents a moderate spectral dynamic range over 10 dB [16]. the fabricated SWFT also faces the problem of large transmission loss, which results in a poor signal-noise ratio, but in contrast, presents a moderate spectral dynamic range over 10 dB [16]. In this study, a weak waist-enlarged fiber taper (WWFT) and its rapid inline fabrica-

bust structure, easy fabrication and flexible mode excitation, etc. However, it needs repeated cleaving and splicing during the fabrication process, and the structural health of the sensor still needs to be improved due to the embedded splicing points in the fiber. Besides, similar to the intermodal interferometers mentioned in the previous paragraph,

*Sensors* **2021**, *21*, x FOR PEER REVIEW 2 of 10

In this study, a weak waist-enlarged fiber taper (WWFT) and its rapid inline fabrication method are further proposed, and the application for multiplexed curvature measurement is studied. The WWFT could be fabricated rapidly with a weak electric arc discharge released on the inline sensing fiber, and the repeated cleaving-splicing process is no longer needed. It has no fusion points and no obvious structural deformation in appearance compared with general single mode fiber (SMF). Therefore, the insertion loss of the WWFTbased interferometer is greatly reduced, which makes it possible to realize multiplex curvature sensing with phase-demodulation scheme. Unfortunately, the WWFT device presents a relatively low spectral dynamic range. Finally, the experimental results of curvature measurement demonstrate that the proposed sensor has diverse responses for the curvature changes in two orthogonal directions, showing its great potential applications in vector curvature sensing. tion method are further proposed, and the application for multiplexed curvature measurement is studied. The WWFT could be fabricated rapidly with a weak electric arc discharge released on the inline sensing fiber, and the repeated cleaving-splicing process is no longer needed. It has no fusion points and no obvious structural deformation in appearance compared with general single mode fiber (SMF). Therefore, the insertion loss of the WWFT-based interferometer is greatly reduced, which makes it possible to realize multiplex curvature sensing with phase-demodulation scheme. Unfortunately, the WWFT device presents a relatively low spectral dynamic range. Finally, the experimental results of curvature measurement demonstrate that the proposed sensor has diverse responses for the curvature changes in two orthogonal directions, showing its great potential applications in vector curvature sensing.

#### **2. Sensor Fabrication 2. Sensor Fabrication**  Figure 1a shows the schematic of WWFT and the constructed intermodal interferom-

Figure 1a shows the schematic of WWFT and the constructed intermodal interferometer, which has no difference in appearance with the SMF. Figure 1b shows the specific inline fabrication process of WWFT. It is fabricated by using the manual splicing mode of the fiber splicer (Fujikura Ltd., FSM 60S, Tokyo, Japan) in our experiment. Firstly, two segments of SMFs with polymer coating stripped and end cleaved are prepared, and they are separately put onto the motor stages of the fiber holder. The motors are impelled and the fiber end faces are aligned, then the splicing process is temporarily suspended. In this step, the motors are adjusted to avoid fiber bending during arc discharge. Secondly, an SMF whose polymer coating has been stripped is used to replace the aligned fibers on the motor stages, and arc discharge is then released when the two motor stages are pushed forward with a designed overlap distance, defined as *L*vp. In this study, three *L*vp values of 10, 20 and 30 µm are selected to fabricate the WWFT in the experiment for comparison. Using this method, the WWFT is rapidly embedded into the sensing SMF fiber. The arc discharge intensity and duration time are set as 'standard mode' and 200 ms respectively. Note that, if the motor can be automatically aligned, the aligned process with two cleaved fibers in the first step is no longer needed and the duration of fabrication process could be shortened to just a few seconds. Therefore, the WWFT is possible to be massively fabricated by a specially designed arc releasing electrical and mechanical structure for rapidly arranging a large-scale sensing network. Figure 1c shows the side views of the WWFT. It can be seen that, different from the SWFT, it has no difference in appearance with SMF, and no splicing point is embedded in the sensing fiber, which gives a strong physical structure and ultralow insertion. eter, which has no difference in appearance with the SMF. Figure 1b shows the specific inline fabrication process of WWFT. It is fabricated by using the manual splicing mode of the fiber splicer (Fujikura Ltd., FSM 60S, Tokyo, Japan) in our experiment. Firstly, two segments of SMFs with polymer coating stripped and end cleaved are prepared, and they are separately put onto the motor stages of the fiber holder. The motors are impelled and the fiber end faces are aligned, then the splicing process is temporarily suspended. In this step, the motors are adjusted to avoid fiber bending during arc discharge. Secondly, an SMF whose polymer coating has been stripped is used to replace the aligned fibers on the motor stages, and arc discharge is then released when the two motor stages are pushed forward with a designed overlap distance, defined as *L*vp. In this study, three *L*vp values of 10, 20 and 30 μm are selected to fabricate the WWFT in the experiment for comparison. Using this method, the WWFT is rapidly embedded into the sensing SMF fiber. The arc discharge intensity and duration time are set as 'standard mode' and 200 ms respectively. Note that, if the motor can be automatically aligned, the aligned process with two cleaved fibers in the first step is no longer needed and the duration of fabrication process could be shortened to just a few seconds. Therefore, the WWFT is possible to be massively fabricated by a specially designed arc releasing electrical and mechanical structure for rapidly arranging a large-scale sensing network. Figure 1c shows the side views of the WWFT. It can be seen that, different from the SWFT, it has no difference in appearance with SMF, and no splicing point is embedded in the sensing fiber, which gives a strong physical structure and ultralow insertion.

**Figure 1.** Schematic of WWFT and its fabrication: (**a**) Schematic of WWFT and the constructed intermodal interferometer; (**b**) Inline fabrication of WWFT; (**c**) XY side views of the fiber taper with overlap of 20 µm. The red circle represents the fiber taper location.

Next, two WWFTs are sequentially embedded into the sensing SMF to construct an intermodal interferometer, as shown in Figure 2. For the experimental test, the WWFTbased interferometer is connected with ASE light source (Golight Tech., ASE-C&L-10, Shenzhen, China) and CCD spectral module (Bayspec Inc., FBGA-1525-1610, Fremont, CA, America) for transmission spectrum acquisitions. The acquired spectra are then analyzed by using a LabVIEW program for Fast Fourier Transform (FFT), and the peak amplitudes/phases of the spectra after FFT are individually discussed. The schematic diagram of the experimental setup is shown in Figure 2. intermodal interferometer, as shown in Figure 2. For the experimental test, the WWFTbased interferometer is connected with ASE light source (Golight Tech., ASE-C&L-10, Shenzhen, China) and CCD spectral module (Bayspec Inc., FBGA-1525-1610, Fremont, CA, merica) for transmission spectrum acquisitions. The acquired spectra are then analyzed by using a LabVIEW program for Fast Fourier Transform (FFT), and the peak amplitudes/phases of the spectra after FFT are individually discussed. The schematic diagram of the experimental setup is shown in Figure 2.

Next, two WWFTs are sequentially embedded into the sensing SMF to construct an

**Figure 1.** Schematic of WWFT and its fabrication: (**a**) Schematic of WWFT and the constructed intermodal interferometer; (**b**) Inline fabrication of WWFT; (**c**) XY side views of the fiber taper with

*Sensors* **2021**, *21*, x FOR PEER REVIEW 3 of 10

overlap of 20 μm. The red circle represents the fiber taper location.

**Figure 2.** Experiment setup of WWFT-based interferometer. ASE: amplified spontaneous emission, CCD: charge-coupled device, WWFT: weak waist-enlarged fiber taper. **Figure 2.** Experiment setup of WWFT-based interferometer. ASE: amplified spontaneous emission, CCD: charge-coupled device, WWFT: weak waist-enlarged fiber taper.

#### **3. Spectral Characterization of WWFT-Based Interferometer 3. Spectral Characterization of WWFT-Based Interferometer**

For the following analyses, sensors with different overlap values *L*vp (10, 20 and 30 μm) are experimentally fabricated for comparison. Besides, various MZ interferometer samples with different interferometer lengths of *L* are also prepared*.* Figure 3 shows the transmission spectra of the WWFT-based interferometers with overlaps of 10 μm, 20 μm and 30 μm, respectively. We can observe that the induced insertion losses could be as low as 0.15 dB when *L*vp = 20 μm. It increases slightly with the overlap distance enlarged since more high-order modes are excited and experience large transmission loss. From the figure, it also can be seen that if the overlap distance remains the same, the fringe contrast decreases as the interferometer length elongates. For example, when *L*vp = 20 μm, the maximum fringe contrasts are 0.68 dB, 0.37 dB, 0.12 dB, respectively, for the interferometer lengths of 30.8 mm, 48.7 mm and 70.1 mm. Actually, when the interferometer length increases, the optical power ratio between sensing arm and reference arm of interferometer gradually decreases since the high-order mode suffers a higher transmission loss. Finally, For the following analyses, sensors with different overlap values *L*vp (10, 20 and 30 µm) are experimentally fabricated for comparison. Besides, various MZ interferometer samples with different interferometer lengths of *L* are also prepared. Figure 3 shows the transmission spectra of the WWFT-based interferometers with overlaps of 10 µm, 20 µm and 30 µm, respectively. We can observe that the induced insertion losses could be as low as 0.15 dB when *L*vp = 20 µm. It increases slightly with the overlap distance enlarged since more high-order modes are excited and experience large transmission loss. From the figure, it also can be seen that if the overlap distance remains the same, the fringe contrast decreases as the interferometer length elongates. For example, when *L*vp = 20 µm, the maximum fringe contrasts are 0.68 dB, 0.37 dB, 0.12 dB, respectively, for the interferometer lengths of 30.8 mm, 48.7 mm and 70.1 mm. Actually, when the interferometer length increases, the optical power ratio between sensing arm and reference arm of interferometer gradually decreases since the high-order mode suffers a higher transmission loss. Finally, it results in the deterioration of fringe contrast. *Sensors* **2021**, *21*, x FOR PEER REVIEW 4 of 10 optical path difference could be also expressed as = ∆, where ∆*n* is the difference of effective refractive index between LP01 and the excited high-order mode. Therefore, associating with the two different expressions of above, the ∆*n* could be determined in experiment, and the value is 2.59 × 10−3, which is very close to the theoretical value of 2.73 × 10−3 between LP01 and LP11 mode using beam propagation method. Hence we can infer that the excited high-order mode is LP11 mode. The electric field distributions of LP01 and LP11 mode are displayed as the inset of Figure 4a.

**Figure 3.** Transmission spectra of WWFT-based interferometers with different overlap values of (**a**) 10 μm, (**b**) 20 μm and (**c**) 30 μm, respectively. **Figure 3.** Transmission spectra of WWFT-based interferometers with different overlap values of (**a**) 10 µm, (**b**) 20 µm and (**c**) 30 µm, respectively.

**Figure 4.** (**a**) Spectra after FFT for the WWFT-based interferometers with different interferometer lengths, inset: Electric field distributions of LP01 and LP11 mode; (**b**) Optical path difference versus interferometer length with overlaps of 10 μm,

> Owing to the low insertion loss of WWFT, the proposed sensor can be applied for curvature multiplex sensing. Figure 5a shows the schematic of the curvature measurement for two cascaded WWFT-based interferometers. The fiber sensor 1# and 2# are with interference lengths *L* = 41.7 mm and 69.5 mm, respectively. Besides, the overlap distances are set as 20 μm and 30 μm, respectively, for the two sensors. For the two cascaded interferometers, the total insertion loss is measured as approximately 2.5 dB, and it can be further reduced with a smaller overlap value. Besides, the amount of multiplexed sensors depends on the measurable spectral range. According to the equation δ = /2∆ in

**4. Results and Analyses**

20 μm and 30 μm, respectively.

(**c**) 30 μm, respectively.

Since WWFT is with non-adiabatic characteristics, part of the energy of LP01 is leaked and coupled into the high-order modes [21]. Next, in order to evaluate the specific mode order excited by WWFT, the Fast Fourier Transform (FFT) is applied to the transmission spectra in Figure 3. Herein, Figure 3b is selected, and another two sensors with interferometer lengths of 41.2 mm and 59.1 mm are supplemented for the following analyses. The results are shown in Figure 4a. For each interferometer with a certain physical length *L*, the optical path difference (OPD) *δ* could be calculated through the relationship of *δ* = *cKLm*/2∆*ν* [22], where *KLm* is the index of peak amplitude location in the FFT spectra, and ∆ν is the whole frequency range of CCD spectral module. It should be noted that *Klm* is a parameter with no dimension, and it varies with the frequency range of CCD spectral module. Hence, the calculated OPD itself is independent with the CCD module. Figure 4b illustrates the calculated OPD versus interferometer length *L* of the WWFT-based interferometer with different overlap distances of 10, 20 and 30 µm. It is well known that the optical path difference could be also expressed as *δ* = ∆*nL*, where ∆*n* is the difference of effective refractive index between LP01 and the excited high-order mode. Therefore, associating with the two different expressions of *δ* above, the ∆*n* could be determined in experiment, and the value is 2.59 × 10−3 , which is very close to the theoretical value of 2.73 × 10−3 between LP01 and LP11 mode using beam propagation method. Hence we can infer that the excited high-order mode is LP11 mode. The electric field distributions of LP01 and LP11 mode are displayed as the inset of Figure 4a. **Figure 3.** Transmission spectra of WWFT-based interferometers with different overlap values of (**a**) 10 μm, (**b**) 20 μm and

optical path difference could be also expressed as = ∆, where ∆*n* is the difference of effective refractive index between LP01 and the excited high-order mode. Therefore, associating with the two different expressions of above, the ∆*n* could be determined in experiment, and the value is 2.59 × 10−3, which is very close to the theoretical value of 2.73 × 10−3 between LP01 and LP11 mode using beam propagation method. Hence we can infer that the excited high-order mode is LP11 mode. The electric field distributions of LP01 and

*Sensors* **2021**, *21*, x FOR PEER REVIEW 4 of 10

LP11 mode are displayed as the inset of Figure 4a.

**Figure 4.** (**a**) Spectra after FFT for the WWFT-based interferometers with different interferometer lengths, inset: Electric field distributions of LP01 and LP11 mode; (**b**) Optical path difference versus interferometer length with overlaps of 10 μm, 20 μm and 30 μm, respectively. **Figure 4.** (**a**) Spectra after FFT for the WWFT-based interferometers with different interferometer lengths, inset: Electric field distributions of LP01 and LP11 mode; (**b**) Optical path difference versus interferometer length with overlaps of 10 µm, 20 µm and 30 µm, respectively.

#### **4. Results and Analyses 4. Results and Analyses**

Owing to the low insertion loss of WWFT, the proposed sensor can be applied for curvature multiplex sensing. Figure 5a shows the schematic of the curvature measurement for two cascaded WWFT-based interferometers. The fiber sensor 1# and 2# are with interference lengths *L* = 41.7 mm and 69.5 mm, respectively. Besides, the overlap distances are set as 20 μm and 30 μm, respectively, for the two sensors. For the two cascaded interferometers, the total insertion loss is measured as approximately 2.5 dB, and it can be further reduced with a smaller overlap value. Besides, the amount of multiplexed sensors depends on the measurable spectral range. According to the equation δ = /2∆ in Owing to the low insertion loss of WWFT, the proposed sensor can be applied for curvature multiplex sensing. Figure 5a shows the schematic of the curvature measurement for two cascaded WWFT-based interferometers. The fiber sensor 1# and 2# are with interference lengths *L* = 41.7 mm and 69.5 mm, respectively. Besides, the overlap distances are set as 20 µm and 30 µm, respectively, for the two sensors. For the two cascaded interferometers, the total insertion loss is measured as approximately 2.5 dB, and it can be further reduced with a smaller overlap value. Besides, the amount of multiplexed sensors depends on the measurable spectral range. According to the equation δ = *cKL*/2∆*ν* in the last paragraph, a larger spectral range results in a shorter optical path difference, and further, a higher resolution in the spectra after FFT, hence more sensors can be multiplexed. During the test, sensor 1# is placed between two supporting points with a separated distance of 115 mm, while sensor 2# keeps still. A one-dimensional translation stage with circular metal head is placed at the center of the sensor 1#. With the micrometer screw thread pushed forward, a fiber arc shape is formed. The curvature can be calculated by the expression of ρ = 1/*R* = *8d*/(4*d* 2 + *D*2 ) [23], where *R* is the curvature radius of the fiber arc induced by pushing the one-dimension fiber translation stage, *D* is the distance of two

supporting points, and *d* is the bending displacement. In our case, *d D*, therefore, the curvature is simplified as ρ = 1/*R* = *8d*/*D*2 . arc induced by pushing the one-dimension fiber translation stage, *D* is the distance of two supporting points, and *d* is the bending displacement. In our case, *d << D*, therefore, the curvature is simplified as ρ *=* 1*/R = 8d/D*2.

the last paragraph, a larger spectral range results in a shorter optical path difference, and further, a higher resolution in the spectra after FFT, hence more sensors can be multiplexed. During the test, sensor 1# is placed between two supporting points with a separated distance of 115 mm, while sensor 2# keeps still. A one-dimensional translation stage with circular metal head is placed at the center of the sensor 1#. With the micrometer screw thread pushed forward, a fiber arc shape is formed. The curvature can be calculated by the expression of ρ *=* 1*/R = 8d/*(4*d*2 *+ D*2) [23], where *R* is the curvature radius of the fiber

*Sensors* **2021**, *21*, x FOR PEER REVIEW 5 of 10

**Figure 5.** (**a**) Schematic of curvature test for the two cascaded WWFT-based interferometers; (**b**) Spectra and (**c**) Spectra after FFT for the WWFT-based interferometer 1# and 2# at 0° direction and 90° direction, respectively. **Figure 5.** (**a**) Schematic of curvature test for the two cascaded WWFT-based interferometers; (**b**) Spectra and (**c**) Spectra after FFT for the WWFT-based interferometer 1# and 2# at 0◦ direction and 90◦ direction, respectively.

Based on the experimental setup in Figure 5a, different external curvatures are applied on the sensor 1#, and the corresponding transmission spectra at the two orthogonal directions (0° and 90°, schematic of fiber rotation angle is indicated in the inset of Figure 5a) are shown in Figure 5b. It should be noted that the axis Y in Figure 5b represents the relative intensity; it is obtained by normalizing the transmission spectrum amplitude with the input power, and then the gravity center of the spectral intensity shifts to zero to eliminate the DC component, hence both positive and negative values are shown. It can be Based on the experimental setup in Figure 5a, different external curvatures are applied on the sensor 1#, and the corresponding transmission spectra at the two orthogonal directions (0◦ and 90◦ , schematic of fiber rotation angle is indicated in the inset of Figure 5a) are shown in Figure 5b. It should be noted that the axis Y in Figure 5b represents the relative intensity; it is obtained by normalizing the transmission spectrum amplitude with the input power, and then the gravity center of the spectral intensity shifts to zero to eliminate the DC component, hence both positive and negative values are shown. It can be seen that the wavelength shifts in Figure 5b are insignificant and irregular, which signifies that the curvature information cannot be directly demodulated. Figure 5c shows the corresponding spectra after FFT; it can be seen that the peak amplitudes of sensor 2# almost keep still, while the peak amplitudes of sensor 1# become lower as the curvature increases in both 0◦ and 90◦ directions. This phenomenon is reasonable since the LP11 mode is more sensitive to the bending effect, which leads to an enlarged transmission loss. Hence, the fringe contrast at the second WWFT point starts to deteriorate and further results in the reduction of peak amplitude.

Next, the continuous data acquisition is conducted, and the dynamic results are shown in Figure 6. The peak phases of spectra after FFT are continuously acquired as curvature varies between 0.48~0.94 m−1 (0◦ direction) ◦and 0.48~0.76 m−1 (90◦ direction). Figure 6a shows the phase shift when the WWFT-based fiber is placed at the initial position (defined as 0◦ direction). It can be seen that, when the curvature increases, the phase shift of sensor 1# is enlarged synchronously. Meanwhile, sensor 2# keeps almost unchanged, indicating that the curvature detections are independent between two sensors. Next, when the fiber sensor 1# is rotated to 90◦ , as shown in Figure 6b, we can see that the phase of sensor 1# shifts to the negative direction and presents different phase shifts when compared with those in 0◦ direction. This phenomenon signifies that the proposed sensor has diverse responses in two orthogonal directions, which can be further potentially applied for vector curvature sensing. Besides, for the 90◦ direction, we can also observe the cross-talk for sensor 2#. Theoretically speaking, no phase shift should occur for sensor 2# since no curvature is applied. This experimental phenomenon may be explained by the resolution ability of spectra after FFT, which is determined by the measurable spectral range of the CCD spectral module. As shown in Figure 5c, the peak amplitude of sensor 2# is easily affected by that of sensor 1# when the spectrum resolution ability is limited, and the cross-talk becomes more significant when the sensitivity decreases. Finally, the sensitivities in both 0◦ and 90◦ directions are shown in Figure 6c. The curvature sensitivities are determined to be 79.1◦/m−1 and −48.0◦/m−1 , respectively. The different curvature sensitivities are mainly due to the anisotropic distribution of the excited high-order mode (LP11) in different axes. It should be explained that most of previous works relating to optical fiber curvature sensors are based on wavelength demodulation in unit of nanometer (nm) or amplitude demodulation in unit of decibel (dB). While in this study, the curvature sensing is realized based on phase demodulation in unit of degree (◦ ). It should be noted that the wavelength and phase are two independent parameters which are inconvenient to convert to each other in our demodulation system. Hence, it is difficult to compare the curvature sensitivity performance presented in this study with the previously reported studies. Here, the sensitivity performance comparison is not discussed.

Finally, the temperature cross-sensitivity of the presented WWFT-based interferometer is discussed. The sensor is heated from 22.0 ◦C to 37.5 ◦C. Figure 7a shows the transmission spectra of the sensor with interferometer length of 41.2 mm and overlap value of 20 µm. We can hardly observe the wavelength shift with wavelength demodulation scheme. Then, based on phase demodulation scheme, the phase shift versus temperature is indicated in Figure 7b, and the temperature cross-sensitivity is determined to be 1.37 ◦/m−1 . When compared with the curvature sensitivities mentioned above, the temperature crosssensitivity is insignificant. Besides, the phase fluctuation is monitored for 4 min at 22 ◦C, and the result is shown in Figure 7c. The recorded standard deviation is determined to be 0.071◦ . Here, we use three times the standard deviation as the phase resolution, i.e., 0.213◦ . Hence, the curvature resolutions of the proposed sensor are determined to be around 0.213/79.1 = 2.69 × 10−3 m−1 at 0◦ direction and 0.213/48.0 = 4.44 × 10−3 m−1 at 90◦ direction.

**Figure 6.** Phase responses of the WWFT-based interferometer 1# and sensor 2# in (**a**) 0° direction and (**b**) 90° direction; (**c**) Phase shift versus curvature in 0° and 90° direction. **Figure 6.** Phase responses of the WWFT-based interferometer 1# and sensor 2# in (**a**) 0◦ direction and (**b**) 90◦ direction; (**c**) Phase shift versus curvature in 0◦ and 90◦ direction.

Finally, the temperature cross-sensitivity of the presented WWFT-based interferometer is discussed. The sensor is heated from 22.0 °C to 37.5 °C. Figure 7a shows the transmission spectra of the sensor with interferometer length of 41.2 mm and overlap value of 20 μm. We can hardly observe the wavelength shift with wavelength demodulation scheme. Then, based on phase demodulation scheme, the phase shift versus temperature is indicated in Figure 7b, and the temperature cross-sensitivity is determined to be 1.37 °/m−1. When compared with the curvature sensitivities mentioned above, the temperature cross-sensitivity is insignificant. Besides, the phase fluctuation is monitored for 4 min at 22 °C, and the result is shown in Figure 7c. The recorded standard deviation is determined to be 0.071°. Here, we use three times the standard deviation as the phase resolution, i.e., 0.213°. Hence, the curvature resolutions of the proposed sensor are determined to be around 0.213/79.1 = 2.69 × 10−3 m−1 at 0° direction and 0.213/48.0 = 4.44 × 10−3 m−1 at 90°

direction.

**Figure 7.** Temperature cross-sensitivity of WWFT-based interferometer: (**a**) Transmission spectra as temperature increases from 22 °C to 37.5 °C; (**b**) Phase shift versus temperature for the sensor with *L* = 41.2 mm and *L*vp = 20 μm; (**c**) Continuous monitoring of phase fluctuation at 22.0 °C. **Figure 7.** Temperature cross-sensitivity of WWFT-based interferometer: (**a**) Transmission spectra as temperature increases from 22 ◦C to 37.5 ◦C; (**b**) Phase shift versus temperature for the sensor with *L* = 41.2 mm and *L*vp = 20 µm; (**c**) Continuous monitoring of phase fluctuation at 22.0 ◦C.

#### **5. Conclusions 5. Conclusions**

This study demonstrates a multiplexed WWFT-based curvature sensor and its rapid inline fabrication. When compared with the other types of fiber taper, the proposed WWFT has no difference in appearance with the SMF, which greatly reduces the insertion loss. The WWFTs could be rapidly embedded into the inline sensing fiber without the repeated cleaving-splicing process and splicing point, which enhances the solidity of sensor. Owing to the low transmission loss which could be as low as 0.15 dB, the WWFTbased interferometer can be further multiplexed based on phase-demodulation scheme. The experimental results verified that curvature variation can be dynamically detected for two cascaded WWFT-based interferometers. Besides, the proposed sensor shows diverse responses for the curvature changes in two orthogonal directions, and the sensitivities are This study demonstrates a multiplexed WWFT-based curvature sensor and its rapid inline fabrication. When compared with the other types of fiber taper, the proposed WWFT has no difference in appearance with the SMF, which greatly reduces the insertion loss. The WWFTs could be rapidly embedded into the inline sensing fiber without the repeated cleaving-splicing process and splicing point, which enhances the solidity of sensor. Owing to the low transmission loss which could be as low as 0.15 dB, the WWFT-based interferometer can be further multiplexed based on phase-demodulation scheme. The experimental results verified that curvature variation can be dynamically detected for two cascaded WWFT-based interferometers. Besides, the proposed sensor shows diverse responses for the curvature changes in two orthogonal directions, and the sensitivities are

determined to be 79.1°/m−1 and −48.0°/m−1. Correspondingly, the curvature resolutions in

determined to be 79.1◦/m−1 and −48.0◦/m−1 . Correspondingly, the curvature resolutions in two orthogonal directions are estimated as 2.69 × 10−3 m−1 and 4.44 × 10−3 m−1 , respectively. This result confirms the feasibility for vector curvature sensing applications in the future.

**Author Contributions:** Conceptualization, Y.G. and D.Y.; methodology, Y.G.; software, Y.G.; validation, L.W., X.L. and X.H.; formal analysis, D.Y.; investigation, L.W.; resources, Y.D.; data curation, D.Y.; writing—original draft preparation, D.Y. and Y.G.; writing—review and editing, D.Y and Y.G. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was supported, in part, by Stabilization Support Program for Higher Education Institutions of Shenzhen, grant number<PHONE_NUMBER>2156001, and, in part, by National Natural Science Foundation of China, grant number 51808347.

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

**Informed Consent Statement:** Not applicable.

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

#### **References**


**Reflecting on PASUC Heritage Initiatives through Time, Positionality, and Place**

#### **Scott Hutson 1,\*, Céline Lamb 1 , Daniel Vallejo-Cáliz 1 and Jacob Welch 2**

Received: 17 March 2020; Accepted: 8 April 2020; Published: 14 April 2020

**Abstract:** This paper reports on heritage initiatives associated with a 12-year-long archaeology project in Yucatan, Mexico. Our work has involved both surprises and setbacks and in the spirit of adding to the repository of useful knowledge, we present these in a frank and transparent manner. Our findings are significant for a number of reasons. First, we show that the possibilities available to a heritage project facilitated by archaeologists depend not just on the form and focus of other stakeholders, but on the gender, sexuality, and class position of the archaeologists. Second, we provide a ground-level view of what approaches work well and which do not in terms of identifying aspects of cultural heritage that are relevant to a broad swath of stakeholders. Finally, we discuss ways in which heritage projects can overcome constraints to expanding community collaboration.

**Keywords:** archaeology; community museums; gender and sexuality; Maya