Abstract

Stress wave which is caused either by an explosion in a borehole or by an accidental explosion in a tunnel is supposed to be considered under certain circumstances when it propagates through the surrounding rock masses which contain holes in cylindrical form. Studying the ground motion induced by the cylindrical wave propagation is of practical significance for underground rock engineering and underground energy exploitation. The current study presents a numerical study on the ground motion caused by cylindrical P-wave propagation across a rock mass with a structural plane using a discrete element numerical method, UDEC. Firstly, the accuracy and validity of the cylindrical wave propagation simulation in UDEC and of the induced ground vibration are confirmed by comparison with the theoretical results for a special case that there is no structural plane in a rock mass. Secondly, cylindrical wave propagation across a rock mass with a structural plane is simulated, and then, the particle velocity on the ground surface is subsequently obtained. Finally, parametric researches are carried out on the influence of the monitoring point’s position, the structural plane stiffness, and the frequency of incident wave on the peak particle velocities (PPVs) of the ground vibrations.

1. Introduction

Dynamic loads caused by engineering blasting usually spread in rock mass in the form of stress waves. Stress wave propagation in rock masses is a crucial topic in rock dynamics and rock engineering. Especially, when stress waves arrive at the ground, vibrations will be caused, which may threaten the safety and stability of buildings and structures on the ground. Therefore, it is of practical significance to pay attention to the ground motion caused by stress waves in the construction and design of underground engineering. An important index used to measure the vibration intensity is the peak value of each physical quantity of the medium particle motion, such as the peak value of displacement, velocity, and acceleration [1]. Generally, the peak particle velocity (PPV) is the most widely used one to quantify the ground vibration [25]. At present, in many countries′ construction and design codes, PPV is one of the important criterions in assessing the structural responses under dynamic loads.

Chronically, blast-induced waves have generally been treated as plane waves in researches or practices, which is reasonable when the source is relatively far. However, the assumption of plane waves is no longer suitable in the following two cases: (1) when blasting operations are carried out near by the existing structures and (2) when the rock damage near the blasting site needs to be monitored and controlled. In those two cases, the curvature of the near-field wave fronts cannot be ignored. Especially, in practical explosion engineering, if the cavity or borehole is long enough for the dynamic loads, the stress wave propagation in the elastic vibration zone near blasting source could not be accurately treated as a plane wave [6, 7]. When a beam of cylindrical wave propagates though a medium, the wave front takes place a form of a cylinder, and the velocity of each particle is constant in the radial direction on a wave front. In addition, cylindrical waves have an obvious attenuation due to geometrical spreading, which is significantly different from plane waves [1]. In recent years, certain researches on cylindrical wave propagation have been carried out. Such as, Li et al. [3] carried out a preliminary estimation induced by an underground explosion according to assuming the blast loading as a cylindrical P-wave by using the time-domain recursive method (TDRM) [8, 9]. This approach was then improved to predict underground tunnel dynamic response [4]. In view of the circular tunnel under cylindrical P-wave action, Lu et al. [10] theoretically studied the particle velocity response of surrounding rocks. However, the existing theoretical researches on cylindrical wave focus mainly on the induced ground motion as wave propagates in an unjointed rock mass. Few researches are involved in the induced ground motion caused by the wave propagation in a rock mass with discontinuities.

Discontinuities are found in rock mass, such as natural fractures, joints, and faults, which affect the propagation of stress waves. Reflection and transmission of stress wave will take place on a rock structural plane. So far, a lot of theoretical analyses have been carried out to study stress wave propagation across jointed rock masses. Several theoretical methods of calculation have been developed, such as the displacement discontinuity method (DDM) [1114], the equivalent modulus method (EMM) [15], and the virtual wave source method (VWSM) [16, 17]. However, most of the researches focus either on the case that the wave front is a plane or on the case that the stress wave lacks of the curved wave front, such as cylindrical wave and spherical wave. When a cylindrical wave propagates in a rock mass with a natural structural plane, multiple reflections will occur between the structural plane and the ground surface. The ground vibration will be superposed by the multiple reflected waves. Obviously, the ground motion is quite complicated; hence, analytical solutions are very difficult to obtain. What is worse, the experimental study seems difficult to conduct due to dangers, complexities, and large expense. Therefore, numerical simulation seems to be a convenient alternative and acceptable method.

