Consistent numerical implementation of hypoelastic constitutive models

In hypoelastic constitutive models, an objective stress rate is related to the rate of deformation through an elasticity tensor. The Truesdell, Jaumann, and Green–Naghdi rates of the Cauchy and Kirchhoff stress tensors are examples of the objective stress rates. The finite element analysis software ABAQUS uses a co-rotational frame which is based on the Jaumann rate for solid elements and on the Green–Naghdi rate for shell and membrane elements. The user subroutine UMAT is the platform to implement a general constitutive model into ABAQUS, but, in order to update the Jacobian matrix DDSDDE in UMAT, the model must be expressed in terms of the Jaumann rate of the Kirchhoff stress tensor. This study aims to formulate and implement various hypoelastic constitutive models into the ABAQUS UMAT subroutine. The developed UMAT subroutine codes are validated using available solutions, and the consequence of using wrong Jacobian matrices is elucidated. The UMAT subroutine codes are provided in the “Electronic Supplementary Material” repository for the user’s consideration.


Introduction
Hypoelasticity is a rate form of elastic material model [1], in which an objective stress rate is linearly related to the rate of deformation by means of a fourth-order elasticity tensor which, in general, is not obtainable from a strain energy density. Originally, Dienes [2] showed that the zero-graded hypoelastic model, i.e. a hypoelastic model with constant isotropic elasticity tensor, exhibits oscillation in simple shear, if it is constructed based on the Jaumann rate of the Cauchy stress. However, zero-graded hypoelastic models which are based on the Truesdell or Green-Naghdi rates do not suffer this problem [3]. We remark that the definition of elements in ABAQUS is based on the Jaumann rate for solid elements [4] and on the Green-Naghdi rate for structural elements (shells, membranes, beams, trusses) [5,6], as mentioned in the ABAQUS Theory Manual, Section 1.5.3 [7].
This study aims to formulate and implement various hypoelastic constitutive models into the ABAQUS UMAT (user material) subroutine. To attain this, it is essential to express the elasticity tensor of the hypoelastic model in terms of the elasticity tensor which relates the Jaumann rate of the Kirchhoff stress tensor to the rate of deformation. According to Pinsky et al. [8], such relations seemed difficult to be constructed for models associated with the Green-Naghdi stress rates; however, the kinematical relations provided in Mehrabadi and Nemat-Nasser [9] enable us to establish such connections.
The study starts with a review of some basic definitions, the concept of objective rate and the structure of hypoelastic constitutive models. Next, the relations between the elasticity tensors of various hypoelastic where the material point X denotes the position of a particle in the reference configuration B and the spatial point x = φ(X, t) is the current placement of point X at time t. The placement or configuration of the body B at time t is denoted The codomain restriction to B t ≡ φ(B, t) of the configuration map φ( ·, t) is required to be invertible, continuous, and differentiable along with its inverse, i.e. a diffeomorphism. The deformation gradient F is a two-point tensor [12], mapping material vectors M attached at point X into spatial vectors m = F M attached at the spatial point x and is defined, in components, as The determinant J = det F of the deformation gradient F describes the local change of volume and, due to the requirement of invertibility and regularity of φ, it must be strictly positive. Cauchy's polar decomposition theorem allows to express the deformation gradient F as where U and V are symmetric and positive definite tensors and R is a proper orthogonal tensor. Moreover, the right stretch tensor U is completely material, the left stretch tensor V is completely spatial and the ZAMP Consistent numerical implementation of hypoelastic Page 3 of 23 156 rotation tensor R is, like F , a two-point tensor. The right and left Cauchy-Green deformation tensors are completely material and completely spatial, respectively, and are defined as The Eulerian spatial velocity v is defined as the vector field such that v(x, t) = v(φ(X, t), t) =φ(X, t), (6) whereφ(X, t) ≡ (∂ t φ)(X, t) is the Lagrangian spatial velocity at x = φ(X, t). The gradient of the spatial velocity field v with respect to the spatial coordinates is called the velocity gradient and has components ij = v i,j . The symmetric and skew-symmetric parts of the velocity gradient are the rate of deformation d and the spin or vorticity tensor w, respectively: The skew-symmetric tensor where R is the two-point rotation tensor of the polar decomposition (4) of F , is called rigid spin [6]. Its skew-symmetry can be easily shown by taking the time derivative of R R T = i, where i is the spatial identity tensor. It can also be shown that the spin tensor w and the rigid spin tensor Ω are related by (see, e.g. [13] or [14]) There are two cases for which w and Ω coincide. The first case is when the motion is rigid and F = R, U = I,U = 0 and = w (or, equivalently, d = 0). The second case is when the (normalised) eigenvectors of U remain constant during the motion, which implies thatU has the same eigenvectors and thus commutes with U . The fundamental measure of stress in continuum mechanics is the Cauchy stress σ, which linearly relates the unit normal vector n at a point on the boundary ∂R t of an arbitrary region R t ≡ φ(R, t), subset of the current configuration B t ≡ φ(B, t), to the corresponding surface traction vector t n , i.e. t n = σ n. (11) The Cauchy stress σ is power-conjugated to the rate of deformation d (or, equivalently, to the velocity gradient ), in the sense that the internal power (or deformation power) in an arbitrary region R t is given by Another measure of stress, often employed in numerical applications, is the Kirchhoff stress τ , which is obtained by pulling the integral (12) on R t back to the referential region R, subset of B by means of the theorem of the change of variables, i.e.
Palizi, Federico and Adeeb ZAMP from which we obtain the relation

