# An improved one-point integration method for large strain elastoplastic analysis

Laurent Stainier, Jean-Philippe Ponthot

► To cite this version:

Laurent Stainier, Jean-Philippe Ponthot. An improved one-point integration method for large strain elastoplastic analysis. Computer Methods in Applied Mechanics and Engineering, 1994, 118 (1-2), pp.163-177. 10.1016/0045-7825(94)90111-2. hal-01005235

HAL Id: hal-01005235

https://hal.science/hal-01005235

Submitted on 2 Nov 2021

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
# An improved one-point integration method for large strain elastoplastic analysis

L. Stainier, J.Ph. Ponthot

L.T.A.S.-Thermomécanique, University of Liège, Rue Ernest Solvay, 21, B-4000 Liège, Belgium

## Abstract

In this paper, we present a new one-point integration method generalizing Flanagan and Belytschko's method (Internat. J. Numer. Methods Engrg. 17 (1981) 679-706) and a modification of Belytschko and Bindeman's method (Comput. Methods Appl. Mech. Engrg. 88 (1991) 311-340), both in the frame of large deformation elastoplastic analysis. These stabilization methods are combined with the radial return method used to integrate the constitutive law. Plane strain problems are first considered, and the method is then generalized to axisymmetrical situations. The explicit time integration scheme with its critical timestep is also considered. A few examples are presented that show the great time savings that can be obtained with reduced integration without any loss of accuracy, and even with a gain in the solution quality, since the underintegrated elements prove to be 'flexurally superconvergent'.

## 1. Introduction

The use of one-point integration rules in modelling large inelastic strains via finite elements provides two advantages: (i) it prevents locking from appearing in nearly incompressible media (NIM), (ii) it reduces the number of operations in the integration of the constitutive laws for each element. Unfortunately, drawbacks are generated: spurious zero energy modes develop that can ruin the quality of the numerical solution. Hence, it is mandatory to add stabilization forces to the nodal internal forces to preserve the advantages of the one-point integration system.

Pioneer work in one-point integration elements and their stabilization was carried out by Flanagan and Belytschko [1, 3, 4]. They presented two-dimensional and three-dimensional small strains elements with stabilization by artificial damping or by artificial stiffness, the latter being preferred. Unfortunately, these stabilization methods necessitate a parameter to control their amplitude. Liu and Belytschko [5] extended the method to a heat conduction element. In the following years, the efforts focused on determining this control parameter and even getting rid of it. For that purpose, Belytschko et al. [6] introduced the idea of using a variational principle while Liu et al. [7, 8] proposed using a Taylor development of the strain rate.

Later, Belytschko and Bachrach [9] presented a quintessential bending element (QBE), with no locking for NIM. This element is based on a first order development of the strain field and on the use of the Hu-Washizu variational principle for determining the stress field. They propose several formulations for the linear shear field and compare them. It is also worth noting the work of Koh and Kikushi [10], who review different full or selective reduced integration methods. Later, Liu et al. [11] developed a QBE using the Hu-Washizu variational principle for non-linear problems. Jetteur and Cescotto [12] also used this variational principle to develop an element for the analysis of large inelastic strains without locking for NIM. Belytschko and Bindeman [2] used the assumed strain method of Simo and Hughes [13] and built a strain field avoiding the locking phenomenon. They also extended their element to non-linear problems.
In the present paper, after a brief review of the hourglass phenomenon (Section 2), we first extend the artificial stiffness method in [1] for small elastic strains to large elastoplastic deformation (Section 3). In Section 4, starting with the strain field established in [2], we also propose a straightforward use of the constitutive law (instead of a variational principle) and a modification of the stabilization forces in order to insure frame-independency. All these developments (Sections 3 and 4) are made in the context of the radial return method of Wilkins [14] for the integration of the elastoplastic constitutive law. In Section 5, we generalize these methods for the previously considered 2-D plane strain situation to the axisymmetrical one. Some attention is paid in Section 6 to the time dependent problems where an explicit scheme with a critical timestep is used. Finally, we present some examples illustrating the excellent behaviour of our one-point integration elements in Section 7.

## 2. The existence of hourglass modes in one-point integration for quadrilateral elements [1,3]

The strain rate $\dot{\epsilon}_{ij}$ can be written as the symmetric part of the material velocity gradient,

$$ \dot{\epsilon}_{ij} = \frac{1}{2}(\dot{u}_{i,j} + \dot{u}_{j,i}), \qquad (2.1) $$

where $\dot{u}_i$ is the velocity component in spatial direction $i$, and where the notation, $i$ denotes a derivative with respect to spatial coordinates $x_i$ and a superposed dot the material rate of change. In the context of a Lagrangian formulation, we can write

$$ \dot{u} = N_i \dot{u}_i \qquad (2.2) $$

since the material derivative of shape functions $N_i$ equals zero. Note that the Einstein convention on summation indices is assumed and that Latin capitals (e.g. $I$) are relative to each of the four vertices. The shape functions of an isoparametric bi-linear element are given by

$$ N_i = \frac{1}{2}(1 + \xi_i \xi)(1 + \eta_i \eta), \qquad (2.3) $$

where $\xi$ and $\eta$ are the reference square coordinates (no summation on the indices in this case).

If we write the strain rate tensor in a vectorial form, we have

