Face gear generating grinding residual model based on the normal cutting depth iterative method

The generating grinding method has the characteristics of high machining accuracy and high efficiency. Therefore, it is widely used in the finishing process of the face gear tooth surface. The physical performance of the grinding process is an important factor that influences the machining accuracy of the face gear tooth surface. Therefore, to improve the manufacturing accuracy of the face gear, it is necessary to accurately simulate the grinding process of the face gear. Owing to the influence of the complex spatial geometric characteristics of the face gear tooth surface and the grinding wheel, establishing a surface residual modeling method for the face gear is a key challenge in simulating the generating grinding process. To address the aforementioned issues, this study proposes a normal cutting depth iterative method to calculate the grinding residual in the face gear generating grinding process. This method considers the complex 3D spatial characteristics of the face gear tooth surface and establishes the spatial residual model of each node of the face gear surface in the process of generating grinding. Compared with other residual algorithms based on 2D truncation or Boolean operation of face gears, this algorithm has higher computational accuracy and efficiency. Thus, it lays the foundation for accurately establishing the complex spatial microscopic surface topography in the face gear generating grinding process.


Introduction
Face gear pairs have the advantages of a large transmission ratio, high bearing capacity, and compact structure. They are the core transmission components of key research and development equipment such as helicopters, top-grade precision CNC machine tools, and precision machining robots [1]. Meeting the development requirements of high efficiency, high precision, and intellectualization for such equipment involves high requirements in terms of the manufacturing accuracy and surface quality of its profiles. The generating grinding method has the characteristics of high machining accuracy and high efficiency. Hence, it is widely used in the finishing process of face gears. To improve the generating grinding manufacturing accuracy of face gears, it is necessary to accurately simulate the grinding process [2]. And the dynamic physical parameters are the crucial factors of the grinding process. Owing to the influence of the complex spatial geometric characteristics of the face gear tooth surface and the grinding wheel, establishing a high precision and high efficiency model to simulate the grinding contact area and the manufacturing surface topography is a key challenge in simulating the generating grinding process. To investigate the surface topography of face gear grinding, previous studies conducted detailed research and exploration of the grinding mechanism, machining technology of the face gear generating grinding process, and methods for modeling the face gear grinding physics and surface topography. These studies are reviewed below.
First, with regard to simple technologies such as surface grinding and cylindrical grinding, several studies have been conducted to investigate the formation mechanism of the surface topography and the generation mechanism of the physical field during the grinding process. Malkin [3] proposed a calculation model for the dynamic effective abrasive particle number and determined the grinding surface topography and grinding force by combining this model with the grinding wheel-workpiece surface contact length algorithm. Durgumahanti [4] conducted a single abrasive particle experiment using a light wheel inlaid with a single particle, tested the model correlation coefficients of the grinding force identification plowing force, chip deformation force, and sliding force, and measured the surface topography of the grinding wheel via electron microscope experiments. Furthermore, they investigated the prediction model of the grinding force in the surface grinding process by integrating the calculation method of the grinding wheel-workpiece actual grinding domain. Thus, they laid the foundation for calculating the grinding surface topography. Hecker [5] established a conical single abrasive particle grinding force model and calculated the expected value of the cutting depth according to the relationship between the chip cross-sectional area of the conical particle and the material removal amount. Based on the Rayleigh probability density model, the cutting depth for each effective abrasive particle was determined. By integrating the actual grinding wheel-workpiece contact area model, the surface roughness and grinding force during surface grinding were calculated. Based on the Rayleigh probability distribution equation, Agarwal [6,7] studied a calculation method for the undeformed chip thickness of spherical particles. According to the relationship between the geometric model of the chip cross-sectional area of spherical particles and the material removal amount, the expected value of the cutting depth was calculated, and the actual cutting depth of each effective abrasive particle was determined. Furthermore, by integrating the actual grinding wheel-workpiece contact length model, the surface roughness and cutting force during surface grinding were calculated. Based on the theory of point contact grinding, Qiu [8] established the chip deformation model and determined the sliding area of point contact grinding. In addition, by integrating the single abrasive particle grinding force model, the grinding force model of point contact grinding was established. This laid the foundation for establishing an accurate generating grinding surface topography model. Leonesio [9] predicted the time-domain curve of the grinding force by establishing a time-domain cutting thickness model during the grinding process. They also used the modal experiment method to identify the modal mass, modal stiffness, and damping ratio of the grinding wheel spindle system. Furthermore, the vibration of the contact position during the grinding process was predicted, and the surface topography of the grinding workpiece was fitted. Ding [10][11][12][13] proposed a method for determining the surface topography of the grinding wheel with obliquely arranged particles and investigated the 3D surface topography of the grinding wheel via electron microscope experiments. The actual height of the particles was measured based on the Johnson transformation method, which proved that the height of the particles fitted based on the Rayleigh distribution is close to the actual height of the particles. On this basis, the radial bounce and motion trajectory of the grinding wheel were measured using a micrometer, and the surface roughness of the workpiece was predicted.
Based on the aforementioned studies of the grinding mechanism, researchers have conducted detailed studies of the grinding process, grinding physical field, and grinding surface topography of face gears. Guo [14] proposed a machining process for face gear grinding using a five-axis CNC machine tool and accurately modeled the surface topography of the face gear during the grinding process. Hui [15] proposed a process scheme for worm-shaped grinding wheel face gears based on a six-axis spur gear grinder and achieved efficient and high-precision machining of face gears. Wang [16] conducted a detailed study of the face gear grinding process. They proposed a new grinding method for face gears based on CBN grinding wheels and accurately calculated the CBN grinding wheel profile. Ming [17] studied the generating grinding process of face gears and established a mathematical model of the face gear surface manufacturing error under the influence of the motion error of each axis of the machine tool. Furthermore, by integrating the principal curvature of the face gear and the grinding wheel and the principal curvature direction, the equivalent diameter of the grinding wheel was obtained and the surface roughness of the face gear grinding was determined [18]. Haifeng [19] investigated the particle trajectory motion model according to the generating grinding process and determined the surface roughness of the involute gear profile using the plane projection mapping method. By integrating the single-particle grinding cross-sectional area model and based on the equivalent diameter of the grinding wheel, Wang [20] determined the surface roughness of the face gear tooth surface. Moreover, based on the plane formed by the meshing point and its normal direction, a calculation model for the grinding residual of the face gear was established for accurately modeling the surface topography of the face gear. On this basis, by integrating the grinding wheel-gear contact area in the grinding process, Wang [21] accurately modeled the grinding force and grinding temperature of the face gear, studied the boundary conditions of the 3D stress-strain differential equation, and accurately determined the participating stress on the surface of the face gear.
Thus, researchers have conducted detailed studies of the surface topography modeling method of face gear grinding. And the grinding residual is one of the important factors for the face gear surface topography. To further improve the modeling accuracy of the surface topography of the face gear, this study proposes a calculation model for generating grinding residuals of face gears based on a normal cutting depth iterative algorithm. The framework of the algorithm is shown in Fig. 1. First, mathematical models of the grinding wheel profile and the face gear tooth surface are established. According to the process scheme, the model of the grinding trajectory is established. On this basis, the mathematical model of the grinding wheel profile in the face gear coordinate system at different times is investigated, and the tooth surface of the face gear is gridded. Furthermore, the initial normal cutting residual matrix of the tooth surface is established according to the cutting allowance. Subsequently, through Taylor series expansion, the tangent plane of the grinding wheel is established. The normal cutting depth of each grid node at each moment is calculated on the basis of iterative approximation of the grinding wheel tangent plane. Moreover, based on the parameters of the meshing point of the grinding wheel and face gear, the initial conditions for the fastest convergence of the normal cutting depth iterative method are established. By taking the meshing point as the initial solution point, the grinding area is calculated and the normal cutting depth of each grid node is determined. Moreover, the normal grinding residual matrix of the whole tooth surface at the previous moment is used as the initial condition to calculate the face gear normal grinding residual matrix at the current moment. The steps mentioned before are repeated to simulate the whole grinding process. Finally, the generating grinding residual model of the face gear is obtained. Compared with the previous methods, the main advantages of this method are as follows: (1) This method combines the grinding trajectory, grinding wheel profile model, and face gear tooth surface model, and proposes a 3D spatial residual calculation method that improves the calculation accuracy of the face gear residual and lays the foundation for calculating the grinding residuals of the face gear with high precision. (2) Compared with the Boolean operation method, this method is a numerical solution method. The solution accuracy of the normal cutting depth does not depend on the grid accuracy, which improves the calculation accuracy. Moreover, this method simplifies the 3D volume grid when solving the grinding residuals into a 2D surface mesh (with 2D grid parameters), which improves the algorithm efficiency.

