#### Abstract

A semianalytical solution of stress and displacement in the strain-softening and plastic flow zones of a salt cavern is presented. The solution is derived by adopting the large deformation theory, considering the nonlinear Hoek–Brown (H-B) strength criterion. The Romberg method is used to carry out numerical calculation, and then, the large deformation law of displacement is analyzed. The results are compared with those obtained by former numerical methods, and the solutions are validated. The results indicate that the displacement of the plastic zone decreases with the increase in distance away from the salt cavern. Similarly, it decreases with an increase in the geological strength index or running pressure, with the running pressure having a more significant effect on the displacement. It increases with the dilation angle, and the impact degree gradually increases. Compared with the softening parameter, , of the plastic zone, the flow parameter, , has little impact on the displacement. The displacement of the plastic zone obviously increased when considering the strain-softening of salt rock. When considering the shear dilation and softening behaviors of salt rock, the analytical solution obtained by employing the experiential regression Hoek–Brown (H-B) criterion, which considers many factors such as the structural characteristics of the salt formation and the rock mass quality, is safer and closer to the actual situation. This study can provide reference for many applications, including but not confined to analyzing the deformation of the surrounding rock of an underground salt cavern storage facility during construction.

#### 1. Introduction

Salt caverns are considered the best places in the world for energy reserves [1–5]. The reason is that salt cavern storage facilities, compared to other underground storage facilities, have many advantages. For example, it has simple hydrological conditions, a wide distribution range, and a complete structure and stratum status. Moreover, the mechanical properties of rock salt, compact structure, low porosity and permeability, and damage self-recovery capability can ensure the safety of salt cavern. An underground salt cavern storage facility can be constructed using a water solution to form underground adits in deep salt formations. However, the redistribution of the stress and strain around the cavity and the damage to the mechanical properties of the surrounding rock are closely linked with the solution process. In recent years, seal failure, volume shrinkage, excessive deformation or heavy collapse, serious ground subsidence, etc. have become among the most arduous challenges to the engineers. For example, in the Hutchison salt cavern gas storage in Kansas, USA, a ground subsidence with 90 m in diameter was observed. Internal pressure caused excessive volume shrinkage, and the Eminence salt cavern storage in Mississippi, USA, lost about 40% of its volume in 2 years [6]. Thus, analysis of the stability of underground salt cavern gas storage is very important for controlling this hazard.

In order to ensure the safety and stability of underground adits, a large number of analytical methods for the stress and deformation of the surrounding rock have been suggested for ideal underground spherical caverns in the past. For stability analysis of deep underground cavity, the Hoek–Brown (H-B) failure criterion has been found to be more suitable than the commonly used Mohr–Coulomb (M-C) criterion. The H-B failure criterion was proposed for hard rock masses and very poor qualities of rock masses. It can reflect the inherent nonlinear damage characteristics of rock masses. Using the nonlinear H-B failure criterion, Brown et al. [7] derived a solution for the stress and displacement of the plastic zone of a tunnel. Using the nonassociated flow rules and considering three cases for the elastic strain of the plastic zone, Park and Kim [8] put forward a solution for the displacement of the plastic zone of a circular tunnel. Employing the H-B criterion, Lee and Pietruszczak [9] analyzed the stress and deformation field of circular adits with strain-softening behavior and presented a solution for the surrounding rock of a tunnel by assuming that the criterion parameters were linearly transformed with the shear strain. Introducing the nonassociated flow rules in the design of a hydropower tunnel, Rojat et al. [10] derived a closed-form solution for the surrounding rock of a tunnel using the H-B criterion. Park [11] dealt with a similarity solution for a spherical or circular adits excavated in elastic-strain-softening rock mass compatible with the M-C or H-B yield criterion. The partial differential equations from stress equilibrium, constitutive law, and consistency equations are simplified with first-order ordinary differential equations that solved via the Runge–Kutta (R-K) method. In China, Liu and Lin [12] derived the elastoplastic deformation of a tunnel by adopting the M-C criterion. Adopting the H-B strength theory and limit equilibrium equation, Yao et al. [13] and Pan et al. [14] conducted elastoplastic analyses of the surrounding rock in an idealized roadway with hydrostatic pressure. Wang et al. [15] proposed a brittle-plastic solution of a spherical cavity, and the solution is used in the strain-softening analysis. They also proposed numerical solutions and studied their application to determining stress–strain curve of strain-softening H-B rock masses. Wang et al. [16] proposed an elastic analytical solution for a spherical salt cavern with a lateral pressure coefficient (*λ*) and used the H-B failure criterion to derive the internal pressure limit when the plastic failure of the gas storage cavern did not appear. In order to study excavation problems for tunnels exhibiting strain-softening behavior, considering the effect of the confining stress, Wen and Yang [17] and Cai et al. [18] conducted numerical calculations by adopting the R-K method. They obtained the deformation around a tunnel and analyzed the displacement of the surrounding rock of a deep underground chamber by using the H-B criterion. Using the elastoplastic theory, Wang et al. [19], Fan et al. [20], and Hu and Zhu [21] presented analytical solutions to the expansion of spherical cavities and analyzed the effects of softening and shear dilation on the spherical expansion. Adopting the M-C criterion and nonassociated flow rule and introducing the large deformation theory, Wang and Zhu [22] obtained an elastoplastic analytical solution to the cavity expansion. Zhang et al. [23] established a multistep brittle-plastic model by simplifying the strain-softening process of the postfailure region, and they proposed analytical solutions of a spherical cavern with considering the deterioration of elastic parameters. Both M-C and H-B criteria are included in their analysis. Zhang et al. [24] analyzed the stress and deformation of the surrounding rock in a spherical salt cavern under ideal elastoplastic conditions, and they predicted the creep deformation of the surrounding rock via theoretical analysis.

