Abstract
Simultaneous multiple fracturing is a key technology to facilitate the production of shale oil/gas. When multiple hydraulic fractures propagate simultaneously, there is an interaction effect among these propagating hydraulic fractures, known as the stress-shadow effect, which has a significant impact on the fracture geometry. Understanding and controlling the propagation of simultaneous multiple hydraulic fractures and the interaction effects between multiple fractures are critical to optimizing oil/gas production. In this paper, the FDEM simulator and a fluid simulator are linked, named FDEM-Fluid, to handle hydromechanical-fracture coupling problems and investigate the simultaneous multiple hydraulic fracturing mechanism. The fractures propagation and the deformation of solid phase are solved by FDEM; meanwhile the fluid flow in the fractures is modeled using the principle of parallel-plate flow model. Several tests are carried out to validate the application of FDEM-Fluid in hydraulic fracturing simulation. Then, this FDEM-Fluid is used to investigate simultaneous multiple fractures treatment. Fractures repel each other when multiple fractures propagate from a single horizontal well, while the nearby fractures in different horizontal wells attract each other when multiple fractures propagate from multiple parallel horizontal wells. The in situ stress also has a significant impact on the fracture geometry.
1. Introduction
In the past two decades, a great number of researchers have been attracted to improving hydraulic fracturing techniques, which provide large surface area contact with formations and facilitate the production of oil and gas, as evident from the experience in North America [1]. Simultaneous multiple fracturing (simul-frac) is a key hydraulic fracturing technology that fractures two or more parallel horizontal wells simultaneously, which can significantly increase the volume and surface area of the flow path to improve the productivity and recovery efficiency [2–4]. A great number of attempts have been made to simultaneously fracture multiple adjacent horizontal wells to generate complex fracture networks. Although field practices have yielded a significantly improvement from simul-frac over stand-along wells [5], microseismic data [6], and numerical simulations [3, 7–12] also show a complex fracturing network produced by simul-frac, the reasons for simul-frac success are still not clear.
Multiple hydraulic fractures treatment is a complex process, as it needs to consider not only the hydraulic fracturing process but also fracture interaction between multiple fractures. The hydraulic fracturing process is a typical hydro-mechanical-fracture (H-M-F) coupling problem, in which the following three basic processes require discussion [13]: (i) the rock deformation induced by the fluid pressure acting on the fracture surface, (ii) the fluid flow in the fracture, and (iii) the fracture growth. The significant consequences of multiple fractures interaction are stress-shadow effects that lead to stress field and fracture geometry changes.
With the development of computers, numerical techniques have become reliable and convenient tools for studying hydraulic fracturing processes. The finite element method (FEM) [14] is a well-established method for rock engineering problems, and it has been utilized to simulate hydraulic fracture propagation for several decades [15]. However, meshes in the FEM and the fracture surfaces need to be coincident with each other. Therefore, advanced remeshing technologies have been developed for fracturing simulation [16, 17]. FEM with advanced remeshing technologies is still difficult to implement for field-scale reservoirs because of the meshing burden. The extended finite element method (XFEM) [18–20] and the generalized finite element method (GFEM) [21, 22] have been proposed to overcome the meshing burden associated with FEM and applied to model hydraulic fracturing. However, XFEM and GFEM still have difficulties in representing multiple complex fracture networks. The boundary element method (BEM) is also used for hydraulic fracturing problems [23–25]. However, it is difficult for BEM to address interfaces when they are in close proximity [26]. The discrete element method (DEM) [27–29] is well suited for the simulation of hydraulic fracture propagation in rock mass because it inherently incorporates fabric features. However, fracture propagation cannot be simulated directly but is often presented as a function of bond breakage between DEM particles. Furthermore, there is initially interspace between particles regardless of the particle arrangement. The combined finite-discrete element method (FDEM) [30] can well model the transition from continuous to discontinuous behavior in rock mass. FDEM inherits the advantages of FEM in describing elastic deformations and the capabilities of DEM in capturing discontinuities. FDEM is widely used for model collision, fracturing, and fragmentation where both discontinuum and continuum are involved [31, 32]. By coupling it with a fluid flow solver, FDEM has also been applied in handling hydromechanical (H-M) coupling problems [33–36].
Although FDEM has been well-developed and applied by many researchers to study the process of hydraulic fracturing, very little literature can be found about the study of stress-shadow effect with FDEM. In this paper, FDEM coupled with a fluid algorithm, referred to as FDEM-Fluid, is developed to handle simultaneous multiple hydraulic fractures propagation. The organization of this paper is as follows. In Section 2, the fundamental principles of FDEM are briefly presented. In Section 3, the fluid algorithm in fracture networks is presented. In Section 4, the coupling process and a flowchart of FDEM-Fluid are proposed. In Section 5, verification examples are presented, with comparisons with analytical or experimental results. In Section 6, simultaneous multiple fractures propagation in single and multiple horizontal wells under the stress-shadow effect are investigated. Finally, conclusions are drawn in Section 7.
2. FDEM for Solid Phase
2.1. Principles of FDEM
As shown in Figure 1, the computational model is discretized with 3-node triangular elements and 4-node joint elements, which are inserted between the triangular elements. The stress state and deformation of the bulk material are captured by the linear-elastic triangular elements and the nonelastic joint elements. Fracture initiation and propagation within the continuum are simulated by the breakage of the joint elements on the basis of the combined single and smeared crack model [39] when the slip or opening distance is significantly large. The failure criterion will be explained in great detail in Section 2.4.
2.2. Governing Equation
The movement of discrete bodies and the nodal coordinates are updated by the following governing equation:where and are the lumped mass and damping diagonal matrices of the system, respectively. is the vector of the nodal coordinates. , , and represent the internal resisting forces, the external loads, and the contact forces, as shown in Figure 2. The internal resisting forces include the contribution from the elastic reaction forces and the interelement bonding forces. The external nodal forces are applied boundary condition. To model quasistatic phenomena by dynamic relaxation or nonlinear material behavior, numerical damping C is used in the governing equation to account for energy dissipation. The matrix C is calculated aswhere represents the constant damping coefficient and I is an identity matrix.
2.3. Contact Detection and Interaction
A highly efficient No Binary Search (NBS) algorithm [40] is used to detect all the contact couples that are in contact. As soon as contacting couples are detected, contact interaction is performed to calculate forces between discrete bodies by a potential function method. The contact function is given as follows:where is the penalty term, is the boundary of the contactor element, and and are position vectors of the points falling on the target and contactor elements, respectively.
A distributed contact force is generated according to the shape and size of the overlap between the contactor and target elements. As the contactor penetrates an area into the target (Figure 3), a contact force () is generated as follows:where and are overlapping points of the contactor and target, respectively; is the corresponding potential function. The contact force only depends on the endpoints and is independent of the path taken.
2.4. Material Failure
Failure of rock material is simulated by the combined single and smeared crack model [39], which is developed from nonlinear elastic fracture mechanics (NLEFM) [41, 42]. According to this theory, when the strength of the material is reached, a fracture process zone (FPZ) characterized by nonlinear behavior starts to form ahead of the fracture tip (as shown in Figure 4(a)). The material in the FPZ, albeit damaged, is still able to transfer load across the fracture surfaces. In FDEM, fracture surfaces are assumed to coincide with the edges of the triangular elements, and the relative fracture displacements are evaluated by joint elements, as illustrated in Figure 4(b).
(a)
(b)
Based on the local stress and deformation field, fractures initiate and propagate according to three different failure modes, which are often denoted as Mode I (opening failure), Mode II (sliding failure), and Mixed Mode I-II. The constitutive behavior of the joint element is illustrated in detail in the following.
Mode I: as shown in Figure 5(a), before the fracture initiation, the stress at the fracture tip is linear with the fracture opening. Then, the fracture initiates when the value of the fracture opening, , reaches a critical value, , which corresponds to the tensile strength . After the fracture propagates, the normal stress does not disappear immediately but decreases with increasing fracture opening. A stress-free surface (a real fracture) is generated when the opening reaches a certain residual value, .
(a)
(b)
(c)
Mode II: as shown in Figure 5(b), the tangential bonding stress between the two fracture surfaces is a function of the slip distance and the normal stress on the fracture surfaces. The critical slip () corresponds to the intrinsic shear strength (), which is calculated according to the Mohr-Coulomb criterion as follows:where is the normal stress acting on the fracture surfaces; and represent the intrinsic intact material friction angle and the internal cohesion, respectively.
Upon undergoing critical slip, , the tangential stress is gradually reduced to a residual value, , which corresponds to a pure friction resistance:where is the residual friction angle.
Mixed Mode I-II is defined by a coupling criterion between Mode I and Mode II (see Figure 5(c)). This mode describes a combination of shear and tensile failure and is defined by the following:where is the opening distance and is the slip distance.
The values of the residual opening and slip depend on Mode I and Mode II fracture energy release rates, and , which are defined as follows:
3. Fluid Flow in Fracture Networks
As shown in Figure 6, two sets of grid systems, namely, the solid grid and the hydraulic grid, are introduced to represent the solid phase and the fluid phase, respectively. The intact rock mass is discretized with solid grids by using 3-node triangular elements (in gray) and 4-node joint elements (in red). The fluid domain is discretized with hydraulic grids by using 2-node hydraulic elements (in blue), which are inserted into the broken joint elements.
3.1. Basic Assumptions of Fluid Flow
The hydraulic fracturing process is complicated, especially in regard to H-M-F coupling. Therefore, the following hypotheses are introduced: the rock mass is impermeable and the fluid flows only in fractures. In the FDEM model, the flow only happens in the broken joint elements; the fluid is assumed to be laminar, viscous, and incompressible; and ) the fracture apertures can be obtained directly through the displacement of broken joint element surfaces.
3.2. Governing Equations
A cubic law based on parallel-plate flow model [43] is used to simulate the fluid flow in hydraulic elements, in which the flow rate can be expressed as where is the equivalent aperture, is the fluid pressure gradient, and is the viscosity of the fluid.
3.3. Update of Fluid Pressure and Saturation in the Fracture Network
Hydraulic elements connected with the hydraulic boundary are defined as the hydraulic fracture network, which is the only flow path for the fluid flow. Taking hydraulic nodes as an example to illustrate the fluid flow in fracture network, the influence domain of node is the union of all hydraulic elements connected to node (see Figure 6). is a hydraulic element connected to node . and are fluid pressures at nodes and , respectively. The pressure gradient between and can be given bywhere and are the vertical coordinates of hydraulic nodes and , respectively. is the length of hydraulic element .
The flow rate from to can be calculated as
Considering the fluid flow saturation, the actual flow rate from to will decrease when the hydraulic node is unsaturated. Equation (11) can be modified aswhere is a function of saturation, and ; is the saturation of the node. It should be noted that if , ; if , saturation should be updated as follows: where is the saturation of node at the previous time step. and are the volumes of influence domain of this node at the current time step and previous time step, respectively. Then, , and. The volume of a hydraulic node is defined as half of the total volume of all hydraulic elements connected to the node.
The total flow rate of node can be calculated as the summation of the flow rates in the influence domain of node . The fluid pressure of hydraulic node can be updated as follows [44]:where is the pressure at the previous time step; is the bulk modulus of the fluid; is the total flow rate; and is the time increase.
Note that there is a convergence criterion herein; that is, the hydraulic time step should be less than a critical value:where is the node volume and is the permeability factor of the th flow path. For all nodes, the minimum value of is used in the algorithm.
4. Hydro-Mechanical-Fracture (H-M-F) Coupling Procedure
The hydro-mechanical-fracture (H-M-F) coupling process is illustrated in Figure 7. Fluid flows in the fracture network, which is the only path of fluid flow. The fluid pressure applied on the fracture surfaces influences the fracture aperture, which in turn affects the fluid flow behavior. Furthermore, the fractures may propagate and connect with each other under the drive of fluid pressure, which improves the conductivity of the rock mass.
4.1. Coupling Procedure of FDEM-Fluid
This FDEM-Fluid consists of the following two parts: the solid solver, which calculates the deformation of rock mass and fracture propagation, and the fluid solver, which calculates the fluid diffusion process and fluid pressure in fracture networks. The schematic diagram of the coupling between the two solvers is presented in Figure 8.
A typical time step in FDEM-Fluid is illustrated in Figure 9. First, nodal forces and elemental deformation are calculated based on the deformation of the triangular elements and joint elements. Second, fracture initiation and propagation are detected and marked, and, then, the hydraulic elements are updated when the joint elements are broken. Third, the current apertures of fractures are calculated. Fourth, the fluid pressure of the nodes is computed and applied on the triangular element as nodal forces. Then, a new loop starts with the updating of the nodal force.
4.2. Fluid Pressure Acting on the Fracture Surfaces
As shown in Figure 10, the hydraulic pressure, simplified as linearly distributed loadings, is acting on the triangular element boundaries along the fracture surfaces. Then, the fluid pressure is converted to the nodal forces of the relevant triangular elements in the context of FDEM.
and are the hydraulic pressure on hydraulic nodes and , respectively. The total force on this surface can be calculated aswhere and represent the normal force and shear force, respectively, is the length of the fracture surface, and is the equivalent aperture.
The normal force on each node is perpendicular to the edge and can be calculated as
The shear force on each node is opposite to the fluid flow and can be calculated as
4.3. Equivalent Aperture
The fracture aperture can be obtained directly through the displacement of the broken joint element surfaces. Considering the H-M coupling, the fracture aperture is assumed to be a linear distribution along the fracture surface, as shown in Figure 10. and represent the aperture of the two endpoints of the broken joint element. According to Iwai [45], the equivalent aperture of the fracture is as follows: where , and. It is obvious that the equivalent aperture is less than.
5. Validation of FDEM-Fluid
Two numerical tests for the present FDEM-Fluid are carried out, and the results are compared with previous analytical and experimental results.
5.1. Validation 1: Unsteady Flow in a Single Fracture
This test investigates the fluid flow in a single fracture by FDEM-Fluid. The mesh sensitivity of fluid flow in FDEM-Fluid is also discussed. The geometry and boundary conditions of fluid flowing through a single fracture are shown in Figure 11(a). The rock sample is 1 m in length and 0.2 m in height. The right end of the sample is an impermeable boundary, and a hydraulic pressure MPa is applied to the left side. The aperture of the fracture m is constant during the simulation. The hydraulic parameters include the bulk modulus of water GPa and the viscosity of water Pa·s.
(a) Sample geometry
(b) Three sets of mesh with mesh sizes of 0.1 m, 0.05 m, and 0.02 m
The hydraulic pressure distribution with time and location is investigated. The analytical solution of this example is as follows [46]:where; is the distance to the left side; and is the fluid pressure.
Three sets of meshes, as shown in Figure 11(b), with mesh sizes = 0.1, 0.05, and 0.02 m, are utilized to investigate the accuracy of FDEM-Fluid and the mesh sensitivity of fluid flow.
The hydraulic pressure distributions at different time with the three sets of meshes are shown in Figure 12. Note that only hydraulic nodes have fluid pressure. Fluid flows from the left to right side, and the hydraulic pressure decreases from left to right. The pressure at the same location increases with the passing of time.
A comparison between the analytical and numerical solution is presented in Figure 13. A nonlinear decrease of fluid pressure along the fracture can be observed. Fluid pressure increases with the passing of time. The numerical results of the three sets of meshes all agree well with the analytical solution; therefore, the accuracy and reliability of FDEM-Fluid for simulating fluid flow in fractures are validated.
The mesh sensitivity of fluid flow is shown in Figure 14. Comparing the error between the analytical solution and the numerical solution for the three sets of mesh sizes, it can be found that the error decreases when the mesh size decreases. In other words, the lower the mesh size is, the more accurate the numerical solution.
(a) Error analysis at = 0.6 ms
(b) Error analysis at = 2.4 ms
5.2. Validation 2: Hydraulic Fracturing of a Cubic Sample with Preset Fractures
An experiment designed by Jiao et al. [37] is used to validate this FDEM-Fluid in simulating hydraulic fracturing. A cube plaster sample with preset fractures and a central injection hole is prepared, which is then fractured by injecting high-pressure water into the borehole with no other loads applied. In addition, Chen et al. conducted a similar physical experiment [38] where an anisotropic confining stress is applied.
The geometry and cross section of the sample are shown in Figure 15. Two preset fractures with a 30° inclination angle are symmetrically set along the diameter direction of the borehole. The length and depth of the preset fractures are 20 and 90 mm, respectively. The numerical mesh with 11502 triangular elements is shown in Figure 15(c).
(a) Sample geometry
(b) Sample cross section
(c) Model mesh
Two models are taken in to account: in model I, the confining stress is ignored; in model II, the compressive confining stresses MPa and MPa are applied in the and directions. A constant pressure, MPa is applied on the wellbore surface uniformly. The mechanical parameters and hydraulic parameters are listed in Table 1.
The fracture geometry evolution and fluid pressure distribution are shown in Figures 16(a) and 17(a). The displacement magnitude contours of the two models are shown in Figures 16(b) and 17(b). Before fractures initiation, fluid flows into the preset fractures, and the fluid pressure increases gradually. With the passing of time, the prefractures start to propagate from the fracture tips, and fluid flows into the newly generated fracture. The plaster samples are eventually cut into two symmetrical blocks. Note that the prefractures start to propagate at step 3000 in model I and at step 4000 in model II. The plaster samples are cut into two blocks at step 12000 and step 15000 in model I and model II, respectively. This means that the fracture initiation and propagation in model II are harder than those in model I.
(a)
(b)
(a)
(b)
Comparing the two fracture geometries, the fractures propagate strictly along the direction of preset fractures in model I. However, in model II, the fractures obviously turn toward direction, which is parallel to the direction of the maximum principle stress.
The numerical result of model I agree well with the experimental result conducted by Jiao et al. [37], as seen in Figure 18. The result of model II is also similar to the existing results (Figure 19), where the preset fractures gradually reorient to the direction of maximum principal stress which is consistent with the observation by Chen et al. [38], Almaguer et al. [47], and Wang [48].
(a) Top view
(b) Sectional view
These two examples validate the ability of FDEM-Fluid to model hydraulic fracture correctly with and without considering the influence of in situ stress.
6. Simultaneous Multifracture Treatments
Simultaneous multiple hydraulic fracturing in combination with horizontal drilling has been the key technology in the exploitation of unconventional gas/oil. This technology can significantly increase the volume and surface area of the flow path, so as to improve the productivity and recovery efficiency. For simultaneous multiple-fracture propagation, fracture interaction (i.e., the stress-shadow effect) plays an important role in fracture geometry. In this section, two examples are studied to investigate fracture interactions during simultaneous fracturing processes.
6.1. Simultaneous Multifracture Propagation from a Horizontal Well
In this section, a single horizontal wellbore with two parallel initial fractures is simulated and compared to the results of Wu et al. [10]. As shown in Figure 20, two preset fractures are orthogonal to the wellbore axes, with the offset between the two parallel preset fractures being 10 m. A constant injection rate m3/s is applied at the middle of the horizontal wellbore. The numerical model with 123466 unstructured meshes is shown in Figure 20(b). The mechanical parameters and hydraulic parameters are listed in Table 2 [10].
(a)
(b)
This model is simulated under isotropic and anisotropic in situ stress conditions. The differential stress for the anisotropic stress is 0.9 MPa, with the maximum horizontal stress oriented in the direction.
Figures 21 and 22 show the fracture geometry and hydraulic pressure distribution of fracturing in the two models. The fluid first flows into the preset fractures from the borehole, and fluid pressure is applied on the fracture surfaces. After step 4000, the two fractures start to propagate from the fracture tips, driven by the fluid pressure. Then, fluid flows into the newly generated fractures.
(a) Step 4000
(b) Step 8000
(c) Step 12000
(d) Step 16000
(a) Step 4000
(b) Step 8000
(c) Step 12000
(d) Step 16000
As shown in Figures 21 and 22, in both cases, two preset fractures propagate away from each other with a side-by-side configuration because of the mechanical fracture interaction (i.e., stress-shadow effect). Compared to the isotropic in situ stress condition, the fractures under the anisotropic in situ condition have smaller curvature because of the confining stress (see Figure 22).
As shown in Figure 23, the fractures under the anisotropic in situ stress condition propagate longer in the direction than those under the isotropic in situ stress condition, while the influence domain in the -direction is smaller. The complex influence domain is critical to the exploitation of unconventional gas/oil because it influences the surface area contact with the formation or fracture network.
(a) Influence domain of fractures in the -direction
(b) Influence domain of fractures in the -direction
The reason fractures turn is because of the mechanical fracture interaction (i.e., the stress-shadow effect), which exerts additional stress on the surrounding rock, resulting in a local change in direction of the maximum horizontal stress and a deviation of the fracture path.
As shown in Figures 24 and 25, there is a domain surrounding the two adjacent hydraulic fractures, where the stress-shadow effect obviously changes the stress field around the fracture tip and the direction of fracture propagation. Under the anisotropic in situ stress condition, the anisotropic in situ stress can reduce the influence of the stress-shadow effect. In other words, anisotropic far-field stresses tend to suppress fracture curving [49], which indicate that the stress shadow and anisotropic stress compete with each other.
(a) Step 4000
(b) Step 8000
(c) Step 12000
(d) Step 16000
(a) Step 4000
(b) Step 8000
(c) Step 12000
(d) Step 16000
As shown in Figure 26, the FDEM-Fluid results agree well with previous numerical results [10] in both isotropic and anisotropic in situ stress conditions. Therefore, this FDEM-Fluid method is able to simulate the multiple hydraulic fracturing, with considering the H-M-F coupling, the stress-shadow effect, and in situ stress.
(a) Under isotropic in situ stress field condition
(b) Under anisotropic in situ stress field condition
6.2. Simultaneous Multifracture Propagation from Multiple Horizontal Wells
In this example, simultaneous propagation of two sets of initial fractures from two horizontal wells is investigated. As shown in Figure 27, the horizontal wells are 4 m apart, and two preset fractures are orthogonal to each wellbore. The offset between the two parallel perforations in the same wellbore is 4 m, and the offset between two fractures in different horizontal wells is 1.5 m. The fluid pressure in the horizontal wellbores is constant at 10 MPa, and the fluid viscosity is taken as 0.001 Pa·s. The properties of the rock mass are the same as in Table 2.
The evolution of the preset fracture geometry and the distribution of fluid pressure at different time steps are shown in Figure 28. The fluid flows into the preset fractures. Then, the fractures start to propagate, driven by the fluid pressure. It can be observed that fractures and which propagate from each horizontal well turn toward each other because the induced stress-shadow effect changes the maximum principle stress at the fracture tips (see Figure 28(c)). The phenomenon that hydraulic fractures from different horizontal wells attract each other is consistent with the observation by Wang [4], Kumar and Ghassemi [50], and Wu [51]. Compared to fracture , fracture has a larger deviation because of the stress-shadow effect induced by both fracture and . However, and propagate away from each other.
(a) Step 4000
(b) Step 8000
(c) Step 12000
(d) Step 16000
Figure 29 shows the distribution of the maximum principle stress and displacement magnitude. Figure 30 shows the direction of the maximum principal stress map (the green line segment denotes the direction) at steps 4000 and 8000. As the fracture is propagating, it alters the local in situ stresses around it, and the direction of the maximum principle stress can be reversed close to the fracture surface. When the hydraulic fracture propagates into the “stress alteration zone” that induced by the other fracture, the fracture reorients itself to propagate along the new direction of the local maximum principal stress.
(a)
(b)
(a)
(b)
7. Conclusions
The interaction (i.e., the stress-shadow effect) during simultaneous multiple-fracture propagation plays an important role in fracture geometry, which can significantly affect the productivity and recovery efficiency. Therefore, understanding and controlling the interaction effects between multiple fractures are critical. In this study, a hydraulic fracturing simulator, named FDEM-Fluid, based on FDEM and the parallel-plate flow model is presented to capture the development of complex hydraulic fracturing propagation. This FDEM-Fluid is verified by comparing it with existing analytical and experimental results, and good agreements are observed. Then, the multifracture treatment in the hydraulic fracturing process is investigated with FDEM-Fluid. The main conclusions from the study include the following:(1)When multiple fractures are orthogonal to the horizontal well and propagate from a single horizontal well simultaneously, fractures repel each other as a result of the stress-shadow effect, resulting in a complex fracture geometry. Furthermore, the anisotropic in situ stress, in which the maximum horizontal stress is orthogonal to the horizontal well, has a great influence on the curvature and length of the fracture geometry because of the competition between the stress-shadow effect and the anisotropic stress.(2)When multiple fractures are orthogonal to the horizontal wells and propagate from multiple parallel horizontal wells simultaneously, the nearby fractures in different horizontal wells attract each other as a result of the stress-shadow effect, which alters the local in situ stresses and the direction of the maximum principle stress. The fracture tips in these fractures tend to reorient themselves to be perpendicular to the main fractures. This leads to the development of a more complex fracture geometry and prevent fractures from penetrating deeper into the opposite wells.
Conflicts of Interest
The authors declare that there are no conflicts of interest related to this paper.
Acknowledgments
This research is supported by the National Basic Research Program of China (973 Program) (Grant no. 2014CB046900) and the National Natural Science Foundation of China (Grant no. 41602296).