# Modeling and Numerical Computation of Necking in Round Bars Using a Total Lagrangian Elastoplastic Formulation

Anh Le Van, Philippe Le Grognec

► To cite this version:

Anh Le Van, Philippe Le Grognec. Modeling and Numerical Computation of Necking in Round Bars Using a Total Lagrangian Elastoplastic Formulation. Computer Modeling in Engineering and Sciences, 2001, 2 (1), pp.63-72. 10.3970/cmes.2001.002.063. hal-01385743

HAL Id: hal-01385743

https://hal.science/hal-01385743

Submitted on 22 Oct 2016

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Modeling and Numerical Computation of Necking in Round Bars Using a Total Lagrangian Elastoplastic Formulation

A. Le van<fnref n="1" />, P. Le Grognec<fnref n="1" />

**Abstract:** Necking is a bifurcation phenomenon observed in round bars under tensile loading and has been investigated in numbers of papers. In the present work, it is modeled within the framework of finite rate-independent plasticity. The theory is based on thermodynamic foundations developed for standard materials and results in a total Lagrangian formulation for finite plasticity, where the total strain is decomposed additively according to [Green and Nagdhi 1965] and the hardening is characterized by a nonlinear isotropic hardening law of the saturation type.

The discretization and consistent linearization of the elastic-plastic equation set using the standard finite element procedure lead to a low-cost algorithm, robust enough to deal with necking problems.

The numerical computations of necking are performed on cylindrical bars with various boundary condition types and the corresponding results compared with those obtained in the literature.

**keyword:** elastoplasticity, necking bifurcation, total Lagrangian formulation, finite element analysis

# 1 Introduction

Necking is an bifurcation phenomenon observed in circular cylindrical bars under tensile loading and has been investigated in numbers of works for its theoretical modeling and numerical computation. The pioneering papers dated back to the 70's and were due to [Needleman (1972)], [Argyris and Doltsinis (1979)]. A lot of significant results were obtained in these works, although the computed necked configurations did not resemble experimental observations. Later, [Besdo (1987)] was among the first to give a realistic necking shape. More recently, studies using other formulations have been developed

by [Simo (1988)], [Simo (1992)], [Simo and Hughes (1998)] and [Brünig (1998)], providing further results on necking. The numerical solution by Simo necessitates mixed finite elements, moreover the line search method is used to get across the limit point.

Whereas most papers considered a cylindrical circular bar with both ends assumed to be shear free, [Needleman (1972)] dealt with two types of boundary conditions, shear free ends and cemented to rigid grips ones. [Brünig (1998)] went further by also considering the whole tensile specimen with shear free stiff heads and then compared the results corresponding to the three cases.

Most of the previous investigations were performed within the framework of finite rate-independent plasticity, assuming a yield criterion in the stress space and an isotropic hardening. Outside this context, analyses using the strain space representation of the finite plasticity were proposed by [Besdo (1987)] and [Simo (1988)]. Other nontraditional approaches can be used to successfully solve the necking problem, as shown by [Atluri and Murakawa (1980)], using a modified complementary energy principle instead of the conventional principle of virtual work. Also, the integral equation method was shown to be capable of solving bifurcation and post-bifurcation problems, see [Okada, Rajiyah and Atluri (1990)].

On the other hand, thermal aspects in necking were taken into account by [Kleiber (1992)] who presented a simple solution scheme for elastoplasticity as well as thermoplastoplasticity. Also, one can use none of finite plasticity but a viscoplastic type constitutive law in small strains in order to describe necking behavior in viscoplasticity, e.g. see [Rappaz, Bellet and Deville (1998)].

In the present work, necking is studied within the framework of finite rate-independent plasticity. The theory is based on thermodynamic foundations established for generalized standard materials [Halphen and Nguyen (1975)], [Germain, Nguyen and Suquet (1983)] and results in a total Lagrangian formulation for finite plas-

<footnote n="1">¹Laboratoire de Mécanique et Matériaux, Division Mécanique des Structures, Ecole Centrale de Nantes. 1, rue de la Noë, BP 92101 - 44321 Nantes cedex 3, France</footnote>
ticity. The total strain is decomposed additively according to [Green and Naghdi (1965)]. The yield function is of von Mises type and depends on the symmetric Piola-Kirchhoff stress tensor. The hardening in the material is characterized by a nonlinear isotropic hardening law of the saturation type, as in [Simo (1988), Part II] and [Simo (1992)]. However, here it relates Lagrangian (material) variables, namely the Lagrangian hardening variable and the Lagrangian effective plastic strain.

The numerical implementation using standard finite element analysis is simple and yet robust enough to accurately reproduce the necking phenomenon. The local integration of the elastic-plastic equation set is carried out using the fully implicit integration scheme (backward Euler difference scheme), which ensures both the stability and the symmetry of the consistent tangent modulus. The numerical computations of necking are carried out both on shear free end bars including an initial small geometrical defect and geometrically perfect bars with gripped ends. The corresponding results are then discussed in some detail and compared with those obtained in the literature.

