Mathematical modelling of power skiving for general profile based on numerical enveloping

Power skiving is an effective generating machining method for internal parts like gears with respect its high productivity. The general mathematic modelling for power skiving is the basis for cutting tools design, machining precision evaluation, and machining process optimization. Currently, mainly studies are focus on the involute gear machining with adopting the analytical enveloping equation. However, these analytical methods have failed to deal with overcutting for general profile skiving tasks. Moreover, little attention has been devoted to investigate the power skiving process with taking variable configuration parameters, which is significant to control the machined surface topography. Herein, we introduce a mathematic modelling method for power skiving with general profile based on the numerical discrete enveloping. Firstly, the basic mathematic model of power skiving is established, in which the center distance is formulated as polynomial of time. By transforming the power skiving into a forming machining of the swept volume of cutting edge, a numerical algorithm is designed to distinguish the machined transverse profile via the discrete enveloping ideology. Especially, the precise instant contact curve is extracted according to the feed motion speed inversely. Finally, simulations for involute gear and cycloid wheel are carried out to verify the effectiveness of this method and investigate the influence of variable radial motions on the machined slot surface topography. The results validate capability of this method on simulate the dynamic power skiving process for general desire profiles and evaluate the machined results.

Rotation angle of cutter φ w Rotation angle of workpiece Ω Section angle between axes of cutter and workpiece E 0 Center distance between cutter and workpiece E c Install eccentricity of cutter f w Section plane of workpiece N t Dividing number of time N p Dividing number of profile N s Dividing number of cutting edge K e (t) Polynomial of time for center distance R c Pitch circle radius of cutter

Introduction
Power skiving is a typical generating machining method, which provides high productive in machining parts with periodic features like gears. With installing as a pair of cross-axis gear meshing, the cutting edges on the end face of cutter remove a layer of material via the relative axial motion for the part [1]. This working principle indicates that power skiving not only adopts continuously generating motion as gear hobbing to ensure high machining efficiency, but also retains the advantages of gear shaping that is capable to manufacture both external and internal part, non-through slots, doublelinked gears, and so on. As invented in 1910 by Pittler [2], however, power skiving underwent slow development with subject to the poor stiffness of the machine tools and the short tool life. In recent years, with the development of tool materials, especially, the spindle technology and numerical control systems, the power skiving has performed superiority in gear production. Relevant power skiving solutions [3][4][5] such as machine tools and cutters have been provided by companies like Gleason and Pittler. Meanwhile, theories on power skiving have been studied in difference aspects to enhance its application.
As the basis of successful machining, mathematical model of power skiving and the cutting edge curve were studied widely [6][7][8][9][10][11][12][13][14][15][16][17][18][19]. Jin [6] investigated the analytical theory of gear skiving and pointed out that the action line of skiving is as same as the action line of spiral gears. Li et al. [7][8][9] analyzed the working principle of skiving and proposed a cutter design approach to error-free gear skiving. Guo et al. [10,11] investigated the skiving tool design and cutting mechanism of cylindrical gears, and introduced a multiple blades taper skiving tool [12]. Radzevich [13] introduced the design and computation principle of the skiving cutter for gear skiving. Tsai [14] established the mathematical model for design and analysis of power skiving tool for involute gear cutting. Stadtfeld [15] introduced the power skiving technology of Gleason, including the generation kinematics, cutter geometry, chip geometry, machine tool configuration, and processing software. Tomokazu et al. [16] established a calculation model for internal gear skiving with a pinion-type cutter having pitch deviation and run-out. Moriwaki et al. [17,18] investigated the cutting tool parameters of cylindrical skiving cutter with sharpened angle for internal gears to optimize skiving cutter design. Shih and Li [19] proposed an error-free conical power skiving cutter design method via meshing theory with considers the variation of center distance. Jia et al. [20] developed a discrete enveloping-assisted cutting edge curve calculation method for skiving cutter. Li et al. [21] studied the design of power skiving tool by considering the interference and minimizing the machining deviations. These researches provided the general theory and mathematical model for power skiving as well as the calculation method of the cutting edge curve design.
In addition, the simulation of power skiving process is significant to estimate and optimize the skiving process and cutting tools. Research works are mainly devoted to geometric precision simulation [22][23][24][25][26][27][28][29][30][31] and the working condition analysis of cutter [32][33][34][35][36][37]. Commonly, the geometric simulations include three kinds of manner. First, the analytical ones adopted the meshing equation for geometry calculation, like Guo et al. [22,23] studied the theoretical tooth profile errors of gears in skiving and investigated the tool setting errors on gear skiving accuracy; further, they [24] studied the cutting edge correction method for conical cutter; Zheng et al. [25] generalized the machine kinematics correction and TCA to gear skiving. Secondly, the numerical methods represent the machining process in a discrete way for concrete calculation, like Jia et al. [26] developed a numerical simulation method for non-involute gear power skiving with adopt approximate external enveloping; Zheng et al. [27] developed a novel z-mapbased numerical method to calculate tooth flank and investigate the influence of eccentricity error on the surface roughness; Inui et al. [28] developed a triple-dexel based geometric simulation method for power skiving process with GPU computing acceleration. Thirdly, the CAD-based methods perform the cutting process as Boolean operation and finished it via commercial CAD package, like Antoniadis et al. [29,30] simulated the kinematics of the cutting process with the aid of commercial CAD software, which allows the precise determination of the non-deformed chips and cutting forces; Tapoglou [31] simulated the chip thickness of skiving with the aid of CAD. In the basis of geometric simulation, cutting force of power skiving is modelled via dividing the cutting edge as a serial of oblique cutting element and integrating the dynamic force of all the elements. McClosky et al. [32] and Onozuka et al. [33] reported that these modelling results are consistent with the experiment cutting force; Li et al. [34,35] further discussed the temperature during power skiving process. Besides, some comprehensive simulations were developed such as Klocke et al. [36] investigated the influences of skiving configuration parameters on the working performance like tool wear and chip welding, and Schulze et al. [37] studied the kinematic process of skiving and investigated the chip formation mechanisms of skiving with adopting 3D-finite element simulation. Most of the researches are focused on the machining of involute gears with adopting analytical mathematic model. However, it is difficult to process the conditions like self-intersection and overcutting, particularly for universal target profiles. The discrete numerical enveloping method proposed in [26] demonstrates excellent performance to deal with these drawbacks in general profile skiving simulation, but it lacks the capability to investigate the entire skiving process with combining time-varying configuration and kinematic parameters.
Aiming to investigate the dynamic machining process for power skiving by overcoming the shortcoming of [26], this work developed a discrete enveloping-based mathematic modelling method for power skiving, which is the effect of the general profile tasks with strong robustness. The remainder is organized as follows. The basic mathematic model for skiving is introduced in Section 2. The numerical simulation method for power skiving is studied in Section 3. Then in Section 4, the deviation estimation is described. In Section 5, several skiving tasks are simulated and concluded at last. 2 The mathematical model for power skiving

