© M. Qi et al., published by IFP Energies nouvelles, 2018

1 Introduction

The glutenite reservoirs are generally characterized by low-porosity, low- permeability and deep burial. The core analysis of glutenite reservoir shows that the pore space is intergranular and the radius of pore throat is much smaller than sandstone reservoir (Meng et al., 2010). With all these reservoir characteristics, hydraulic fracturing is inevitable in glutenite field development. Horizontal well fracturing process can effectively utilize the rock characteristics of glutenite reservoirs and the interference effect between micro-fractures generated from each perforation, forming a more complex fracture network, thus increasing the stimulated volume of the reservoir and improving the field oil production. According to the field application experience of the glutenite block in Shengli Oilfield, China, after the horizontal fracturing, the daily well production is 13.7t, much more than 1.2t with conventional straight well fracturing treatment, the utilization rate of reserves is increased to 96.7% as well (Zhao, 2014), which shows that the horizontal well fracturing technology has a good field application prospect in reconstruction of glutenite reservoirs.

In the past several years, a lot of effort has been put in research on the numerical simulation of hydraulic fracture propagation, and several models for predicting the fracture morphology, including PK, PKN, KGD, P3D, and PL3D, etc. have been put forward (Detournay, 1990). In the research on hydraulic fracture propagation, many coupling terms are considered, including the fracture plane deformation caused by fluids pressure, fluid-solid coupling in fluids flow inside the fracture, influence of reservoir heterogeneity, etc. (Guo et al., 2015b, 2016; Olson and Taleghani, 2009; Rui et al., 2018). Currently, numerical simulation of hydraulic fracturing for formation with heterogeneous media and random fracture propagation directions is the hotspot of research (Delorme et al., 2013a, 2013b; Li et al., 2011; Weng et al., 2011). The numerical simulation of complex fracture propagation mainly involves the following five methods: (1) conventional finite element method: fracture propagation is simulated by the principle of damage mechanics. One is cohesive element method (Paul et al., 2017), which does not apply to random fracture propagation. In another method, all elements are assumed to have the same type. According to the judgment criteria, the dynamic fracture propagation can be simulated by changing the material properties of the elements (Faivre et al., 2016). Although this method is affected by the size of elements with lower precision, it is easy to simulate the random fracture propagation; (2) extended finite element method (Belytschko and Black, 2015): this method has a high precision in simulating fracture propagation. When there are a large amount of other heterogeneous media (natural fractures), the results may be difficult to converge; (3) boundary element method (Olson, 2008): it has certain advantages in simulating random fracture propagation, but may have problems in treating the load and flow on fracture surface; (4) meshless method (Moës et al., 2015): this is a new approach to the simulation of fracture propagation, and easy to use for random fracture propagation. However, theoretical maturity and calculation efficiency should be further improved; (5) analytic model (Weng et al., 2011): the simulation precision is lowered due to the use of analytic model and formula for determining the direction of fracture propagation.

A lot of effort has been put in research on the hydraulic fracture propagation in the glutenite reservoirs (Zhao et al., 2007, 2015). Li et al. adopted the finite element method (FEM) to simulate the dynamic process of rock failure in the glutenite reservoirs, conducted research on influence of horizontal geostress difference, gravel size and content on the fracture morphology, and identified five fracturing modes (termination, branching, deflection, penetration, and attraction) of fracture encountering the gravel based on numerical simulation, obtaining some valuable conclusions (Li et al., 2013). Ma et al. conducted the physical experiment in six natural glutenite rock specimens, and quantified the influence of gravel size, geostress difference, and fracturing fluids type on the fracture morphology, width, and breakdown pressure (Ma et al., 2017). Meng et al. applied the large-size true tri-axial hydraulic fracturing system to simulate the hydraulic fracture propagation in the artificial specimens, simulated the fracture propagation in the conditions of presence or non-presence of gravel, different geostress difference, and different gravel sizes, and quantified the influence of gravel on fracture propagation by combining with fracturing treatment curves (Meng et al., 2010). Liu et al. applied the technology of high resolution micro-computed tomography in the glutenite specimen from the hydraulic fracturing true tri-axial physical experiment, obtaining the visual description of fracture propagation in the glutenite specimen (Liu et al., 2004). Based on CDEM method, and referring to the result of the physical experiment by Liu et al, Ju et al. established the 3D model of the glutenite specimen, and conducted the influence of heterogeneity and principal stress difference on the position of fracture initiation and propagation (Ju et al., 2016). In the above study, the rule of fracture propagation in the glutenite reservoirs was simulated in terms of experiment and numerical simulation, but the scale and duration of simulation were limited, and there was little research on influence of gravel and matrix rock mechanical parameters, fracturing fluids displacement and viscosity, gravel size, etc. on the propagation of hydraulic fracture.

