Global and Local Mechanical Responses for Necking of Rectangular Bars Using Updated and Total Lagrangian Finite Element Formulations

In simulations of forged and stamping processes using the finite element method, load displacement paths and three-dimensional stress and strains states should be well and reliably represented.The simple tension test is a suitable and economical tool to calibrate constitutive equations with finite strains and plasticity for those simulations. A complex three-dimensional stress and strain states are developed when this test is done on rectangular bars and the necking phenomenon appears. In this work, global and local numerical results of the mechanical response of rectangular bars subjected to simple tension test obtained from two different finite element formulations are compared and discussed. To this end, Updated and Total Lagrangian formulations are used in order to get the three-dimensional stress and strain states. Geometric changes together with strain and stress distributions at the cross section where necking occurs are assessed. In particular, a detailed analysis of the effective plastic strain, stress components in axial and transverse directions and pressure, and deviatoric stress components is presented. Specific numerical results are also validated with experimental measurements comparing, in turn, the performance of the two numerical approaches used in this study.


Introduction
Currently there are many technological applications where computational modeling with the finite element method is extremely useful when, in particular, the problem studied leads to complex triaxial stress states.Problems such as metal forming, impact, and behavior of energy dissipation devices, among others, are typical examples in which it is necessary to consider large strains and elastoplasticity in the finite element models.In all of these cases, different finite element types and formulations are used.However, it is not usual to assess their predicting capabilities of the global and local mechanical responses of the problem.
Calibration of the constitutive equations with experimental results in the large strain range is a requisite for reliable finite element modelling.For this purpose, the simple tension test is one of the most appropriate ways in practice to characterize the mechanical response of metals at large plastic strain ranges and triaxial stress states.
The simple tension test has a first stage where small or moderate strains are observed.In this case, the material behavior is practically linear.In a second stage, strains increase with little changes on the applied load.Finally, after the maximum load is reached, strains are concentrated in a zone where large geometric changes produce the necking in the specimen.In the necking zone, large plastic strains and a triaxial state of stresses can be found.
Bridgman [1] and Davidenkov and Spiridonova [2] proposed analytical expressions in order to describe the relationships between stress components and yield stress at the necking section.Early numerical simulations of the simple tension test are due to Wilkins [3], Chen [4], and Needleman [5], among others.Norris et al. [6] and Goicolea [7] compared numerical and experimental results for steel and aluminum specimens, respectively.Later, Hallquist [8], Simó [9], Ponthot [10], Cabezas and Celentano [11], and García Garino et al. [12] discussed this problem as well.From an application's point of view, Safari et al. [13], Paul [14], and Kamaya and Kawakubo [15] studied the flow of equivalent stress in terms of effective plastic strains using the simple tension test and different experimental procedures.Moreover, the works by Wang et al. [16], Coppieters and Kuwabara [17], Kim et al. [18], and Yazzie et al. [19] are dedicated to characterizing the global postnecking behavior in rectangular bars.It is important to remark that all of these works are mainly devoted to predicting the global response of the material during the tensile test.Therefore, the analysis of the local response is a relevant aspect that needs further research.
Both the Updated and Total Lagrangian formulations (see Bathe's textbook [20] and hereafter, resp., denoted as ULF and TLF) have been used for necking simulation in the tensile test.In practice, ULF has been preferred because the constitutive laws obtained from experimental results are usually written in terms of spatial variables [9,10,12].However, Cabezas and Celentano [11] successfully simulated the necking problem using a TLF based approach.Although there is a lot of information in the literature regarding ULF and TLF, it should be noticed that, up to the authors' knowledge, a comprehensive assessment of both formulations in problems involving metals subjected to large plastic strains has not been previously addressed.
The aim of this paper is to present a detailed study of the global and local mechanical responses of rectangular samples subjected to simple tension conducted in order to obtain strain and stress distributions at the minimum cross section.In particular, effective plastic strain, stress components in axial and transverse directions, and pressure and deviatoric stress components are presented and discussed.To this end, finite element simulations are carried out using the ULF approach with H1/P0 elements [21,22] and the TLF with Bbar elements [11] considering in both cases the same constitutive model.Thus, the detailed analysis of the strain and stress states at the necking section of rectangular bars using these two numerical approaches is the main original contribution of the present work.
The large strain kinematics and constitutive model are provided in Sections 2 and 3, respectively.The finite element implementation of the ULF and the TLF is presented in Section 4. In Section 5, the numerical simulation of the simple tension test using a SAE1045 steel rectangular sample is studied in order to compare the results obtained with both formulations.Global and local characteristic variables (such as applied load in terms of true strain and geometric changes at the necking section) obtained with such approaches considering the same finite element mesh are analyzed in detail.Furthermore, distributions of effective plastic strain, stresses in the longitudinal and transverse directions, and pressure and deviatoric stresses are discussed at the loading stages showing large plastic strains and a triaxial stress state.Finally, Section 6 summarizes the main conclusions of this work.

Large Strain Kinematics
The large strain kinematics is defined in terms of deformation gradient tensor F = x/X, where X and x are the coordinates of the material configuration 0 Ω and the current configuration  Ω, respectively, in a Cartesian coordinates system, as can be seen in Figure 1.
For the elastoplastic case, the multiplicative decomposition of the deformation gradient tensor F in its elastic F  and plastic F  components [23,24] can be written: In the original configuration 0 Ω the Green Lagrange strain tensor E and its elastic E  and plastic E  components are defined as where C = F  F and C  = F  F  are, respectively, the right Cauchy Green tensor and its plastic component.Finally, 1 is the identity tensor.
In the current configuration of Almansi strain tensor e and its elastic e  and plastic e  components are defined as where the Finger tensor is defined as b It is convenient to introduce the rate of deformation tensor d and its elastic and plastic components, d  and d  , respectively, that are defined as d =  V (e) ; d  =  V (e  ) ; where  V denotes the Lie derivative.The Lie derivative of a tensor z is defined [25] as  V (z) =  * [(/)( * z)], where  * and  * are the push-forward and pull-back operators [25], respectively.It is immediately recognized that ( It is important to point out that, from the multiplicative decomposition, additive decompositions for E, e, and d in their elastic and plastic components can be obtained.

Constitutive Model
The constitutive law is written in the context of irreversible thermodynamics of solids by using the free energy , decomposed additively in its elastic   and plastic components   [21]: being  a set of internal variables.For metals under large strains, the elastic strains are negligible, even for the case of large plastic deformations.The elastic component of the gradient tensor F  approaches the unity tensor 1, and the elastic component on the Finger tensor b −1 tends to the metric tensor g.In this case, the distinction between intermediate and current configurations has no meaning.Then it is possible to write the elastic component of free energy function as a quadratic function of elastic component of Almansi strain tensor e  : where  and  are constants of the material.The Cauchy stress tensor can be obtained from (7) as [26]  =   (e  ) e  = tr (e  ) 1 Plasticity is taken into account by a classical J2 model: where (s,   ) is the yield function, s is the deviatoric Cauchy stress tensor, and  is the actual yield stress written in terms of the effective plastic strain   whose increment is [(2/3)d  : d  ] 1/2 .The associated flow rule is given by as γ is the plastic multiplier and n is outward normal to the yield function.In addition to that, the loading/unloading or Kuhn-Tucker conditions [21,27,28] are given by γ ≥ 0; The actual yield stress has a power hardening law given by where  and  are material constants and  is a parameter obtained by constraining (12) to have an initial yield strength of  0 .Furthermore, the spatial tensor a  is defined in the usual manner for a linear elastic isotropic material at the spatial as with where  is the Krönecker delta.Finally, the continuum tangent elastoplastic tensor a  is defined as where the hardening parameter is  = /  .

Finite Element Implementation
The finite strain constitutive elastoplastic model presented previously in Section 3 has been implemented in both Total and Updated Lagrangian formulations.In the first case H1/P0 finite element is used while a so-called B-bar finite element is considered for the Total Lagrangian formulation.In order to integrate the constitutive model, in both cases a predictorcorrector scheme [9,21] is used.The following summarizes the most important characteristics of these implementation processes.

Incremental Iterative Equilibrium Solution.
The global discretized equilibrium equations, for load level at the present quasistatic case, can be written in matrix form for a certain time  as [20] where R is the residual vector, F int is the internal force vector, and F ext denotes the incremental external force vector.For a given load increment this equilibrium equation is solved at time  + Δ by an incremental iterative Newton-Raphson scheme: where K  is the tangent stiffness matrix and Δu is the displacement vector increment.As usual in the finite element computations, the tangent stiffness matrix is split into the material stiffness matrix K  contribution related to the elastoplastic constitutive behavior and geometric stiffness matrix K  contribution related to the nonlinear effects of the strain measure: The displacements and the spatial coordinates are updated to fulfill a convenient convergence criterion.

Updated Lagrangian Formulation with H1/P0 Finite Element.
In order to compute the elastic predictor problem the elastic Finger tensor b −1 is updated in a multiplicative way as where the incremental deformation gradient tensor f is given by The elastic part of the Finger tensor plays the role of an internal variable in this model.Details of the derivation of ( 18) can be found in the work of García Garino [21].
If necessary the plastic corrector is based on a backward Euler integration scheme, leading to where the corrector term  +Δ n is obtained by the radial return algorithm [9].Then we obtain the elastic Almansi strain tensor from Then, substituting (20) in above equation gives Replacing ( 22) in ( 8), the following expression of Cauchy stress tensor is obtained: in which the trial stress is given by To avoid numerical locking due to incompressible plastic deformation, an H1/P0 mixed finite element is used.The H1/P0 element is the extension to the three-dimensional case of the well known Q1/P0 element proposed by Nagtegaal et al. [29].The key ingredient for this element is the smoothing of the pressure field averaging its value at the element level.
The internal force vector, in the current configuration  Ω, is computed as where  is the mean pressure and   is the volume in the current configuration.The incremental external force vector is also computed in  Ω.
The tangent stiffness matrix  K  at the deformed configuration  Ω is written in terms of the so-called material and geometrical stiffness matrices  K  and  K  , respectively, as where B is the standard matrix of the derivative of the shape functions /  x and G is the strain-displacement matrix for large strains derived by linearization of B. To compute the contribution of K  to the tangent stiffness matrix the elastoplastic consistent tangent modulus [30] in the spatial configuration is used.

Total Lagrangian Formulation with B-Bar Finite Element.
In the elastic predictor problem the trial elastic Green Lagrange tensor is computed first as Then the trial Second Piola Kirchhoff stress [20] results: where +Δ A  is the fourth-order elasticity tensor written in the original configuration and  * is the pull-back operator [25].Finally the trial Cauchy stress tensor is computed as where  * is the push-forward operator [25].
If necessary in the plastic, corrector problem, the plastic multiplier  is computed using a radial algorithm in a similar way to the one discussed for ULF formulation in Section 4.2.Then the increment of the plastic Green Lagrange tensor is given by Consequently the elastic Green Lagrange tensor is computed as Finally, the Second Piola Kirchhoff stress tensor results: In this case the internal forces vector is computed in the original configuration.In order to avoid numerical locking due to incompressible plastic deformation, a B-bar finite element and an improved strain-displacement matrix B [31] are used (see Appendix A for details).On this basis, the internal force vector is given by where 0  is the volume in the material configuration.Furthermore, in the TLF the incremental external force vector is computed in the original configuration.
Finally, the contributions of  K  and  K  to  K  are given by where A  =  * a  , B is the improved strain-displacement matrix for large strains [31], and H is the strain-displacement matrix for large strains derived by linearization of B.

Finite Element Mesh and Material
Properties.The geometry of the rectangular bar used in the numerical simulation of the simple tension test is shown in Figure 2(a).The initial geometry has 2 0 = 50mm length (initial extensometer length) and a rectangular cross section of 2 0 = 6mm thickness and 2 0 = 12.5 mm width.Due to the symmetry of the specimen only one-eighth of it has been modeled.Thus the length, thickness, and width in the finite element model are  0 ,  0 , and  0 , respectively.The 3D finite element mesh is shown in Figure 2(b).The total number of the mesh elements is 3440 hexahedral elements with 4266 nodes.A refined mesh in the central section is used to correctly describe the deformation gradients expected in the necking zone.
To ensure location of necking at the central part of the specimen, a geometric imperfection is imposed as a linear width variation along its length.The width is reduced to 1.376% in the central area of the specimen.
To simulate the tensile test, displacements   are imposed on the  plane (at the  = 25 coordinates) until it reaches 5 mm (i.e., 20% of final elongation or engineering strain).

Numerical Results and Discussion.
In what follows the results obtained with the two implementation processes presented in previous section are compared.

Evolution of Axial
Load.The evolution of the tensile force  applied to the specimen as a function of the logarithmic strain in the neck obtained with both implementation processes is shown in Figure 3.The logarithmic strain or true strain is obtained as the natural logarithm of the ratio of the initial area  0 of the neck section of the specimen and the area  of the same section, for each value of applied load.The numerical results are in good agreement with the experimental ones reported by Cabezas and Celentano [11].
An overall good agreement between the two finite element implementations can be observed in these curves.In both cases is obtained a linear behavior for loads of less than 33950 N, which is the initial yield load, presenting deformations very close to zero.On the other hand, a nonlinear relationship between deformation and load is obtained for greater values than 33950 N.
We can state that both implementation processes have a behavior very close to each other for the global parameter behavior of the tensile test given by applied load-true strain ratio.

Evolution of Necking.
As characteristic parameter of the necking phenomenon it is possible to consider the evolution of the specimen cross-sectional dimensions in the necking zone.The symmetry plane  = 0 is taken as a characteristic section of the necking.The initial width  0 is the largest dimension of the cross section and it is in the coordinate direction .The initial thickness  0 is the smallest dimension of the cross section and it is in the coordinate direction .
Figure 4 shows the evolution of the ratio between the current width  of the cross section of the specimen in the necking zone and the initial width  0 , as a function of the engineering strain Δ/ 0 .The engineering strain is defined as the ratio between the total elongation of the specimen (Δ =  −  0 ) and the initial length ( 0 ).
The results obtained numerically in this work and experimentally by Cabezas and Celentano [11]  experimental values.It may be noted that both implementation processes deform more than the experimental reported values.
The evolution of the ratio between the current thickness  of the cross section of the specimen (in the necking zone) and initial thickness  0 in terms of engineering strain Δ/ 0 is plotted in Figure 5.The experimental results obtained by Cabezas and Celentano [11] are also plotted.It can be observed that both implementation processes generally conform well to the experimental results.Up to engineering strain values of 0.12 the agreement is excellent, while for values greater than 0.12 both implementation processes deform more than the reported values experimentally.It may be noted that TLF deforms slightly more than ULF.
From comparing Figures 4 and 5 for the same engineering strain value greater than 0.12, it can be stated that the reduction in thickness is greater than the reduction in width of the cross section of the specimen in relation to their initial values  0 and  0 .That is, the smaller dimension of the cross section has a greater reduction in size.In Figure 6 deformed configurations in the central cross section of the specimen are shown for the two implementation processes with the original configuration.The central cross section is placed at plane of symmetry  = 0.These deformed configurations correspond to a fixed displacement of   = 5 mm, that is, 20% of final elongation.It should be noted that the width and thickness reductions are consistent with those shown in Figures 4 and  5.
It can be seen in Figure 6 that with TLF larger necking is obtained than with Updated Lagrangian one.Moreover, for both implementation processes it appears that reducing the thickness  (coordinate direction ) is greater than the reduction of the width  (coordinate direction ).
The two implementation processes represent very well the evolution of necking in the cross section of the specimen corresponding to the plane of symmetry.It should be noted that the cross section in the necking zone has a greater geometric change in the direction of smaller dimension.

Contours of Effective Plastic Strain.
The effective plastic strain contours at the end of the simulation, corresponding to the fracture stage for elongation of 20.0%, are shown in Figures 7(a) and 7(b) for the TLF and the ULF formulations, respectively.As can be seen, the global distributions of the effective plastic strain are very similar between the two implementation processes, although the maximum values are slightly different.These maximum values are placed in the transverse plane of symmetry ( = 0).The maximum effective plastic strain value occurs in the central point of the cross section of symmetry of the specimen, where maximum necking is developed.The effective plastic strain obtained with TLF is 1.027.It is slightly higher than that obtained with ULF which is 0.906.

Effective Plastic Strain Distribution in the Cross
Section of the Neck.Since the maximum effective plastic strains are slightly different it is of interest to analyze the distribution of these in the central section of the necking.This is done through analyzing the strain values on Gauss points of the used finite element mesh.The location of these Gauss points at the  section close to the transversal plane of symmetry ( = 0) can be seen in Figure 8.
The effective plastic strain level curves plotted at the deformed configuration are shown in Figures 9 and 10 for the TLF and ULF, respectively.It can be seen that in both cases the shapes of the level curves are similar and they approach circumferences dominated by the smaller dimension of the specimen.However, the maximum values are different.
In order to analyze the effective plastic strain development, values greater than 0.9 are compared in Figures 9 and  10.While Figure 9 shows an important zone of the section with values greater than 0.9, in Figure 10 these values are restricted to a setting close to the section center.Thus, TLF developed a largest area of effective plastic strain, greater than 0.9.It seems that the gradient of the effective plastic strain is greater in the TLF results than the ones in the ULF.
Two interesting straight lines are shown in Figure 8.One of them is placed close to the free surface of the specimen and it corresponds to the Gauss points location at undeformed coordinate equal to 2.9032.The other one is placed close to the  plane of symmetry of the specimen and it corresponds to the Gauss points location at the undeformed -coordinate equal to 0.0788.Hereafter, the former is called free surface line (FS line), while the latter is the plane of symmetry line (PS line).
The effective plastic strain in the necking section as a function of the deformed -coordinate is shown in Figure 11 for both implementation processes and for the same load compared to the previous effective plastic strain level curves.The deformed -coordinate is the variable along the width  of the deformed section.As it can be seen in Figure 11, the TLF has higher values than the ULF one along both lines (FS line and the PS line).This is more enhanced from the center of the section (where more plasticity effects have been developed) to the edges.
Although not shown, according to the hardening law given by ( 12), the equivalent stress distributions have the same profiles compared to those of the effective plastic strain shown in the previous figure.

Axial and Transverse Stress Distributions.
It is also important to analyze the behavior of some of the Cauchy stress tensor components.The more representative ones are the axial component   , in the applied load direction, and the one in transverse direction   .Smoothed values for the stress components are presented below.At each element, the smoothed value for each stress component considered is obtained computing the average of the corresponding Gauss point values at FS line and PS line, respectively.
Figure 12 shows the axial stress   (MPa) in the longitudinal direction of the specimen in terms of the deformed -coordinate, evaluated at the FS line and the PS line.The values given by both implementation processes are very close together in both lines, presenting maximum values at center of the section.Higher values are obtained in the PS line than in the FS line.Both implementation processes provide higher values of   (MPa) in the center of the section than in the edges.This corresponds to the plastic flow evolution.
Figure 13 shows the transverse stress   (MPa) in terms of the deformed -coordinate, evaluated in the same points as the previous figure.The   stress corresponds to the transverse axial stress in the smaller direction of the specimen cross section.As in the previous figure, the values for both implementation processes are very close together in both lines (FS Line and PS Line).Both implementation processes approach quite well   equal zero boundary condition for the free edge along the FS line.Similar   variations along the PS line are given by both implementation processes.While at the center of the cross section the   stress has traction values, near the free edge compression values appear.5.2.6.Pressure and Deviatoric Stress Distributions.Once again, smoothed values for the pressure and deviatoric normal stress components are presented below, using the procedure previously described and discussed.As has been suggested by Gabaldón [32], the distributions of these variables at the necking zone are shown in Figures 14 and 15, applied to a rectangular bar in this work.Overall, nearly coincident results are obtained for both implementation processes.In particular, very similar pressure distributions along the PS line and FS line can be appreciated where larger gradients along the -direction are observed in the PS line than in the FS line.Moreover, it is seen that the axial deviatoric stress   distribution is practically constant in the whole section.In addition, the transverse deviatoric stress components   and   present a small linear variation along the PS line and FS line.
In these figures, the analytical approximations for the pressure and the axial deviatoric stress for cylindrical specimens (see Appendix B for details), that is,  =   −2/3 and   = 2/3, are also shown.The pressure results obtained with both implementation processes have similar behavior to those obtained with the analytical approximation of cylindrical specimens.Furthermore, the axial deviatoric stress for rectangular specimen has a constant distribution as in the cylindrical case.

Conclusions
The response of the simple tension test of a rectangular specimen of SAE1045 steel into a postnecking behavior has been obtained using two well known large strain finite element implementation processes.An Updated Lagrangian formulation with a H1/P0 element (ULF implementation) and a Total  Lagrangian formulation with a B-bar element (TLF implementation) both using hexahedral isoparametric elements have been used for the numerical simulations.Although the plastic constitutive model is the same for both approaches, the implementation in each formulation is quite different, and their performances on the prediction of complex threedimensional stress and strain states such as those occurring at the necking section of a tensile sample are not usually discussed in literature.
Global and local mechanical results have been obtained and compared using the two finite element implementation processes.Moreover, some specific numerical results at the necking section have been also validated with experimental measurements.As a global comparison parameter the applied load as a function of the true strain has been used.These curves obtained with both finite element implementation processes are very close to the experimental results.
At the central cross section where the necking is more evident, some local parameters have been evaluated and discussed.Particularly, the deformed shape of the central cross section and local constitutive parameters, such as the effective plastic strain, stress components, pressure, and deviatoric stresses have been compared and analyzed.
Both numerical solutions represent very well the geometric evolution of the necking in the cross section of the specimen corresponding to the plane of symmetry.Very well fitting has been found regarding the changes in thickness  and width of the cross section in terms of engineering strain compared with the experimental results.
It should be noted that the cross section of the specimen in the necking zone presents greater geometric change in the direction of smaller dimensions.Both implementation processes predict in excess these deformations compared with the experimental results.Particularly, the TLF results are a bit higher than the ULF.

B. Pressure and Axial Deviatoric Stress for Cylindrical Specimens
The pressure in a cylindrical specimen can be defined as where   and   are, respectively, the radial and circumferential components of the Cauchy stress tensor.The pressure at a point of the solid for a stress state present in the neck [33] can be expressed as

Figure 3 :
Figure 3: Tensile force as a function of logarithmic strain in the neck.

Figure 4 :
Figure 4: Ratio of current to initial width versus engineering strain.

Figure 6 :
Figure 6: Deformed configurations in the necking zone for a final elongation of 20%.

Figure 7 :
Figure 7: Comparison of effective plastic strain contours for a final elongation of 20%.

Figure 8 :
Figure 8: Gauss points in the undeformed configuration.

Figure 9 :
Figure 9: Level curves of effective plastic strain for TLF and 20% of final elongation.

Figure 10 :
Figure 10: Level curves of effective plastic strain for ULF and 20% of final elongation.

Figure 11 :
Figure 11: Effective plastic strain in the neck section along coordinate axis  for a final elongation of 20%.

Figure 12 :Figure 13 :
Figure 12: Axial stress   (MPa) in the neck section along coordinate axis  for a final elongation of 20%.

Figure 14 :Figure 15 :
Figure 14: Pressure (MPa) in the neck section along coordinate axis  for a final elongation of 20%.

2 )
Considering that the J2 yield criterion for metals in the neck section can be expressed as  −   =  (B.3)and substituting the above equation in (B.2), the pressure gives =   − in (B.4)   is decomposed into spherical and deviatoric parts, the axial deviatoric stress can be expressed as are plotted in Figure 4.It can be observed that both implementation processes, generally, have a good agreement with the experimental results.Up to engineering strain values of 0.11 the agreement is excellent, while, for values greater than 0.11 the adjustment of both implementation processes is good compared to 0 Figure 5: Ratio of current to initial thickness versus engineering strain.