Mathematical Model for Tapping Simulation to Predict Radial Pitch Diameter Difference of threads

As the most important component of parts, thread has a great inﬂuence on the mechanical properties and service performance of parts. In order to ensure the quality of the thread, the thread quality inspection standard involves 11 main thread characteristic-s, of which the surface roughness of the thread is the most studied, and the research on Radial Pitch Diameter Diﬀerence (RPDD) is still blank. In this paper, a quasi-static model of the tapping process is developed based on the roundness error mechanism of the hole, which includes cutting force and cutting damping force. Due to the regenerative nature of cutting, the force on each cutting edge depends on both the tool’s current position and previous position. According to the eigenvalues and eigenvectors of the discrete state-transition matrix, RPDD is ﬁnally determined, and the inﬂuence of the chamfer length and the spindle speed on RPDD is simulated by this model. The results demonstrate that the chamfer length and spindle speed will aﬀect RPDD, and the RPDD is the smallest when the chamfer length is 2 threads and the spindle speed is 1400 rev/min. The development of this model not only provides a cheap and eﬀective method for the study of RPDD, but also lays a foundation for further experimental research.


Introduction
In modern industry, threaded connections are often used when parts and pipe connections need to be assembled in a non-destructive manner. The performance of threaded connections is critical to many fields such as petroleum, aerospace, shipbuilding, high-speed rail, nuclear energy, and automobile manufacturing [1][2][3].According to statistics, threaded connections generally account for 60% of the total mechanical components in each machinery and equipment [4].Therefore, it is of great significance to improve thread processing quality [5].
Poor thread quality will affect the mechanical properties and service performance of the components, such as tensile strength, torsion strength, vibration resistance, connection reliability, pipeline tightness, etc. [2,3,6].The factor that has the greatest influence on the mechanical properties and service performance of threaded components is the geometry and dimensional accuracy of the thread [2].Dong [2] and Leon [7] studied the influence of thread dimensional conformance on the selfloosening resistance and static strength of the threaded connection which will be greatly reduced due to the unqualified outer diameter, middle diameter or inner diameter of bolts and nuts. Nassar et al. [8] showed that unqualified thread root radius will adversely affect the fatigue performance of pre-tightened threaded fasteners. In order to ensure the quality of thread, the thread quality inspection standard involves 11 main thread characteristics, and RPDD is one of them. RPDD is defined as the maximum difference among the pitch diameters in all radial directions within a lead. When a thread with an excessive RPDD is matched with a qualified thread, the sharply changing contact area between the internal and external threads may cause uneven load distribution on the thread surface, and severe stress concentration occurs at the position with a smaller pitch diameter, which accelerates wear and tear, shorten the service life, and severely cause "looseness" or even "slip buckle" [3].
Compared with external threads, it is more difficult to ensure the quality of internal threads, especially small diameter internal threads. The machining of internal threads is very complicated, which is usually the last process of workpiece manufacturing. Any machining failure or reduced accuracy can not achieve the perfect assembly of components without gaps, and even lead to huge economic losses [9,10]. Tapping is a process widely used in the manufacturing of internal threads. Although in the past few years, other processing methods have been used to manufacture internal threads with great success, such as thread milling and turning [11][12][13][14][15][16], tapping is almost the only method to manufacture small diameter internal threads [17]. However, there are relatively few literatures on thread quality in tapping. Fromentin et al. [18] showed that the use of suitable oil-based lubricants in tapping can improve the mechanical properties and surface integrity of internal threads. Piska et al. [19] compared taps with PVD and TiN composite coatings and uncoated taps to process C45, and the results showed that the use of coatings can effectively reduce tool damage and improve the surface finish of threads. Hsu et al. [20] used the Taguchi method to test and analyze the influence of tool parameters and cutting conditions on thread quality during tapping, and found that the larger the helix angle, the worse the quality of the thread. Bratan et al. [21] developed a tool to improve the threading process considering the combined deforming-cutting tap. The results have shown that utilising proposed combined deforming-cutting taps in the processing of internal threads with a small diameter for M3-M6 in aluminium alloys enhanced the accuracy and surface quality of threads. It can be seen that the current literatures on the thread quality obtained by tapping is mostly focused on the influence of various factors on the thread surface quality, and the research on the RPDD of the thread is almost blank. Therefore, it is very necessary to analyze the factors that affect the dimensional of R-PDD by establishing a mathematical model of RPDD.
For a standard internal thread, its pitch diameter can be determined by the nominal diameter and pitch of the thread, that is, pitch diameter = nominal diameterpitch × 0.6495.Therefore, it can be considered that R-PDD of the thread is the result of the combined effect of the roundness error of the root circle and the pitch error. The mechanism of roundness error of the root cycle is similar to that of the hole. From the literature [22], it is known that the roundness error of the hole is largely caused by the vibration of the drill. The tap usually shows transverse vibration, torsional vibration and axial vibration during tapping, but there is no coupling relationship between the transverse vibration and axial/torsional vibration of the tap [23]. Because the torsional stiffness and axial stiffness are much greater than the bending stiffness, the torsion and axial displacement of a cantilever tool (such as a tap) with a large length-to-diameter ratio is much smaller than the transverse displacement under a given force. Therefore, torsional vibration and axial vibration are ignored, that is, the pitch error generated in tapping is ignored. The pitch diameter of the thread is only related to the nominal diameter of the thread, which means that the R-PDD is caused by the roundness error of the root cycle. Bayly et al. [24] established a quasi-static reaming model to explain the vibration of the tool during the cutting process and the resulting roundness errors in reamed holes, and process parameters such as the direction and the angular frequency of rotation of the tool will affect the dimensional of the roundness error. Deng et al. [25] believed that the dynamic excitation force will cause the deflection and the vibration of the tool, leading the roundness error of the hole, and developed a system control equation composed of exciting force to study mechanism of roundness error of the hole. This will provide great help for this research.
Since the previous research in this area is still blank, the blind cutting experiments are not only difficult to determine the factors that affect the dimensional of R-PDD in tapping, but the cost of the experiment is also very high. Based on previous studies, this paper will establish a quasi-static model of vibration with reference to the analysis method of roundness error of the hole, and finally determine the dimensional of RPDD, laying a solid foundation for the next step of experimental research.