In this paper, a damage mechanics method is introduced to establish a 2D coupled Flow-Stress-Damage (FSD) model of hydraulic fracture propagation in the glutenite reservoirs, and the model is resolved through programming. The influence of reservoir principal stress difference, hydraulic fracturing treatment parameters, gravel content, tensile strength, etc. on the hydraulic fracture propagation, is studied, aiming at capturing the profound understanding of horizontal well hydraulic fracture network propagation in the glutenite reservoirs.

2 Mathematical model of fracture propagation

Compared with other hydraulic fracturing numerical simulation models, the FSD coupled model can simulate the progressive failure process of materials. Its calculation method is based on the finite element and statistical damage theory, taking the heterogeneity of the material properties into consideration, and reflecting in the simulation result (Li et al., 2012; Yang et al., 2011). In the process of simulating fracturing, failure treatment was carried out for units that satisfy a given strength criterion, thus achieving the simulation of heterogeneous material fracturing. The FSD model can be applied not only to simulate the dynamic fracturing process of various reservoirs (including but not limited to shale, Guo et al., 2017; glutenite, hard rock, Zhu et al., 2015 and brittle rock, Li et al., 2011; etc.), but also to the study of fracturing in deviated (Zhu et al., 2015) and horizontal wells (Guo et al., 2015a).

Based on the classical Biot fluid-solid coupling equation, the damage mechanics method was used to establish the 2D fracture propagation model with seepage-stress-damage coupling (Biot, 1941, 2004). By self-compiled FORTRAN program, the finite element simulation of hydraulic fracturing was carried out. Through the redefinition of the mechanical parameters of the meshes in the finite element model, the effect of gravel in the sand conglomerate reservoir on the fracture expansion is simulated. The ABAQUS software was used to open the grid file and select grids that need to be set to gravel, and the independent parameter assignment software was used to assign the gravel grid separately. Finally, the other grids were assigned in the computation program. The software “TECPLOT” was used for data analysis and visual processing.

Rock deformation is modeled with the theory of linear elasticity. Coupling of fluid flow and geomechanical deformations are based on the poroelastic theory developed by Biot. The fluid flow inside the fracture is usually simplified to flow along a channel with a very narrow opening, which is represented by a nonlinear partial differential equation that relates the fluid flow velocity with the fracture width and pressure gradient along the fracture. The fracture propagation process is usually considered in the framework of Linear Elastic Fracture Mechanics (LEFM) theory (Tian et al., 2001). To simplify the processes, we used different type of element with their own properties in our model. Fluid element represents the induced hydraulic fracture. It has very high permeability and very low value of Young’s modulus and Poisson’s ratio. Gravel element represents the gravel inside the reservoir which has higher tensile strength and Young’s modulus than normal formation element and has lower permeability and porosity. We have used formation element with formation properties.

2.1 Seepage equation (Chen et al., 1995)

(1)where P is the pore pressure (MPa); μ is the fluid viscosity; C _t is the total compression factor (C _f  + C _l ), C _f is the rock compressibility, C _l is the fluid compressibility; t is the time; K is the permeability (×10⁻³ μm²). This equation is the mathematical equation for unstable seepage of elastic, porous, single-phase micro-compressive fluid.

2.2 Deformation field equation

(1) Equilibrium equation (Yang et al., 2002)

Suppose the deformation is linear elastic. According to the effective stress principle, the equation characterizing the coupling between the stress field of solid and pressure field of fluid can be obtained:(2)where is the tensor of effective stress; P is the pore pressure (MPa); f _i is the external load (MPa); δ _ij is the Kronecker constant; α is the effective stress coefficient ranging between 0 and 1.

(2) Geometric equation(3)where ε _ij is strain tensor; u is displacement.

(3) Constitutive equation (Rahman et al., 2009)

The stress and strain of elastic material are linearly related:(4)