Several numerical modeling methods have been used to study wave propagation across discontinuous rock mass. For example, Banadaki and Mohanty [18] used ANSYS AUTODYN, one of the finite element methods, to conduct numerical simulation of rock fracture caused by explosive stress wave, while the Johnson–Holmquist model was confirmed to be applicable to the analysis of the rock fracture mode due to stress wave. By using the FLAC3D code which employs an explicit finite difference method, Sainoki and Mitri [19] conducted a simulation on rock mass vibrations induced by stress waves arising from production blasts. At the same time, they presented a methodology to calibrate blast vibration attenuation. Vorobiev and Antoun [20] proposed the equivalent continuum method to simulate nonlinear wave propagation in jointed media. Fan et al. [21] presented a numerical manifold method (NMM) simulation of stress wave propagation through fractured rock mass; later, they also verified the validity and accuracy of using the NMM in modeling fractured rock mass. Shao et al. [22] conducted a numerical study on isotropic elastic waves in the layered medium by the discontinuous Galerkin method. De Basabe et al. [23] applied the discontinuous Galerkin method for elastic wave propagation through fractures with linear-slip model, and their simulation results were calibrated with analytic solutions. Based on the detailed analysis of the damping, block size, and joint property in the discontinuous deformation analysis (DDA) method, Fu et al. [24] conducted a simulation on stress wave propagation across the rock mass. However, few researches are based on the abovementioned methods to study the influence of rock discontinuity on wave propagation and on the induced ground motion associated with.

In addition to those above numerical methods, the discrete element method (DEM) is more widely used and more efficient to analyze wave propagation in jointed rock masses, in comparison with the theoretical solution [25]. By using UDEC, a kind of code based on DEM, Brady et al. [26] conducted a numerical study on the slip of a single joint under an explosive loading. Fan et al. [27] performed UDEC modeling on stress wave propagation across a single rock joint. They also discussed the influence of velocity boundary and stress boundary on their simulation results. Zhao et al. [28] studied numerically the stress wave propagation though a set of joints with linear and nonlinear characters. Zhu et al. [29] conducted a simulation study on wave propagation characteristics between multiple intersecting joint sets. As for the cylindrical wave propagation across rock mass, Deng et al. [30] studied the damage mechanism of a circular tunnel under the action of cylindrical waves caused by blasting in jointed rock masses. Their studies combined two approaches: numerical simulation and field monitoring. Based on the above, Deng et al. [31] adopted AUTODYN to model the explosion process as well as to simulate shock wave propagation in jointed rock masses surrounding the explosion chamber by coupling UDEC. In comparison with the test data of an actual large-scale explosion case, the UDEC-AUTODYN hybrid modeling method was proved to be accurate and efficient. Yoo et al. [32] adopted a pseudostatic discrete element to analyze the circular tunnel dynamic response and to evaluate the effect of deformation under seismic wave. Among those above, relatively few researches focus on the influence of the structural plane on the ground vibration caused by the cylindrical wave. If a structural plane exists in the rock mass, it will affect the stress wave propagation as well as the induced ground motion caused by engineering blasting. On the one hand, the structural plane will attenuate the stress waves passing through; on the other hand, the superposition of multiple reflected waves generated between the structural plane and the ground surface will strengthen the ground vibration response. The influences of a structural plane on the cylindrical wave inducing ground vibration need further study.

The current study numerically simulates the ground motion induced by cylindrical wave propagation across a rock mass with a structural plane by UDEC. Firstly, theoretical analysis and numerical simulation are conducted, respectively, on the ground motion for a special rock mass without any structural plane. By comparing with the theoretical method for the special case, UDEC is proved to be efficient to model cylindrical wave propagation across rock masses and the caused ground motion. Then, the cylindrical wave propagation across the rock mass with a structural plane is simulated and the ground motion is calculated numerically. Finally, the parametric studies such as the structural plane stiffness, the frequency of incident wave, and the position of monitoring point of the ground are carried out.

2. Ground Vibration of a Rock Mass with No Structural Plane

2.1. Theoretical Formulations

Supposing the ground is a smooth and horizontal free surface, the ground motion caused by cylindrical stress waves contains two parts [3], namely, vibration caused by the body wave and vibration caused by the surface wave. In Figure 1, IP denotes the incident wave and RP and RS denote P- and S-wave reflected by the ground surface, respectively.

