Numerical implementation of the Murnaghan material model in ABAQUS / Standard

The paper presents a numerical implementation of the Murnaghan material model (M) [1] in the finite element method software ABAQUS / Standard v. 6.14 [2]. The UHYPER user subroutine is employed, which is suitable for the class of isotropic hyperelastic models [3]. As a special case of the M model, the Saint Venant-Kirchhoff (SVK) model is considered [4]. Formal verification on the basis of elementary tests is performed. Among others, a special attention is paid to a simple shear deformation. In all tested types of deformation, analytical values confirms results based on the finite element procedure within assumed numerical precision and accuracy. It should be noted that the stored-energy function of the M and SVK models do not meet any requirements of the mathematical theory of non-linear elasticity [4, 5]. Therefore, these models are suitable for relatively small deformations, while there are no restrictions on finite rotations. As an example of applications, a tube under axial compression is considered in two cases. Various starting parameters for the Riks procedure [6, 7] are adopted to obtain different solutions of corresponding boundary value problem. Material parameters of steel are considered according to Lurie [8].


Introduction
The Murnaghan (M) and Saint Venanta-Kirchhoff (SVK) hyperelastic material models are widely discussed in the literature, see [1,4,5,8].These models are of great theoretical and practical importance, e.g. in the stability analysis or studies of waves propagation [9].The aim of the article is to implement the M model and its special case, i.e. the SVK model in the ABAQUS / Standard software.To this end, user subroutine UHYPER is employed.
In order to verify the implementation of the models, elementary tests based on homogeneous deformations are performed.In the paper, results only for a simple shear test are presented.As an example of applications of discussed models, a cylindrical tube under compression is analyzed.The boundary value problem is considered in two variants, i.e. with and without a circular cutout in the central part of the tube.The imperfection is incorporated to invoke locally significant inhomogeneous deformation of the body.

Stored-energy potential of the Murnaghan model
A stored-energy function of the Murnaghan model in terms of the invariants of the Lagrangian strain tensor E is given by, cf. [1, 5, 10] where 0  , 0  are Lamé constants (Murnaghan's first order constants) and stand for Murnaghan's second order constants.The function (1) defines elastic potential of SVK material model if the second order constants are omitted.Table 1.Material parameters of the M model, cf.[8].A constitutive equation, which links the second Piola-Kirchhoff tensor T and Lagrangian strain tensor E , is defined by the partial derivative of the stored-energy function with respect to E , see [5, 10]: Relationship between the first Piola-Kirchhoff stress tensor S and the deformation gradient F is given by  S FT .It is important due to the fact that the stress tensor S appears in the differential equilibrium equation with respect to the initial configuration of a body.On the other hand, the spatial description of equilibrium requires the Cauchy stress tensor.It is related to the first Piola-Kirchhoff tensor by / T J  σ SF .

Implementation of the M model in ABAQUS/Standard
In order to implement the M model in the ABAQUS/Standard software, UHYPER user subroutine is employed, see [3].Firstly, the energy potential needs to be expressed as a function of invariants with the volume change eliminated.It is convenient to choose decomposition of the deformation gradient as follows: ).The UHYPER procedure requires to implement a stored-energy function per unit volume in the form and its corresponding derivatives with respect to the invariants, see [3].Then, hypoelastic objective stress rate with the symmetric fourth-order tangent operator C is obtained as where "  " means the Zaremba-Jaumann derivative , and " ."denotes dot product between a fourth-order and a second-order symmetric tensors, cf.[11][12][13].
In the case of the Murnghan material, the energy potential in a proper form for UHYPER is given as , , 1 .
Because of relations only five material parameters in (5) are independent.These are expressed in terms of the Lamé and Murnaghan's constants as follows: First and second order derivatives of ( 5) with respect to the invariants 1 2 , , I I J have to be defined in the procedure.Moreover, in the case of the hybrid formulation third order derivatives with respect to J are required.A subroutine written in FORTRAN/C++ should be link to an input file as described in the ABAQUS/Standard documentation, see [3].

Formal verification of the implementation
The implementation of the M model is verified in simple tests with homogeneous deformations using 3D elements and 2D (plane strain) as well, e.g.C3D8 and CPE4.The results are compared with the values obtained on the basis of analytical formulas.In addition, the procedure is tested in the case of the Saint Venanta-Kirchhoff (SVK) model.Tests are performed on the basis of material parameters of steel according to Table 1.
Basic numerical tasks concerned uniaxial and biaxial deformations.For such homogeneous states, stress tensor components are easily derived from the constitutive equation.Therefore, the accuracy of numerical results may be straightforward verified.
In all tested types of deformations, analytical values confirms results based on the finite element procedure within assumed precision of calculations.Nevertheless, a boundary value problem based on the constitutive equation of the M model is in general ill-posed.It arises from the fact that the stored-energy function is not a convex one with respect to principal deformations i  .What is more, the potential does not fulfill the requirement of proper growth, see [5].Hence, volumetric changes are not described properly.In the case of significant range of i  , the function fails to meet 0 W  for chosen parameters of steel, cf.[8].The situation is physically unacceptable.  (dashed).
One of the most important homogeneous deformation is a simple shear.It is defined as an isochoric plane deformation, i.e. 1 J  .Corresponding deformation gradient F is represented by with  being called the amount of shear.
In contrast to the uniaxial or biaxial deformation, the state associated with simple shear is also characterized by the rotation.The rotation tensor R may be found using polar decomposition of   F RU VR .The right U and the left V Cauchy-Green tensors describe stretches in material and spatial sense, respectively.Unlike the linear theory of elasticity, normal stresses appears in the case of finite strains (the Poynting effect).It is especially relevant type of deformation for nearly incompressible materials since shearing is the major behavior of these in applications, e.g.rubber.