$$ \begin{bmatrix} \dot{\epsilon}_{11} \\ \dot{\epsilon}_{22} \\ 2\dot{\epsilon}_{12} \end{bmatrix} = \begin{bmatrix} b_1^1 & 0 \\ 0 & b_2^1 \\ b_2^1 & b_1^1 \end{bmatrix} \begin{bmatrix} \dot{d}_1 \\ \dot{d}_2 \end{bmatrix} = B \cdot \dot{d}, \qquad (2.4) $$

where the shape functions derivatives are denoted by $b_{il} = N_{l,i}$. Working with one integration point only, all values will be calculated at the centroid of the element, i.e. at $\xi = \eta = 0$. The central shape functions gradients can then be written as

$$ \begin{aligned} b_1^1 &= \frac{1}{2A} [y_{24} \quad y_{31} \quad y_{42} \quad y_{13}], \\ b_2^1 &= \frac{1}{2A} [x_{42} \quad x_{13} \quad x_{24} \quad x_{31}], \end{aligned} \qquad (2.5) $$

where $A$ is the element area

$$ A = \frac{1}{2}(x_{31}y_{42} + x_{24}y_{31}), \qquad (2.6) $$

and where $x_{ij}$ and $y_{ij}$ are written for $x_i - x_j$ and $y_i - y_j$, respectively.

The null-space of the $B$ matrix $NS(B)$ is the set of all non-zero displacement vectors $d$ that do not produce any strain, i.e. $B \cdot d = 0$. Since $B$ is of size $3 \times 8$ and of rank 3 and there are 8 degrees of freedom, this null-space is of size 5 and is given by

$$ NS(B) = \begin{bmatrix} t & 0 & -x_2 & h & 0 \\ 0 & t & x_1 & 0 & h \end{bmatrix}, \qquad (2.7) $$

where
<figure><img src="image_2.png" /><figcaption>Fig. 1. The null-space of B.</figcaption></figure>

$$t^i = [1 \ 1 \ 1 \ 1], \quad h^i = [1 \ -1 \ 1 \ -1], \qquad (2.8)$$

and where the nodal coordinates vectors are respectively

$$x_1^i = [x_1 \ x_2 \ x_3 \ x_4], \quad x_2^i = [y_1 \ y_2 \ y_3 \ y_4]. \qquad (2.9)$$

The first three columns of NS(B) ($n=1,2,3$) are rigid body displacements (see Fig. 1), but the two additional ones ($n=4,5$) are spurious modes. They are called hourglass modes (because of their trapezoidal shape) and imply a non-zero strain except at the centroid, precisely where it is measured. If no care is taken, these modes become unbounded and can develop until the solution loses any physical significance. Thus, our goal will be to determine stabilization forces that prevent these modes developing without hampering the effective displacement modes.

Note that these vectors have the following fundamental properties:

$$b_i^t \cdot t = 0, \quad b_i^t \cdot h = 0, \quad t^i \cdot h = 0, \quad b_i^t \cdot x_j = \delta_{ij}. \qquad (2.10)$$

The first two of these relations signify that neither translation nor hourglass modes produce any strain at the centroid of the element while the third one shows the orthogonality between these two displacement modes. The fourth one is readily demonstrated, recalling that the $b_i$ are the shape functions derivatives.

### 3. Generalization of Flanagan and Belytschko's method to large elastoplastic strains

#### 3.1. Flanagan and Belytschko's artificial stiffness method [1, 3]

To control the hourglass modes, these authors propose introducing stabilization forces proportional to the amplitude of these modes for the element considered by the mean of an artificial stiffness. This leads to the following expression for the stabilization forces:
$$f_i^{HG} = \frac{\kappa}{2} A(b_j^t \cdot b_j) Q_i \gamma, \quad Q_i = (K + \frac{4}{3}\mu)q_i, \qquad (3.1)$$

where $\kappa$ is a user defined control parameter, $K$ is the material bulk modulus, $\mu$ its shear modulus, $A$ and $b_i$ are respectively defined in (2.6) and (2.5). The $Q_i$ can be considered as hourglass stresses, and the size $q_i$ of the hourglass mode in direction $i$ is given by

$$u_{ii}^{HG} = q_i h_i, \quad q_i = u_{ii} \gamma_i, \qquad (3.2)$$

with the hourglass shape vector

$$\gamma_i = \frac{1}{4}[h_i - (h_j x_{ij}) b_{ij}]. \qquad (3.3)$$

This vector corresponds to the hourglass mode for an element of arbitrary shape, while $h$ is relative to the reference square element.

The main disadvantage of this method is the presence of the control parameter $\kappa$ in (3.1). That is the reason why many authors [2,6-12] have tried to develop stabilization methods without any parameter, most of these methods being based on a linearized expression of the displacement field.

## 3.2. Extension of the method to large elastoplastic strains

In a non-linear situation, the hourglass stresses in (3.1) are incrementally computed, but there is no problem of objectivity, since the expression (3.1) is frame independent.

In the case of deviatoric plasticity (Von Mises yield surface for example), a simple and efficient method for the integration of the constitutive law over a load increment is the radial return of Wilkins [14]. This method consists of finding an approximate plastic stress along the radius issued from an elastic predictor (see Fig. 2 for the case of the Von Mises yield criterion with isotropic hardening). In practice, we first compute the equivalent elastic stress at the centroid of the element