Configuration of power skiving
In a general power skiving system, the cutter and workpiece is set up as a pair of cross-axis in Fig. 1. Workpiece coordinate system S w : O w -X w Y w Z w is attached to the workpiece, and its Z w -axis is coincided with the cylindrical workpiece axis. Similarly, cutter coordinate system S c : O c -X c Y c Z c is attached to the cutter, and its Z c -axis is coincided with the cutter axis. In initial, we define that X w -axis is going through the origin O c and X c -axis is coinciding with X w -axis. Meanwhile, with respect to the meshing of cutter and workpiece, the shaft angle Ω between their axes and the nearest distance E 0 between their axes are given as follows: where j w denotes the helix direction of workpice (j w = 0 indicates spur, j w = 1 indicates right hand, and j w = − 1 indicates left hand), and j c denotes the helix direction of cutter (j c = 0 indicates spur, j c = 1 indicates right hand, and j c = − 1 indicates left hand); β w and β c are the helix angles on the pitch circles, respectively, for the workpiece and the cutter; and k io denotes the type of skiving (k io = 1 denotes external skiving, and k io = − 1 denotes internal skiving). Moreover, in practical applications, the cutter takes an eccentricity E c along the Z c -axis from the origin O c , looking to improve the working condition of cutting edges and to avoid the interferences during the machining process.

