Abstract

Critical for the seismic safety of immersed tunnels is the magnitude of deformation and force developing in the segment joints. To investigate the seismic performance of segment joint in immersed tunnel, this paper presents a series of shaking table arrays tests that were performed on a microconcrete tunnel model embedded in the soil. The tests take account of the uniform and wave passage effect in different apparent wave velocity of longitudinal seismic excitation. The result showed that the wave passage effect had a great impact on axial force, bending moment and deformation of joints. The comparison showed that structural response under nonuniform earthquake excitation is larger than that under uniform excitation. The simplified model was established in ABAQUS for numerical analysis. The soil around the tunnel was simplified as spring-damper, the tunnel was simplified as beam element, and the joint was simulated by nonlinear spring. The numerical simulation results were in good agreement with the experimental data. In addition, the model was analyzed by changing input apparent wave velocity, joint stiffness and joint number. The results showed that the deformation of joints was smaller and the deformation of flexible joints was greater under high apparent wave velocity.

1. Introduction

Many immersed tunnels have been constructed in soft ground at port areas, where the response of both the soft ground and the immersed tunnels are amplified during earthquakes. Since China is located in one of the most active seismic zones in the world, the earthquake resistance of these immersed tunnels must be checked from various points of view, especially the magnitude of deformations developing in the segment joints. The seismic excitation may affect the overall stability of the immersed tunnel, leading to decompression of the joint gaskets, jeopardizing the water tightness of the tunnel. Presently, no uniform criteria have been available for the seismic design of segmental joints. The safety of segmental joints is especially important and worth of extraordinary attention [1].

Recently, more and more researchers have got to investigate the immersed tunnel. Anastasopoulos [2] applied finite element method (FEM) to investigate the seismic response of a very-deep immersed tunnel in Greece, under the simultaneous action of longitudinal, transversal, and vertical seismic excitation. The joints between the tunnel segments are modeled realistically with special nonlinear hyperelastic elements. Oorsouw [3] investigated the behavior of segmental immersed tunnel subjected to seismic loading and especially on the sensitive segment joint. Jun-Hong Ding[4] used the finite-element code LS-DYNA to analysis the large-scale seismic response of immersed tunnel. The behavior of nonlinear material such as soil, nonreflecting boundary definition, and soil-tunnel interaction were taken account. Yang Ding [4] and Tayloy [5] established numerical models of immersed George Massey tunnel for nonlinear dynamic geotechnical analyses and the centrifuge models was used to verify and calibrate the numerical models. However, most of the literature carries out theoretical and numerical simulation analysis on immersed tunnel, and few of the actual experimental studies are carried out. Asynchronous ground motion does not significantly affect buildings, but they may have significant impact on the seismic response of extended structures such as bridges and tunnels. Consideration of spatially varied seismic input motion for such structures is challenging owing to its complexity [6]. Haitao Yu [7, 8] carried out multipoint shaking table tests on a long immersed tunnel. Jun Chen [9, 10] investigated the utility tunnel under nonuniform earthquake wave excitation. However, the force of the joint was not measured in the test. Immersed tunnel is a long linear structure, it is very important to consider the nonuniform earthquake excitation for seismic response analysis. Seismic performance of segmental joints is particularly important for tunnel safety. Therefore, it is necessary to carry out shaking table tests to research the seismic performance of segment joints under wave passage effect excitation.

The Zhoutouzui tunnel is a newly built immersed tunnel in Guangzhou. The whole length is 2200 m, including immersed tunnel (about 340 m) under the Pearl River [11], as is shown in Figure 1. The first and last segments are variable cross section. The geotechnical engineering investigation report on the Zhoutouzui variable cross section immersed tunnel indicates that, the tunnel is laid on rock-soil layers with variable physical mechanics properties, and lithologies in this area puts up quite complex, in some parts intensely weathered zones and weak weathered zones are in alternating layers, which makes the mechanics analysis of this tunnel under earthquake quite difficult.

To investigate the mechanical properties of immersed tunnel joint under nonuniform earthquake, this paper presents a series of shaking table arrays tests that were performed on an immersed tunnel model embedded in the soil. The wave passage effect was considered. The main results obtained from the tests were summarized and discussed. Further, a FEM of the test immersed tunnel was established. The results obtained from the suggested numerical model are compared with experimental measurements in terms of displacement. The comparisons show that the numerical results match the experimental measurements quite well. For further analysis, the different stiffness of joints, different apparent wave velocity and multi joints were taken into account.

