Application of the FEM Method to Modeling and Analysis of the Cold Thread Rolling Process—Part 1: Conditions for Ensuring a Plane State of Deformations

The article concerns the application of the FEM method for the prediction of stress and deformation states in a workpiece during the thread rolling process (TR). The analysis covered a new kinematic variant of the TR process in which the basket of the head rotates and is torque-driven, while the workpiece is stationary and the head with the rollers moves axially relative to the workpiece. The TR process was considered as a geometrical and physical non-linear initial and boundary problem with non-linear, moving, and variable in time and space boundary conditions. The boundary conditions in the contact areas of the tool with the workpiece were unknown. An updated Lagrange (UL) description was used to describe the physical phenomena at a typical incremental step. The states of strain and strain rate were described by non-linear relationships without linearization. New discrete systems of motion and deformation equations of the object in the TR were introduced, which take into account the change in the stiffness of the system during the TR process. This equation was solved by the central differences method (explicit). The material parameters were estimated during tensile tests to determine the characteristics of the C45 steel, and a new semi-empirical method was used to determine the relationship yield stress, effective true strain, and effective true strain rate in the thread rolling process. A modified Cowper–Symonds material model was also used to model the displacement process of the wedge on an elastic/visco-plastic body reflecting the TR process. A non-linear dependency of material hardening module depending on strain and strain rate was introduced. To confirm the plane state of deformation and spatial state of stress, an experimental investigation was carried out. The computer models were validated, and a good convergence of the results was obtained. Applications in the ANSYS/LS-Dyna program were developed to simulate the TR process. The developed applications enable a comprehensive time analysis of the states of displacement, strain, and stress occurring in an object consisting of a workpiece (shaft) and a tool (roller) for the case of a plane strain state and a spatial stress state. Exemplary results of numerical analyzes are presented to explain the influence of the friction coefficient on the condition of the thread quality, and the state of deformations and stresses were shown.


Introduction
The first attempts to model the process of thread rolling on a two-roller thread rolling machine with a steady rest were carried out in [1]. Formulas were developed for the required output diameter of the blank for various types of external threads (metric, trapezoidal) and the recommended bevel angles for the shaft, assuming that the volume of the material remains unchanged before and after rolling. The permissible deviations of the semi-finished product diameter were determined, depending on the required tolerance of the thread execution. Recommendations and nomograms for rolling parameters were also developed, including roll pressure force, rolling speed, and roll feed speed, as well as a selection of cooling and lubricating liquids. The formulas for forces, moments, and 2 of 40 rolling work were given, as well as the allowable deviations of the roll profile half-angle and the allowable deviations of the roll profile pitch and tolerances of the outer diameter depending on the type of thread and its dimensions.
Olszak W. et al.'s team developed a new method of designing angular head rolls for the rolling of external threads [2]. The rollers have rolling rings with gradually increasing diameters. The rational design of these tools consists, among other things, of obtaining the desired or acceptable load conditions of individual rings, as well as of entire rollers. The dividing methods of the machining allowance between the rolling rings were developed, and the geometrical parameters of the rolls and formulas for the components of the force during the rolling of deep profiles were determined, taking into account the elastic and plastic deformation of the part material: where: F y , F z -components (radial and circumferential) of the rolling force acting on the rolling ring, a-width of the layer pressed by the ring, h-hollow of the ring in the material, b 1 , b 2 -empirical coefficients depending on other rolling conditions. The relationship structure (1a and 1b) allowed for the determination of the forces acting on the finishing rings (when a = 0). The investigations were carried out using three rollers, where each of them had either one, two, or three rings [3]. Those variations allowed for different operating conditions of the outer and inner roller rings. The force acting on the outer ring is defined as half of the force acting on the roller with two rings. The inner ring was loaded with a force equal to the difference of forces acting on the rollers with three and two rings or with two and one ring. By summing up the forces acting on the active rings of the roller, the total load on the roller was obtained.
Ivanov V. [4] developed equations describing the required axial profile of rolls rolling long threads by the axial method [5,6] and gave empirical relationships for the rolling forces and moments [7,8], as well as equations of the contact surface of the part material with the active surface of the rolls in different phases of the process. Łyczko K. [9] developed an experimental and theoretical methodology for determining the curvilinear axial profile of a tap forming a thread in (internal) holes based on the analysis of allowances measured by the volume of the material displaced through an individual root. Łyczko derived formulas for the volume of the material deformed by an individual root depending on the deformation depth. The obtained outlines were approximated by a circle, a parabola, a logarithmic curve, or by the fourth degree polynomial. Additionally, models were derived in the form of a fourth-degree polynomial for the increase in the microhardness of the thread, being formed depending on the number of roots. Łyczko derived the formulas for the recommended diameters and phases of the thread hole depending on the design and dimensions of the forming screw-tap and the type and size of the thread.
Since 2002, the application of computer modeling to the thread rolling process has been performed using the finite element method (FEM). Domblesky J. and Feng F. developed computer applications in the DEFORM-PRO v. 3.0 [10] and DEFORM-3D v. 2.3 [11] systems to simulate the rolling process of metric threads. However, the applications did not take into account the sensitivity of the material to the rate of strain and the history of the process. It was also not possible to obtain results for the complete immersion of the rolling rolls in the workpiece. This was due to the excessively large amount of non-linearities occurring during rolling and the instability of the numerical calculation process. Kinematic diagrams (a-c) of three methods of rolling external threads and the rolling heads according to these methods (d-f): (a)-stationary basket of the head and rollers driven and a rotating non-driven workpiece that moves axially relative to the head; (b)-stationary head basket, rotating and non-driven rollers and rotating and driven workpiece; the head with the rollers moves axially relative to the workpiece; (c)-stationary object, rotating and driven basket of the head and rotating and non-powered rollers; the head with the rollers moves axially relative to the workpiece: 1-head basket, 2-rollers, 3-workpiece, v -rolling speed, v -peripheral speed of the head.
The first kinematic variant is that the head basket (1) is stationary (n 0), while the rollers (2) are driven by the torque M (n 0), which sets the workpiece (3) in rotation (n 0) [1,2,4,10,11,14,15,16,17] (Figure 1a). In addition, the workpiece (3) during the machining process moves axially relative to the head with the rollers. This method is used for the rolling of threads of any length and is also called the through centerless rolling method. The force F 0 acting on the object causes the material to move in the tangential direction, causing spatial states of displacements, strains, and stresses. The second method is that the head basket (1) is stationary (n 0), while the workpiece is driven by the torque M (n 0), which sets the rollers in rotation (n 0) [2,15,18] (Figure 1b). In addition, the head (1) with the rollers (2) moves axially relative to the workpiece (3) during the machining process. This method is used for rolling threads of limited length. The force F 0 acting on the rollers causes displacement of the object material in the tangential direction, causing a spatial state of displacements and strains.
In order to create plane states of displacements and strain in the workpiece, as well as a spatial state of stresses, the Author used an innovative method of thread rolling ( Figure 1c) [15,16,17,19]. This method is based on the fact that the head basket (1) is driven Figure 1. Kinematic diagrams (a-c) of three methods of rolling external threads and the rolling heads according to these methods (d-f): (a)-stationary basket of the head and rollers driven and a rotating non-driven workpiece that moves axially relative to the head; (b)-stationary head basket, rotating and non-driven rollers and rotating and driven workpiece; the head with the rollers moves axially relative to the workpiece; (c)-stationary object, rotating and driven basket of the head and rotating and non-powered rollers; the head with the rollers moves axially relative to the workpiece: 1-head basket, 2-rollers, 3-workpiece, v r -rolling speed, v h -peripheral speed of the head.
The first kinematic variant is that the head basket (1) is stationary (n h = 0), while the rollers (2) are driven by the torque M r (n r > 0), which sets the workpiece (3) in rotation (n w > 0) [1,2,4,10,11,[14][15][16][17] (Figure 1a). In addition, the workpiece (3) during the machining process moves axially relative to the head with the rollers. This method is used for the rolling of threads of any length and is also called the through centerless rolling method. The force F w > 0 acting on the object causes the material to move in the tangential direction, causing spatial states of displacements, strains, and stresses.
The second method is that the head basket (1) is stationary (n h = 0), while the workpiece is driven by the torque M w (n w > 0), which sets the rollers in rotation (n r > 0) [2,15,18] (Figure 1b). In addition, the head (1) with the rollers (2) moves axially relative to the workpiece (3) during the machining process. This method is used for rolling threads of limited length. The force F r > 0 acting on the rollers causes displacement of the object material in the tangential direction, causing a spatial state of displacements and strains.
In order to create plane states of displacements and strain in the workpiece, as well as a spatial state of stresses, the Author used an innovative method of thread rolling (Figure 1c) [15][16][17]19]. This method is based on the fact that the head basket (1) is driven by the torque M h and rotates at a speed of n h > 0, while the object (3) is stationary (n w = 0). The rollers (2) are not driven but rotate at speed n r > 0. In addition, the head (1) with Materials 2023, 16, 4647 4 of 40 the rollers (2) during the machining process moves axially relative to the workpiece (3). A characteristic feature of this kinematic variant is that tangential forces do not act on the workpiece (F w = 0 and F r = 0), but only normal forces in planes parallel to the axial plane act. Thus, with such a system of forces, there is no displacement of the material in the tangential direction. This method is used to rolling threads with increased accuracy but limited lengths.
All three kinematic variants of thread rolling are implemented in the Center of New Technologies at the Faculty of Mechanical Engineering of the Koszalin University of Technology. The rolling station according to the first method is shown in the Figure 2a In addition, in the third method, in order to ensure plane states of displacements and strains, it is necessary to maintain the tip clearance between the rollers and the shaped thread outline shown in Figure 3. Ensuring plane states of displacements and strain of the material during thread rolling makes it possible to shape the thread outline with greater dimensional accuracy and to obtain a thread of unprecedented quality.   Other factors that adversely affect the quality of the thread are caused by insufficient accuracy of modeling of the rolling process. Thus far, the rolling process has been described by equations of motion that do not take into account the dynamics of the process or do not take into account the variable stiffness and variable damping of the system during TR process. Therefore, it is assumed that the stiffness matrix and the damping matrix are constant (K = const and C = const) [10,11,14,20].
In addition, the mechanical properties of materials based on the uniaxial tensile test of cylindrical samples on a testing machine are described in accordance with the standard on a global scale without taking into account the deformation states in the neck [27], which means that the values of true plastic deformations are underestimated by up to three times. This is due to the fact that it is impossible to determine the elongation value of the sample at which a separation crack occurs inside the neck, and it is impossible to determine the actual values of spatial displacements of the material inside the neck, which are then used to calculate the true values of material strain.
The adoption of inappropriate measures of physical quantities and inappropriate methodologies for calculating deformation and stress states causes that the mechanical properties of materials are described inadequately to the requirements. The semi-empirical method proposed by the Author does not have these disadvantages. For example, for C45 steel, the true limit deformation (the so-called failure strain) calculated according to the standard is ϕ f = 0.51 [-](ε f = 0.66 [-]), while the true failure strain calculated taking into account the true ones in the neck, according to the proposed methodology, Geometric and physical non-linearity of the system is also omitted by using inadequate material models. The use of Johnson-Cook [30] or Cowper-Symonds [14] constitutive equations means that the variability of the material-hardening modulus during machining is not taken into account and it is treated as a body with linear hardening, the so-called bilinear hardening material model. In addition, the invariability of the plastic strain rate of the material is assumed, assuming a constant value of the effective strain rate in the rolling process. Meanwhile, in the thread rolling process, the maximum effective strain rate occurs in the initial phase of contact of the rollers with the workpiece and amounts to . ϕ e = 10 4 s −1 , and then decreases to . ϕ e = 0 s −1 in the axial plane. The real impact of the strain rate on the hardening of the material is large and increases the yield stresses by up to 20%. However, in the assumed material models, the efficiency of the strain rate is logarithmized, so, for example, for . ϕ e = 1 s −1 , ln( . ϕ e ) = 0, while for . ϕ e = 10 4 s −1 , ln( . ϕ e ) = 9.21, and for . ϕ e = 10 5 s −1 , ln( . ϕ e ) = 11.51. As a result, the calculated sensitivity of the material to strain rate is underestimated. It is also not possible to calculate the function K . ϕ , which depends on the effective strain rate, for . ϕ e = 0 s −1 . In our own works [15][16][17][18][19]28,29], the incremental updated Lagrange description, variational formulation, dynamic theory of metal materials with velocity-type memory, contact mechanics, non-linear theory of continuous medium, and the Finite Element Method are used.
Our own research shows [15][16][17][18][19]28,29] that in typical thread rolling tools, when the profile is completely filled, the deformed material flows not only in planes parallel to the axial plane, but also in the tangential direction (perpendicular to the axial plane), and then there are spatial states of displacements, strains, and stresses ( Figure 3).
In accordance with our work, it is shown by the outline of the rolling of the thread roll that when increasing the height of the working ring, the tip clearance ensures that the top of the thread profile does not make contact with the roller. Then, with some approximation, it can be assumed that plane states of displacements and strains, as well as a spatial state of the stresses, appear in the deformed material, and there is no unfavorable jamming of the shaft material during rolling. In order to test the hypothesis, model tests, numerical simulations using the FEM method, and experimental verification were carried out ( Figure 4).
axial plane, but also in the tangential direction (perpendicular to the axial plane), and the there are spatial states of displacements, strains, and stresses ( Figure 3).
In accordance with our work, it is shown by the outline of the rolling of the threa roll that when increasing the height of the working ring, the tip clearance ensures that th top of the thread profile does not make contact with the roller. Then, with some approx mation, it can be assumed that plane states of displacements and strains, as well as a sp tial state of the stresses, appear in the deformed material, and there is no unfavorab jamming of the shaft material during rolling. In order to test the hypothesis, model test numerical simulations using the FEM method, and experimental verification were carrie out ( Figure 4). Determining the states of displacements and deformations of the material in th thread-rolling process is a fundamental issue because it determines the possible simplif cations not only of the numerical model but also of the model used for model investig tions. In addition, reducing the complexity of the model reduces the time-consuming ca culations and the demand for the computing power of computers, which is of significa practical importance when designing the rolling process. Knowledge of the mechanism plastic material flow during thread rolling, and mainly of the states of displacement of th workpiece in various phases of the process and after unloading-the so-called sprin back-are necessary for the correct design of the outline of the active surface of the knead ing rolls.
The article deals with the determination of thread rolling processing conditions ensure flat states of displacement and deformation of the workpiece material. To predi stress and strain states and the displacement and spring back of the workpiece materi during and after the thread rolling process, the FEM method was used. The updated L grange description was used to describe the physical phenomena at a typical increment step. Non-linear incremental models without any linearization were used to describ physical phenomena. The results of the numerical simulation were experimentally ver fied on the example of selected parameters of the thread. Figure 4 shows an exemplary roller for rolling of a trapezoidal thread Tr22×4, who geometry of the active (working) rings was modified so that during the rolling proces  Determining the states of displacements and deformations of the material in the threadrolling process is a fundamental issue because it determines the possible simplifications not only of the numerical model but also of the model used for model investigations. In addition, reducing the complexity of the model reduces the time-consuming calculations and the demand for the computing power of computers, which is of significant practical importance when designing the rolling process. Knowledge of the mechanism of plastic material flow during thread rolling, and mainly of the states of displacement of the workpiece in various phases of the process and after unloading-the so-called spring back-are necessary for the correct design of the outline of the active surface of the kneading rolls.
The article deals with the determination of thread rolling processing conditions to ensure flat states of displacement and deformation of the workpiece material. To predict stress and strain states and the displacement and spring back of the workpiece material during and after the thread rolling process, the FEM method was used. The updated Lagrange description was used to describe the physical phenomena at a typical incremental step. Non-linear incremental models without any linearization were used to describe physical phenomena. The results of the numerical simulation were experimentally verified on the example of selected parameters of the thread. Figure 4 shows an exemplary roller for rolling of a trapezoidal thread Tr22×4, whose geometry of the active (working) rings was modified so that during the rolling process, plane states of material displacement and deformation occur in planes parallel to the axial plane of the formed thread.