Motions of power skiving
Power skiving consists of two kinds of coupled motion as illustrated in Fig. 1, i.e., the meshing motion and the feeding motion. The cutter and workpiece rotate φ c and φ w around Z caxis and Z w -axis, respectively, with constant transmission ratio, which is termed as main meshing motion, contributing to material cut off. Meanwhile, a differential rotation Δφ c is implemented on the cutter to ensure the meshing motion when the cutter performs a synchronous linear feed motion f along Z w -axis, which contributes to produce the complete slot. The rotation angles of the power skiving system are given as below: where z w and z c are the slot number of workpiece and teeth number of cutter, respectively; Δφ c is the differential rotation angle of the skiving cutter, and it is determined by both the axial feeding Δf of the skiving cutter along Z w -axis and the configuration of power skiving system as follows: Aiming to satisfy the diverse requirements in practical applications of skiving, such as tooth profile crown correction and machining error sensitivity analysis, we involve an additional feed motion K e to enhance this virtual kinematic model of power skiving as in [38], which expresses the constant radial distance E 0 between the cutter and the part. The variances of feed motion K e is formulated as polynomial of time as follows: where q is the polynomial order of the radial feed, c e is the polynomial coefficients of radial feed, and ΔE is the variance in radial direction. Furthermore, proper rotation directions of the cutter and the workpiece are crucial to ensure successful material removal during skiving. The rotation of cutter must ensure that the cutting edge is sliding inside the slot of workpiece along the feeding direction.
3 Mathematical modelling of power skiving process by numerical enveloping

Equivalent expression of power skiving as form machining
The aforementioned working model of power skiving indicates the engagement of each tooth cut off a layer of material in slot with the help of both meshing motion and feeding motion, while the succession of engagements produces the whole slot via the feeding motion. This revels that the skiving process can be taken as a forming machining process for helical or spur gear like drawn in Fig. 2, in which the swept volume of one cutting edge (SVC) relative to the gear slot during engagement works as the wheel and its instant contact curve generates the slot surface following the feed motion.
Consequently, one can deduce the generating process of power skiving as follows: each cutting edge generates one or several points on the desired surface at every engaging moment, and these points further develop a forming curve on the desired slot surface during each engagement of one cutting tooth.

Numerical enveloping of power skiving
Commonly, the generating points during power skiving process are determined by enveloping equations. It tells that the external normal of generating points on the SVC surface is perpendicular to the velocity vector of their feed motion. However, the analytical method might invalid in cases like singular points, overcutting, and multiple cutting edge. For this reason, this work models the skiving process as a forming machining process in a numerical way. The key points are approximating SVC by a serial of cutting edges, and then obtaining the transversal profile of workpiece and the contact curve via identifying the external enveloping curve of the intersection curves between SVC and the workpiece transversal section.
As shown in Fig. 3, the detailed numerical modelling and investigating for the skiving process are arranged in 7 steps.
Step. 1: Dividing the cutting edge curve r(u) = [x(u),y(u),z(u),1] T into continuous points for the purpose of maintaining the generality of this method, where u is the radial parameter of the cutting edge.
Step. 2: Determining the rotation interval [φ s , φ e ] of one tooth engaging cutting by figuring out the rotation angles of cutter touching the addendum circle of worpiece.
Step. 3: Generating the cutting edge curve r w (u,t) in workpiece system S w based on the mathematical model of power skiving as given as below: where M w c (t) is the transform matrix from cutter system S c to workpiece system S w at t moment. It is given by basic homogeneous coordinate transform matrices according to the kinematic model of skiving.
Step. 4: As shown in Fig. 3a, establishing the swept volume of one cutting edge SVC through assemble N t cutting edges r w (u,t i ) with rotation angle uniformly within where t i is the time moment of the i th rotation angle. It is given by the angular speed of cutter ω c as follows: Step. 5: As shown in Fig. 3b, taking SVC to perform a forming machining on the transversal section of workpiece f w (z = z 0 ). Essentially, it is intersecting all the cutting edge curve r w (u,t i ) of SVC to the transversal section f w along the feed motion f(τ), which is formulated as an independent motion as the function of τ. In fact, additional rotation Δφ c (τ) is involved to ensure that the intersection curves are symmetrically distributed along the X w -axis. In this basis, one can figure out the intersect point E i,j for each i th cutting edge with specified radial parameter u j which satisfies Eq. (10).
Obviously, the proper τ 0 can be ascertained by numerical searching methods, since Eq. (10) is a function of variable τ, in which t i specifics a constant locus of cutting edge curve. Consequently, the point E i,j can be directly given as: Step. 6: Distinguishing the external envelop profile of all the intersect curves on f w . The external enveloping profile can be derived numerically as a serial of points along the X w -axis since the intersection curves are distributed around X w -axis. At each radius within the specified radial range [r f ,r a ], setting a line perpendicular to X-axis, then its intersection point for every intersect curve on f w , is easy obtained as the intersection of two straight lines segment in X-O-Y plane. After all the intersection points between the intersect curves and current line have been identified, the point that taking the maximum/minimum y-component is specified as the enveloping profile points C k as shown in Fig. 3c. Concurrently, for each profile point, the time serial number of corresponding cutting edge indicates the cutting occur time of machined profile point C k .
Step. 7: Extracting the instant contact curve on the machined slot surface as in Fig. 3d. At first, transforming each machined profile point C k inversely to the rake flank S R that defined by the cutting occur time t k according to the feed motion and additional rotation Δφ c (t). The intersection point Q k that cutting occurs on the cutting edge can be solved via searching proper variable τ p satisfying Eq. (12).
Subsequently, the instant contact point P k on the machining surface can be determined by translating Q k from cutter coordinate system to workpiece coordinate system based on its corresponding time moment t k as follows: In final, all the instant contact points on the slot surface compose a spatial curve, which is the generating line of one cutting edge on the slot surface during one engagement of skiving process.