2. Experiment Setup

2.1. Shaking Table Array

In this experimental test, the shaking table array system at Beijing University of Technology (BJUT, Beijing, China), which was able to support multisupport inputs, was used for the dynamic experimental test of the immersed tunnel under asynchronous excitations. Each table consisted of a 1 × 1 countertop, actuators, link rods, and a base, as shown in Figure 2. The performance parameters of the table array are shown in Table 1 [12]. In this experiment, four independent shaking tables were used, and the space between each table was 1 m.

2.2. Test Model

For the test model under dynamic loads, the physical parameters of the structural dynamic characteristics were as follows: (1) structural geometry (l); (2) material properties, including elastic modulus (E), mass density (ρ), Poisson’s ratio (µ), stress (σ), and strain (ε); (3) loads, including force (F) and moment (M); and (4) dynamic indicators, including time (t), frequency (ω) damping ratio (ξ), acceleration (a), stiffness (k), and mass (m). Thus, the solution equation for structural dynamics problems can be expressed as [13]

First, because the aim of the experiment was to investigate the response of an immersed tunnel under seismic excitation, the influence of gravity on the model was ignored. Second, a strength model with the same materials as the prototype was adopted, e.g., the materials and the strain of the tunnel model were the same as that of the prototype. Thus, it was assumed that (1) vertical force did not affect the transverse stiffness of the structure, i.e., the stiffness similitude parameter Sk = SESL; (2) the strain of the tunnel model was the same as that of the prototype, i.e., strain similitude parameter Sε = 1; and (3) the damping ratio of the tunnel model was the same as that of the prototype, i.e., the damping ratio similitude parameter Sξ = 1. According to equation (1), the similitude relations of the tunnel model could be obtained by dimensional analysis.