Modeling of Thread Rolling Process
The thread rolling process was considered as a geometrically and physically non-linear boundary-initial problem, in which there are non-linear, moving and variable in time and Mathematical description of non-linear phenomena requires the use of other than in linear problems, principles of formulating boundary-initial problems and more complex methods of solving them.
Geometric non-linearity results from the non-linear relationship between strain and displacement. It is caused by a change in the geometry of the workpiece during the machining process.
Physical (material) non-linearity is caused by non-linear mechanical properties of the workpiece material, leading to a non-linear stress-strain relationship.
For the modeling and analysis of the thread rolling process as a real object, which consists of a three-rollers angle head and a workpiece (see Section 2.5.6), a modern algorithm was used, with variational method and FEMs ( Figure 5) [20,[31][32][33][34][35].
plane states of material displacement and deformation occur in planes parallel to the axial 230 plane of the formed thread. Physical (material) non-linearity is caused by non-linear mechanical properties of the 244 workpiece material, leading to a non-linear stress-strain relationship. 245 For the modeling and analysis of the thread rolling process as a real object, which 246 consists of a three-rollers angle head and a workpiece (see section 2.5.6.), a modern algo-247 rithm was used, with variational method and FEMs ( Figure 5) [20,[31][32][33][34][35]. 248 Modelling of the process Real object

PHYSICAL MODEL OF CONTINUOUS OBJECT
Mathematical modelling FEM discretization

The Concept of Incremental Description
Due to the strong geometrical and physical nonlinearities occurring in the thread rolling process, the updated Lagrange description [31][32][33][34][35] (Figure 6), known in non-linear mechanics, was used. Using this description, for a typical time step τ = t + ∆t, where ∆t is a small increment of time, calculated on the basis of the natural frequency of the system, it is assumed that solutions for time steps in the interval [0, t] are known, and the solutions for the next moments are sought. As a result, for successive discrete moments of time τ = 0, ∆t, 2∆t, . . . the geometry of the body and the states of incremental displacements, displacement velocities, accelerations, strains, strain rates, and stresses are determined. rolling process, the updated Lagrange description [31][32][33][34][35] (Figure 6), known in non-linear 253 mechanics, was used. Using this description, for a typical time step τ t Δt, where Δt 254 is a small increment of time, calculated on the basis of the natural frequency of the system, 255 it is assumed that solutions for time steps in the interval [0, t] are known, and the solutions 256 for the next moments are sought. As a result, for successive discrete moments of time τ 257 0, Δt, 2Δt, … the geometry of the body and the states of incremental displacements, dis-258 placement velocities, accelerations, strains, strain rates, and stresses are determined. The concept of incremental description has already been verified by the author and 262 was used to describe other plastic-forming processes: cutting [36][37][38][39][40][41][42], thread rolling [15-263 17,28,29], burnishing rolling [43][44][45][46][47][48][49], drawing [50][51][52], grinding with a single abrasive 264 grain [53], calculating the fatigue strength of products [54][55][56], and diamond sliding bur-265 nishing [56]. 266 It is assumed that at the moment t, the initial (original) configuration of the body C 267 and the current configuration C are known. On the other hand, the dynamic equilibrium 268 configuration C at the next moment τ = t + Δt is sought ( Figure 6). 269 The application of the updated Lagrange description also requires the adoption of The concept of incremental description has already been verified by the author and was used to describe other plastic-forming processes: cutting [36][37][38][39][40][41][42], thread rolling [15][16][17]28,29], burnishing rolling [43][44][45][46][47][48][49], drawing [50][51][52], grinding with a single abrasive grain [53], calculating the fatigue strength of products [54][55][56], and diamond sliding burnishing [56].
It is assumed that at the moment t, the initial (original) configuration of the body 0 C and the current configuration t C are known. On the other hand, the dynamic equilibrium configuration t+∆t C at the next moment τ = t + ∆t is sought ( Figure 6).
The application of the updated Lagrange description also requires the adoption of appropriate coordinate systems, the definition of adequate measures of increments of displacements, strains, and stresses, as well as procedures for their accumulation at the incremental step. Three systems of orthogonal Cartesian coordinates have been adopted in this paper: {x} movable object-related, {y} movable tool-related, and {z} fixed reference system ( Figure 6).

Measures of Increments of Physical Quantities
The tool was treated as a rigid or elastic body, and the object as a body in which elastic deformations (in the range of the reversible zone) and viscous and plastic (in the range of the irreversible zone) may occur, assuming mixed (kinematic and isotropic) non-linear strain hardening. For the build a material model, abbreviated as E/VP, the non-linear Huber-Mises-Hencky (HMH) plasticity condition, plastic potential, and flow law associated with the yield surface was used. The initial states of displacements, stresses, and strains after previous operations were also introduced.
The strain state of the workpiece is described by the symmetric Green-Lagrange strain tensor t t T ε , while the strain rate state is described by the Green-Lagrange strain rate tensor In a similar manner, the following formula for the increments of the components τ t ∆ . ε ij of the strain rate incremental tensor τ t ∆T . ε are obtained [49,50]: The introduction of non-linear terms means that exact formulas will be used, without any linearization.
The stress state and stress increment were described by the second (symmetric) Pioli-Kirchhoff stress tensor t t T σ and its increment τ t ∆T σ , respectively. The components ∆σ ij of the stress increase tensor τ t ∆T σ for the body E/VP were calculated from the following formula: where ψ is the load factor and is ψ = 1 for loading and neutral and ψ = 0 for unloading processes, A kl are the tensor components and B is a scalar (hardening function), and for mixed hardening were calculated according to the following Equations: In the case of kinematic hardening, formula (16) takes the following form: while for isotropic hardening-the following form is taken: In Equations (4)- (8), ∆ϕ kl are the components of the increment of the tensor of total true strain (elastic, plastic, and viscous), C ijkl is the fourth-order tensor characterizing the elastic properties of the object material, ) takes the following form: where K ϕ and K . ϕ are one-parameter functions of yield stresses depending on strains and strain rates, respectively, with the boundary condition that for . ϕ (VP) e = 0, K . ϕ = 1, must be met. The functions K ϕ and K . ϕ , in the general case, have the following forms: The material constants a − h, n 1 − n 4 , K and ϕ 0 in the yield stress models (9), (10), and (11) were determined in experimental studies, and for C45 steel, are listed in Section 2.5.4.
The increment of the yield stress at a typical step t → t + ∆t is defined by the following equation [44]: where ∆ϕ from the data obtained in a series of tensile tests at different true strain rates using virgin material specimens (see Section 2.5.1). These partial derivatives are calculated using a geometric interpretation of the derivative [44]. Then, the partial derivative at a given time t is equal to the tangent of the angle of inclination of the tangent to the corresponding hardening curve at time t: where K φ and K φ are one-parameter functions of yield stresses depending on strains 314 and strain rates, respectively, with the boundary condition that for φ 0, K φ 1, 315 must be met. The functions K φ and K φ , in the general case, have the following forms: The material constants a h, n n , K and φ in the yield stress models (9), (10), 317 and (11) were determined in experimental studies, and for C45 steel, are listed in section 318 2.5.4.