The constitutive relation of damage mechanics is expressed as follows:(5)where D is the damage variable:(6)where E _c is the elastic modulus of the damage material; E is the elastic modulus of the damage material.

(4) Heterogeneity of rock medium

In order to describe the heterogeneity of the reservoir rock, suppose its mechanical properties obey the Weibull distribution, which is defined by the following distribution density function (Weibull, 1951):(7)where σ _c is the mechanical parameter for the elements of rock medium (elastic modulus, Poisson’s ratio). σ ₀ is the parameter related to the average of all element parameters; m is uniformity coefficient, which reflects the uniformity of rock medium.

(5) Permeability-damage coupling equation

When the maximum tensile principal stress of the element reaches the tensile strength of the rock (σ ₃ ≤ −σ _t ), tensile damage occurs (see criteria for tensile failure). The constitutive equation for rocks with elastic brittleness and residual stress is shown as follows:(8)

Permeability k is given by:(9)where σ _t is uniaxial tensile strength; σ _tr is residual strength at initial tensile damage; ε _t0 is the threshold of strain under tensile damage according to criteria for uniaxial tension. When the strain under uniaxial tension reaches ε _t0, the element begins to undergo damage, but does not totally lose stress bearing ability. When the maximum tensile principle strain reaches the specific ultimate strain ε _tu, the element is considered to completely lose the stress bearing ability. The element is completely damaged and undergoes tensile failure, i.e. D = 1. At this time, in the permeability equation, the effective stress coefficient α is equal to 1. ξ (ξ > 1) is the coefficient of sudden jump of permeability. It characterizes the times by which the permeability is increased before and after damage of element under the same stress state. β is coupling parameter, which can be determined by stress-strain permeability test.

Under the action of compressive or shear stress, the Mohr-Coulomb yield criterion can be used as the second criterion of damage threshold (criterion of shear failure):(10)where σ ₁ and σ ₃ are maximum and minimum principal stress; φ is inner frictional angle; f _c is uniaxial compressive strength.

The constitutive equation under uniaxial compressive stress is written as follows:(11)

Permeability k is expressed as(12)where f _cr is residual stress; ε _c0 is the maximum compressive principal strain when the maximum compressive principal stress reaches the level of uniaxial compressive strength.

Combining equations (1)–(12) with initial conditions and boundary conditions, the overall control equation is obtained. The discrete control equations in finite element method (Eqs. (1) and (2)) are used to iteratively solve the two processes of fluid seepage and rock deformation.

3 Model validation

To validate the rationality of numerical model, a comparison of numerical simulation with the result of physical simulation of hydraulic fracturing in the natural glutenite specimens from Ma et al. (2017) is conducted, with same experimental parameters, and similar gravel sizes and distribution morphology. The accuracy of the model is validated by determining whether the fracture propagation morphology and propagation rule of the fracture tip encountering the gravel in the numerical simulation is the same as those in the physical simulation. According to the parameters listed in Table 1, the numerical simulation is operated with Specimen I with larger gravel grain and Specimen II with smaller gravel grain.

As shown in Figures 1a and 1b, in Specimen I with higher content of pebble gravel with sizes <20 mm, the hydraulic fracture morphology is dominated by the single-wing fracture along the horizontal maximum principal stress, with fracture propagation tendency along the interface between gravel and sandstone, but the gravel has little influence on the global fracture morphology, which is due to small gravel size. The result of numerical simulation in Figure 2 (Specimen I) shows that though the fracture deflects locally after encountering the gravel, the fracture generally propagates along the maximum horizontal geostress, forming the relatively straight fracture plane.

In Specimen II with higher content of gravel with sizes between 40 and 100 mm (Figs. 1c and 1d), the fluids are injected where the gravel is located, and the result of fracturing is shown in Figure 1d. Compared with Specimen I, Specimen II has more complex fracture propagation morphology, forming a series of fracturing morphologies such as fracture reorientation and branched fractures, which is affected by large-size gravel. The fracture propagation trajectory in the numerical simulation is shown in Figure 2. Specimen II displays more complex morphology than Specimen I, and after encountering the large-size gravel, the branched fractures are generated inside the gravel, forming complex fracture morphologies, such as deflection penetrating the gravel, and propagation around the gravel.