Quasi-static model
During the tapping process, as the tap enters the predrilled hole at an axial rate synchronized with the thread lead, each tooth cuts a layer of material from the surface formed by the previous tooth. In order to ensure the smooth execution of the cutting process, the thread end of the tap is made into a tapered surface with κ r as a chamfer to truncated the full thread, and a number of straight flutes for chip evacuation evenly distributed over the entire thread length divide the thread into several continuously distributed cutting teeth. The major cutting edges of the chamfered section of the tap are formed by the intersection of the tapered surface and the straight flute surfaces, and the minor cutting edges lie on the flank surfaces of the tap threads. As can be seen from the above description, the geometry of the tap is very complicated. In order to convenience the description of the tapping process, some coordinate systems first need to be established. The global coordinate system X-Y-Z whose origin O 1 is at the center of the hole is fixed on the workpiece surface, and the Z axis coincides with the hole axis. The cutting edge coordinate system is U-V-W with the origin O 2 at any designated point of the elementary cutting edge, and the W axis coincides with the major cutting edge. The local rotating coordinate system T-R-A is the transition coordinate system between the global coordinate system and the cutting edge coordinate system. Rotating the global coordinate system around the Z axis through φ i can be converted into the T-R-A coordinate system that changes with the position of the elementary cutting edge, and then rotating this coordinate system around the X and Y axes through λ and κ r respectively to convert to cutting edge coordinates system, as shown in Figures 1a and 2a. The three coordinate systems can be transformed into each other through the transformation matrix. The relationship between them can be expressed as: where As shown in Figure 1b, the dynamic model of the tap only considers the two orthogonal degrees of freedom of the tool in the radial direction. The differential equation of the two-degree-of-freedom transverse vibration of the dynamic tapping system in the global coordinate system is: among them, M, C, K are the mass matrix, damping matrix and stiffness matrix of the tap respectively, and they are all second-order matrices. Whitehead [26] and Towfighian [27] mentioned that the vibrations to cause roundness errors are divided into two categories: chatter and low frequency vibration. In tapping process, the speed of the tap is usually low, and the frequency is much lower than the lowest natural frequency of the tool. At this time, the stiffness term in the equation plays a leading role, while the inertia term and damping term can be ignored. Therefore the above equation becomes the following form: Equation 4 eliminated the inertia and damping terms in the transverse vibration differential equation is called a quasi-static model. The exciting force F includes regenerative cutting force and cutting damping force due to the waviness on the pre-drilled surface. For the convenience of research, let the number of tap flutes be N f , and the tap are evenly dispersed into N z slices with length dz along the axial direction, then the force of the entire tap is distributed to the N f × N z elementary cutting edges. According to the cutting principle of the tap, each elementary cutting edge is regarded as an oblique cutting model.

