A numerical investigation on Contact Mechanics applications using eight-node hexahedral elements with underintegration techniques

Abstract

In the present work the performance of finite element formulations with different reduced integration strategies is evaluated for Contact Mechanics applications. One-point quadrature and selective reduced integration are utilized here using hourglass control to suppress volumetric and shear locking for materials with incompressible plastic behavior and bending-dominated problems. A corotational formulation is adopted to deal with physically and geometrically nonlinear analysis and the generalized-α method is employed for time integration in the nonlinear dynamic range. The contact formulation is based on the penalty method, where the classical Coulomb’s law is used considering a convected coordinate system for three-dimensional friction with large deformation and finite sliding. Contact problems involving deformable and rigid bodies, as well as static and dynamic analysis, are investigated and results are analyzed considering the different underintegration formulations proposed here.

Keywords:
Contact Mechanics; Finite Element Method; Reduced Integration; Nonlinear Analysis; Elastoplasticity

Graphical Abstract

1 INTRODUCTION

Element technology is mainly concerned with the development of accurate and efficient finite element formulations, especially for large-scale and nonlinear problems. In the field of Contact Mechanics, applications frequently require high computational efforts due to the inherently nonlinear nature of the contact problem, taking into account the contact search algorithms and procedures for numerical evaluation of contact forces and stress tensor components. In this context, underintegration and stabilization techniques may be utilized to obtain computationally efficient algorithms and element formulations free of hourglass and locking instabilities.

Among the existing approaches to solve numerically the contact problem, the penalty and Lagrange multiplier methods are the most popular. In the penalty method, the impenetrability restriction for two contacting bodies is only satisfied approximately using a penalty parameter and a gap function to evaluate small interpenetrations between the so-called master and slave bodies. Sliding may be also permitted considering an interface constitutive formulation with elastoplastic characteristics (Michalowski and Mroz (1978)), such as the Coulomb friction model, and a computationally efficient algorithm for contact search is required, especially for large deformation problems. In this case, implicit methods are very convenient for the solution of low-frequency dynamic applications, although a careful procedure for the linearization of the virtual work principle is necessary. The highly nonlinear nature of the contact problem is enhanced when materials with inelastic behavior are analyzed.

A continuum framework for frictional contact with large deformation was introduced by Laursen and Simo (1993), where a convected coordinate system was developed in order to obtain a friction law with frame indifference. However, some limitations were observed when sliding movements extend over the boundaries of adjacent elements, considering that the corresponding convected coordinate systems are locally defined and, therefore, discontinuous. A frictional contact formulation for large slip was presented by Agelet de Saracibar (1997) and Chen et al. (1999) proposed a numerical formulation similar to Laursen and Simo (1993), where the convected coordinate system is maintained, but the sliding term of the contact constitutive equation is evaluated using the reference configuration to include finite sliding. Detailed information on Computational Contact Mechanics may be found in Wriggers (2006) and Laursen (2010).

A low-order finite element such as the eight-node hexahedron is a natural choice for three-dimensional problems. Unfortunately, finite element formulations with full quadrature are very inefficient for large-scale systems and may be subject to numerical instabilities due to locking phenomena for incompressible materials and bending-dominated problems. Significant improvements are obtained in terms of efficiency and accuracy by using low-order elements with underintegration numerical techniques, where reduced integration is applied on the stress/strain tensor components uniformly or selectively (see, for instance, Belytschko et al. (2014)). However, stabilization procedures are required to avoid excitation of spurious modes, which can be suppressed using hourglass control techniques, although hourglass control cannot fully remove these kinematic modes, particularly for coarse meshes. These hourglass modes are found when displacement fields are associated with spurious zero strain energy modes.

In order to deal with the incompressibility constraint, finite elements with mixed formulation may be adopted, but improvements in the element performance are restricted to constrained problems. Finite element models based on displacement formulation suffer from volumetric locking, which can be remedied using higher order elements. However, higher order elements are not recommended for large-scale computations and low-order elements are also subject to volumetric locking for incompressible or nearly incompressible materials. In this case, reduced integration is generally utilized considering selective-reduced integration or uniform-reduced integration (one-point quadrature). For eight-node hexahedral elements, selective-reduced integration refers to one-point integration over the volumetric terms and 2x2x2 quadrature points for the deviatoric terms. While selective-reduced integration leads robust solid elements, this aspect is not observed for shell elements (Belytschko et al. (2014)).

Although the incompressible behavior is unusual in linear analysis, some materials are incompressible in the plastic range. Nagtegaal et al. (1974) identified that displacement predictions are often underestimated when the finite element method is applied to elastoplastic analysis, especially for perfectly plastic materials. A mean-dilatation formulation was then proposed to improve the finite element solution in the fully plastic range. Hughes (1980) developed the so-called B-bar method using a projection method for displacement formulations, where the strain tensor is broken up into deviatoric and dilatational parts, with only the later being underintegrated (selective integration). The F-bar method was proposed by Souza Neto et al. (1996) to deal with large-strain incompressibility and Elguedj et al. (2008) presented a nonlinear F-bar formulation for large deformation employing a modified minimum potential energy principle and the B-bar method.

One of the first works on hourglass control for hexahedral elements with one-point quadrature is due to Kosloff and Frazier (1978), where h-stabilization was adopted. Later, Flanagan and Belytschko (1981) developed a stabilization scheme using a user-specified parameter to define an hourglass control based on artificial damping and artificial stiffness. In their element formulation, the stabilization matrix and hourglass-resisting force vector were obtained using anti-hourglass modes. Volumetric and shear locking were not investigated. Liu et al. (1985) proposed a stabilization matrix without user-specified parameters utilizing anti-hourglass vectors obtained from partial derivatives of the strain components with respect to the local natural coordinates. Element stabilization was accomplished by expanding the strain tensor components in Taylor series up to bilinear terms. Nevertheless, shear locking was not eliminated.

Schultz (1985) utilized an hourglass control scheme where the stress tensor components are expanded in Taylor series about the element center. The hourglass modes were suppressed by retaining the first and second order derivative terms in the Taylor series expansion. A stabilization approach based on the assumed strain method and corotational coordinate system was presented by Belytschko and Bindeman (1993) for eight-node hexahedral elements with one-point quadrature. Liu et al. (1994) proposed an underintegrated eight-node hexahedron formulation based on the model introduced by Liu et al. (1985), where shear locking instabilities were eliminated by describing the generalized strain components in a local corotational system and omitting certain nonconstant terms of the shear strain tensor.

Stainier and Ponthot (1994) presented a one-point integration formulation for large strain elastoplastic analysis using quadrilateral elements, where the stabilization scheme proposed by Flanagan and Belytschko (1981) is generalized and a modification of the assumed strain stabilization method introduced by Belytschko and Bindeman (1991) is presented. A finite element formulation involving hexahedral elements with one-point quadrature applied to elastoplastic analysis may be found in Liu et al. (1998). Reese (2005) also investigated stabilization techniques for one-point integration formulations applied to three-dimensional elastoplastic analysis with finite strains, where an optimized parameter was proposed for the stability matrix. A sensitivity study with respect to the optimized parameter for perfectly elastoplastic materials was carried out by Reese (2007) using an eight-node solid-shell finite element formulation based on reduced integration and hourglass stabilization.

A simple formulation for hexahedral elements with one-point quadrature was presented by Hu and Nagy (1997) using the multi-quadrature formulation proposed by Liu et al. (1994). The corotational coordinate system is adopted to eliminate singular modes associated with shear locking and volumetric locking was suppressed evaluating the dilatational part of the gradient matrix at the center of the element. A similar formulation was adopted by Duarte Filho and Awruch (2004) for geometrically nonlinear static and dynamic analysis of shells, plates and beams. Braun and Awruch (2008) extended that formulation to deal with highly nonlinear dynamic problems, where the Generalized-α method was utilized for time integration. A numerical model for elastoplastic analysis was implemented by Braun and Awruch (2013a, 2013b) to simulate the mechanical behavior of soils using explicit stress integration and the Modified Cam-Clay model. Extensions of the numerical model presented by Duarte Filho and Awruch (2004) were also successfully applied to composite materials (see Andrade et al. (2007)) and fluid-structure interaction (see Braun and Awruch (2009)).

To the authors’ knowledge, papers dedicated to the study of solid elements with one-point quadrature applied to three-dimensional contact problems and elastoplastic materials are scarce. Therefore, the main objective of the present work is to investigate the performance of an eight-node hexahedral element formulation using one-point quadrature, developed by Braun and Awruch (2008, 2013a, 2013b), when it is applied to contact problems with elastoplastic materials and compare its efficiency and robustness with the B-bar method, which is a selective-reduced integration technique that is widely used in the literature and also implemented into some finite element software.

The numerical model utilized here is an extension of the eight-node hexahedral finite element formulation developed by Braun and Awruch (2008, 2013a, 2013b). The contact formulation is built using the algorithms proposed by Bittencourt and Creus (1998) and Wriggers (2006), where a three-dimensional frictional contact model with finite sliding, large deformation and elastoplastic materials is considered. A corotational coordinate system is adopted to deal with geometrically and physically nonlinear applications and the contact problem is kinematically described utilizing a convected coordinate system. The penalty method is adopted to approximately enforce the impenetrability condition using an implicit scheme, where the Coulomb constitutive model and an iterative search algorithm are employed. A consistent linearization is adopted for the weak form of the equations of motion and time integration for nonlinear dynamic analysis is performed using the generalized-α method. The sparse system of equations obtained is solved using the PARDISO solver (De Coninck et al. (2016); Verbosio et al. (2017); Kourounis et al. (2018)). Contact problems involving deformable and rigid bodies, static and dynamic analyses are investigated taking into account the influence of the underintegration approaches presented here on the numerical predictions.

The present work is organized as follows: in section 2, the contact mechanics formulation utilized here is presented; in section 3, the incremental analysis and the finite element implementation are discussed; in section 4, the underintegration and stabilization techniques proposed in this work are detailed; in section 5, numerical analyses are performed and results obtained with numerical formulations proposed here are compared with predictions obtained by other authors, including comparisons with ANSYS commercial software. Finally, conclusions on the present investigation are summarized in section 6.

2 CONTACT MECHANICS FORMULATION

2.1 Contact kinematics

Consider the motion of two deformable bodies with reference configurations at t = 0 denoted by Ω^α (α = 1,2), as indicated in Figure 1, which are approaching each other and eventually come into contact under a finite deformation process. Each body is defined here as a bounded domain Ω^α ⊂ R³, with their respective boundaries Γ^α (α = 1,2) decomposed into three parts: ${\Gamma }_{\sigma }^{\alpha }$ denoting a boundary region with prescribed surface loads, ${\Gamma }_u^{\alpha }$ indicating a boundary region with prescribed displacements and ${\Gamma }_c^{\alpha }$ is the boundary region where contact forces are developed during a time interval t _i ≤ t ≤ t _f , when the bodies come into contact.

Figure 1
Definitions on contact kinematics applied to finite deformation problems.

Assuming that X ^α (α = 1,2) represents material points localized on their respective boundary surfaces Γ^α in the reference configuration and x ^α (α = 1,2) are their corresponding positions in the current configuration, the impenetrability constraint is expressed here using the closest point projection method, which leads to the following inequality equation:

$g_N=\left(x^2-{\overline{x}}^1\right)\cdot {\overline{n}}^1\ge 0$ (1)

where $x^{\alpha }=X^{\alpha }+u^{\alpha }$, uα (α = 1,2) are the displacement vectors associated with the material points $X^{\alpha }$ and g _N is the normal gap. The outward unit normal vector ${\overline{n}}^1$ is obtained considering convective coordinates and the minimization of a distance function for a master-slave pair of points, i.e.:

${\overline{n}}^1=\frac{{\overline{a}}_1^1\times {\overline{a}}_2^1}{\Vert {\overline{a}}_1^1\times {\overline{a}}_2^1\Vert }$ (2)