Objective stress rates
The principle of material frame indifference states that the constitutive equations must be form-invariant under changes of the frame of reference, i.e. under arbitrary rototranslations [1,16]. The regular substantial time derivative of a spatial measure of stress, such as the Cauchy stress, is not frame-indifferent, as we shall briefly show now. Let us denote the substantial time derivative by a superposed dot, i.e.
Under a spatial rotation Q, the stress transforms as Therefore, the rateσ transforms as (see, e.g. [13,14]) which does not preserve the form of the transformation (16). In contrast, it can be shown [13,14] that the rate of deformation d is frame-indifferent. Objective stress rates have been introduced precisely to overcome the problem suffered by stress rates and are all essentially based on the use of Lie derivatives, as elegantly shown by Marsden and Hughes [12]: the Oldroyd rate is precisely a Lie derivative, the Jaumann rate is a linear combination of Lie derivatives, the Green-Naghdi rate is modelled after a linear combination of Lie derivatives, and finally, the Truesdell rate (of the Cauchy stress) is a Lie derivative involving the volume form (or volume element). In this work, we use the objective rates below, in the notation of Bonet and Wood [13]. Truesdell Rate of the Cauchy Stress: Truesdell Rate of the Kirchhoff Stress (also called Oldroyd rate [17] and coincident with the Lie derivative of the Kirchhoff stress [13]): Jaumann Rate of the Cauchy Stress: Jaumann Rate of the Kirchhoff Stress: Green-Naghdi Rate of the Cauchy Stress: Green-Naghdi Rate of the Kirchhoff Stress: The proof of objectivity and non-objectivity of various kinematic and stress variables can be found in Sections 4.15 and 5.6 of Bonet and Wood [13] or Bonet et al. [14]. A more extensive discussion is provided in Chapter 1, Box 6.1 of the book by Marsden and Hughes [12], who also make the distinction between objectivity with respect to isometries (i.e. rototranslations), which coincides with frame indifference, and objectivity with respect to diffeomorphisms, which coincides with the condition of covariance.

Co-rotational frames
As discussed in Sect. 2.2, the definition of objective rates such as Truesdell, Jaumann, and Green-Naghdi is based on the notion of Lie derivative. Here, we show that, in an appropriate co-rotational frame, the components of the Jaumann and Green-Naghdi stress rates can be expressed as the time rate of the components of the stress (see also Section 1.5.3 of ABAQUS Theory Manual [7]). A time-dependent basis {e α } 3 α=1 is co-rotational with respect to the spin tensor w if it transforms following Poisson's Theorem 3 locally, i.e. according to the value of the spin tensor w at the point considered:ė α = w e α = e α w T = −e α w. (24) Considering that w = w μν e μ ⊗ e ν and using the properties of the tensor product, Eq. (24) readṡ In this w-co-rotational frame, the Cauchy stress σ (exactly the same considerations can be made for the Kirchhoff stress τ ) reads and thus its substantial time derivative iṡ where we used Eq. (25). Switching indices μ and α in the second term and ν and β in the third term, we have and, solving forσ αβ e α ⊗ e β , we obtaiṅ σ αβ e α ⊗ e β = (σ αβ e α ⊗ e β )˙+ σ αν w νβ e α ⊗ e β − w αμ σ μβ e α ⊗ e β .
On the right-hand side of Eq. (29), we recognise the Jaumann rate of the Cauchy stress, i.e.
Therefore, comparing Eqs. (29) and (30), we finally obtain that, in this w-co-rotational basis, the components of the Jaumann rate are equal to the (substantial) time derivatives of the components of the stress, i.e.
Analogously, if we define a basis {e α } 3 α=1 that is co-rotational with respect to the rigid spin tensor Ω, in the sense thatė we can show that, in this Ω-co-rotational basis, the components of the Green-Naghdi rate equal the (substantial) time derivatives of the components of the stress, i.e.