## 2 Theoretical background

Let us consider a three-dimensional continuum body which at initial time $t_o$ occupies a region $\Omega_o$ bounded by surface $S_o$ and assume that the body undergoes a quasi-static elastic-plastic transformation. The initial configuration $\Omega_o$ will be chosen as the reference one and the equations of the problem expressed using Lagrangian variables relative to this configuration.

The motion of the body in the Lagrangian description is defined by the displacement field denoted by $\vec{U}$. The deformation inside the body is measured by the Green strain tensor defined in terms of the displacement

$$ \vec{E} = \frac{1}{2} (\overline{\overline{grad}} \vec{U} + \overline{\overline{grad}}^T \vec{U} + \overline{\overline{grad}}^T \vec{U} \cdot \overline{\overline{grad}} \vec{U}) \quad (1) $$

According to [Green and Naghdi (1965)], [Green A. E.; Naghdi P. M. (1971)], the strain tensor $\vec{E}$ is decomposed additively as

$$ \vec{E} = \vec{E}^e + \vec{E}^p \quad (2) $$

where the so-called plastic strain $\vec{E}^p$ is assumed to be symmetric and objective, so that the elastic strain $\vec{E}^e$ is

symmetric too, see [Green and Naghdi (1965)], [Casey and Naghdi (1980)].

In a quasi-static process, the equilibrium equation reads

$$ \mathrm{div}(\vec{F} \cdot \vec{\Sigma}) + \rho_0 \vec{f}_0 = \vec{0} \quad (3) $$

where $\vec{F}$ is the deformation gradient, $\vec{\Sigma}$ the second Piola-Kirchhoff stress tensor (symmetric), $\rho_0$ the reference mass density and $\vec{f}_0$ the body force. The product $\vec{\Pi} = \vec{F} \cdot \vec{\Sigma}$ is the first Piola-Kirchhoff stress tensor, it is not symmetric whereas $\vec{\Pi} \cdot \vec{F}^T$ is.

The total Lagrangian formulation presented herein is based on thermodynamic foundations developed for the so-called generalized standard materials [Halphen and Nguyen (1975)], [Germain, Nguyen and Suquet (1983)]. Within this framework, we consider that the state of any particle and at any time is characterized by the following state variables : the total strain $\vec{E}$, the plastic strain $\vec{E}^p$ and another internal variable to account for the deformation history, represented by the scalar hardening variable denoted by $\alpha$.

To complete the equation set for the time-independent elastic-plastic problem, we shall provide the state laws together with the evolution laws. The state laws give the symmetric stress $\vec{\Sigma}$ and the thermodynamic conjugate $A$, dual of the hardening variable $\alpha$, by differentiating the free energy per unit volume $w$

$$ \vec{\Sigma} = \frac{\partial w}{\partial \vec{E}^e} \qquad A = \frac{\partial w}{\partial \alpha} \quad (4) $$

As usual, the stored-energy $w$ is assumed to be the sum of the elastic energy $w^e$ -function of the elastic strain $\vec{E}^e$ only - and the hardening energy $w^\alpha$ - function of the hardening variable $\alpha$

$$ w(\vec{E}^e, \alpha) = w^e(\vec{E}^e) + w^\alpha(\alpha) \quad (5) $$

For metals, we can assume small elastic strains and choose the elastic energy as a quadratic function of the elastic strain

$$ w^e(\vec{E}^e) = \frac{1}{2} \vec{E}^e : \tilde{D} : \vec{E}^e = \frac{1}{2} E_{ji}^e D_{ijkl} E_{lk}^e \quad (6) $$

where $\tilde{D} = 2\mu\vec{I} + \lambda\vec{J}$ denotes the isotropic elasticity tensor ($\vec{I}$ is the fourth-order identity tensor, $\vec{e}_i = \vec{e}_i \otimes \vec{e}_j \otimes \vec{e}_i$, $\vec{J} = \vec{I} \otimes \vec{I}$, $\vec{I}$ the second order identity tensor). In order to
describe nonlinear hardening effects, the hardening energy is chosen as a function containing exponential terms

$$w^{\alpha}(\alpha) = c_0\alpha + \frac{1}{2}c_2\alpha^2 + (c_1 - c_0)\left(\alpha + \frac{\exp(-c_3\alpha)}{c_3}\right) \quad (7)$$

where the coefficients $c_i$ are material-dependent constants.

The state laws then yield the stress-strain relation and the nonlinear relation between hardening variables

$$\bar{\Sigma} = \bar{D} : \bar{E}^e = \bar{D} : (\bar{E} - \bar{E}^p) \quad (8a)$$

$$A = c_0 + c_2\alpha + (c_1 - c_0)(1 - \exp(-c_3\alpha)) \quad (8b)$$