319
The increment of the yield stress at a typical step t → t ∆t is defined by the follow-320 ing equation [44]: where Δφ and Δφ are the increments of the visco-plastic effective true strain and 322 the effective true strain rate, respectively. 323 In order to further evaluate the above expression for Δσ φ , φ from (12), it 324 is necessary to obtain partial derivatives ∂σ φ , φ / ∂φ and ∂σ φ , φ / 325 ∂φ from the data obtained in a series of tensile tests at different true strain rates using 326 virgin material specimens (see section 2.5.1.). These partial derivatives are calculated us-327 ing a geometric interpretation of the derivative [44]. Then, the partial derivative at a given

The Equation of the Motion of a Discrete Object on a Typical Step Time
The variational, non-linear equation of motion and deformation of the object was derived using the stationarity condition of the functional, which is the incremental of the total energy of the system, at a typical step time. This equation has an entangled form, and in order to solve it, a spatial approximation appropriate to the finite element method was used, obtaining discrete systems of an equation of motion and the deformation of a typical finite element and for a typical time step t → τ = t + ∆t in the thread rolling process. By writing the specified equations of motion for all finite elements separated from the tool and workpiece after their transformation to a common global coordinate system {z} and after expanding and summing over all finite elements, the equation of motion of a discrete object in the updated Lagrangian description in a typical step incremental was obtained with the following form [44]: where: M-temporary mass matrix; C-temporary damping matrix; K, ∆K-respectively, temporary stiffness matrix and its increment; F, ∆F-respectively, temporary column vector of the internal loads of nodes and its increment; R, ∆R-respectively, temporary column vector surface loads and its increment; ∆r-column vector of node displacement increments; ∆ r-column vector of node acceleration increment.
The System of Equation (13) contains N equations differential of the second order with constant (incremental step) coefficients in which 2N known elements of the vector of internal forces F and external forces R and 3N 2 elements of the matrix M, C and K occur, but it contain s4N + N 2 unknowns, i.e., the components of the vectors: an increment of displacements of nodes ∆r, an increment of velocity of nodes ∆ . r, an increment of node accelerations ∆ .. r, an increment of the internal loads of the object ∆F, and N 2 unknown elements of the object stiffness increment matrix ∆K. In this equation, it is also a part of the components of the vector of increase of external loads ∆R in contact areas that is unknown.
Equation (13) was solved using partial linearization and the method of central differences (explicit), in which the differential approximation of the derivatives of partial displacements is assumed according to the following Equations [32]: ..
where a 1 = 1 2·∆t and a 2 = 1 ∆t 2 are the integration constants, then Equation (13) takes the following form: where ∼ M is the effective mass matrix and ∼ Q the effective load vector: From Equation (15), the unknown column vector of the displacements of the nodes of the discrete object at the end of the considered time step r t+∆t is calculated for given initial r 0 and boundary values for the case of the thread rolling process.
In the already formulated constitutive relationships (material models) Equations (2)-(8), the Cauchy index notation was used to facilitate mathematical operations. In order to use these relationships for the Object Motion Equation (16) derived and to develop computer applications for a numerical simulation of the thread rolling process, it is convenient to write them in a matrix form.
Relationship (4) now takes a matrix form. The increment of the stress tensor ∆T σ , for mixed hardening, is calculated from the following formula: where C (E) is the Hook matrix, ∆ϕ is the column vector of Green-Lagrange of total true strain increment, A is the matrix, and B is the scalar, according to the formulas: The increment of the total strain tensor ∆T ε is calculated from the following formula: where is the scalar, ∆ϕ * * are is the column vector of Green-Lagrange of total true strain, and G is a matrix, according to the following formula: In numerical calculations, the matrix notation of Formula (18) obtained after its discretization by the finite element method is also used: where − B and ∼ B are the temporary matrices of the linear and non-linear dependence of the strain increment ∆ε on the displacement increment of the nodal points ∆r, respectively; is the stress matrix of the discrete system; [∆q] and {∆q} are, respectively, a matrix and a column vector of the displacement increment of the nodal points of the system on a step time.
In the calculation process, it is necessary to accumulate incremental quantities, while the problem of accumulation of the components of the vectors of the displacement increase ∆u i , and the components of the true strain increment tensors ∆ϕ ij and the components of the true strain rate increment tensors ∆ . ϕ ij at any time τ is trivial, and it consists in adding their respective components derived in the previous steps.
The accumulation of the effective yield stress increments ∆σ Y and the increments of stress tensors ∆T σ and their components ∆σ ij is not trivial and it requires the use of appropriate algorithms named "Yield stress accumulation" and "Stress accumulation".