Linearisation of the deformation
Linearisation is essential in finite element formulations [13]. The linearisation of the deformation about a specific configuration map φ entails the evaluation of the deformation gradient F and of all derived quantities, after a perturbation is applied to φ. We call this perturbation an infinitesimal displacement, 4 which we denote δu when seen as a function of the material point X and δu when seen as a function of the spatial point x, i.e.
The perturbed configuration map is thus The deformation gradient of the perturbed configuration mapφ is given by (see Section 4.2 in the book by Marsden and Hughes [12])F where F is the deformation gradient of the unperturbed φ and we recall that the large nabla, ∇, denotes the gradient performed with respect to the referential coordinates X I . It is helpful to express (36) as a multiplicative decomposition [19]. Indeed, by definition of inverse, we havȇ Then, the transformation rule for the gradient states that where we recall that the small nabla, ∇, denotes the gradient with respect to the spatial coordinates x i . Therefore, we can write (37) asF where is the deformation gradient mapping from the configuration B t = φ(B, t) to the perturbed configuration B t =φ(B, t). Finally, the multiplicative decomposition equivalent to (36) is given by (39) and (40) as In the following, we shall need to perform the linearisation in time, i.e. considering the relative infinitesimal displacement between the configuration B t = φ(B, t) at time t and the configuration B s = φ(B, s) at time s > t. In this case, the perturbed configuration mapφ and the perturbed deformation gradient F are replaced by the configuration map φ( ·, s) and deformation gradient F ( ·, s) at time s > t. Thus, we have and where the subscript s emphasises that the relative displacement (δu) s (X, t) = φ(X, s) − φ(X, t) points to φ(X, s). Passages analogous to those seen above yield where the relative deformation gradient (δF ) s (x, t) is given by We note that this treatment is based on the concept of relative deformation, for which we refer the Reader to the treatise by Eringen [20].

Hypoelasticity
Hypoelasticity describes a class of elastic materials defined in rate-form. Truesdell and Noll [1] proved that the general forms of hypoelastic constitutive equation must necessarily be given by where s is an objective rate of a stress tensor s, which could be either the Cauchy stress σ or the Kirchhoff stress τ , d is the rate of deformation and C is the fourth-order elasticity tensor, which, in general, depends on the stress s, i.e. C =Ĉ(s). As mentioned in Sect. 1, in this study, the hypoelastic constitutive equations are assumed to be zerograded, which means that the fourth-order elasticity tensor C is isotropic and independent of the stress, i.e. can be expressed as where the superscript zero stands for zero-graded material, λ and μ are Lamé's constants, i is the spatial identity tensor with components δ ij (the Kronecker delta), and the special tensor product ⊗ is defined by Curnier et al. [21] and is such that the component representation of (47) is Given the Young's modulus E and the Poisson's ratio ν, the Lamé's constants can be calculated via the well-known relations

Hypoelasticity in ABAQUS user subroutine UMAT
A clear understanding of variables provided by the UMAT subroutine is essential to properly implement constitutive models into ABAQUS. In this section, we shall discuss the essential theoretical and numerical aspects of the UMAT-subroutine variables. It should be noted that, also for the case of a hyperelastic material, which is described by a strain energy density, the Jacobian matrix of the UMAT subroutine DDSDDE should be updated based on a particular rate-form equation. In this sense, the treatment below holds for both hyperelastic and hypoelastic models.

