Cutterhead approximation machining method of line contact spiral bevel gear pairs based on controlling topological deviations

To improve the meshing performance and increase the bearing capacity and service life of spiral gear pairs, a cutterhead approximation machining method based on controlling geometric topological deviations was proposed to solve the problem where line contact spiral bevel gears with tapered teeth depth cannot be machined by cutterheads. First, a mathematical model of the line contact conjugate flank was established, and the geometric topological deviations model of the comparison between the machining tooth flank and the theoretical tooth flank was obtained. Second, the deviation of each tooth flank point was calculated by the topological differential calculus analysis. Finally, with the machining tooth flank approaching the theoretical tooth flank as the modification objective, the additional cutting motions and machining compensation parameters of cutterheads were obtained to control the machining tooth flank deviations and reduce them to the allowable deviations of the theoretical tooth flank, and the cutterhead approximation machining of line contact spiral bevel gear pairs with tapered teeth depth was realized. The contact simulation analysis and rolling test verified the correctness of the line contact conjugate flank model and the feasibility of the cutterhead approximation machining method.


Introduction
Spiral bevel gears can realize power and motion transmission for intersecting axes, which are the key parts of mechanical transmission and are widely applied in many fields. According to the form of teeth depth, spiral bevel gears can be divided into constant teeth depth and tapered teeth depth. Spiral bevel gears with constant teeth depth can be directly machined by cutterheads to form conjugate flanks based on the generating gear principle. However, when spiral bevel gears with tapered teeth depth are machined in the same way, the cutterhead axis should be perpendicular to the root cone instead of the pitch cone. As a result, the axes of the two cutterheads are not parallel when machining the gear and pinion, that is, both cutter surfaces cannot match each other, which makes the machined gear pairs unable to achieve conjugate meshing in theory.
Based on the local conjugate theory [1], with the aim of local contact for the mating flank, traditional Gleason's technology realized the local conjugate meshing of spiral bevel gears. However, problems such as difficulty in adjusting the contact area and poor bearing capacity of gear pairs remained. For the above problems, Gleason Corporation successively proposed the helix form method [2] and duplex helical method [3], which adjusted the position of the contact area by adding the helical feed motions to the gear and the pinion, respectively, during cutting. Litvin et al. [4][5][6] presented the local synthesis, which overcame the problem where local conjugate theory did not control the second-order contact parameter of the tooth flank by presetting three second-order contact parameters at the reference point.
In recent years, with the development of gear machining technology, many effective tooth flank design and machining theories have been proposed. Shih et al. [7,8] developed a novel ease-off methodology for flank modification, which constructed the desired ease-off topographies by calculating and synthesizing the ease-off along the contact path from the predesigned transmission error and along the contact line from the predesigned bearing ratio. To identify the values of the machine tool settings required to obtain flank modifications in hypoid gears, Artoni et al. [9][10][11][12] solved the nonlinear least-squares formulation by the Levenberg-Marquardt method with a trust-region strategy. On this research basis, an efficient methodology was presented to restore the designed functional properties of hypoid gears, that is, the gear deviations can be mapped into equivalent pinion deviations, added to those of the pinion itself, and cumulatively compensated for by applying corrective machine tool settings to the pinion. Nie et al. [13,14] established the pinion target tooth flank and modification mathematical model by adjusting the factors of tooth flank mismatch topography along five directions, which improved the meshing performance and correct contact area of spiral bevel gears. Mu et al. [15] modified the tooth surface by using an arc blade instead of cutterheads, which reduced the maximum contact stress on the tooth surface and avoided tooth edge contact under misalignment or heavy load. Based on the highly flexible universal five-axis milling machine tools, Alvare et al. [16] studied the complete machining process of the large-sized spiral bevel gear and found the optimal cutting strategy in all strategies with different milling paths and scallop heights. Zhou et al. [17] reduced the machining error of the contact area through the optimization method of cutter position and orientation based on the tooth flank model and NC programming technology. Sun et al. [18] analyzed the reason for the bias of tooth contact in the traditional generating method and effectively improved the quality of tooth flank by adopting the spread-out helix modified roll. Efstathiou et al. [19] established the kinematic simulation model based on the milling process and studied the effect of generating feed rate and cutting tool geometry.
However, the tooth flank modification and manufacturing technology of spiral bevel gears based on the local conjugate principle are still limited by the small contact area for improving the bearing capacity and service life. Therefore, some experts have begun to break through the limitations of traditional manufacturing equipment and machining theory and carried out meshing theory research on line contact spiral bevel gear pairs with larger contact areas.
Taking the tooth flank design of spiral bevel gears with spherical involute as an example, Hong et al. [20] proposed a new cutting method based on the theory for designing the tooth profile of spiral bevel gear, which can process spherical involute tooth profile, but the tooth profile error is large in actual machining. García-García et al. [21] obtained bevel gears with exact spherical involute based on the additive manufacturing technology, but the bearing capacity is limited and machining efficiency is low. Ailibao et al. [22] proposed a generating method for spherical involute tooth flanks based on linear approximation technology but lack machining test verification. Sun et al. [23] realized the dual machining of spiral bevel gears by using the gear machining tooth flank as the cutter surface of the pinion, but the machining efficiency is too low to be industrialized.
To solve the problem where line contact spiral bevel gear pairs cannot be machined by cutterheads, a cutterhead milling approximation machining method was proposed in this paper, that is, by controlling the geometric topological deviations between the machining tooth flank and theoretical tooth flank, the modification objective where the machining tooth flank approaches the theoretical tooth flank will be realized.
The plan for the remainder of this paper is as follows. In Sect. 2, the mathematical modeling of line contact conjugate flank was studied, and meshing equations and conjugate flank equations of bevel gear pairs were obtained. Section 3 discussed the modeling process of a geometric topological deviations model between the machining tooth flank, and the theoretical tooth flank was established. Then, the additional cutting motions and machining compensation parameters of cutterheads were obtained to control the machining tooth flank deviations in Sect. 4. Finally, Sect. 5 carried out the test verification of theoretical modeling and machining methods.