Comparison shows that the fracture propagation from the result of numerical simulation is basically consistent with that of physical simulation not only in terms of fracture global propagation morphology but also in concrete details, e.g., propagation of fracture encountering gravel (penetrating gravel, propagation around gravel, and attraction of gravel to the fracture) is displayed similarly in both simulations. Thus, the numerical model has higher reliability, and it is applicable in the simulation of hydraulic fracture propagation morphology in the glutenite reservoirs.

4 Results and analysis

The schematic of the simulating model is shown in Figure 3. Both the X and Y direction of the model is of the length of 2 m, where the X direction is the horizontal minimum stress direction, and the Y direction is the horizontal maximum stress direction. The basic model uses a perforation density of 16 perforations per meter, and to minimize the effect of the boundary effect, the perforations were placed symmetrically in the middle of the model. The gravel size is 40–100 mm, evenly distributed in the model. The basic model parameters are mainly derived from logging and core analysis of glutenite reservoirs in Yong’an block of Shengli Oilfield, as partly shown in Table 2. Other mechanical parameters of matrix and gravel used in the simulation are the same as those in the core experiments, as shown in Table 1.

4.1 Influence of horizontal geostress difference

Hydraulic fracture always propagates along the maximum horizontal geostress (σ _ν > σ _H > σ _h), thus the straight double wings more likely occur in the formation with a higher difference between maximum horizontal geostress and minimum horizontal geostress (Δσ = σ _H − σ _h) (Guo et al., 2015b). Under the influence of the gravels, the hydraulic fracture pattern shows greater complexity than other reservoirs especially for those with a lower geostress difference. For the glutenite reservoir with horizontal well fracturing, the stress state of the reservoir leads to different fracture network morphology, as shown in Figure 4.

In Model 4(a) with a stress difference of 5 MPa, main fractures are generated at perforations that located at the edge of the section. In the middle of the perforation section, the microfractures generated from the tip of the perforations generally extend along the X direction and converge to each other. The morphology of the fracture network nearby the well is like square mesh, in this paper we call this kind of fracture network square mesh-like network. The main fractures that extended from the perforations at edge also produced substantial steering under the influence of gravels and the interference effect of multi-fracture. In Model 4(b) with a stress difference of 10 MPa, although there are still some micro-fractures generated from the middle of perforation section propagate towards X direction and intersect with others, the outwardly extending fractures are not confined to those generated by the perforations at the edge. Main fractures that extended outwardly will still occur deflection under the influence of gravel and intersect with each other. The stress difference in Model 4(c) is 15MPa, which can be seen different from the other two models. Although square mesh-like fracture network still exists in the middle of the perforation section, the fractures generated by the upper perforations on the left side of the model are converged into the main fracture with a bigger width and extend along the oblique 45-degree angle, and the lower part form into the main fracture that extends straightly. The area covered by the fracture network in Model 4(c) is relatively smaller than the other two models.

According to the previous analysis, when the stress difference is small, fractures generated from the middle of perforation section will expand along the direction of the minimum principal stress under the influence of the multi-fracture disturbance and the influence of the gravel. When the geostress difference is large, the deflection effect of the main fracture is not that obvious. With the comparison of these models, the number of main fractures and the stimulated reservoir area in Model 3(b) with the stress difference of 10 MPa is the largest, so the application effect of horizontal well fracturing is better in reservoirs that have an appropriate geostress difference.

4.2 Influence of fracturing fluids injection displacement

The fracturing fluids injection displacement is a treatment parameter requiring deliberate concern during the hydraulic fracturing plan design in the glutenite reservoir. In the process of hydraulic fracturing, the fracturing fluids injection displacement influences the hydrostatic pressure inside the fracture and the penetration force on gravel by fracture tips reaching the gravel, and thus it influences the gravel failure, and propagation morphology and length of the fracture. Also, high displacement is conducive to the formation of more complex hydraulic fracturing network. The study of Olson has shown that the complexity of the fracture path is closely related to the high relative net pressure (low horizontal stress difference or high net pressure) (Olson 2008; Olson and Taleghani, 2009). The high net pressure is more related to the induced stress and the interaction between multiple fractures, conducting to the bending of the fracture expansion path, thus forming a complex fracture network.