Mathematical model related to face gear generating grinding
The tooth surface of the involute face gear is formed by the generating motion of the spur gear shaper or the grinding wheel. Hence, to study the grinding process of the involute face gear, it is necessary to first establish the model of the grinding wheel profile. Its mathematical model r s (θ s ) is given by Eq. (1), where r bs is the radius of the base circle, θ s is the involute unfold angle, and θ os is the initial involute unfold angle.
The truncation of the grinding wheel is shown in Fig. 2. The circle O' denotes the spur gear meshing with the face gear, and its tooth shape is the same as the truncation of the grinding wheel. The bold solid line denotes the involute truncation of the grinding wheel (spur gear). In the grinding wheel coordinate system O g -X g Y g Z g , the mathematical model for the grinding wheel profile S g (θ s , β) is given by Eq. (2), where R is the distance between the center O' of the spur gear and the center of the grinding wheel, and β is the angle between the plane where the grinding wheel truncation is located and the x-y plane.
In this study, the generating grinding process of the grinding wheel feeding along the tooth thickness direction is investigated. The grinding process is shown in Fig. 2. When the grinding wheel is not in contact with the face gear, the grinding wheel and the face gear workpiece first simulate the meshing process between the spur gear and the face gear to rotate by an angle and then move to the specified position such that the distance between the center of the spur gear center O' and the tooth tip of the face gear is equal to the tooth root circle radius r f of the spur gear. The transmission ratio of the face gear and the spur gear is set to n, the rotation step of the grinding wheel is set to φ, and the meshing angle of the spur gear is set to φ ts . The grinding wheel-face gear contact trace at the meshing position is set as the first contact trace. Then, when grinding the i-th contact trace, the turning angle of the grinding wheel φ si and the turning angle of the face gear φ 2i are given by Eq. (3), where n g is the number of teeth of the spur gear and n w is the number of teeth of the face gear.
(1) s s = r bs sin s + os − s cos s + os −r bs cos s + os + s sin s + os To simplify the calculation, the grinding wheel coordinate system O g -X g Y g Z g is set to coincide with the spur gear coordinate system O'-X'Y'Z'. The feed rate of the grinding wheel along the tooth thickness direction is set to f. Based on the parameters, the profile model S Ogi (t,θ s ,β) of the grinding wheel in the face gear coordinate system O w -X w Y w Z w at different times t when grinding the i-th contact trace can be obtained as Eq. (4) shows. Figure 3 shows the coordinate system relationship between the grinding wheel and the face gear. When the tooth tip of the grinding wheel meshes with the tooth root of the face gear, the coordinate axis z g of the grinding wheel coordinate system is in the same direction as the coordinate axis x 2 of the face gear coordinate system. Further, the coordinate axis x g of the grinding wheel coordinate system is in the same direction as the coordinate axis y 2 of the face gear coordinate system. In addition, the coordinate axis y g of the grinding wheel coordinate system is in the same direction as the coordinate axis z 2 of the face gear coordinate system. Based on the principle of gear meshing, by integrating the truncation of the grinding wheel, the face gear tooth surface can be modeled. The tooth surface model S 2 (φ s ,θ s ) is expressed as Eq. (5) shows.

