# Modeling and computation of necking instability 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 computation of necking instability in round bars using a total lagrangian elastoplastic formulation. European Congress on Computational Methods in Applied Sciences and Engineering, Sep 2000, Barcelone, Spain. hal-01381102

HAL Id: hal-01381102

https://hal.science/hal-01381102

Submitted on 14 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 COMPUTATION OF NECKING INSTABILITY IN ROUND BARS
## USING A TOTAL LAGRANGIAN ELASTOPLASTIC FORMULATION

Anh Le van and Philippe Le Grognec

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.
e-mail: anh.le-van@ec-nantes.fr, philippe.le-grognec@ec-nantes.fr

**Key words:** necking instability, plasticity, total Lagrangian formulation.

**Abstract.** Necking is an instability phenomenon observed in round bars under tensile loading and has been investigated in numbers of works. In the present work, it is modeled within the framework of finite rate-independent plasticity. The theory is based on thermodynamical foundations for standard materials¹ and results in a total Lagrangian formulation for finite plasticity. The total strain is decomposed additively according to Green and Nagdhi². The yield function is of von Mises type and depends on the symmetric Piola-Kirchhoff stress tensor. Considering nonlinear isotropic hardening, the hardening law is written in the form of an exponential function. However, here it relates Lagrangian (material) variables, namely the Lagrangian hardening variable and the Lagrangian effective plastic strain. The local integration 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 implementation is simple as it works with classical finite elements. For the necking problem considered herein the proposed algorithm is robust enough to pass through the first critical point with no use of line-search.

The numerical computations are performed both on shear free end bars including an initial small geometrical defect and geometrically perfect bars with gripped ends. The derived load-displacement curve compares very well with those obtained in the literature, although the material considered and particularly the formulations adopted are quite different.
# 1 INTRODUCTION

Necking is an instability phenomenon observed in round bars under tensile loading and has been investigated in numbers of works for its modeling and numerical computation. The pionnering papers dated back to the 70's and were due to Needleman³, Mc Meeking and Rice⁴, Argyris and Doltsinis⁵. The obtained results were rather limited, in particular the computed deformed meshes did not resemble experimental observations. Besdo⁶ gave first a realistic necking configuration, he also seems to be the only one to solve the problem with a total Lagrangian formulation in the strain space for finite plasticity. More recent studies using other formulations have been developed by Simo⁷⁻⁹ and Brünig¹⁰, 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 critical point.

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

In the present work, necking is studied within the framework of finite rate-independent plasticity. The theory is based on thermodynamical foundations for standard materials¹ and results in a total Lagrangian formulation for finite plasticity. The total strain is decomposed additively according to Green and Nagdhi². The yield function is of von Mises type and depends on the symmetric Piola-Kirchhoff stress tensor. Considering nonlinear isotropic hardening, the hardening law is written in the form of an exponential function as in References 7 to 9. However, here it relates Lagrangian (material) variables, namely the Lagrangian hardening variable and the Lagrangian effective plastic strain. The local integration 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.

# 2 THEORETICAL BACKGROUND

Below we summarize the total Lagrangian formulation which is derived from thermodynamical foundations for standard materials¹. The initial configuration is chosen as the reference one and the equations of the problem will be expressed using Lagrangian variables relative to this configuration.

The Green strain tensor defined in terms of the displacement $\vec{U}$

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

is decomposed additively as
$$ \overline{\overline{E}} = \overline{\overline{E}}^e + \overline{\overline{E}}^p \quad (2) $$

According to Naghdi et al²,¹¹, the so-called plastic strain $\overline{\overline{E}}$ in Equation (2) is assumed to be symmetric and objective, with the result that the elastic strain $\overline{\overline{E}}$ is symmetric too.

The equilibrium equation reads

$$ \overline{\mathrm{div}} \overline{\overline{\Pi}} + \rho_o \overline{\overline{f}}_o = \overline{\overline{0}} \quad (3) $$

where $\overline{\overline{\Pi}}$ is the first Piola-Kirchhoff stress tensor. The stress $\overline{\overline{\Pi}}$ is not symmetric contrary to the second Piola-Kirchhoff stress tensor $\overline{\overline{\Sigma}}$ defined by $\overline{\overline{\Pi}} = \overline{\overline{F}} \cdot \overline{\overline{\Sigma}}$, where $\overline{\overline{F}}$ is the deformation gradient.

The free energy per unit volume $w$ is assumed to be the sum of the elastic energy $w^e$ - function of the elastic strain $\overline{\overline{E}}$ - and the hardening energy $w^\alpha$ - function of one scalar hardening variable denoted by $\alpha$

$$ w(\overline{\overline{E}}, \alpha) = w^e(\overline{\overline{E}}) + w^\alpha(\alpha) \quad (4) $$