Mathematical model of line contact conjugate flanks
The mathematical model of line contact conjugate flanks was established. As shown in Fig. 1, coordinate systems S 1 (x 1 , y 1 , z 1 ) and S 2 (x 2 , y 2 , z 2 ) are rigidly connected to gear 1 and gear 2, respectively. Auxiliary coordinate systems S p (x p , y p , z p ) and S g (x g , y g , z g ) are used to determine the positional relationship between gear 1 and gear 2. According to the positional relationship between each coordinate system shown in Fig. 1, the transformation matrix from coordinate system S 2 to S 1 can be obtained: where l 1 and l 2 are the distances from the pitch apex to the crossing point of axes for gear 1 and gear 2, respectively. e is the offset. φ 1 and φ 2 are the rotation angles of gear 1 and gear 2, respectively. The rotation angle φ 1 can be represented as φ 1 = φ 2 ·z 2 /z 1 , where z 1 and z 2 are the number of teeth for gear 1 and gear 2, respectively. Let the position vector of any meshing point P on the tooth flank Σ 2 be expressed in coordinate system S 2 as follows: where v and w are the tooth flank parameters used to determine the position vector of point P.
The unit normal vector n (2) 2 at the point P is expressed in coordinate system S 2 as follows: where , and N (2) z2 are the components of normal vector N (2) 2 on axes x 2 , y 2 , and z 2 . The relative motion velocity v (2) 12 at the meshing point P in coordinate system S 2 can be represented as follows: At meshing point P, the two conjugate flanks satisfy the following meshing equation:

Equations (2) and (3) show that n (2)
2 is a function with parameters v and w, and v (2) 12 is a function with parameters v, w, and φ 2 . Therefore, Eq. (5) can be expressed as scalar functions: The equation for the instantaneous contact line on tooth flank Σ 2 can be obtained as follows: (2) (2) 2 ⋅ (2) 12 = 0 Equation (7) is expressed as a Cartesian coordinate equation: The equation for conjugate flank Σ 1 can be obtained by transforming Eq. (7) in coordinate system S 2 to S 1 : Equation (9) is expressed as a Cartesian coordinate equation:

Geometric topological deviations model
Due to the large number of teeth and longer cutting time for the gear, to ensure machining efficiency, the gear tooth flank should be set as the datum tooth flank for priority machining. First, the machine tool cutting model was established, and the gear machining tooth flank equation based on the cutterhead generating method was solved. Second, the pinion theoretical tooth flank equation, which is conjugate meshing with the gear machining tooth flank, was obtained by applying the mathematical model of line contact conjugate flanks. Then, based on the pinion theoretical tooth flank equation and the pinion machining tooth flank equation obtained by the generating gear principle, the geometric topological deviations model of the comparison between the machining tooth flank and theoretical tooth flank can be established.