Samples Tested and Chemical Compositions
The tensile strength tests were carried out at an ambient temperature of 18 • C for 5 cylindrical samples made of drawn steel bars in accordance with the PN-EN 10002-1:2004 standard [58] (Figures 8a and 9a) and for 5 three-stage cylindrical samples (Figures 8b and 9b). The material parameters of C45 steel (1.0503) used in numerical calculations were determined in a uniaxial tensile test, on a Zwick Roell Z400E testing machine. The deformation of the samples was measured using the advanced optical measuring system ARAMIS Ad-jusTable 12M, which uses the technique of digital image correlation, coupled with testing machine and HBM MGCplus data acquisition system. The system allows non-contact measurement of displacements and strains of test samples of any shape and size from 20 × 15 × 4 mm to 2500 × 1800 × 480 mm. The measurement results are obtained in real time in the full research cycle and the entire field of observation. The data-recording rate is from 25 fr/s to 100 fr/s, the resolution is 12 million pixels (pixel size 3.45 µm). The system is at the Rapid Prototyping Center in Technical University of Koszalin ( Figure 10). The influence of strain and strain rate on the yield stress obtained from experiment was compared with the basis of catalog data [59]. The chemical composition was tested in the BATORY Research Laboratory in Chorzów, Poland, by optical emission spectrometry using an optical emission spectrometer type ARL 3460 with spark excitation, manufactured by Thermo Electron in accordance with the PN-H-04045:1997 standard [60]. The tests were carried out at an ambient temperature of   In the plastic zone-deformation uniform ( Figure 11-section AB), i.e., flow curves for materials insensitive to changing the strain rate, the relationship between stress and strain were described by the following non-linear empirical models: Hollomon's J.H. [61], Ludwig's P. [21], and Swift's H.W. [22].  In the plastic zone-deformation uniform ( Figure 11-section AB), i.e., flow curves 441 for materials insensitive to changing the strain rate, the relationship between stress and 442 strain were described by the following non-linear empirical models: Hollomon's J.H. [61], 443 Ludwig's P. [21], and Swift's H.W. [22].  Figure 11. Stress-strain engineering approximate, and true (hardening) curves: ε is the effective engineering plastic strain, ϕ is the effective true plastic strain, σ Y0 is the initial yield stress.
A fundamental problem is to determine the state of stress and strain of the material beyond the tensile strength at the time of loss of stability of the specimen (i.e., at the time of formation of the neck- Figure 11, point B). The formation of the characteristic throat uniform interferes with uniaxial stress state. In the smallest cross-section of the specimen, radial and circumferential stresses also appear. Thus, there is a three-dimensional stress state. To determine the state of true stress in the neck ( Figure 11-section BD), the most known methods used so far were: Bridgman P.W. [23], Davidenkov N.N. and Spiridonova N.E. [24], Siebel E. [25], and Szczepiński W. [26]. The effective true stress σ e in the specimen at any time of the tensile test in the section BD and the local stress σ z in the z-axis in the smallest cross-section of the specimen, depending on the distance r from the center, are calculated using the Bridgman, Davidenkov-Spiridonova, Siebel, or Szczepiński methods [62]. The use of material parameters calculated with these methods for numerical calculations of physical phenomena occurring in the thread rolling process led to a significant overestimation of force parameters (see Section 2.5.3) and, consequently, to product quality inconsistencies. This can also lead to the destruction of the surface layer of the workpiece. Calculation errors are mainly caused by difficulties in determining the value of the neck radius.
The paper uses a more accurate method of determining stress and strain states in a tensile test at any time and place of cylindrical samples, using the finite element method (FEM) [62,63]. This allows the development of the required material-hardening curve in the form of dependencies between effective true stress and effective true strain, and to determine the required mechanical parameters of the workpiece.

Calculation the Mechanical Material Parameters Using Semi-Empirical Method
The application Numerical tensile test of cylindrical specimen is developed in APDL in the ANSYS/LS-Dyna system [64]. The application allows for a time analysis of the states of displacements, strains, and stresses in any point of the specimen and at any time t in the duration of the tension test, with the history taken into consideration, for the following data: the material-related factors of the specimen: Young's module, Poisson ratio, initial yield stress, strain hardening modulus, sensitivity to strain rate, strain rate hardening modulus, true failure strain, boundary conditions, -the geometrical factors of the specimen, -dynamic or static load. -the initial state of strains and stresses in the specimen, The essence of this method is to evaluate the required mechanical parameters of the material over a few stages. The first stage is experimental research on the universal testing machine ( Figure 10). A tensile test and measurement of force F D are carried out, as well as elongation ∆L. Then, specimen geometry and displacement (grid) (Figures 12a and 13a) or deformation (specks) of mapped coordination points (Figures 12b and 13b) were measured. Then, the surface area A = π·r 2 of the minimal section of the specimen and the fracture average axial stress σ f at rupture also are calculated as the ratio of the tensile force F D to the surface area F D − σ f = F D /A. N.E. [24], Siebel E. [25], and Szczepiński W. [26]. The effective true stress σ in the speci-  Grid node displacements or speck deformation during the tensile testing are the basis for evaluation of local true strain. The basic mechanical properties of the investigated steel are also determined. Then, in the developed application Numerical tensile test of cylindrical specimen in the ANSYS/LS-Dyna ver. 2019R3 system, numerical analysis of the tensile test for the cylindrical specimen is performed.
The ability to analyze the tensile test of the specimen step by step not only allows for the determination of stress and strain at anytime and anywhere in the specimen, but also allows for the construction of effective true stress-effective true strain curves σ Y − ϕ e . In Figure 14, curves were developed on the ground of the proposed semi-empirical method are presented, compared and calculated based on the existing methods. Grid node displacements or speck deformation during the tensile testing are the basi for evaluation of local true strain. The basic mechanical properties of the investigated stee are also determined. Then, in the developed application Numerical tensile test of cylindrica specimen in the ANSYS/LS-Dyna ver. 2019R3 system, numerical analysis of the tensile tes for the cylindrical specimen is performed.
The ability to analyze the tensile test of the specimen step by step not only allows fo the determination of stress and strain at anytime and anywhere in the specimen, but als allows for the construction of effective true stress-effective true strain curves σ φ . I Figure 14, curves σ φ for C45 steel non-linear hardening ∂σ φ ; φ ∂φ ≅ tgα const. were developed on the ground of the proposed semi-empirica method are presented, compared and calculated based on the existing methods. Additionally, as a result of the conducted research, the values of constants in one parameter functions K φ in the model (9) and K φ in the model (11) were determined φ 0.049, a 0; b 0; c 0; d 1; e 0; f 0; g 0; h 0.02; n 0.39; n 0.6 which finally have the following forms: Additionally, as a result of the conducted research, the values of constants in oneparameter functions K ϕ in the model (9) and K . ϕ in the model (11) were determined: ϕ 0 = 0.049, a = 0; b = 0; c = 0; d = 1; e = 0; f = 0; g = 0; h = 0.02; n 2 = 0.39; n 3 = 0.6, which finally have the following forms: while the graphs of these functions are shown in Figure 15.  (24) and (25), respectively.

Effective True Stress-Effective True Strain Models
In numerical calculations, in order to calculate the true strain and true stress states at each incremental time step, it is necessary to know the yield stress model of the deformed material. Two material models with non-linear hardening were used in the paper: the first, according to Formulas (9)- (11), and the second, the Cowper-Symonds model, which is normally used in the ANYSYS program [64]. This model takes into account mixed isotropic-kinematic or kinematic plastic hardening and the effect of effective plastic true strain and an effective of plastic true strain rate, according to the Power Relation: where -K = 1; K ϕ = R e + β·E p (ϕ ), according to the following formula: where modulus E T (ϕ (VP) e ) dependent on the instantaneous effective true strain: In order to determine the remaining constants, K, n, andn 1 in model (9), it was first linearized by logarithmization, obtaining the following Regression Equation: ln σ Y (ϕ (VP) e , . ϕ (VP) e ) = lnK + n·lnK ϕ + n 1 ·lnK . ϕ (30) where K ϕ and K .
Equation (30), in the system of real variables, can be written in the general form: wherê However, after coding the input variables x 1 and x 2 according to formula [63]: Equation (31) takes the form:Ŷ where k o , k 1 , and k 2 -unknown coefficients of the Regression Equation, x i -input coded variables,x i ∈ [−1; +1]; and i = 1, 2. Then, the experiment was performed according to the two-level plan and with the actual and coded values of the parameters listed in Table 3. Table 3. The two-level experiment plan.

Real Values
Coded Values Five-fold repeatability of the tests was used. This task required the following steps to be conducted [65,66]:

1.
Determination of variability range of the studied parameters.

2.
Choice of the class of the mathematical model.

4.
Gathering the experiment results.
Calculating the inter-row variance and standard deviation. 7.
Checking the homogeneity of variance. 8.
Calculating the coefficients of regression function.