UMAT-subroutine variable DDSDDE
Based on Section 1.1.44 of ABAQUS User Subroutines Reference Guide [7], in order to implement a hyperelastic material, for which the Cauchy stress tensor is explicitly expressed in terms of the deformation gradient and other kinematic variables, it is necessary to evaluate the consistent Jacobian matrix DDSDDE. This is defined as the matrix representing the fourth-order tensor D featuring in the expression 5 Equation (50) can be promptly justified by using the definition (14) of Kirchhoff stress τ to eliminate the Cauchy stress σ and the definition (21) of Jaumann rate τ of the Kirchhoff stress, i.e.
If we define C τ as the fourth-order tensor relating τ to d, i.e.
comparison of Eqs. (51) and (52) yields A similar formulation is employed for the Jacobian matrix in hypoelastic materials, i.e. for rate-form constitutive models. As we shall show later in Sect. 4, the employment of the consistent Jacobian matrix DDSDDE based on Eq. (53) results in convergence, while other Jacobian matrices may give rise to slow convergence [6,22] or even divergence.
In order to implement a hypoelastic model into UMAT subroutine, it is essential to relate the elasticity tensor of the hypoelastic model to the fourth-order tensor D of Eq. (53), represented by the consistent Jacobian matrix DDSDDE. Below, we report the tensor D that must be used for each choice of hypoelastic formulation (proofs provided in "Appendix A"). Truesdell Rate of the Cauchy Stress: Truesdell Rate of the Kirchhoff Stress: Jaumann Rate of the Cauchy Stress: Jaumann Rate of the Kirchhoff Stress, i.e. the point of departure, Eq. (53): Green-Naghdi Rate of the Cauchy Stress: Green-Naghdi Rate of the Kirchhoff Stress: In the expressions above, we used where λ α are the principal stretches (eigenvalues of the left-stretch tensor V , so that λ 2 α are the eigenvalues of the left Cauchy-Green deformation b), and are the eigenprojections of the left Cauchy-Green deformation b, with i being the spatial second-order identity tensor, as seen earlier. The proof of Eqs. (54), (55), and (56) can also be found in Wooseok et al. [23] and the proof of Eqs. (58) and (59) is based on the relation.
the derivation of which can be found in Mehrabadi and Nemat-Nasser [9], and also Zhou and Tamma [24].

UMAT-subroutine variables DSTRAN and DROT
The UMAT-subroutine variables DSTRAN and DROT contain the components of the incremental strain and incremental rotation, respectively, between the configuration B n = φ(B, t n ) at time t n and the The incremental deformation is evaluated following the procedure outlined in Sect. 2.4 and employs a midpoint formula [25] considering the configuration where t n is the time at the beginning of the increment and (δt) n = t n+1 − t n is the increment from t n to t n+1 . Our goal here is to determine the relative displacement gradient, the midpoint displacement gradient, the midpoint velocity, and the midpoint velocity gradient. From the midpoint deformation rate, symmetric part of the midpoint velocity gradient, we shall derive the incremental strain DSTRAN.
The deformation gradient at time t n+1 is obtained as a function of the relative deformation gradient and of the deformation gradient at time t n , following Eq. (44), as with the relative deformation gradient as in Eq. (45), i.e.
where we use the subscript n + 1 in place of t n+1 in (δF ) n+1 (x, t n ) and (δu) n+1 (x, t n ). Now, we approximate the displacement between t n and t n+1/2 as half the displacement between t n and t n+1 : Using again Eq. (45) for the relative deformation gradient between t n to t n+1/2 with the displacement (66), we have which, together with Eq. (65), yields the alternative expression Similarly, we approximate the displacement between t n+1/2 and t n+1 as where x = φ(X, t n+1/2 ) and it is understood that the displacement 1 2 (δu) n+1 (x, t n ) must be paralleltranslated from x = φ(X, t n ) to x = φ(X, t n+1/2 ). The corresponding incremental deformation gradient is where ∇ denotes the gradient operator at x = φ(X, t n+1/2 ), as opposed to the gradient operator ∇ at x = φ(X, t n ). Using the incremental deformation gradients (δF ) n+1 (x, t n ) and (δF ) n+1 (x , t n+1/2 ), we can write the multiplicative decomposition which, comparing with Eq. (64) and using Eqs. (67) and (70), yields The midpoint displacement gradient ∇ (δu) n+1 (x, t n ) is evaluated from Eq. (72), considering Eq. (68), as The velocity in the configuration B n+1/2 is approximated by the midpoint discrete derivative where, again, the parallel translation from x = φ(X, t n ) to x = φ(X, t n+1/2 ) is understood. Based on this approximation, the velocity gradient in the midpoint configuration B n+1/2 is given by with ∇ (δu) n+1 (x, t n ) given by Eq. (73). The corresponding deformation rate and spin are The infinitesimal strain increment δ n is the tensor corresponding to the UMAT-subroutine variable DSTRAN and is obtained from the midpoint deformation rate (76a) by multiplying by the time increment (δt) n : The midpoint spin (76b) is instead used to calculate the incremental rotation tensor Q n corresponding to the UMAT-subroutine variable DROT, according to the Hughes-Winget algorithm (see Hughes and Winget [26] and Section 14.10.6 of Neto et al. [25]), as the usage of which will be explained in the Sect. 3.3.