The expansion of hydraulic fracturing network under three sets of injection displacement is shown in Figure 5. The injection displacement of Model 5(a) is 0.0005 m³ s⁻¹. The fracture width in this model is relatively short, and micro-fractures that extending from different perforations are easy to converge and communicate with each other. But the formation of the main fracture is difficult, and the stimulated area is significantly smaller than the other two models. The injection displacement of Model 4(b) is 0.0009 m³ s⁻¹, which leads to a much more complicated fracture network due to the rise of the net pressure. The confluence effect of the micro-fractures generated from the perforations is weakened concerning smaller displacement, and 2–3 micro-fractures are converged to form an independent main fracture easily, and six of these main fractures can be found in this model. The injection displacement in Model 5(c) is 0.0012 m³ s⁻¹. In Model 5(c), the fracture tip frequently deflects after encountering gravels, resulting in a higher tortuosity of the fracture than the other two sets of models. The fracture network morphology in the middle of the perforation no longer appears in a square mesh-like shape, and the deviation of the main fractures is much more obvious, especially which generated from the leftmost perforation. But the number of main fractures in this model is five, one less than Model 5(b), thus affecting the model’s stimulated area. Therefore, with the increase in injection displacement, the fracture deflection is more obvious, but the reservoir stimulation effect is not necessarily better.

4.3 Influence of fracturing fluids viscosity

The viscosity of the fracturing fluid and the injection displacement affect the relative static pressure of the tip of the hydraulic fracture, which affects the mutual interference between multi-fractures and the expansion after encountering gravels. Also, the high viscosity fracturing fluids benefit in forming higher and wider fracture and reduce the fracturing fluids leak off effectively. The viscosity of fracturing fluids changes obviously as the formation temperature and pressure vary after the injection into the reservoir. Therefore, it is of significance to study the influence of fracturing fluids viscosity on the fracture propagation (Fig. 6).

In Model 6(a), the viscosity of the fracturing fluid is 20 MPa s. Compared with the other two models, the fractures are more prone to propagate around the gravel instead of penetrating them. Due to low fracture net pressure, no main fracture was formed in the middle of the perforation section, and the fracture deviation is not that significant as the other two models. The fracture network morphology of Model 6(b) and Model 6(c) are roughly similar. Fractures of Model 6(c) with a higher fluid viscosity have better connectivity and deflection degree. But the main fracture generated from the left third perforation converged with other fractures faster than that in Model 6(b), thus affecting the stimulated area. Therefore, fracturing fluid viscosity should not be too low when conducting horizontal well fracturing. But high fracturing fluid viscosity will accelerate the gathering of fractures, which is not conducive to develop a high stimulated area.

4.4 Influence of perforation density

When the distance between the perforations is close, the stress concentration and interference effect caused by other perforations and micro-fractures make the stress distribution more complicated, thus influencing the extension of the fracture network. The interaction effect of perforations is conducive to communication and gathering of the microfractures, which makes it much easier to generate main fractures. To study the influence of perforation density on the fracture network in the glutenite reservoir, three models with perforation density of 8, 16 and 24 perforations per meter are established. The simulation results are shown in Figure 7.

In Model 7(a) with eight perforations, the induced stress and interference effects generated by perforations are relatively small, which result in the independence of each microfracture. Unlike the other two models, micro-fractures expand radically along the direction of horizontal maximum principal stress (Y direction) at first. While even if communication and gathering occur under the influence of gravels and the interference effect from other fractures during the process of fracture expansion, still the complexity and deviation effect of the fracture network is not compared with the other two models. As the perforation density increases, the interaction between perforations and micro-fractures gradually increases, which leads to a faster communication and aggregation rate. Micro-fractures form a square mesh-like network as shown in the Model 7(b). Model 7(c) reflects the expansion of the fracture network under high perforation density. The speed of communication and gathering of micro-fractures is much faster, and a square mesh-like network was formed immediately after the injection. The net pressure of the main fracture in Model 7(c) is greater than the others due to the aggregation of more micro-fractures, therefore result in a more obvious deflection.

4.5 Influence of gravel content

Compared with the conventional sandstone reservoir, the glutenite reservoirs have more complex lithology. For example, rock specimens sampling from the glutenite reservoirs in Yong’an Block, Shengli Oilfield show great heterogeneity, and the gravel sizes in the specimen vary between 1 cm and 20 cm, and with gravel content maximum of 64% and the minimum of 2%. A fracturing model of a horizontal well in sandstone reservoirs was established as a control group, and two models of gravel content 10% and 40% were set up to study the effect of gravel content on the formation of fracture network as shown in Figure 8.