where ${\overline{a}}_{\beta }^1={\overline{x}}_{,\beta }^1\left({\overline{\xi }}^{\beta }\right)$ (β = 1,2) are tangent vectors defined in a local convective coordinate system related to the current configuration (see Figure 2), with boundary ${\gamma }^1$ and parametric coordinates ${\overline{\xi }}^{\beta }$, which correspond to the orthogonal projection of the slave point x2 onto the current master surface, i.e. ${\overline{x}}^1$.

In the tangential direction, the contact problem is evaluated taking into account the stick and slip conditions. In the first case, the slave point sticks to the master body and the convective coordinates ${\overline{\xi }}^{\beta }$ are invariant during the motion. Considering that g _T denotes the relative tangential displacement vector evaluated on the boundary surface of the master body, the following stick condition can be defined:

$g_T=g_{T_{\alpha }}{\overline{a}}^{\alpha }=0(\alpha =\mathrm{1,2})$ (3)

where the tangential gap components are given by:

$g_{T_{\alpha }}=\left(x^2-{\overline{x}}^1\right)\cdot {\overline{a}}_{\alpha }^1$ (4)

which implies that $x^2-{\overline{x}}^1=0$ for the stick case.

Figure 2
The convective coordinate system and the orthogonal projection method.

On the other hand, when sliding occurs in frictional contact problems, the relative velocities of the slave point x ² on the contact surface have to be integrated over the time interval t _i ≤ t ≤ t _f regarding the slip motion. The path of x ² relative to the master surface is then obtained considering that the parametric coordinates ${\overline{\xi }}^{\beta }$ of the projection point ${\overline{x}}^1$ are constantly varying. The path of x ² on the contact surface can be computed using the differential of Eq. (3) as follows:

$d{g}_T={\overline{a}}_{\alpha }^{1}d{\xi }^{\alpha }$ (5)

where the tangent vectors ${\overline{a}}_{\alpha }^1={\overline{x}}_{,\alpha }^1$ are evaluated at the instantaneous projection point ${\overline{\xi }}^{\alpha }$ and $d{\xi }^{\alpha }={\dot{\xi }}^{\alpha }dt$. The time derivative of parametric coordinates ${\dot{\overline{\xi }}}^{\alpha }$ can be obtained considering the time derivative of Eq. (4) and taking into account that Eq. (4) is null at the contact point due to orthogonality between $x^2-{\overline{x}}^1$ and ${\overline{a}}_{\alpha }^1$. This procedure leads to the following system of equations:

${\dot{\overline{\xi }}}^{\beta }={\overline{H}}^{\alpha \beta }\left[\left(v^2-{\overline{v}}^1\right)\cdot {\overline{a}}_{\alpha }^1+g_N{\overline{n}}^1\cdot {\overline{v}}_{,\alpha }^1\right](\alpha ,\beta =\mathrm{1,2})$ (6)

where $v^2={\dot{x}}^2$, ${\overline{v}}^1={\dot{\overline{x}}}^1$ and ${\overline{H}}^{\alpha \beta }={\left({\overline{H}}_{\alpha \beta }\right)}^{-1}$, with ${\overline{H}}_{\alpha \beta }$ defined as:

${\overline{H}}_{\alpha \beta }=\left[{\overline{a}}_{\alpha \beta }-g_N{\overline{b}}_{\alpha \beta }\right]$ (7)

where ${\overline{a}}_{\alpha \beta }={\overline{x}}_{,\alpha }^1\cdot {\overline{x}}_{,\beta }^1$ are the components of the metric tensor and ${\overline{b}}_{\alpha \beta }={\overline{x}}_{,\alpha \beta }^1\cdot {\overline{n}}^1$ are the components of the curvature tensor.

A Lie derivative is finally obtained, which is utilized to compute objective updates of the tangential slip vector g _T , i.e.:

$\Im g_T={\dot{\overline{\xi }}}^{\alpha }{\overline{a}}_{\alpha }^1$ (8)

2.2 Constitutive equations

When contact takes place between two bodies, contact tractions are developed on their respective contact areas, which are calculated here considering constitutive equations relating normal and tangential components with the corresponding gap functions. In the current configuration, the stress vector evaluated in the master surface may be expressed as:

$t^1\left({\overline{\xi }}_1,{\overline{\xi }}_2\right)=t_N^1{\overline{n}}^1+t_{T\beta }^1{\overline{a}}^{1\beta }(\beta =\mathrm{1,2})$ (9)

where $t_N^1$ is the normal contact pressure and $t_{T\beta }^1$ are the contact tangential components, both defined with respect to the projection point $\xi =\left({\overline{\xi }}_1,{\overline{\xi }}_2\right)$ in the local contact surface of the master object. In order to satisfy the action-reaction principle, the traction vector acting on the slave body is determined as $t^2=-t^1\left({\overline{\xi }}_1,{\overline{\xi }}_2\right)$. Notice that $t_N^1=t_N^2<0$ and the tangential components $t_{T\beta }^1$ are null for frictionless contact.

For frictional contact, the Coulomb’ law is utilized in an elastoplastic framework, where the tangential slip is decomposed into elastic (stick) and plastic (slip) parts, i.e.:$g_T=g_T^e+g_T^p$. The Coulomb plastic function may be presented as:

$\Phi =\Vert t_T^1\Vert -\mu t_N^1$ (10)

where Φ > 0 indicates slip and Φ ≤ 0 denotes stick. The tangential part of the friction traction vector is evaluated here considering a convected basis owing to finite deformation, i.e. $t_T^1=t_{T\beta }^1{\overline{a}}^{\beta }$. A flow rule for the tangential slip rate is also defined using the maximum dissipation principle, which leads to:

${\dot{g}}_T^p=\dot{\gamma }\frac{\partial \Psi }{\partial t_T^1}=\dot{\gamma }{n}_T^1$ (11)

where $\dot{\gamma }$ is the plastic parameter and $\Psi =\Vert t_T^1\Vert$ is the plastic potential, which leads to $n_T^1=t_T^1/\Vert t_T^1\Vert$.

In the present work, the penalty method is utilized to remove limitations imposed by the contact constraints to the boundary value problem. The algorithmic treatment of the contact problems is then possible by specifying constitutive equations relating contact traction vector components to the corresponding contact kinematic quantities. In this case, the penalized constitutive relations may be formulated as follows:

$t_N^1=k_N{g}_N$ (12)

where k _N is normal penalty parameter. The contact tangential traction during the stick condition is obtained as:

$t_T^1=k_T{g}_T^e$ (13)

where k _T is the tangential penalty parameter. By using Eqs. (11), (13) and the decomposition of the tangential slip rate vector, one obtains:

${\dot{g}}_T-\dot{\gamma }\frac{\left(t_T^1\right)}{\Vert t_T^1\Vert }=\frac{1}{k_T}\Im t_T^1$ (14)

where $\Im t_T^1$ is the Lie derivative of the tangential traction vector $t_T^1$.

2.3 Weak form

The Hamilton’s principle is usually adopted to describe the dynamic behavior of a continuum, which may be expressed as follows:

${\pi }_H={\displaystyle \underset{t_0}{\overset{t_f}{\int }}\delta \left(K-\pi \right)dt+}{\displaystyle \underset{t_0}{\overset{t_f}{\int }}\delta W_{d}dt}$ (15)

$K=\frac{1}{2}{\displaystyle \int \rho {\dot{u}}^T\dot{u}dv};\delta K=\frac{1}{2}{\displaystyle \int \rho \delta {\dot{u}}^T\dot{u}dv}$ (16)

$\pi ={\displaystyle \int u\left(\epsilon \right)dv}-{\displaystyle \int u^{T}b}dv-{\displaystyle \int u^{T}t}ds;\delta \pi ={\displaystyle \int \delta {\epsilon }^{{\rm T}}\sigma dv}-{\displaystyle \int \delta u^{T}b}dv-{\displaystyle \int \delta u^{T}t}ds$ (17)

$W_d=-{\displaystyle \int u^T{f}_{d}dv};\delta W_d=-{\displaystyle \int \delta u^T{f}_{d}dv}$ (18)

where K and π are the kinetic and total potential energy, respectively, U is the potential energy, V is the external force potential, W _d is the work done by the non-conserving forces fd, including damping, and $u\left(\epsilon \right)={\displaystyle \int \sigma d\epsilon }$ is the specific strain energy. The displacement variations δ u must vanish at the time limits t ₀ and t _f and also on boundary s _u , where displacements are imposed. The material is characterized by the specific mass ρ and constitutive relations must be established between the Cauchy stress tensor σ and the small strain tensor ε. In order to deal with finite deformation problems, a special care is required with respect to the kinematical description of the continuum. In this work a corotational formulation is adopted, where stress and strain components are measured according to local reference systems at the instantaneous configuration.

The semi discrete system of momentum equations is obtained taking into account that they are discrete in space but continuous in time. The space discretization is performed here considering the Bubnov-Galerkin method applied into the context of the Finite Element analysis, where displacement variations associated with the variational form assume the role of weight functions. By integrating by parts the kinetic energy variation presented in Equation (16) and taking into account the restrictions imposed on the displacement variations δ u at the time limits t ₀ and t _f , the following expression is obtained:

${\displaystyle \underset{t_0}{\overset{t_f}{\int }}\left({\displaystyle \int \delta u^T\rho \ddot{u}dv}+{\displaystyle \int \delta u^{T}c\dot{u}dv}+{\displaystyle \int {\left(\nabla \delta u\right)}^T\sigma dv-{\displaystyle \int \delta u^{T}bdv}-{\displaystyle \int \delta u^{T}tds}}\right)dt}=0$ (19)

considering that $f_d=c\dot{u}$, where c is the damping coefficient.

In addition, the solution of the weak form is subjected to boundary conditions $\overline{t}$ on s _σ and $\overline{u}$ on s _u , as well as initial condition given by:

${\displaystyle \int \delta u^T\left[u\left(0\right)-{\overline{u}}_0\right]dv}=0;{\displaystyle \int \delta u^T\left[\dot{u}\left(0\right)-{\overline{v}}_0\right]dv}=0$ (20)

where ${\overline{u}}_0$ and ${\overline{v}}_0$ are the initial displacement and velocity fields, respectively.

Contact is introduced into the variational formulation taking into account a two body contact problem, where the addition of two weak forms is considered as follows:

$\begin{array}{l}{\displaystyle \sum _{i=1}^2{\displaystyle \underset{t_0}{\overset{t_f}{\int }}\{{\displaystyle \underset{v^{(i)}}{\int }{\left(\delta u^{(i)}\right)}^T{\rho }^{(i)}{\ddot{u}}^{(i)}d{v}^{(i)}}+{\displaystyle \underset{v^{(i)}}{\int }{\left(\delta u^{(i)}\right)}^T{c}^{(i)}{\dot{u}}^{(i)}d{v}^{(i)}+{\displaystyle \underset{v^{(i)}}{\int }{\left(\nabla \delta u^{(i)}\right)}^T{\sigma }^{(i)}d{v}^{(i)}}}}}-\\ {\displaystyle \underset{v^{(i)}}{\int }{\left(\delta u^{(i)}\right)}^T{b}^{(i)}d{v}^{(i)}}-{\displaystyle \underset{s_{\sigma }^{(i)}}{\int }{\left(\delta u^{(i)}\right)}^T{\overline{t}}^{(i)}d{s}^{(i)}}-{\displaystyle \underset{s_c^{(i)}}{\int }{\left(\delta u^{(i)}\right)}^T{t}^{(i)}d{s}^{(i)}}\}dt=0\end{array}$ (21)

The equation above is also subjected to initial conditions ${\overline{u}}_0^{(i)}$ and ${\overline{v}}_0^{(i)}$ and boundary conditions ${\overline{u}}^{(i)}$ are imposed on $s_u^{(i)}$. Notice that the virtual work due to traction is decomposed into prescribed tractions ${\overline{t}}^{(i)}$ on $s_{\sigma }^{(i)}$ and tractions $t^{(i)}$ due to contact restrictions on $s_c^{(i)}$.

The contact virtual work δWc can be reduced to a single term over $s_c^{(i)}$ by considering that $t^{(2)}$ = - $t^{(1)}$, i.e.:

$\delta W_c=-{\displaystyle \underset{s_c^{(1)}}{\int }{\left(\delta u^{(1)}-\delta u^{(2)}\right)}^T{t}^{(1)}d{s}^{(i)}}$ (22)

In the frictionless case, Eq. (22) leads to:

$\delta W_c={\displaystyle \underset{s_c^{(1)}}{\int }{t}_N{g}_{N}d{s}^{(i)}}$ (23)

On the other hand, when friction contact is considered, one obtains:

$\delta W_c={\displaystyle \underset{s_c^{(1)}}{\int }\left(t_N\delta g_N+t_{T_{\alpha }}\delta {\overline{\xi }}^{\alpha }\right)d{s}^{(i)}}\left(\alpha =\mathrm{1,2}\right)$ (24)

where $\delta {\overline{\xi }}^{\alpha }$ are the contravariant components of the tangential virtual displacement vector $\delta g_T$, which is evaluated over $s_c^{(i)}$ in the surface coordinate system.

For the solution of the system of nonlinear equations given by Eq. (21), the Newton-Raphson method is adopted here, where linearization procedures must be carried out considering an incremental-iterative approach. Taking into account that a current equilibrated configuration x ^α (α = 1,2) is known at a time instant t and a solution is required at t + Δt, the following linearized problem is proposed:

$\delta W_{\mathrm{int},ext}+\delta W_c+\Delta \left(\delta W_{\mathrm{int},ext}\right)+\Delta \left(\delta W_c\right)=0$ (25)

where δW _int,ext is the internal and external virtual work, excluding contact forces, and Δ(⋅) denotes directional derivative associated with directions of Δu ^α , which is utilized to update the geometric configurations as ${\left(x^{\alpha }\right)}_{t+\Delta t}^{i+1}={\left(x^{\alpha }\right)}_{t+\Delta t}^i+{\left(\Delta u^{\alpha }\right)}^i,$with i indicating iterations within a time increment Δt.

A special attention is given here to the contact term Δ(δWc), whereas Δ(δWint,ext) is obtained considering usual procedures. The directional derivative of δWc is obtained as follows:

$\Delta \left(\delta W_c\right)=\Delta {\displaystyle \underset{s_c^{(1)}}{\int }\left(t_N\delta g_N+t_{T_{\alpha }}\delta {\overline{\xi }}^{\alpha }\right)d{s}^{(i)}}={\displaystyle \underset{{\Gamma }^{(1)}}{\int }\Delta \left(t_N\delta g_N+t_{T_{\alpha }}\delta {\overline{\xi }}^{\alpha }\right)d{\Gamma }^{(i)}}\left(\alpha =\mathrm{1,2}\right)$ (26)

where derivation of $\Delta t_N$, $\Delta \delta g_N$, $\Delta t_{T_{\alpha }}$ and $\Delta \delta {\overline{\xi }}^{\alpha }$ must be performed, which are not presented here for sake of brevity. A detailed description on linearization of the contact virtual work may be found in Wriggers (2006) and Laursen (2010).

3 FINITE ELEMENT MODEL

3.1 Discretized weak form

The finite element model is obtained here considering the semidiscrete approach, where spatial discretization is performed over the continuous weak form given by Eq. (21) using one-point quadrature hexahedral finite elements. In addition, nodal displacements are supposed to be continuous functions of time. The set of nonlinear ordinary differential equations is temporally discretized employing the generalized-α method (Chung and Hulbert (1993); Kuhl and Crisfield (1999)) in order to maintain numerical stabilization over time integration and incremental-iterative procedures are adopted taking into account the Newton-Raphson method. The discretized weak form is then presented as:

${\displaystyle \underset{t_0}{\overset{t_f}{\int }}{\displaystyle \underset{e=1}{\overset{nel}{\cup }}({\displaystyle \underset{v_e}{\int }\rho {\left(\delta u\right)}^T\ddot{u}dv}+{\displaystyle \underset{v_e}{\int }{\left(\delta u\right)}^{T}c\dot{u}dv+{\displaystyle \underset{v_e}{\int }{\left(\delta \epsilon \right)}^T\sigma dv}}}}-{\displaystyle \underset{v_e}{\int }{\left(\delta u\right)}^{T}bdv}-{\displaystyle \underset{s_{{\sigma }_e}}{\int }{\left(\delta u\right)}^T\overline{t}ds}-{\displaystyle \underset{s_{c_e}}{\int }{\left(\delta u\right)}^{T}tds})dt=0$ (27)

Finite element approximations are adopted for geometry and displacement field using hexahedral finite element interpolation functions in the semidiscrete sense, i.e.:

$x\left(X,t\right)=N\left(\xi \right)x_e\left(t\right)u\left(x,t\right)=N\left(\xi \right)u_e\left(t\right)$ (28)

where x _e and u _e are the nodal coordinates and nodal displacement vectors, respectively, which are associated with element e, and N is a row matrix containing the interpolation functions of the eight-node hexahedral element. The strain tensor is approximated at element level using$\epsilon =B{u}_e$, where $B$ is the gradient matrix.

Introducing the finite element approximations for displacement and strain fields at element level into Eq. (27), a set of nonlinear ordinary differential equations is obtained, which may be expressed in the following form:

${\displaystyle \underset{e=1}{\overset{nel}{\cup }}{M}_e{\ddot{u}}_e}+{\displaystyle \underset{e=1}{\overset{nel}{\cup }}{C}_e{\dot{u}}_e}+{\displaystyle \underset{e=1}{\overset{nel}{\cup }}{F}_e^{\mathrm{int}}\left(u_e\right)}+{\displaystyle \underset{e=1}{\overset{nel}{\cup }}{F}_e^c\left(u_e\right)}={\displaystyle \underset{e=1}{\overset{nel}{\cup }}{F}_e^{ext}}$ (29)

where M _e and C _e are the element mass and damping matrices, $F_e^{\mathrm{int}}$ is the internal force vector, $F_e^c$ is the contact force vector and $F_e^{ext}$ is the external force vector, including surface and body forces.

The time discretization is then applied to Eq. (29), where time continuity is broken into time increments Δt. In order to take into account the nonlinear behavior, linearization is adopted according to the Newton-Raphson method, which leads to the following global linearized form:

$\frac{\partial }{\partial u}\left[M{\ddot{u}}_{n+1}^i+C{\dot{u}}_{n+1}^i+F^{\mathrm{int}}\left(u_{n+1}^i\right)+F^c\left(u_{n+1}^i\right)\right]\Delta u=F_{n+1}^{ext}-M{\ddot{u}}_{n+1}^i-F^{\mathrm{int}}\left(u_{n+1}^i\right)-F^c\left(u_{n+1}^i\right)$ (30)

where the solution within the time step Δt is updated iteratively using $u_{n+1}^{i+1}=u_{n+1}^i+\Delta u$. The dynamic analysis is carried out employing an implicit algorithm, considering that displacements $u_n$ and velocities ${\dot{u}}_n$ are available at the beginning of the time step. Although the Newmark’s method is unconditionally stable for linear problems, numerical instabilities are usually observed for applications in the nonlinear range. Therefore, the generalized-α method is adopted here for time integration. Static analysis can be also performed with Eq. (30) by omitting the inertial and damping terms. In this case, time is considered in a fictitious sense associated with the load steps. Additional details on the time integration algorithm utilized in this work may be found in Braun and Awruch (2008, 2013a, 2013b).

3.2 Nonlinear analysis with the corotational reference system

In the present work, nonlinear problems are analyzed considering the corotational approach, where the strain and stress tensors are evaluated in a reference system locally defined at the element quadrature points. In a corotational kinematical description, the nonlinear problem can be presented in rate form taking into account the small strain hypothesis and an objective rate of the Cauchy stress tensor. Assuming large rotations, the motion of a continuum can be decomposed into rigid body motion and pure deformation. Once the spatial discretization is fine enough, this motion decomposition can be performed at element level, where the pure deformation component is a small quantity with respect to the element dimensions.

Strains are update in the corotational system using the following expression:

${\widehat{\epsilon }}^{n+1}={\widehat{\epsilon }}^n+\Delta \widehat{\epsilon }$ (31)

where $\widehat{\epsilon }$ is the corotational strain tensor and $\Delta \widehat{\epsilon }$ is the corotational strain increment, which is obtained at element level using mid-point integration (Hughes and Winget (1980)) as follows:

$\Delta \widehat{\epsilon }={\displaystyle \underset{n}{\overset{n+1}{\int }}\widehat{d}dt}\doteq \frac{1}{2}\left[\frac{\partial \Delta {\widehat{u}}^{def}}{\partial {\widehat{x}}_{n+1/2}}+{\left(\frac{\partial \Delta {\widehat{u}}^{def}}{\partial {\widehat{x}}_{n+1/2}}\right)}^T\right]$ (32)

where $\widehat{d}$ is the corotational strain rate tensor, $\Delta {\widehat{u}}^{def}$ is the corotational displacement increment associated with the deformation component of the total displacement increment $\Delta \widehat{u}=\Delta {\widehat{u}}^{def}+\Delta {\widehat{u}}^{rot}$, with $\Delta {\widehat{u}}^{rot}$ indicating the displacement component related to pure rotation. The increment of deformation displacements is calculated using $\Delta {\widehat{u}}^{def}=R_{n+1/2}\Delta u^{def}={\widehat{x}}_{n+1}-{\widehat{x}}_n$, where Δudef is the global deformation displacement increment, and the geometric configurations ${\widehat{x}}_n$, ${\widehat{x}}_{n+1/2}$ and ${\widehat{x}}_{n+1}$, which are related to a specific finite element e, are obtained from the following transformations:

${\widehat{x}}_n=R_n{x}_n;{\widehat{x}}_{n+1/2}=R_{n+1/2}{x}_{n+1/2};{\widehat{x}}_{n+1}=R_{n+1}{x}_{n+1}$ (33)

where R _n , R _n+1/2 and R _n+1 are orthogonal transformation matrices performing rotations from the global coordinate system to the corotational coordinate system, while x _n , x _n+1/2 and x _n+1 are element geometric configurations defined in the global coordinate system.

Although the corotational Cauchy stress tensor σ is frame invariant, an objective rate tensor is required for stress updates during the incremental-iterative integration procedures. In this work, the Truesdell rate tensor is adopted, which may be expressed as:

${\dot{\widehat{\sigma }}}^{TR}=\dot{\widehat{\sigma }}-\widehat{L}\widehat{\sigma }-\widehat{\sigma }{\widehat{L}}^T+\widehat{\sigma }\cdot tr\widehat{d}$ (34)

where $\widehat{L}$ is the spatial velocity gradient tensor evaluated in the corotational system.

The corotational Cauchy stress tensor is updated using an incremental form of Eq. (34), where the stress increment is calculated considering an incremental constitutive formulation for elastoplastic materials, i.e.:

$\Delta \widehat{\sigma }=D^{ep}\Delta \widehat{\epsilon }$ (35)

where $\Delta \widehat{\epsilon }$ is given by Eq. (32) and the constitutive tensor is defined as:

$D^{ep}=D^e-\frac{D^e{a}_g{a}_f^T{D}^e}{A+a_f^T{D}^e{a}_g}\Delta \widehat{\epsilon }$ (36)

where D ^e is the fourth-order linear elastic tensor, A is the hardening parameter and $a_f$ and $a_g$ are flow vectors associated with the yield function f and plastic potential g, respectively. The constitutive equation is integrated over time considering an explicit algorithm proposed by Owen and Hinton (1980).

In the present formulation, the tangent stiffness matrix and the internal force vector are evaluated in the corotational system as follows:

${\widehat{K}}_e^{\tan }={\displaystyle \underset{{\widehat{v}}_e}{\int }{\widehat{B}}^T\left(D^{ep}+{\widehat{D}}^{geo}\right)\widehat{B}dv};{\widehat{F}}_e^{\mathrm{int}}={\displaystyle \underset{{\widehat{v}}_e}{\int }{\widehat{B}}^T\widehat{\sigma }dv}$ (37)