$$\sigma_{v,i}^{e} = \sqrt{\frac{3}{2} s_{ij}^{e} s_{ij}^{e}} \qquad (3.4)$$

where $s_{ij}^e$ is the elastic deviatoric stress. If this equivalent elastic stress exceeds the yield stress $\sigma_v^n$ (the superscript $n$ recalls that it is the value obtained in the reference state, i.e. at the end of timestep $n$), then we need to calculate a plastic corrector. In the radial return method, the plastic deviatoric stress is supposed to lie along the radius issued from the elastic predictor

$$s_{ij} = \beta s_{ij}^{e} \quad (\beta \in ]0,1]), \qquad (3.5)$$

<figure><img src="image_3.png" /><figcaption>Fig. 2. The radial return method.</figcaption></figure>
where $\beta$ is computed such that the plastic stress be on the updated yield surface. For a linear isotropic hardening, we have

$$ \beta = 1 - \frac{\sigma_{\text{eq}}^{e} - \sqrt{\frac{2}{3}} \sigma_{\text{V}}^{n}}{\sigma_{\text{eq}}^{e} \left(1 + \frac{H}{3\mu}\right)}, \qquad (3.6) $$

where $H$ is the hardening coefficient.

This radial return method is nearly equivalent to using an 'effective shear modulus' $\beta\mu$ throughout the current step. This suggests to compute the hourglass stresses as follows:

$$ Q_i^{n+1} = Q_i^n + (K + \frac{4}{3}\beta\mu)q_i, \qquad (3.7) $$

the hourglass strains $q_i$ being computed from the displacements relative to the reference configuration (the previous equilibrium one since we work in updated Lagrangian formalism, ULF):

$$ q_i = \gamma_i^t \cdot d_i. \qquad (3.8) $$

# 4. An improvement of Belytschko and Bindeman's stabilization forces [2]

## 4.1. Small elastic strains case

In [2], it is shown that the following expression for the strain field leads to a stabilization without shear or volumetric locking:

$$ \epsilon = \epsilon^0 + \epsilon^{\text{HG}} = \begin{bmatrix} \epsilon_{11}^0 + eq_1 h_{,1} - eq_2 h_{,2} \\ \epsilon_{22}^0 - eq_1 h_{,1} + eq_2 h_{,2} \\ 2\epsilon_{12}^0 \end{bmatrix}, \qquad (4.1) $$

where $\epsilon_{ij}^0$ are the strains computed at the center, the $q_i$ are again the hourglass modes amplitudes and $h = \xi\eta$. The best value for the $e$ parameter is shown to be one half in [2]; this will be confirmed in the examples in Section 7.

Let us associate with (4.1) a stress field through a plane strain elastic constitutive law represented by the C matrix (Hooke matrix),

$$ \sigma = C \cdot \epsilon, \qquad (4.2) $$

with

$$ C = \begin{bmatrix} \lambda + 2\mu & \lambda & 0 \\ \lambda & \lambda + 2\mu & 0 \\ 0 & 0 & \mu \end{bmatrix}, \qquad (4.3) $$

where $\lambda$ and $\mu$ are the Lamé's constants.

The stress field is then written

$$ \begin{bmatrix} \sigma_{11} \\ \sigma_{22} \\ \tau_{12} \end{bmatrix} = \begin{bmatrix} \sigma_{11}^0 + Q_1 h_{,1} - Q_2 h_{,2} \\ \sigma_{22}^0 - Q_1 h_{,1} + Q_2 h_{,2} \\ \tau_{12}^0 \end{bmatrix}, \qquad (4.4) $$

where the $Q_i = 2e\mu q_i$ are the hourglass stresses associated with the hourglass strains $q_i$.

Now, let us define the $\tilde{B}$ matrix as the matrix relating the strains and the nodal displacements:
$$ \bar{B} = \begin{bmatrix} b_1^t + eh_{1,1}\gamma^t & -eh_{1,2}\gamma^t \\ -eh_{1,1}\gamma^t & b_2^t + eh_{2,1}\gamma^t \\ b_2^t & b_1^t \end{bmatrix}, \qquad (4.5) $$

Following the $\bar{B}$ method of Simo and Hughes [13], the nodal forces are then obtained by

$$ f = \int_V \bar{B}^t \cdot \sigma \, dV = f^0 + f^{HG}, \qquad (4.6) $$

where

$$ \begin{aligned} f^0 &= \int_V \bar{B}^{0t} \cdot \sigma^0 \, dV \\ f^{HG} &= \int_V \bar{B}^{HGt} \cdot \sigma^{HG} \, dV. \end{aligned} \qquad (4.7) $$

The uncoupling between the central and hourglass parts of the forces is due to the following property of $h$, which can be easily checked:

$$ \int_V h_{,i} \, dV = 0. \qquad (4.8) $$

If we introduce in (4.7) the expressions (4.4) and (4.5), we obtain for the nodal forces

$$ f^0 = A \begin{bmatrix} b_1 \sigma_{11}^0 + b_2 \tau_{12}^0 \\ b_2 \sigma_{22}^0 + b_1 \tau_{12}^0 \end{bmatrix} \qquad (4.9) $$

and

$$ f^{HG} = \begin{bmatrix} 2e(Q_1 H_{11} - Q_2 H_{12})\gamma \\ 2e(-Q_1 H_{12} + Q_2 H_{22})\gamma \end{bmatrix}, \qquad (4.10) $$

where