Model 8(a) reflects the expansion of the hydraulic fracture network under the influence of interaction between multiple fractures during the horizontal fracturing. Not that all the hydraulic fractures developed at each perforation can expand into the main fracture, but the action of deflection frequently occurs under the effect of interference, finally gathering with other fractures. With the absence of gravel, the main fractures in the model are mostly developed under the direction of maximum principal stress, but deflection caused by the border effect also exists in some fractures. Fracture network shows much more deflection and tortuosity in glutenite reservoir compared with the sandstone group. From the simulation results of Models 8(b) and 8(c), it can be seen that the effect of gravel content on the expansion of the fracture network is mainly reflected in the following tips: the communication of microfractures in Model 8(b) is relatively poor; for two symmetrical perforations, fracture expansion and communication from one side hinder that from the other, therefore no square mesh-like network was formed; Model 8(c) shows a better fracture network complexity and deflection, which lead to a larger stimulated area.

4.6 Influence of gravel tensile strength

During the process of hydraulic fracturing, the most important rock failure form is the tensile failure. The fracture generation and propagation during hydraulic fracturing are dominated by the tensile failure in the rock matrix, thus the tensile strength of gravel is an important element influencing the propagation of hydraulic fracture encountering the gravel. A simulation of propagation trajectory of hydraulic fracture network in the glutenite reservoir is operated with the tensile strength of gravel (σ _tg) 15 MPa, 20 MPa, and 25 MPa.

It is shown in Figure 9 that the influence of the gravel tensile strength on the fracture network is mainly reflected in fracture deflection. The fracture is more prone to deflect in the model with higher strength gravel. The tensile strength of Model 9(b) is 20 MPa, the overall fracture morphology is similar to that in Model 9(a), but the fracture is more prone to deflect when encountering gravels, which results in a fracture network with higher tortuosity. Model 9(c) with gravel tensile of 25 MPa generates a fracture network that propagates toward the direction of X integral. But it is worth noting that the general fracture length, number of main fractures and stimulation area in Model 9(c) are the least of these three models. Therefore, the effect of horizontal well fracturing is better in the reservoir with lower gravel tensile strength.

5 Conclusion

In this paper, the influencing factors of fracturing development of horizontal well in the glutenite reservoir are studied by numerical simulation. The hydraulic fracturing performance was evaluated based on the complexity of the fracture network, fracture deviation, and the stimulated reservoir area, and the following conclusions were obtained.

1. Based on the damage mechanics method, a coupled FSD model of hydraulic fracture propagation in the glutenite reservoirs is established, and the accuracy of the model is validated through comparison between results of the physical experiment and numerical simulation.

2. Square mesh-like fracture network was formed in the middle of the perforation section at a geostress difference of 5 MPa. Fracturing deviation effect is not that obvious with a high geostress difference, and the stimulated area is also relatively lower.

3. The injection displacement raises the relative static pressure of the tip of the hydraulic fracture directly, as well as the viscosity of the fracturing fluid, which affects the complexity of the fracture network. Microfractures generated from the perforations are more easily to converge into main fractures with a higher displacement, and fracturing deflection after penetrating gravel occurs more often. Although fracturing fluid viscosity should not be too low when conducting horizontal well fracturing, high fracturing fluid viscosity will accelerate the gathering of fractures, which is also not conducive to develop the high stimulated area.

4. The increase of perforation density mainly enforces the interaction between the perforations and the microfractures, which can effectively promote the communication and connection. In the case of 24 perforations per meter, the aggregation rate of microfractures is obviously faster, and the net pressure of the main fracture is greater than the others, therefore result in a more obvious deflection.

5. The main fractures in the model are mostly developed along the direction of maximum principal stress with the absence of gravel, but deflection caused by interference and border effect also exists. With the gravel content increases, the communication speed of microfractures and complexity of the network increase, as well as the fracture network diversion. The fracture shows great deflection and tortuosity in the reservoir with high gravel tensile strength and is more prone to deflect when encountering gravels, but the stimulated area is smaller instead.

Acknowledgments

The authors would like to acknowledge the financial support of the National Natural Science Foundation of China (Grant No. 51404288), the Fundamental Research Funds for the Central Universities (Grant No. 17CX02077) and the Applied Basic Research Project for Qingdao (Grant No. 17-1-1-20-jch).

References

All Tables

All Figures