UMAT-subroutine variable STRESS
The stress array STRESS is the key UMAT-subroutine variable that provides the user with the components of the Cauchy stress tensor at the beginning of an increment. The user is, then, required to update the STRESS array with the components of the Cauchy stress tensor at the end of the increment. In some cases, such as that of hyperelastic constitutive models, the user can easily utilise the deformation gradient at the end of the increment, i.e. UMAT-subroutine variable DFGRD1, to calculate the Cauchy stress components. However, for rate-dependent constitutive equations, updating the STRESS array requires careful consideration on part of the user. ABAQUS uses the Hughes-Winget algorithm to integrate the rate-form constitutive equations and update the STRESS array (see Section 3.2.2 of ABAQUS Theory Manual [7]). Based on the Hughes-Winget algorithm, the stress update reads [26] where σ n and σ n+1 are the Cauchy stress tensor at time t n and t n+1 , respectively, Q n is the incremental rotation tensor and C σ n is the elasticity tensor relating the Jaumann rate σ of the Cauchy stress to the rate of deformation d (see Eq. (A.14)). By defining we can write Eq. (79) in the form The tensorσ n is the quantity corresponding to the UMAT-subroutine variable STRESS (see Section 1.1.44 of the ABAQUS User Subroutines Reference Guide [7]). Based on Eq. (56), relating C σ to the consistent Jacobian D, we have where D n and σ n are the consistent Jacobian and the stress at the beginning of step n, respectively. Using Eq. (82), we can write (81) in terms of D n as The form of the consistent Jacobian D n depends on the selected hypoelastic model, see Eqs. (54)-(59).