Cutting force
The cutting force of oblique cutting is proportional to the cross-sectional area of the chip to be removed. There-fore, the cutting force F c can be calculated by the following formula: where k c is the cutting force coefficient, h is the cutting thickness, and b is the cutting width. When the tool's actual trajectory (shown by the dashed line in Figure 2b) deviates from the nominal one due to transverse vibration (shown by the solid line in Figure 2b), an additional dynamic uncut chip thickness will attach to the desired static uncut chip thickness. Since the static uncut chip thickness does not cause vibration, only the load caused by the dynamic uncut chip thickness remains in the cutting force. Therefore, the radial, tangential and axial components of the cutting force of elementary cutting edge in the cutting edge coordinate system U-V-W can be expressed as: where ∆h ij represents the dynamic uncut chip thickness of the j th elementary cutting edge corresponding to the i th flute. If the transverse displacement of the tap axis at time t is represented by x(t) and y(t), the transverse displacement ∆R ij of the elementary cutting edge can be determined by the displacement of the tap axis and the tooth angle φ i (Figure 1c): where τ 0 is the time interval between adjacent teeth. The dynamic uncut chip thickness can be expressed as: As the tapping depth increases, each cutting edge will experience partial and full engagement status. Therefore, the cutting force is constantly changing with the effective cutting edge length involved in cutting during the continuous tapping process. It can be seen from the geometric relationship in Figure 3 that the length of the cutting edge can be determined by its start and end positions in the axial direction. The start and end positions of each cutting edge in the axial direction can be expressed as: where k represents the serial number of the cutting edge, h ks is the start position of the k th cutting edge, h ke is the end position of the k th cutting edge, and P is the pitch. Since the tap cutting edge is not continuously distributed, the window function g(z) is introduced to identify whether the elementary cutting edge corresponding to a certain flute engages with the workpiece at time t. The elementary cutting force can be expressed in the U-V-W coordinate system as: The elementary cutting force in the T-R-A coordinate system can be expressed as follows based on coordinate transformation: k cu = τ s sin ψ n cos(β n − γ n ) + tan λ tan η sin β n cos 2 (ψ n + β n − γ n ) + tan 2 η sin 2 β n k cv = τ s sin ψ n cos λ sin(β n − γ n ) cos 2 (ψ n + β n − γ n ) + tan 2 η sin 2 β n k cw = τ s sin ψ n cos(β n − γ n ) tan λ + tan η sin β n cos 2 (ψ n + β n − γ n ) + tan 2 η sin 2 β n (14) where τ s is the shear stress, ψ n is the normal shear angle, λ is the inclination angle, β n is the normal friction angle, γ n is the normal rake angle, and η is the chip flow angle, and these values can be obtained according to the method in literature [28,29].In the T-R-A coordinate system, the cutting force of all elementary cutting edges corresponding to a certain flute can be expressed as: where l 0 is the start cutting position, l is the current cutting position, they can be expressed as: where D h is the diameter of the pre-drilled hole with a value of 8.5mm, d is the tap nominal diameter, and d 0 is the chamfer point diameter.
By converting the cutting forces of all flutes from the local coordinate system to the global coordinate system and then adding them, the cutting forces ∆F cx , ∆F cy and ∆F cz in the tangential, radial and axial directions of the tap can be calculated as follows: through simplification, Equation 18 can be expressed as: where