$$ \begin{aligned} H_{ij} &= \int_V h_{,i} h_{,j} \, dV \\ H_{11} &= \frac{2}{3} \frac{y_{31}^2 + y_{42}^2}{A}, \quad H_{22} = \frac{2}{3} \frac{x_{31}^2 + x_{42}^2}{A}, \quad H_{12} = -\frac{2}{3} \frac{x_{31}y_{31} + x_{42}y_{42}}{A}. \end{aligned} \qquad (4.11) $$

Note that we have assumed the Jacobian matrix determinant to be constant over the element when computing the integrals. An interesting feature of these stabilization forces is that they do not depend on the bulk modulus $K$, thus allowing them to be used in the analysis of incompressible materials (e.g. rubber).

## 4.2. Extension to the case of large strains and rotations

Let us now check if these forces are objective, i.e. independent of the orientation of the reference frame, by assuming a rigid rotation of an angle $\theta$ described by

$$ \begin{bmatrix} \tilde{x} \\ \tilde{y} \end{bmatrix} = \begin{bmatrix} \cos \theta & -\sin \theta \\ \sin \theta & \cos \theta \end{bmatrix} \begin{bmatrix} x \\ y \end{bmatrix}, \quad \tilde{x} = R \cdot x. \qquad (4.12) $$

Replacing the old nodal positions by the new ones given by (4.12), we can see that the forces (4.10) are frame dependent. Thus, we have to modify their expression. Starting from the general expression:
$$ \begin{bmatrix} f_1^{\text{HG}} \\ f_2^{\text{HG}} \end{bmatrix} = \begin{bmatrix} 2e((AH_{11} + BH_{12} + CH_{22})Q_1 + (DH_{11} + EH_{12} + FH_{22})Q_2)\gamma \\ 2e((DH_{22} + EH_{12} + FH_{11})Q_1 + (AH_{22} + BH_{12} + CH_{11})Q_2)\gamma \end{bmatrix} \quad (4.13) $$

and stating that $\tilde{f}^{\text{HG}} = R \cdot f^{\text{HG}}$, we find that the only possible frame invariant expression of the hourglass nodal forces is

$$ \begin{bmatrix} f_1^{\text{HG}} \\ f_2^{\text{HG}} \end{bmatrix} = \begin{bmatrix} 2e(H_{11} + H_{22})Q_1\gamma \\ 2e(H_{11} + H_{22})Q_2\gamma \end{bmatrix}. \quad (4.14) $$

By comparison to (4.10), it means that we take $H_{12} = 0$ and replace both $H_{11}$ and $H_{22}$ by a global value $H_{11} + H_{22}$.

Thus, we shall use the following nodal forces, where we have replaced $H_{11} + H_{22}$ by its expression in function of the $b_i$:

$$ f^0 = A \begin{bmatrix} b_1 \sigma_{11}^0 + b_2 \tau_{12}^0 \\ b_2 \sigma_{22}^0 + b_1 \tau_{12}^0 \end{bmatrix}, \quad (4.9) $$

$$ f^{\text{HG}} = \frac{8}{3} e A (b_1^t \cdot b_1 + b_2^t \cdot b_2) \begin{bmatrix} Q_1 \gamma \\ Q_2 \gamma \end{bmatrix}. \quad (4.15) $$

As already stated, the formalism we use is the ULF, in which the reference state for each step is the latest equilibrated configuration available. In order to guarantee incremental objectivity, we work with corotational strains and stresses and with the instantaneous final rotation. This consists of splitting the total strain in a constant strain without rotation giving the corotational stresses and a pure rotation occurring instantaneously at the end of the time increment, leading to the Cauchy stresses (for more details, see [15]).

If we used this method with the hourglass quantities, an hourglass shear stress would appear because of the final rigid rotation. As we have designed our hourglass strain field without shear to avoid locking, we cannot accept a shear stemming from the rotation. That is why we use the above expression for the hourglass forces in which we introduce the nodal displacements during the considered step. These forces being frame invariant, objectivity is preserved. Thus, we shall work with the following incremental scheme:

(i) the hourglass strains are calculated from the nodal displacements relative to the reference configuration

$$ q_i = \gamma^t \cdot d_i; \quad (4.16) $$

(ii) the hourglass stresses are obtained from the hourglass strains; the elastic predictor is given by

$$ \Delta Q_i = 2e\mu q_i, \quad Q_i^{n+1,e} = Q_i^n + \Delta Q_i, \quad (4.17) $$

and the corrected plastic stresses are (the radial return method is still used)

$$ Q_i^{n+1} = \beta Q_i^{n+1,e}, \quad (4.18) $$

where $\beta$ is given by (3.6);

(iii) the stabilization forces are calculated from the updated stresses,

$$ f_i^{\text{HG}} = \frac{8e}{3} A (b_j^t \cdot b_j) Q_i \gamma. \quad (4.19) $$

## 5. Axisymmetrical geometry

In the case of an axisymmetrical geometry, the analysis is limited to a one radian sector with radius-varying thickness. The geometrical properties (mass, stiffness, etc.) of a surface element are thus dependent on its distance with respect to the revolution axis. But if we work with one-point integration
elements, the properties of such an element are considered to be those at its center in the radial plane, thus leading to an error in the repartition of the nodal forces, and this error increases when the element approaches the central axis. That is why we have to correct these forces to take into account the variation of the radius across the element.