Results and discussion
Four simulations are performed in ABAQUS to validate the update procedure of the Cauchy stress array STRESS and the Jacobian matrix DDSDDE for the hypoelastic constitutive equations. Each model consists of only one eight-node brick element C3D8 with unit dimensions L = 1 (Fig. 1). The models are divided into displacement-based models and force-based models. The displacement-based models are used to check the correctness of the updated Cauchy stress tensor. In these models, the iterative procedure is not involved and, thus, the Jacobian matrix does not affect the material response. In contrast, in the forcebased models, the implementation of the correct Jacobian matrix plays a crucial part in the convergence of the analysis. Employing the force-based models serves to validate the Jacobian matrices used in the hypoelastic constitutive equations and to show the consequences of using a Jacobian matrix that does not correspond to that required by the element type (for solid elements, that related to the Jaumann rate). In ABAQUS, the Static-General step employs the Newton-Raphson method for the iteration. For the force-based models, the incrementation type is set to automatic, and, for the displacement-based models, the fixed incrementation type is used. For automatic incrementation, the initial increment size is   selected 100 times smaller than the maximum increment size and the minimum increment size is chosen 100 times smaller than the initial increment size. The hypoelastic constitutive equations used in these four models are based on: (i) the Jaumann rate of the Cauchy stress, (ii) the Jaumann rate of the Kirchhoff stress, (iii) the Truesdell rate of the Cauchy stress, (iv) the Truesdell rate of the Kirchhoff stress, (v) the Green-Naghdi rate of the Cauchy stress, and (vi) the Green-Naghdi rate of the Kirchhoff stress. The constitutive equations are all zero-graded. The Young's modulus E and Poisson's ratio ν are, respectively, set to 20 (with consistent units) and 0.2. We verified analytically that, for the cases studied, the material parameters E and ν appear linearly in the expressions of the Cauchy stress and, therefore, varying their magnitude does not introduce any further  nonlinearity. Therefore, the conclusions drawn in this section are valid for arbitrary (and, naturally, physically admissible) values of the Young's modulus E and the Poisson's ratio ν. In the first model (Fig. 1a), the nodes located on the bottom surface of the element are fixed in all directions, the nodes on the upper surface are fixed in x 2 and x 3 directions and a displacement equal to 5 is applied at the upper nodes in x 1 direction. This model is a displacement-based model subjected to simple shear. Using the zero-graded hypoelastic constitutive equations in this model, the numerical Cauchy stress tensors can be checked with the analytical counterparts provided in [2,3] (see Fig. 2 for the analytical solutions of the zero-graded hypoelastic constitutive equations under simple shear).
In simple shear, the local volume does not change (i.e. J = 1) and, thus, the Kirchhoff stress tensor coincides with the Cauchy stress tensor. Therefore, the hypoelastic constitutive equations based on the Kirchhoff stress tensor are equivalent to those based on the Cauchy stress tensor. Fig. 3 compares the analytical solutions (CASE 4) with the numerical ones, considering the fixed increment size as 0.1 (CASE 1), 0.01 (CASE 2), and 0.001 (CASE 3). For all hypoelastic constitutive equations, the numerical solutions corresponding to the increment size of 0.001 lie on the analytical solutions. For the hypoelastic models associated with the Jaumann and Green-Naghdi rates, the convergence occurs at the increment size of 0.01 (CASE 2). In the second model (Fig. 1b), the nodes on the back surface are fixed in all directions; the nodes on the front surface are fixed in x 2 and x 3 directions and a displacement equal to 5 is applied in x 1 direction. This model is a displacement-based model subjected to uniaxial extension, which we use to double-check the update procedure of the Cauchy stress tensor for the zero-graded hypoelastic constitutive equations. In uniaxial extension, the eigenvectors of the right stretch tensor U remain constant during the motion and thus the eigenvectors ofU coincide with those of U . Consequently, U andU commute and the spin tensor w reduces to the rigid spin Ω, as can be seen from Eq. (10). Accordingly, in uniaxial extension, the hypoelastic constitutive equations based on the Jaumann rates show the same mechanical behaviour as those based on the Green-Naghdi rates. Fig. 4  In the third model (Fig. 1c), the bottom nodes are fixed in all directions; by defining a rigid body constraint, the upper nodes are rigidly constrained to a reference point located at the centre of the upper surface and a concentrated force in the x 1 -direction is applied to the reference point. The reason behind the employment of the rigid constraint is that the application of equal concentrated forces at the upper nodes does not produce simple shear deformation. This is a force-based model, used to check the correctness of the Jacobian matrix for the zero-graded hypoelastic constitutive equations.
In this model, the applied force is in equilibrium with the shear Cauchy stress component, i.e. σ 12 (Fig. 1c). Noting that the iterative procedure, i.e. the Newton-Raphson method, is unable to pass the critical points (e.g. snap-through), for each hypoelastic constitutive equation, the value of the applied force must be selected based on the shear stress response in the displacement-based model under simple shear. Based on Fig. 3b c, for the hypoelastic constitutive equations associated with the Truesdell and Green-Naghdi rates, the shear stress is increasing in the course of the simple shear motion and its maximum value occurs at u 1 = 5. Nevertheless, for the hypoelastic constitutive equations associated with the Jaumann rates, the maximum shear stress occurs at a critical point located at u 1 = 1.57 (Fig. 3a). Employing the zero-graded hypoelastic constitutive equations and the corresponding applied force values into the model, in Fig. 5, a comparison is provided between the analytical solutions (CASE 4) and the numerical solutions using the automatic incrementation with maximum increment sizes as: 0.1 (CASE 1), 0.01 (CASE 2), and 0.001 (CASE 3). For all hypoelastic constitutive equations, the maximum increment size of 0.001 (CASE 3) provides numerical solutions which are located on the analytical solutions. For the hypoelastic constitutive equations based on the Jaumann and Green-Naghdi rates, even the numerical solutions associated with the maximum increment size of 0.1 lie on the analytical counterparts (Fig. 5a,  b). For the fourth model (Fig. 1d), the nodes in the back surface are fixed in all directions, the nodes in the front surface are restricted in x 2 and x 3 directions and, in the remaining degrees of freedom (DOFs), i.e. DOFs in x 1 direction at the nodes in the front, equal concentrated forces f 1 are applied. This model is a force-based model subjected to uniaxial extension, used to double-check the correctness of the Jacobian matrix for the zero-graded hypoelastic constitutive equations.
In this model, the sum of the concentrated forces, i.e. 4f 1 , is in equilibrium with the normal Cauchy stress component σ 11 (Fig. 1d). Thus, for each hypoelastic constitutive equation, the value of f 1 must be selected based on the response of the stress component σ 11 in the displacement-based model subjected to uniaxial extension (Fig. 4). As is clear from Fig. 4b, for the hypoelastic constitutive equations based on the Jaumann or Green-Naghdi rates of the Kirchhoff stress tensor, a critical point (snap-through) exists at u 1 = 1.71, whereas the stress component σ 11 is increasing during the uniaxial extension motion in other constitutive equations (Fig. 4a, c, d). In Fig. 6 The elasticity tensor C σ n describes the relation between the Jaumann rate of the Cauchy stress and the rate of deformation (see Eq. (A.14)) and is related to the consistent Jacobian by C σ n = D n − σ n ⊗ i, via Eq. (56). Since C σ n is used in the Hughes-Winget algorithm [26] to update the stress (Eq. (79)), a possible coding error is to use C σ n in place of the consistent Jacobian D n (Eq. (57)). Figure 7 compares  the numerical solutions based on the correct Jacobian D n (CASE 1) and the incorrect Jacobian C σ n = D n − σ n ⊗ i (CASE 2) for the force-based model subjected to uniaxial extension. Figure 7 shows how using the wrong Jacobian matrices results in non-convergence of the analyses at the early stages of the loading. Nevertheless, before the simulation fails to converge and stops (which is indicated by the red dots in the CASE 2 plots in Fig. 7), the Cauchy stress components are identical in both cases of the correct and incorrect Jacobian matrices. Therefore, the correctness of the Jacobian matrix is crucial for the convergence of the analyses, but has no effect on the correctness of the updated Cauchy stress tensor.

