Abstract
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.
1. 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 B-bar 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.
2. Large Strain Kinematics
The large strain kinematics is defined in terms of deformation gradient tensor , where and are the coordinates of the material configuration 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 in its elastic and plastic components [23, 24] can be written:
In the original configuration the Green Lagrange strain tensor and its elastic and plastic components are defined aswhere and are, respectively, the right Cauchy Green tensor and its plastic component. Finally, is the identity tensor.
In the current configuration of Almansi strain tensor and its elastic and plastic components are defined aswhere the Finger tensor is defined as and its elastic component by .
It is convenient to introduce the rate of deformation tensor and its elastic and plastic components, and , respectively, that are defined as where denotes the Lie derivative. The Lie derivative of a tensor is defined [25] as , 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 , , and in their elastic and plastic components can be obtained.
3. 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 approaches the unity tensor , and the elastic component on the Finger tensor tends to the metric tensor . 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 :where and are constants of the material.
The Cauchy stress tensor can be obtained from (7) as [26]
Plasticity is taken into account by a classical J2 model:where is the yield function, is the deviatoric Cauchy stress tensor, and is the actual yield stress written in terms of the effective plastic strain whose increment is . The associated flow rule is given byas is the plastic multiplier and is outward normal to the yield function. In addition to that, the loading/unloading or Kuhn-Tucker conditions [21, 27, 28] are given by
The actual yield stress has a power hardening law given bywhere and are material constants and is a parameter obtained by constraining (12) to have an initial yield strength of .
Furthermore, the spatial tensor is defined in the usual manner for a linear elastic isotropic material at the spatial aswith , where is the Krönecker delta. Finally, the continuum tangent elastoplastic tensor is defined aswhere the hardening parameter is .
4. 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 predictor-corrector scheme [9, 21] is used. The following summarizes the most important characteristics of these implementation processes.
4.1. 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 is the residual vector, is the internal force vector, and 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 is the tangent stiffness matrix and is the displacement vector increment. As usual in the finite element computations, the tangent stiffness matrix is split into the material stiffness matrix contribution related to the elastoplastic constitutive behavior and geometric stiffness matrix contribution related to the nonlinear effects of the strain measure:
The displacements and the spatial coordinates are updated to fulfill a convenient convergence criterion.
4.2. Updated Lagrangian Formulation with H1/P0 Finite Element
In order to compute the elastic predictor problem the elastic Finger tensor is updated in a multiplicative way aswhere the incremental deformation gradient tensor 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 towhere the corrector term 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 aswhere 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 at the deformed configuration is written in terms of the so-called material and geometrical stiffness matrices and , respectively, aswhere is the standard matrix of the derivative of the shape functions and is the strain-displacement matrix for large strains derived by linearization of . To compute the contribution of to the tangent stiffness matrix the elastoplastic consistent tangent modulus [30] in the spatial configuration is used.
4.3. 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 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 aswhere 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 byConsequently 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 [31] are used (see Appendix A for details). On this basis, the internal force vector is given bywhere 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 and to are given bywhere , is the improved strain-displacement matrix for large strains [31], and is the strain-displacement matrix for large strains derived by linearization of .
5. Numerical Simulation of the Simple Tension Test
5.1. 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 mm length (initial extensometer length) and a rectangular cross section of mm thickness and mm width.
(a) Geometry of the specimen
(b) Finite element mesh
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 , , and , 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.
Boundary conditions are imposed to satisfy symmetry conditions. Displacements (coordinate ), (coordinate ), and (coordinate ) are fixed at the planes , and , respectively.
To simulate the tensile test, displacements are imposed on the plane (at the coordinates) until it reaches 5 mm (i.e., 20% of final elongation or engineering strain).
The material considered in the numerical analysis is SAE 1045 steel with the following mechanical properties: elastic modulus MPa, Poisson’s ratio , MPa, , and MPa.
5.2. Numerical Results and Discussion
In what follows the results obtained with the two implementation processes presented in previous section are compared.
5.2.1. 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 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.
5.2.2. 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 is taken as a characteristic section of the necking. The initial width is the largest dimension of the cross section and it is in the coordinate direction . The initial thickness 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 , as a function of the engineering strain . The engineering strain is defined as the ratio between the total elongation of the specimen () and the initial length ().
The results obtained numerically in this work and experimentally by Cabezas and Celentano [11] 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 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 in terms of engineering strain 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 and . 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 . These deformed configurations correspond to a fixed displacement of 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.
5.2.3. 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 ().
(a) TLF
(b) ULF
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.
5.2.4. 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 () 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.
5.2.5. 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, and , 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.
6. 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 three-dimensional 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.
The effective plastic strain distributions obtained with both implementation processes do not lead to significant differences in the results. Particularly in the cross section of the specimen in the central area of the neck (plane of symmetry ), the contours of effective plastic strain are similar in shape but have slightly different maximum values as a consequence of the slightly different geometric changes in the direction of smaller thickness of the section.
Both axial stress components analyzed, the longitudinal and transverse ones, have quite similar responses, with maximum values in the central area and minimum values in the width of the cross section where necking occurs. Compression values in the width of the specimen are obtained for the transverse stress component.
Pressure and deviatoric stress distributions have quite similar behaviors for both implementation processes. Pressure distributions have larger gradients along the -direction in the PS line than in the FS line. Axial deviatoric stress and transverse deviatoric ones have nearly constant and linear behavior, respectively. It can be noted that pressure and the axial deviatoric stress distributions can be approximated for the analytical expressions of cylindrical specimens.
In summary, the simple tension test in large strain range, in which the necking phenomenon appears, is an appropriate alternative to calibrate the constitutive plasticity J2 models in both ULF and TLF implementation processes. The J2 constitutive model implemented in such formulations was found to provide a realistic response during the tensile test of rectangular bars. Finally, the results for the global and local variables obtained with the two implementation processes exhibit a good agreement.
Appendix
A. B-Bar Finite Element and Matrix
By means of an improved strain-displacement matrix [31] the numerical locking is avoided in the TLF implementation. Based on the deformation gradient multiplicative standard decomposition into deviatoric and volumetric parts and assuming a selective (or averaged or smoothed) numerical integration for the volumetric part of , the derived strain-displacement matrix, directly obtained by linearization of the Green-Lagrange strain tensor, can be classified within the so-called B-bar methods.
In 3D problems the expression of for a given node iswhere with such that is deviatoric part of , is volumetric part of numerically integrated in a selective (of averaged or smoothed) form, is standard strain-displacement matrix for large strain analysis (see [20]), and is improved contribution.
For 3D case iswhere are components of the right Cauchy Green deformation tensor and is as follows: as , and are components of the nodal displacement vector .
B. Pressure and Axial Deviatoric Stress for Cylindrical Specimens
The pressure in a cylindrical specimen can be defined aswhere 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
Considering that the J2 yield criterion for metals in the neck section can be expressed as and substituting the above equation in (B.2), the pressure gives
Finally, if in (B.4) is decomposed into spherical and deviatoric parts, the axial deviatoric stress can be expressed as
Competing Interests
The authors declare that they have no competing interests.
Acknowledgments
The financial support provided by SeCTyP-Universidad Nacional de Cuyo projects 06B/253, 06B/308, and B008 is gratefully acknowledged. Diego Celentano acknowledges the financial support provided by FONDECYT (Project no. 1130404).