In the axisymmetrical case, the nodal forces are written (the stresses being measured at the centroid, they are constant)

$$f = \left( \int_{\Omega} \bar{B}^t r \, d\Omega \right) \cdot \sigma, \qquad (5.1)$$

where

$$\bar{B}^t = \begin{bmatrix} N_{,1} & 0 & \frac{1}{r}N & N_{,2} \\ 0 & N_{,2} & 0 & N_{,1} \end{bmatrix}, \quad \sigma^t = [\sigma_r \ \sigma_z \ \sigma_{\theta} \ \tau_{rz}], \qquad (5.2)$$

with $N$ the shape functions vector and $r$ the radial coordinate. Assuming, as before, a constant Jacobian determinant, we obtain

$$f_{\text{int}} = Ar \begin{bmatrix} b_1 + HRh & 0 & \frac{1}{4r}t & b_2 + HZh \\ 0 & b_2 + HZh & 0 & b_1 + HRh \end{bmatrix} \cdot \sigma, \qquad (5.3)$$

where $r$ is now (and until the end of the section) the radius at the center of the element, $t$ and $h$ are the vectors defined in (2.8), and

$$HR = \frac{z_{31}r_{42} + z_{42}r_{31}}{24Ar}, \quad HZ = -\frac{r_{31}z_{42}}{12Ar}. \qquad (5.4)$$

In the case of a square element located just next to the symmetry axis, the nodal forces obtained by (5.3) are the same as those obtained by an exact integration (with constant stresses).

If we calculate the null-space of the axisymmetrical $B$ matrix, we find hourglass modes slightly different from the plane ones:

$$g = [HR(t'r) + HZ(t'z)]t - 4HRr - 4HZz + [1 + HR(h'r) + HZ(h'z)]h, \qquad (5.5)$$

where $r$ and $z$ are the nodal axisymmetrical coordinates vectors (for more details, see [16]). But we will use the plane ones for the stabilization, since the difference is very small (see Fig. 3) and since we have noticed no improvement in the solution by using the effective cylindrical modes.

<figure><img src="image_4.png" /><figcaption>Fig. 3. The plane and the cylindrical hourglass modes.</figcaption></figure>
Note that a purely axisymmetrical study of the hourglass phenomenon can be found in Matejovic and Adamik [17] and in Matejovic [18], but it leads to a heavier and less general formulation.

## 6. An explicit scheme for time integration

To analyse dynamic problems, we have chosen an explicit integration scheme because of its simplicity of implementation in a developing computation code. But its major drawback is to be conditionally stable: the time step that can be used is bounded by a critical value.

The working sequence is as follows:

(i) the nodal accelerations are given by the equilibrium equation

$$ M \cdot \ddot{d}^n = f_{\text{ext}}^n - f_{\text{int}}^n, \qquad (6.1) $$

where the constant mass matrix is given by

$$ M = \int_{V(t)} \rho N^t \cdot N \, dV, \qquad (6.2) $$

with N the shape functions matrix.

(ii) the nodal velocities and positions are then obtained by the central difference scheme

$$ \dot{d}^{n+1/2} = \dot{d}^{n-1/2} + \Delta t^n \ddot{d}^n, \quad d^{n+1} = d^n + \Delta t^n \dot{d}^{n+1/2}. \qquad (6.3) $$

Note that the nodal accelerations are computed with a lumped mass matrix (by a row summation). This operation avoids a matrix inversion and also has the advantage of lowering the natural frequencies of the element, which are raised by the explicit integration (for more details on explicit schemes, see [19]).

The critical timestep $\Delta t_{\text{crit}}$ can be obtained from the Courant-Friedrichs-Levy condition that states that a strain wave cannot pass through the entire element during one single timestep. The different waves that are encountered in elastoplastic problems are the elastic pressure wave, the elastic shear wave and the plastic wave, whose speeds in a unidimensional medium are respectively given by (see [19])

$$ C_p = \sqrt{\frac{E}{\rho}}, \quad C_s = \sqrt{\frac{\mu}{\rho}}, \quad C_{pl} = \sqrt{\frac{E_t}{\rho}}, \qquad (6.4) $$

where $\rho$ is the current density, $E$ is Young's modulus and $E_t$ is the plastic tangent stiffness. So we see from (6.4) that the larger wave speed is $C_p$, and we can express the critical timestep as

$$ \Delta t_{\text{crit}}^n = \gamma \frac{L^n}{C_p}, \qquad (6.5) $$

where $\gamma$ is a security coefficient and $L^n$ a characteristic length (the maximum distance between two adjacent nodes, for example) of the elements at the time $t^n$.

Another approach is given by the theory of vibrations, which states that the critical timestep is

$$ \Delta t_{\text{crit}} = \frac{2}{\omega_{\text{max}}}, \qquad (6.6) $$

where $\omega_{\text{max}}$ is the maximal natural pulsation of the element.

Using the bounds established by Flanagan and Belytschko [4], we obtain the following expression:

$$ 2C_p^2 b_{il} b_{il} \le \omega_{\text{max}}^2 \le 4C_p^2 b_{il} b_{il}, \quad \Delta t_{\text{crit}} = \gamma \sqrt{\frac{\rho}{(\lambda + 2\mu)b_{il}b_{il}}}, \qquad (6.7) $$

where $\gamma$ is again a security coefficient, for which a value of 0.9 seems to be well suited.