Normal cutting depth iterative algorithm and its initial condition solution
As the grinding wheel profile is a sculpture surface, obtaining the normal cutting depth is a nonlinear equation system problem, which is difficult to solve. In view of this problem, by combining the mathematical models of the grinding wheel profile, tooth surface, and grinding wheel motion trajectory, this study proposes a normal cutting depth iterative algorithm. The grinding residuals of the face gear surface are established by calculating the normal cutting depth of each grid node of the face gear. This algorithm is a numerical solution algorithm for obtaining the normal cutting depth by approximating the grinding wheel tangent plane. It converts the nonlinear equation system problem into a linear equation system problem, which can be solved to obtain the normal cutting depth accurately and efficiently. The principle of the algorithm is shown in Fig. 4, where P is the grid node on the tooth surface, and its normal intersects the grinding wheel profile at a point A. PA is the normal cutting depth. A 0 is the initial condition of the iteration. P t is the meshing point of the grinding wheel and the face gear at time t. A k is the iteration point of the k-th step. S 0 ' and S k ' are the tangent O g x g y g Fig. 2 The generating grinding wheel of the face gear  (1) First, the initial conditions of the algorithm must be determined. As shown in Eq. (2), the parameters of the grinding wheel profile model are θ s and β. As the meshing points of the grinding wheel and the face gear at each grinding moment as Fig. 5 shows are all on the x-y plane of the grinding wheel coordinate system, the grinding wheel-face gear contact area is near the x-y plane. For the initial condition of the iterative calculation, the parameter β is taken as 0. Further, the involute parameter θ s of the meshing point of the face gear and grinding wheel is set to be the initial condition of the iteration.
(2) The tangent plane S0' of the initial iteration point A0 of the grinding wheel is determined. The tangent plane S0' and the normal PA of the grid node P intersect at a point A0'. The parameter coordinates (θs0', β0') of the intersection point A0' of the tangent plane and the normal line are calculated. The calculated parameter coordinates (θs0', β0') are mapped to the face gear tooth surface to obtain the tooth surface point A1. Its parameter coordinates (θs1, β1) are the same as the parameter coordinates (θs0', β0') of A0'.
(3) The tangent plane S 1 ' continues to be generated. These steps are repeated until the iteration converges.
The mathematical models required for each step of the algorithm are established. The normal vector of the face gear must be calculated first. The normal direction of the tooth surface could be established based on Eq. (5). As shown in Eq. (6), 2 s is the partial derivative of the face gear tooth surface S 2 with respect to the parameter φ s ; and 2 s is the partial derivative of S 2 with respect to the parameter θ s . Further, n x2 , n y2 , and n z2 are the components of the normal direction of the face gear tooth surface in the x w , y w , and z w directions, respectively.
According to the normal direction and the grid node coordinates, the normal equation is listed as shown in Eq. (7), where u is the line parameter.
Subsequently, based on the Taylor expansion, the tangential model S k ' of the k-th iteration point A k on the grinding in the x w , y w , and z w directions of the face gear coordinate system, respectively, while n βgx , n βgy , and n βgz are the components of k in the x w , y w , and z w directions, respectively.
Then, the A k+1 parameter coordinates (θs k+1 , β k+1 ) are given by Eq. (9), where Δ sk and Δ k are the increments of the parameters θ s and β of the k-th iteration step.
By combining Eqs. (7) and (8), the incremental parameters Δθ s and Δβ in Eq. (9) can be obtained as shown in Eq. (10), where x g , y g , and z g are the coordinates of the grinding wheel profile iteration point.
The iterative process is conducted according to the method. When the solution condition conforms to Eq. (11), the iteration is considered to converge, where ε is the preset iterative convergence condition.
The normal cutting depth is given by Eq. (12), where a ij is the normal cutting depth corresponding to the j-th grid node of the i-th contact trace; a p is the normal cutting allowance; t, sk , k is the coordinate of the k-th iteration point on the grinding wheel profile; and 2 s , s is the coordinate of the grid node on the tooth surface.
n tgz n y2 x g −n tgy n z2 x g −n tgz n y2 x 2 +n tgy n z2 x 2 −n tgz n x2 y g +n tgx n z2 y g n gy n tgz n x2 −n gz n tgy n x2 −n gx n tgz n y2 +n gz n tgx n y2 +n gx n tgy n z2 −n gy n tgx n z2 + +n tgz n x2 y 2 −n tgx n z2 y 2 +n tgy n x2 z g −n tgx n y2 z g −n tgy n x2 z 2 +n tgx n y2 z 2 n gy n tgz n x2 −n gz n tgy n x2 −n gx n tgz n y2 +n gz n tgx n y2 +n gx n tgy n z2 −n gy n tgx n z2 Δ s = − n gz n y2 x g −n gy n z2 x g −n gz n y2 x 2 +n gy n z2 x 2 −n gz n x2 y g +n gx n z2 y g n gy n tgz n x2 −n gz n tgy n x2 −n gx n tgz n y2 +n gz n tgx n y2 +n gx n tgy n z2 −n gy n tgx n z2 − +n gz n x2 y 2 −n gx n z2 y 2 +n gy n x2 z g −n gx n y2 z g −n gy n x2 z 2 +n gx n y2 z 2 n gy n tgz n x2 −n gz n tgy n x2 −n gx n tgz n y2 +n gz n tgx n y2 +n gx n tgy n z2 −n gy n tgx n z2