The analytical elastoplastic solution for a space problem accounting for large deformation requires a complicated mathematical deduction, and, even when it is combined with a numerical method, there are a limited number of studies that can be referenced. Wang and Qian [25] proposed a finite difference method. Adopting the large deformation theory, Guan et al. [26] proposed a numerical procedure that reflects the strain-softening behavior and studied the instability problem of a plane strain axially symmetric tunnel. The above studies mainly involved two-dimensional research on circular tunnels. Jiang et al. [27] provided key information for predicting deformational performance for a tunnel undergoing large deformation problem. Zhang et al. [28] presented a large strain numerical solution of a circular tunnel for the ground reaction curve in strain-softening rock masses. In consideration of the tunnelling-induced large deformation, Mo et al. [29] analyzed the cavity contraction in M-C and H-B media.

In the previous results for salt cavern storage projects, there has been little information on the large deformation [30] properties of salt rock, and these results were obtained mainly by employing the M-C criterion. Accounting for the erosion effect of solution mining on the cavity wall, the H-B criterion, which is based on empirical regression, is especially practical with salt cavern gas storage. Hence, this paper considers the large deformation around a salt cavern by introducing the H-B criterion and taking the strain-softening property of salt rock into account. This study can have significant implications for the design and construction of a deep hollow spherical domain.

#### 2. Mechanics Model for Spherical Salt Cavern Gas Storage

In situ investigations indicated that salt caverns have all irregular shapes; they are different from each other. Most of the salt storage caverns are pear-like in the design stage. However, in order to facilitate the theoretical analysis and derivation, we can establish a spherical symmetry model for the stability analysis of a salt cavern under hydrostatic stress (side pressure coefficient is 1). Salt cavern gas storage is here regarded as a deep hollow spherical domain subjected to internal and external pressures, the radius of the salt cavern is designated as . The gas pressure (approximately homogeneous running pressure), , and hydrostatic pressure, , are the inner and outer boundaries of the model, respectively. Gravity and side pressure coefficient are not considered in the model. In order to meet the demands of spherical symmetry, the coordinate origin of the salt cavern mechanics model is designated in the center of the spherical adits. It agrees with the hypothesis of isotropy during the construction of a salt cavern. There exist elastic and plastic zones around the salt cavern. The far filed of the cavern () enters into an elastic state. The stress of the elastic zone is smaller than the strength of the salt rock; as a result, the deformations can recover after unloading. At , the radial stress and tangential strain are designated as and , respectively. The near-zone of the cavern () enters into a strain-softening state. The rock in the strain-softening zone will undergo microcrack propagation and plastic deformation. At , the radial stress and tangential strain are designated as and , respectively. Likewise, the near-zone of the cavern () enters into a plastic flow state. The rock in the plastic flow zone undergoes ultimate failure. What is more, the salt rock around the cavity satisfies the continuity conditions both in the elastoplastic interface, , and in the interface of the strain-softening and plastic flow zones, . Figure 1 shows the elastoplastic model of the rock surrounding spherical salt cavity.

