Research Article | Open Access
Xiang-qian Xu, Hao-bin Zhou, "Study on the Influence of Pulse Current Cathodic Protection Parameters of Oil Well Casing", Advances in Materials Science and Engineering, vol. 2019, Article ID 2847345, 9 pages, 2019. https://doi.org/10.1155/2019/2847345
Study on the Influence of Pulse Current Cathodic Protection Parameters of Oil Well Casing
Pulse current cathodic protection had obvious advantages over traditional DC cathodic protection, but when square-wave pulse current was used for cathodic protection of oil well casing, the influence of pulse protection parameters on potential distribution was not very clear. Therefore, according to the current situation of oil fields in China, the typical formation characteristic structure of oil wells was selected and the cathodic protection model of pulse current for oil well casing was established. Based on the established model, the effects of pulse current parameters (frequency and duty cycle) and anode distance on the cathodic protection effect of square-wave pulse current were studied. The results showed that better protection effect could be obtained by choosing appropriate anode distance, high frequency, and moderate duty cycle.
The damage of oil well casing affects not only its normal production but also the normal work of adjacent wells, resulting in damage to the normal production of the whole oilfield . Cathodic protection technology is widely used abroad to protect casing from corrosion . External current cathodic protection of oil well casing is a standard method to restrain corrosion of casing outer wall . This technology is a very economical and effective method to prevent corrosion of casing outer wall . Pulse current cathodic protection has been successfully applied to casing protection of cluster wells in 1968 . Some researchers have proved that pulse current cathodic protection system of oil well casing has more uniform current density, deeper protection depth, and better protection effect than sacrificial anode cathodic protection system and direct current cathodic protection system [6, 7]. However, there are many variable electrical parameters of pulse current, such as frequency, amplitude, duty cycle, and pulse shape, which have no clear influence on the protection effect. Furthermore, the geological environment and stratigraphic structure of the casing are complex, and then the dielectric properties of the cathodic protection current flowing through the casing are different. So at present, only a few oil fields use this technology in a small scale.
The field application results show that the pulse current cathodic protection has more uniform current density, deeper penetration, smaller current, and more uniform potential distribution on the casing surface, which makes the casing have deeper protection depth than the traditional DC cathodic . For example, in the DC cathodic protection of deep well casing, the protection depth can only reach about 1500 m. If the pulse current cathodic protection technology is used, the protection depth can be extended to more than 3000 meters . Metwally introduced a theoretical investigation of the pulse cathodic protection systems to show how they behave under different operating conditions by using a 3D field approach software package. The effectiveness of pulsed current cathodic protection in desert and coastal areas was compared in Metwally’s paper. Metwally gave many suggestions on the application of pulse current cathodic protection. However, there was no suggestion on the application of pulse current cathodic protection for oil well casing with multilayer rock structure . George studied the corrosion rate and anode life of pipeline under steady state and pulsed DC output conditions . At present, the design of pulsed current cathodic protection for casing in oil and gas wells is often empirical or can only reach the level of semiempirical and semitheoretical with the help of cathodic protection theory for long-distance pipelines. Although it is known that the application of pulse current cathodic protection to casing in production wells can extend the protection depth, it is often unclear how the cathodic protection potential is distributed at different depths and how deep each well can be protected. Lack of sufficient understanding of the protection mechanism hinders the application and promotion of this technology. According to the current situation of oil fields in China and the typical formation characteristics, the cathodic protection model of pulsed current for oil well casing was established, and the effects of wave pulse current and other factors on the cathodic protection effect of square-wave pulse current were studied.
According to the structure of oil well casing, the model of pulse current cathodic protection for casing outer wall was established. The boundary conditions were used to solve the differential equations in the model. The cathodic protection potential and current density at any depth on the outer surface of casing were analyzed theoretically. Furthermore, the model was applied to typical stratum structure in China. According to the characteristics of stratum structure, the calculation methods of protective potential of casing in each stratum were obtained, respectively. The characteristics of different formations were analyzed, the laboratory model of casing was established, and the pulse current cathodic protection test was carried out.
The influence of pulse current parameters (frequency and duty cycle) on casing protection potential was analyzed through the test results, and the rule of influence was obtained. With this experimental model, the effect of anode distance on casing protection potential was further tested. The influence of basic parameters of pulse current cathodic protection for simulated oil well casing was discussed through experiments.
Through the simulation of the established test model, it is found that in such multilayer geological structure, the oil well casing can be well protected when the frequency of pulse power supply parameters is 3500 Hz and the duty cycle is 50%, and there is an optimum anode distance. This can then be of a great help to design engineers and field engineers responsible for design and installations of oil well casing in the field.
2. Establishment of Pulse Current Cathodic Protection Model for Casing
2.1. Theoretical Analysis of Cathodic Protection Potential Distribution of Casing
The schematic diagram of the external pulse current model of the casing is shown in Figure 1. The cathodic protection current on the outer surface of the casing can be expressed as follows:where is a small section on the casing, E is the casing potential of the small section, is the transition resistance, and I is the casing current of the small section. Here, the earth potential is zero.
In addition, when the current flows through the section, there will be voltage drop due to the resistance of the pipeline itself. If the average current passing through the section pipeline is I and the resistance per unit length of the pipeline is , the voltage generated by the small pipeline can be expressed by the following equation:
The general solutions of the above differential equations are as follows:where A1, B1, A2, and B2 are constants.
Boundary conditions: if x = 0, then E = E0, I = I0; if , then E = 0, I = 0.
Formula (7) shows that the applied potential on the casing path varies exponentially. Its characteristic is that the potential near the confluence point decreases rapidly, and the farther away from the confluence point, the slower the decline. The potential E depends mainly on the transition resistance RT. According to formulas (1) and (6), the protective potential at any point of casing can be calculated:
Formula (8) shows that the decline rate of the applied potential curve of the casing (the gradient of potential change) depends on the attenuation coefficient and the pipeline transition resistance RT. If the effect of pipeline spontaneous potential is considered, the protective potential E of the casing can be expressed as follows:where U0 is the pipeline spontaneous potential.
2.2. Potential Distribution of Cathodic Protection in Casing
In engineering application, the stratum structure through which oil casing crosses is complex, including soil surface layer, aquifer, various rock layers, and oil-bearing layers. Oil wells in the Shaanxi-Gansu-Ningxia Basin with wide distribution and stable lithology are mainly distributed at the junction of Shaanxi, Gansu, and Ningxia. The lower part of the stratigraphic structure is mainly purple-red mudstone, and the bottom is yellow medium-to coarse-grained feldspar sandstone with huge oblique bedding. Mudstone and argillaceous siltstone with purple-red predominance in the middle of the basin are intercalated with fine-grained feldspar sandstones with oblique bedding. The upper part is red fine-to-coarse-grained feldspar sandstone with large oblique bedding, which contains fine gravel and mud gravel. The schematic diagram of the underlying structure model is shown in Figure 2. Based on the pulse current model of Figure 1 and formulas (4) and (5), the current field model of the casing pulse current cathodic protection system including the outer wall interface of the casing is established. The relationship between casing surface protection potential and depth x at any depth along well depth is studied by the model.
According to the theoretical analysis in the previous section and Figure 2, the calculation method of casing protection potential in the topsoil layer is as follows:where is the potential of the topsoil layer, RTx1 is the resistivity of the topsoil layer, is the attenuation coefficient of the topsoil layer, Lx1 is the depth of the topsoil layer, and RX1 is the resistance at x meters from the wellhead in the soil layer.
The first rock layer casing protection potential is as follows:where is the potential of the first rock layer, RTx2 is the resistivity of the first rock layer, is the attenuation coefficient of the first rock layer, Lx2 is the depth of the first rock layer, and RX2 is the resistance at x meters from the wellhead in the first rock layer.
The second rock layer casing protection potential is as follows:where is the potential of the second rock layer, RTx3 is the resistivity of the second rock layer, is the attenuation coefficient of the second rock layer, Lx3 is the depth of the second rock layer, and RX3 is the resistance at x meters from the wellhead in the second rock layer. Using the resistivity and depth data of water, rock, clay, and oil layers, the protective potential law of casing in each layer can be obtained. The following experiments are based on the above theoretical analysis and the resistivity and depth of each layer obtained from the data.
3.1. Test Model
Based on the above model and the resistivity of the strata layer, the test model was established, as shown in Figure 3. The cathodic protection system for simulated vertical oil well casing was built in the 5 × 0.5 × 0.5 m PVC flume. The test medium was salt water, and the cathode was 40 angle steel. The cathode was placed 20 cm below the liquid level of the tank horizontally. The water entry depth of the anode was the same as that of the cathode. The power supply was a self-made pulse power supply. The cathodic protection potential was collected by the oscilloscope.
The formation resistivity of the model in Figure 2 was calculated for the protection potential of casing in a reasonable range. According to formulas (10)–(12) and the requirement of cathodic protection potential, the resistivity of each layer can be calculated separately. But the formation resistivity varies greatly in each oil production area. And it is also affected by formation structure density and porosity, pore fluid content, and liquid properties. Sun et al. measured the resistivity of multilayer stratum structure in a certain area, and the variation range of resistivity was between 0 and 30 . Zhang simulated the effect of different resistivity of different strata on cathodic protection potential . Ghimire et al. studied the change of formation resistivity around wellbore . Liu et al. studied the formation resistivity changes of oil wells in Zhidan area, northern Shaanxi . They found that the average resistivity was 5.2 and the resistivity was not higher than 21 when using acoustic time difference to measure formation resistivity. By analyzing the data in the above literatures, it can be concluded that the formation resistivity varies greatly, and it is difficult to find the law of formation resistance variation.
The oil wells in the literature  and the oil well stratigraphic structure models in the paper are all located in the Shaanxi-Gansu-Ningxia Basin. Therefore, the resistivity of each layer in the simulation test is measured between 0 and 21 based on Liu’s research. And considering the requirement of the cathodic protection potential, different resistivity values are calculated for different strata.
According to formula (10), the resistivity of the soil layer is calculated as follows:
When the boundary conditions and the resistivity of casing are known, the resistivity of soil layer can be calculated by setting the appropriate value of the protective potential of soil layer. In the same way, the resistivity of the first rock layer and the second rock layer can be calculated according to formulas (11) and (12). When the protective potentials , , and are −0.85 V, the boundary condition is 0.3 V, is 0.180 A in the model, and the resistivity of the casing is , the calculation results are shown in the table below.
According to the formation resistivity in Table 1, the formation resistance in the model of Figure 3 is calculated. Assuming that, in Figure 3, resistance RX1 is located 180 meters from the wellhead, the resistance RX2 is located 1000 meters from the wellhead, resistance RX3 is located 1500 meters from the wellhead, and resistance is located 1800 meters from the wellhead, according to Figure 2, those resistance RX1, RX2, RX3, and RX3' are the resistances of the soil layer, the first layer of rock layer, the upper part of the second layer of rock layer, and the lower part of the second layer of rock layer, respectively. According to the stratum structure of Figure 2, the resistance values of RX1, RX2, RX3, and are calculated to be 896 , 7942 , 12183 , and 15704 by using the resistivity of the water-bearing stratum 1 and the resistivity in Table 1. Therefore, the calculated resistance values were brought into the simulation test of Figure 3. Experiments were carried out in accordance with Figure 3, as shown in Figure 4.
3.2. Test Process
The test process included the preparation of test instruments and equipment, sample preparation, sample loading, determination of test time, selection of power supply parameters, etc. This experiment included the following processes.
3.2.1. Test Equipment Preparation
(1)The self-made pulse power supply was selected as the power supply.(2)The saturated copper sulfate electrode was selected as the reference electrode.(3)The graphite anode was selected as the anode.(4)The corrosion test tank was customized according to size. The white PVC plate trough of 5 × 0.5 × 0.5 m was custom-made to insulate and prevent stray current interference.(5)TDS1002 oscilloscope was used to monitor the potential, amplitude, and frequency of pulse waveform.
3.2.2. Sample Preparation
The national standard 40 angle steel was selected as the sample with a length of 4.5 m. The sample was polished and rust removed before the test. After polishing, the sample was degreased and dried with acetone.
3.2.3. Sample Loading Method
In this experiment, the shelf type was selected and the samples were placed on the self-made PVC plastic shelf. In order to place the sample correctly, a bracket was placed at the bottom of the flume. The sample was placed on the bracket. The following principles were followed in the test process.(1)During the test, the sample was properly placed and insulated from the metal in the test system, and the stray current interference was prevented.(2)The sample intruded into the environmental medium as a whole, and the depth of immersion was 20 cm. The depth of immersion of the auxiliary anode and the reference electrode was consistent with the depth of the sample.(3)During the test, the sample was easy to take and place.
3.2.4. Planned Test Time
According to the relevant research results, the annual corrosion rate of oil well casing in the stratum structure of Shaanxi-Gansu-Ningxia Basin varies from 0.1 mm/a to 1.0 mm/a. Table 2 shows the time required for corrosion experiments. According to the data in Table 2, the simulation experiment time was 72–168 hours.
3.2.5. Selection of the Orthogonal Test Method
In order to investigate the protective effect of oil well casing under impulse current, six levels of frequency (F) and duty cycle (P) and three levels of duty cycle (P) were selected. Table 3 lists the experimental factors and levels.
4. Results and Discussion
In order to study the relationship between the outer surface protection potential E and depth x at any depth along well depth, DC current cathodic protection test, pulse current cathodic protection test, and the influence of anode distance on the effect of pulse current cathodic protection were carried out according to the established model. According to the established cathodic protection model of casing pulsed current and the resistivity of different formations, the potential E at depth x can be calculated. Here, the soil layer with reference electrode 1 was set up at 180 meters, the rock layer with reference electrode 2 at 1000 meters, the rock layer with reference electrode 3 at 1500 meters, and the rock layer with reference electrode 4 at 1800 meters. The transition resistances at each stratum were calculated, respectively, and the resistors were connected to the test device for testing. According to the criteria for judging the effect of cathodic protection, the effective protection can be achieved when the polarization potential falls between −0.77 and −l.03 V.
4.1. DC Current Cathodic Protection Test
Figure 5 is a DC current cathodic protection potential diagram. After 72 hours of establishment and operation of the simulation system, data acquisition was carried out for the experiment. Data were recorded every 60 minutes, as shown in Figure 5.
From the analysis in Figure 5, the results can be obtained. The protective effect of the casing at the reference electrode 1 is better, and the protective potential is within the range of cathodic protection potential. However, the protective potentials at reference electrodes at 2, 3, and 4 sites are not within the range of cathodic protection potential and are not protected by cathodic protection. The overall protection depth of the specimen is not enough.
4.2. Pulse Current Cathodic Protection Experiment
Detailed data of square-wave stable protection potential under different duty cycles were recorded, respectively. The protective potentials at four locations were recorded in six groups, each with an interval of 1 hour. The average value was obtained when analyzing the data. Figures 6–8 are experimental results at different frequencies and duty cycles.
Figure 6 shows the effect of pulse current cathodic protection for bushing with duty cycle of 20% at different frequencies. By analyzing Figure 6, the following results can be obtained. The protective potential at reference electrode 1 exceeds the range of protective potential from 1500 Hz to 4500 Hz, with a maximum of −1.03 V, and is in an overprotective state. At very low or very high frequencies such as 500 Hz or 5500 Hz, the potential at reference electrode 1 exceeds the protective potential less. From 1500 Hz to 5500 Hz, the protective potential at reference electrode 2 is within the range of the cathodic protection potential, and the protective potentials at reference electrode 3 and reference electrode 4 are lower than the range of cathodic protection potential, without cathodic protection. In DC current cathodic protection, the protective potentials at reference electrode 1 and reference electrode 2 are in the range of −0.77 V to 1.03 V, and the protective potentials at reference electrode 3 and reference electrode 4 are lower than cathodic protection potential, which have not been effectively protected. Compared with DC current cathodic protection and pulse current cathodic protection, it is found that with the increase of casing depth, the protection potential decreases sharply, and the trend of DC current cathodic protection is more obvious. When the frequency of pulse current cathodic protection is lower than 3500 Hz or lower, the potential value of pulse current cathodic protection is lower than that of DC current cathodic protection. When the frequency of pulse current cathodic protection is between 3500 Hz and 5500 Hz or higher, the potential value of pulse current cathodic protection is higher than that of DC current cathodic protection.
Figure 7 shows the effect of pulse current cathodic protection for bushing with duty cycle of 50% at different frequencies. By analyzing Figure 7, the following results can be obtained. When the frequency is 2500 Hz or 3500 Hz, the protective potentials at reference electrode 1 and reference electrode 2 are in the range of protective potential, and they are well protected. From 1500 Hz to 5500 Hz, the protective potentials at reference electrodes 3 and 4 are still lower than cathodic protection potential, but the protective potential value is higher than that of 20% duty cycle time, which are closer to the protective potential range. In DC current cathodic protection, the protective potentials at reference electrodes 1 and 2 are in the range of −0.77 V to −1.03 V, and the protective potentials at reference electrodes 3 and 4 are lower than cathodic protection potential, which have not been effectively protected. Compared with DC current cathodic protection and pulse current cathodic protection, it is found that DC current cathodic protection has a steep descent trend at reference electrode 1 and reference electrode 2 than pulse current cathodic protection. When the frequency is 3500 Hz, the protective potential at the reference electrode 3 is within the range of cathodic protection potential, and cathodic protection is obtained. At other frequencies, the reference electrode 3 is not cathodically protected. When the frequency of pulse current cathodic protection is higher than 500 Hz, the potential values of pulse current cathodic protection at reference electrode 3 and reference electrode 4 are higher than those of DC current cathodic protection.
Figure 8 shows the effect of pulse current cathodic protection for bushing with duty cycle of 70% at different frequencies. By analyzing Figure 8, the following results can be obtained. The protective potential at reference electrode 1 is in overprotection state because of the higher protective potential. The protective potential at reference electrode 2 is in the protective potential range, and it is effectively protected. The protective potentials at reference electrode 3 and reference electrode 4 are lower than that of cathodic protection potential range, which does not reach cathodic protection. In DC current cathodic protection, the protective potentials of reference electrode 1 and reference electrode 2 are in the range of −0.77 V to −1.03 V, and the protective potentials of reference electrode 3 and reference electrode 4 are lower than cathodic protection potential, which have not been effectively protected. The experimental results of DC current cathodic protection and pulsed current cathodic protection show that the protection potential of pulsed current cathodic protection is higher than that of DC current cathodic protection at different frequencies.
By combining Figures 6–8, the difference of data was calculated. A set of differences was obtained by subtracting the cathodic protection potential data in Figure 7 from the cathodic protection potential data in Figure 6. The specific method was to subtract the cathodic protection potential measured at the same position and frequency. The calculation results are shown in Table 4.
In the same way, another set of differences was obtained by subtracting the cathodic protection potential data in Figure 8 from the cathodic protection potential data in Figure 7. The calculation results are shown in Table 5.
By combining Tables 4 and 5, it can be seen that if the duty cycle of the pulse power supply is 50%, the difference between reference electrode 2, reference electrode 3, and reference electrode 4 is positive, so the protection effect is better than that at duty cycle 20%. When the duty cycle of the pulse power supply is 70%, the difference between the reference electrode 2, reference electrode 3, and reference electrode 4 appears to be negative. Therefore, the optimal duty cycle parameter is 50%. Figure 7 shows that when the frequency is 3500 Hz, the cathodic protection potential is the largest at reference electrode 2, reference electrode 3, and reference electrode 4, so the optimal frequency is 3500 Hz. So it can be seen that when duty cycle is 50% and frequency is 3500 Hz, the effect of pulse current cathodic protection is better than that of DC current cathodic protection, and the protection effect is the best. Therefore, it can be concluded that the protection effect is better when the frequency is appropriate in the middle and the duty cycle is 50% in the test. The above simulation results provide a reference for the application of pulse current cathodic protection technology in actual oil well casing system.
4.3. Effect of Anode Distance on Pulse Current Cathodic Protection
According to the optimum parameters (F = 3500 Hz, ) obtained from the above simulation tests, the cathodic protection potential distribution of each measuring point is shown in Figure 9.
Figure 9 shows the potential distribution curve on the casing surface at different anode distance D. When the anode distance D is 1 cm, the distribution of points on the surface of the casing is between 698 and 1061 mV. The upper part of the casing (at reference electrode 1) has begun to appear overprotection phenomenon, while the lower part (at reference electrode 3) has not yet reached effective protection. When D is 5 cm, the distribution range of points on the surface of casing is 705–1067 mV. The upper part of casing (at reference electrode 1) is still protected, but it is not obvious. The lower part (at reference electrode 4) is still not effective, but it is improved when D is 1 cm. When D is 10 cm, the distribution of points on casing surface is between 735 and 1014 mV. There is no overprotection at the upper part (at reference electrode 1). The lower part (at reference electrode 4) of the casing is still not effective. When D is 20 cm, the point distribution on the casing surface is between 666 and 1059 mV. The upper part of the tube (at reference electrode 1) has begun to appear overprotection phenomenon, while the lower part (at reference electrode 4) has not yet reached effective protection.
From the analysis of the above results, we can draw a very important conclusion: There is an optimum anode distance in this experiment. However, only by adjusting the anode distance, the whole casing cannot be fully protected; that is, from the analysis of the experimental data, when the anode distance is 10 cm, the protection effect of the whole casing is better. Therefore, it is very difficult to achieve effective protection for the whole line by adjusting D value only. At this time, the parameters of pulse current should be adjusted appropriately.
When using pulse current to cathodic protection of oil well casings distributed in Shaanxi-Gansu-Ningxia Basin, the following conclusions can be drawn from the above theoretical analysis and experiments.(1)The cathode will be better protected when the parameters of pulse current cathodic protection are selected at a higher frequency and a moderate duty cycle(2)There is an optimum value of the anode distance, which makes the pulse cathodic protection system provide better protection effect
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.
This research was financially supported by the Shaanxi Science and Technology Coordination Innovation Project (Grant no. 2013KTCL-04).
- M. H. Parsa, S. R. Allahkaram, and A. H. Ghobadi, “Simulation of cathodic protection potential distributions on oil well casings,” Journal of Petroleum Science and Engineering, vol. 72, no. 3-4, pp. 215–219, 2010.
- P. Miltiadou and L. C. Wrobel, “Optimization of cathodic protection systems using boundary elements and genetic algorithms,” Corrosion, vol. 58, no. 11, pp. 912–921, 2002.
- NACE International, Application of Cathodic Protection for External Surfaces of Steel Well Casings HD-1994, NACE International, Houston, TX, USA, 1994.
- R. Brousseau, S. Qian, B. Gummow et al., “Investigating the effect of pulse cathodic protection underneath a disbanded coating,” Journal of Materials Performance, vol. 32, no. 6, pp. 21–26, 1993.
- T. M. Doniguian, American Gas Association Operating Section Proceedings, MGA, Arlington, VA, USA, 1982.
- J. Bauman, “Well casing cathodic protection utilizing pulse current technology,” in Proceedings of the NACE Corrosion 2004 Conference, New Orleans, LA, USA, March-April 2004.
- N. Bich and J. Bauman, “Pulsed current cathodicprotection of well casings,” Journal of Materials Performance, vol. 34, no. 4, pp. 17–21, 1995.
- Y. Qiu and H. Wang, “Simulation study on square wave pulse current cathodic protection of oil well casing,” Oil-Gasfield Surface Engineering, vol. 19, no. 6, pp. 35–37, 2000.
- Y. Qiu, Z. Guo, H. Lin et al., “Pulsed current cathodic protection technology,” Corrosion Science and Technology Protection, vol. 13, no. 4, pp. 226–229, 2001.
- I. A. Metwally, “Factors affecting pulsed-cathodic protection effectiveness for deep well casings,” Anti-Corrosion Methods and Materials, vol. 56, no. 4, pp. 195–205, 2009.
- J. George, Experimental Setup to Study the Effect of Pulse Width Modulated Signal on Impressed Current Cathodic System, National Association of Corrosion Engineers International, Houston, TX, USA, 2018.
- H. Sun, W. Ni, and L. Xin, “Signal detection of lateral resistivity LWD based onC4D technique with oil-based drilling fluid. Hangzhou,” Journal of Central South University (Science and Technology), vol. 49, no. 5, pp. 47–51, 2018.
- Z. Zhang, Research on Resistivity Logging Through Steel Casing, Huazhong University of Science & Technology, Wuhan, China, 2009.
- B. Ghimire, H. Daigle, and K. E. Gray, “Modeling formation resistivity changes due to invasion and deformation during initial leak-off test build-up,” Journal of Petroleum Science and Engineering, vol. 169, pp. 69–80, 2018.
- X. Liu, S. Zhang, and Q. Qian, “Logging identification method research of chang 612 low resistivity oil layer in Zhidan Fanchuan Region. Northern Shaanxi,” Journal of Yanan University (Natural Science Edition), vol. 36, no. 2, pp. 1281–1288, 2018.
Copyright © 2019 Xiang-qian Xu and Hao-bin Zhou. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.