Considering the limitations of the laboratory space and the actuation capacity of the shake tables, the test immersed tunnel was scaled to 1/60 of the prototype tunnel in geometric dimensions. The tunnel was made out of microconcrete, and galvanized steel wire gauze was used for the reinforcement. The microconcrete mixture ratio of the experiment was: cement (42#): sand: lime: water = 1 : 6.0 : 0.6 : 0.5. The modulus of elasticity was 7410 N/mm2. The compressive strength of the concrete cube was 5.679 N/mm2. The density of the microconcrete was similar to that of the prototype concrete, so the similitude parameter of density was set as 1, and the elastic modulus similitude parameter was set as 1/4. According to the scale factor, the model of two middle segments corresponded to a 503 mm × 160 mm rectangle tunnel having an equivalent concrete lining thickness of about 20 mm. The first and the end models corresponded to a variable cross section, the dimensions of which are shown in Figure 3.

In this experiment, the most important requirement for the similitude relations of the dynamic test model was to determine the similitude parameters of acceleration and time. The influence of the gravity soil model was also ignored. Because it was difficult to meet the requirements of similarity relation between the model and prototype, the predominant period of similarity relation was used for the design of the model soil, ensuring the predominant period of prototype soil and model soil similarity. The shear wave velocity was changed with depth and the property of the soil. The predominant period of the soil can be calculated by the equivalent shear wave velocity:

Thus,due to

The shear wave velocity of soil iswhere m is the model and P is the prototype.

In this model soil, according to the similarity ratio.

According to the geotechnical engineering investigation report, the equivalent shear wave velocity was 250 m/s. According to (6), the equivalent shear wave velocity of the model soil was set as 88.75 m/s. The results of the similitude relation are presented in Table 2.

To obtain an appropriate equivalent shear wave velocity, clay was chosen, and different proportions of sawdust were mixed in as test samples. Three groups of samples with different proportions were selected for the resonant column test. The proportions were: 1 : 2.5, 1 : 3 and 1 : 3.5. The max shear modulus and shear wave velocity are shown in Table 3. Finally, sawdust: prototype: water = 1 : 3:2.7 was chosen for the experiment tests [14]. The height of soil on the prototype tunnel was 2.29 m, and the density was 2 g/cm3. The depth of water was 6 m, and the density was 1 g/cm3. The water that was converted to the depth of soil was equal to 3 m. According to the similitude relations, the embedded soil depth was 14 cm.

2.3. Joint

Limited by the test conditions, the similarity ratio of the model is set as 1/60. As a result, many parts of the joint cannot be made according to the prototype structure. Considering the joint failure caused by excessive tensile displacement at the joint, the design focus of the joint is to design the tensile stiffness of the model. The shear key at the joint is designed to make the joint work normally, without considering the similarity of the shear key.

The longitudinal deformation of the tunnel depends significantly on the gasket and prestressed tendon. In this work, 38 sets of tendons were used in the prototype tunnel. The whole tension stiffness was 83,406.2 KN/m. The compressive stiffness of the gasket was simplified to the bilinear model by the result of the finite element analysis, as illustrated in Figure 4. The Gina gasket was simplified to a rubber ring with a section equal to the cross section of the tunnel model, and the thickness was 2 cm. The angle steel was embedded at the end of the tunnel model. The angle steel embedded on the side of the tunnel was cut, and the steel bar welded with angle steel was bound with steel wire gauze. The end steel shell was welded to the angle steel, and the rubber ring was glued to the side of the end steel shell. The horizontal and vertical shear key was into the rectangular ring, and its thickness was 2 mm. One end of the shear key was welded to the end steel shell, the other end was free, so that the free end could insert into a moveable middle steel shell with a section larger than the section of the shear key. The gap in the shear keys was to prevent collisions during the test, as is shown in Figure 5(a). The prestressed cables were simplified to six bolts with a diameter of 4 mm. The effective tensile length of the cables was 55 mm, as is shown in Figure 5(b). The cables were prestressed by tightening the nuts. The bolt with a large diameter in the middle of the joint was used to protect the joint from damage during installation. When the position of the tunnel was ensured, the bolt was removed.

Four boxes were fabricated and installed on each shaking table to load the soil in the test. The overall dimension of the box was 7.3 m × 3.2 m × 1.2 m (length × width × height), as is shown in Figure 6(a). The first and the last boxes were 1.5 m × 3.2 m (length × width), and the two middle boxes were 2 m × 3.2 m (length × width). The box frames were welded by angle steel, with a dimension of 70 mm × 70 mm × 5 mm. The bottom of the box was a steel plate, which had a thickness of 10 mm. For the purpose of achieving nonuniform seismic excitation, the gap between each box was 100 mm. The bottom of the two adjacent boxes was lapped with steel plates, and bolts were fixed between the steel plates. Butter was smeared between the contact areas to decrease friction. Before conducting nonuniform seismic excitation tests, the bolts were removed. Both sides of the two adjacent boxes were inserted into square steel tubes, with a dimension of 100 mm × 100 mm. The square steel tubes and the boxes were connected by bolts, as shown in Figure 6(b). When the nonuniform excitation was applied, the square steel tubes were removed, as shown in Figure 6(c). To avoid large deformation at the bottom of the box, joist steel and stiffener was welded under the plate. Rubber sheet and polystyrene foam board was used to decrease the boundary effect. The thickness of the rubber sheet and foam board was 15 mm and 200 mm, respectively. The rubber sheets and box walls were fixed by multiple small bolts. The foam boards and rubber sheets were filled with spray foam.

The whole model of the test was a system composed of model boxes, soil and a tunnel embedded in soil. To ensure that the vibration of the model box did not affect the dynamic response of the model soil, the natural frequency of the model box was far away from the natural frequency of the model soil; that is, the natural frequency of the model box was much higher or lower than that of the model soil. Moreover, the natural frequency of the model should be less than the maximum working frequency of the shaking table.

To verify the applicability of boxes, numerical analysis of the box was carried out by ABAQUS. Beam and shell elements were used to simulate the frame and floor of the box, respectively. The simulation results showed that the fundamental frequency of boxes at both ends was 33.9 Hz. The fundamental frequency of the middle boxes was 22.6 Hz. The four boxes were tied together, and the overall fundamental frequency was 23.8 Hz. When the soil was taken into consideration, the frequency of the box-soil was 11.6 Hz. The frequency difference between the whole box and box-soil was nearly two times, so the box could meet the test requirements.

3. Shaking Table Test Program

3.1. Model Installation

To avoid relative sliding between the soil and box bottom, gravel concrete was prepoured on the box bottom and scratched with a broom. The model soil was compacted manually every 20 cm. When the free-field test was completed, the soil around the model tunnel was excavated, and the tunnel was buried in the soil. The top of the soil was then covered with a plastic sheet to prevent evaporation, and the soil was compacted with adhesive weight, as shown in Figure 7.

3.2. Instrumentation

The layout of sensors is depicted in Figure 8. A stands for the accelerometer, D stands for displacement and F stands for force sensor. The accelerometers were applied to the soil and tunnel. A1∼A4 were arranged on the tunnel to measure the acceleration response of the tunnel, A5–A7 were arranged on the joints and A8–A11 were buried in different heights of the soil. Laser displacement sensors (D1∼D3) were arranged on the joints. Twelve force sensors were applied on the joints, and each joint was arranged with four force sensors. The positions of the sensors are shown in Figure 8(c). To prevent damage to the accelerometers, they were wrapped in balloons before being buried into the soil, as shown in Figure 9(a). The force sensors were connected with the steel shell through the screw, and the two ends of the screw were fixed by bolts on both sides of the steel shell, as shown in Figure 9(b). The laser displacement sensor was fixed on a horizontal shear key on one side of the joint, and a baffle was placed on the other side of the joint to receive the laser signal, as is shown in Figure 9(c) and 9(d).

3.3. Test Case

To simulate time lag in the arrival of the waveform at each box bottom, considering the wave passage effect, it was assumed that the wave motion propagated from shaking table 1 to 4, meaning that the time difference in the wave motion between adjacent shaking tables is [15]where D is the distance between the adjacent shaking tables along the direction from shaking tables 1 to 4, and Vp is the apparent wave velocity. Chi-Chi ground motion acceleration was applied to wave velocities of 100 m/s, 200 m/s, 300 m/s, 600 m/s, and infinity (uniform). The magnitude was gradually increased from 0.06 g to 0.4 g. According to the scaling law, the frequency and dynamic time in the shaking table tests were n time and 1/n those of the prototype, respectively. The sampling frequency of the data acquisition was set at 1000 Hz.

4. Test Results and Analysis

4.1. Acceleration Response

Figure 10 shows the variation of peak value of A10 with the change of input peak ground acceleration (PGA). It was seen that for PGA< 0.19 g the soil had good linearity. When PGA >0.19 g, the peak value exhibited distinctly nonlinear behavior. Figure 11 shows the amplification factor (AF) of soil and tunnel under different apparent wave velocity in PGA = 0.12 g and 0.26 g. As can be observed that when PGA = 0.12 g, AF of the soil and tunnel were 1.63 and 1.65 under uniform excitation. However, for PGA = 0.26 g, AF of those were 1.42 and 1.56, the difference was larger than those in PGA = 0.12 g. It is also found that the AF of tunnel decrease with the high input PGA. The reason is that the higher of the input PGA, the more obvious of the soil nonlinearity is shown. When considering the wave passage effect, the AF of tunnel was larger than that of soil under input PGA = 0.12 g. However, for input PGA = 0.26 g, the AF of the tunnel was smaller than that of soil (except for the case of Vp = 300 m/s). This is caused by the nonlinear characteristics of soil, and the wave passage effect, which will lead to the sliding between soil and tunnel.

Due to space limitations, only input PGA = 0.12 g is presented in this paper. It is interesting to compare the actual input excitation of the shaking table with that developed in the soil and tunnel. For this purpose, the time history and corresponding autospectrum of sensors A4 and A10 are depicted in Figure 12. It is seen from these figures that the acceleration time-history curves of A4 and A10 had a large amount of synchronism under uniform seismic excitation. When considering the wave passage effect, the movement of the tunnel (A4) was slower than that of soil (A10). The autospectrum figures show that the curve of A4 was similar to that of A10 under uniform seismic excitation. However, autospectrum figures present multiple peak values under wave excitation, making the spectral composition much richer than those in uniform excitation.

4.2. Axial Force of Joints

Because the initial pretension could not be measured in dynamic experiment, so the initial value of each force sensor was set to zero. Under uniform seismic excitation, the joint 1 (J1) mainly undertake the pressure, and the J3 always undertake the tension. Compared to the J1 and J3, J2 experience less force, as is shown in Figure 13 The range of force varies within 120 N. It is explained that the tunnel essentially “follows” the uniform excitation. However, when considering the wave passage effect, the axial force was amplified about 10 times. Central joint experienced less force than the terminals and the maximum value of tension was occurred in J1. J3 experienced the biggest pressure in all the case. The maximum of tension was occurred in the case of Vp = 100 m/s, the value of which was 1436 N. According to the similitude relations, the maximum axial force of the prototype was 20822 KN. The mean tensile stress of prestress tendon was 20.9 MPa, adding the initial stress 2 MPa. The stress was within acceptable limits in all cases.

4.3. Bending Moment of Joints

Figure 14 illustrates the bending moments Mz and MY in different apparent wave velocity cases. It is shown that J1 and J3 had identical bend direction in uniform excitation. When considering the wave passage effect, the bending moment decreased along the wave propagation direction. Mz was larger than MY. Changing the Vp led to a slight difference of the bending moment.

4.4. Deformation of Joints

Deformation time histories of the immersion joins are portrayed in Figure 15 and the maximum value of deformation is shown in Figure 16. The gaskets experience a slightly initial compression by the bolt. We set the initial compression is zero. The seismic oscillation caused successive decompression and recompression of the gaskets. Considering the wave passage effect, joints experienced more relative displacement than the ones in uniform excitation. While for Vp = 200 m/s, the decompression Δmax = 0.047 mm. Converting into the prototype tunnel, the deformation was 28.2 mm, less than precompressed value (50 mm), so the deformation within acceptable limit.

5. Finite Element Analysis

5.1. Finite Element Model of the Immersed Tunnel

The finite element model (FEM) was used to perform nonlinear dynamic transient analysis of the tunnel by ABAQUS. The model layout is depicted in Figure 17. Tunnel segments were simulated using beam elements. Each immersion joint was modeled with six nodes, which were rigidly connected to the single end beam with special transitional rigid elements [2, 16]. Adjacent nodes were connected to each other with single-degree-of-freedom nonlinear springs, representing the stiffness of the joint. All beams were connected to the soil through interaction springs and dashpots. The first and last segments were also connected with springs and dashpots.

The analysis was conducted in two stages. First, static pressure was applied to the end of each segment to simulate the initial hydrostatic longitudinal compression. At the second stage, the model was subjected longitudinal dynamically to earthquake shaking. The acceleration time histories, which were measured by A10, were applied to the supports of springs and dashpots with the experimental time lag.

5.2. Soil-Tunnel Interaction Parameters

To obtain proper values of the longitudinal (x) and vertical (z) supporting spring and dashpot constants, a rigid long rectangular foundation on half-space was utilized [15]. For the surface of a homogeneous half-space, the vertical static stiffness is determined according to

Longitudinal stiffness iswhere ; G is the shear modulus that get from Table 3; 2L = 5.75 m is the length of the model tunnel; 2B = 0.6 m is the width of the tunnel; ν is the Poisson’s ratio of soil and the Ab is the area of tunnel.

The dynamic stiffness iswhere K is the static stiffness of “spring”; is times the dynamic stiffness coefficient. is plotted in Figure 17.

The vertical radiation dashpot coefficient is determined bywhere , Vs is the shear wave velocity of soil, and ρ is the density of soil. The value of Vs and ρ were selected from table 3. is plotted in Figure 18.

The longitudinal radiation dashpot coefficient is

The total dashpot is determined using

Dynamic stiffness and radiation dashpot coefficient are shown in Figures 18 and 19.

For the fully embedded tunnel in homogeneous half-space, the static stiffness iswhere D = 0.16 m is height of the model tunnel and is actual sidewall-soil contact area.

The vertical dynamic stiffness is obtained by

Longitudinal static stiffness is

Dynamic stiffness  = 1 is estimated in terms of L/B, D/B, and d/B.

The vertical radiation dashpot coefficient is

The longitudinal radiation dashpot coefficient is

5.3. Joint Spring Parameters

The joints between the tunnel segments were modeled with nonlinear springs. When they compress, the springs refer to the Gina gasket. When the springs are in tension, the stiffness of springs is equivalent to the rigidity of prestressed tendons. The force-deformation is shown in Figure 20. Assuming that the stiffness of the Gina gasket is denoted by k0, different Gina gasket stiffness (1/4 k0, 1/2 k0, 3/4 k0) is considered in this paper.

5.4. Dynamic Analysis Results
5.4.1. Comparison of Experimental and Numerical Results

To verify the rationality and reliability of the model, the data of A10 in the test was regarded as the input ground motion in FEM. The baseline of the data was corrected, the data was filtered, and 0–50 Hz was reserved. The wave was input along the longitudinal direction of the model, and the time lag was the same as the experiment. Figure 21 shows the deformation of the test and FEM in the case of PGA = 0.12 g. The trend of the two curves was consistent, and the maximum tensile value of the two curves was close. However, there was permanent compression deformation in the test curve, and it was shifted downward after 2 s. If the effect of permanent deformation was ignored, the curve of the experiment and FEM were consistent with each other.

5.4.2. Results of Different Apparent Velocities

Different apparent velocities from 300 m/s∼600 m/s were input into the model. The maximum deformation of J1 is shown in Figure 22. When the apparent wave velocity was less than 1000 m/s, the value of tension was obviously greater than that of compression. When the velocity was greater than 1000 m/s, the difference between tension and compression was small. Besides, with the increase of wave velocity, the maximum value decreased.

5.4.3. Joint Stiffness

To consider the deformation of joints under different stiffness, four kinds of stiffness were considered in the FEM. The apparent wave velocity was 600 m/s. The results are shown in Figure 23. When the joint stiffness was k0, the maximum compression deformation was at the J1, and the minimum compression deformation was at the J3. When the stiffness of the joint changed to 1/4 k0, the maximum compression deformation appeared at the J2. As the stiffness decreased, the compression and tensile deformation of joints increased gradually, indicating that the flexible joints could withstand greater deformation.

5.4.4. Multijoints

To investigate the effect of the number of joints on the seismic performance of the immersed tunnel, three joint and four-joint cases were considered in the numerical analysis. The same ground motion was input to the model, and the apparent wave velocity was 600 m/s. Figure 24 illustrates the deformation time history curves of each joint. The time history curves indicate that the shape of curves with three joints and four joints are quite different. When there are three joints in the tunnel, the maximum compression deformation is 0.1702 mm, and the maximum tension is 0.252 mm. When there are four joints, the maximum compression is 0.166 mm, and the maximum decompression is 0.184 mm. Therefore, increasing the number of joints in the same length immersed tunnel can improve the seismic performance of the joint.

6. Conclusion

This paper presented a series of shaking table arrays tests that were performed on a microconcrete tunnel model embedded in the soil. The prototype tunnel was constructed under the Pearl River in China. The soil was made of clay mixed with sawdust. The joints between the tunnel segments were simplified by a rubber ring, and a box was designed to contain the soil and tunnel. The tests took into account the uniform and wave passage effect in different apparent wave velocity of longitudinal seismic excitation. Compared to the existing research, the joints were redesigned in this test, so as to research the response of segment joint force and displacement under nonuniform seismic excitation. The results showed that the tunnel and the soil maintained synchronous motion under uniform seismic excitation. However, when the wave passage effect was considered, the tunnel and soil experienced sliding at the interface. The longitudinal wave passage effect and its input direction can, therefore, have a great impact on the axial force, bending moment and deformation of joints. The wave passage effect will make the joint deformation tend to nonuniformity. The comparison shows that structural response under nonuniform seismic excitation is larger than that under uniform excitation. Therefore, the effect of nonuniform seismic excitation should be considered in the design of immersed tunnels.

The simplified model was established in ABAQUS for numerical analysis. The soil around the tunnel was simplified as spring damper, and the tunnel was taken as the dynamic foundation to calculate the dynamic stiffness coefficient of the spring. The tunnel was simplified as a beam element, and the joint was simulated by a nonlinear spring. The numerical simulation results were in good agreement with the experimental data. The immersed tunnel was analyzed by changing the input apparent wave velocity, joint stiffness and joint number. The results showed that the deformation of joints was smaller, and the deformation of flexible joints was greater under high apparent wave velocity. The increase in the number of joints reduced the deformation of joints. It is suggested that the response of different segment lengths under earthquake should be considered in the design of immersed tunnel, so as to determine the best segment length.

Although this research was prompted by the needs of a specific project, many of the conclusions in this study are sufficiently general and may apply to the design of other immersed tunnels.

Data Availability

Data are available upon request.

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Acknowledgments

This project was supported by the National Natural Science Foundation of China (52178392), the S & T program of Hebei Province in China (21375403D, 20375410D), the Natural Science Foundation of Hebei Province in China (E2020210017), and Beijing Natural Science Foundation (8192047).