The particle velocity of the incident wave is denoted as and the particle velocities of the reflected P- and S-wave are denoted as and , respectively. Since the geometric attenuation occurs during cylindrical wave propagation across elastic rock mass, in addition, the amplitude of the cylindrical wave is proportional to ; thus, the incident wave reaching the ground can be expressed aswhere is the radius of the incident wave front, i.e., the distance between the point A and the wave source O. According to the reflection theory of elastic waves on a free surface [33], the velocities of reflected waves can be expressed aswhere and are the wave impedance of P- and S-wave, there are and , and and are the phase velocities of P- and S-wave in the intact rock; and are the emergence angles of the reflected P- and S-wave and satisfy Snell′s law, i.e., ; and is Poisson’s ratio in the intact rock.

According to the propagation characteristics and emergence angles of the incident and reflection waves, the horizontal and vertical components of the body wave caused vibration velocity of the ground surface particle, i.e., and , can be, respectively, derived by combining equations (1) and (2) and are expressed in the matrix form as

Since the intensity, direction, and time delay of the waves arriving at the adjacent particles are not equal, the differential velocity between arbitrary two adjacent points on the ground and will cause surface wave, which is expressed as [3]where and are the horizontal and vertical components of the surface wave at ground point . Surface wave amplitude attenuation conforms to , where d is the propagation distance. The surface wave velocity can be obtained by the approximate equation given by Achenbach [34]. Therefore, the vibration of point A, as shown in Figure 1, caused by the propagation of surface wave generated at any point denoted as xi can be expressed aswhere is the distance between the ground point and point A (the horizontal coordinate is ). To sum up, the horizontal and vertical components of the vibration velocity of the ground surface point A can be expressed as [3]where is the interval distance between two adjacent points on the ground.

2.2. Numerical Simulation by UDEC

The UDEC model of ground vibration caused by cylindrical wave propagation in a rock mass without any structural plane is shown in Figure 2(a). The model is semisymmetric, and its dimension is 150 m × 40 m. The bottom boundary is set to a symmetrical boundary, while the top boundary is taken as a free boundary to simulate the ground surface. In order to eliminate the effect of reflection wave superposition, the left and right boundaries are viscous boundaries [35]. The mesh size of the block is 1 m, which is much smaller than the wavelength.

Because of the complexity of the cylindrical charge explosive process, the initial wave is simplified into a cylindrical elastic semisinusoidal stress wave acting evenly on the concave surface of the semicircular cylindrical cavity, as shown in Figure 1. The initial wave is applied by uniform distributed velocity loading. It is a half-cycle sinusoidal stress wave along the radial direction, with a frequency of 50 Hz, amplitude of 1 m/s, and duration of 0.01 s. The rock media in the model are assumed to be homogeneous and elastically intact. The mechanics parameters of the rock media are shown in Table 1, which are consistent with Zhao’s study [36]. In order to focus on the analysis of the influence of the structural plane on wave propagation, the material damping of rock is ignored and set as η = 0.

As shown in Figure 2(a), a set of monitoring points Pn are arranged on the ground surface. The horizontal distance between the monitoring point and O′ is h. If the horizontal coordinate of O′ is 0, the coordinates of other monitoring points can be calculated. The horizontal and vertical vibration velocities of each monitoring point are recorded as and , respectively.

Figure 3 shows the waveforms of and at four monitoring points on the ground surface, where the horizontal coordinates of each point are X0 = 0, 20, 40, and 60 m, respectively. For any monitoring point, the vibration waveform of either vertical or horizontal components is similar to that of the incident wave, which is roughly half sinusoidal. The maximum value of the waveform curve is defined as the peak particle vibration velocity, i.e., PPV. As shown in Figure 3, each vibration wave curve has a certain time delay, which is the required time of that the wave travels from the initial wave front to the ground surface. Therefore, with the increasing X0, the wave propagation distance increases, which leads to a longer time delay. Thus, for a certain point except that X0 = 0 m, the vertical and horizontal vibration wave curves have a same time delay. Since the stress wave is vertically incident upon the ground surface at the point X0 = 0 m, only vertical vibration is generated, and the horizontal vibration velocity should be 0. However, when the cylindrical wave reaches the ground surface, due to the velocity difference between two adjacent particles, the surface wave will be generated, as described in Section 2.1. Thus, the horizontal velocity waveform at X0 = 0 m shown in Figure 3, in all probability is caused by the surface wave. The time delay of horizontal velocity curve is longer than that of vertical velocity because surface waves occur later than body waves do. For other vibration wave curves shown in the figure, there is a jumping point in the descent section after the wave reaches the peak, forming a second small crest, which is mainly the result of the vibration superposition caused by body wave and surface wave. Moreover, as X0 increases, the time interval between surface wave and body wave acting on the point increases.