Face gear surface grinding residual model based on the normal cutting depth iterative algorithm
Based on the model mentioned before, the normal cutting depth of each grid node can be obtained. As shown in Fig. 6, the normal cutting depth of all the grid nodes can be determined; thus, the face gear machining surface residual can be obtained. In Fig. 6, the red surface denotes the tooth surface of the face gear, the blue line denotes the normal of the face gear, and the black surface denotes the actual processing morphology. And the machining surface residual of the face gear can be modeled at each time. The  Fig. 7 The flow diagram of the face gear grinding residual prediction algorithm meshing point P t is set to be on the i-th contact trace. The modeling process is shown in Fig. 7. The specific steps are as follows.
(1) According to the feed speed f, the pose of the grinding wheel in the face gear coordinate system O w -X w Y w Z w at time t is calculated, and the model for the grinding wheel profile in the face gear coordinate system O w -X w Y w Z w is established.
(2) According to the contact trace of the face gear, the grid of the tooth surface is established for calibrating the face gear tooth surface residual model. Moreover, the tooth surface residual matrix G t is established in correspondence to the tooth surface grid node, as shown in Eq. (13), where a pij is the grinding normal residual of the j-th grid node of the i-th contact trace at time t, m is the number of grinding contact traces, and n is the number of grid nodes of each contact trace.
(3) The solution for the grinding contact area is determined. First, according to the grinding time t, the meshing point P t parameters (φ st , θ st ) on the i-th contact trace are obtained, as shown in Eq. (14). Further, as shown in Fig. 8, the grinding contact area is obtained on the tooth surface. Then the normal cutting depth of each point of the contact trace on the tooth surface grid, where the meshing point P t is located, is calculated. When the nor- p12 a p13 ... a p1j ... a p1n  a p21 a p22 a p23 ... a p2j ... a p2n  a p31 a p32 a p33 ... a p3j ... a mal cutting depth adjacent to the meshing point P t reaches the conditional expression given by Eq. (15), it means that the point cannot be on the ground at time t. The solution in this direction of the contact is ended, and it is commenced in another direction. The numbers next to the arrows in Fig. 8 represent the solution sequence. When there are grid nodes in the second direction on the contact trace that satisfy Eq. (15), the solution of the peripheral contact trace is commenced. The solution is commenced from the point closest to P t on the peripheral contact trace until the two outermost contact traces are obtained. No grid node is on the ground for the two contact traces, and the solution ends.
(4) The normal residual apij,t of each grid node in the grinding area is compared with the normal residual apij,t-1 corresponding to the grid node at Gt-1 at time t-1. When apij,t < apij,t-1, the element at this position of Gt-1 is updated until all the elements are updated to obtain the face gear tooth surface residual matrix Gt at time t. (5) The aforementioned steps are repeated continuously. The face gear tooth surface residual model at each time is simulated until the end of the simulation. Thus, the final tooth surface residual model is obtained.