9.
Statistical analysis of the regression function. 10. Examination of the significance level of the correlation coefficient. 11. Checking the adequacy of the mathematical model. 12. Decoding the regression function.
Using matrix calculus, the column vector {k} = [k 0 k 1 k 2 ] T of the unknown coefficients in Equation (32) was calculated from the matrix formula: where: -X -input variable matrix of dimension N×L, for data N = 4 and L = 3, -covariance matrix, which reduces to the form: Y -column vector of the average values of the experimental results.
The boundaries of the confidence region for Regression Function (30) were determined from the following formula: form) and its transposition: 1 4 x T x = 1 4 1 + x 1 2 + x 2 2 ŷ i -average values of model outputs for plan points calculated from Equation (30), y i -average values of experimental results. The test results after statistical processing according to the algorithm presented above were used to develop regression equations in the form (30), and after decoding and delarithmization, the following power form of the function for the yield stress of steel C45 was finally obtained: where K = 980[MPa]; n = 0.2322; n 1 = 0.1222 weredetermined to be the values of the parameters of the Regression Equation.
The model, according to formula (36), does not require the calculation of material hardening modules at each incremental step, as in the case of the Cowper-Symonds model, because the influence of hardening is already included in the values of true strain. The material parameters in the visco-plastic yield stress model (9) were determined on the basis of the results of experimental tests and it was proven that they have a significant impact on the accuracy of the simulation.
The material constant values K and n in the yield stress model (9) were also verified using the Heyer method [67]: In an analogous way, as with model (36), the values of the constants in the Cowper-Symonds model (27) were calculated and the following final form was obtained: After substituting the dependence (40) into Formula (39), the equation that is convenient to use in numerical calculations was obtained: All the material parameters in the visco-plastic yield stress models (9)- (11) and (27), necessary for numerical calculations, are summarized in Table 4.
Due to the material non-linearity of the thread rolling process, in the process of step-614 by-step numerical calculations, the use of the yield stress model (39)  for C45 steel are presented graphically in Figure 16.  (28) and (27), 623 respectively, and from approximation Equations (40) and (41), respectively.  (28) and (27), respectively, and from approximation Equations (40) and (41), respectively. where: -mass density of the object material; E-Young's modulus; ν-Poisson's ratio; K-hardening parameter; Re-initial yield stress; ϕ 0 -initial strains; n, n 2 -hardening coefficients depending on the strain; n 1 , n 3 -hardening coefficients depending on the strain rate, m; P-material constants determining the sensitivity of the material to the rate of plastic deformation; C-parameter dependent on the strain rate; a-d coefficients depending on the strain; e-h-coefficients depending on the strain rate; ϕ f -true failure strain.
Dependence of yielding stressesσ Y ϕ (VP) e , . ϕ (VP) e for C45 steel in the whole range of variability of true effective strain and true effective strain rate in the thread rolling process determined by the semi-empirical method are graphically presented in Figure 17.  Figure 17) was compared with the relationships based on Equations (36) and (42). It was found that functions (36) and (42) describe this relationship very well, and the correlation coefficient is close to unity (R = 0.985). For example, in Figure 18, the graphs obtained on the basis of Equations (36) and (42)

Numerical Evaluation of Mechanical Parameters
The verification of the determined material parameters was carried out by numerical analysis of the tensile strength of a standardized cylindrical sample and a comparison of the sample geometry obtained numerically with those measured experimentally. For this purpose, an original application called "Stretching a cylindrical sample" in the APDL language in the ANSYS/LS-Dyna program was developed. Figure 19 presents a discrete model of cylindrical specimen applied in the numerical simulation procedure. The model was divided into NE = 59, 200 finite elements and NN = 59, 749 nodes with the aid of eight-node elements of the SOLID164 type. As input data for the simulation, material parameters previously determined from the experiment for the power model were introduced. Calculations were performed in 180 steps. In order to control the convergence of the solution, relative values relative to the stress and strain tolerances were established. Figure 20 shows maps of effective stress and effective strain in the specimen at steps 178, 179, and 180. The highest effective stresses occur in the neck before the fracture inside the specimen and are σ e = 755 [MPa] (Figure 20a). However, the plastic strain of the specimen is the largest at the center of the neck and is ϕ      Figure 20 shows selected stages of numerical analysis of the sample stretching process. Figure 20a,b already shows the stage of formation of the so-called "necks". It is the next stage after the stage of elastic deformation. The material flowed and the neck formed. This was followed by plastic hardening of the material. When the narrowing of the so-called "neck" occurs, there is a stress concentration in this area. The strengthening capacity of the material is no longer high enough to overcome the loss in diameter of the material. After exceeding the permissible deformations, the sample breaks in its smallest cross-section (Figure 20c,d). The greatest stresses occur under the surface of the sample at the fracture site. It is in this area that cracking occurs first, and then it moves outward from the sample.
The ability to analyze the tensile test of the specimen step by step not only allows the determination of stress and strain at anytime and anywhere in the specimen, but also for the construction of curves true stress-true strain σ e − ϕ (VP) e . In the Figure 17 were presented curve for C45 steel and non-linear hardening developed on the basis of the proposed method. In the calculation, the geometrical dependences on the neck of the specimen were used for step 178 (crack initiation inside specimen), step 179 (crack growth), and step 180 (final fracture) ( Figure 20). The obtained values of the effective stresses during sample fracture were compared with the values of stresses calculated on the basis of existing methods (Table 5). During the thread rolling process, the rolled shaft copies the shape of the profile of the rolls, so some of it moves into the free space between the flanks of the thread profile on the rolls. The degree of space filling (total or partial) depends on many factors, including, among other things, to the outside diameter of the semi-finished product and the rolling interference, and leads to obtaining a thread with a full or incomplete profile. Each rolling ring of the roll occupies a specific position relative to the theoretical outline of the thread to rolled and the surface of the threaded shaft with an outer diameter d 2 determining the rolling allowance ( Figure 21). During rolling, individual rings shape a thread profile suitable for their penetration into the material. At a given moment, this outline is shaped only on a specific part of the circumference, and the loads affect the metal with a frequency depending on the number of rolls and its rotational speed. When the ring of rollers has moved relative to the workpiece, elastic deformation occurs at the beginning, consisting of temporary changes in the distance between atoms in the volume of the crystal lattice and intra-crystalline and inter-crystalline shifts. As a result of the further impact, irreversible displacements of atomic lattice defects occur, causing plastic deformation with the material squeezed out to the sides. To a large extent, the material moves in the radial direction (into the free spaces between the roll rings), gradually increasing the height of the thread. On the other hand, the remaining part of the material is displaced in the axial direction by the action of the feed and the side surfaces of the rolling rings. Until then, it is assumed that the state of displacement and deformation is plane, i.e., it occurs in planes parallel to the axial plane of the object. When analyzing the shaping of the thread profile during rolling, it can be assumed that the thread profile is symmetrical in various phases of the process. In the initial phase of rolling ( Figure 22), a clear difference is noticeable between the displacement of the extreme volumes compared to the middle part of the outline. This is due to the relatively small penetration of the first rolling rings into the workpiece, resulting in an insignificant impact on the central part of the profile. As the next rings penetrated, the entire thread outline was gradually enlarged. The change in the shape and width of the valley on the top results from the increase in the filling of the free space between the roll rings. At the same time, the influence of the lubricant was noticed, which facilitates the movement of the material along the through surfaces of the rolls, causing a slight increase in the height of the thread profile (Figure 22 top) concerning dry rolling, during which the movement of the material is inhibited (Figure 22 bottom). The overlapping of successive outlines obtained after the passage of the first and subsequent rings allows for a schematic representation of the formation of the thread (Figure 21). The height H i of each of the obtained outlines consists of two quantities: 1. penetration w i of each tip of the roll rolling ring into the shaft, measured from the diameter d 2 of the semi-finished product, 2.
radial outflow h i -displacement of the deformed material in the direction opposite to the axis of the thread.
There is a positive correlation between the ring penetration w i and the height of the flash h i for each group of workpieces with similar physical and mechanical properties [9].
Since in the case of rolled threads, an incomplete profile is most often formed, while a full profile is a borderline case, it can assume that the thread is in a plane state of deformation of the material, which was displaced parallel to the axial plane (Ox 2 x 3 , see Section 3.1.1) of the object. To verify this assumption, experimental investigations were carried out on samples in the form of cylindrical samples made of C45 steel.
In the semi-finished product, as a grinded shaft, longitudinal axial grooves were created, which were filled by flat steel C45 (Figure 23a) and copper M1E (Figure 23b). On the bars, the thread was made by the axial method on the stand with an angular head made by FETTE (Figure 24), which contains three rolls ( Figure 25). View samples carried out in the thread are shown in Figure 26.   (a) (b) Figure 25. View of the processing system of the FETTE rolling mill for thread-rolling using the axial method (a) and view of the rolling head using the axial method (b): 1-roller, 2-hydraulic holder, 3-roller axis, 4-thread diameter adjustment mechanism, 5-workpiece, 6-hydraulic vice, 7three-roller head, 8-head closing and opening system, 9-cooling and lubricating system.   Figure 25. View of the processing system of the FETTE rolling mill for thread-rolling using the axial method (a) and view of the rolling head using the axial method (b): 1-roller, 2-hydraulic holder, 3-roller axis, 4-thread diameter adjustment mechanism, 5-workpiece, 6-hydraulic vice, 7-three-roller head, 8-head closing and opening system, 9-cooling and lubricating system.

