To predict fragment separation during rock cutting, previous studies on rock cutting interactions using simulation approaches, experimental tests, and theoretical methods were considered in detail. This study used the numerical code LS-DYNA (3D) to numerically simulate fragment separation. In the simulations, a damage material model and erosion criteria were used for the base rock, and the conical pick was designated a rigid material. The conical pick moved at varying linear speeds to cut the fixed base rock. For a given linear speed of the conical pick, numerical studies were performed for various cutting depths and mechanical properties of rock. The numerical simulation results demonstrated that the cutting forces and sizes of the separated fragments increased significantly with increasing cutting depth, compressive strength, and elastic modulus of the base rock. A strong linear relationship was observed between the mean peak cutting forces obtained from the numerical, theoretical, and experimental studies with correlation coefficients of 0.698, 0.8111, 0.868, and 0.768. The simulation results also showed an exponential relationship between the specific energy and cutting depth and a linear relationship between the specific energy and compressive strength. Overall, LS-DYNA (3D) is effective and reliable for predicting the cutting performance of a conical pick.

1. Introduction

Rock cutting is frequently encountered in some industries, for example, coal mining, tunnel excavation, and oil exploitation, and is the major function of roadheaders and drilling machines. For a given rock formation, the variation in rock morphology and cutting forces in the progress of rock cutting is very important for designing cutting tools [1].

Many scholars have conducted experimental and theoretical research on mechanical performance and rock behavior. Experimental procedures were performed with a linear cutting machine by Biligin et al. [1] and Kel et al. [2]. The results of these studies were used to verify the numerical simulation results from Jiang et al. [3] and Huang et al. [4]. Using a polycrystalline diamond compact cutting test bench, Munoz et al. [5] studied the relationships among cutter inclinations, cutting forces, and rock properties. A series of linear cutting tests for observing the progress of chip formation was performed by Che et al. [6]. The classic and most authoritative theory was proposed by Evans [79] and has been accepted by scholars and used widely. The proposed improvements on this theory consider the friction angle and the compressive strength of the rock [10, 11]. Moreover, considerable numerical simulations of rock cutting have also been performed to research chip separation and have obtained scientific results [3, 4, 1239].

Generally, experimental research is the most reliable and effective method but is too costly in terms of both time and money [17]. Theoretical research is the most systematic method but is not intuitive. In contrast, the numerical simulation method can obtain accurate and available results when the conditions are specified. Furthermore, this method is faster and more detailed than the experimental method. As summarized from the references [3, 4, 1239], both two-dimensional (2D) and three-dimensional (3D) software have their own characteristics. The 2D analysis code is an effective method for simulating crack formation and fragment separation and is frequently used by scholars, while the 3D code is good for obtaining the cutting force. Therefore, few scholars have successfully simulated crack propagation and fragment separation using the 3D analysis code.

In this paper, the explicit finite element method (FEM) code LS-DYNA (3D), a computational modeling method that is good at simulating impact, blast, and penetration, was employed to model the interaction between a conical pick and rock. A damage material and erosion criteria were used in the code to dominate the rock failure and to delete the elements. The influences of the speed of the conical pick, the depth of the cut, and the mechanical properties of the base rock on the crack extent, fragment formation, and cutting force were examined. The cutting forces were verified by experimental and theoretical studies, and the relationships among specific energy, compressive strength, and cutting depth were analyzed.

2. Theoretical Considerations and Previous Simulation Studies

2.1. Theory Studies on Cutting Force

The point attack pick theory was first proposed by Evans [7]. Evans proposed that the cutting process of the point attack pick is 3D, unlike the chisel pick with a 2D attack. He demonstrated that the dominant properties of the cutting force are the tensile and compressive strength of rock, as formulated inwhere is the mean peak cutting force; , the subscript for FC, is the name of researcher; is the tensile strength of the rock; is the compressive strength of the rock; and is the tip angle of the conical pick, as shown in Figure 1.

Evans’ theory provided a reference for later research. Roxborough and Liu [10] proposed a modification theory based on Evans’ linear cutting theory, as given in (2). They demonstrated that tensile strength was the dominant rock property that is responsible for rock failure. The friction angle between the rock and cutting pick was taken into account. With the same parameters, the cutting force values calculated by (2) are closer to the experimental values than those calculated by (1). whereis the friction angle between the pick and rock and the other notations are the same as those in (1). In this equation, is limited between 16° and 30°.