The elastic and plastic ranges at any particle in the body are defined by the yield function $f$ of von Mises type, depending on the symmetric stress $\bar{\Sigma}$ and the hardening variable $A$

$$f(\bar{\Sigma}, A) = \sqrt{\frac{3}{2} \|\bar{S}\|} - \sigma_0 - A \quad (9)$$

where $\bar{S}$ is the deviator of stress $\bar{\Sigma}$ and $\sigma_0$ the initial yield stress.

The evolution laws result from the generalized maximum plastic dissipation rule and provide the plastic strain and hardening variable rates by differentiating the yield function with respect to its variables

$$\dot{\bar{E}}^p = \dot{\lambda} \frac{\partial f}{\partial \bar{\Sigma}} - \dot{\alpha} = \dot{\lambda} \frac{\partial f}{\partial A} \quad (10)$$

where the plastic multiplier $\dot{\lambda}$ satisfies $\dot{\lambda} \ge 0$ if $f(\bar{\Sigma}, A) = 0$ and $\frac{\partial f}{\partial \bar{\Sigma}} : \bar{\Sigma} \ge 0$, otherwise $\dot{\lambda}=0$. Of course, there is no relation between the Lamé constant $\dot{\lambda}$ and the plastic multiplier $\dot{\lambda}$. Relation (10a) is the plastic flow rule and relation (10b) the hardening law. They identically fulfill the second law of thermodynamics for continua, written in terms of objective Lagrangian quantities.

Since the plastic strain $\dot{\bar{E}}^p$ is a Lagrangian variable, its rate is obtained by ordinary differentiation with respect to kinematic time $t$. There follows from relations (9)-(10) that the hardening rate is equal to the Lagrangian equivalent plastic strain rate

$$\dot{\alpha} = \dot{P} = \sqrt{\frac{2}{3} \dot{\bar{E}}^p : \dot{\bar{E}}^p} \quad (11)$$

Relations (1)-(3), (8)-(10) form a coupled system of 30 scalar partial differential equations with 30 scalar unknowns $\vec{U}, \vec{E}, \vec{E}^e, \vec{E}^p, \bar{\Sigma}, A, \alpha$, and $\dot{\lambda}$.

### 3 Numerical implementation

In order to numerically solve the above stated elastoplastic problem in the context of the finite element method, it is more convenient to replace the equilibrium equation (3) by the virtual power principle

$$\forall \vec{U}^*, \quad \int_{\Omega_o} \vec{\Pi}^T : \overline{\overline{\mathrm{grad}} \vec{U}}^* d\Omega_o - \int_{\Omega_o} \vec{f}_o \vec{U}^* d\Omega_o - \int_{S_o} \vec{U}^* \cdot \vec{\Pi} \cdot \vec{N} dS_o = 0 \quad (12)$$

The nominal stress vector $\vec{T} = \vec{\Pi} \cdot \vec{N}$ related to the reference normal $\vec{N}$ is known on the part of boundary where dead loads are prescribed.

The discrete formulation is obtained using the following vector and matrix notations. A second-order tensor $\vec{T}$ (symmetric or not) is represented by a column vector $\{T\}$ of nine components as follows

$$\{T\}^T = \{T_1, T_2, T_3, T_4, T_5, T_6, T_7, T_8, T_9\}^T \\ = \{T_{11}, T_{21}, T_{31}, T_{12}, T_{22}, T_{32}, T_{13}, T_{23}, T_{33}\}^T \quad (13)$$

Thus, the components of vector $\{T\}$ are related to those of tensor $\vec{T}$ by $\{T\}_k = T_{pk} \leftrightarrow k=3(r-1)+p$. It should be noted that second order tensors are represented by vectors with nine components instead of six as usually done. In fact, this representation mode, which is necessary for unsymmetric matrices such as the displacement gradient $\vec{H}$ or the stress $\vec{\Pi}$, does not entail extra computational cost, as we shall see below. Similarly, in matrix notation, a fourth-order tensor $\vec{T}$ is represented by a $(9 \times 9)$ matrix $[T]$ as follows

$$[T] = \begin{bmatrix} T_{1111} & T_{1112} & T_{1113} & T_{1121} & T_{1122} & T_{1123} & T_{1131} & T_{1132} & T_{1133} \\ T_{2111} & T_{2112} & T_{2113} & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ T_{3111} & T_{3112} & T_{3113} & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ T_{1211} & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ T_{2211} & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ T_{3211} & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ T_{1311} & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ T_{2311} & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ T_{3311} & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \end{bmatrix} \quad (14)$$

Explicitly, the components of matrix $[T]$ are related to those of tensor $\vec{T}$ by $T_{kl} = T_{prqs} \leftrightarrow k=3(r-1)+p, l=3(s-1)+q$.
Within an element with nne nodes (nne = node number of the element), the displacement gradient $\bar{H} \equiv \overline{grad} \bar{U}$ is approximated by