Here we choose the elastic energy as a quadratic function of the elastic strain $\overline{\overline{E}}$, and the hardening energy containing exponential terms

$$ w^e(\overline{\overline{E}}) = \frac{1}{2} \overline{\overline{E}} : \tilde{D} : \overline{\overline{E}}^e \quad (5) $$

$$ w^{\alpha}(\alpha) = c_{0}\alpha + \frac{1}{2}c_{2}\alpha^{2} + (c_{1} - c_{0})(\alpha + \frac{\exp(-c_{3}\alpha)}{c_{3}}) \quad (6) $$

where $\tilde{D} = 2\mu\tilde{I} + \lambda\tilde{J}$ ($\tilde{I}$ is the fourth-order identity tensor, $\tilde{I} = \vec{e}_i \otimes \vec{e}_j \otimes \vec{e}_j \otimes \vec{e}_i$, $\tilde{J} = \overline{\overline{I}} \otimes \overline{\overline{I}}$, $\overline{\overline{I}}$ the second order identity tensor) denotes the isotropic elasticity tensor and the coefficients $c_i$ are material constants.

The state laws then yield the stress-strain relation and the nonlinear hardening law relating the hardening variable $A$ to its conjugate $\alpha$

$$ \overline{\overline{\Sigma}} = \frac{\partial w}{\partial \overline{\overline{E}}} = \tilde{D} : \overline{\overline{E}}^e = \tilde{D} : (\overline{\overline{E}} - \overline{\overline{E}}^p) \quad (7) $$

$$ A = \frac{\partial w}{\partial \alpha} = c_0 + c_2 \alpha + (c_1 - c_0)(1 - \exp(-c_3 \alpha)) \quad (8) $$

The yield function $f$ depending on the second Piola-Kirchhoff $\overline{\overline{\Sigma}}$ and the hardening variable $A$ is of von Mises type
$$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 provide the plastic strain and hardening variable rates

$$\dot{\bar{E}}^p = \dot{\lambda} \frac{\partial f}{\partial \bar{\Sigma}} \qquad -\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$.

Since the plastic strain $\dot{\bar{E}}^p$ is a Lagrangian variable, its rate is obtained by ordinary differentiation with respect to 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), (7)-(10) form a system of 30 scalar equations with 30 scalar unknowns $\vec{U}, \bar{E}, \bar{E}^e, \bar{E}^p, \bar{\Sigma}, A, \alpha$, and $\dot{\lambda}$.

### 3 NUMERICAL IMPLEMENTATION

The above stated elastoplastic problem is now numerically solved using the weak formulation, where the equilibrium equation (3) is replaced by the virtual power principle

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

where $\Omega_o$ and $S_o$ denote the reference initial domain and boundary of the body, respectively.

The nonlinear matrix system resulting from the discretization of the virtual power principle is solved by a time-stepping incremental procedure. The solution is computed at every time step by using the standard finite element method combined with a Newton type iterative scheme. Within each iteration, the local integration of the elastoplastic equations set (7)-(10) is carried out in order to compute the stress $\overline{\overline{\Pi}}$ (or $\overline{\Sigma}$) for a given strain increment $\Delta E$ and to build up the internal force vector. The integration procedure is 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 each local iteration

$$\Delta\Delta\lambda = -\Delta\lambda + \frac{\sqrt{\frac{3}{2}} \frac{\bar{S}}{\bar{\Sigma}} : \frac{\bar{S}}{\bar{\Sigma}} - \sigma_0 - A_{n-1}}{\|\bar{S}\|} \quad (13)$$
where h stands for the material property c₂ in (6), μ is the shear modulus, S the deviator of the elastic trial stress $\Sigma \equiv \Sigma_{n-1} + \tilde{D} : \Delta E$. The subscript (n-1) denotes known quantities at the previous time step. Then, the plastic multiplier increment Δλ obtained at convergence of the local integration allows us to update state variables in the classical way for the herein considered nonlinear isotropic hardening model.

We also compute the consistent tangent modulus $\partial \Sigma / \partial \bar{E}$ at each iteration by

$$ \frac{\partial \Sigma}{\partial \bar{E}} = \tilde{D} - \frac{6\mu^2}{3\mu + h} \frac{\Sigma \otimes \Sigma}{\|\Sigma\|^2} - 4\mu^2 \sqrt{\frac{3}{2} \frac{\Delta\lambda}{\|\Sigma\|^2} \left( \tilde{I} - \frac{\tilde{J}}{3} - \frac{\Sigma \otimes \Sigma}{\|\Sigma\|^2} \right)} \quad (14) $$

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

**Proposition.** By using an appropriate change of variables, the following symmetry condition can always be satisfied for expression of $\Sigma$:

$$ \forall \bar{T} \text{ (not necessarily symmetric), } \Sigma(\bar{T}, \text{state at } t_{n-1}) = \Sigma(\bar{T}, \text{state at } t_{n-1}) \quad (15) $$