Goktan [11] also completed a modification on (1), as formulated in (3). According to (1), when decreases to zero, the cutting force is greater than zero. Goktan examined and overcame the above deficiency and showed that the friction angle is one of the most important factors in the peak cutting force. Goktan considered that the main failure form of rock is tensile damage, so the compressive strength was removed in his theory.where is set as 8.5° and the other parameters are the same as those in (1) and (2).

One of the most accepted methods for predicting the cutting rate uses specific energy, which is the energy expenditure for cutting per rock volume [40]. The specific energy calculation method in this paper is given in the following equation:where is specific energy in kWh/m3, is the mean cutting force in , is cutting distance in , and is the number of the elements that separate from the base rock.

2.2. Previous Numerical Simulations of Rock Cutting

In general, three modeling methods are frequently used to investigate the tool-rock interaction: FEM, discrete element method (DEM), and finite difference method (FDM) [3, 4, 1239]. Various numerical methods and their simulation results are summarized in Table 1.

FEM is a very practical numerical procedure for simulating engineering physics problems. Both 2D and 3D FEM are widely used to solve 2D and 3D dimensional problems. In the field of rock mechanics, various numerical simulation codes, such as LS-DYNA [3, 4, 12, 13, 1518], ABAQUS [14], and rock failure process analysis (RFPA) [1929], have been used for modeling. As shown in Table 1, LS-DYNA was widely used by scholars to simulate rock-tool interaction. Jiang et al. [3], Huang et al. [4], and Li et al. [12, 13] conducted research using LS-DYNA (3D). In all of these papers, the cutting forces of the conical pick were calculated, but no scholar used this method to obtain the fragments. In conclusion, LS-DYNA (3D) is difficult to use for simulating crack propagation, as shown in Figure 2(a). Jaime et al. [15, 16] and Menezes et al. [17, 18] employed LS-DYNA (2D) to simulate the rock cutting process. The simulation results indicated that LS-DYNA (2D) is an effective approach for numerical crack generation and fragment separation, as shown in Figure 2(b). RFPA, which was developed by the Northeastern University of China, is a professional rock failure analysis tool and has been widely used to analyze 2D rock failure problems in recent years. As a scientific research team, Tang and his partners have performed many studies on rock fragmentation mechanism, simulating the progress of rock cutting using RFPA (2D) [1929].

DEM is a numerical method for simulating a discontinuous medium, and it is widely used in slope stability, tunnel excavation, and rock dynamic behavior simulations. The representative codes of DEM are EDEM and PFC, and the method includes both 2D and 3D modes. Dai et al. [30] and Su and Akcin [31] employed DEM 3D to simulate the rock cutting process. These two studies are similar with regard to the rock cutting process, such that the particles were separated out; however, every particle was independent due to the software characteristics, as shown in Figure 2(c). Lei et al. [32], Rojek [33], Huang et al. [34], and so on [35] conducted numerical experiments using DEM (2D) and achieved satisfactory results. Crack propagation is well implemented by PFC (2D), but no chip separation occurred, as shown in Figure 2(d).

FLAC is a fast Lagrangian analysis program based on FDM and is an international geotechnical engineering analysis software. Stavropoulou [36], Innaurato et al. [37], and Fang and Harrison [38, 39] employed FLAC (2D) to investigate the influence factors of chip formation.

3. Numerical Modeling Details

All the simulations in this study were carried out using the code LS-DYNA (3D). In the simulation, a conical pick and a cuboid rock were modeled, as shown in Figure 3. The conical pick moved at a specified speed, and the base rock was stationary. The rock that was cut was modeled as a cuboid block with a length of 100 mm, in Figure 3, a height of 40 mm, in Figure 3, and a width of 80 mm. The 3D 8-node hexahedron elements were used to account for both the conical pick and the rock. There were 320,000 elements and 335,421 nodes for the rock, and there were 1512 elements and 1737 nodes for the conical pick.