Analysis of the Thread Rolling Process in Real Conditions
After the thread rolling process, the cross-sections of transverse to the axis of the shaft were made in three planes passing respectively through the threadless surface (sections I-I), the top (sections II-II), and the bottom of the thread (sections III-III) (Figure 27).
introduced discontinuity of the material, and not from the spatial state of deformation. 781 The proportion of material at upper part of the profile is only 3.5% of the total thread 782 profile. Thus, with sufficient accuracy for engineering practice, it can be assumed that in 783 the rolled thread material there is a spatial state of stress and plane states of displacement 784 and strain in planes parallel to the axial plane (Ox x , see section 3.1.1) of the object. 785 Then, the displacements and deformations of the material during the rolling process 786 can be treated as a series of single identical processes which are characterised by a small 787 contact area (compared to the dimensions of the object and tool) between the active sur-788 face of the roll. 789 Initially, the contact of the tool with the workpiece occurs over small lengths and 790 deformations occur in small volumes. Then, as the process progresses, these zones en-791 large, and in the final phase, they cover the entire outline of the thread. Since the length 792 of contact and the volume of the deformed material for a single tooth are relatively small 793 compared to the dimensions of the object and the cylinder, the profile-shaping process 794 can be considered as the contact of a rigid tool with a deformed semi-infinite body. It has been found that the rectilinear groove planes (1) and (2), which are parallel to the axial plane (Ox 2 x 3 , Section 3.1.1), also remain parallel after thread rolling. It can be assumed that plane states of displacement and strain occur during thread rolling. Additionally, in cross-sections, the side planes of the grooves remain almost linear. The visible slight curvature of the plane in the upper part of the thread profile results rather from the introduced discontinuity of the material, and not from the spatial state of deformation. The proportion of material at upper part of the profile is only 3.5% of the total thread profile. Thus, with sufficient accuracy for engineering practice, it can be assumed that in the rolled thread material there is a spatial state of stress and plane states of displacement and strain in planes parallel to the axial plane (Ox 2 x 3 , Section 3.1.1) of the object.
Then, the displacements and deformations of the material during the rolling process can be treated as a series of single identical processes which are characterised by a small contact area (compared to the dimensions of the object and tool) between the active surface of the roll.
Initially, the contact of the tool with the workpiece occurs over small lengths and deformations occur in small volumes. Then, as the process progresses, these zones enlarge, and in the final phase, they cover the entire outline of the thread. Since the length of contact and the volume of the deformed material for a single tooth are relatively small compared to the dimensions of the object and the cylinder, the profile-shaping process can be considered as the contact of a rigid tool with a deformed semi-infinite body.
Under such assumptions, i.e., for a plane strain state, in the Cartesian coordinate system (0xyz) as in Figure 28, the radial strains ε x and axial strains ε z are significant, while the circumferential (tangential) strain is negligible ε y ≈ 0. This condition, assuming the incompressibility of the object material, leads to the following relationship: ε z = −ε x .
where H is the distance of the axis of the rollers from the axis of the object. 814 The length of the base of the isolated solid is equal to the pitch of the thread P. 815 The planes π and π also limit the part of the active surface of the roller ring with 816 the diameter D in contact with the considered profile ( Figure 29). As in the case of an 817 object, any k-th ring of a cylinder is treated as a set of elementary volumes (solids) with a 818 very small width Δx Δx → 0 . Then, the curvature of the top of the ring (with radius 819 D /2) and its variable penetration into the object over the width can be neglected. 820 The separated part of the cylinder, penetrating into the object to a depth equal to: where w is the depression of the roller in the axial section, causes the displacement of its 822 material parallel to the plane x x , shaping the outline of the thread. Any k-th thread profile is treated as a set of elementary volumes (solids) with a very small width ∆x 1 (∆x 1 → 0). Then, the inclination of the thread profile (resulting from the thread pitch), the outline of its top (with radius d/2), and the variable depression in the workpiece of the cylinder ring at the width ∆x 1 can be neglected. Each such volume is separated by two planes, π 1 and π 2 , parallel to the axial plane x 2 x 3 and spaced from it by the values of x 1 and x 1 + ∆x 1 , respectively ( Figure 28). The height h 3 of the thread profile in this section is: where H is the distance of the axis of the rollers from the axis of the object. The length of the base of the isolated solid is equal to the pitch of the thread P. The planes π 1 and π 2 also limit the part of the active surface of the roller ring with the diameter D w in contact with the considered profile ( Figure 29). As in the case of an object, any k-th ring of a cylinder is treated as a set of elementary volumes (solids) with a very small width ∆x 1 (∆x 1 → 0). Then, the curvature of the top of the ring (with radius D w /2) and its variable penetration into the object over the width ∆x 1 can be neglected.
The separated part of the cylinder, penetrating into the object to a depth w x equal to: where w is the depression of the roller in the axial section, causes the displacement of its material parallel to the plane x 2 x 3 , shaping the outline of the thread.  Figure 29. Scheme of thread rolling process using four rolls by axial method (real object) (a); schem of the division of the thread profile and roller rings into solids with a width of dx1 (b); cross-secti of the k-th thread profile and the k-th roller ring with the plane π (c); numerical simulation mod implementation (d).

Model Tests of the Thread Forming Process
Assuming a plane deformation state and a spatial state of stress, in order to bett understand the physical phenomena occurring in the thread rolling process and to dete mine the mechanism of displacement and strain of the material, model tests were carri out, i.e., on an appropriately enlarged geometric model and with the use of model mat rial.
On the basis made in Section 3.1.1 of the analysis of thread rolling in real condition in which there is a spatial state of stress and a flat state of deformation, a stand for mod testing of the object (Figure 29d) was developed in which an identical state of stresses an strains of the model material is produced. The purpose of the model tests carried out, i. on a properly enlarged geometric model and with the use of model material, was to bett understand the physical phenomena occurring in the process of thread rolling and to d termine the mechanism of displacement and deformation of the material.
In order to ensure the condition of rheological similarity, yield stresses were dete mined for the model material-plasticine. This is due to its availability, low price, ea machinability, plasticity (yielding stresses are 100-1000 times lower than the correspon ing metal stresses), the possibility of its modification and regeneration, and thus, its r peated use. The basic component of plasticine is calcium carbonate CaCO3 and bindin agents such as water, mineral, vegetable and animal fats, and other softening and har ening agents and dyes. By selecting the appropriate plasticine composition, it is possib to create various rheological models of bodies (i.e., perfectly plastic, rigid-plast strengthening, etc.) as well as multiple versions of the same rheological model wi Figure 29. Scheme of thread rolling process using four rolls by axial method (real object) (a); scheme of the division of the thread profile and roller rings into solids with a width of dx 1 (b); cross-section of the k-th thread profile and the k-th roller ring with the plane π 1 (c); numerical simulation model implementation (d).