Summary
In this study, we implement various hypoelastic constitutive models into the finite element analysis software ABAQUS through the user subroutine UMAT. For the formulation of the consistent Jacobian, i.e. the matrix DDSDDE, ABAQUS uses the elasticity tensor relating the Jaumann rate of the Kirchhoff stress to the rate of deformation for solid elements and the elasticity tensor relating the Green-Naghdi rate of the Kirchhoff stress to the rate of deformation for shell elements. Therefore, it is essential to relate the elasticity tensor of various hypoelastic constitutive models to the elasticity tensor associated with the consistent Jacobian. In regard to the importance of the consistent Jacobian, it is shown that the usage of wrong Jacobian matrices would give rise to the non-convergence of the analyses in the early stages of loading. Additionally, in order to update the Cauchy stress in the various hypoelastic models presented, the comprehension of ABAQUS co-rotational framework and UMAT-subroutine variables such as STRESS, DSTRAN, and DROT is essential, and this is why they were all described in detail. The correctness of the stress array STRESS and Jacobian matrix DDSDDE in the zero-graded hypoelastic constitutive equations is checked using displacement-based and force-based models subjected to simple shear and uniaxial extension. This work is aimed at providing a step-by-step guide to the implementation of hypoelastic materials in ABAQUS, but the procedures shown can be adapted to the modelling of hyperelastic materials as well.
The constitutive equation of a hypoelastic model based on the Truesdell rate of the Cauchy stress is in which the elasticity tensor C σ • is given. Considering the expression of the time derivative of the volume ratio (a proof can be found, e.g. in Section 4.14 of [14]), i.e. the following expression will be achieved: Now, based on Eq. (53), we can finally write and, therefore: Now, employing the above relation into Eq. (A.8) leads to, Finally, considering Eq. (53), the elasticity tensor C σ of the hypoelastic model, which is known, could be related to the tensorial version of the consistent Jacobian, i.e. Based on Eq. (62), we haveσ from which, rearranging the terms, Now, using Eqs. (20) and (A.14) and index notation, we obtain Employing Eq. (60) and eliminating d kl from both sides, we get