Simulation result validation
To verify the accuracy of the algorithm, the involute face gear grinding process is simulated using the face gear generating grinding residual algorithm based on the normal cutting depth iterative method. In addition, a three-coordinate measuring instrument is used to measure the grinding residuals of the face gear processed under the same process parameters. The parameters of the face gear are summarized in Table 1. First, the relative positions of the grinding wheel and the tooth surface of the face gear are determined. The relative positions of the grinding wheel and the tooth surface of the face gear at different times are shown in Fig. 9.
According to the algorithm, the tooth grinding residuals in the grinding process at different times are determined. Figure 10 shows the iterative convergence process for obtaining the normal cutting residual of a certain grid node. As can be seen, the algorithm requires only three iteration steps to converge when obtaining the normal cutting residual, because the initial condition of the iteration is close to the numerical solution. The grinding residuals of the face gear grinding process are shown in Fig. 11, where the red circle denotes the position where the grinding wheel is grinding. The colored region in Fig. 11 represents the normal cutting residual. During the grinding process, the grinding wheel is in point contact with the tooth surface. Thus, the contact area between the grinding wheel and the rough surface of the face gear is a spot with a large normal cutting depth in the middle and a small normal cutting depth at the periphery. On the face gear surface, the normal cutting residual near the meshing point is small and the normal cutting residual far from the meshing point is large.
The tooth grinding residuals of the face gears processed with three generation angles of 0.914°, 1.14°, and 1.31° are simulated. As shown in Fig. 12, when the generation angle is 1.14° or 1.31°, there is a certain yellow area between the contact traces of the tooth tip area, indicating that the machining residual is large. However, when the generation angle is 0.914°, there are blue areas between the contact traces of the tooth tip area, indicating that the machining residual is smaller and the grinding track is denser. As can be seen from the results, the larger the generation angle spacing, the larger is the machining residual, which is consistent  . 9 The relative position and relative orientation of the grinding wheel and the face gear with the real machining situation. Moreover, the machining residuals of the tooth surface ground under an unfold angle of 0.914° or 1.14°are less than 2 μm, while the machining residuals of the tooth surface ground under an unfold angle of 1.31°are less than 4 μm. According to the four aforementioned sets of process parameters, the Makino machine tool is used to grind the face gear, and the German Wenzel three-coordinate measuring instrument is employed to measure the face gear. The positioning accuracy of the three-coordinate measuring machine is 0.0001 mm. The machining measurement processes are shown in Fig. 13. Figure 13(a) shows the generating grinding process of the face gear, and Fig. 13(b) shows the three-coordinate measurement process of the face gear tooth surface. In Fig. 13(a), the rotation angle of the axis B is equal to the nφ s and the rotation angle of the axis C is equal to the φ s . The grinding wheel is inlaid on the tool holder with threads. The face gear is fixed on the tooling and then clamped to the worktable of the Makino machine tool. Taking the B-axis turning angle (90°) of the machine tool as the initial position, the face is ground. The face gear tooth surface is divided into nine areas, i.e., three areas for the position of the tooth surface near the tooth tip, three areas for the middle position of the tooth surface, and three areas for the position of the tooth surface near the tooth root. The tooth surface machining error of each area is measured repeatedly using a three-coordinate measuring instrument. The number of each measurement area is shown in Fig. 14.
The simulation results of the top, middle, and root positions of the face gear tooth surface are compared with the measured results, and the results are shown in Table 2. As can be seen, the machining residuals of the tooth tip are larger, whereas the machining residuals of the middle and root parts are smaller, which is consistent with the trend of the simulation results. Furthermore, it can be seen that most of the errors between the simulation results and the actual measured results of the tooth surface in terms of the generation angle ground face gear are within 20%. This verifies that     the algorithm is basically accurate, although there is a certain error between some measured results and the predicted results. The main reasons influencing the measurement errors are grinding wheel topography [20] and positioning accuracy of the three-coordinate measuring instrument. However, these factors are not considered when predicting the surface residual in this study.

