An axisymmetric nodal averaged finite element

A nodal averaging technique which was earlier used for plane strain and three-dimensional problems is extended to include the axisymmetric one. Based on the virtual work principle, an expression for nodal force is found. In turn, a nodal force variation yields a stiffness matrix that proves to be non-symmetrical. But, cumbersome non-symmetrical terms can be rejected without the loss of Newton-Raphson iterations convergence. An approximate formula of volume for a ring of triangular profile is exploited in order to simplify program codes and also to accelerate calculations. The proposed finite element is intended primarily for quasistatic problems and large irreversible strain i.e. for metal forming analysis. As a test problem, deep rolling of a steel rod is studied.


INTRODUCTION
The finite element method (FEM) is widely used in science and engineering for numerical analysis of problems which beggar analytical methods.The specific character of solid mechanics tasks is that (except for porous media) a body shape can be varied greatly while body's volume alters insignificantly.But, the approximating functions proposed by classical FEM (see for example the textbook Zienkiewicz (1971)) do not comply with the incompressibility requirement.At first sight, this does not induce any drawback because the convergence theorems have been proven for many FEM calculation schemes.Nevertheless, in practice (mainly when a personal computer is used) an approximate solution can be found to be too far from an exact one and such tasks were really discussed in last decades.This circumstance forced to revise all the classical approximations of tensor fields.Many new approaches have been developed but they require a special review.Plane strain, three-dimensional, and shell problems have been involved but the axisymmetric one remains out of consideration so far and the present work aim is elimination of this gap.
The first attempt to overcome the volumetric locking difficulties was the application of selective reduced integration schemes to displacement-based nonlinear finite elements (bi-linear, tri-linear, second order, and third order).The deviatoric items of a variational equation are integrated by larger number of Gauss points P1,…,Pn while the volumetric items -by less number of points Q1,…,Qm (m<n).In general, it is impossible to meet the incompressibility requirements at all the points Pi but it is possible to do it at points Qi.This method suffers from many shortcomings: sizable computer memory is needed to store the state parameters such as stress, strain, damage, etc. at Gauss points; problems appear when transferring data from one finite element mesh to another during the remeshing procedure; and so on.Nevertheless, the reduced-integration technique came into further development in Reese (2002Reese ( , 2005) ) but is not popular now.
Mixed formulations (see for example Simo et al. (1985)) are numerous and widespread.They are employed in many commercial software products: ANSYS, NASTRAN, ABAQUS, SUPERFORGE, etc. Degrees of freedom comprise both nodal coordinates and nodal stresses.As a variant, only the hydrostatic component of the stress tensor is used in Sussman and Bathe (1987).A finite element matrix appears to be sign-indefinite and may be ill-conditioned such that some preconditioner may be required.The convergence analysis of mixed formulations is presented in monograph by Brenner and Scott (2008).The comparative analysis of a mixed constant pressure hexahedral element Q1P0 and a stabilized nodally integrated tetrahedral by Puso and Solberg (2006) had shown that the results were identical in some tests while the benchmark "Cook membrane" demonstrated the advantage of element by Puso and Solberg.The Q1P0 element converges slowly in this problem due to poor bending performance.
The mixed methods came into further development.A second order stabilized Petrov-Galerkin finite element framework was introduced by Lee et al. (2014).This formulation, written as a system of conservation laws for the momentum and the deformation gradient, yields equal order of convergence for displacements and stresses.In Gil et al. (2014), the formulation is enhanced for nearly and truly incompressible materials with some novelties.Unfortunately, that work is restricted by hyper-elastic scenarios.
Enhanced-assumed strain approach has produced many treatments -see for example Simo and Armero (1992), Simo et al. (1993), Puso (2000), and Areias et al. (2003).An attractive variant was presented by Broccardo et al. (2009) where the assumed strain technique is combined with the nodal average strain formulation but only elastic and hyper-elastic scenarios were under consideration.
Nodal integration, or nodal averaging, yields a particle-type method where stress and material history are located exclusively at the nodes and can be employed when using meshless or finite element shape functions.This particle feature is desirable for large deformation settings because it avoids the remapping or advection of the state parameters required in other methods.The nodal integration relies on fewer stress point evaluation than most other methods (Puso et al. (2008)).Three classical nodal integration schemes are known: nodal strain method by Beissel and Belytschko (1996), stabilized conforming nodal integration by Chen et al. (2001), and nodal averaging by Dohrmann et al. (2000).The enhanced versions of these schemes are given in Puso et al. (2008) where the stability of nodal integration for both meshless and finite element approximations was investigated.
Meshfree, or meshless, methods do not employ a finite element base functions therefore these methods do not suffer from the mesh distortion.More complex (in comparison to finite element) meshfree base functions provide with smooth solutions but computational cost may appear to be too high.A comparative study of the three meshless approaches was carried out in Quak et al. (2011) and a comprehensive review is given by Cueto and Chinesta (2015).
Taking into account all the aforecited, one can conclude that there exist many ways to design an axial symmetric finite element involving super-large strain and high gradients of tensor fields.Therefore, to make an appropriate choice, some special requirements are needed.In this work, the following requirements are determinative: 1) the volumetric locking is damped; 2) the locking due to bending is damped; 3) spurious modes are absent; 4) a finite element passes the "patch test"; 5) a theoretical framework is simple; 6) a program code is short.
A finite element by Puso and Solberg (2006) meet all the requirements and its formulation (based on the nodal integration technique -see above) will be extended on the axial symmetric event here.Unfortunately, two shortcomings are inherent to the formulation.First, mesh distortion causes bad performance and a remeshing procedure together with a data transfer procedure (from an old mesh to a new one) are needed is such case.This problem is studied and resolved in Zhang and Dolbow (2017) specially for the nodal averaged stabilized element by Puso and Solberg (2006).Second, despite exhibiting very good behavior in terms of displacements the formulation tends to exhibit non-physical hydrostatic pressure fluctuations -see Pires et al. (2004) and Gee et al. (2009).In order to overcome this drawback, a smoothing procedure is applied at the final stage of numerical analysis -see Puso and Solberg (2006).Smoothed pressures are evaluated by averaging nodal pressures at element centers and then reaveraging these element pressures back out to the nodes.If it is necessary, this procedure is repeated several times in order to smooth the isoline pattern.So, both the shortcomings of the nodal averaged stabilized element can be rectified.
An approximation based on the simplest finite element (the linear triangle) and the nodal averaging procedure over triangles of a nodal complex is proposed in this work.Due to many reasons, finite elements of low order prove to be preferable for non-linear problems and large deformation.Besides, mesh generators provides 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. (2000) and further analyzed by Bonet et al. (2001).As it turns out, nodal averaging provides weakened constraints such that coarse mesh accuracy can often be achieved for near incompressible materials.It was pointed out also that the formulation was prone to spurious lowenergy modes.In order to suppress them, Puso and Solberg (2006) proposed to enter stabilizing terms into a discrete form of 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. (2000) alleviates both the volumetric locking and locking due to bending.However, some modern methods give the same result as in Puso and Solberg (2006) -see for example Puso et al. (2008).But, stabilized nodal averaging is the simplest among them from the point of view of both a theoretical framework and a program code.The present work extends that formulation on the axisymmetric problem.
This work has one more aspect associated with an impotent industrial application -numerical simulation of the deep rolling working.Such the treatment creates a near-surface hardened layer and compressive residual stress preventing from crack generation in metallic parts.Traditionally, this simulation is performed in three dimensional setting.But, in some events it is possible to simplify a model to an axisymmetric one and to reduce the computational cost drastically -see Section 3.