Replacing the $b_{il}$ by their expression (2.5), we see that $b_{il}b_{il}$ is proportional to the sum of the
<figure><img src="image_5.png" /><figcaption>Fig. 4. Locus of optimal element shapes.</figcaption></figure>

element diagonals' square lengths. On that basis, it is easily shown that the square element is optimal, in the sense that it gives the larger timestep for a given area (see Fig. 4).

The Flanagan and Belytschko [4] bounds have been obtained from the one integration point quadrilateral without stabilization. We can check that the eigenvalues of the stabilized element are still comprised between those bounds, what would allow us to keep the expression (6.7) for the critical timestep. Considering a linear situation, the eigenvalues problem for the stabilized element is

$$ (K - \omega^2 M)d = 0, \quad K_{ij,j} d_{ij} = \omega^2 \frac{\rho A}{4} d_{ii}, \qquad (6.8) $$

with the linear stiffness matrix given by

$$ \begin{align} K_{ij,j} &= K_{ij,j}^{0} + K_{ij,j}^{HG}, \notag \\ K_{ij,j}^{0} &= A[\lambda b_{ij}b_{jj} + \mu(\delta_{ij}b_{kl}b_{kj} + b_{jj}b_{ij})], \tag{6.9} \\ K_{ij,j}^{HG} &= \frac{16}{3}e^2 A \mu (b_{k\kappa}b_{k\kappa}) \delta_{ij}\gamma_i\gamma_j. \notag \end{align} $$

As in [4], we can develop the eigenmodes in a linear combination of the $b_i$, but also of $h$:

$$ d_{il} = \alpha_{ik}b_{kl} + \beta_{i}h_{i}, \qquad (6.10) $$

Let us now define the following matrix: $a_{ij} = b_i^t b_j$. In the principal system of $a$, the set of equations (6.8) can be written as (values expressed in this system are noted by a superposed circumflex):

$$ \begin{gather} \frac{\rho\omega^2}{4} \alpha_{ij} = \lambda\delta_{ij}\alpha_{km}\hat{a}_{km} + \mu(\alpha_{im}\hat{a}_{jm} + \alpha_{jm}\hat{a}_{im}) + \frac{16e^2}{3}\mu\hat{a}_{nn}\hat{f}_j(\alpha_{im}\hat{a}_{km}\hat{f}_k - \beta_i), \tag{6.11} \\ \frac{\rho\omega^2}{4}\beta_i = \frac{4e^2}{3}\mu\hat{a}_{nn}(\beta_i - \alpha_{im}\hat{a}_{km}\hat{f}_k), \notag \end{gather} $$

where

$$ f_i = h_i x_{il}. \qquad (6.12) $$

For rectangular or diamond-shaped elements, the $f_i$ are zero and the two equation sets become uncoupled. The first set is the same as in [4], while the second one gives the following hourglass natural pulsation:

$$ \omega^2 = \frac{8e^2}{3} C_s^2 (a_{11} + a_{22}). \qquad (6.13) $$

Since this pulsation is proportional to the elastic shear waves speed $C_s$ given by (6.4), it is smaller than the value (6.7) obtained in [4]. Thus, the critical eigenvalue is still the one associated with the pressure mode and the critical timestep (6.7) is to be considered.
For elements of arbitrary shape, the eigenproblem (6.8) can be numerically solved. For the geometries we have tried (e.g. trapezoidal), the maximal natural pulsations found are still between the bounds given in [4] (see [16]). Thus, we keep the critical timestep (6.7) in our explicit scheme. Since the supplementary terms in the $\vec{B}$ matrix are stabilizing ones, one can consider the expression (6.7) as conservative.

## 7. Examples

The previous theoretical developments have been implemented in METAFOR [20], which is the developing large deformation analysis code of SAMCEF [21], a general purpose finite element program developed by Samtech and the L.T.A.S.

### 7.1. Uniformly loaded beam

We take here the example of a beam simply supported at both ends, as in [1]. This beam is 0.8 m long, 0.1 m high and is uniformly loaded with a $18 \times 10^5$ N/m pressure. We will use midspan symmetry. The load is applied in $t=0$ and the computation spans over 15 ms. We will consider an elastoplastic material (linear isotropic hardening) with the following physical properties:

$$ \rho = 1000 \text{ kg/m}^3, \quad E = 10^9 \text{ N/m}^2, \quad \nu = 0, \quad \sigma_v^0 = 2 \times 10^5 \text{ N/m}^2, \quad H = 1.11 \times 10^8 \text{ N/m}^2, $$

Fig. 5 shows the deformed geometry (scale 1:1) at 15 ms for one integration point elements with 'incompressible' stabilization ($e=0.5$). We see that the element keeps good behaviour (some hourglass modes tend to appear, but remain well bounded), even confronted with very large strains.

We can also compare the midpoint displacement curves obtained with different methods and meshes. We have considered two meshes: a 4 × 8 elements mesh and a finer 8 × 16 elements mesh. In Fig. 6, we see the results obtained by different methods with the 4 × 8 mesh. We see that the elements with

<figure><img src="image_6.png" /><figcaption>Fig. 5. Deformed geometry ($t=15$ ms) with reduced integration ('incompressible' method with $e=0.5$).</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>Fig. 6. Midpoint displacement curves for different methods with 4 × 8 mesh.</figcaption></figure>

