Identification of the critical depth-of-cut through a 2D image of the cutting region resulting from taper cutting of brittle materials

Optical materials play a key role in both imaging and non-imaging optical systems. However, a large volume of the materials suffers from low fracture toughness and appears to be brittle in nature, including optical crystals, infrared optical glasses, and transparent ceramics, to mention a few [1–3]. For practical applications, ultra-precision diamond turning and diamond milling are commonly adopted to directly generate optical surfaces with complicated shapes on the optical materials. For instance, a freeform Alvarez lens in the mid-wave infrared was generated on single crystal germanium through ductile milling [4], and a discontinuous Fresnel structure was successfully obtained on germanium by means of diamond turning [5]. Moreover, taking advantage of tool-servo-assisted diamond turning, more complicated freeform surfaces [6], microstructured lens arrays [7, 8], and hierarchical microstructures [9] were also generated on optical brittle materials.

To facilitate the diamond cutting process, efforts have been devoted to studies on the underlying mechanism for material removal, mainly focusing on the brittle-ductile transition behavior [2, 7, 10, 11]. The critical depth-of-cut (DoC), at which the brittle-ductile transition occurs, directly reflects the machinability of the material under the specified cutting conditions, attracting intensive investigations from both theoretical [12] and experimental aspects [1, 2, 11]. No matter which aspect is focused upon, experimental determination of the critical DoC is commonly employed based on taper cutting featuring linearly varying DoC. After cutting, a microscopic image of the cutting region is accordingly taken for the determination of the DoC. In general, there are two sorts of imaging methods for surface characterization, namely the optical-microscopy-based two-dimensional (2D) imaging method, and the more predominant optical-surface-profiler-based three-dimensional (3D) imaging method [1, 13].

Although the 3D imaging method provides more detailed information of surface topography in the cutting region [1, 12, 14, 15], micro-cracks in the 2D image are more easily identified due to the high contrast between the cracked and optical surfaces; further, the facilities required for the 2D image is much cheaper and more widespread to access [13, 16]. With the 2D imaging method, the process starts with the measurement of the groove width at the onset point of the micro-cracks in the image and accordingly calculates the critical DoC from the simple geometrical relationship of the employed tool [13]. However, measurement of the groove width in the image is really subjective, with large uncertainty, and the influence of the spatial angular discrepancy between the imaging plane and the workpiece surface on the measured width of the groove is ignored. In addition, determination of the onset point of the micro-cracks for both the 2D and 3D imaging methods is also very subjective and without consistency, inevitably leading to low estimation accuracy with poor convergence for the critical DoC. Moreover, manual handling in both the 2D and 3D imaging methods is laborious and time-consuming, especially when processing a large volume of captured images.

To overcome the aforementioned defects in determining the critical DoC, an automatic identification method is proposed based on 2D image analysis. For determining the critical DoC, searching for the onset point of the micro-cracks in conventional methods is replaced by extracting the envelope curve of the whole cracking region, which may significantly improve the estimation uncertainty. Furthermore, by best-fitting the envelope curve of the micro-cracks with the projected theoretical profile of the cutting points with a specified DoC, the critical DoC can be effectively identified.

To identify the critical DoC for brittle materials, a schematic of the commonly adopted taper cutting is illustrated in figure 1(a). During cutting, the diamond tool moves straight with a constant velocity $v_{\rm t}$. Meanwhile, the sample is fixed with an inclination angle $\gamma$ related to the moving direction of the diamond tool, leading to cutting with linearly varying DoC. Thereby, a cutting region with both smooth and crack features is obtained, respectively suffering from ductile- and brittle-mode material removal with the practical DoC being smaller and larger than the critical DoC.

Zoom In Zoom Out Reset image size

As shown in figure 1(a), taking the horizontal plane as the reference, there are two inevitable installation angles of the cutting tool around the D₁ and D₂ directions. Since the tool edge is round, influences of rotation around the D₂ direction on features of the cutting region can be ignored. On the other hand, the rotation around the D₁ direction can be regarded as a part of the inclination angle $\gamma$.