Simulation frame of power skiving
Simulating the complete slot surface of workpiece undergoing power skiving is significant to analyzing the machining result like gear crowning. For this purpose, the workpiece is assembled by a serial of disks with same thickness along axial direction. Meanwhile, only one slot is simulated with respect to the symmetric of workpiece. The resulted generating line of cutting edge that works on each axial plane of workpiece composes the machined slot surface.
The flowchart of this numerical simulation for power skiving is provided as in Fig. 4; the concrete seven steps for calculating the generating line are demonstrated as a submodule. The simulation frame is implemented to following case studies.
In practical machining, the cutting edge and the slot of workpiece are meshing periodically; consequently, the axial displacement of every cutting engagement must comply with this meshing motion. For simplicity, we are taking the assumption that the axial position of f w can be given arbitrarily in this study, without considering the cusps on the machined slot surface that are produced by the inherent interrupted meshing cutting mechanism.

Machined error evaluation
In order to evaluate the simulation precision of proposed method, the deviation between the simulated profile and the theoretical one is calculated. As shown in Fig. 5, for each generating curve of every cross section on the workpiece, every generating point decides a cross section. The shortest distance from the generating point to the point on the theoretical profile on this cross section is specified as machining error, instead of the minimum distance from the simulated point to the theoretical surfaces along its normal vector. Therefore, the machining error is expressed along the generating curve on the slot surface.

Numerical examples
Aiming to demonstrate the effectiveness of this method, simulations of power skiving with various configurations for an internal involute gear and a cycloid wheel are performed, and the machining errors are investigated.

Simulation for internal gear power skiving
For the convenience of evaluating the accuracy of this method, an internal helical involute gear power skiving was simulated with the parameters listed in Table 1, since the cutter and the gear blank are meshing as a pair of standard involute gear with the same normal module and normal pressure angle. At first, the standard involute curve like cutting edges for two types of rake flanks was calculated by a method in [20]. Then, the cutting edge curve and the enveloping curves on the gear transversal section for case 1 are demonstrated as in Fig. 6(a), in which the distinguished enveloping curves were also involute curves in Fig. 6b. As shown in Fig. 6c, the cutters for case 1 and case 2 adopted different rake flanks. The swept volumes of cutting edge relative to the gear blank for case 1 and case 2, which are demonstrated respectively in Fig. 6d and Fig. 6e, show difference as well as the geometries of their instant contact curve.
Compared with the standard involute curve, the deviations of machined profile for these two cases are illustrated in Fig. 7, in which the maximum deviations were no more than 0.15  μm. The accuracy indicates that the developed method is capable to simulate the machining error.
In practical machining process, various teeth surface like crowning and conical might be adopted to avoid the interferences on the part gear shoulder and to obtain specific working performance. In common, the teeth surface can be corrected by properly changing the polynomial coefficient of cutting path along the full axial feeding range. Unlike traditional analytical methods, the developed numerical enveloping simulation method was used to study the machined teeth surface by power skiving with two types of center distance K e (t) for both the two cutting edge configurations.
A linear radial motion K e (t) = E 0 +0.0067t-0.04 and a quadratic radial motion K e (t) = E 0 -0.0025t 2 + 0.03t-0.09 are carried out on the power skiving cutter as shown in Fig. 8a and b, respectively. A constant axial feed speed 1 mm/s was performed for the full axial length that started at L = − 6 mm and exited at L = 6 mm. In Fig. 8, the theoretical tooth surface is represented by the gray grid, and the machined tooth surface is defined by the contact curve as blue grid. One can find out the teeth surface performed a taper shape after linear radial motion, in which the deviations were about 24 μm at L = 6 mm, and the deviations were nearly − 23 μm at L = − 6mm, and both sides of the tooth surface showed the similar topography but following their contact curve, respectively. The simulation for the quadratic radial motion is provided in Fig. 8b, where the machined tooth surface performed a parabolic variation over the full length of tooth. The deviations were about − 34 μm at the two ends of the tooth surface and approached to the minimum nearly 0 μm at the middle. Besides, the small differences of the corresponding grid nodes of tooth surface between case 1 and case 2 are reasonable, since the distributions of contact curve are different in these two cases are different. In all, these results were consistent with the anticipated correction for teeth; i.e., linear polynomial radial motion generates a cone-shaped correction, and quadratic polynomial radial motion produces a drum-like correction.
Simulations implemented by commercial software VERICUT are provided to proof the effectiveness of developed numerical simulation method. As shown in Fig. 9, a raw gear blank is fixed for machining with specific cutting edge curve, and the linear radial motion and quadric radial motion with cutting edge as case 1 are performed, respectively. Through the color map of machining residual of the slot surface, without considering the machining residual cusps, one can see that the top section curve shows undercut about 25 μm, while the bottom section curve shows overcut about 20 μm, consisting with the simulation in Fig. 8. Meanwhile, the quadric radial motion demonstrates overcut about 25-30 μm on the two end section while zero-machining error in the middle section of the slot. The same trend as provided in Fig. 8 indicates that the developed model is capable to simulate the power skiving machining process with radial motion like gear crowning.