The conical pick had an impact angle of 57°, and the tip angle of the cutter was 80°. In addition, the friction angle between the conical pick and the base rock was set as 8.5°. To improve the calculation efficiency, the upper part of the cutter was removed, and the cutter was modeled as a rigid body. Nonreflecting boundaries were implemented for all of the rock surfaces except the front. The constraint conditions of the surface nodes for the model were detailed as follows. The nodes distributed on each surface were treated as a group. The node group displacement and rotation conditions were the following: (a) the bottom node group of the rock was constrained in all degrees of freedom, including in the -, -, and -directions and -, -, and -axis rotations, (b) both the right surface and the left surface node group of the rock were constrained in the -direction, (c) the back surface of the rock was constrained in the -direction, and (d) the conical pick was constrained in the - and -directions. The conical pick moved at various straight-line speeds of 1, 2, 3, 5, and 10 m/s in the -direction. For the given straight-line speed, the simulations were carried out for cutting depths of 3, 6, 9, and 12 mm. For the given cutting speed and cutting depth, the simulations were performed for four types of rock with different mechanical properties, as shown in Table 2 [1]. These mechanical properties were adopted to match the experimental results. With the various combinations of the previous parameters, 80 groups of simulations were performed. The conical pick and the base rock were assigned different materials. The conical pick was assumed to be a rigid material, and the rock used a damage material model [41, 42]. The damage constitutive law and erosion criteria were used to enable the conical pick erosion into the rock, to extend cracks, and to form fragments. The ERODING_SURFACE_TO_SURFACE and AUTOMATIC_GENERAL contact types were applied between the cutter and the base rock. Considering a rock with high hardness, the elements could be distorted in the cutting process, and then the hourglass energy would occupy a large proportion of the total energy. The situation greatly influenced the computational accuracy and the cutting force results. To overcome the negative volume effect, the fully integrated calculation method was adopted. The calculation time is automatically specified by the code based on the time step. Mass scaling measures were not employed to ensure that the calculation time did not influence the real time or calculation results.

4. Numerical Results and Discussion

The numerical simulation results of the chipping progress between the conical pick and the rock for a cutting depth of 9 mm, a speed 3 m/s, and a base rock of sandstone-1 are shown in Figure 4. During the cutting process, the conical pick exerted force on the rock, an initial crack formed and propagated, and then the fragments separated from the base rock. As shown in Figure 4(a), the cracks generated on the rock around the cutter due to some elements failing under tensile stress and shear stress. Most of the elements in front of the cutter were reserved and did not fail primarily because of the elastic deformation and the compressive stress. Cracks were randomly distributed because the rock is modeled as a homogeneous material; therefore, there was no preferential direction of crack propagation. As shown in Figure 4(b), the cracks propagated outward by tensile stress along with increasing displacement of the conical pick in the -direction [26]. Around the cutting pick, some rock elements failed due to shear stress and compressive stress. As shown in Figure 4(c), more cracks appeared without regularity, and some fracture regions formed. Cracks propagated on the surface and inside of the rock, which is expected to form larger cracks and chips before the fracture zone forms. Then, the cracks connect to each other on the surface and inside of the rock, and a small quantity of fragments separate from the base rock. The size of the fracture region around the pick increased due to the high-pressure region breaking away from the high confining pressure constraints. As shown in Figure 4(d), with increasing cutting distance, the fragments with random amounts and morphology were separated from the base rock. Meanwhile, a considerable number of elements between the fragments and the base rock failed due to tensile stress. Figure 4(d) also shows that fracture regions appeared at the positions where the cracks connected rather than concentrating around the cutting pick. As shown in Figure 4(e), the shape of the cutting groove was irregular because of the fragment formation. The roughness of the groove was closely related to the fragment size. Specifically, the size of the fragment increased with increasing roughness.

The variation in cutting force with distance in this numerical simulation is shown in Figure 4(f). There are three peak forces in the cutting process, and each peak force represents one cutting cycle of stages I, II, and III. Stage I is the ascending stage of the cutting force, corresponding to Figures 4(a), 4(b), and 4(c), mainly consisting of crack propagation. Stage II is the descending stage of the cutting force, matching Figure 4(d), mainly consisting of fragment separation. The sizes of the fragments have an upward trend corresponding to the increasing peak cutting force. In stage III, the cutter withstood a cutting force much smaller than the peak force due to the fragment formation. Specifically, the rock in front of the conical pick was removed in the form of a massive fragment. Hence, only a few rocks are left in the cutter speed direction to resist the conical pick.