$$ \{H\} = [G]\{U\}^e \quad (15) $$

In the above relation, vector {$H$} has 9 components, {$U$}^e is the vector of element nodal displacements with 3nne components, matrix [G] is of dimension (9×3nne) and defined by

$$ [G] \equiv \begin{bmatrix} N_{1,1} \mathbf{I}_3 & N_{2,1} \mathbf{I}_3 & & & N_{nne,1} \mathbf{I}_3 \\ N_{1,2} \mathbf{I}_3 & N_{2,2} \mathbf{I}_3 & & & N_{nne,2} \mathbf{I}_3 \\ N_{1,3} \mathbf{I}_3 & N_{2,3} \mathbf{I}_3 & & & N_{nne,3} \mathbf{I}_3 \end{bmatrix} \quad (16) $$

where $N_i$ are the usual shape functions and $\mathbf{I}_3$ denotes the unit matrix of the 3rd order.

Note that matrix [G] does *not* depend on the nodal displacements. Moreover, it has the advantage of being sparse and *regularly* populated. This allows us to compute the internal force vector as well as the tangent matrix in a simple and efficient way, as shown in proposition 2 below.

The discretization and consistent linearization of the elastic-plastic equation set lead to a nonlinear matrix system which can be solved by an incremental procedure. The convergence at every time step requires a Newton type iterative scheme. Within each iteration, the local integration of the elastoplastic equations set (8)-(10) is carried out in order to compute the stress $\overline{\Pi}$ and the consistent tangent tensor $\partial \overline{\Pi}/\partial \overline{F}$ corresponding to a given strain increment $\Delta \overline{E}$. Then we can build up the element internal force vector and the element tangent matrix defined as

$$ \{\Psi\}^e = \int_e [G]^T \{\Pi\} d\Omega_o \quad (17) $$

$$ [K]^e = \frac{\partial \{\Psi\}^e}{\partial \{U\}^e} = \int_e [G]^T \left[ \frac{\partial \Pi}{\partial \overline{F}} \right] [G] d\Omega_o \quad (18) $$

The local integration procedure is formally the same as in small strains, the rate equations are discretized using the fully implicit integration scheme (backward Euler scheme) and we get the correction for the plastic multiplier at the k-th local iteration [Simo and Taylor (1985)]

$$ \delta\lambda^{(k)} = -\Delta\lambda^{(k-1)} + \frac{\sqrt{\frac{3}{2} \|\overline{S}^E\|} - \sigma_o - A_{n-1}}{3\mu + h} \quad (19) $$

where $h=\partial^2 w^\alpha / \partial \alpha^2$, $\mu$ is the shear modulus, $\overline{S}^E$ the deviator of the elastic trial stress, $\overline{S}^E \equiv 2\mu (\mathrm{dev}(\overline{E}) - \overline{E}_{n-1}^p)$. The subscript (n-1) denotes known quantities at the previous time step. The plastic multiplier increment $\Delta\lambda$ obtained at convergence allows us to update the unsymmetric stress $\overline{\Pi}$ by using relation $\overline{\Pi} = \overline{F} \cdot \overline{\Sigma}$ and

$$ \overline{\Sigma} = \overline{S}^E - 2\mu \sqrt{\frac{3}{2} \Delta\lambda \frac{\overline{S}^E}{\|\overline{S}^E\|}} + \frac{1}{3} (3\lambda + 2\mu) (\mathrm{tr}\overline{E}) \overline{I} \quad (20) $$

On the other hand, the solution convergence at each iteration requires computing the consistent tangent modulus $\partial \overline{\Sigma} / \partial \overline{E}$

$$ \frac{\partial \overline{\Sigma}}{\partial \overline{E}} = \overline{D} - \frac{6\mu^2}{3\mu+h} \frac{\overline{S}^E \otimes \overline{S}^E}{\|\overline{S}^E\|^2} - 4\mu^2 \sqrt{\frac{3}{2}} \frac{\Delta\lambda}{\|\overline{S}^E\|} \left( \overline{I} - \frac{\overline{f}}{3} - \frac{\overline{S}^E \otimes \overline{S}^E}{\|\overline{S}^E\|^2} \right) \quad (21) $$

Eventually, the consistent tangent tensor $\partial \overline{\Pi}/\partial \overline{F}$ entering in the construction of global structural tangent matrix is derived from the above tensor $\partial \overline{\Sigma}/\partial \overline{E}$ by means of the following proposition.

**Proposition 1**

By using an appropriate change of variables, the following symmetry condition can always be satisfied for the expression of stress $\overline{\Sigma} : \forall \overline{T}$ (not necessarily symmetric),

$$ \overline{\Sigma}(\overline{T}, \text{state at } t_{n-1}) = \overline{\Sigma}(\overline{T}^T, \text{state at } t_{n-1}) \quad (22) $$