**(a)**

**(b)**

The rock salt formation where salt cavern gas storage is located has a depth arranging from 1000 to 2000 m, so surrounding rock of the salt cavern is situated in high geostress levels. In situ investigations and compression tests indicated that the large plastic deformation played an unneglectable role during the analysis of the mechanical properties of rock salt under triaxial stress. Thus, the shear dilatation and strain-softening properties of rock salt [30–33] are considered in this study, and a strain-softening mechanical model is used.

#### 3. Strength Criterion

The generalized H-B criterion is an experiential regression strength criterion. The H-B failure criterion was proposed in 1980s originally as follows [34]:where is the uniaxial compressive strength of the intact rock (MPa) and is a dimensionless material constant, with the value of varying from 0.001 to 25.0.

To eliminate some of the deficiencies with the application of the original H-B failure criterion, a modified rock mass failure criterion has been developed and is presented by Hoek et al. [35]. Hoek et al. suggested a specific form of the strength criterion:where is the uniaxial compressive strength of the intact rock (MPa) and , , and are dimensionless material constants.

The H-B failure criterion was updated by discussing the use of the geological strength index (GSI) to describe the structure and condition of a wide range of types of rock masses and the selection of the disturbance factor *D* [36, 37]. The dimensionless material constants , , and can be calculated by (2) [36, 37]. The parameter varies over a range from 0.5 to 0.666, where . Considering the difficulty to obtain theoretical solutions of nonlinear numerical problems, to simplify the calculation, the commonly used H-B strength criterion equation (1b) is used and the parameter is given a value of 0.5 in this paper:where the superscripts “peak” and “res” represent the peak parameter and residual parameter, respectively. Hoek and Brown proposed a nonlinear relation between the major and minor principal stresses at failure, assuming independence of the intermediate principal stress [29]. In recent years, many studies on generalized three-dimensional H-B strength criterion have been carried out [38–41]. The commonly used H-B criterion is modified continuously, and the generalized three-dimensional model is the future direction in the research and application.

#### 4. Solutions in the Elastic Zone

For the symmetry of a spherical salt cavern subjected to homogeneous internal and external pressures, the corresponding relationships between both the stresses and the strains can be written as and . As a result, the salt rock around the cavern satisfies the force equilibrium equation:

The compressive stress is assumed to be positive, the initial in situ stress is not considered, and the stresses and strains in the elastic zone can be obtained, respectively, as follows:

Then, the elastic strains within the elastic zone at can be written as follows:where and .

What is more, the displacement caused by the salt cavern’s construction can be derived as follows:

Assuming that , represents the radial stress at . The values of and are not yet known. When elastic–perfectly plastic model is used, and strength parameters and in the H-B yield criterion are constant values, the values of and can be obtained [24]. However, if strain-softening occurs in the plastic zone and if and in the H-B criterion change after peak strength, numerical methods have to be used to get the values of and/or . Besides, the stresses and strains, on the interface between elastic zone and softening zone and on the interface between softening zone and residual zone, have to be used in the solving process; the detail determination method of the elastoplastic interface and its corresponding stress will be analyzed later. Meanwhile, the stresses and displacements at the elasto-softening interface have to be used to obtain the stresses and displacements on the softening zone and the residual zone.

#### 5. Stress Analysis in the Plastic Zone

##### 5.1. Stress and Radius of the Strain-Softening Zone

Within the near-zone of a salt cavern, , to build the H-B softening constitutive model, it is of critical importance to define a variable as the softening coefficient, which can reflect the linear degradation law of each softening model parameter during the strain-softening process (Figure 2). The coefficient is given as , where is the tangential plastic strain in the strain-softening zone. The functional relation between and is established as follows:where can be used to describe the H-B parameter or ; and are the corresponding peak parameter and residue parameter, respectively, and can be obtained by substituting and into (2), respectively; and represents the residue value of .

According to the definition of the softening coefficient, we can obtain the following:where, the parameter is the curve slope for the strain-softening process. When *K* = 0, the elastic-brittle-plasticity model degenerated to an ideal elastoplastic model [42].