Cutting damping force
The actual cutting edge is not sharp, but is a small arc, as shown in Figure 4b.Therefore, the edge force and the cutting damping force are generated because the material under the flank face of the cutting edge is extruded [30].Since the edge force does not cause vibration, only the cutting damping force remains. The cutting damping force is caused by the interference between the uneven workpiece surface caused by the tool vibration and the tool flank surface, which causes the effective flank angle α ef f (as shown in Figure 4a) to change. Literature [31] mentioned that the cutting damping force is composed of the normal force F dv and the friction force F du , and the normal force F dv is proportional to the volume V of the workpiece material extruded under the flank face of the cutting edge: where K sp is the specific contact force, which depends on the material and the geometry of the cutting edge. The volume V is related to the geometry of the cutting edge, the vibration speedċ perpendicular to the machined surface and the cutting speed v c (Figure 4a): where l w is the interference length between the workpiece surface and the flank face of the tool. According to the model proposed by Ahmadi et al. [30], the interference length l w can be written as: where r h is the hone radius, α is the flank angle and θ is the separation angle which defines the position where the workpiece material is no longer removed as chips through the rake face, but is extruded under the flank face of the tool. Therefore, the normal component F dv of the cutting damping force can be written as: where The friction component F du of the cutting damping force is considered to be proportional to the normal component F dv , so it can be written as: where µ is the Coulomb friction coefficient, which is taken as 0.3 according to literature [32]. Therefore, the radial and tangential components of the cutting damping force of the elementary cutting edge in the U-V-W coordinate system can be expressed as: where n is the spindle speed, v c is the cutting speed at different axial positions.
v c = [πd 0 + nP t tan κ r ]n Substituting the dynamic uncut chip thickness into Equation 28, and then based on the coordinate transformation, the elementary cutting damping force in the T-R-A coordinate system can be determined as follows: Like the method of obtaining the cutting force, the tangential ∆F dx , radial ∆F dy and axial ∆F dz components of the cutting damping force can be expressed as: ∆x ∆y (31) through simplification, Equation 31 can be expressed as:

Elastic force
In the quasi-static equation, KX is the elastic force, where K is the stiffness matrix of the tool. For a twodegree-of-freedom system, it can be expressed as: Due to the small bending deformation of the tool, the vibration of the tool can be considered to be restricted to a plane perpendicular to the axis of the undeflected tool. Therefore, the relationship between the elastic force and the axis offset can be expressed as: It is assumed that the tap is a cantilever beam whose force is concentrated on the tip of the tool in this analysis. According to literature [27],the stiffness matrix can be written as: where EI is the bending stiffness of the tool material,E is the elastic modulus of material of the tap,and L is the total length of the tap.

Solutions of quasi-static equations
Now substituting Equations 19 and 32 into Equation  4, the quasi-static equilibrium equation is rewritten as follows according to the displacement of the tap axis: The current value of X is expressed as a function of the previous value, and the friction time delay τ is approximately a fraction of the cutting time delay τ 0 , that is, τ = τ 0 /m. Therefore, the successive tool axis position equation can be expressed as: where Equation 38 is a difference equation with variable coefficients. In order to obtain the solution of the equation, based on the freezing coefficient method [33], the successive cutting process is discretized into several steps with time τ , each of which is approximated by a constant coefficient difference equation. For the convenience of solving, Equation 38 is written as a large matrix form: . . .
this is where the matrix Q t−τ is the state-transition matrix at time t−τ ,I is a second-order identity matrix. According to the literature [34], Equation 41 can be written as follows: where λ is a characteristic exponent, whose imaginary and real parts represent the oscillation frequency and growth or decay rate respectively. Equation 42 is a characteristic equation with m eigenvalues. The eigenvalues and eigenvectors satisfying the matrix Q t−τ are µ t−τ and q t−τ respectively, which can be determined by solving the characteristic equation. Therefore, the instantaneous position vector X i t−τ of m modes at time t − τ can be extracted from the m eigenvectors respectively. The displacement vector X i of the i th mode at time t can be written as: At time t, the displacement vector X of the tool axis can be obtained by superimposing the displacement vectors of m modes. The pitch diameter of the thread will also change with the axis offset. The pitch diameter at time t can be expressed as: where D 2 (t) is the pitch diameter of the thread at time t, and D is the nominal diameter of the thread which is equal to the nominal diameter d of the tap,the offset ∆D(t) at time t can be expressed as: where φ(t) is the angle of rotation in time t.  RPDD of the thread can be expressed as the difference between the maximum and minimum pitch diameters: where D 2max and D 2min are the maximum and minimum pitch diameters of the thread, respectively.