Fig. 2. Components of the Cauchy stress tensor -simple shear test -parameters of steel
Plots of ij  and corresponding results obtained on the basis of the UHYPER subroutine are presented in Fig. 2 5 Examples

A tube under axial compression
The first numerical test is an axial compression of a circular tube with characteristic dimensions: height 20 m h R  , internal radius 0.95 at the other one according to the coordinate system as in Fig. 3.

Fig. 3. Plots of averaged stress 33
 -displacement control -M and SVK model.https://doi.org/10.1051/matecconf/201819601042XXVII R-S-P Seminar 2018, Theoretical Foundation of Civil Engineering Fig. 3 presents plots of averaged stress 33  at the loaded end of the tube as a function of displacement of the end of the tube.Characteristic deformations of the body corresponding to marked points on equilibrium paths are shown in the figure.The main difference between these two results is a range of the displacement load obtained during the process.In the case of the SVK model, computations are carried out in the full scope, i.e. , while for the M model, the extreme displacement value is 3 2 .
It results directly from the significantly different form of the stored-energy function for material parameters of the models.In both cases, the local form of buckling is observed, followed by a very rapid reduction of the stiffness of the tube.

A tube with an imperfection
The second numerical test concerns an axial compression of a tube with a circular cutout, which is located in the central part of the element.Radius of the hole equals 0.2 m R .The introduction of the imperfection is intended to induce relatively large deformations in contrast to the previous example.This type of problem is important from the point of view of applications of cylindrical shells in engineering, see [14] and the references given there.Boundary conditions is specified in the same manner as in the previous task.In addition, calculations are carried out by adopting two different variants of starting parameters in the Riks algorithm, see [2].The results for these variants are marked on the graphs M-1, SVK-1 and M-2, SVK-2 respectively, see Fig. 4-6.both the M and SVK models.In both cases, the buckling form is a global one.Again, results in the case of the SVK model are obtained in the significantly larger displacement range compare to the M model.

Remarks and conclusions
Regardless of the selection of starting parameters, equilibrium paths are convergent in the post-critical behavior.This is especially evident in the case of the SVK model.It should be noted here that the various values of these parameters significantly affect the equilibrium path in the entire scope of the solution.
An influence of a finite element type or a mesh quality on the solution is not investigated in the paper.In the case of the tube with the cutout, where significant local deformations are obtained, a minimization of errors related to the mesh would probably lead to performance a greater range of the displacement load in the M model.For the problem based on the SVK model, clearly better conditioned due to the form of the energy potential, a significant improvement of the obtained solution is doubtful.
The presented solutions of the boundary value problems indicate that the application of the M model in practice is justified for relatively small deformations.Standard metallic materials, including steel, are characterized by significantly lower values of the yield stress than these obtained even for small deformations in the examples.Thus, solution and results should be analyzed in the base of a hyperelasto-plastic model, see [11].A scope of applicability of the M model in the context of constraints on material parameters is discussed in [5,10].
model is built of 6002 finite elements (3 elements by thickness) of type C3D8 with linear interpolation functions.In order to compare results, the problem is solved for the Murnaghan model (M) as well as for Saint-Venant Kirchhoff one (SVK).The loading is realized by displacement control.Thus, following boundary conditions are applied:

Fig. 4 .Fig. 5 .
Fig. 4. Plots of averaged stress 33  -displacement control -M model.The relevant plots are divided into fragments preceding the first critical point (black color) and the further range of the equilibrium path (blue color).Given points on the paths indicate sequence of obtained results in the process.Indeed, the non-smooth equilibrium path for both models confirms the significant impact of local deformations induced by the imperfection on the behavior of the analyzed tube.In the case of the M model, the buckling form is global, while in the SVK model local instability is observed associated with the cutout.Similarly to the previous example, the

Fig. 6
Fig. 6 presents a comparison of plots of averaged stress 33  at the loaded end of the tube as a function of displacement of the end of the tube in the case of second variant of starting parameters in the Riks procedure.Characteristic deformations of the body corresponding to marked points on equilibrium paths are shown in the figure.In compare to the first variant, see Figs. 4 and 5, a significantly smoother equilibrium path is obtained for : MATEC Web of Conferences 196, 01042 (2018) https://doi.org/10.1051/matecconf/201819601042XXVII R-S-P Seminar 2018, Theoretical Foundation of Civil Engineering