FINITE ELEMENT FORMULATION
Write the discretized variational equation taking into account (1): where the first sum corresponds with nodal complexes while the second one -with triangles; f is an external load; C is the isotropic elasticity tensor and C  is the stabilization tensor (on optimal choice of C  see Puso and Solberg (2006)).
It is clear that the influence of the stabilization tensor drops as a finite element mesh is refined but for a coarse mesh tensor C  becomes to be necessary.Formally, one can consider that the tessellation 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.Turn to the axisymmetric formulation itself.In cylindrical coordinates = r y 1 , = h y 2 , , the Cauchy stress tensor and the small strain one take the form (here, the small displacement field It is necessary to rewrite the variational principle (2) taking into account large irreversible (plastic) strain and the cpecific form (3) of tensors σ and ε (the Euler approach is adopted throughout the work i.e. all tensor fields refer to the current configuration rather than the initial one): where -average strain at node n A ; -approximate volume for a ring of triangular profile m T ; -average radius of a ring of triangular profile m T .Equation (4) has to be supplemented by a constitutive equation.As far as metals and alloys refer to near rigidplastic 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 , and circular where the summation over repeated indexes is implied ; the differentiation sign is separated by a comma; ) , (   is calculated according to where vector Due to the approximate formula for volume m v , expression (9) yields is the well-known in FEM matrix of coordinate differences of triangle vertices (three nodes are numerated by 1,2,3 counter clockwise and index m is omitted): As for expression (11), taking into account the one-point integration scheme for the classical linear triangle element, function T where this function equals to Now, return to variational principle (4).Substituting the displacement gradient approximation (8) and the approximation of circular deformation (10), find a nodal force at node p A : 6/11 where ij  is Kronecker's delta-symbol.Substituting ( 12) and ( 14) into ( 15) and taking into account Point out that this expression is simpler than (15) 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 ( 4)).Finally, a stiffness matrix is obtained by means of the nodal force variation.To this goal, find the stress tensor variation from the constitutive equation ( 5) with Yaumann's rate: is the vortex tensor.According to ( 17), the nodal stress variation at node n A is where n denotes "no summation" over index n.Besides, matricies n b (see ( 12), ( 13)), matrix b ~ (see ( 14) with comments), and scalars m r , m s are also dependent variables and subjected to variation.Therefore, the total expression for the nodal force ( 16) variation is too cumbersome and consists of two items: symmetrical and nonsymmetrical.In practice of metal forming analysis, a non-symmetrical item can be rejected without the loss of convergence.A symmetrical positively defined item of the stiffness matrix of the nodal finite element is obtained by the first item in the right hand of (19) and is of the form (p and q are node numbers of complex n C ) while if 1 i  or 1 k  then some special items should be added to the stiffness matrix respectively: and if 1 i k   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 in Zienkiewicz (1971) and constitutive equation ( 18)).This completes the nodal averaged axisymmetric finite element formulation.