To study the base rock failure mode, two test points at different positions were selected to analyze the variation in the stress, as shown in Figures 4(a) and 4(b). Positive values indicate compressive stress, while negative values indicate tensile stress. The point did not fail when the shear stress and compressive stress were at a maximum, as shown in Figure 5, and test point 1 bore a maximum compressive stress of 31.56 MPa and a maximum shear stress of 54.42 MPa at 0.01299 s. The test point failure when the tensile stress reached its maximum of 11.37 MPa indicated that the test point failure mode is tensile failure. This result is consistent with the assumption in the model. The tensile stress did not reach 11.6 MPa at the point failure due to the influence of the sampling frequency. Test point 2 also experienced tensile failure at a maximum tensile stress of 11.51 MPa, as shown in Figure 5(b). The fluctuation in the stress is due to the propagation of the stress wave. Therefore, the assumption made for the base rock failure mode was valid, reliable, and consistent with the assumption in the theoretical model [7, 10, 11].

The cutting results for the simulation that were carried out at a cutting depth of 9 mm and cutter speeds of 1 m/s and 2 m/s with sandstone-1 are shown in Figure 6. The cutting results for these parameters were similar to the cutting results for a cutting speed of 3 m/s, as shown in Figure 4(e). The damage area was largest at a cutting speed of 1 m/s and was smallest at a cutting speed of 2 m/s, with that for a cutting speed of 3 m/s in between. The cutting speed has little influence on the peak cutting force. However, the peak cutting force has a strong consistency with the fragment size. Therefore, chip formation did not significantly vary with changing cutting speed.

The cutting process and the variation in the cutting forces were simulated at a cutting depth of 3 mm and a cutting speed of 3 m/s with sandstone-1, as shown in Figures 7(a) and 7(b). However, cracks and fragments are difficult to obtain at a cutting depth of 3 mm. Accordingly, the cutting depth can be considered the most important factor that influences fragment formation by comparing the results at cutting depths of 3 mm to those at 9 mm, as shown in Figures 4 and 6.

To obtain the cutting depths at which the rock transitions from ductile to brittle regimes, simulations were carried out at cutting depths of 4 mm and 5 mm and at a cutting speed of 3 m/s for the four types of rock material. The cutting results at a cutting depth of 4 mm of sandstone-1 as shown in Figure 7(c). Rock fragments began to appear, though the quantity was very small. Furthermore, at the cutting depth of 5 mm, as shown in Figure 7(d), cracks and fragments were easy to obtain. Therefore, the cutting depth of 4 mm is the transition depth from the ductile to brittle regime for sandstone-1. Other simulation results showed that the cutting depth of 5 mm was the transition depth from the ductile to brittle regime for sandstone-2, 6 mm was the transition depth for sandstone-3, and 4 mm was the transition depth for limestone. Thus, the transition cutting depth is closely related to the compressive strength and elastic modulus of the base rock. Specifically, the transition cutting depth decreased with increasing compressive strength and elastic modulus. Thus, a greater cutting depth is associated with a larger crushing area and tends to generate larger fragments. In addition, the size and number of fragments determine the shape and size of the cutting groove. Thus, the cutting groove can be used to analyze fragment separation.

The variation in the forces is influenced by the cutting depth, as shown by comparing Figures 7(b) and 4(f). The mean peak cutting force increases with increasing cutting depth, which is in agreement with theoretical studies [7, 8, 11] and an experimental study [5]. The frequency of the cutting force increases with increasing cutting speed because the formation of finer fragments leads to high frequency and low peak cutting force.

Analyzing the shape of the cutting groove shows that compared with the other three rock materials shown in Figures 8(a), 8(e), and 8(g), sandstone-2, shown in Figure 8(c), obtained the fewest fragments. In contrast, analyzing the cutting groove shows that limestone obtained the most fragments. The result tends to be related to the elastic modulus of the rock: a higher elastic modulus indicates a stiffer rock and a higher crack spacing due to fracture energy release [43, 44]. Specifically, a rock with a higher elastic modulus has a tendency to separate into larger fragments. Although the elastic modulus of sandstone-3 is greater than that of sandstone-1, the cutting results showed that sandstone-1 obtained more fragments than sandstone-3 because the fragment size increased with increasing compressive strength of the rock [45].