Then one has

$$ \frac{\partial \overline{\Pi}}{\partial \overline{F}} = (\overline{I} \cdot \overline{\Sigma})^T + \overline{F} \cdot \frac{\partial \overline{\Sigma}}{\partial \overline{E}} \overline{F}^T \quad (23) $$

It can be shown that the above discretization procedure leads to existence of a discrete potential $\tilde{w}(\overline{F}$, state at $t_{n-1}$) such that $\overline{\Pi}^T = \partial \tilde{w}(\overline{F}) / \partial \overline{F}$. We deduce that tensor $\partial \overline{\Pi}/\partial \overline{F}$ possesses themirror symmetry : $(\partial \overline{\Pi}/\partial \overline{F})_{ijkl} = (\partial \overline{\Pi}/\partial \overline{F})_{lkji}$ and that its matrix representation $[\partial \overline{\Pi}/\partial \overline{F}]$ is symmetric. From (18) there follows the symmetry of the tangent consistent modulus.

To end this section, let us note that one can take advantage of the sparse and regular nature of matrix [G], Relation (16), in order to prove the following proposition and
remarkably reduce the operation number for the computation of the element internal force vector {$\Psi$}^e and the element tangent stiffness matrix [K]^e, Relations (17) and (18).

**Proposition 2**

∀i, j∈{1, ..., 3nne},

$$ \Psi_i^e = \int_e N_{a,r} \Pi_{3(r-1)+p} N_{b,s} d\Omega_o \quad (24) $$

$$ K_{ij}^e = \int_e N_{a,r} \frac{\partial \Pi_{pr}}{\partial F_{qs}} N_{b,s} d\Omega_o \quad (25) $$

where implicit summations are made over r, s∈{1,2,3}, whereas integers a, b∈{1, ..., 3nne}, and p, q∈{1,2,3} are expressed in terms of integers i, j by $a=\text{Int}(\frac{\pm 2}{3})$, $b=\text{Int}(\frac{\pm 2}{3})$, $p=i-3(a-1)$, $q=j-3(b-1)$

(Int=integer part of)

Relations (24) and (25) show that computing the element internal force vector and the element tangent stiffness matrix is not so tremendous a task.

# 4 Numerical results

## 4.1 Geometry and boundary conditions

The numerical computations are performed on a circular cylindrical bar of initial length $2l_0=48$mm and of initial radius $r_0=4$mm.

Concerning the boundary conditions, the bar is modeled in two ways:

(i) The first one is most commonly adopted in the literature. During the loading process, the ends of the bar remain shear free so that the radial displacement thereat is not prevented. To trigger the initiation of necking, an axisymmetric geometric imperfection is introduced in shape of a full cosine wave with an amplitude of 1% of the initial radius.

(ii) Another way to compute the necking behavior is to consider a geometrically perfect bar with its ends cemented to rigid grips, as in [Needleman (1972)] and [Brünig (1998)]. Note that in this case there is no need of a geometrical defect for the onset of necking.

## 4.2 Material properties

The Young modulus is $E=2.10^{11}$Pa, the Poisson ratio $\nu=0.3$ and the initial yield stress $\sigma_0=400$MPa. The coefficients of the hardening law (8b) are so adjusted that

<figure><img src="image_2.png" /><figcaption>**Figure 1:** Perfect bar with gripped ends. Initial and final meshes.</figcaption></figure>

the numerically derived load-displacement curve best approximates the experimental curve of a given standard steel (XC48C): $c_0=0$ (no initial hardening), $c_1=220$MPa, $c_2=-560$MPa and $c_3=15$.

## 4.3 Meshes

Because of obvious symmetry, we only have to consider one quarter of the specimen the length of which is $l_0=24$mm. For each boundary condition type, the numerical computations are carried out using two axisymmetric meshes (Figure 1):

(i) A coarse mesh with 125 axisymmetric eight-node quadrilaterals.

(ii) In order to check the sensitivity of the numerical results to mesh refinement, another computation is performed over a finer mesh with 400 elements, figure 1b.

## 4.4 Numerical results in the case of the shear free end bar with imperfection

The bar is subjected to axial end displacements corresponding to a maximum relative elongation of 23% attained after 25 time steps and Figure 1 shows the undeformed and final configurations.

With the 125- and 400-element meshes, on average 4 iterations are necessary to satisfy the convergence tolerance of 1.e-6, with a peak of 8 iterations just after going through the limit point on the load-displacement curve
<figure><img src="image_3.png" /><figcaption>Figure 2: Load-displacement curves</figcaption></figure>

shown in Figure 2. As a matter of fact, the plots corresponding to the coarse and the finer meshes are found to be indistinguishable, thus validating the insensitivity of the numerical results to mesh refinement.