'incompressible' stabilization and $e = 0.5$ and those with generalized Flanagan and Belytschko's method give a less rigid solution than the elements with $e = 1$ and the four integration points elements (these elements in fact have only one Gauss point for their volumetric part in order to avoid locking) which present little difference. If we took the four Gauss points elements as a reference, we could think that the optimal value for $e$ is 1, but let us remember that in flexion, the beam behaviour is not perfectly represented by four elements across its height. When we compare the previous results with those obtained with a 8 × 16 mesh (Fig. 7), we see that the 'incompressible' underintegrated elements with $e = 0.5$ behave even better than the four integration points elements.

So, we can say that our 'incompressible' underintegrated element (with $e = 0.5$) is 'flexurally superconvergent'. This is due to the fact that the shear is totally underintegrated, thus avoiding locking.

Now, let us compare the CPU time and number of timesteps needed by each method for the 4 × 8 mesh (see Table 1). Computations were performed on a Vax Station 3100/M76.

So, we see that the number of timesteps is quite the same for each case, but the time needed by reduced integration elements for the computation of one step is about half the time needed by the classical element. Thus, we can say that we have an element giving a better solution with less calculation time. In the following, we will fix the value of $e$ to 0.5.

<table><caption>Table 1
Uniformly loaded beam</caption><thead><tr><th>Method</th><th>No. of timesteps</th><th>CPU time (s)</th><th>Relative CPU time</th></tr></thead><tbody><tr><td>Generalized Flanagan and Belytschko (κ = 0.1)</td><td>926</td><td>49.32</td><td>0.60</td></tr><tr><td>'Incompressible' method (e = 1.0)</td><td>921</td><td>47.77</td><td>0.58</td></tr><tr><td>'Incompressible' method (e = 0.5)</td><td>921</td><td>43.71</td><td>0.53</td></tr><tr><td>Four integration points</td><td>918</td><td>82.04</td><td>1</td></tr></tbody></table>
<figure><img src="image_8.png" /><figcaption>Fig. 7. Midpoint displacement curves for 4 × 8 and 8 × 16 meshes.</figcaption></figure>

## 7.2. Taylor bar impact

We will now consider a classical benchmark test: a cylindrical copper bar (32.4 mm high and of radius 3.2 mm) impacting a rigid wall with an initial speed of 227 m/s. The bar behaviour will be computed over 80 µs with a 5 × 50 axisymmetrical elements mesh. The copper physical properties are the following:

$$ \rho = 8930 \text{ kg/m}^3, \quad E = 117 \times 10^9 \text{ GPa}, \quad \nu = 0.35, \quad \sigma_v^0 = 400 \times 10^6 \text{ MPa}, \\ H = 100.08 \times 10^6 \text{ MPa}. $$

First, we can compare the isocurves of the equivalent plastic strain over the deformed geometry obtained with the four Gauss points elements and the two types of reduced integration elements (see Fig. 8(a) and (b)). We see that the curves are corresponding quite well. If we look at Table 2, we see that the maximum equivalent plastic strain obtained with 'full' integration is a little higher than those obtained with reduced integration. This can be explained by the fact that this maximum strain is located at the bottom center of the bar and thus is better approached with four Gauss points.

We can also compare geometrical quantities as the final height and radius (see Table 2). We see that the results are quite the same for the different methods. Note that we have compared with the DYNA2D [22] solution and with the underintegrated elements without any stabilization. Here, we see that the stabilization is not really necessary, although it leads to a slightly irregular mesh.

If we compare the number of timesteps and CPU time needed by each method, we see that underintegration leads to an important time saving (more than half the 'full' integration time). We also see that DYNA2D [22] is the quickest of all, and by far, but we must note that this code is entirely
<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /><figcaption>Fig. 8. Equivalent plastic strains in the final geometry.</figcaption></figure>

<table><caption>Table 2
Taylor bar impact</caption><thead><tr><th>Method</th><th>No. of timesteps</th><th>CPU time (min:s)</th><th>Relative CPU time</th><th>Final foot radius (mm)</th><th>Final height (mm)</th><th>ε<sup>pl</sup><sub>max</sub></th></tr></thead><tbody><tr><td>Four integration points</td><td>8419</td><td>75:39</td><td>1.00</td><td>7.13</td><td>21.43</td><td>3.114</td></tr><tr><td>'Incompressible' method (e = 0.5)</td><td>7736</td><td>35:41</td><td>0.47</td><td>7.06</td><td>21.04</td><td>2.849</td></tr><tr><td>Generalized</td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Flanagan and Belytschko (κ = 0.1)</td><td>8564</td><td>39:35</td><td>0.52</td><td>6.92</td><td>21.05</td><td>2.938</td></tr><tr><td>Underintegrated without stabilization</td><td>7914</td><td>29:18</td><td>0.39</td><td>7.08</td><td>21.10</td><td>2.868</td></tr><tr><td>DYNA2D</td><td>9506</td><td>16:04</td><td>0.21</td><td>7.16</td><td>21.43</td><td>3.000</td></tr></tbody></table>

devoted to explicit resolution of dynamic problems, whereas METAFOR is more general purpose (explicit dynamic or implicit quasistatic problems). We must also say that we have encountered other problems that DYNA2D [22] cannot satisfactorily solve (insufficient stabilization) whereas METAFOR proved to be able to do so.

