Simulation of surface hardening in the deep rolling process by means of an axial symmetric nodal averaged finite element

Surface hardening by deep rolling can be considered as the axial symmetric problem in some special events (namely, when large R and small r radii of the deforming roller meet the requirement R>> r). An axisymmetric nodal averaged stabilized finite element is formulated. The formulation is based on a variational principle with a penalty (stabilizing) item in order to involve large elastic-plastic strain and near to incompressible materials. The deep rolling process for a steel rod is analyzed. Axial residual stress, yield stress, and Odkvist’s parameter are calculated. The residual stress is compared with the data obtained by other authors using a three-dimensional statement of the problem. The results obtained demonstrate essential advantages of the newly developed finite element.


Introducton
Plastic deformation is used not only in material forming processes. This physical property of metals and alloys is also applied for improving the surface behaviour under mechanical and chemical influences. As is known, surface affects such characteristics of metallic components as resistance to plastic deformation (surface hardness), resistance to corrosion, fatigue fracture, fretting fatigue, and wear, because most cases of material failure originate from the exterior layers of the workpiece [1]. When surface experiences mechanical contact loads, the role of those characteristics grows dramatically. Mechanical working is one of the ways to improve the surface behavior. Such working has been employed for a long time and the three most popular techniques are: shot peening, burnishing, and deep rolling. Each of them has its own advantages and disadvantages but the deep rolling stands out against a background of other (both mechanical and many non-mechanical) techniques owing to large values of hardness, residual stress, and hardening depth reached. By combining deep rolling and burnishing processes one can elevate the surface hardness by several times for some materials [2,3]. This combined process is shown in Figure 1 and Figure 2 where the burnishing roller follows the deforming one and all the process can be realized by means of a turning lathe.
The present work is dedicated to the construction of a simple mathematical model such that the model can be studied by means of a personal computer. Mechanical surface treatment has attracted considerable scientific interest in the last few decades due to both spread of modern experimental equipment for microstructure investigations and accelerated possibilities of computers. Note that the process of multiple local deforming is the most complicated one in solid mechanics because large gradients of tensor fields superimpose on material non-uniformities induced by plastic flowing. Classical finite elements fail in this case because they are not consistent with the incompressibility requirement. Many improved finite element schemes were developed but in spite of essential progress most of them suffered from volume locking, bending locking, and spurious low-energy modes. Efforts to eliminate these drawbacks have resulted in several new approaches and the simplest of them is the stabilized nodal averaging method [4][5][6]. In the present work, simplicity is the main criterion of choice of an appropriate modern approach. Due to many reasons, finite elements of low order prove to be preferable for non-linear problems and large plastic deformation. Besides, mesh generators produce triangles and tetrahedrons as a rule; more complex figures especially for irregular meshes cause serious problems. The procedure of nodal averaging was first formulated by Dohrmann et al. [4] and further analyzed by Bonet et al. [5]. As it turns out, nodal averaging provides weakened constraints such that coarse mesh accuracy can often be achieved for near incompressible materials. It was also pointed out that the formulation was prone to spurious low-energy modes. In order to suppress them, Pusso and Solberg [6] proposed to enter stabilizing terms into a discrete form of a variational equation. A number of different formulations have been developed to rectify poor performance of low-order elements under incompressible and near incompressible conditions but the approach by Dohrmann et al. [4] alleviates both the volumetric locking and locking due to bending. However, some modern methods give the same result as in Pusso and Solberg [6] -see for example Puso et al. [7]. But, stabilized nodal averaging is the simplest among them owing to a trivial theoretical basis and a short program code.
One more problem is that periodic local deformation requires three-dimensional modelling despite the fact that a workpiece looks like an axisymmetric body. However, 3D simulation is too timeconsuming therefore the task has to be reduced to a simpler one (see Section 3). In contrast to this approach, the Authors propose to consider the deep rolling problem as the axial symmetric one in some special events (namely, when large R and small r radii of the deforming roller meet the requirement R>> r). Thus, if an effective axial symmetric finite element will be developed then modelling of the process can be simplified considerably. In this work, an appropriate nodal (i.e. attached to a mesh node) stabilized finite element is formulated. An approximation scheme based on the simplest triangle element and the averaging procedure over triangles of a nodal complex is proposed and the methodology [6] is extended on axial symmetric problems.