In addition, the peak cutting force decreases with decreasing elastic modulus of the rock, as shown by comparing the results in Figure 8(b) to those in Figures 8(d), 8(f), and 8(h). The cutting force is more stable for a lower elastic modulus rock due to the formation of finer fragments. With increasing elastic modulus, the cutting force frequency decreased, but the amplitude increased. Hence, the elastic modulus of the rock greatly influences fragmentation behavior. In addition, the compressive strength has the same influence on the cutting performance.

The simulation results of the cutting forces at different cutting depths, speeds, and rock properties are given in Table 3. The relationships between the cutting force and cutting speed at different rock properties for sandstone-1, sandstone-2, sandstone-3, and limestone are shown in Figure 9. The mean peak cutting force was defined as follows. First, mean cutting force was calculated, and all the peak cutting forces were determined. Second, the values of the peak cutting force and mean cutting force were compared. Finally, the peak cutting forces that were greater than the mean cutting force were accumulated and averaged. Thus, the mean peak cutting force was calculated.

The mean peak cutting forces vary nonsignificantly with all the simulated speeds for a cutting depth of 3 mm. However, the mean peak forces did not vary significantly with cutting speed up to 3 m/s at cutting depths of 6 mm, 9 mm, and 12 mm. The influence of the conical pick speed on the mean peak cutting force was significant only at higher cutting speeds, particularly 10 m/s, and at cutting depths of 6 mm, 9 mm, and 12 mm. Specifically, the influence of conical pick speed on the mean peak cutting force is significant only at higher cutting depths and higher cutting speeds. This outcome can be attributed to the strain rate, which is the strain variation over time. The dynamic strength of the rock had a significant dependence on the strain rate. The change in the rock strength is not obvious at a lower strain rate, but when the strain rate increased to a certain value, the strength of the rock increased significantly [46, 47]. In addition, the strain rate increased significantly with increasing loading rate [48, 49]. Thus, the strain rate of the rock increased at higher cutting speeds, leading to an increase in the rock strength. Therefore, lower cutting speeds had almost no effect on the mean peak cutting force; however, when the speed value increased to a certain value, the cutting force increased with increasing cutting speed.

The relationships between the mean peak cutting force and cutting depth for the different rock types of sandstone-1, sandstone-2, sandstone-3, and limestone are shown in Figure 10. The mean peak cutting force increases almost linearly with increasing cutting depth for a given cutting speed. With increasing cutting depth, the resistance to the conical pick and cutting force both increase.

On the basis of the rock properties presented in Table 2 and the parameters given in the modeling section, the mean peak cutting forces on the conical pick were calculated using (1), (2), and (3). The results of the theoretical and experimental values at depths of 3 mm, 6 mm, and 9 mm are summarized in Table 4 [1]. The theoretical values at a depth of 12 mm are not included in this table because reference [1] did not perform an experimental study at this depth. A cutting speed of 0.127 m/s was applied to the experiment, which was less than the cutting speed used in the simulation in the paper. Because a lower speed had almost no influence on the cutting force, the experimental values were available.

The numerical, theoretical, and experimental values were significantly different from the results shown in Tables 3 and 4. However, as presented in Figure 11, a significant correlation was found between the numerical values at cutting speeds of 1 m/s and the other study values since a cutting speed of 1 m/s is closer to real excavation values.

A linear correlation coefficient of 0.698 was obtained from Figure 11(a). The regular pattern of the cutting force is consistent between the numerical study and Evans’ theory. As demonstrated in Figures 11(b) and 11(c), the numerical study and the theoretical research have a strong coincidence at linear correlation coefficients of 0.811 and 0.868, respectively. Moreover, a linear correlation coefficient of 0.768 is obtained from the relationship between the numerical modeling study and experimental research, as presented in Figure 11(d).

The relationship between the numerical and theoretical results and the relationship between the numerical and experimental results was verified by linear regression analysis. Accordingly, analysis of variance was carried out with a confidence level of 0.95. The cutting forces obtained from LS-DYNA (3D) were taken as the independent variable, while the cutting forces obtained from the theoretical and experimental method were taken as the dependent variable. Since the P values obtained from the regression analysis were less than 0.05, as shown in Table 5, the relationships are valid [50].

The linear correlation coefficient between the numerical study and Evans’ theory is lower than the coefficient from the theories of the other two researchers. The result occurs because the mechanical properties of compressive strength and tensile strength mainly lead to rock failure in Evans’ cutting theory; in contrast, the friction angle was taken into account in Roxborough’s and Goktan’s cutting theories, and only tensile strength determined rock failure. Moreover, the highest linear correlation coefficient was obtained between the numerical study and Goktan’s theory, in agreement with the studies performed by Biligin et al. [1] and Su and Akcin [31].