2.3. Verification of the Cylindrical Wave Propagation

From the numerical simulation, the horizontal and vertical components of the peak vibration velocity of several chosen points on the ground surface are obtained and plotted in Figure 4. The results for each case calculated based on the theoretical method as described in Section 2.1 are also drawn in Figure 4. Obviously, the PPVs obtained from UDEC are very close to the ones calculated by the theoretical method. The possible reason to cause the difference between the theoretical and numerical methods is that the models of these two methods are not completely consistent. The theoretical model is designed in semi-infinite space, which is not fully implemented by UDEC. There are some small errors caused by grid division, boundary setting, and other factors in the UDEC model. However, as shown in Figure 4, the error is small enough to be ignored. It indicates that the UDEC can be used to simulate the ground motion of rock mass caused by cylindrical wave propagation.

Moreover, the figure shows that the horizontal and vertical components of PPVs vary with different locations of the ground points. The horizontal component of PPV, Vτ, increases at first and then decreases with the increasing of the horizontal coordinate of the monitoring point. Differently, with the increase of the distance on the ground point, the vertical component of PPV, Vn, almost monotonously decreases, and the decrease trend is not linear. Additionally, obvious inflection points can be observed in both curves, with x coordinate approximately 35 m.

3. Ground Vibration of a Rock Mass with a Structural Plane

The numerical model to study the ground vibration caused by the cylindrical wave propagation in rock mass with a structural plane is shown in Figure 2(b). In the model, the structural plane is supposed to be a fictitious plane, parallel to the ground surface. The normal distance between the structural plane with ground surface is denoted as S (0 < S < 35 m). Other conditions of the model are the same as those in Figure 2(a). The normal stiffness of the structural plane is denoted as kn, while the shear stiffness is denoted as ks. Similarly, histories of the vibration velocity components at the monitoring point, arranged on the upper boundary of the model, are monitored and denoted as and .

The ground vibration response varies with changed parameters, such as the stiffness of the structural plane, distance between the structure and the ground, incident wave frequency, and the location of monitoring points. Therefore, based on the calculation results of numerical simulation, the results of the parametric analysis on the effects of these factors contributed to the ground motion are discussed in this section.

3.1. Effect of the Distance between the Structural Plane with Ground Surface

Figure 5 shows the relationship between Vτ and Vn at 3 points on the ground with S. The coordinates of the 3 points, X0, are 0, 20, and 40 m, respectively. In the calculation, kn = 3.5 GPa/m, ks = 1.0 GPa/m, respectively. It can be seen from the figure that the PPVs vary with different distances between ground surface and the structural plane. The variation tendencies of the peak particle vibration velocities with S of the three points are very similar. With the increase of the distance between ground surface and the structural plane, the horizontal component of PPV Vτ slightly reduces, while the vertical component of PPV Vn increases firstly and then decreases. The situation of S = 0 m can be regarded as the special case that there is no structural plane in the rock mass.

When S = 0 m, i.e., there is no structural plane in the rock mass, and the body waves reach directly at the ground and cause ground motion. The surface waves induced from the differential particle velocities can also superimpose the ground vibration. Since a structural plane exists in the rock mass, e.g., S = 5 m, multiple reflected waves are induced between the structural plane and the ground, which make contribution to the ground vibration as well. When the stiffness of the structural plane is large, the strength attenuation of the transmitted waves through the structural plane is small. The increases of ground particle velocities caused by the multiple reflected waves will be greater than the decreases of ground vibration caused by the decreasing strength of the transmitted body waves. Thus, taking the multiple reflections of stress wave between the structural plane and ground into account, the value of Vn for S = 5 m is naturally greater than that for S = 0 m.