Materials and methods
Let an arbitrary triangulation (i.e. N nodes  Figure 3 illustrates these denotations.  The idea of stabilized nodal averaging can be demonstrated by the simplest example -plane strain problem with small elastic strain. Let ) ( be a FEM-approximation of the displacement field. By n ε denote the average strain at node procedure follows the meshing one and new finite elements appear while each of them is a one-third of a nodal complex -see Figure 2 where one-thirds of triangles (each of triangles is dissected into three parts of equal area) form a new (nodal averaged) finite element. Note that all the nodes except for the central one are located outside the new finite element. As for the axisymmetric event, a volume formula for a ring of triangular profile is not so trivial as for triangular or tetrahedron. This results in too cumbersome expressions for both nodal force and stiffness matrix. In the present work, an approximate formula of volume for such a ring is exploited in order to simplify a program code and to accelerate calculations. Let us proceed to the axisymmetric formulation itself. In cylindrical coordinates = r y 1 , = h y 2 ,  = y 3 , the Cauchy stress tensor and the small strain one take the form It is necessary to rewrite the variational principle (1) taking into account large irreversible (plastic) strain (the Euler approach is adopted throughout the work i.e. all tensor fields refer to the current configuration rather the initial one). Contracting by a pair of indexes, we obtain where:  Equation (2) has to be supplemented by a constitutive equation. As far as metals and alloys refer to near rigid-plastic materials, it is convenient to write that equation by means of any corotational rate "r" (for example, Yaumann's rate "j"): where d is the deformation rate tensor. Therefore, two stress tensors take part in step-by-step numerical analysis: real stress σ and stabilizing stress σ  . Introduce a FEM-approximation for the displacement field u, its gradient Due to the approximate formula for volume Substituting the displacement gradient approximation in variational principle (2), find the nodal force at node p A : where ij  is Kronecker's delta-symbol. Substituting (5) and (6) into (7) and taking into account   . (8) Note that this expression is simpler than (7) because all volumes (both nodal and of rings of triangular profile) disappeared due to the successful choice of the approximate formula for volume m v (see the elucidations to equation (2)). Finally, the stiffness matrix is obtained by means of the nodal force variation. In order to deduce a proper expression, find the stress tensor variation from the constitutive equation with Yaumann's rate: is the vortex tensor. Then the nodal stress variation at node n A is s are also dependent variables and subjected to variation. Therefore, the total expression for the nodal force variation is too cumbersome and consists of two items: symmetrical and non-symmetrical. In practice of metal forming analysis, a non-symmetrical item can be rejected without the loss of convergence. The symmetrical positive-definite item of the stiffness matrix of the nodal finite element is of the form (p and q are node numbers of complex n C ) then some special items should be added to the stiffness matrix respectively: then one more item should be added: When assembling the global stiffness matrix, the stabilizing matrices are also taking into account (they are calculated according to classical expressions [8]). This completes the nodal averaged axisymmetric finite element formulation.

Results and discussion
The problem of simulation for multiple local deforming (including the deep rolling treatment) is very difficult from the computational point of view due to the fact that this process may be accompanied not only by high gradient values of the tensor fields but also by essential transformations of the material microstructure leading to high material non-homogeneity. So, the superposition of the material non-homogeneity with high gradients of tensor's components leads to an extremely complicated isoline pattern for stress, strain and other variables. As a rule, such problems are considered in a three-dimensional statement (though a workpiece looks like an axisymmetric body) when a small "representative" part is cut out the entire workpiece and some "reasonable" boundary conditions are imposed at the cut surface. The developed axial symmetric finite element reduces the total computational time radically and allows considering the whole workpiece rather than its small "representative" part (this may be necessary in some cases, for example in the event of deep penetration of the residual stress). Further, consider a numerical example. A steel rod of radius 20   mm and of length 100  l mm is rolled by two rollers of transverse radius 908 . 0  r mm and of longitudinal one r R  while the rollers move in opposite directions ( Figure 5) . The process shown in Figure 5 is not realistic but this process possesses the bilateral symmetry therefore the total number of nodes may reduced twice. In practice, the working is carried out by a single roller or by two-three rollers rotating in a common plane. The finite element mesh is presented in Figure 6 where the finest discretisation is marked out by black; the size of this region is 5×1mm and the splitting numbers here are 50×10. A contact problem statement conforms to [9,10]. A similar problem in 3D statement was investigated in [11] with the only difference: longitudinal roller's radius was equal to mm. The results of both works (the present and [11]) are compared in Figure 7 (the finite elements meshes correspond to each other), where the axial component of residual stresses is presented (two other components are small).
The negative sign demonstrates the compressive character of this stress. The peak values are almost equal. The difference between the two plots is explained by different values of radius R. It is obvious that if r R  then the contact zone is of the form of a very oblong oval and is much larger than in case 5 . 2  R mm. Correspondingly, the hardening depth is larger too when r R  . It should be pointed out that the three-dimensional model [11] was calculated up to depth 3 mm only and in the present work an extrapolation of those results up to depth 4 mm is given. The developed finite element permits to calculate the model without any restriction in depth and the proper result is presented in Figure 8. It is seen that at depth 3.5 mm the axial component changes its sign and reaches a maximum at depth 6mm. The fractional increase of yield stress (here, it is a measure of hardening) in different cross   sections of the rod is shown in Figure 9, where L is the distance from the rod middle. It is seen that the hardening profiles are not identical hence the deep rolling process is not stabilized. The reason is surface buckling in front of the roller, and this local distortion grows during the process of working. Finally, the isolines of Odkvist's parameter are presented in Figure 10 in the range 0.078-0.78. It is seen, they are not parallel to the surface. This confirms that the process is not stabilized. Thus, it is necessary to model a whole workpiece in order to obtain an adequate result but 3D modelling of periodic local deformation is too time-consuming for modern computers. In contrast to 3D approach, it is proposed to consider the deep rolling problem as the axisymmetric one in some special events (the restriction in radii r R  is required) and to simulate the behaviour of the whole workpiece. In this work, only the axial component of the residual stress tensor is analyzed. As for the other components, one can turn to the results [11,12] where all the components were calculated. Those results indicate that the stress component in the radial direction is not as significant and may be ignored altogether. However, the circular stress component in the direction which coincides with rolling direction is considerable but still is significantly lower than the stress component in the axial direction which is perpendicular to the direction of rolling. Therefore, only the axial stress component is considered for evaluation of the effects of the rolling parameters on the residual stress distribution induced by deep rolling.
Discuss some practical aspects of the work. To fulfil the deep rolling treatment, it is necessary to design an appropriate technological process (numerical modelling is an essential stage of this design) and to manufacture a proper rig (see for example Figure 1 and Figure 2). Therefore, the total cost of a produced detail will depend on the sum of the design cost and the rig manufacturing one. At first sight, a roller of small radius is preferable because its cost is lower in comparison to a roller of large radius. On the other hand, the design will require a three dimensional numerical model in such case because a two dimensional model leads to great errors. The results on ball burnishing modelling [13] indicate that the difference between calculated and measured residual stress values vary by 40% for 2D simulation and only by 7% for 3D one. According to that work, the plane strain assumption for 2D model prohibits the material from flowing along the out-of-plane tangential direction, and thus, compared to 3D simulations, the plastic strain decreases clearly. Consequently, the maximum residual stress introduced is affected. As for deep rolling simulation, the results [13] should be corrected. Indeed, the contact zone for the deforming roller in the deep rolling process is of the oval form and plastic flowing in tangential direction is far less than in axial one if this oval is strongly oblong. Hence, the error of axial stress 2D calculations may be not as large as 40% for a sufficiently large roller. One can conclude from this resume: 3D model is necessary if a small roller is employed in the treatment. But, this requires a supercomputer application. In turn, each supercomputer has its own operating  system and its own programming peculiarities. Therefore, the program debugging cost is added to the detail one. Note that a 3D model is too time-consuming and that supercomputer time cost is much greater than the time cost of a personal computer. A simplified 3D simulation may bring essential mistakes as it was shown in this Section hence the initial model has to be considered and the application of supercomputers appears to be necessary. On the other hand, the cost of a big roller is greater but all the design of the technological process can be carried out by means of a personal computer because the axisymmetric model proves to be sufficient. The comparison of these two variants (small roller or large roller) is the task to be resolved at the initial stage of the process design. It is obvious that if the material of which the roller will be made is cheap then large roller's size is preferable. And if the material is of a high price then the opportunity of compound roller may be considered when the roller consists of a work part (material of a high price) and an inside (a cheap material). One more issue to be considered is how to determine a concrete value of roll's radius R. Strictly speaking; a series of three dimensional numerical experiments has to be performed for increasing values of R as R→∞. Taking into account a great cost of such investigation, it may be reasonable to be confined to initial stages of the treatment. An appropriate approximation (that is adequate for a given concrete event) by the axial symmetric model will give the optimal value of radius R. In this case, both the design and the control of the process will be an easy task. If it is not the case (for example if R ≈ r) then serious problems may occur.

Conclusion
A stabilized nodal averaged axisymmetric finite element was developed. The formulation is based on the virtual work principle for large irreversible strain together with a nodal averaged approximation and an approximate formula for a ring of triangular profile. An additional stabilizing stress tensor appears to be necessary for incremental numerical analysis. Simple expressions for nodal force and stiffness matrix permit to obtain fast algorithms and short program codes. The finite element can be exploited successfully for numerical simulation of the most complicated task in solid mechanicsmultiple periodic local deformations as it takes place in the deep rolling process. It becomes possible to analyze a whole detail under working rather than its small three-dimensional part. The necessary condition is that longitudinal roller's radius must be sufficiently large. If this is not the case then the task cannot be considered as the axisymmetric one hence both a mathematical simulation and a process control may be too problematic.
A steel rod under deep rolling treatment is simulated. Profiles of residual stress (its axial component) and the yield stress are obtained. The residual stress is calculated in all the cross section rather than in a thin near-surface layer as this occurs in three dimensional models. The Odkvist's parameter isolines are not parallel to the surface wrought hence the process is not stabilized. A bulged surface forms in front of the roller while this bulging grows in the process of working. Taking into account that deep rolling may require many passes, three dimensional simulations of such processes seem to be ineffective.