After the rock around the cavern enters the strain-softening state, and vary with the softening coefficient and conform to the H-B criterion. To determine the stress within the strain-softening zone, the classical H-B criterion can be rewritten as follows:

Combining (3) and (10), the radial stress can be derived as follows:

Substituting into (11) at , we can obtain the radial stress of the strain-softening zone as follows:where ; ; and vary with the variation of softening coefficient ; and .

According to (10), the tangential stress can be determined as follows:

Accounting for the boundary condition at , we can obtain the following:where the parameters and are given the residual values.

Because (14) is a nonlinear equation, , the iteration method can be employed to determine the value of in the strain-softening zone at . To determine the radius of the strain-softening zone, adding the stress–strain curve in Figure 1(b) and assuming that the elastic strain in the plastic zone abides by Hooke’s law, we can express the stress in the strain-softening zone as follows:where .

Considering the boundary conditions at and combining (13) and (15), assuming , we can obtain the following:where and are obtained by and at , respectively.

##### 5.2. Stress and Radius in the Plastic Flow Zone

Within the near-zone of a salt cavern, , the stress satisfies (1b) and (3). The H-B parameters and are residual values, which can be derived by combining (2) and parameter . Combining (1b) and (3) results, we obtain the radial stress of the plastic flow zone as

Substituting into (17) at results in the radial stress:where and .

Adding (1b), we can determine the stress of the plastic flow zone, , as follows:

Accounting for the continuous conditions, at , the stress of the plastic flow zone should equal that of the strain-softening zone . Combining (12) and (18), we can obtain the radial stress as

Adding (20), we can determine the range of the plastic flow zone aswhere .

#### 6. Deformation Analysis in the Plastic Zone

According to the plastic theory, the strain in the plastic zone can be divided into the elastic and plastic strains. For salt rock, the large strain theory is more suitable for the plastic strain. Generally when the degree of deformation is less than 10%, it can be expressed by engineering strain; otherwise, logarithmic strain should be used. The deformation of salt rock under triaxial stress state indicates the large deformation characteristics of salt rock under high confining pressure. Using the logarithmic strain to analyze the large deformation characteristics of salt rock is more suitable [30]. Therefore, the strain can be expressed in the logarithmic form as , where and are variables in the calculations. Then, the geometric equation can be rewritten as (22). Although the elastic strain in the plastic zone has various definitions [43], it is assumed that it conforms to Hooke’s law in this study:

##### 6.1. Displacement Analysis in the Strain-Softening Zone

In the zone of , the strains can be stated as follows:where the superscripts “*e*” and “*p*” denote the elastic and plastic strains of the strain-softening zone, respectively.

Employing the flow rule establishes the relationship between the components of the plastic strain. According to the nonassociated flow law, we can get the relationship between the plastic strain increments as follows:where *h* is a function of , , and is the angle of dilation obtained by a triaxial compressive test.

Combining (23) and (24) results in the following:

Based on the integration of (25) and by adding (5), we can derive the following:where can be given as .

After substituting the radial stress of (12) and tangential stress of (13) into the above form for , we can obtainwhere , , and .

Combining (22) and (26) leads to the following equations:

Definite integral calculation can be performed on the intervals and , where represents the displacement around a salt cavern in the elastic zone at . Based on (7), we can obtain .

The integration of (28) leads to the displacement for a large deformation in the strain-softening zone:

Using the method of definite integration by substitution, we can deduce (29) as follows:

Furthermore, the solution of (30) can be derived based on Romberg integration. At , according to (28), the tangential strain at the interface of the strain-softening zone and plastic flow zone can be given as follows:where .

Combining (16) and (31), we can obtain the nonlinear differential equation for . However, this equation has no closed-form solution. To find the values of and , the numerical iteration method is employed with the aid of MATLAB. Next, substituting into (20) produces . Then, substituting into (21), we can obtain and .

##### 6.2. Displacement Analysis in the Plastic Flow Zone

Likewise, in the zone of , the strains can be expressed as follows:

Based on the nonassociated flow law, we can determine the plastic strain increments as follows:where is a function of ; ; and can be got through a triaxial compressive test.

Combining (32) and (33), we will have the following:

Based on the integration of (34) and by adding (5), we can derive the following:where .

After substituting the radial stress of (18) and tangential stress of (19) into , we can obtainwhere ; ; ; and .