NUMERICAL EXAMPLES
The following two example problems demonstrate the abilities of the new finite element.The first test is a benchmark involving large strain, elastoplasticity, and specific "stone wall" boundary conditions.The second example shows a possible industrial application when a three dimensional model is reduced to an axisymmetric one with considerable cost saving.As for the stabilizing tensor C  , the natural choice C C  is not optimal.Better results may be obtained by letting (see Puso and Solberg (2006)) where p is a small penalty parameter (p=0.05 is recommended) and Lame's constants   and   are defined in the following way.For elastic material with Lame's constants  and  : For plastic materials, a shear modulus based on the plastic tangent modulus ET is used: A contact problem statement conforms to Morrev (2007Morrev ( , 2009Morrev ( , and 2011)).Only the symmetrical positively defined items ( 20), ( 21), ( 22) of the finite element matrix are evaluated.In this example, which is adopted in Puso and Solberg (2006), performance on a large deformation, confined plasticity problem is evaluated.The undeformed steel (E =2.987×10 4 ksi, ν=0.3) billet is 2 in in height and diameter and quarter symmetry triangular mesh model is used for simulation as seen in Figure 3.The billet is squashed by two parallel rigid plates to thickness 0.5 in.A contact surface of type "stone wall" (i.e. with infinite frictional coefficient) is used to vertically constrain the outer boundary as it expands (bulges) horizontally.The material is elastoplastic with the power law hardening of the form

Cylindrical billet upsetting
In this example, the tangent modulus ET , which is used in stabilizing tensor C  , was the tangent of the power law equation evaluated at εp=0: The only difference from the original work by Puso and Solberg (2006) in the problem statement is the optimal penalty parameter value p=0.01 instead of p=0.05.That is, the axisymmetric finite element and the tetrahedral one are not identical in spite of identical meshes and common theoretical framework!The results are presented in Figure 3 where the smoothing procedure (see Section 1) was applied twice to both the isoline patterns.It is difficult to compare these isolines with the results by Puso and Solberg (2006) because in that work color pictures were used to visualize the tensor fields.Nevertheless, the qualitative similarity is evident.

Deep rolling of a steel rod
A surface hardening process of a steel rod by means of deep rolling (DR) is studied.DR is an extremely complex process (see for example Radchenko et al. (2015) and Gryadunov et al. (2015)) because of high gradient values of tensor fields and due to the induced non-uniformity of a material.As a rule, such problems are solved in a threedimensional statement (despite a workpiece looks like an axisymmetric body) when a small representative part of a detail is considered together with some reasonable boundary conditions.In contrast to this approach, the proposed axisymmetric finite element permits to consider the entire detail in some events.The point is that a shape of a contact zone depends on ratio R/r where R is longitudinal roller's radius and r is transversal one.If r R  then the contact zone is, roughly speaking, round and the material flows in all directions from under the roller.But, if r R  then the contact surface is of the form of a very oblong oval hence the material flow is negligible in the direction of rolling and predominates in the direction of the workpiece axis, i.e. the flow character is the same as in axial sym-   The operation shown in Figure 4 is not realistic but this process possesses the bilateral symmetry therefore the total number of nodes may be reduced twice.In practice, the working is carried out by a single roller or by twothree rollers rotating in a common plane.The finite element mesh is presented in Figure 5 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; the total number of nodes is 2349.The feature of this example is too high gradient values of tensor field (due to local character of deformation) such that a great value 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 penetration depth of plastic deformation is larger too when r R  .It should be pointed out that the threedimensional model in Majzoobi et al. (2010) was calculated up to depth 3 mm only and in the present work an extrapolation of those results up to depth 4 mm is used.The developed finite element permits to evaluate the residual stress without any restriction in depth and a proper result is presented in Figure 7.A stabilized nodal averaged axisymmetric finite element was developed.It 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 solving the most complicated task of solid mechanics -multiple periodic local deformations as it takes place in the DR process.The axisymmetric finite element and the tetrahedral one are not completely identical in spite of identical meshes and common theoretical framework.In some cases (when the restriction in radii r R  holds), it becomes possible to calculate an entire Latin American Journal of Solids and Structures, 2018, 15(2), e14 10/11 detail under working rather than its small three-dimensional part with considerable reduction in computation time and input data preparing.
plane region is given.Let a nodal complex n C be a set of triangles sharing node n A .Denote by   n a set of triangle numbers in complex n C ; by   np -a set of triangle numbers in intersection of complexes n C and p C (i.e.] [ m -a set of node numbers of triangular m T .Figure 1 illustrates these denotations.
in the last but one sum runs over the intersection of nodal complexes p C and n C (this intersection can be empty, consists of one or two triangles or coincides with the entire nodal complex if p at the center of triangle m

Figure 3 :
Figure 3: Upsetting of cylindrical billet.Effective plastic strain is shown for three different thicknesses: (a) 2 in; (b) 0.82 in; and (c) 0.5 in.
metric models.Consequently, a three dimensional model can be reduced to axial symmetric one if radius R is sufficiently large.The computation time diminishes drastically in such case.Turn to an appropriate 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 4).Note once again that if r R  then the problem cannot be considered as an axial symmetric one.Rolling depth 15 .0  h mm; rolling step d=0.5 mm; number of passes is equal to 1. Material is steel AISI 4340; hardening law is

Figure 5 :
Figure 5: Finite element mesh ), (24), (25) appears to be necessary instead of Latin American Journal of Solids and Structures, 2018, 15(2parameter value p=0.05 instead of p=0.01 in the previous example.The residual stress component along the rod axis is to be found.A like problem was investigated inMajzoobi et al. (2010) with the only difference: longitudinal radius was equal results of both the works are compared in Figure6.The peak values are almost equal.

Figure 6 :
Figure 6: Axial residual stress up to deph 4mm

Figure 7 :
Figure 7: Axial residual stress along the entire deph