Time domain simulation is performed based on the tool motion equation of Equation 38
, which is also quasistatic. The static equilibrium position of the tool, sub- ject to regenerative cutting and cutting damping forces, is found at every time step in the simulation, and then the dimensional of RPDD is determined by the Equation 47.The tool used in the simulation is a standard straight flute high-speed steel tap, and its parameters are shown in Table 1.The workpiece is an AISI1045 plate with pre-drilled holes. In the cutting simulation, 6 groups of tests whose parameters are listed in Table  2 were performed to study the influence of the chamfer length and the spindle speed on RPDD.  In order to analyze the influence of the chamfer length and the spindle speed on RPDD, the eigenvalues and eigenvectors of all modes of each group of tests are obtained by simulation. The vibrations of first four modes which are selected from all modes according to the rate of decay or growth of each mode because of their importance are considered here, and RPDD obtained is a combination of these four modes. Figure 5 shows the change of the pitch diameter of the first full thread and the movement trajectory of the outermost cutting edge of the tap in the global coordinate system during the entire cutting process when the chamfer length is 2 threads, 4 threads and 8 threads. From the comparison of the three figures, it is found that the change of the chamfer length has an effect on RPDD.As can be seen in Figure 7a, the dimensional of RPDD shows a trend of first increasing and then decreasing with the increase of the chamfer length and the smallest RPDD is obtained when the chamfer length is 2 threads. The possible reason is that fewer cutting edges are involved in cutting, causing the vibration amplitude to increase briefly and then become a steady periodic vibration. Figure 6 shows the change of the pitch diameter of the first full thread and the movement trajectory of the outermost cutting edge of the tap in the global coordinate system during the entire cutting process when the spindle speed is 700 rev/min, 1050 rev/min, 1400 rev/min and 1750 rev/min. From the comparison of the three figures, it is found that the change of the spindle speed has an effect on RPDD.As can be seen in Figure 7b, the dimensional of RPDD shows a trend of first decreasing and then increasing with the increase of the spindle speed and the smallest RPDD is obtained when the chamfer length is 1400 rev/min. The possible reason is that more cutting edges quickly participate in cutting when the spindle speed increases, resulting in a rapid increase in damping and a decrease in amplitude. But at the same time, the increase of the spindle speed will also lead to the decrease of the damping change rate, so a higher speed will restrain the amplitude reduction.

Summary
A quasi-static tapping model, including delayed cutting force and cutting damping force was established. A similar cutting force model and cutting damping force were developed based on the work of previous researchers. The force equilibrium equation is expressed as a discretetime matrix equation describing the successive state of the tool. According to the eigenvalues and eigenvectors of the state transition matrix, the characteristic solution of the equilibrium equation is found, and then R- This model was used to simulate the tapping of AISI1045 plates using standard straight flute high-speed steel taps, and the dimensional of RPDD under different chamfer lengths and spindle speeds was obtained. By comparison, it is found that the chamfer length and spindle speed will affect RPDD, and RPDD is the smallest when the chamfer length is 2 threads and the spindle speed is 1400 rev/min. Since the research on R-PDD is still blank so far, the development of this model not only provides a cheap and effective method for the study of RPDD, but also lays a foundation for further experimental research.
This model provides a method for the research of RPDD and the research on the stability of tapping in the literature [23] also confirms the conclusion of this article, but it needs to be verified by experiments in order to better illustrate its effectiveness. In addition, in the actual tapping process, the distributed force acts on the entire engagement length of the tool, but this model applies concentrated force on the tip of the tool, which will have a certain impact on the accuracy of the model. These issues will be studied in the future.

Funding information
The research is financially supported by the National Natural Science Foundation of China (No. 51275333)

Conflict of interest
The authors have no relevant financial or non-financial interests to disclose.

Availability of data and material
The data sets supporting the results of this article are included within the article and its additional files.

Authors' contributions
Jie Ren developed a mathematical model to predict the radial diameter difference of threads, analyzed the simulation results, and was a major contributor in writing the manuscript. Xianguo Yan provided guidance for the writing of manuscript. All the authors read and approved the final manuscript.

Ethics approval
All analyses in this paper are based on previously published research and this paper does not involve animal and human testing, so this item is not applicable to this paper.

Consent to participate
All analyses in this paper are based on previously published research and this paper does not involve animal and human testing, so this item is not applicable to this paper.

Consent for publication
The Author confirms: that the work described has not been published before (except in the form of an abstract or as part of a published lecture, review or thesis); that it is not under consideration for publication elsewhere; that its publication has been approved by all co-authors ,if any; that its publication has been approved (tacitly or explicitly) by the responsible authorities at the institution where the work is carried out. The author agrees to publication in the Journal indicated below and also to publication of the article in English by Springer in Springer's corresponding Englishlanguage journal. The copyright of the English article is transferred to Springer effective if and when the article is accepted for publication.   Tool path and cutting damping force geometry Change of the pitch diameter of the first full thread and the movement trajectory of the outermost cutting edge of the tap in the global coordinate system during the entire cutting process when the chamfer length is 2 threads, 4 threads and 8 threads Change of the pitch diameter of the first full thread and the movement trajectory of the outermost cutting edge of the tap in the global coordinate system during the entire cutting process when the spindle speed is 700rev/min, 1050rev/min, 1400rev/min and 1750rev/min