Machine tool cutting model
The gear of a spiral bevel gear pair is usually machined by the forming method or generating method, and the pinion is often processed by the roll ratio modification method and tilting method. The roll ratio modification method requires the roll ratio modification mechanism to realize the regular speed changing of the roll ratio. Its machine tool cutting model is similar to the generating method, which only needs to change the roll ratio from constant to variable. The tilting method needs a machine tool with the tilt and swivel mechanisms to realize the changes in tilt angle and cutting position for the cutterhead, with the addition of tilt and swivel coordinate systems relative to the generating method.
Considering the generality of the tooth flank model, when establishing the machine tool cutting model, the gear and the pinion will be machined by the generating method and tilting method, respectively. The generating method can be regarded as a special machine tool cutting model of the tilting method with the tilt and swivel mechanisms at the initial zero position. In other words, the tilt and swivel angles are both zeros when the gear is machined by the generating method. Therefore, the cutting model of the cradle-type bevel gear milling machine with the tilt and swivel mechanisms was established, as shown in Fig. 2.

Gear machining tooth flank equation
Based on the above cutting model, the machining coordinate systems for the gear were established in Fig. 3, coordinate systems S t (x t , y t , z t ) and S 2 (x 2 , y 2 , z 2 ) are rigidly connected to the cutterhead and the work gear, respectively. S c represents the coordinate system of the generating gear, and S 0 represents the machine coordinate system. S a , S b , and S g are auxiliary coordinate systems used to determine the positional relationship between the cutterhead and the work gear. Transforming position vector r (t) t of cutter surface Σ t in coordinate system S t to S 2 , position vector r (2) 2 of gear machining tooth flank Σ 2 can be obtained in coordinate system S 2 : where n (t) t ( 2 ) and v (t) 2t ( 2 , 2 , 2 ) represent the unit normal vector of the cutter surface and the relative motion velocity between the cutterhead and the work gear in coordinate system S t , respectively. M 2t ( 2 , 2 ) represents the transformation matrix from coordinate system S t to S 2 . ξ 2 = (u 2 ,   Fig. 3 Machining coordinate systems for the gear θ 2 , r 02 , α 2 , W 2 ) are the cutter parameters of the gear. Here, the position parameters are represented by u 2 and θ 2 , the profile angle is represented by α 2 , the nominal radius by r 02 , and the cutter point width by W 2 . λ 2 = (φ c2 , S 2 , q 2 , i 2c , δ M2 , E 2 , X 2 , X B2 ) represent the machining parameters of the gear tooth flank, including the rotation angle of generating gear φ c2 , the radial setting S 2 , the initial cradle angle setting q 2 , the roll ratio i 2c , the machine root angle δ M2 , the vertical offset E 2 , the increment of machine center to back X 2 , the sliding base feed setting X B2 , and the rotation angle of the gear φ 2 , which can be represented as φ 2 = i 2c ·φ c2 . The above parameters are constant except for u 2 , θ 2 , and φ c2 . Therefore, when solving the tooth flank points, the gear machining tooth flank equation can be regarded as an expression of the variables u 2 , θ 2 , and φ c2 for iterative solution.

Pinion theoretical tooth flank equation
Transforming position vector r (2) 2 in coordinate system S 2 to S 1 , position vector r (1) 0 of pinion theoretical tooth flank Σ 0 , which is conjugated with gear machining tooth flank Σ 2 can be expressed as follows: where n (2) 2 ( 2 , 2 , 2 ) and v (2) 12 ( 2 , 2 , 2 ) represent the unit normal vector of the gear machining tooth flank and the relative motion velocity between the gear and pinion in coordinate system S 2 , respectively. M 12 ( 1 , , 2 ) represents the transformation matrix from coordinate system S 2 to S 1 . Since the transmission ratio of the spiral bevel gear pair is fixed, the pinion theoretical tooth flank equation can also be regarded as an expression of the variables u 2 , θ 2 , and φ c2 .