Conclusion
To solve the nonlinear solution problem of surface residual prediction in the face gear grinding process, this study proposed a face gear tooth surface residual model based on a normal cutting depth iterative algorithm. The proposed model can accurately determine the tooth surface residuals of face gear generating grinding. Its advantages are as follows.
1) The algorithm combines the real grinding wheel grinding trajectory, the grinding wheel profile, and the face gear tooth surface model, and it employs a 3D spatial residual calculation method that improves the calculation accuracy of the face gear residual.
2) Compared with the Boolean operation method, this algorithm is a numerical solution method. The solution accuracy of the normal cutting depth does not depend on the grid accuracy, which improves the calculation accuracy. Moreover, this method simplifies the 3D grid when determining the grinding residuals into a 2D grid, which improves the algorithm efficiency.
As the algorithm does not consider the impact of the surface topography of the grinding wheel when determining the machining residual of the face gear, there is a certain error in the solution. In the future, by combining this algorithm with the mathematical model of the abrasive particle motion trajectory cluster based on the measured abrasive particle position, the real mapping relationship between the grinding wheel and the face gear will be studied to improve the prediction accuracy of the surface topography of the face gear. Data availability All data generated during this study are included in this published article.
Code availability Some parts of the code about the mathematical model during the current study are available from the corresponding author on reasonable request.

Declarations
Ethics approval Not applicable.
Consent to participate Sijie Cai, Zhiqin Cai, Bin Yao, Zhihuang Shen, Jianchun Liu, Haipeng Huang, Bingjing Lin, Jianchun Lin, and Haibin Huang are willing to be listed as the author of this essay. We further confirm that all authors have checked the manuscript.

Consent to publication
This manuscript is our original work and has not been published nor has it been submitted simultaneously elsewhere, in whole or in part. We further confirm that all authors have checked the manuscript and have agreed to the publication.

Competing interests
The authors declare no competing interests.