The slopes of the fitting equation between the numerical and theoretical results are 0.513, 0.754, and 0.938. The regression slopes suggest that the cutting forces obtained with the numerical simulation are always greater than the values calculated by theoretical models, which is consistent with experimental results. The intercepts of 0.907, 1.452, and 1.88 were acceptable when compared with the numerical results. The numerical results underestimated the experimental results by a factor of 3.5 due to the instability of the experimental cutting forces. The cutting forces were significantly influenced by the joints, bedding planes, discontinuities, and hardness inclusion of rock specimens under the experimental conditions, but these factors were not considered in the numerical simulation. This assumption caused most of the difference between the cutting force obtained from numerical results and that obtained from experimental results [3, 31]. Furthermore, since the cutting forces are also significantly influenced by the round pick tip [51] and the conical pick having approximately the same wear in the experimental conditions, the experimental cutting force was larger. However, the wear was simplified and idealized in the numerical study. Hence, the significantly different performance conditions result in different cutting forces between the model and experimental study. The rock-tool interaction was also influenced by a range of physical, mechanical, and technological parameters. Therefore, a certain dispersion of results is always expected, and trends due to the main variables were recognized at both the experimental and the numerical stage. To obtain results that are closer to the experimental studies, the joints, bedding planes, and hardness inclusion of the base rock should be considered in the numerical simulation. In addition, the cutting wear is also an important simulation aspect that needs further research.

The estimation of the specific energy is important for predicting the cutting efficiency, as explained in (4). The parameters in the equation are shown in Table 6. The relationship between the uniaxial compressive strength and specific energy and the relationship between the cutting depth and specific energy are shown in Figure 12. The specific energy decreased with increasing cutting depth and increased with increasing rock uniaxial compressive strength. A linear relationship is shown between uniaxial compressive strength and specific energy.

The minimum value of the correlation coefficient is 0.917, indicating strong agreement between the uniaxial compressive strength and specific energy. An exponential relationship exists between the uniaxial compressive strength and specific energy, consistent with the experimental results [52]. The minimum value of the correlation coefficient is 0.943, which indicates that the fitting results are ideal. The regression analysis results that were obtained at a confidence level of 0.95 are shown in Table 7. All of the values are less than 0.05, indicating that the relationships are valid and reliable.

5. Conclusion

In this research, the explicit finite element code LS-DYNA (3D) was used to model the interaction between a conical pick and base rock. In the simulation, damaged material and erosion criteria were combined to investigate the behavior during rock cutting. The crack propagation and fragment formation were studied with varying cutting depth, speed of the conical pick, and mechanical properties of the rock. The main conclusions are as follows:(1)The morphology and number of fragments obtained in the cutting process varied significantly with varying cutting depth and mechanical properties of the base rock. However, the morphology of the rock varied nonsignificantly at a lower cutting depth, since only a few fragments were generated at the lower cutting depth.(2)The cutting forces during the cutting process also varied significantly with the cutting depth and rock properties. With decreasing cutting depth, compressive strength, and elastic modulus of the base rock, the frequency of the cutting force increased, but the peak cutting force decreased.(3)The mean peak cutting forces were verified by theoretical and experimental results. A strong linear correlation coefficient was achieved between the numerical and theoretical results.(4)The linear correlation coefficient obtained between the numerical and experimental studies was lower than that obtained between the numerical and theoretical studies, since there were many uncertain factors in the experimental process that influenced the cutting force.(5)The specific energy increased linearly with increasing uniaxial compressive strength and decreased exponentially with increasing cutting depth.(6)This paper successfully proves that the explicit FEM code LS-DYNA (3D) can be used to simulate fragment separation.

Conflicts of Interest

The authors declare that there are no conflicts of interest regarding the publication of this paper.


This work was supported by projects of the National Natural Science Foundation of China (Grant no. 51375282), the Natural Science Foundation of Shandong Province (Grant no. ZR2014EEM021), the science and technology development program of Shandong province (Grant no. 2014GGX103043), and the Qingdao postdoctoral researcher applied research project (Grant no. 2016120).