Pinion machining tooth flank equation
The solution process of the pinion machining tooth flank equation is similar to that of the gear, but the tilt and swivel coordinate systems are added based on the generating method, for which it is necessary to reestablish the machining coordinate systems to solve it.
Based on the cutting model in Fig. 2, machining coordinate systems for the pinion were established in Fig. 4, coordinate systems S t (x t , y t , z t ) and S 2 (x 1 , y 1 , z 1 ) are rigidly connected to the cutterhead and the work pinion, respectively.
S 0 represents the machine coordinate system. Coordinate systems S Q , S β , S J , and S I are rigidly connected to the cradle, eccentric, swivel, and tilt drum, respectively. S a , S b , S c , S d , S e , and S g are auxiliary coordinate systems used to determine the positional relationship between the cutterhead and the work pinion. Transforming position vector r (t) t of cutter surface Σ t in coordinate system S t to S 1 , position vector r (1) 1 of pinion machining tooth flank Σ 1 in coordinate system S 1 can be obtained: where ( 1 , 1 , 1 ) represent the unit normal vector of the cutter surface and the relative motion velocity between the cutterhead and work pinion in coordinate system S t , respectively. M 1t ( 1 , 1 ) represents the transformation matrix from coordinate system S t to S 1 . ξ 1 = (u 1 , θ 1 , r 1 , α 1 ) are the cutter parameters of the pinion. Here, the position parameters are represented by u 1 and θ 1 , the profile angle is represented by α 1 , and the cutter point radius by r 1 . λ 2 = (Q, β, I, J, δ M1 , E 1 , X 1 , X B1 , i 1c ) represent the machining parameters of the pinion tooth flank, including the cradle angle Q, which is the sum of the initial cradle angle Q 0 and the rotation angle of the generating gear φ c1 , the eccentric angle β, the swivel angle J, the tilt angle I, the machine root angle δ M1 , the vertical offset E 1 , the increment of the machine center to the back X 1 , the sliding base feed setting X B1 , the roll ratio i 1c , and the rotation angle of pinion φ 1 , which satisfies the relation φ 1 = i 1c ·φ c1 . Similar to the gear machining tooth flank equation, the pinion machining tooth flank equation can be regarded as an expression of the variables u 1 , θ 1 , and φ c1 .

Geometric topological deviations model
According to the obtained pinion theoretical tooth flank Σ 0 and machining tooth flank Σ 1 , a 3D model of both tooth flanks on the same pinion root cone was established. As seen from Eq. (13), pinion theoretical tooth flank Σ 0 is the secondary enveloping surface of cutter surface Σ t (cutter surface Σ t -gear machining tooth flank Σ 2 -pinion theoretical tooth flank Σ 0 ), and pinion machining tooth flank Σ 1 was directly enveloped by cutter surface Σ t (cutter surface Σ t -pinion machining tooth flank Σ 1 ) so that the variable parameters u, θ, and φ at each tooth flank point are generally not the same. As a result, the 3D model of both tooth flanks on the same pinion root cone will have an angular offset, as shown in Fig. 5.
To solve the geometric topological deviations between the two tooth flanks, pinion machining tooth flank Σ 1 should be rotated by ∆φ around its axis z 1 so that its midpoint will coincide with the midpoint of the theoretical tooth flank Σ 0 . The pinion machining tooth flank Σ 1 after rotation can be expressed as follows: here, pinion theoretical tooth flank can be solved. Then, the deviation ∆d between any point on pinion machining tooth flank Σ 1 and its corresponding point on theoretical tooth flank Σ 0 can be expressed as follows: Equation (18) is transformed into vector form:

Additional cutting motions
To enable the machined pinion tooth flank to more accurately approximate its line contact theoretical tooth flank, additional cutting motions should be added to the cutterhead in the process of pinion machining flank modification. Based on the obtained pinion machining tooth flank equation of the tilting method, the roll ratio can be modified to make the relationship between the rotation angle of the generating gear and the work pinion satisfy the modification polynomial as follows: where C 0 and D 0 represent the second-order ratio and thirdorder ratio, respectively, and their initial values can be set to zero. Substituting Eq. (20) into Eq. (13), the pinion machining tooth flank equation can be reestablished: The geometric topology deviations model and the additional cutting motions of the cutterhead were obtained. Then, with the pinion machining tooth flank approaching its theoretical tooth flank as the goal, the machining compensation parameters of the cutterhead were obtained to control the machining tooth flank deviations and reduce them to the allowable deviations.
Since the nominal cutter radius and profile angle have been standardized, cutter parameters ξ 1 should be kept unchanged without increasing the cost of cutter manufacturing, and only machining parameters λ 1 should be modified. Therefore, the deviations should be simplified from the expression d (ξ 1 , λ 1 ) to expression d (λ 1 ) containing only parameters λ 1 , and the deviations of all the tooth flank points can be expressed as a vector: Taking machining parameters λ 1 as the optimization variable and the minimum square sum of the deviations d (λ 1 ) as the goal, the least square method optimization model can be established: The Levenberg-Marquardt method with a trust-region strategy is used to solve the optimization model. Generally, the initial deviations between the pinion machining tooth flank and its theoretical tooth flank are relatively large. Directly using the iterative method to modify machining parameters λ 1 will cause a large offset of the cutter relative to the initial position, that is, interference or deviation with the work pinion. To make the modified machining parameters have practical engineering significance, the initial deviations can be reduced by setting a deviation default value in advance, that is, some of the machining parameters are modified in advance so that the cutterhead and the work pinion do not offset along the z axis of the machine.
The influence of each tooth flank machining parameter on the relative position between the cutter and work pinion is analyzed. As shown in Fig. 7, the parameters that can prevent the cutter and the work pinion from being offset along the z axis of the machine include Q, β, E 1 , i 1c , C 0 , and D 0 . Therefore, the machining parameters to be modified in advance are ρ 1 = [Q, β, E 1 , i 1c , C 0 , D 0 ].
Then, with premodified machining parameters, ρ 1 as the optimization variable and the minimum square sum of tooth flank deviations d (ρ 1 ) as the goal, the least square optimization model can be established as follows: The premodified machining parameters that reduce the deviation below the default value are solved, and the premodified deviations are obtained as d ρ (λ 1 ). Then, the minimum sum of squares of d ρ (λ 1 ) is taken as the goal and substituted into Eq. (23) to reestablish the optimization model: with such step-by-step modification by the above method, the machining compensation parameters that can satisfy the requirement of controlling the deviations and reducing them to the allowable deviations can be obtained to realize the effective approximation of the pinion machining tooth flank to its theoretical tooth flank. In addition, the approximation process is shown in Fig. 8. where f a (φ c ), f b (φ c ), f x (φ c ), f y (φ c ), and f z (φ c ) represent the motion functions of φ a , φ b , D x , D y , and D z , respectively, with the rotation angle of generating gear φ c . By setting the value of φ c and substituting it into Eq. (29), the instantaneous machine point position corresponding to each axis of the NC machine tool can be obtained, and then the cutting motion planning of the line contact spiral bevel gear pair can be realized.

3D model of the spiral bevel gear pair
Taking a pair of spiral bevel gear pairs with tapered teeth as the example pairs, the case analysis and cutting test were carried out on the drive side (gear convex and pinion concave). Table 1 shows the basic parameters of the example pair, from which the initial machining parameters of the (29) example gear and pinion are determined and are listed in Tables 2 and 3, respectively.
Substituting the parameters in Tables 2 and 3 into Eqs. (12) and (13), respectively, generates the 3D model of the pinion machining tooth flank and its theoretical tooth flank by discretizing the tooth flank. The geometric topology deviations model of the initial pinion machining tooth flank and its theoretical tooth flank was established, as shown in Fig. 10. Calculating the deviation value at each effective tooth flank point, the results in Fig. 10 show that the average value of the tooth flank deviations is 396.9 μm, the sum of squares is 11,566,311.1 μm 2 , and the absolute value of the maximum is 874.5 μm.
The data show that the deviation between the initial pinion machining tooth flank and its theoretical tooth flank is large, so it needs to be modified in advance based on some of the tooth flank machining parameters. The premodified machining parameters are listed in Table 4.
According to the machining parameters in Table 4, the geometric topology deviations model of the premodified   pinion machining tooth flank and its theoretical tooth flank was established. As shown in Fig. 11, the average value of the tooth flank deviation is 18.1 μm, the sum of squares is 32,131.5 μm 2 , and the absolute value of the maximum is 101.7 μm. The data show that through previous modification, the tooth flank deviations were effectively reduced on the premise of ensuring no offset between the cutterhead and the work pinion. Then, all the machining parameters of the pinion tooth flank were modified, and the modified machining parameters are listed in Table 5.
Based on the parameters in Table 5, the geometric topology deviations model of the modified pinion machining tooth flank and its theoretical tooth flank was established. As Fig. 12 shows, the average value of the deviations is 0.6 μm, the sum of squares is 31.7 μm 2 , and the absolute value of the maximum is 2.8 μm.  The results show that, based on the cutterhead approximation machining method mentioned in this article, the pinion machining tooth flank can be effectively modified to the allowable deviations of its theoretical tooth flank by adjusting the cutting motions of the cutterhead so that the pinion machining tooth flank can effectively approach its line contact theoretical tooth flank.