Model Tests of the Thread Forming Process
Assuming a plane deformation state and a spatial state of stress, in order to better understand the physical phenomena occurring in the thread rolling process and to determine the mechanism of displacement and strain of the material, model tests were carried out, i.e., on an appropriately enlarged geometric model and with the use of model material.
On the basis made in Section 3.1.1 of the analysis of thread rolling in real conditions, in which there is a spatial state of stress and a flat state of deformation, a stand for model testing of the object (Figure 29d) was developed in which an identical state of stresses and strains of the model material is produced. The purpose of the model tests carried out, i.e., on a properly enlarged geometric model and with the use of model material, was to better understand the physical phenomena occurring in the process of thread rolling and to determine the mechanism of displacement and deformation of the material.
In order to ensure the condition of rheological similarity, yield stresses were determined for the model material-plasticine. This is due to its availability, low price, easy machinability, plasticity (yielding stresses are 100-1000 times lower than the corresponding metal stresses), the possibility of its modification and regeneration, and thus, its repeated use. The basic component of plasticine is calcium carbonate CaCO 3 and binding agents such as water, mineral, vegetable and animal fats, and other softening and hardening agents and dyes. By selecting the appropriate plasticine composition, it is possible to create various rheological models of bodies (i.e., perfectly plastic, rigid-plastic, strengthening, etc.) as well as multiple versions of the same rheological model with different plastic resistances [27]. The material model was developed based on the static compression test.
In order to compare the obtained material characteristics of the model material with the characteristics of C45 steel, they were presented in the diagram (Figure 30). Identical inclination angles of the lines in the diagram in the logarithmic coordinate system indicate that the obtained material characteristics for the model material are similar to the material characteristics of C45 steel. Based on the values of the initial yield point determined in strength tests for the model material and steel, which are R * e = 45. 63 [kPa] (plasticine), R e = 440 [MPa] (C45 steel), the value of the rheological scale factor is k r = 440, 000/45.63 = 9643. different plastic resistances [27]. The material model was developed based on the sta compression test.
In order to compare the obtained material characteristics of the model material wi the characteristics of C45 steel, they were presented in the diagram (Figure 30). Identic inclination angles of the lines in the diagram in the logarithmic coordinate system indica that the obtained material characteristics for the model material are similar to the materi characteristics of C45 steel. Based on the values of the initial yield point determined strength tests for the model material and steel, which are R * 45.63 kPa (plasticin R 440 MPa (C45 steel), the value of the rheological scale factor is k 440,00 45.63 9643. Flat states of displacement and strain as well as spatial states of stress during mod tests of the thread rolling process were obtained by appropriately shaping the tool, punc and compressed sample (Figure 31a). The two-part sample, shown in Figure 31b, w made of plasticine. On the dividing surface of each half, a grid with shapes, mesh size and location of line families was applied, as in Figure 31b. Intersections of grid lines crea material points, from whose displacements the kinematics of the process are determine Grid lines were applied with a special device.   Flat states of displacement and strain as well as spatial states of stress during model tests of the thread rolling process were obtained by appropriately shaping the tool, punch, and compressed sample (Figure 31a). The two-part sample, shown in Figure 31b, was made of plasticine. On the dividing surface of each half, a grid with shapes, mesh sizes, and location of line families was applied, as in Figure 31b. Intersections of grid lines create material points, from whose displacements the kinematics of the process are determined. Grid lines were applied with a special device.  [27]. The material model was developed based on the sta compression test. In order to compare the obtained material characteristics of the model material w the characteristics of C45 steel, they were presented in the diagram (Figure 30). Identi inclination angles of the lines in the diagram in the logarithmic coordinate system indic that the obtained material characteristics for the model material are similar to the mater characteristics of C45 steel. Based on the values of the initial yield point determined strength tests for the model material and steel, which are R * 45.63 kPa (plasticin R 440 MPa (C45 steel), the value of the rheological scale factor is k 440,00 45.63 9643. Flat states of displacement and strain as well as spatial states of stress during mod tests of the thread rolling process were obtained by appropriately shaping the tool, pun and compressed sample (Figure 31a). The two-part sample, shown in Figure 31b, w made of plasticine. On the dividing surface of each half, a grid with shapes, mesh siz and location of line families was applied, as in Figure 31b. Intersections of grid lines cre material points, from whose displacements the kinematics of the process are determine Grid lines were applied with a special device.   The coordinates of the nodes were determined in the Cartesian reference system associated with the sample. In the subsequent stages of the process ("step by step") ( Figure 32), i.e., for different values of the punch depth, deformations of the applied mesh were observed and recorded (photographed), and then measured using a WERTH optical microscope with a digital readout with accuracy to 0.001 mm.
The coordinates of the nodes were determined in the Cartesian reference system as-874 sociated with the sample. In the subsequent stages of the process ("step by step") ( Figure   875 32), i.e., for different values of the punch depth, deformations of the applied mesh were 876 observed and recorded (photographed), and then measured using a WERTH optical mi-877 croscope with a digital readout with accuracy to 0.001 mm. In model investigations, by changing the surface roughness of the punch, it is also possible to produce different friction conditions, e.g., different values of the friction coefficient (µ ≈ 0 ÷ 0.39) in the contact areas of the punch with the model material. In Figure 33, the results of model investigations of the process of shaping the trapezoidal thread profile for various friction coefficients are shown. The coordinates of the nodes were determined in the Cartesian reference system as-874 sociated with the sample. In the subsequent stages of the process ("step by step") (  Figure 33 shows a significant influence of the friction coefficient on the state of displacements of nodal points, and thus on the state of material deformation. In the absence of friction in the contact area of the tool with the workpiece (Figure 33a), when forming the thread profile, the material is not inhibited by the punch and slides freely along the contact plane. The lack of curvature of the vertical lines of the finite element mesh is visible. An increase in the coefficient of friction results in an increase in the inhibition of material movement. For high values of the friction coefficient (Figure 33b,c), the material is strongly braked in the contact area. There are also areas of material adhesion. There is then greater movement of material in areas further away from the contact zone. This is shown by the curvature of the mesh vertical lines towards the bottom of the thread, mainly in areas close to the contact zone.
A numerical application in the ANSYS/LS-Dyna system was also developed, which enables a comprehensive time analysis of deformation states (displacements and strains) occurring in a model object consisting of a rectangular sample and a punch, for spatial states of stress and plane states of displacement and strain. The application allows for the simulation of the process using two methodologies. The first method consists of introducing given boundary conditions in the contact area for displacements determined in model visioplasticity tests (Figure 33). In the second methodology (unknown boundary conditions in the contact area), the incremental theory was used, performing step-bystep calculations. In this case, appropriate modeling of the contact area was performed. The second methodology additionally allows for the determination of the distribution of pressures in the contact area, which was not possible with the first methodology. In this case, the contact area was modeled using CONTA172 and TARGE169 finite elements. The CONTA172 element defines an element that is in contact and slip with the deformed surface defined by the TARGE169 element. The results of the computer simulation were compared with the results of the model investigation, obtaining a very good agreement ( Figure 34). strongly braked in the contact area. There are also areas of material adhesion. There is 897 then greater movement of material in areas further away from the contact zone. This is 898 shown by the curvature of the mesh vertical lines towards the bottom of the thread, 899 mainly in areas close to the contact zone.

Sensitivity Analysis
System sensitivity analysis is defined as a measure of the change in the system response caused by a change in a selected parameter called the decision variable. In the case of FEM calculations of the TR process, an important step is to determine the sensitivity of the effective maximum strains and stresses at discrete points of the workpiece to: the change in the dimensions of the finite element-shape factor (SF), the number of finite elements (NE), and the shape element function (SEF).

The Influence of the Shape Factor of Finite Elements (SF)
The shape factor of the finite element SF was defined as the ratio of the height B to the width A of the element: SF = B/A. It is desirable that the aspect ratio be close to unity. Due to the profile of the thread (strong geometric non-linearity, especially at the bottom and top of the thread), discrete models are used, built of finite elements with a large aspect ratio. In order to reduce the number of degrees of freedom of the model, the number of finite elements is minimized by increasing the height of the elements. In order to determine the rational shape of the FE, computer simulations were carried out to determine the impact of SF on the distribution of effective stress and strain and the accuracy of the thread profile mapping. Due to the symmetry of the model, half of the outline was considered. The analysis was carried out using a numerical application developed in the APDL language in the ANSYS/LS-Dyna system. An eight-node element with a non-linear shape function of the PLANE183 type was used to discretize the object. In order to determine the maximum equivalent stresses and strains according to the non-linear Huber-Mises-Hencky (HMH) hypothesis in a trapezoidal thread as a function of the SF, model variants with different division grids were prepared. The following shape factor values were adopted: SF = 0.11; 0.25; 0.43; 0.66; 1; 1.5; 2.33; 4; and 9. Example results of the simulation of the impact of SF on the accuracy of the tool mapping in the workpiece for the friction coefficient µ ≈ 0 on the contact surface, for the case of a trapezoidal thread, are shown in Figure 35. It was found that the most favorable results of the calculated effective stresses and strains were obtained for SF = 1, and therefore, such a value of SF was adopted for further analyses.
impact of SF on the distribution of effective stress and strain and the accuracy of the thread 936 profile mapping. Due to the symmetry of the model, half of the outline was considered. 937 The analysis was carried out using a numerical application developed in the APDL lan-938 guage in the ANSYS/LS-Dyna system. An eight-node element with a non-linear shape 939 function of the PLANE183 type was used to discretize the object. In order to determine 940 the maximum equivalent stresses and strains according to the non-linear Huber-Mises-941 Hencky (HMH) hypothesis in a trapezoidal thread as a function of the SF, model variants 942 with different division grids were prepared. The following shape factor values were 943 adopted: SF 0.11; 0.25; 0.43; 0.66; 1; 1.5; 2.33; 4; and 9. Example results of the simula-944 tion of the impact of SF on the accuracy of the tool mapping in the workpiece for the fric-945 tion coefficient μ 0 on the contact surface, for the case of a trapezoidal thread, are 946 shown in Figure 35. It was found that the most favorable results of the calculated effective 947 stresses and strains were obtained for SF 1, and therefore, such a value of SF was 948 adopted for further analyses.

The Influence of the Number of Finite Elements (NE)
The next stage of the sensitivity analysis was to determine the impact of the finite element mesh density on the calculation results. In order to determine the mapping of the tool in the workpiece and to determine the maximum effective stresses and strains according to the non-linear HMH hypothesis in the trapezoidal thread as a function of the mesh density, variants of models with different numbers of finite elements were prepared. Models containing finite elements were adopted for the numerical calculations: NE = 220; 920; 3680; 8280; 22, 280; 32, 880; and 45, 600, for SF = 1 and µ ≈ 0.
The impact of the number of finite elements NE on the accuracy of the tool outline mapping is shown in Figure 36. It can be concluded that the use of a minimum NE = 22, 280 of elements already ensures the correct mapping of the active surface of the tool on the shaft surface. In addition, the increase in the number of finite elements results in a more accurate mapping of the areas of contact between the object material and the tool.
impact of SF on the distribution of effective stress and strain and the accuracy of the thread 936 profile mapping. Due to the symmetry of the model, half of the outline was considered. 937 The analysis was carried out using a numerical application developed in the APDL lan-938 guage in the ANSYS/LS-Dyna system. An eight-node element with a non-linear shape 939 function of the PLANE183 type was used to discretize the object. In order to determine 940 the maximum equivalent stresses and strains according to the non-linear Huber-Mises-941 Hencky (HMH) hypothesis in a trapezoidal thread as a function of the SF, model variants 942 with different division grids were prepared. The following shape factor values were 943 adopted: SF 0.11; 0.25; 0.43; 0.66; 1; 1.5; 2.33; 4; and 9. Example results of the simula-944 tion of the impact of SF on the accuracy of the tool mapping in the workpiece for the fric-945 tion coefficient μ 0 on the contact surface, for the case of a trapezoidal thread, are 946 shown in Figure 35. It was found that the most favorable results of the calculated effective 947 stresses and strains were obtained for SF 1, and therefore, such a value of SF was 948 adopted for further analyses.

The Influence of the Shape Element Function (SEF)
In practice, increasing the accuracy of calculations while reducing the number of finite elements used for discretization can be obtained by increasing the degree of the polynomial, the so-called shape element function (SEF). In order to achieve such an effect, higher-order elements were used in the analyses (elements with eight nodes with a square shape element function-PLANE 183 type) and for comparison, lower-order elements were used (elements with four nodes with a linear shape element function-PLANE42). Example results of calculations of equivalent stresses and strains in a trapezoidal thread are shown in Figure 37.
nite elements used for discretization can be obtained by increasing the degree of the pol-970 ynomial, the so-called shape element function (SEF). In order to achieve such an effect, 971 higher-order elements were used in the analyses (elements with eight nodes with a square 972 shape element function-PLANE 183 type) and for comparison, lower-order elements 973 were used (elements with four nodes with a linear shape element function-PLANE42). 974 Example results of calculations of equivalent stresses and strains in a trapezoidal thread 975 are shown in Figure 37.  In the developed application in the ANSYS system, the object was modeled as a plane 990 (2D) state of deformation and a spatial state of stress [10,11,16,18]. Due to the symmetry, 991 half of the profile was considered. The tool was treated as a perfectly rigid body (E→∞) 992 ( Figure 38), while the shaft was treated as a non-linear elastic/plastic body with non-linear 993 hardening, which was approximated by the Cowper-Simonds power model (37)

Effective Discrete Model for Trapezoidal Threads
On the basis of the sensitivity analysis carried out and additionally taking into account the friction conditions in the contact areas of the tool with the workpiece (µ = 0), the following effective model was adopted for further numerical calculations: shape function SF = 1; number of finite elements NE = 32, 280; number of nodes NN = 99, 433; number of degrees of freedom DOF = 198, 866 and non-linear of the shape element function-PLANE183 type ( Table 6). In the developed application in the ANSYS system, the object was modeled as a plane (2D) state of deformation and a spatial state of stress [10,11,16,18]. Due to the symmetry, half of the profile was considered. The tool was treated as a perfectly rigid body (E→∞) (Figure 38), while the shaft was treated as a non-linear elastic/plastic body with nonlinear hardening, which was approximated by the Cowper-Simonds power model (37) ( Table 4). The object was discretized with eight-node PLANE183 finite elements with a non-linear shape function. In the contact area, the boundary conditions for displacements are unknown. The application allows for the simulation of the process using the two methodologies presented in Section 3.1.2. In the first method, boundary conditions for the tool and the shaft were selected in accordance with the methodology presented in Section 3.1.1 shown in Figure 29, where the thread rolling depth was selected on the basis of equation (39). The displacement in the y direction was loaded on the tool to a depth equal to half the height for the Tr12×3 thread with trapezoidal outline w y = u y = 0.75 mm ( Figure 38, Table 7). In the second method, the contact of bodies was modelled with the finite elements of the TARGE169 and CONTA172 types. On the bottom of the model, the translational and rotational degrees of freedom were taken off. The calculation was carried out for the material properties from Table 4. Three values of the friction coefficient on the contact surfaces µ 1 = 0, µ 2 = 0.2, and µ 3 = 0.39 were assumed and the shape coefficient SF = 1 ( Figure 38, Table 7). The finite elements grid contains NE = 32,880 of finite elements.
shown in Figure 29, where the thread rolling depth was selected on the basis of equation 1000 (39). The displacement in the y direction was loaded on the tool to a depth equal to half 1001 the height for the Tr12×3 thread with trapezoidal outline w u 0.75 mm (Figure 38, 1002 Table 7). In the second method, the contact of bodies was modelled with the finite ele-1003 ments of the TARGE169 and CONTA172 types. On the bottom of the model, the transla-1004 tional and rotational degrees of freedom were taken off. The calculation was carried out 1005 for the material properties from Table 4. Three values of the friction coefficient on the 1006 contact surfaces μ 0, μ 0.2, and μ 0.39 were assumed and the shape coefficient 1007 SF 1 (Figure 38, Table 7  The coefficient of friction has a significant impact on the value and distribution of 1017 equivalent stresses (Table 8)

Simulation Parameters Value
Process temperature 20 • C Type of thread Tr12×3 Thread-rolling depth w y = u y = 0.75 mm Shaft Material C45 Roller (tool material) a perfectly rigid body, E→∞ Type of Finite Elements for shaft and roller PLANE 183 Type of Finite Elements for contact TARGE169 and CONTA172 Variety of the friction coefficient µ 1 = 0, µ 2 = 0.2, and µ 3 = 0.39 Shape coefficient SF = 1

Results of Numerical Analysis
Exemplary numerical simulation results for trapezoidal threads are shown in Figures 39 and 40. Analysing the distribution of strain intensity, stress, and deformation of the finite element mesh, a brief influence of lubrication conditions was observed. The coefficient of friction also has a significant impact on the value and distribution 1030 of deformations (Table 8) The coefficient of friction has a significant impact on the value and distribution of equivalent stresses (Table 8). For µ = 0, the highest effective stresses are σ e = 2390 [MPa] and occur at the bottom of the thread on the fillet radius, and based on the thread profile, in the zone adjacent to the contact surface (MX, Figure 39a). With the increase of the friction coefficient, it moves upwards along the contour, taking values from σ e = 2410 [MPa] for µ = 0.2 (Figure 39b) to σ e = 3610 [MPa] for µ = 0.39 (Figure 39c). For µ > 0, a local area of minimum stresses (MN) appears at the bottom of the thread on the axis of symmetry, which assumes a characteristic wedge shape. This minimum, along with the increase of the friction coefficient, moves deeper into the material, covering a larger and larger area. The coefficient of friction also has a significant impact on the value and distribution of deformations (Table 8). For µ = 0, the largest effective deformations are ε e = 2.058 and occur (as in the case of stresses) on the fillet radius and the basis of the thread profile (MX, Figure 40a). With the increase of the coefficient of friction, this area moves upwards along the contour, taking values from ε e = 2.803 for µ = 0.2 (Figure 40b) up to ε e = 3.514 for µ = 0.39 (Figure 40c). For µ > 0, the area of the local minimum (MN) appears at the bottom of the thread on the axis of symmetry, which takes a characteristic wedge shape and is the result of material adhesion. The area of this minimum moves further and further away from the contact surface, going deeper into the material, and covering a larger area as the value of the friction coefficient increases.

Conclusions
External cold thread rolling with a trapezoidal profile is a complex process in technological terms. This process was analysed as a geometrically and physically non-linear initial and boundary problem with variable stiffness of the system during the rolling process. Boundary conditions are variable in time and space, while in areas of contact between the tool and the workpiece, they are unknown. The measurement of process parameters, such as a displacement zone, temperature, stress, structural change, etc., during the thread rolling process using known techniques of measurement is impossible. Conclusions can only be reached about the properties of the product after the thread rolling.
The currently used methodologies for calculating stress states during the stretching of cylindrical samples result in an overestimation of bursting stresses by about 5%, while the true strains are underestimated by about three times. The applied new semi-empirical methodology allows for accurate determination of the relationship between true stress and true strain and strain rate. This allowed for an increase in the accuracy of the calculations of strain states, strain rates, and stress in the workpiece for the entire range of strain variability and strain rate during the thread rolling process.
The numerical application of the semi-empirical method proposed allows for more accurate descriptions of the material mechanical properties. It also allows us to determine the moments of cracking (invisible) inside the neck of the sample (Figure 20a,b) and the states of true strain and true stress in the neck.
New constitutive equations were introduced, taking into account the effect of strain and strain rate on material hardening and the non-linear variation of the strain-dependent modulus E T (ϕ e ) and the second strain-rate-dependent modulus . E T ( . ϕ e ) ( Figure 17). The proposed method also makes it possible to accurately determine the failure strain and the value of the spring back value after breaking the sample.
Of the three presented kinematic variants of the rolling process, only the new variant proposed by the Author with a rotating and driven basket ensures the generation of flat states of displacement and strain as well as spatial states of stress in the deformed material.
Model and experimental studies have shown the validity of assuming the presence of a spatial state of stress and plane states of displacement and strain in the rolled thread shaft, occurring in planes parallel to the axial plane. Those made it possible to significantly simplify the numerical model of the process and reduce the calculation time multiple times over.
An application of modern numerical methods and computing systems allows for an analysis of complex physical phenomena occurring in the process under investigation. The application developed in the ANSYS system in the APDL language enables a time analysis of the rolling process with the consideration of the changeability of the roller ring's profile.
Based on the course of physical phenomena in the working zone, the technological quality of the thread can be predicted.
Exemplary numerical analyses show that the friction coefficient significantly affects the states of displacements, deformations, and stresses in the surface layer of the screw and is one of the factors determining its technological and functional quality. The best operational quality of the screw was obtained when rolling with abundant lubrication (µ ≈ 0).
The outline of the rollers and the pattern of their penetration into the material are also important. These factors are decisive for the flow mechanics of the object material on which the thread outline and displacement states depend, as well as for the deformations. The results of the numerical simulation can be the basis for the selection of the outline of the working surface of the rolls and the kinematics of the process.
The results of the numerical calculations confirm the possibility of a correct analysis of the material deformation process during thread rolling, both according to method I (set boundary conditions) and method II, without knowing the boundary conditions in the contact area. For the first method, it was necessary to perform model tests using plasticine. However, the limitations of defining the boundary conditions imposed in the computational program are laborious. The disadvantage of this method is also the inability to determine the pressures in the contact area. These inconveniences are eliminated by the second calculation methodology consisting of the introduction of contact finite elements (CONTA and TARGE). An additional benefit of using the second methodology is the development of the pressure distribution in the contact area, which is not possible experimentally.
The results obtained from the numerical simulation can be used to optimise the thread rolling process, reduce design cycles and costs at the beginning of production, reduce overall production costs, and increase product quality.
The large number of factors affecting the quality of the thread and non-linearity of physical phenomena occurring in the rolling process mean that in order to better understand these factors and determine the mechanism of displacement and deformation of the material, it is advisable to conduct model tests, i.e., on a properly enlarged geometric model and with the use of model material. The developed methodology of model research allows for a better understanding of the physical phenomena occurring in the process of thread rolling. This knowledge is essential for the proper design and control of this complex plastic-forming process.