Combining (32) and (34) results in the following:

Definite integral calculation can be performed on and , where is based on (30) and represents the displacement in the strain-softening zone at . The integration of (37) leads to the displacement for a large deformation in the plastic flow zone:

Using the definite integration by substitution, we can deduce (38) as follows:

Then, (39) can be solved by the method of Romberg integration.

#### 7. Discussion

When , there is no softening behavior, and the zoned phenomenon does not remain in the plastic zone. In addition, the constitutive relationship of the salt rock is the ideal elastoplastic one. The radial stress and tangential stress in the plastic zone can be written as (40) and (41), respectively:where , , and and are stationary values obtained by substituting the parameter into (2):

Considering the boundary condition at , we can express (14) as follows:

The radial stress at , , can be derived through numerical iteration. Then, plastic radius can be determined according to (40):

When calculating the displacement in the plastic zone, (26) is rewritten as follows:where , , , and .

Combining (22) and (44), we can obtain the following:

Moreover, definite integral calculation can be performed on and , where is based on (7) and represents the displacement in the elastic zone at the interface of the plastic zone and elastic zone. The integration of (45) leads to the displacement for a large deformation in the plastic zone with an ideal elastoplastic rock mass:

Using the method of the definite integration by substitution, we can deduce (46) as follows:

Furthermore, the solution of (47) can be determined based on Romberg integration.

#### 8. Case Study

Large-scale underground salt cavern gas storage construction has been carried out in Huaian, Jiangsu, Pingdingshan, Henan, Wanzhou, Chongqing, Yunying, Hubei, etc., those storage groups will be used as oil and gas reserves. At present, there is just one underground gas storage salt cavern (Jintan gas storage) in China at running condition [44]. To show the accuracy of the proposed solutions, the results of the displacement, the radial and circumferential stresses, and the variation of are compared using the dataset as shown in Table 1. Table 1 lists the salt cavity structural and rock mechanical parameters. , , and denote the elastic modulus, poison ratio, and dilation angle of salt rock, respectively. Table 2 lists the H-B parameters of the salt rock. After the water encroachment of salt cavity storage, the damaged quality and structure of the rock mass decrease the GSI value. Accounting for the empirical regression H-B parameter, GSI affects the displacement, with different GSI values of 60 and 70.

Primarily, through numerical iteration, is calculated on the basis of (14). Furthermore, the nonlinear differential equations (48) are utilized according to (16) and (31). Then, the unknown variables and are determined with the aid of the MATLAB software. Next, , and are calculated one by one. Table 3 lists the related data:where .

First, it is particularly necessary to validate the analytical result. Therefore, the elastic-brittle-plasticity model is simplified to an ideal elastoplastic model. The solution of (47) is derived and compared with the numerical solutions that are carried out with the use of the H-B constitutive model in FLAC3D. Figure 3 shows the principle stress and displacement of the rock masses around the cavity, respectively. It can be seen from Figure 3 that the solutions obtained with different H-B criterion parameters are in accord with the distribution law of deep spherical adits and agree with the published outcomes [24]. The variation in GSI affects the principle stress but not obviously. Without considering softening behavior, the analytical solution is closer to the numerical solution. However, when taking the softening behavior into account, the plastic zone is enlarged, and a great difference exists between the analytical and numerical results, particularly at the cavity wall. The numerical solution for the tangential stress is significantly larger than the analytical solution. It can be seen from Figure 3(b) that GSI markedly affects both the displacement at the cavity wall and the plastic radius of the surrounding rock. For example, the GSI increases by 10 and, as a result, the displacement at the cavity wall decreases from 0.72 m to 0.51 m, while the plastic radii *R*_{1} and *R*_{2} decrease by 6% and 12%, respectively. The influence of GSI on *R*_{2} is more obvious. The numerical result for the displacement obtained when adopting the ideal elastoplastic model is evidently smaller than the analytical result obtained when considering the softening behavior, but it should be noted that the numerical result is fundamentally in accord with the analytical result obtained without taking the softening behavior into account. It can be seen from the results that the variation of GSI and/or the softening behavior greatly affects the results. Thus, if the effects of the two factors are not considered, the conclusions will have large deviations.

**(a)**

**(b)**

