An efficiency comparison of numerical implementation approaches of hyperelastic constitutive model in ABAQUS/Standard

The main goal of the article is to compare a computational efficiency of different implementation of a hyperelastic material model in the ABAQUS/Standard v. 6.14 [1]. The software offers basically two ways to proceed, namely UMAT and UHYPER user subroutines [2]. The procedures are employed to implement an isotropic, compressible neoHookean material model [3]. Corresponding stored-energy function, constitutive equations and the consistent tangent operator are presented. These are essential to program the subroutines. Some theoretical and numerical aspects of different implementation approaches are discussed. As examples, a tube under axial compression and a contact problem of disc are considered. On the basis of obtained results, selected aspects of computational efficiency and quality of solutions are compared.


Introduction
The main goal of the article is to compare a computational efficiency of different implementation of a hyperelastic material model. The ABAQUS/Standard FEM software [1] offers basically two ways to proceed. The first approach is to employ user subroutine UHYPER [2], which can be used to define an isotropic model. It requires that a storedenergy function and its derivatives be defined with respect to the strain invariants.
The other way, more general one, is to define a material's constitutive model using subroutine UMAT [2]. Basically, the procedure requires definition of the Jacobian matrix necessary in the FEM and the Cauchy stress tensor in proper vector notation. Similar to the UHYPER, a user may specify material behavior dependent on field variables or states ones. Moreover, the solver allows to proceed with non-symmetric matrices.
Followed by a formal verification of the implemented procedures, more complex boundary value problems are considered, i.e. a tube under axial compression and a contact problem of disc in a plane strain deformation. On the basis of obtained results, selected aspects of computational efficiency and quality of solutions are compared.
where S is the first Piola-Kirchhoff stress tensor and F defines a 'deformation gradient' as D ( ) defines an actual configuration of the body with respect to an initial one.
, the mapping  is oriented-preserving and locally invertible. In order to avoid interpenetration of matter, it is necessary to put some restrictions on ( ) W F , see [5]: From now on, we make the assumption that stored-energy function is isotropic, i.e.
where 1 The relationship between the Cauchy stress tensor σ , the first and second Piola-Kirchhoff ones, namely The Kirchhoff stress tensor J  τ σ is work conjugated to the symmetric rate of deformation tensor 1 sym( )   D FF with respect to the initial volume so that

A simple compressible material model
The one the simplest class of compressible materials is defined by which refers to compressible neo-Hooekan model, cf. [3]. The invariant 1 I is define as can be interpreted as an initial shear modulus. Note that the stored-energy function does not depend on the second invariant 2 I . Generally, the potential (7) is not polyconvex one. Nevertheless, it satisfies a rank-one convexity condition if   J  is convex, cf. [6]. From now on, we assume that A general form of the constitutive equation based on (7) in material and spatial description respectively reads as The function  needs to fulfill proper growth conditions, see [5]. We consider the form defined by where 0  is the classical Lamé constant found in linear elasticity, cf. [3]. Evidently, Corresponding constitutive equations are given by (10)

UHYPER
The main goal of the article is to compare a computational efficiency of different implementation of a hyperelastic material model. The ABAQUS /Standard FEM software offers basically two ways of to do so. The first approach is to employ user subroutine UHYPER, which can be used to define an isotropic model. It requires that the storedenergy function and its derivatives be defined with respect to the invariants of the form the Zaremba-Jaumann objective stress rate, see [4,7,8], can be shown to be The derivatives of the stored-energy function defined by (7) and (9) are presented in Table 1. Note that third other derivatives do not result from (13). Thus, in the case of pure displacement formulation these can be omitted in a UHYPER user subroutine. Table 1. Derivatives of the stored-energy function.

First order derivatives
Second order derivatives Third order derivatives

UMAT
User subroutine UMAT allows a general way to define a material's constitutive model, see [2]. Basically, the routine requires a definition of Jacobian matrix (DDSSDE) necessary in the FEM procedure and an array STRESS containing components of the Cauchy stress tensor. Similar to the UHYPER, a user may define material behavior dependent on field variables or states ones, see [2]. In order to implement a hyperelastic material model using UMAT, it is convenient to define a symmetric fourth-order tensor as which is known as the material (Lagrangian) elasticity tensor, cf. [4,8]. It is conjugated to the second Piola-Kirchhoff stress tensor as .  T E C . Since we are interested in the spatial description, the push forward operation is applied to give the Eulerian elasticity tensor     : : Finally, the Zaremba-Jaumann rate of the the Kirchhoff stress is given by The material elasticity tensor corresponding to the compressive neo-Hookean material (7) is of the form Therefore, the elasticity tensor included in (16), cf. [7], can be expressed as The UMAT procedure requires the explicit definition of a material's constitutive relation. Because of the fact that constitutive description of a material in the ABAQUS/Standard is based on the Cauchy stress, the tangent moduli must be transformed as follows, cf. [2]: which is obtained after incorporating an integration procedure in (16). The way leads to a consistent Jacobian, which ensure rapid convergence in the Newton-Raphson procedure. Thus, on the basis of the forward Euler method, the Jacobian is given as stated in (20).
The ABAQUS/Standard solver stores the stress tensor in a one-dimensional array of the size equals to the number of non-zero components NTENS (6 in 3-dimensional case, 4 in plane cases), and the Jacobian is stored as a two-dimensional array of the size NTENS x NTENS . Shear strains components are defined as the engineering ones  It can be noted that, despite the general similarity of the presented measures, the solutions of tasks using UHYPER subroutine take 5% less CPU time and require around 2-5% fewer iterations than UMAT one. All of the 2D tasks yield exactly the same results for both subroutines. However, the convergence of the jobs for model C differs-the UHYPER model is generally more stable than UMAT one (although that depends on the time increment assumedthe larger the time increment, the larger the difference between methods). All tasks, static compression and tension of the tube and the compression solved via Riks method, are aborted in earlier phase for UMAT model. The reaction and deformations of the model are shown in fig. 5  and fig. 6. It is important to note that the load capacity of the tube is the same for both Riks and the static standard analyses obtained on the basis of UMAT procedure.

Conclusions
Two different approaches to the definition of hyperelastic user material in the ABAQUS software are analyzed. UHYPER and UMAT subroutines implementing the compressible neo-Hookean material are created and multiple FE jobs are compared.
Although the jobs yield similar qualitative results for two-dimensional problems, the UHYPER subroutine shows better convergence, as well as requires less CPU time. The three-dimensional problem, on the other hand, shows great divergence between the obtained results. Clearly, these are different solutions of nonlinear boundary value problem as for such class there is no unique solution. UMAT model is generally less stable, aborting the jobs due to too many cutbacks in incrementation in a row.
Authors would like to thank Professor Stanisław Jemioło for his encouragement and useful suggestions on the subject.