where $\widehat{B}$ is the gradient matrix evaluated in the corotational system and ${\widehat{D}}^{geo}$ is the geometric constitutive tensor related to Truesdell rate tensor components. In order to solve the system of linearized equations of motion (Eq. 30), the tangent stiffness matrix and the internal force vector must be obtained in the global coordinate system as follows:

$K_e^{\tan }=R^T{\widehat{K}}_e^{\tan }R;F_e^{\mathrm{int}}=R^T{\widehat{F}}_e^{\mathrm{int}}$ (38)

An elastoplastic framework is also adopted for integration of the interface constitutive equation in the case of friction contact, where stick and slip (slide) conditions are identified considering the classical Coulomb’s friction law (see Eq. 10). The slide condition is associated with yield, where Φ > 0 and irreversible slip occurs. On the other hand, the stick condition is obtained when Φ < 0, where a small and reversible (elastic) slip is permitted. The problem is solved taking into account a rate formulation, which is integrated over time steps considering elastic trial and return mapping.

An elastic predictor is first considered, where a linear path and reversible behavior are assumed within the time increment for nodes on the contact surface, i.e.:

${\left(t_N\right)}_{n+1}=k_N{\left(g_N\right)}_{n+1};{\left(t_{T_{\alpha }}^{trial}\right)}_{n+1}={\left(t_{T_{\alpha }}\right)}_n+k_T{\left({\overline{a}}_{\alpha \beta }\right)}_{n+1}\left[{\overline{\xi }}_{n+1}^{\beta }-{\overline{\xi }}_n^{\beta }\right](\alpha ,\beta =\mathrm{1,2})$ (39)

where ${\overline{a}}_{\alpha \beta }$ denotes the metric tensor components while ${\overline{\xi }}_n^{\beta }$ and ${\overline{\xi }}_{n+1}^{\beta }$ indicate the convective coordinates at t _n and t _n+1 .

Notice that the tangential traction in the previous equation cannot be evaluated when a contact point slides over the element boundaries, considering that a convective coordinate system is locally defined for every finite element. Consequently, the load steps must be chosen carefully to avoid this condition and the target domain must be extended in order to allow that the slave node projection to be slightly outside the element domain. After that, it is assumed that the slave node is projected on the adjacent element and the numerical process continues. Another solution is presented by Chen et al. (1999), where the position vector of the contact point is mapped to the reference configuration from its parametric coordinate system, such that the increment of relative displacement can be computed.

The slip function is then checked and the tangential traction components are defined according to the following conditions:

${\Phi }_{n+1}^{trial}=\Vert {\left(t_T^{trial}\right)}_{n+1}\Vert -\mu \left|{\left(t_N\right)}_{n+1}\right|\Rightarrow \{\begin{array}{l}{\left(t_{T_{\alpha }}\right)}_{n+1}={\left(t_{T_{\alpha }}^{trial}\right)}_{n+1}if{\Phi }_{n+1}^{trial}\le 0\hfill \\ {\left(t_{T_{\alpha }}\right)}_{n+1}=\mu \left|{\left(t_N\right)}_{n+1}\right|\frac{{\left(t_{T_{\alpha }}^{trial}\right)}_{n+1}}{\Vert {\left(t_T^{trial}\right)}_{n+1}\Vert }if{\Phi }_{n+1}^{trial}>0\hfill \end{array}$ (40)

It is important to notice that the normal contact traction component is given according to the global coordinate system, while the tangential contact traction components are written in terms of local convective coordinate axes.

3.3 Contact stiffness matrix and contact residual vector

A contact element formulation is adopted here to deal with node-to-surface contacts, where four-node master elements are defined considering every hexahedral element face subject to contact with a slave node defined by x ². A five-node contact element is then formulated using the following interpolation equation for the master surface:

$x^1=\stackrel{⌣}{N}\left(\xi \right)x_e^1$ (41)

where $\stackrel{⌣}{N}$ is a row matrix containing the interpolation functions of the four-node quadrilateral element expressed in terms of local surface coordinates ξ ₁ and ξ ₂. Vector $x_e^1$ contains the nodal coordinates of a master element e.

In order to obtain the normal gap g _N (see Eq. 1), ${\overline{x}}^1$ and ${\overline{n}}^1$ must be determined locally for every slave node x ² using closest point projection, where the convective coordinates ${\overline{\xi }}_1$ and ${\overline{\xi }}_2$ associated with x ² are computed considering the following system of nonlinear equations:

$\left[x^2-x^1\left(\xi \right)\right]\cdot x_{,\alpha }^1=\left[x^2-\stackrel{⌣}{N}\left(\xi \right)x_e^1\right]\cdot {\stackrel{⌣}{N}}_{,\alpha }\left(\xi \right)x_e^1=0(\alpha =\mathrm{1,2})$ (42)

which is solved iteratively using the Newton-Raphson method (see Wriggers (2006)). The closest point $\overline{\xi }=({\overline{\xi }}_1,{\overline{\xi }}_2)$ is then obtained and the gap function (Eq. 1) is evaluated. When contact is identified (g _N ≤ 0), the corresponding slave node generates contributions to the weak form.

The contact virtual work δW _c can be discretized using the general form given by Eq. (24) as follows:

$\delta W_c(\overline{\xi })={\displaystyle \underset{s_{ce}}{\int }\left[t_N\delta g_N+{\left(\delta g\right)}^T{t}_T\right]ds}\approx {\displaystyle \sum _{se}^{ns}\left[t_{Nse}(\overline{\xi })\delta g_{Nse}(\overline{\xi })+{\left(\delta g_{Tse}(\overline{\xi })\right)}^T{t}_{Tse}\right]A_{se}}$ (43)

which leads to:

$\delta W_c(\overline{\xi })={\displaystyle \sum _{se}^{ns}{A}_{se}{\left(\delta g(\overline{\xi })\right)}^T{F}_{se}^c}$ (44)

where ns is the number of slave nodes, A _se is the area surrounding the slave node, which is assumed to be included in the normal penalty parameter k _N , δ g is the displacement variation related to the slave nodes on the master surface and $F_{se}^c$ is the contact residual vector, expressed by:

$F_{se}^c=k_N{g}_{Nse}{\stackrel{⌣}{N}}_{se}^T+t_{Tse\alpha }{\left({\stackrel{⌣}{D}}_{se}^{\alpha }\right)}^T(\alpha =\mathrm{1,2})$ (45)

with:

$\begin{array}{l}{\stackrel{⌣}{N}}_{se}=\left[\begin{array}{ccccc}{\overline{n}}^1& -{\stackrel{⌣}{N}}_1(\overline{\xi }){\overline{n}}^1& -{\stackrel{⌣}{N}}_2(\overline{\xi }){\overline{n}}^1& -{\stackrel{⌣}{N}}_3(\overline{\xi }){\overline{n}}^1& -{\stackrel{⌣}{N}}_4(\overline{\xi }){\overline{n}}^1\end{array}\right];\\ {\stackrel{⌣}{D}}_{se}^{\alpha }={\overline{H}}^{\alpha \beta }\left[{\stackrel{⌣}{T}}_{\beta }-g_{Nse}{\stackrel{⌣}{N}}_{\beta }\right](\alpha ,\beta =\mathrm{1,2})\end{array}$ (46)

where:

$\begin{array}{l}{\stackrel{⌣}{T}}_{\beta }=\left[\begin{array}{ccccc}{\overline{a}}_{\beta }& -{\stackrel{⌣}{N}}_1(\overline{\xi }){\overline{a}}_{\beta }& -{\stackrel{⌣}{N}}_2(\overline{\xi }){\overline{a}}_{\beta }& -{\stackrel{⌣}{N}}_3(\overline{\xi }){\overline{a}}_{\beta }& -{\stackrel{⌣}{N}}_4(\overline{\xi }){\overline{a}}_{\beta }\end{array}\right];\\ {\stackrel{⌣}{N}}_{\beta }=\left[\begin{array}{ccccc}0& -{\stackrel{⌣}{N}}_{\mathrm{1,}\beta }(\overline{\xi }){\overline{n}}^1& -{\stackrel{⌣}{N}}_{\mathrm{2,}\beta }(\overline{\xi }){\overline{n}}^1& -{\stackrel{⌣}{N}}_{\mathrm{3,}\beta }(\overline{\xi }){\overline{n}}^1& -{\stackrel{⌣}{N}}_{\mathrm{4,}\beta }(\overline{\xi }){\overline{n}}^1\end{array}\right]\end{array}$ (47)

The tensor components ${\overline{H}}^{\alpha \beta }$ are calculated using Eq. (7) and the friction traction components $t_{Tse\alpha }$ are obtained considering the incremental procedure defined by Eqs. (39) and (40).

The contact contribution to the global tangent stiffness is obtained at element level using linearization of Eq. (43) based on the Newton-Raphson method applied in the context of the penalty formulation (see also Eq. 26), i.e.:

$\Delta \left(\delta W_c\right)={\displaystyle \sum _{se}^{nelc}{A}_{se}{\left(\delta g_{se}\right)}^T{K}_{se}^c\Delta u_{se}}$ (48)

which leads to the contact stiffness matrix:

$K_{se}^c=K_{Nse}^c+K_{Tse}^c$ (49)

with:

$K_{Nse}^c=k_N{\stackrel{⌣}{N}}_{se}^T{\stackrel{⌣}{N}}_{se}+k_N{g}_N\left[{\stackrel{⌣}{N}}_{\alpha }^T{D}_{se}^{\alpha }+{\overline{a}}^{\beta \alpha }{\stackrel{⌣}{T}}_{\alpha }^T\left({\stackrel{⌣}{N}}_{\beta }-{\stackrel{⌣}{D}}_{se}^{\gamma }\left({\overline{n}}^1\cdot {\overline{a}}_{\beta ,\gamma }\right)\right)\right](\alpha ,\beta ,\gamma =\mathrm{1,2})$ (50)

where ${\overline{a}}^{\beta \alpha }$ are components of the inverse of the metric tensor, with components denoted by ${\overline{a}}_{\beta \alpha }$. Notice that the normal gap can be quite large during the first iterations, when a large incremental step is employed, which may lead to divergence in the Newton-Raphson process. Therefore, only the first term in the previous equation is employed in the first iterations and the terms multiplied by gN are added later.

The tangential contribution to the tangent stiffness matrix ($K_{Tse}^c$) may be decomposed into two parts, where the first tangent stiffness matrix is written as:

$\begin{array}{l}{K}_{Tse1}^c=t_{T\alpha }{H}^{\alpha \eta }[\left({\stackrel{⌣}{T}}_{\eta \beta }^T+{\stackrel{⌣}{T}}_{\beta \eta }^T+{\widehat{T}}_{\eta \beta }^T\right)D_{se}^{\beta }+{\left(D_{se}^{\beta }\right)}^T\left({\stackrel{⌣}{T}}_{\eta \beta }+{\stackrel{⌣}{T}}_{\beta \eta }+{\widehat{T}}_{\eta \beta }\right)-\\ \left({\overline{a}}_{\eta }\cdot {\overline{a}}_{\beta ,\gamma }+{\overline{a}}_{\beta }\cdot {\overline{a}}_{\eta ,\gamma }+{\overline{a}}_{\eta ,\beta }\cdot {\overline{a}}_{\gamma }-g_N{\overline{n}}^1\cdot {\overline{a}}_{\eta ,\beta \gamma }\right){\left({\stackrel{⌣}{D}}_{se}^{\beta }\right)}^T{\stackrel{⌣}{D}}_{se}^{\gamma }-\\ {\stackrel{⌣}{E}}^T{\stackrel{⌣}{E}}_{,\eta }-{\stackrel{⌣}{E}}_{,\eta }^T\stackrel{⌣}{E}-g_N{\widehat{N}}_{\eta \beta }^T{\stackrel{⌣}{D}}_{se}^{\beta }-g_N{\left({\stackrel{⌣}{D}}_{se}^{\beta }\right)}^T{\widehat{N}}_{\eta \beta }](\alpha ,\beta ,\gamma =\mathrm{1,2})\end{array}$ (51)

with the following vectors and matrices:

$\begin{array}{l}{\stackrel{⌣}{T}}_{\alpha \beta }=\left[\begin{array}{ccccc}0& -{\stackrel{⌣}{N}}_{\mathrm{1,}\beta }(\overline{\xi }){\overline{a}}_{\alpha }& -{\stackrel{⌣}{N}}_{\mathrm{2,}\beta }(\overline{\xi }){\overline{a}}_{\alpha }& -{\stackrel{⌣}{N}}_{\mathrm{3,}\beta }(\overline{\xi }){\overline{a}}_{\alpha }& -{\stackrel{⌣}{N}}_{\mathrm{4,}\beta }(\overline{\xi }){\overline{a}}_{\alpha }\end{array}\right];\\ {\widehat{T}}_{\alpha \beta }=\left[\begin{array}{ccccc}{\overline{a}}_{\alpha ,\beta }& -{\stackrel{⌣}{N}}_1(\overline{\xi }){\overline{a}}_{\alpha ,\beta }& -{\stackrel{⌣}{N}}_2(\overline{\xi }){\overline{a}}_{\alpha ,\beta }& -{\stackrel{⌣}{N}}_3(\overline{\xi }){\overline{a}}_{\alpha ,\beta }& -{\stackrel{⌣}{N}}_4(\overline{\xi }){\overline{a}}_{\alpha ,\beta }\end{array}\right];\\ {\widehat{N}}_{\alpha \beta }=\left[\begin{array}{ccccc}0& -{\stackrel{⌣}{N}}_{\mathrm{1,}\alpha \beta }(\overline{\xi }){\overline{n}}^1& -{\stackrel{⌣}{N}}_{\mathrm{2,}\alpha \beta }(\overline{\xi }){\overline{n}}^1& -{\stackrel{⌣}{N}}_{\mathrm{3,}\alpha \beta }(\overline{\xi }){\overline{n}}^1& -{\stackrel{⌣}{N}}_{\mathrm{4,}\alpha \beta }(\overline{\xi }){\overline{n}}^1\end{array}\right];\\ \stackrel{⌣}{E}=\left[\begin{array}{ccccc}I& -{\stackrel{⌣}{N}}_1(\overline{\xi })I& -{\stackrel{⌣}{N}}_2(\overline{\xi })I& -{\stackrel{⌣}{N}}_3(\overline{\xi })I& -{\stackrel{⌣}{N}}_4(\overline{\xi })I\end{array}\right];\\ {\stackrel{⌣}{E}}_{,\alpha }=\left[\begin{array}{ccccc}0& -{\stackrel{⌣}{N}}_{\mathrm{1,}\alpha }(\overline{\xi })I& -{\stackrel{⌣}{N}}_{\mathrm{2,}\alpha }(\overline{\xi })I& -{\stackrel{⌣}{N}}_{\mathrm{3,}\alpha }(\overline{\xi })I& -{\stackrel{⌣}{N}}_{\mathrm{4,}\alpha }(\overline{\xi })I\end{array}\right]\end{array}$ (52)

where I is a 3x3 identity matrix and 0 is a 3x3 null matrix.

The second tangent stiffness matrix is defined for each frictional state, i.e. stick or slip condition. By using $\Delta {\xi }^{\beta }={\overline{\xi }}_{n+1}^{\beta }-{\overline{\xi }}_n^{\beta }$, the second stiffness for the stick condition is given by:

$K_{Tse2}^c=k_T{\left(D_{se}^{\alpha }\right)}^T\left\{{\overline{a}}_{\beta \alpha }{D}_{se}^{\beta }-\Delta {\xi }^{\beta }\left[{\stackrel{⌣}{T}}_{\beta \alpha }+{\stackrel{⌣}{T}}_{\alpha \beta }-({\overline{a}}_{\alpha ,\zeta }\cdot {\overline{a}}_{\beta }+{\overline{a}}_{\alpha }\cdot {\overline{a}}_{\beta ,\zeta })D_{se}^{\zeta }\right]\right\}(\alpha ,\beta ,\zeta =\mathrm{1,2})$ (53)

while for the slip condition:

$\begin{array}{l}{K}_{Tse2}^c=\mu {\left(D_{se}^{\alpha }\right)}^T\{-k_N\frac{t_{T_{\alpha }}^{trial}}{\Vert t_T^{trial}\Vert }{\stackrel{⌣}{N}}_{se}+\frac{\left|t_N\right|}{\Vert t_T^{trial}\Vert }[k_T{\overline{a}}_{\beta \alpha }{D}_{se}^{\beta }-\\ k_T\Delta {\xi }^{\beta }\left({\stackrel{⌣}{T}}_{\beta \alpha }+{\stackrel{⌣}{T}}_{\alpha \beta }-({\overline{a}}_{\alpha ,\zeta }\cdot {\overline{a}}_{\beta }+{\overline{a}}_{\alpha }\cdot {\overline{a}}_{\beta ,\zeta })D_{se}^{\zeta }\right)]+\\ \frac{\left|t_N\right|}{{\Vert t_T^{trial}\Vert }^3}{t}_{T_{\alpha }}^{trial}{t}_T^{tria{l}^{\beta }}[-{\stackrel{⌣}{P}}_{\beta }+k_T\Delta {\xi }^{\gamma }\left({\stackrel{⌣}{T}}_{\beta \gamma }+{\widehat{T}}_{\gamma \beta }\right)+\\ D_{se}^{\zeta }\left(t_T^{trial}\cdot {\overline{a}}_{\beta ,\zeta }-k_T{\overline{a}}_{\zeta \beta }-k_T\Delta {\xi }^{\gamma }({\overline{a}}_{\gamma ,\zeta }\cdot {\overline{a}}_{\beta }+{\overline{a}}_{\gamma }\cdot {\overline{a}}_{\beta ,\zeta })\right)]\}(\alpha ,\beta ,\gamma ,\zeta =\mathrm{1,2})\end{array}$ (54)

where:

${\stackrel{⌣}{P}}_{\alpha }=\left[\begin{array}{ccccc}0& -{\stackrel{⌣}{N}}_{\mathrm{1,}\alpha }(\overline{\xi })t_T^{trial}& -{\stackrel{⌣}{N}}_{\mathrm{2,}\alpha }(\overline{\xi })t_T^{trial}& -{\stackrel{⌣}{N}}_{\mathrm{3,}\alpha }(\overline{\xi })t_T^{trial}& -{\stackrel{⌣}{N}}_{\mathrm{4,}\alpha }(\overline{\xi })t_T^{trial}\end{array}\right]$ (55)

Notice that the first tangent stiffness matrix due to frictional contact is symmetric, while the second part is no longer symmetric in both frictional states.

The contact between a deformable body and a rigid surface is obtained as a special case of the previous development. The rigid surface is assumed to be the master surface, described by a parametric equation, and its motion is considered to be prescribed. Thus, the degrees of freedom of the master surface vanish in the previous equations and the stiffness matrix turn out to be the same presented by Bittencourt and Creus (1998). In other words, only the rows and columns related to the degrees of freedom of the slave nodes are computed.

4 UNDERINTEGRATION AND STABILIZATION TECHNIQUES

In the present work, two underintegration techniques are utilized to suppress shear and volumetric locking effects from the eight-node hexahedral element formulation, especially for materials with incompressible or nearly incompressible behavior. A uniform reduced integration procedure with hourglass control is adopted using one-point quadrature for both the deviatoric and volumetric parts of the strain tensor. In addition, the B-bar method is also employed here to avoid volumetric locking for elastoplastic analysis involving materials with the incompressibility constraint.

One-point quadrature may be used for reduced integration considering that shape functions and derivatives are evaluated at the center of the element. Unfortunately, the so-called hourglass modes may be excited and some numerical procedures are necessary to stabilize the element formulation. In order to alleviate volumetric locking, the gradient matrix B is decomposed into two parts as follows:

$B\left(\xi ,\eta ,\zeta \right)=\tilde{B}\left(0\right)+\widehat{B}\left(\xi ,\eta ,\zeta \right)$ (56)

where $\tilde{B}\left(0\right)$ is the gradient matrix corresponding to the volumetric terms of the strain tensor, which is evaluated at the element center, and $\widehat{B}\left(\xi ,\eta ,\zeta \right)$ is the gradient matrix corresponding to the deviatoric terms, which must be expanded using Taylor series at the center of the element up to the bilinear terms. Equation (56) is then rewritten as:

$B\left(\xi ,\eta ,\zeta \right)=\overline{B}(0)+{\widehat{B}}_{,\xi }(0)\xi +{\widehat{B}}_{,\eta }(0)\eta +{\widehat{B}}_{,\zeta }(0)\zeta +2{\widehat{B}}_{,\xi \eta }(0)\xi \eta +2{\widehat{B}}_{,\eta \zeta }(0)\eta \zeta +2{\widehat{B}}_{,\xi \zeta }(0)\xi \zeta$ (57)

where $\overline{B}(0)$ includes the sum of the volumetric and deviatoric parts of the gradient matrix, both obtained from one-point quadrature. By performing the same procedure with the stress tensor σ, one obtains:

$\sigma \left(\xi ,\eta ,\zeta \right)=\overline{\sigma }(0)+{\widehat{\sigma }}_{,\xi }(0)\xi +{\widehat{\sigma }}_{,\eta }(0)\eta +{\widehat{\sigma }}_{,\zeta }(0)\zeta +2{\widehat{\sigma }}_{,\xi \eta }(0)\xi \eta +2{\widehat{\sigma }}_{,\eta \zeta }(0)\eta \zeta +2{\widehat{\sigma }}_{,\xi \zeta }(0)\xi \zeta$ (58)

The internal virtual work at element level is evaluated as follows:

$\delta W_e^{\mathrm{int}}=\delta u_e^T{\displaystyle \underset{{\Omega }_e}{\int }{B}^T(\xi ,\eta ,\zeta )}\sigma (\xi ,\eta ,\zeta )dv$ (59)

By substituting Eqs. (57) and (58) into Eq. (59), the following expression is obtained:

$\begin{array}{l}\delta W_e^{\mathrm{int}}=\delta u_e^T[{\overline{B}}^T(0)\overline{\sigma }(0)+\frac{1}{3}{\widehat{B}}_{,\xi }^T(0){\widehat{\sigma }}_{,\xi }(0)+\frac{1}{3}{\widehat{B}}_{,\eta }^T(0){\widehat{\sigma }}_{,\eta }(0)\eta +\frac{1}{3}{\widehat{B}}_{,\zeta }^T(0){\widehat{\sigma }}_{,\zeta }(0)+\\ \frac{1}{9}{\widehat{B}}_{,\xi \eta }^T(0){\widehat{\sigma }}_{,\xi \eta }(0)+\frac{1}{9}{\widehat{B}}_{,\eta \zeta }^T(0){\widehat{\sigma }}_{,\eta \zeta }(0)+\frac{1}{9}{\widehat{B}}_{,\xi \zeta }^T(0){\widehat{\sigma }}_{,\xi \zeta }(0)]{\Omega }_e\end{array}$ (60)

where Ω_e is the current element volume. Notice that the first term on the right-hand side of Eq. (60) is identified with the one-point quadrature internal virtual work, while the remaining terms, which are also evaluated at the element center, are associated with the element stabilization.

In order to eliminate the spurious singular modes, Hu and Nagy (1997) recommend adding an hourglass-resisting force f ^hg to the element internal force vector f ^int:

$f^{\mathrm{int}}=f^{o-p}+f^{hg}=\left({{\rm K}}^{o-p}+{{\rm K}}^{stab}\right)u$ (61)

where f ^o-p is the internal force vector referring to one-point quadrature, which is obtained using the constitutive relation, i.e:

$f^{o-p}={\overline{{\rm B}}}^{{\rm T}}(0)\overline{\sigma }(0){\Omega }_e=\left[{\overline{{\rm B}}}^{{\rm T}}(0)D\overline{{\rm B}}(0){\Omega }_e\right]u={{\rm K}}^{o-p}u$ (62)

where D is the constitutive tensor (D ^e or D ^ep ). On the other hand, the hourglass-resisting force vector f ^hg is expressed as:

$\begin{array}{l}{f}^{hg}=[\frac{1}{3}{\widehat{{\rm B}}}_{,\xi }^{{\rm T}}(0){\widehat{\sigma }}_{,\xi }(0)+\frac{1}{3}{\widehat{{\rm B}}}_{,\eta }^{{\rm T}}(0){\widehat{\sigma }}_{,\eta }(0)+\frac{1}{3}{\widehat{{\rm B}}}_{,\zeta }^{{\rm T}}(0){\widehat{\sigma }}_{,\zeta }(0)+\\ \frac{1}{9}{\widehat{{\rm B}}}_{,\xi \eta }^{{\rm T}}(0){\widehat{\sigma }}_{,\xi \eta }(0)+\frac{1}{9}{\widehat{{\rm B}}}_{,\eta \zeta }^{{\rm T}}(0){\widehat{\sigma }}_{,\eta \zeta }(0)+\frac{1}{9}{\widehat{{\rm B}}}_{,\xi \zeta }^{{\rm T}}(0){\widehat{\sigma }}_{,\xi \zeta }(0)]{\Omega }_e={{\rm K}}^{stab}U\end{array}$ (63)

where K ^stab is the stabilization stiffness matrix, which is obtained considering the stabilization matrix E proposed by Hu and Nagy (1997) in order to determine the stress tensor derivatives, i.e.:

${\widehat{\sigma }}_{,\xi }={\rm E}{\widehat{\epsilon }}_{,\xi };{\widehat{\sigma }}_{,\eta }={\rm E}{\widehat{\epsilon }}_{,\eta };{\widehat{\sigma }}_{,\zeta }={\rm E}{\widehat{\epsilon }}_{,\zeta };{\widehat{\sigma }}_{,\xi \eta }={\rm E}{\widehat{\epsilon }}_{,\xi \eta };{\widehat{\sigma }}_{,\eta \zeta }={\rm E}{\widehat{\epsilon }}_{,\eta \zeta };{\widehat{\sigma }}_{,\xi \zeta }={\rm E}{\widehat{\epsilon }}_{,\xi \zeta }$ (64)

with E given as:

${\rm E}={\left[\begin{array}{cc}e& 0\\ 0& e\end{array}\right]}_{(6x6)}withe=\left[\begin{array}{ccc}2{\mu }^{*}& 0& 0\\ 0& 2{\mu }^{*}& 0\\ 0& 0& 2{\mu }^{*}\end{array}\right]$ (65)

where μ* = μ for elastic materials, considering that μ is the Lamé constant associated with shear. For elastoplastic materials, this alternative leads to loss of element stiffness when a certain stress limit is reached. In order to improve the element performance, Reese (2005) proposed an optimized parameter (μ* = μ ^opt) for the stabilization matrix. This parameter (μ ^opt) leads to the smallest value that provides sufficient strength to inhibit hourglass modes and may be expressed as:

${\mu }^{opt}=\mu \frac{H}{E+H}$ (66)

where E and H are the Young and hardening moduli, respectively, considering that H is evaluated at the onset of plastification. Reese (2005) mentions that for situations where μ* is much smaller than μ ^opt, a too soft behavior may be observed, while much higher values may lead to locking effects. In addition, the author also states that the influence of μ* on the solution decreases noticeably when the number of elements is increased. However, μ ^opt is null for perfectly elastoplastic materials, considering that H = 0 in this case. Consequently, the stabilization matrix is removed. As a solution, a percentage of 𝜇 can be adopted for perfectly elastoplastic materials, i.e.:

${\mu }^{*}=k\mu 0\le k\le 1$ (67)

Reese (2007) indicates that μ* should be chosen to be as small as possible when problems with perfectly elastoplastic materials are analyzed, e.g. 0.001% of the Young’s modulus. Nevertheless, an optimum value for μ* is not clearly defined in this case.

The element developed so far is only free of volumetric locking. In order to remove shear locking, shear components of the strain tensor must be described using an orthogonal corotational coordinate system and linearly interpolated in a single direction of the reference system as follows:

$\begin{array}{l}{\epsilon }_{xy}\left(\xi ,\eta ,\zeta \right)={\epsilon }_{xy}\left(0\right)+{\widehat{\epsilon }}_{xy,\zeta }\left(0\right)\zeta \\ {\epsilon }_{yz}\left(\xi ,\eta ,\zeta \right)={\epsilon }_{yz}\left(0\right)+{\widehat{\epsilon }}_{yz,\xi }\left(0\right)\xi \\ {\epsilon }_{xz}\left(\xi ,\eta ,\zeta \right)={\epsilon }_{xz}\left(0\right)+{\widehat{\epsilon }}_{xz,\eta }\left(0\right)\eta \end{array}$ (68)

Consequently, one can observe that ${\widehat{B}}_{xy,\zeta }$, ${\widehat{B}}_{yz,\xi }$ and ${\widehat{B}}_{xz,\eta }$, which correspond to the deviatoric strain components ${\widehat{\epsilon }}_{xy}$, ${\widehat{\epsilon }}_{yz}$ and ${\widehat{\epsilon }}_{xz}$, are the only gradient matrix derivatives that are not null.

Finally, the gradient matrix obtained with reduced integration B(0) is replaced by the uniform gradient matrix B’(0) proposed by Flanagan and Belytschko (1981):

${B^{\prime }}_i=\frac{1}{{\Omega }_e}{\displaystyle \underset{{\Omega }_e}{\int }{B}_i}\left(\xi ,\eta ,\zeta \right)d{\Omega }_e\left(i=\mathrm{1,2,...,8}\right)$ (69)

The selective reduced integration technique can be also adopted to overcome the incompressibility constraint, where only the volumetric terms of the stiffness matrix are evaluated using a lower order integration rule. Although this simple procedure is very effective for isotropic elastic materials, it is not the same when elastoplastic materials are considered. An alternative procedure was proposed by Hughes (1980) using a decomposition of the gradient matrix into dilatational and deviatoric terms. This method is usually referred as the B-bar method.

In the B-bar method, the gradient sub-matrices B _i are first decomposed into dilatational and deviatoric parts as follows:

$B_i=B_i^{dil}+B_i^{dev}\left(i=\mathrm{1,2,...,8}\right)$ (70)

where:

$B_i^{dil}=\frac{1}{3}\left[\begin{array}{ccc}{B}_1& B_2& B_3\\ B_1& B_2& B_3\\ B_1& B_2& B_3\\ 0& 0& 0\\ 0& 0& 0\\ 0& 0& 0\end{array}\right]$ (71)

This dilatational part must be improved in order to obtain an effective formulation for incompressible applications. The sub-matrix $B_i^{dil}$ is then replaced by ${\overline{B}}_i^{dil}$, which is defined as:

${\overline{B}}_i^{dil}=\frac{1}{3}\left[\begin{array}{ccc}{\overline{B}}_1& {\overline{B}}_2& {\overline{B}}_3\\ {\overline{B}}_1& {\overline{B}}_2& {\overline{B}}_3\\ {\overline{B}}_1& {\overline{B}}_2& {\overline{B}}_3\\ 0& 0& 0\\ 0& 0& 0\\ 0& 0& 0\end{array}\right]$ (72)

Equation (70) is also rewritten as follows:

${\overline{B}}_i={\overline{B}}_i^{dil}+B_i^{dev}\left(i=\mathrm{1,2,...,8}\right)$ (73)

where $B_i^{dev}=B_i-B_i^{dil}$ are fully integrated and ${\overline{B}}_i$ are given explicitly by:

${\overline{B}}_i=\left[\begin{array}{ccc}{B}_5& B_6& B_8\\ B_4& B_7& B_8\\ B_4& B_6& B_9\\ B_2& B_1& 0\\ 0& B_3& B_2\\ B_3& 0& B_1\end{array}\right]$ (74)

where:

$\begin{array}{l}{B}_4=\left({\overline{B}}_1-B_1\right)/3;B_5=B_1+B_4;B_6=\left({\overline{B}}_2-B_2\right)/3\\ B_7=B_2+B_6;B_8=\left({\overline{B}}_3-B_3\right)/3;B_9=B_3+B_8\end{array}$ (75)

In the present work, the sub-matrix components ${\overline{B}}_i$ are locally evaluated using:

${\overline{B}}_i={\displaystyle \underset{{\Omega }_e}{\int }{\overline{B}}_{i}d\Omega }/{\displaystyle \underset{{\Omega }_e}{\int }d\Omega }$ (76)

5 NUMERICAL APPLICATIONS

5.1 Upsetting of a cylindrical billet

This example is a classic problem presented by Taylor and Becker (1983), where a three-dimensional cylindrical billet is compressed by a perfectly rough rigid surface in order to reduce its length to 64% of the initial length. The billet is characterized with a diameter of 20 mm and length of 30 mm. Only 1/8 of the billet is modeled due to symmetry, with 480 hexahedral elements and 671 nodes, as shown in Figure 3.

The problem is simulated assuming von Mises yield criterion for the material and 200 load steps were performed using the Newton-Raphson method. Material and numerical parameters utilized in this example are presented in Table 1.

Figure 3
Initial mesh for the cylindrical billet.

Table 1
Material and contact properties for billet upsetting analysis.

Young’s modulus - E	200 kN/mm2
Linear hardening modulus - H	0.3 kN/mm2
Yield stress - σ y	0.7 kN/mm2
Poisson coefficient - ν	0.3
Normal penalty parameter - k N	3.5 x 105 kN/mm
Tangential penalty parameter - k T	3.5 x 104 kN/mm

Load-deflection curves computed for both element formulations are presented in Figure 4, where results obtained by Taylor and Becker (1983), Bittencourt (1994) and Simo and Ju (1989) are also presented. One can observe that by using B-bar formulation, the numerical analysis was unable to finish, reaching a 61.12% upset and showing some jumps when the side of the billet folded onto the rigid surface. On the other hand, when the one-point quadrature formulation was applied, the analysis was unable to converge, even for the first load steps, where many elements reached the yield stress. Thus, although the present material is not perfectly elastoplastic, convergence was only obtained when k = 0.12 was utilized in Eq. (67) to compute the stabilization matrix. Notice that, in general, the computed load-deflection curves showed a good agreement with results obtained by Taylor and Becker (1983) and Simo and Ju (1989), while Bittencourt (1994) presented a smaller load due to the finite element mesh employed.

Figure 4
Load-displacement curves for billet upsetting analysis.

Deformed mesh configurations on the symmetry plane, corresponding to the last converged load step, are shown in Figure 5 according to the underintegration formulations utilized in this work.

Figure 5
Deformed billet mesh configurations for the last converged load step.

It is observed that both underintegration techniques utilized in this work showed numerical difficulties due to excessive element distortion. Other mesh configurations were tested and even with a higher number of elements it was not possible to improve the present solution considerably. This drawback may be circumvented using an arbitrary Lagrangian-Eulerian kinematic formulation, which can describe large deformation problems with reduced finite element distortion.

5.2 Copper cylinder impact on rigid surface

This benchmark example is frequently known as the ‘Taylor’s bar problem’ and consists of a small cylinder with initial length of 32.4 mm and initial radius of 3.2 mm, which impacts against a flat surface with an initial velocity of 227 m/s. The material properties of the copper cylinder are presented in Table 2.

Table 2
Material properties of the copper cylinder.

Young’s modulus - E	117 x 109 N/m2
Linear hardening modulus - H	100 x 106 N/m2
Yield stress - σ y	400 x 106 N/m2
Poisson coefficient - ν	0.35
Specific mass - ρ	8930 kg/m3

The dynamic analysis is performed during 80 μs with Δt = 1.6 x 10^-8 s and spectral radius r_α = 1.0. Due to symmetry, a quarter of the bar is modeled with 1080 hexahedral elements and 1517 nodes. A frictionless contact is assumed with normal penalty parameter equal to 10⁹ N/m.

The last converged configuration obtained by each element formulation is shown in Figure 6. The numerical models with full-integration and B-bar formulation completed the 80 μs analysis, while the analysis using elements with one-point quadrature was unable to converge, even for the initial time stages. It is possible to see that the standard hexahedral element, i.e. with full integration, showed volumetric locking (Figure 6a). In addition, the hexahedral element formulation with one-point quadrature was unable to control the numerical instabilities that arose in this case (Figure 6c).

Figure 6
Converged configurations for the Taylor impact problem.

The final cylinder radius and final cylinder length obtained here are compared with those obtained by other authors in Table 3. Notice that results obtained using standard formulation showed wrong solutions due to locking effects, while results computed here with the B-bar formulation indicated a good agreement with those obtained by other authors.

Table 3
Taylor impact problem - final cylinder radius and length.

Author	Element type	Formulation	Final radius	Final length
Present work	8-node hex.	Full integration	4.94 mm	20.29 mm
Present work	8-node hex.	B-bar	7.25 mm	21.44 mm
Present work	8-node hex.	One-point quadrature	-	-
Bonet and Burton (1998)	8-node hex.	-	6.95 mm	-
Bonet and Burton (1998)	4-node tet.	Full integration	5.55 mm	-
Aymone et al. (2001)	8-node hex.	B-bar	7.10 mm	21.44 mm
Zienkiewicz et al. (1998)	8-node hex.	Mixed	7.11 mm	21.48 mm
Zienkiewicz et al. (1998)	4-node quad.	Mixed	7.11 mm	21.48 mm

5.3 Longitudinal impact between two bars

This numerical example consists of a frontal impact between two elastic bars. Note that the contact is the only nonlinearity present in the problem. Both bars have the same dimensions, with a cross section of 1 x 1 cm and 5 cm length. It is initially considered that bar A is in motion with a speed of 1 cm/s and bar B is stationary. For the discretization of the bars, 10 hexahedral eight-node finite elements are used for each bar along the longitudinal direction, as shown in Figure 7. Physical properties and computational parameters adopted in this problem are shown in Table 4.

Figure 7
Impact problem between two bars: geometry and mesh configurations.

It is well known from the classical theory of elasticity (see Timoshenko and Goodier (1970)) that the velocity at the moment of impact between the two bars is equal to 0.5 cm/s, producing waves that propagate towards the end of the bars and later reflect towards the contact surface. The impact duration is equal to 1 second, when the energy of bar A is transferred to bar B.

Table 4
Impact problem between two bars: physical properties and computational parameters.

Young’s modulus - E	1 N/cm2
Poisson coefficient - ν	0.0
Specific mass - ρ	0.01 kg/cm3
Normal penalty parameter - k N	4 N/cm
Time step - Δt	5 x 10-3 s
Total time - T	5 s
Spectral radius - r α	1.0

Figure 8 shows the contact surface displacements over time for each element formulation. One can see that both formulations are in agreement and consistent with the predictions obtained by the classical theory of elasticity.

Figure 8
Two bars under frontal impact: contact surface displacements over time.

5.4 Contact between two beams

This example is presented by Chen et al. (1999), which consists of two cantilever beams in contact undergoing large deformations. Both beams are 20 mm long, 2 mm wide and 1 mm thick, where the left end of the upper beam and the right end of the lower beam are fixed, as shown in Figure 9. A prescribed displacement of 30 mm is applied at the fixed end of the upper beam, while both beams are prevented from moving along the depth direction.

Figure 9
Geometrical and boundary conditions for the contact problem between two beams.

It is assumed that both beams have an elastic behavior with Young’s modulus (E) equal to 200 GPa and Poisson coefficient (ν) equal to 0.30. The contact is performed with friction coefficients μ = 0.0 and 0.5. Some oscillations were observed in the present results when small load steps were employed, especially for the frictional case. Thus, the analysis is performed with 500 loads steps. The normal and tangential penalty parameters employed here are equal to 10⁴ and 10³ N/mm, respectively.

The problem is initially discretized using a finite element mesh of 20 x 2 x 1 (length x height x width) hexahedral elements with full integration in order to compare the present predictions with the results obtained by Chen et al. (1999), who also employed a standard hexahedral finite element formulation. Comparisons were also carried out using ANSYS with SOLID45 elements and no extra displacement shapes.

The relationship between the prescribed displacement and the reaction force at the fixed end is presented in Figure 10. One can observe that the present results are very similar to the ones presented by Chen et al. (1999) and by ANSYS. Note that the reaction forces for the frictional case were larger than those observed for the frictionless case, considering that the beams lost contact before reaching the total prescribed displacement δ = 30 mm.

Analyses using hexahedral elements with B-bar formulation and one-point quadrature were also performed but they presented smaller reaction forces than those shown by the case with full integration. Results computed here are shown in Figure 11 and compared with those obtained using the software ANSYS with SOLID185 elements, which also presents the B-bar formulation and uniform reduced integration.

Figure 10
Prescribed displacement vs. reaction force for the contact problem between two beams: results with full integration elements.

The results presented in Figures 10 and 11 probably indicate the occurrence of shear locking for the element models with full integration and B-bar formulation. It is also observed that the ANSYS element SOLID185, using uniform reduced integration, presented a softer behavior when compared with the results obtained by the element formulation using one-point quadrature, which is employed in the present work.

Figure 11
Prescribed displacement vs. reaction force for the contact problem between two beams: results with underintegrated elements.

Therefore, the mesh configuration was modified to suppress shear locking effects and the problem was finally discretized using a mesh of 60 x 4 x 1 (length x height x width) hexahedral elements with B-bar formulation and one-point quadrature. Figure 12 shows the reaction force at the fixed end as a function of the prescribed displacement for both underintegration techniques, where results obtained using the software ANSYS are also presented. One can see that results presented here were very similar and the shear locking was successfully suppressed. Notice that the ANSYS element SOLID185 with uniform reduced integration showed again a softer behavior, although it is in a good agreement with the predictions obtained with the one-point quadrature element proposed in the present work.

Figure 12
Prescribed displacement vs. reaction force for the contact problem between two beams: results with underintegrated elements and refined mesh.

The deformed configuration for μ = 0.5, when the prescribed vertical displacement (δ) is equal to 15 mm, is shown in Figure 13.

Figure 13
Deformed configurations for the contact problem between two beams: results with underintegrated elements and refined mesh.

5.5 Soil-pile interaction

A single pile surrounded by a three-dimensional soil layer and submitted to axial or lateral loads is analyzed in this example. It is assumed that the pile and the soil have elastic behavior with Young’s moduli equal to E_p = 2 x 10⁷ kPa and E_s = 2 x 10⁴ kPa, respectively, and Poisson coefficients equal to ν_p = 0.30 and ν_s = 0.45, respectively. The square pile is 10 m long and 0.5 m in width, which is embedded in a 17 m deep soil layer.

The computational domain is presented in Figure 14, where only half of the spatial domain is modeled using 6760 hexahedral elements and 8004 nodes due to symmetry. Nodes located on the base of the computational domain are fixed in all directions, while nodes on the lateral walls are constrained only in the normal direction.

A rough frictional contact is assumed for the interface elements, where soil and pile can open a gap between them but no sliding is allowed. The normal and tangential penalty parameters employed here are both equal to 1.2 x 10⁵ kN/m.

Results referring to head settlement, when the pile is subjected to an axial load of 2,500 kN, and lateral deflection, when the pile is subjected to a lateral load of 216 kN at its head, are presented in Table 5. Note that due to symmetry, half of the previous load values are applied in the numerical model. Results computed here are compared with numerical results presented by Trochanis et al. (1991) and elastic solutions provided by Poulos and Davis (1980).

Figure 14
Mesh configuration and geometry for the soil-pile interaction problem.

Table 5
Pile head settlement and lateral deflection.

Author	Head settlement [cm]	Lateral deflection [cm]
Present work with B-bar formulation	1.73	0.62
Present work with one-point quadrature	1.73	0.71
Trochanis et al. (1991)	1.66	0.63
Poulos and Davis (1980)	1.70	0.80

Notice that the computed head settlements are quite similar, while the lateral deflection results showed a good agreement between the present work with B-bar formulation and Trochanis et al. (1991). On the other hand, numerical models showed a stiffer behavior for lateral deflection when compared to the analytical solution proposed by Poulos and Davis (1980).

The final deformed shape for the pile subjected to a lateral load is shown in Figure 15, where one can observe the soil-pile separation near the pile head for both underintegration techniques.

Figure 15
Final deformed shape for laterally loaded pile (amplification factor of 100).

The soil behavior away from the pile, when subjected to axial or lateral load, is also analyzed in this problem. Therefore, the soil surface displacements, u for the axial load model and w for the lateral load model, are normalized with respect to the head pile displacements U _p and W _p, respectively. The distance from the pile position, s, is normalized with respect to the pile width, D. Figure 16 shows a comparison considering the vertical displacements computed on the soil surface nodes as a function of the distance from the pile, while Figure 17 shows a comparison on the horizontal displacements along the loading line and normal to the loading line.

Figure 16
Soil surface settlements for the soil-pile interaction.

Figure 17
Soil-pile interaction - horizontal displacements on the soil surface: (a) along the loading line and (b) normal to the loading line.

Figures 16 and 17 show a good agreement between results obtained here and those obtained by Trochanis et al. (1991) and by Poulos and Davis (1980). Note that the horizontal displacements normal to the loading line are smaller and decrease faster than those computed along the loading line.

5.6 Soil-structure interaction

The present example seeks to evaluate the dynamic interaction of a shallow foundation resting on a three-dimensional soil layer with a cyclic load applied at the column free end. Geometrical characteristics of the present application are presented in Figure 18a and due to symmetry, only half of the spatial domain is modeled with 1948 eight-node hexahedral elements and 2561 nodes, as shown in Figure 18b.

Figure 18
Soil-structure interaction problem: (a) geometrical characteristics and (b) initial mesh configuration.

A rough frictional contact between the soil-structure interface is assumed and, therefore, the interface elements can open a gap but no sliding is allowed. The normal and tangential penalty parameters employed here are both equal to 4.0 x 10³ kN/m.

The soil is modeled with the Drucker-Prager yield criterion, assuming coincidence at the outer edges with the Mohr-Coulomb surface, while the structure is considered as an elastic material (see Table 6). Notice that the Young’s modulus for the soil is a function of the initial confining pressure (σ ₃ in kN/m²), which is evaluated at the center of the finite element using a soil specific weight equal to 16 kN/m².

Table 6
Physical parameters adopted for the soil-structure interaction problem.

Soil	Young’s modulus - E	40000 x (σ3/100)0.86 kN/m2
	Poisson’s ratio - ν	0.33
	Angle of friction - ϕ	42°
	Cohesion - c	1.2 kN/m2
	Specific mass - ρ	1.631 x 103 kg/m3
Structure	Young’s modulus - E	30 x 106 kN/m2
	Poisson’s ratio - ν	0.20
	Specific mass - ρ	2.548 x 103 kg/m3

Nodes located on the base of the computational domain are fixed in all directions, while nodes on the lateral walls are constrained only in the normal direction, just as the structure nodes that are on the symmetry plane.

A dynamic analysis is performed over 10 s with Δt = 5.0 x 10^-3 s, where the axial load description is shown in Figure 19 and the horizontal load is defined as F_H = 2sin(2πt) kN.

Figure 19
Load description adopted for the soil-structure interaction problem.

Two points are used for displacement measurements: point A, corresponding to the point where the loads are applied and point B, which corresponds to a point on the soil-structure interface, which is located at a central point underneath the footing in the symmetry axis.

Results computed here are compared with those obtained using the finite element commercial software ANSYS. The soil is modeled using SOLID65 elements with extra displacement shapes included, which support the classic Drucker-Prager model, while the structure is modeled using SOLID185 elements with the B-bar formulation. Both elements are defined by eight nodes having three degrees of freedom at each node.

Horizontal and vertical displacements evaluated at point B over time are plotted in Figures 20 and 21, respectively, while the horizontal displacement response at point A is shown in Figure 22. It is important to mention that ANSYS did not converge using the Newmark’s method and it was necessary to use the Generalized HHT-α method, where a numerical dissipation is introduced, similar to the Generalized-α method adopted in the present algorithm.

Figure 20
Horizontal displacement response at point B for the soil-structure problem.

Figure 21
Vertical displacement response at point B for the soil-structure problem.

Figures 20, 21 and 22 show that it was necessary to employ the Generalized-α method to stabilize the time integration process, with r _α = 0.8 for the analysis with B-bar formulation and r _α = 0.7 for the analysis with one-point quadrature. It can be seen that results obtained with ANSYS were similar to those measured in the present work with a small difference after 5 s. One can see that the ANSYS formulation showed a large amount of numerical dissipation, which strongly dissipates vertical displacement oscillations, especially after t = 4 s.

In order to evaluate the structure response resting on a rigid base, two other configurations were tested: structure resting on rigid contact surface and structure with base nodes fixed in all directions. The horizontal displacements measured at point A for both models are presented in Figure 23, where one can notice that predictions with B-bar formulation are in agreement with ANSYS results, while the predictions with one-point quadrature showed slightly larger displacements when compared with ANSYS and B-bar results.

As expected, horizontal displacement values obtained with these base models were very small when compared with those obtained considering a soil base. The maximum values obtained for the structure resting on a rigid base were at least 200 times smaller than those obtained with the ANSYS soil-structure model, as shown in Figures 22 and 23. This difference demonstrates the importance of considering soil-structure interaction when the structure is resting on a flexible soil.

Figure 22
Horizontal displacement response at point A for the soil-structure problem.

Figure 23
Horizontal displacement response at point A for the structure resting on rigid base.

6 CONCLUSIONS

Numerical simulations employing underintegration techniques associated with the eight-node hexahedral finite element were performed in this work to evaluate their performance when applied to contact problems. It was observed that both techniques showed, in general, very similar responses and presented a good agreement with numerical and analytical predictions presented by other authors. The B-bar element formulation adopted here obtained a good performance for all the numerical problems tested; however, numerical difficulties were observed when the finite element mesh was subject to high distortions, such as those presented in the billet upsetting analysis. The B-bar element formulation obtained good results for a cylinder impacting on rigid wall, considering that it is a volumetric locking-free element. However, the shear locking issue was identified when the formulation was applied to bending problems, such as the contact problem between two beams. On the other hand, the one-point quadrature element formulation presented numerical difficulties when employed to some elastoplastic problems, which was unable to control instabilities that appeared in the Taylor’s bar impact test. It was observed that the uniform reduced integration element formulation requires finer meshes when applied to elastoplastic applications in order to capture the plastic front. Nevertheless, the one-point quadrature element showed a numerical response very similar to that obtained with the B-bar formulation for elastic or elastoplastic problems with low element distortion. It was also observed that the Generalized-α method was important to maintain a stable time integration process in highly nonlinear dynamic analyses. Although results concerning numerical efficiency were not presented in this work, hexahedral element formulations with one-point quadrature usually lead to important reductions in the processing time, which is a very important feature for highly nonlinear 3D contact problems.

Acknowledgments

The authors would like to thank the National Council for Scientific and Technological Development (CNPq, Brazil) and Brazilian Federal Agency for Support and Evaluation of Graduate Education for the financial support.

References

• Agelet de Saracibar, C. (1997). A new frictional time integration algorithm for large slip multi-body frictional contact problems. Computer Methods in Applied Mechanics and Engineering 142: 303-334.
• Andrade, L. G., Awruch, A. M., Morsch, I. B. (2007). Geometrically nonlinear analysis of laminate composite plates and shells using the eight-node hexahedral element with one point integration. Composites Structures 79: 571-580, 2007.
• Aymone, J. L. F., Bittencourt, E., Creus, G. J. (2001). Simulation of 3D metal-forming using an arbitrary Lagrangian-Eulerian finite element method. Journal of Materials Processing Technology 110(2): 218-232.
• Belytschko, T., Bindeman, L. P. (1991). Assumed strain stabilization of the 4-node quadrilateral with 1-point quadrature for nonlinear problems. Computer Methods in Applied Mechanics and Engineering 88: 311-340.
• Belytschko, T., Bindeman, L. P. (1993). Assumed strain stabilization of the eight node hexahedral element. Computer Methods in Applied Mechanics and Engineering 88: 311-340.
• Belytschko, T., Liu, W. K., Moran, B., Elkhodary, K. I. (2014). Nonlinear Finite Elements for Continua and Structures, 2nd ed., John Wiley & Sons, Chichester.
• Bittencourt, E. (1994). Treatment of contact impact problem in large deformation by the finite element method. Ph.D. Thesis (in Portuguese), Universidade Federal do Rio Grande do Sul, Brazil.
• Bittencourt, E., Creus, G. J. (1998). Finite element analysis of three-dimensional contact and impact in large deformation problems. Computers and Structures 69: 219-234.
• Bonet, J., Burton, A. J. (1998). A simple average nodal pressure tetrahedral element for incompressible and nearly incompressible dynamic explicit applications. Communications in Numerical Methods in Engineering 14(5): 437-449.
• Braun, A. L., Awruch, A. M. (2008). Geometrically non-linear analysis in elastodynamics using the eight-node finite element with one-point quadrature and the generalized-α method. Latin American Journal of Solids and Structures 5: 17-45.
• Braun, A. L., Awruch, A. M. (2009). A partitioned model for fluid-structure interaction problems using hexahedral finite elements with one-point quadrature. International Journal for Numerical Methods in Engineering 79: 505-549.
• Braun, A. L., Awruch, A. M. (2013a). An efficient model for numerical simulation of the mechanical behavior of soils. Part 1: theory and numerical algorithm. Soils and Rocks 36: 159-169.
• Braun, A. L., Awruch, A. M. (2013b). An efficient model for numerical simulation of the mechanical behavior of soils. Part 2: applications. Soils and Rocks 36: 171-182.
• Chen, X., Nakamura, K., Mori, M., Hisada, T. (1999). Finite element analysis for large deformation frictional contact problems with finite sliding. JSME International Journal Series A 42(2): 201-208.
• Chung, J., Hulbert, G. M. (1993). A time integration algorithm for structural dynamics with improved numerical dissipation: the generalized-α method. Journal of Applied Mechanics 60: 371-375.
• De Coninck, A., De Baets, B., Kourounis, D., Verbosio, F., Schenk, O., Maenhout, S., Fostier, J. (2016). Needles: toward large-scale genomic prediction with marker-by-environment interaction. Genetics 203(1): 543-555.
• Duarte Filho, L. A., Awruch, A. M. (2004). Geometrically nonlinear static and dynamic analysis of shells and plates using the eight-node hexahedral element with one-point quadrature. Finite Elements in Analysis and Design 40: 1297-1315.
• Elguedj, T., Bazilevs, Y., Calo, V. M., Hughes, T. J. R. (2008). B and F projection methods for nearly incompressible linear and non-linear elasticity and plasticity using higher-order NURBS elements. Computer Methods in Applied Mechanics and Engineering 197: 2732-2762.
• Flanagan, D. P., Belytschko, T. (1981). A uniform strain hexahedron and quadrilateral with orthogonal hourglass control. International Journal for Numerical Methods in Engineering 17: 679-706.
• Hu, Y. K., Nagy, L. I. (1997). A one-point quadrature eight-node brick element with hourglass control. Computers and Structures 65: 893-902.
• Hughes, T. J. R. (1980). Generalization of selective integration procedures to anisotropic and nonlinear media. International Journal for Numerical Methods in Engineering 15: 1413-1418.
• Hughes, T., Winget, J. (1980). Finite rotations effects in numerical integration of rate constructive equations arising in large deformation analysis. International Journal for Numerical Methods in Engineering 15: 1862-1867.
• Kosloff, D., Frazier, G. A. (1978). Treatment of hourglass patterns in low order finite elements codes. International Journal for Numerical and Analytical Methods in Geomechanics 2: 57-72.
• Kourounis, D., Fuchs, A., Schenk, O. (2018). Toward the next generation of multiperiod optimal power flow solvers. IEEE Transactions on Power Systems 33(4): 4005-4014.
• Kuhl, D., Crisfield, M. A. (1999). Energy-conserving and decaying algorithms in non-linear structural dynamics. International Journal for Numerical Methods in Engineering 45: 569-599.
• Laursen, T. A., Simo, J. C. (1993). A continuum-based finite element formulation for the implicit solution of multibody, large deformation frictional contact problems. International Journal for Numerical Methods in Engineering 36: 3451-3485.
• Laursen. T. A. (2010). Computational Contact and Impact Mechanics, 2nd ed., Springer, Berlin.
• Liu, W. K., Guo, Y., Tang, S., Belytschko, T. (1998). A multiple-quadrature eight-node hexahedral finite element for large deformation elastoplastic analysis. Computer Methods in Applied Mechanics and Engineering 154: 69-132.
• Liu, W. K., Hu, Y. K., Belytschko, T. (1994). Multiple quadrature underintegrated finite elements. International Journal for Numerical Methods in Engineering 37: 3263-3289.
• Liu, W. K., Ong, J. S. S., Uras, R. A. (1985). Finite element stabilization matrices - a unification approach. Computer Methods in Applied Mechanics and Engineering 53: 13-46.
• Michalowski, R., Mroz, Z. (1978). Associated and non-associated sliding rules in contact friction problems. Archives of Mechanics 30: 259-276.
• Nagtegaal, J. C., Parks, D. M., Rice, J. R. (1974). On numerical accurate finite element solutions in the fully plastic range. Computer Methods in Applied Mechanics and Engineering 4: 153-177.
• Owen, D. R. J., Hinton, E. (1980). Finite Elements in Plasticity: Theory and Practice, Pineridge Press, Swansea.
• Poulos, H. G., Davis, E. H. (1980). Pile foundation analysis and design. John Wiley and Sons, New York.
• Reese, S. (2005). On a physically stabilized one point finite element formulation for three-dimensional finite elasto-plasticity. Computer Methods in Applied Mechanics and Engineering 194: 4685-4715.
• Reese, S. (2007). A large deformation solid-shell concept based on reduced integration with hourglass stabilization. International Journal for Numerical Methods in Engineering 69: 1671-1716.
• Schultz, J. C. (1985). Finite element hourglassing control. International Journal for Numerical Methods in Engineering 21: 1039-1048.
• Simo, J. C., Ju, J. W. (1989). On continuum damage-elastoplasticity at finite strains. Computational mechanics 5(5): 375-400.
• Souza Neto, E. A., Peric, D., Dutko, M., Owen, D. R. J. (1996). Design of simple low order finite elements for large strain analysis of nearly incompressible solids. International Journal of Solids and Structures 33: 3277-3296.
• Stainier, L., Ponthot, J. P. (1994). An improved one-point integration method for large strain elastoplastic analysis. Computer Methods in Applied Mechanics and Engineering 118: 163-177.
• Taylor, L. M., Becker, E. B, (1983). Some computational aspects of large deformation, rate-dependent plasticity problems. Computer methods in applied mechanics and engineering 41(3): 251-277.
• Timoshenko, S. P., Goodier, J. N. (1970). Theory of Elasticity. McGraw-Hill.
• Trochanis, A. M., Bielak, J., Christiano, P. (1991). Three-dimensional nonlinear study of piles. Journal of Geotechnical Engineering 117(3): 429-447.
• Verbosio, F., De Coninck, A., Kourounis, D., Schenk, O. (2017). Enhancing the scalability of selected inversion factorization algorithms in genomic prediction. Journal of Computational Science 22: 99-108.
• Wriggers, P. (2006). Computational Contact Mechanics, 2nd ed., Springer, Berlin.
• Zienkiewicz, O. C., Rojek, J., Taylor, R. L., Pastor, M. (1998). Triangles and tetrahedra in explicit dynamic codes for solids. International Journal for Numerical Methods in Engineering 43(3): 565-583.

Edited by

Editor:

Marco L. Bittencourt

Publication Dates

• Publication in this collection
  18 June 2021
• Date of issue
  2021

History

• Received
  08 Feb 2021
• Reviewed
  06 Apr 2021
• Accepted
  06 May 2021
• Published
  10 May 2021