As shown in Figure 5, the vertical component of the PPV value in the case with the presence of a structural plane is somewhat larger than that in the case with the absence of the structural plane in a certain range. Take the case X0 = 0 m as an example, is approximately 0.59 m/s for S = 0 m, is approximately 0.82 m/s for S being about 6 m, while is approximately 0.58 m/s for S = 20 m. The maximum Vn is recorded near S = 6 m, and it is 1.4 times the value of Vn in the case with the presence of a structural plane. For the other given cases, i.e., X0 = 20 m and X0 = 40 m, the curves of Vn look very similar to that for S = 0 m. As the two curves show, Vn increases firstly and then decreases with the increase of S, and finally it reaches maximum when S is also about 6 m. Similarly, within the range of 0 < S < 20 m, the value of Vn is greater than that in the case with the presence of a structural plane (i.e., S = 0 m). The maximum magnification factors (MMFs) of PPV (X0 = 20 m; X0 = 40 m) are approximately 1.5. The values of MMF over 1 (; ) indicate that the structural plane reflects more stress waves when S is relatively small (0 < S < 20), and those reflected waves arrive at the ground again and again, finally result in the superposition of ground motion. As shown in Figure 5, the higher the S value is, the lower the value of vertical component of PPV is. Especially when S > 20 m, the vertical component of PPV with structural plane, Vn, is smaller than that without the structural plane. The reason is that the number of multiple reflected waves generated between the ground and the structural planes within a same period will be reduced with the increasing S. The decreases of multiple reflected waves make the ground vibration reduce. When S = 20 m, the shortest time required for the transmission wave from the structure surface to reach the ground surface once again after multiple reflections is about 0.103 s (t = 3S/cp = 3 × 20/5830 s). It is larger than the duration of the incident wave, i.e., 0.01 s. Thus, the multiple reflected waves contribute rarely to the ground vibration (MMF < 1). To sum up, when S > 20 m (MMF > 1), the superposition effect of multiple reflected waves on surface vibration can be ignored.

3.2. Effect of the Position of Monitoring Point

In order to describe the location of the monitoring point, the position angle is defined as θ = tan−1(h/l), shown in Figure 2(b). Figure 6 shows the relation between PPVs and θ at the case that S = 20 m. For the convenience of comparison, the curves of the PPVs in the case without the presence of the structural plane are shown as the dashed line in Figure 6.

It can be seen from the figure that the horizontal component of PPV, Vτ, at θ = 0° for both the given cases (with/without the presence of structural plane) is minimum and is approximately equal to zero. That means the transmitted S wave will not generate under this condition (θ = 0°) that the incident P-wave impinges vertically on the structural plane. Besides, with the increasing of the absolute value of the position angle (i.e., incident angle), Vτ firstly increases until |θ| reaches about 45°, and then it decreases progressively. As mentioned before, the PPV components along the horizontal and vertical direction are the combinations (along the horizontal and vertical directions) of the components of P- and S-wave occurred on the ground. Since the position angle, θ, increases, the incident angle of the wave impinging on the ground increases. Meanwhile, the emergent angles of reflected waves increase accordingly, which leads to the decreasing velocity in the vertical direction and the increase velocity in the horizontal direction. Furthermore, the attenuation increases with the increasing position angle. That will result in the reduction of each wave’s intensity. To sum up, the tendency of the curves is driven by the combination of the above two aspects.

Moreover, for the same monitoring point, the horizontal component of PPV for the situation without structural plane is higher than that with a structural plane. This difference occurs because that the reflection of the structural plane on the stress wave weakens the strength of the body waves reaching the ground surface. As shown in Figure 6, Vn reaches the maximum value when the position angle is zero and decreases with the increase of the absolute value of the position angle. In the case without structural plane, Vn almost monotonously decreases with the increase of |θ|. While in the case with a single structural plane, small fluctuations are found along the Vn curve. The main reason of those fluctuations is that the multiple reflected waves between structural plane and ground surface will reach the ground and strengthen the particle vibration of the ground.

3.3. Effect of Structural Plane Stiffness

Figure 7 illustrates the peak particle velocities varying with the normalized stiffness of structural plane at the points X0 = 0, 20, and 40 m on the ground, where S = 20 m. The symbol Kn denotes the normalized stiffness of the structural plane, and the normal stiffness and shear stiffness of the structural plane are and , where . The normalized stiffness is adopted from 0–4.