Contact simulation analysis
Contact simulation analysis was carried out on the example pair, and the main parameters of the finite element analysis model were set as follows: (1) the material properties of the gear and pinion are the same, with an elastic modulus of 2.06 × 10 5 MPa, a Poisson's ratio of 0.3, and a density of 7.85 × 10 −9 t/mm 3 . (2) The angular velocity of the pinion  13 Instantaneous contact line between pinion concave and gear convex is 1 rad/s, and the resistance torque applied to the gear is 1 N·m. As shown in Fig. 13, the analysis results show that the gear and the pinion experience in line contact during the instantaneous engagement.
Based on the finite element results, the instantaneous contact area of the pinion concave from the mesh-in tooth to the mesh-out tooth was obtained. And the total theoretical contact area of this pinion concave was obtained by projecting Step time / s  each instantaneous contact area into the same shaft section of the pinion. As shown in Fig. 14, it can be seen that the example pair is in full-tooth flank contact.

Cutting and rolling tests
To further verify the correctness of the tooth flank mathematical model and machining method, the cutting and rolling tests were carried out. Figure 15 shows the cutting test processes of the pinion and gear based on a JCB32 five-axis CNC machine. Figure 16 shows the actual contact area of the example pair obtained by a Y9550 bevel roll tester. The result of the rolling test in Fig. 16 shows that the machined example pair is in full-tooth flank contact, which is consistent with the contact simulation result. The test results show that the cutterhead approximation machining method based on controlling topological deviation can machine line contact spiral bevel gears under the conditions of existing machining equipment.

Conclusion
To solve the problem where line contact spiral bevel gear pairs cannot be machined by cutterheads, this paper proposed the cutterhead milling approximation machining method based on controlling the geometric topology deviations, which allows the conjugate flanks machined by cutterheads to effectively approach the line contact theoretical tooth flanks. The conclusions are as follows: 1. A geometric topological deviations model between the pinion machining tooth flank and its theoretical tooth flank was established, which realized the comparison between the line contact theoretical tooth flank and the machining tooth flank. First, based on the mathematical model of line contact conjugate flanks, the gear tooth flank equation machined by the cutterhead generating method and the theoretical tooth flank equation of the pinion conjugated to the gear were obtained. Then, the geometric topological deviations model between the obtained pinion theoretical tooth flank and the pinion machining tooth flank machined by the tilting method was established. Finally, the deviation of each tooth flank point was calculated by topological differential calculus analysis. 2. Cutterhead milling approximation machining method of line contact spiral bevel gear pairs based on controlling topology deviations was proposed. That is, with the pinion machining tooth flank approaching its theoretical tooth flank as the modification, by adjusting the cutting motions of the cutterhead, the pinion machining tooth flank was modified to the allowable deviations range of its theoretical tooth flank, and the cutterhead milling approximation machining of line contact spiral bevel gear pairs was realized. 3. The contact simulation results show that the example pair obtained by the cutterhead milling approximation machining method is in line contact. The cutting and rolling tests obtained the gear pair with full-tooth flank contact. The test results show that the cutterhead approximation machining method based on controlling topological deviations can machine line contact spiral bevel gears.
Author contribution All authors contributed to the study's conception and design. Methodology, data analysis, and validation were performed by Mingyang Wang. Resources and supervision were performed by Yuehai Sun. The first draft of the manuscript was written by Mingyang Wang, and all authors commented on previous versions of the manuscript. All authors read and approved the final manuscript.
Funding This work was supported by the National Natural Science Foundation of China (Grant Numbers 51875395).
Availability of data and material All data generated or analyzed during this study are included in this article.

Competing interests
The authors declare no competing interests.