The load-displacement curve representing the tensile load F versus the axial end displacements U<sub>B</sub> is shown in Figure 2. The tensile load, after reaching its maximum value equal to 1.48 times the elastic yield force, decreases when the necking takes place. For a relative elongation of 23%, it falls down to 76% of its maximum value.

For comparison, we also perform the computations over a *perfectly cylindrical* bar. As expected, this bar remains cylindrical during the whole deformation process and the necking never takes place. The corresponding tensile force plotted in Figure 2 attains its maximum value later than the imperfect bar and very slightly decreases afterwards. We also find that bifurcation occurs about the maximum load point, as in [Needleman (1972)] and [Brünig (1998)], although the precise location of bifurcation point is not known. It is noteworthy that the maximum load point (1.48 times the elastic yield force) and the bifurcation point (engineering strain of 16% approximately) are very close to those obtained by [Okada, Rajiyah and Atluri (1990)] for the necking of a tensile plate in plane strain.

Now let us look at the change in the radius of the central cross-section of the bar which is an essential indicator of the necking behavior. At the beginning, the deformation

<figure><img src="image_4.png" /><figcaption>Figure 3: Necking radius vs. axial elongation</figcaption></figure>

is almost homogeneous and the bar remains cylindrical. However, beyond a certain elongation the radial displacement becomes suddenly more pronounced in the central region of the bar, and eventually the high concentration of strain there leads to necking phenomena.

The curve in Figure 3 giving the central cross-section radius versus elongation is tangent as expected to that of the geometrically perfect bar, which is virtually a straight line.

Beyond an engineering axial strain of 15% the former rapidly deviates from the latter, and at an elongation of 23% the current radius drops to about 36% of its initial value. This behaviour is very close to the numerical and experimental results depicted by [Simo (1988)] although the material considered and particularly the formulations adopted in the previous reference and the present work are quite different. The curve shown in Figure 3 can also be compared with that computed by [Brünig (1998)], yet noting that the last Reference gives the change in cross-section *area* instead of the change in cross-section radius.

Figure 4 shows the evolution of the plastic region at different deformation stages of the imperfect shear free end bar. The plastic deformation begins in the central region when the axial elongation U<sub>B</sub>/l<sub>o</sub> is equal to σ<sub>o</sub>/E=0.2% (Figure 4a), and it rapidly sweeps through the whole bar afterwards. Up to an elongation of 15% - i.e. until the maximum load level (see Figure 1) - the bar still remains entirely plastic (Figure 4b). Yet immediately after - at 16% elongation - there suddenly occurs a remarkable de-
<figure><img src="image_5.png" /><figcaption>**Figure 4:** Geometrically imperfect bar with shear free ends. Evolution of plastic regions (shaded area)</figcaption></figure>

<caption>**Figure 5:** Perfect bar with gripped ends. Initial and final meshes.</caption>

crease of the plastic zone (Figure 4c). For larger extensions of the bar, bifurcation has occurred and unloading proceeds rapidly. At the final elongation of 23% only the central necked region of the bar remains plastic whereas all other parts are unloaded elastically (Figure 4d).

The changes in the shape of the plastic zone observed above agree qualitatively well with those depicted in [Brünig (1998)]. However, in the last Reference the plastic zone decreases more rapidly, for instance at 20% elongation of the bar it is virtually concentrated in the central 1/10th of the bar length. Elastically unloaded regions are also detected in the necking of a tensile plate in plane strain, as shown in [Atluri and Murakawa (1980)], and their spread is quite similar to that depicted in Figure 4.

The maximum value of the Lagrangian equivalent plastic strain P (see Equation (11)) at the final elongation of 23% amounts approximately to 87% and is found on the lateral surface, in the immediate vicinity of the central cross-section of the bar. At the same time, the minimum value of P is about 13%, at a point on the circumference of the end sections.

## 4.5 Numerical results in the case of the gripped end perfect bar

So far the ends of the bar are shear free and the onset of necking is due to an initial small geometrical defect, as most commonly modeled in the literature. It should

be noted that [Needleman (1972)] and [Brnig (1998)] also considered another manner to simulate the necking, by taking a geometrically perfect bar with its ends cemented to rigid grips. For comparison, we take up again the numerical computations with these new geometry and boundary conditions.

With the coarse mesh (125 elements) the maximal engineering strain of 23% requires 25 time steps, whereas with the finer mesh (400 elements) it requires one more step, just to get across the limit point when the load reaches its maximum value. Apart from this slight difference in computational effort, all the curves resulting from the two meshes are quite identical.

Figure 5 shows the undeformed and final configurations. Here the gripped end boundary conditions act as an alternative geometrical imperfection and allows the necking bifurcation from the uniform tensile stress state. Despite the evident difference in the deformed shapes, the load-displacement curves shown in Figure 2 are almost the same for the gripped and free end bars.