From Figure 7, with the increasing of Kn, both Vτ and Vn show an increasing tendency. The rising rate is relatively high for Kn less than 0.5. At a certain monitoring point, Vn is greater than Vτ when the normal stiffness is given. Meanwhile, Figure 7 shows that the PPVs tend to be stable with the increase of structural plane stiffness when Kn is greater than 1; in other words, the continuous increase of stiffness has little effect on the peak particle vibration velocities.

When the structural plane is very soft, i.e., Kn near to zero, the horizontal and vertical components of the PPV are zero which means all the waves are reflected from the structural plane. As theoretically identified and described by Li [9], there are only reflected waves produced when stress waves impinge on the kind of weak discontinuities acting as vacuum-free surfaces.

3.4. Effect of the Incident Wave Frequency

The relationship between the PPVs and the frequency of the incident cylindrical wave at the points X0 = 0 and 20 m is illustrated in Figure 8, given that S = 20 m. In calculation, the normal stiffness and shear stiffness of the structural plane are adopted as kn = 3.5 GPa/m and ks = 1.0 GPa/m, respectively. It can be observed from Figure 8 that Vτ and Vn vary with increasing frequency. Specifically, Vn monotonously decreases with f. The curve looks like an inverted “S.” Vn decreases smoothly when the frequency is relatively small (f < 20 Hz) or relatively large (f > 100 Hz) but drops sharply in the range of f = 20–100 Hz. With the increasing wave frequency, Vτ decreases rapidly when the frequency is less than 30 Hz, then decreases slowly, and finally tends to be constant except that small fluctuations exist in the two curves at approximately 30 Hz.

By direct comparison, we can see that the value of Vn at X0 = 0 m is greater than that at X0 = 20 m for a given frequency of incident wave, while the value of Vτ at X0 = 0 m is smaller than that at X0 = 20 m. What is more, when f is relatively high, Vτ and Vn are lower. It reflects that the propagation of high-frequency waves on the structural plane or on the ground is more difficult. In other words, the PPVs have a low value when the frequency is relative high which reflects the structural plane filtering effect to high frequency.

3.5. Assumptions and Limitations of the Numerical Method

It should be noted that the proposed numerical approach has some assumptions. For example, the natural rock mass contains a large number of microcracks and joints which make the rock mass to present nonlinear characteristics. In this simulation, the rock mass is assumed to be linearly elastic. In addition, the structural plane in the rock mass is not generally parallel to the ground surface, which is assumed to be parallel to the surface in this paper. Moreover, the mechanical properties of natural structural plane are complex; however, it is assumed to be a linear elastic model in this paper. In view of these assumptions, the simulation has limitations and can only be applied to the special case of the abovementioned rock mass. In spite of those assumptions, the results in this paper still reflect a certain rule of ground vibrations caused by cylindrical wave propagation across a rock mass with a structural plane. Furthermore, this method can also provide some references for the further study of wave propagation and the caused vibration in complex rock mass with multiple structural planes.

4. Conclusions

The ground vibrations caused by cylindrical wave propagation across a rock mass with a structural plane are investigated numerically by using the UDEC method. Certain conclusions can be obtained. Based on the simulation results, the ground vibrations are mainly affected by three factors: the position of monitoring point, stiffness of structural plane, and frequency of the incident wave. For the cylindrical wave propagation across a rock mass with a structural plane, the PPVs vary at different ground points. The vertical component of PPV Vn is maximum when the absolute value of the position angle θ = 0°, while the maximum point for the Vτ curve is at |θ| = 45°. In addition, the PPVs are sensitive with the normalized stiffness of the structural plane, Kn. The horizontal and vertical components of PPV increase rapidly at low values of the normalized stiffness, and then both of them tend to be stable at a relatively higher value of structural plane stiffness. Moreover, the higher frequency of incident wave causes the greater the intensity of the reflected waves from the structural plane, which reflects the filtering effect of the structural plane on high frequency.

Data Availability

The data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest

The authors declare that they have no conflicts of interest regarding the publication of this paper.

Acknowledgments

The study was supported by the State Key Laboratory for GeoMechanics and Deep Underground Engineering, China University of Mining and Technology (grant no. SKLGDUEK1811), the National Natural Science Foundation of China (grant nos. 41902277 and 51579239), and the Natural Science Basic Research Plan in Shaanxi Province (grant no. 2019-JQ689).