Simulation for cycloid wheel power skiving
Looking to prove the generality of developed method, power skiving for cycloid wheel that is adopted in RV reducer was simulated. The parameters for the short epicycloid equation, the power skiving configuration, and the numerical simulation are listed in Table. 2.
The whole simulation process is demonstrated as in Fig.  10. Firstly, according to the equidistance line of short epicycloid shown in Fig. 10a, the cutting edge curve is calculated by [20], and the external power skiving is set up as in Fig. 10b. Then, the machined profile is distinguished through enveloping all the intersection curves on transversal of wheel as in Fig. 10c. Through inversely tracing the intersection point on the rake flank specified by the recorded time serial number of each profile point, the instant contact point is extracted in work coordinate system. Consequently, the complete contact curve and the swept volume of cutting edge are obtained as in Fig. 10d.
Additionally, to verify the robustness of this method, a complex machining configuration with cubic radial motion K e (t) = 0.000125t 3 -0.0045t 2 + 0.054t-0.216 for t = [0, 12] and K e (t) = − 0.000125t 3 + 0.0045t 2 -0.054t + 0.216 for t = (12,24) is applied on the power skiving cutter for the cycloid wheel machining as shown in Fig. 11. The theoretical tooth surface was represented by gray grid, and the machined surface is illustrated by contact curves in blue grid by sampling nodes. The deviations on the two ends of tooth surface were maximum and were reduced to the minimum near to the middle at each axial path (from A to K). It is consistent with the desired drum-like tooth correction via cubic polynomial radial motion.
As the presented simulations shown, the proposed mathematic model for general profile power skiving is capable to simulate the skiving motion and to evaluate the results of power skiving, and it may be useful for power skiving CNC machine developers as well as the gear manufacturers. Besides, properly changing the polynomial coefficient of feed motions is a solution of gear correction for power skiving.

Conclusions
This work developed a mathematical model of power skiving with strong robustness to simulate the instant contact curve for dynamic machining process and general profile tasks via adopting numerical enveloping ideology. This model is verified and applied by simulations with polynomial radial motions. The main conclusions are drawn as follows: (i) The proposed mathematic model for power skiving includes not only the general configuration parameters but also the center distance formulated as polynomial of time. It is effective for diverse dynamic skiving process modelling rather than the traditional constant kinematic model. Besides, this method takes a strong robustness thanks to the numerical discrete enveloping manner, which ensures that this method is successful in dealing with general profiles like cycloid and in overcoming the overcuts in analytical methods.   the contact curves in axial distribution since the numerical calculation principle. The simulations for involute gear and cycloid wheel machining validate that the machining deviations on teeth surfaces are consistent with the assigned polynomial radial motion on the center distance as well as the simulation results via VERICUT, i.e., linear, parabolic, and cubic, which indicate that properly changing the polynomial coefficient of radial motions is a solution of gear correction for power skiving. (iv) As a numerical method, its accuracy will be affected by the number of profile point and the number of locus curves. Moreover, this method takes the assumption that cutting tooth performs continuous axial feeding motion along one slot, leading to it lack to identify the residual machining height on teeth surface. Further study will be devoted to involving the periodical meshing of cutting teeth during skiving process.