Alejano and Alonso [45] suggested a definition of dilatancy; it is a volume change caused by the shear distortion which could not recover. Hoek and Brown [46] proposed that the dilation angle of rocks may be extremely small. A dilation angle of zero shows that the volume is preserved. Three cases, , , and , are considered here. The relationships between the dilation angles and the displacements within the plastic zone of a salt cavern are illustrated in Figure 4. The figure shows that the displacement at the same point increases with the dilation angle, but the degree of this effect decreases with an increase in the distance away from the salt cavern. The influence of the dilation angle on the displacement at the cavity wall is the most striking. Taking the softening behavior into account, when the angle of dilation is , , or , the displacement at the cavity wall is 0.72 m, 0.84 m, or 0.97 m, respectively. Thus, the extent of the influence increases with the dilation angle. We next consider the displacements of the plastic flow zone for cases of , , and when *h* = 0.7699 () and . Figure 5 illustrates that the displacement of the plastic flow zone decreases with an increase in . In other words, the displacement increases with dilation angle , but the variations are minor. Compared to , the effect of on the displacement is relatively small.

Figure 6 indicates the relationships between the operating pressure and displacement at the cavity wall. As we can see from Figure 6, the displacement decreases with an increase in the operating pressure. When the operating pressure is increased six times, both the analytical results and numerical results are decreased by 80%, and the reduction rates are essentially the same. Under a pressure of 2 MPa, an analytical solution of 2.4 m is attained, but when the softening behavior is not considered, the analytical solution and numerical solution are merely 0.55 m and 0.75 m, respectively. Thus, the difference is evident. The variations are caused by the fact that the numerical solution does not take the softening behavior into account, and the definition of the elastic strain within the plastic zone is different. The results demonstrate that the operating pressure has a significant effect on the displacement of the salt cavern storage at the cavity wall. In order to ensure safe operation of the reservoir, the minimum running pressure should not be too low. Thus, it is safer to design the operating pressure based on the analytical solution.

#### 9. Conclusions

(1)Under the condition of three-dimensional model, considering the large strain characteristics of salt rock, the logarithmic strain is introduced to analyze the theoretical solution for the large deformation of the surrounding rock of salt cavity. And a semianalytical solution is obtained via Romberg numerical method. The stress, strain, and large deformation displacement of the strain-softening and plastic flow zones of a salt cavern in H-B rock masses were obtained.(2)The H-B parameter, GSI, was found to have a significant impact on the radius and displacement of the plastic zone. Both the plastic radius and displacement decreased with an increase in the GSI parameter, which indicated that the plastic zone and displacement decreased with an increase in the rock mass.(3)The dilation angle also affected the distribution of the displacement in the salt cavern. The displacement was positively correlated with the angle of dilation. Compared with the impact of the strain-softening parameter, , the flow parameter, , had little effect on the displacement. The dilation angle and GSI affected the displacement of the salt cavern most significantly at the cavity wall.(4)The displacement and plastic zone became obviously greater when the strain-softening behavior was considered, especially for a poor rock mass and lower operating pressure. Therefore, the strain-softening properties cannot be ignored when designing an actual project.

#### Notations

: | Radius of salt cavern |

: | Gas pressure |

: | Hydrostatic pressure |

: | Plastic radius of strain-softening zone |

: | Radial stress at the elastoplastic boundary |

: | Tangential strain at the elastoplastic boundary |

: | Plastic radius of plastic flow zone |

: | Radial stress at the interface of the strain-softening and plastic flow zones |

: | Tangential strain at the interface of the strain-softening and plastic flow zones |

, : | Uniaxial compressive strength of the intact rock |

, , , : | Constants for Hoek–Brown material |

: | Geological strength index |

: | Disturbance factor |

, : | Radial stress and tangential stress of elastic zone |

, : | Radial strain and tangential strain of elastic zone |

, : | Radial stress and tangential stress of strain-softening zone |

, : | Radial strain and tangential strain of strain-softening zone |

, : | Radial stress and tangential stress of plastic flow zone |

, : | Radial strain and tangential strain of plastic flow zone |

: | Young’s modulus |

: | Poisson’s ratio |

: | Shear modulus |

: | Dilation angle |

: | Softening coefficient |

, , : | Displacement in the elastic zone, strain-softening zone, and plastic flow zone, respectively. |

#### 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

This study was supported by the National Natural Science Foundation of China (Grant no. 51504124).