In contrast, Figure 3 shows that the change in the cross-sectional radius in the necked region is rather different: in the case of the perfect bar with gripped ends, the change in the central cross-section radius is slower on the end part of the curve. At an elongation of 23%, the radius is equal to 44% of its initial value instead of 36% as in the case of shear free ends.
The necking radius curves in Figure 3 are different from those given in [Needleman (1972)]. However, they compare very well with those in [Brünig (1998)]. Note that in the last Reference another geometry - more realistic - was also studied, which is a whole cylindrical tensile specimen with two shear free thick ends. This geometry is not considered in the present work, but according to [Brünig (1998)] the results of the griped end perfect bar must be identical to those of the whole specimen.

Figure 6 shows that the evolution of the plastic zones for the gripped end perfect bar differs notably from that of the shear free imperfect bar, especially at the beginning. When the axial elongation $U_B/l_o$ is equal to $\sigma_o/E=0.2\%$, there exist separate plastic zones of different sizes in the gripped end bar. The small plastic zones, hardly visible in Figure 6a, are situated on the end cross-sections and near the lateral surface. The much larger plastic zones lie around the axis, between the central cross-section and the ends of the bar. As the overall extension increases, these plastic zones grow up together and eventually becomes connected. When the maximum load is reached - i.e. at about 15% elongation (see Figure 1) - the entire bar is plastically loaded, except for two small regions near the ends that will remain elastic throughout the whole range of straining (Figure 6b). Just after, with a 1% additional elongation only, unloading takes place and then regularly spreads out (Figure 6c). At the final elongation of 23%, the plastic zone is reduced to the central region of the bar, over a volume slightly smaller than that of the shear free end imperfect bar (Figure 6d).

The development of plastic zones shown above is very close to that described in [Needleman (1972)]. In particular, as in the quoted Reference, we also observe separate plastic zones which become connected, and the existence of small regions near the ends which remain continuously elastic. In the case of the gripped end perfect bar considered here, less details are given in [Brünig (1998)] for the beginning of the plastic zones evolution, however for the final stage all the results are again in quite good comparison.

The maximum value of the Lagrangian equivalent plastic strain $P$ at the final elongation of 23% amounts approximately to 83% and is found on the lateral surface, a little bit farther from the central cross-section of the bar than in the case of the shear free end imperfect bar. Of course, the equivalent plastic strain $P$ remains zero at the few points near the ends that remain elastic throughout.

<figure><img src="image_6.png" /><figcaption>**Figure 6:** Perfect bar with gripped ends. Evolution of plastic regions (shaded area)</figcaption></figure>

## 5 Conclusions

It is the purpose of this work to show the ability of the proposed formulation of finite elastoplasticity to deal with necking bifurcation problems. Below are summarized the main features of the theoretical and numerical formulations. Besides, further comments are also added on the results obtained in the preceding sections.

(i) The proposed theory is based on thermodynamic foundations established for generalized standard materials and results in a total Lagrangian formulation for finite plasticity. Since all relations are expressed in terms of Lagrangian variables, the material frame indifference is automatically satisfied within this formulation. The deformation history in the material has been characterized by a nonlinear relation of the saturation type between the isotropic hardening variable and its conjugate, Equation (8b).

(ii) The numerical implementation using standard finite element analysis is simple and yet robust enough to accurately reproduce the necking phenomenon. Written in terms of the Green strain tensor $\bar{E}$ and the symmetric Piola-Kirchhoff stress tensor $\bar{\Sigma}$, the local integration procedure is formally the same as in small strains. The rate equations are discretized using the fully implicit integration scheme (backward Euler scheme), which ensures both the stability and the symmetry of the consistent tangent modulus. In computing the global tangent matrix,
the update of matrix [∂Π/∂F] is made by means of Relations (21) and (23).

On the other hand, Relations (24) and (25) show that computing the element internal force vector and the element tangent stiffness matrix is not a highly time-consuming operation.

(iii) The numerical computations of necking are carried out both on imperfect bars with shear free ends and geometrically perfect bars with rigid gripped ends.

For each boundary condition type, the numerical computations carried out using a coarse and a finer meshes are found to yield indistinguishable results with the same rate of solution convergence, thus validating the insensitivity of the numerical results to mesh refinement.

The robustness of the implemented solution procedure is successfully assessed since 23% elongation can be attained in only 25 time steps, giving sensible necking simulation of the tensile bar.

All the obtained results have been quantitatively and qualitatively analyzed in some detail. Comparisons are also made with other formulations in the literature, showing agreements or differences in the load-displacement relationship, the necking radius versus axial elongation, and the evolution of plastic zones during the straining history.

(iv) It should be noticed that necking is usually referred to as a bifurcation phenomenon in the literature, whereas a true bifurcation computation of bar necking - as done in [Needleman (1972)] - is rare. [Needleman (1972)] considered a geometrically perfect bar with shear free ends and performed a real bifurcation analysis by discretizing a variational principle of Hill with the finite element method. As commonly done in other works, here we have solved the necking as a standard deformation problem by either introducing an initial imperfection of 1% in the radius or preventing the radial end displacement.