then one has

$$ \frac{\partial \bar{\Pi}}{\partial \bar{F}} = (\tilde{I} \cdot \Sigma)^T + \bar{F} \cdot \frac{\partial \bar{\Sigma}}{\partial \bar{E}} \cdot \bar{F}^T \quad (16a) $$

In Cartesian co-ordinates, relation (16a) reads:

$$ \left( \frac{\partial \bar{\Pi}}{\partial \bar{F}} \right)_{ijkl} = \delta_{il} \Sigma_{kj} + F_{im} \left( \frac{\partial \bar{\Sigma}}{\partial \bar{E}} \right)_{mkn} F_{ln} \quad (16b) $$

# 4 NUMERICAL RESULTS

The numerical computations are performed on a circular cylindrical bar of initial length 2l₀=48mm and of initial radius r₀=4mm. The Young modulus is E=2.10¹¹Pa, the Poisson ratio v=0.3 and the initial yield stress σ₀=400MPa. The coefficients of the hardening law (8) are so adjusted that the numerically derived load-displacement curve best approximates the experimental curve of a given standard steel (XC48C) : c₀=0 (no initial hardening), c₁=220MPa, c₂=-560MPa and c₃=15.

The ends of the bar remain shear free so that the radial displacement there is not prevented, and the onset of necking is due to an initial small geometrical defect. Here we consider an
axisymmetric imperfection in shape of a full cosine wave with an amplitude of 1% of the initial radius.

Because of obvious symmetry, we only have to consider one quarter of the specimen the length of which is $l_0=24\text{mm}$. The numerical computations are carried out using two axisymmetric meshes:

(i) A coarse mesh with 125 eight-node isoparametric quadrilaterals, figure 1a.

(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.

The results corresponding to the two meshes are found to be indistinguishable.

The bar is subjected to axial end displacements corresponding to a maximum relative elongation of 23% attained after 25 time steps. 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 shown in Figure 2.

<figure><img src="image_2.png" /><figcaption>a. Mesh no. 1.</figcaption></figure>

<caption>b. Mesh no. 2.</caption>

<caption>Figure 1. Undeformed and deformed meshes.</caption>
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.

<figure><img src="image_3.png" /><figcaption>Figure 2. Load-displacement curve.</figcaption></figure>

In order to emphasize the role of the geometrical imperfection we take up again the computations over a *perfectly cylindrical* bar. Obviously, the bar remains now cylindrical along 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.

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. At the beginning, the deformation 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 10%, the former rapidly deviates from the latter : 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⁷, although the material considered and particularly the formulations adopted in the quoted reference and the present work are quite different. However, the curve shown in Figure 3 is rather diferent from that computed by Brünig¹⁰ (note that the change in cross-section *area* is reported in Reference 10 instead of the change in cross-section radius).

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

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³ and Brünig¹⁰ also considered another manner to simulate the necking, by taking a *geometrically perfect* bar with its ends cemented to rigid grips. For comparison, we resume the numerical computations with the new geometry and boundary conditions by applying the
presented formulation herein. The corresponding results show that the deformed configurations and the load-displacement curves are analogous to those of the imperfect bar. In contrast, 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.

# 5 CONCLUSIONS

It is the purpose of this work to show the abilities of the proposed formulation to deal with necking instability problems. Thoroughly expressed in terms of Lagrangian variables, it automatically satisfies the material frame indifference and requires no objective rates for the constitutive equations. The numerical implementation is simple as it works with classical finite elements. The same computer program can be used to solve different types of problems in finite displacements or finite plastic strains, with or without instability. In particular, it allows to study buckling and advanced post-buckling of structures of beam, plate or shell type. For these problems as for the necking one the proposed algorithm is robust enough to pass through the first critical point with no use of line-search.

# REFERENCES

[1] Halphen B., Nguyen Q. S., "Sur les matériaux standard généralisés", *J. de Mécanique*, **14**, 1, 39-63, (1975).

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

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

[4] Mc Meeking R. M., Rice J. R., "Finite-element formulations for problems of large elastic-plastic deformation", *Int. J. Solids & Struct.*, **11**, 601-616, (1975).

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

[6] Besdo D., "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, (1987).

[7] Simo J. C., "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**, 199-219, (1988); "Part II : Computational aspects", *Comp. Meth. Appl. Mech. Eng.*, **68**, 1-31, (1988).
[8] Simo J. C., "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, (1992).

[9] Simo J. C., Hughes T. H. R., *Computational Inelasticity*, Springer, (1998).

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

[11] Casey J., Naghdi P. M., "A remark on the use of the decomposition $F=F^cF^p$ in plasticity", *J. Appl. Mech.*, **47**, 672-675, (1980).