To determine the exact cost of the stabilization methods, we have made the same bar impact computation, but with an imposed number of 8000 equal timesteps. The results are given in Table 3. We

<table><caption>Table 3
The stabilization costs</caption><thead><tr><th>Method</th><th>Four integration points</th><th>'Incompressible' (e = 0.5)</th><th>Generalized Flanagan and Belytschko (κ = 0.1)</th><th>Underintegrated without stabilization</th></tr></thead><tbody><tr><td>CPU time (min:s)</td><td>73:25</td><td>37:07</td><td>39:44</td><td>34:30</td></tr></tbody></table>
see that underintegration allows a time saving of about 50%, and that the stabilization costs about 7.5.% of the 'underintegrated computation time'.

## 8. Conclusion

We have seen how to develop an underintegrated element without locking and with a frame invariant stabilization. We also have shown how to generalize Flanagan and Belytschko's method. This element works in modelling large elastoplastic deformations for plane or axisymmetrical structures. It allows important time savings and moreover has a better flexural behaviour than the fully integrated (in its deviatoric part) element while it gives comparable results in other loading cases.

## References

[1] D.P. Flanagan and T. Belytschko, A uniform strain hexadron and quadrilateral with orthogonal hourglass control, Internat. J. Numer. Methods Engrg. 17 (1981) 679-706.

[2] T. Belytschko and L.P. Bindeman, Assumed strain stabilization of the 4-node quadrilateral with 1-point quadrature for nonlinear problems, Comput. Methods Appl. Mech. Engrg. 88 (1991) 311-340.

[3] D.P. Flanagan and T. Belytschko, Correction of article by D.P. Flanagan and T. Belytschko, Internat. J. Numer. Methods Engrg. 19 (1983) 467-468.

[4] D.P. Flanagan and T. Belytschko, Eigenvalues and stable timesteps for the uniform strain hexhedron and quadrilateral, J. Appl. Mech. 51 (1981) 35-40.

[5] W.K. Liu and T. Belytschko, Efficient linear and nonlinear heat conduction with a quadrilateral element, J. Appl. Mech. 20 (1984) 931-948.

[6] T. Belytschko, J.S.-J. Ong, W.K. Liu and J.M. Kennedy, Hourglass control in linear and nonlinear problems, Comput. Methods Appl. Mech. Engrg. 43 (1984) 251-276.

[7] W.K. Liu, J.S.-J. Ong and R.A. Uras, Finite element stabilization matrices - a unification approach, Comput. Methods Appl. Mech. Engrg. 53 (1985) 13-46.

[8] W.K. Liu, T. Belytschko, J.S.J. Ong and S.E. Law, Use of stabilization matrices in nonlinear analysis, Engrg. Comput. 2 (1985) 47-55.

[9] T. Belytschko and W.E. Bachrach, Efficient implementation of quadrilaterals with high coarse-mesh accuracy, Comput. Methods Appl. Mech. Engrg. 54 (1986) 279-301.

[10] B.C. Koh and N. Kikushi, New improved hourglass control for bilinear and trilinear elements in anisotropic linear elasticity, Comput. Methods Appl. Mech. Engrg. 65 (1987) 1-46.

[11] W.K. Liu, T. Belytschko and J.S. Chen, Nonlinear versions of flexurally superconvergent elements, Comput. Methods Appl. Mech. Engrg. 71 (1988) 241-258.

[12] Ph. Jetteur and S. Cescotto, A mixed finite element for the analysis of large inelastic strains, Internat. J. Numer. Methods Engrg. 31 (1991) 229-239.

[13] J.C. Simo and T.J.R. Hughes, On the variational foundations of assumed strain methods, J. Appl. Mech. 53(1) (1986) 51-54.

[14] M.L. Wilkins, Calculation of elastic-plastic flow, in: B. Alder et al., eds., Methods of Computational Physics, Vol. 3 (Academic Press, New York, 1964).

[15] J.C. Nagtegaal and F.E. Veldpaus, On the implementation of finite strain plasticity equations in a numerical model, in: Pittman et al., eds., Numerical Analysis of Forming Processes (Wiley, New York, 1984) 351-371.

[16] L. Stainier, Calcul optimal des matrices d'éléments finis pur des milieux incompressibles en grandes déformations, U.Lg.-F.S.A., Travail de fin d'études, 1992.

[17] P. Matejovic and V. Adamik, A one-point integration quadrilateral with hourglass control in axisymmetric geometry, Comput. Methods. Appl. Engrg. 70 (1988) 301-320.

[18] P. Matejovic, Quadrilateral with high coarse-mesh accuracy for solid mechanics in axisymmetric geometry, Comput. Methods Appl. Mech. Engrg. 88 (1991) 241-258.

[19] T. Belytschko, An overview of semidiscretization and time integration procedures, in: Computational Methods for Transient Analysis (Elsevier, Amsterdam, 1983) 1-66.

[20] J.P. Ponthot, Mode d'emploi pour la version pilote de METAFOR, Module de caucul en grandes déformations, Report TF-24, University of Liège, LTAS-Thermomécanique, 1992.

[21] SAMTECH S.A., SAMCEF - Manuels de l'utilisateur, Liège, 1992.

[22] J.O. Hallquist, User's manual for D'YNA2D - An explicit two-dimensional hydrodynamic finite element code with interactive rezoning, Lawrence Livermore Laboratory, 1984.