In connection with bifurcation aspects, let us note an interesting issue not often considered, which concerns the sensitivity to geometric or material parameters. In the present work, the numerical computations have been carried out with a unique set of bar dimensions and material values. The sensitivity study requires more extended discussions and is beyond the scope of this work.

(v) The computer program used here can be improved with path-following techniques such as arc-length and branch-switching methods, in order to deal with more

general problems in finite displacements and plastic strains, including limit point instabilities and bifurcation points. First results due to these enhancements have been obtained for the pinched cylinder modeled with shell elements in finite plasticity [Le Grognec and Le van (2000)]. Further investigations to treat static plastic buckling of shell-type structures are being developed.

## References

Argyris J. H.; Doltsinis J. St. (1979) : On the large strain inelastic analysis in natural formulation. Part I : quasistatic problems, *Comp. Meth. Appl. Mech. Eng.*, 20, pp. 213-251.

Atluri S. N. ; Murakawa H. (1980) : New general and complementary energy theorems, finite strain, rate sensitive inelasticity and finite elements : some computational studies, pp. 28-48, in *Nonlinear Finite Element Analysis in Structural Mechanics*. Proceedings of the Europe-US Workshop, ed. by W. Wunderlich et al, July 38-31, 1980, Ruhr-Universitat Bochum, Germany.

Besdo D. (1987) : Total Lagrangian strain-space-representation of the elasto-plasticity of metals, pp. 1357-1364, in Proceedings of the 2nd Int. Conf. on Constitutive Laws for Eng. Materials : Theory and Applications, vol. II, ed. by C.S. Desai et al, January 5-8, 1987, Tucson, Arizona, USA.

Brünig M. (1998) : Numerical analysis and modeling of large deformation and necking behavior of tensile specimens, *Finite Elements in Analysis and Design*, 28, pp. 303-319.

Casey J., Naghdi P. M. (1980) : A remark on the use of the decomposition $F=F^e F^p$ in plasticity, *J. Appl. Mech.*, 47, pp. 672-675.

Germain P.; Nguyen Q. S., Suquet P. (1983) : Continuum thermodynamics. J. Appl. Mech., 50, pp. 1010-1020.

Green A. E.; Naghdi P. M. (1965) : A general theory of an elastic-plastic continuum", *Archive for Rational Mechanics and Analysis*, 18, pp. 251-281.

Green A. E.; Naghdi P. M. (1971) : Some remarks on elastic-plastic deformation at finite strain. *Int. J. Engng. Sci.*, 9, pp. 1219-1229.
**Halphen B.; Nguyen Q. S. (1975)** : Sur les matériaux standard généralisés, *J. de Mécanique*, 14, 1, pp. 39-63.

**Kleiber M. (1992)** : On solving problems of advanced plastic flow, pp. 349-361, in *Finite Inelastic Deformations - Theory and Applications, IUTAM Symposium Hannover/Germany 1991*, D. Besdo, E. Stein (Eds.), Springer Verlag.

**Le Grognec P.; Le van A. (2000)** : Buckling and post-buckling of thin elastoplastic cylindrical shells with finite rotations. 2000 ASME International Mechanical Engineering Congress, November 5-10, 2000, Orlando, USA.

**Needleman A. (1972)** : A numerical study of necking in circular cylindrical bars, *J. Mech. Phys. Solids*, 20, pp. 111-127.

**Okada H.; Rajiyah H., Atluri S. N. (1990)** : A full tangent stiffness field-boundary-element formulation for geometric and material non-linear problems of solid mechanics, *Int. J. Num. Meth. Eng.*, 29, pp. 15-35.

**Rappaz M.; Bellet M.; Deville M. (1998)** : *Modlisation Numrique en Science et Gnie des Matriaux*. Volume 10 of the collection "Trait des Matriaux". Presses Polytechniques et Universitaires Romandes.

**Simo J. C. (1988)** : A framework for finite strain elastoplasticity based on maximum plastic dissipation and the multiplicative decomposition. Part I : Continuum formulation, *Comp. Meth. Appl. Mech. Eng.*, 66, pp. 199-219; Part II : Computational aspects, *Comp. Meth. Appl. Mech. Eng.*, 68, pp. 1-31.

**Simo J. C. (1992)** : Algorithms for static and dynamic multiplicative plasticity that preserve the classical return mapping schemes of the infinitesimal theory, *Comp. Meth. Appl. Mech. & Engng.*, 99, 61-112.

**Simo J. C.; Hughes T. H. R. (1998)** : *Computational Inelasticity*, Springer.

**Simo J. C.; Taylor R. L. (1985)** : Consistent tangent operator for rate-independent plasticity. *Comp. Meth. in Appl. Mech. & Engng*, 48, pp. 101-118.