To give a mathematical description of the cutting region, a fixed Cartesian coordinate system $O_s-X_sY_sZ_s$ is firstly defined, where the $O_s-X_sY_s$ plane and the to-be-machined surface of the sample without inclination are coplanar and set as the observation plane for modeling. The origin point O_s is coincident with the initial cutting point, and the $O_sY_s$ axis and the direction vector of the cutting velocity are in the same plane. From the projection relationship shown in figure 1(c), it is shown that the practical DoC may vary with respect to the cutter contact point (CCP) at the round tool edge, resulting in both ductile- and brittle-mode removals of the material for any one cutter location point (CLP). After cutting, there are two characteristic curves in the cutting region, namely the curve C₁ and C₂ as shown in figure 1(b). The curve C₁ is the intersection between the moving tool and the sample surface, and the curve C₂ is the envelope of the cracked region. In general, the practical DoC for an arbitrary point in the envelope may be equal to the critical DoC.

For the point $p(x_{\rm p}, y_{\rm p}, 0)$ at the curve C₁, the coordinate values can be derived from the geometric relationship shown in figure 1(c) as $x_{\rm p}=\sqrt{R_{\rm t}^2-\left(R_{\rm t}-d\right) ^2}$, and $y_{\rm p}=s(\cos{\gamma}){\hspace{0pt}}^{-1}$, where s is the moving distance of the tool, $R_{\rm t}$ is the nose radius of the diamond tool, and d is the distance from the tool tip to the free surface. Therefore, taking s as the parametric variable, the profile of the curve C₁ can be expressed by

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \left\{\begin{array}{@{}ll@{}} x_1=\pm\sqrt{R_{\rm t}^2-\left(R_{\rm t}-s\tan{\gamma}\right) ^2} \nonumber \\ y_1=s(\cos{\gamma})^{-1} \nonumber \\ z_1=0. \end{array} \right. \label{1} \nonumber \end{align} \tag{ 1 }$$

With the envelope curve C₂, it has the same shape as the curve C₁ and can be derived by translating curve C₁ along the $O_sY_s$ and $O_sZ_s$ directions with distances of $\Delta y_{\rm s}=h_{\rm c}(\tan{\gamma}){\hspace{0pt}}^{-1}$ and $\Delta z_{\rm s}=-h_{\rm c}$ as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \left\{\begin{array}{@{}ll@{}} x_2=\pm\sqrt{R_{\rm t}^2-\left(R_{\rm t}-s\tan{\gamma}\right) ^2} \nonumber \\ y_2=s(\cos{\gamma})^{-1}+h_{\rm c}(\tan{\gamma})^{-1} \nonumber \\ z_1=-h_{\rm c} \end{array} \right. \label{2} \nonumber \end{align} \tag{ 2 }$$

where $h_{\rm c}$ denotes the critical DoC.

When capturing the cutting region through optical microscopy, the observation plane $O_s-X_sY_s$ is parallel with the imaging plane of the microscopy. As for the installation of the sample in the microscopy system, there will be inevitable angular discrepancy around both the $O_sX_s$ and $O_sY_s$ axes, as illustrated in figure 2. Considering an arbitrary point $Q(x_2, y_2, z_2)$ at the curve C₂ in the $O_s-X_sY_sZ_s$ system, after rotating around the $O_sX_s$ axis with an angle of $\alpha$, the new position of the point Q in the $O_s-X_sY_sZ_s$ system can be expressed by

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \left\{\begin{array}{@{}ll@{}} x_{\rm r1}=x_2 \nonumber \\ y_{\rm r1}=\sqrt{y_2^2+z_2^2}\cos\theta_1 \nonumber \\ z_{\rm r1}=-\sqrt{y_2^2+z_2^2}\sin\theta_1 \end{array} \right. \label{3} \nonumber \end{align} \tag{ 3 }$$

with

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \theta_1=\arctan\left(-\frac{z_2}{y_2}\right)-\alpha. \label{4} \nonumber \end{align} \tag{ 4 }$$

Zoom In Zoom Out Reset image size

Similarly, after rotating around the symmetric axis of the curve C₁ with an angle of $\beta$, the point Q will reach a new position as expressed by

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \left\{\begin{array}{@{}ll@{}} x_{\rm r2}=-\sqrt{x_2^2+z_2^2}\sin\theta_2 \nonumber \\ y_{\rm r2}=y_{r1}+\left(\sqrt{x_2^2+z_2^2}\cos\theta_2-\left\vert z_2\right\vert\right) \sin\alpha \nonumber \\ z_{\rm r2}=z_{r1}-\left(\sqrt{x_2^2+z_2^2}\cos\theta_2-\left\vert z_2\right\vert\right) \cos\alpha \end{array} \right. \label{5} \nonumber \end{align} \tag{ 5 }$$

with

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \theta_2=\arctan\left(\frac{x_2}{z_2}\right)+\beta. \label{6} \nonumber \end{align} \tag{ 6 }$$

Therefore, the projection profile of the envelope curve C₂ in the imaging plane, namely the envelope of the cracked region in the image taken, can be obtained as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \left\{\begin{array}{@{}ll@{}} x_{\rm e}=-\sqrt{x_2^2+z_2^2}\sin\theta_2 \nonumber \\ y_{\rm e}=\sqrt{y_2^2+z_2^2}\cos\theta_1+\left(\sqrt{x_2^2+z_2^2}\cos\theta_2-\left\vert z_2\right\vert\right) \sin\alpha. \end{array} \right. \label{7} \nonumber \end{align} \tag{ 7 }$$

3.1. Problem statement

As discussed in section 2, there are four parameters governing the envelope curve C₂, namely $P_1=\left[\alpha, \beta, \gamma, h_{\rm c} \right]$. Practically, the margin of the cracked region can be automatically extracted through image analysis of the microscopic image taken. By comparing the theoretical curve and the practical margin, it will be very promising to get an accurate estimation of the critical DoC. In general, there are two key problems to be addressed before conducting a comparison between the theoretical and practical curves: (i) the coordinate system mismatching; and (ii) one-one correspondence for points in the theoretical and practical curves.

• As for the coordinate system mismatching, the theoretical curve in equation (7) is described in the $O_s-X_sY_sZ_s$ system, whereas the extracted margin of the cracked region is described in its local coordinate system. With an arbitrary point $P_{\rm m}(x_{\rm m}, y_{\rm m})$ at the margin, its new coordinates after planar position transition can be expressed by
  $$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \left[\begin{array}{@{}cc@{}} \hat{x}_{\rm {m}} \nonumber \\ \hat{y}_{\rm {m}} \end{array} \right]= \left[\begin{array}{@{}cc@{}} \cos\kappa&-\sin\kappa \nonumber \\ \sin\kappa&\cos\kappa \end{array} \right] \left[\begin{array}{@{}cc@{}} x_{\rm m} \nonumber \\y_{\rm m} \end{array} \right]- \left[\begin{array}{@{}cc@{}} x_{\rm o} \nonumber \\y_{\rm o} \end{array} \right] \label{8} \nonumber \end{align} \tag{ 8 }$$
  where $\left[x_{\rm o}, y_{\rm o} \right]^{\rm T}$ and $\kappa$ are the planar position translation vector and the rotational angle around the $O_sZ_s$ axis, respectively. The position transition, as presented in equation (8), is to align the local coordinate system of the captured image to be consistent with the pre-defined $O_s-X_sY_sZ_s$ system.
• The description of the theoretical curve is based on a parametric equation, as given in equations (2) and (7), for which it is difficult to derive an explicit expression directly relating the two variables $x_{\rm e}$ and $y_{\rm e}$. By adopting $S(\cdot)$ to denote the cubic spline interpolation function, the curve can then be directly described through the two variables as $x_{\rm e}=S\left(y_{\rm e}\right)$. With respect to the point at the margin $\left[\hat{x}_{\rm m}, \hat{y}_{\rm m}\right]$, the coordinates of the corresponding point at the theoretical curve C₂ yield $\left[S(\hat{y}_{\rm m}), \hat{y}_{\rm m} \right]^{\rm T}$. Thereby, the distance between the two points at the margin and the theoretical curve can be accordingly determined by $d=\left\vert S(\hat{y}_{\rm m})-\hat{x}_{\rm m}\right\vert$.

Therefore, the identification of the critical DoC lies in the determination of a set of proper parameters to minimize the deviation between positions of points at the margin of the cracked region and the corresponding points at the theoretical curve C₂. Mathematically, it takes the form of

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \min E=\sqrt{\frac{\sum\nolimits_{n=1}^Nd_n^2}{N}}=\sqrt{\frac{\sum\nolimits_{n=1}^N\left[S(\hat{y}_{\rm m}^{(n)})-\hat{x}_{\rm m}^{(n)}\right]^2}{N}} \label{9} \nonumber \end{align} \tag{ 9 }$$

where N denotes the number of the points extracted from the margin of the cracked region.

From the objective function shown in equation (9), except for the four governing the curve C₂, there are another three unknown parameters $P_2=\left[x_{\rm o}, y_{\rm o}, \kappa \right]$ for the planar position transition to be identified. Obviously, it is a multi-dimensional minimization problem with strong nonlinearity. The differential evolution (DE) algorithm which is a simple but very powerful tool to solve global optimization problems is employed in this study to minimize the objective function to identify the critical DoC $h_{\rm c}$.

3.2. Differential evolution-based minimization

The DE algorithm, which belongs to the evolutionary computation, uses differential computation as the basic operation. It mainly consists of four phases, namely initialization, mutation, crossover, and selection [17]. With the problem as stated above, there are seven parameters to be determined to minimize the objective function as shown in equation (9). The parameters can be expressed in vector form as $\boldsymbol{P}=\left[\alpha, \beta, \gamma, h_{\rm c}, x_{\rm o}, y_{\rm o}, \kappa \right]$.

At the initialization phase, the initial population is constructed through uniform randomization within the searching space. Mathematically, the population is obtained as [18]

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \boldsymbol{P}_{\rm o}^{{(K,L)}}=\boldsymbol{P}_{\rm min}^{(K,L)}+\left(\boldsymbol{P}_{\rm max}^{(K,L)}-\boldsymbol{P}_{\rm min}^{(K,L)}\right)\cdot \boldsymbol{\Omega}^{{(L,L)}} \label{10} \nonumber \end{align} \tag{ 10 }$$

where K and L are the population size and the numbers of parameters (i.e. L  =  7), $\boldsymbol{P}_{\rm min}^{{(K, L)}}$ and $\boldsymbol{P}_{\rm max}^{{(K, L)}}$ are the lower and upper boundaries of the searching space, and each boundary is a K by L matrix. $\boldsymbol{\Omega}^{{(L, L)}}$ is the uniformly distributed random matrix with a size of L by L, and each of the random numbers in the matrix ranges from 0 to 1.

After the initialization, a mutation operation based on differential computation is adopted to generate a mutant population, leading to the target vector $\boldsymbol{\hat{P}}_{\rm o}^{{(k, L)}}$ (the kth column in the mutant population $\boldsymbol{\hat{P}}_{\rm o}^{{(K, L)}}$) as [17]

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \boldsymbol{\hat{P}}_{\rm o}^{{(k,L)}}=\boldsymbol{{P}}_{\rm o}^{{(k_1,L)}}+F_k \cdot \left(\boldsymbol{{P}}_{\rm o}^{{(k_2,L)}}-\boldsymbol{{P}}_{\rm o}^{{(k_3,L)}}\right) \label{11} \nonumber \end{align} \tag{ 11 }$$

where k_i, i  =  1, 2, 3 are mutually exclusive integers ranging from 1 to K, and $F_k \in [0, 2]$ is the mutation factor for scaling the difference vector.

To increase the diversity of the population, a crossover operation is further introduced to form a trial vector ${{{\bf P}}}_{\rm tr}^{{(k, l)}}$, where the number in the kth column and lth row can be expressed by [17, 18]

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {{P}}_{\rm tr}^{{(k,l)}}= \left\lbrace \begin{array}{@{}ll@{}} {\hat{P}}_{\rm o}^{{(k,l)}},~{\rm if}~\varPsi_k\leqslant CR~{\rm or}~k=\zeta_l \nonumber \\ {{P}}_{\rm o}^{{(k,l)}},~{\rm otherwise} \end{array} \right. \label{12} \nonumber \end{align} \tag{ 12 }$$

where ${\Psi}_k\in [0, 1]$ is a uniform random number, and $CR\in[0, 1]$ is the crossover factor. $\zeta_l$ is a random integer ranging from 1 to K.

After the crossover operation, a selection process following the greedy criterion is conducted between the two populations ${{{\bf P}}}_{\rm tr}^{{(K, L)}}$ and $\boldsymbol{P}_{\rm o}^{{(K, L)}}$ to form the next generation population ${{{\bf P}}}_{\rm new}^{{(K, L)}}$ [17]. For the vector in the kth column, the selection process can be described by

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {{{\bf P}}}_{\rm new}^{{(k,L)}}=\left\lbrace \begin{array}{@{}ll@{}} {{{\bf P}}}_{\rm o}^{{(k,L)}},~{\rm if}~E\left({{{\bf P}}}_{\rm o}^{{(k,L)}} \right)\leqslant E\left({{{\bf P}}}_{\rm tr}^{{(k,L)}}\right) \nonumber \\ {{{\bf P}}}_{\rm tr}^{{(k,L)}},~{\rm otherwise} \end{array} \right. \label{13} \nonumber \end{align} \tag{ 13 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle g_{{\rm best}}=\min\left[E\left({{{\bf P}}}_{\rm o}^{{(k,L)}} \right),E\left({{{\bf P}}}_{\rm new}^{{(k,L)}} \right)\right],~\forall~k \nonumber \end{align} \tag{ 14 }$$

where $E(\cdot)$ is the objective function as defined in equation (9), and $g_{{\rm best}}$ is the minimized objective value, directly reflecting the deviation between the margin of the micro-cracks and the envelope curve C₂. Moreover, the 1 by L vector ${\bf p}_{{\rm best}}$ corresponding to the value $g_{{\rm best}}$ will be the best estimated unknown parameters in this step.

By replacing ${{{\bf P}}}_{\rm o}^{{(K, L)}}$ with ${{{\bf P}}}_{\rm new}^{{(K, L)}}$ and repeating the mutation, crossover, and selection phases, the minimization will iteratively run until the stop criteria is satisfied. After finishing the iteration, the final vector ${\bf p}_{{\rm best}}=\left[\hat{\alpha}, \hat{\beta}, \hat{\gamma}, \hat{h}_{\rm c}, \hat{x}_{\rm o}, \hat{y}_{\rm o}, \hat{\kappa} \right]$ will be the most proper parameters, and $\hat{h}_{\rm c}$ is the identified critical DoC.

Numerical simulation is conducted to verify the feasibility of the proposed identification method. A curve governed by randomly chosen parameters $\boldsymbol{P}_1= \left[0.5^{\circ}, 0.2^{\circ}, 0.09^{\circ}, 140~{\rm nm}, \right.$$\left.10~\mu{\rm m}, 5~\mu{\rm m}, 2^{\circ} \right]$ is employed to mimic the envelope curve of the cracked region. For the DE minimization, the population size and iteration step are set as 30 and 100, respectively, and the crossover factor is 0.7. To investigate the practical working performance, uniformly randomized noise with a mean value of 6 μm and amplitude of 12 μm is superimposed on the envelope curve to simulate the extracted margin of the micro-cracks.

The simulation for each case was conducted 20 times, and the envelope curves and the margins of the micro-cracks after position transition are illustrated in figures 3(a) and (b) with good fitness. As for the simulated margin without noise, a mean value of 139.90 nm with a standard deviation of 0.21 nm is achieved for the critical DoC. Compared with the pre-set true value($h_{\rm c}=140~{\rm nm}$), the deviation is only about 0.7‰. Similarly, as for the simulated margin with noise, a mean value of 139.51 nm with a standard deviation of 0.26 nm is achieved for the critical DoC, and the deviation of the mean value from the true one is about 3.5‰. It is slightly larger than the result obtained without noise but maintained as relatively small. The numerical simulation result suggests that the proposed identification method has a very high identification accuracy with sub-nanometric uncertainty.

Zoom In Zoom Out Reset image size

To experimentally investigate the feasibility of the proposed method, taper cutting of single crystal silicon was conducted on an ultra-precision machine tool (Moore nanotech 350FG, USA), adopting a natural single crystal diamond tool (Contour Fine Tooling, UK) with a nose radius of 1.549 mm and rake angle of $-20^{\circ}$. The cutting speed is 0.5 mm s⁻¹. After cutting, the features of the machined surface are captured through an optical microscope (BX60 Olympus, Japan).

First of all, the cutting velocity was oriented along the $\langle1\, 1\, 0\rangle$ direction. The obtained optical microscopic image is illustrated in figure 4(a), and the corresponding gray and binary images are generated through image analysis, as shown in figures 4(b) and (c), respectively. By using the binary image, features in the cracked region will be much clearer to facilitate extraction of the margin of the micro-cracks. The extracted margin and the identified envelope curve are given in figure 4(d), showing good agreement between each other. After conducting the calculation for 20 times, the obtained mean value of the critical DoC is about 140.22 nm with a standard deviation of about 0.11 nm.

Zoom In Zoom Out Reset image size

Moreover, taper cutting along the $\langle\bar{1}\, 1\, 0\rangle$ and $\langle\bar{1\, 1}\, 0\rangle$ directions is conducted to further verify the universality of the proposed method. To avoid repetition, only the obtained gray images of the surfaces resulting in cutting along the two directions are illustrated in figures 5(a) and (c), and the extracted margins and the corresponding identified envelope curves are shown in figures 5(b) and (d), respectively. Similarly, the envelope curves fit well with the margins of the micro-cracks. The identified critical DoCs for the two directions are 125.87 nm and 125.61 nm, and the corresponding standard deviations are about 0.42 nm and 0.52 nm. From the identified results, sub-nanometric uncertainties are obtained for all three obtained surfaces, well demonstrating the effectiveness of the proposed method for identifying the critical DoC.

Zoom In Zoom Out Reset image size

It is noteworthy that there are intermittent micro-cracks at the onset stage of the formation of the cracked region, which is especially evident for cutting using small inclination angles, as illustrated in figures 4 and 5(c). This chaotic region may contribute to the brittle-ductile transition as reported before [13]. From the identified envelope curves shown in figures 4(d), 5(b) and (d), the margins of the generated micro-cracks are well within the envelope curve for any positions in the cutting region, suggesting that even with the growth of lateral micro-cracks, the critical DoC that triggers the onset of the micro-cracks is almost unaffected. Moreover, there is also a hybrid brittle-ductile region inside and close to the envelope curve, suggesting that the brittle-ductile transition region also occurs at CCPs when the practical DoC approaches the critical DoC.

An identification method for the critical DoC is proposed based on a 2D microscopic image of the cutting region resulting from taper cutting of brittle materials. It mainly consists of three stages: (i) extraction of the margin of micro-cracks through image analysis; (ii) theoretical formulation of the projected curve of the cutting points on the imaging plane with respect to a specified DoC during the whole cutting process; (iii) minimization of the deviation between the theoretical curve and the margin of the micro-cracks through differential evolution algorithm.

To assess the proposed identification method, numerical simulation is firstly conducted by using an ideal profile as the crack margin, and the uniformly randomized noise is then superimposed on the profile to simulate the micro-cracks. For both cases, the obtained relative identification errors are less than 3.5‰ with uncertainty less than 0.3 nm. Practical cutting of single crystal silicon is then carried out along the $\langle 1\, 1\, 0\rangle$, $\langle\bar{1}\, 1\, 0\rangle$, and $\langle\bar{1\, 1\,}0\rangle$ directions. The identified critical DoCs are about 140.22 nm, 125.87 nm, and 125.61 nm, respectively. The corresponding uncertainties are about 0.11 nm, 0.42 nm, and 0.52 nm. Both numerical simulation and experimental results well demonstrate the effectiveness of the proposed identification method.

This work was jointly supported by the National Natural Science Foundation of China (51705254, 51675455), the Natural Science Foundation of Jiangsu Province (BK20170836), and the Fundamental Research Funds for the Central Universities (30917011301, 309171A8804).