© AFM, EDP Sciences 2018

1 Introduction

Over the last past years, numerous studies have been concentrated on dynamic mechanical properties of materials. Constitutive laws, the relation between internal forces and deformations, have also received much attention. Under large deformations and high deformation rates such as forming, machining processes and impact analyses, elastoplastic constitutive laws play a significant role in predicting the mechanical behavior of materials.

Although many kinds of constitutive laws have been implemented into the commercial non-linear finite element software Abaqus [1], in some conditions they cannot fit well some complex behaviors of materials such as the deformation process involving large strain, high strain rates, temperature softening, etc. Therefore, Abaqus provides the ability for users to implement only the constitutive flow, or the complete constitutive behavior law by themselves via FORTRAN user subroutines: UHARD or VUHARD for the hardening flow law and UMAT or VUMAT for the complete behavior law. While the user hardening implementation is quite simple and easy, the user material implementation needs more expertise. The major difference from VUMAT to UMAT is that VUMAT subroutines do not need the evaluation of the so-called consistent Jacobian matrix, but only the evaluation of the stress tensor σ₁ at the end of the increment. The usual approach used to write a VUMAT subroutine, because the time increment is very small due to the explicit time integration scheme, is to use a straightforward evaluation of the equivalent plastic strain increment without requiring any local iterative algorithm [2,3].

In this paper, a different approach has been chosen and consist in evaluating, in the VUMAT subroutine, the plastic strain increment using the implicit radial return algorithm presented in Section 2. In order to illustrate the proposed approach and to be able to validate it, the Johnson–Cook constitutive flow law [4], presented in Section 3, has been selected to be implemented into the Abaqus/Explicit code through both a VUHARD and a VUMAT subroutines. The complete implementation of the Johnson–Cook elastoplastic constitutive law through a VUMAT subroutine is presented in details in Section 4. Numerical efficiency and precision of the proposed implementation are finally investigated, in Section 5, through some dynamic benchmarks tests and have shown the efficiency and robustness of the proposed integration scheme. Section 5 also present some alternative constitutive laws results for the Taylor impact test.

2 Time integration algorithm of J₂ plasticity

Time integration algorithms based on finite difference methods are widely used in finite element solutions of mechanics. In these algorithms, time is discretized on a finite grid, and the distance between consecutive points on the grid is defined as the time step Δt. With the positions and some of the time derivatives at time t, the integration scheme calculates the same quantities at a later time t + Δt. Through the iteration of the procedure, the time evolution of those quantities can be traced.

2.1 J₂ plasticity model

The elastic stress–strain equation is usually written in the following incremental form: (1) where is an objective stress rate in order to take into account the objectivity in large deformations, : is the double dot product, is the linear isotropic fourth order elastic tensor given by: (2)with the second order identity tensor, the fourth order identity tensor, ⊗ the Dyadic product, G the shear modulus and K the Bulk modulus linked to the Young’s modulus E and Poisson’s ratio v by the following: (3)and is the so-called elastic part of the rate of deformation tensor generally assumed [5] for hypo-elastic material to conform to the following additive decomposition equation: (4)where is the plastic part of the rate of deformation tensor . According to the literature, J₂ plasticity, also called von Mises yield criterion, was proposed by Huber in 1904 and von Mises in 1913 [6]. It is widely applied in mechanical engineering and metal forming. It assumes the existence of a scalar yield function f given by: (5)where is the von Mises equivalent stress, or effective stress, defined from the deviatoric stress as: (6)and σ^y is the so-called current yield stress of the material, depending in our kind of applications on the the equivalent plastic strain , the equivalent plastic strain rate and the temperature T given by: (7) (8)where η is the Taylor–Quinney [7] latent coefficient defining the amount of plastic work converted into heat energy, C_p is the specific heat coefficient and ρ is the density of the material. The so-called yield function f defines the yield surface (when f = 0), the elastic domain (when f < 0) and the no-access domain (when f > 0). Assuming an associative plastic flow rule, the plastic strain rate can be expressed with the following equation: (9)where γ is a scalar representing the flow intensity and is a second order tensor (the unit normal to the flow stress determined exclusively in terms of the trial elastic stress [8]) representing the flow direction given by: (10)

From equation (7) one can easily obtain: (11)

The plasticity model described here above has to be integrated in time with respect to an incremental objective algorithm. Different expressions of the time derivative of the Cauchy stress tensor used in equation (1) can be used, leading to some noticeable differences in the final results. Abaqus/Explicit for example in case of solid elements uses a Jaumann stress rate for all Built-In constitutive models and a Green–Naghdi stress rate in case of VUMAT user subroutines [1]. This can lead to some difficulties when one wants to compare VUMAT results with Built-In models as we will see further.

2.2 Radial return mapping algorithm

A widely used method for the time integration of J₂ plasticity with isotropic hardening is the return mapping algorithm introduced by Wilkins [9] and Maenchen and Sack [10]. This algorithm is now extensively used and detailed in numerous books such as [6,11] or papers like Simo et al. [8] or Ponthot [12]. We assume that, thanks to the finite element formulation, the strain increment between time t₀ (the last equilibrated configuration) and t₁ (the new equilibrated configuration with t₁ = t₀ + Δt) is and the deviatoric stress tensor at time t₀ is . The trial deviatoric part of the stress tensor and the final hydrostatic pressure p₁ are calculated from the strain increment assuming that the whole step is fully elastic in a first time, so that: (12) and (13)

Conforming to equation (5), at time t₁, the yield function f is therefore: (14) with: (15)and the yield stress of the last equilibrated configuration at the beginning of the current increment (t = t₀). In order to know if the predicted trial stress is admissible or not, a test has to be performed on the sign of the yield function f leading to the two options here after:

• if f ≤ 0, the whole step is fully elastic, which means that the predicted stress is admissible and we conclude that ;

• if f > 0, the trial stress is not admissible and a plastic correction has to be performed in order to compute the final value of .

The plastic correction is computed enforcing the (discrete) generalized consistency parameter f(Γ) = 0 at the end of the current increment (so that time t = t₁): (16) with (17)

Combination of equations (16) and (17) leads to the following final form of the consistency parameter f(Γ): (18) where the only unknown value is the scalar Γ, already defined in equation (11). If one assume here a linear hardening of the material during the time-step, as usually done thanks to the explicit integration scheme used in Abaqus/Explicit, the following analytical solution to equation (16) can be written, as proposed for example by Gao et al. [2] or the Abaqus manual [3]: (19)with , the slope of the hardening law at the current point, assumed to be constant during the time increment Δt. This approach is usually adopted in VUMAT implementations because of its simplicity as this later does not require any iteration to solve the proposed problem. Unfortunately, as it will be presented further, this lead to many instabilities because of the approximation proposed by this approach when the non-linear terms of the constitutive flow law becomes important. In our approach, we have therefore chosen to solve equation (18) using an approach similar to the one presented in Zaera et al. [13] the so called root finding methods used to solve equation (11) will be presented in Section 2.3.

Once the correct final value of the plastic corrector Γ has been obtained, the final values of the deviatoric part of the stress tensor and the updated internal variables are computed thanks to equations (17). At the end of the proposed algorithm, the final stress at the end of the increment is computed using the following equation: (20)

2.3 Root-finding of the non-linear equation

Several root-finding methods can be used to find the numerical solution of equation (18). In our kind of applications, after a correct selection of the bounds of the searching interval for the solution, we know that the requested root is bracketed in a given interval. Once we know that the proposed interval contains a root, several classical procedures are available to refine it. Among the proposed methods, some of them are known to have fast convergence rate, precision and/or robustness [14,15]. Unfortunately, the methods that are guaranteed to converge are known to be the slower ones, while those that exhibit a fast convergence rate can also dash rapidly to infinity without warning if measures are not taken to avoid such behavior. Among the large list of available methods (Chord, Bisection, Regula-Falsi, Ridders’, Newton–Raphson), only two of them have finally been selected for their efficiency and robustness and are presented here after: the bisection method, a slow and sure method and the Newton–Raphson method a fast but sometimes failing method.

2.3.1 The bisection method

The first root-finding method presented in this paper is the so-called bisection method. This method is one that cannot fail, but with a slow rate of convergence. It’s an iterative root-finding method where the predicted interval of the root is successively halved until it becomes sufficiently small, which is also called interval halving method. This process is repeated until the interval is small enough, which satisfies: (21) where ϵ_bis is the error tolerance of the bisection method. Assuming the root to be found is initially bracketed in an interval [x₀,  x₁] satisfying f(x₀)f(x₁) < 0, the method converges linearly, which is comparatively slow, but the main advantage of this method is that it only requires the evaluation of the function f(x) itself and not its derivative f^′(x) as in other methods such as the Newton–Raphson. Therefore, this method can become competitive when the computation of the derivative of the function f(x) becomes long to compute; i.e. when the expression of f(x) is complex.

2.3.2 The safe version of Newton–Raphson method

The Newton–Raphson method [14,15] is the best known method of finding roots because of its simplicity and efficiency (very high rate of convergence). The only apparent drawback of this method is that it requires the evaluation of the derivative f^′(x) of the function f(x), so it’s only usable in problems where f^′(x) can be readily computed, or numerically evaluated as we will see further. From Taylor series expansion of f(x) about x, we obtain the following expression: (22)

The Newton–Raphson is an iterative process starting with a first guess x₀ for a root of the function f(x) and the process is repeated until the convergence criterion given by: (23) where is the error tolerance of the Newton–Raphson method, is reached.

In some situations, the Newton–Raphson method appears to have poor global convergence, because the tangent line is not always an acceptable approximation of the function. But more important, when f(x) is a piecewise function defined by two different expressions on the left and right side of a given point x_s, f(x) will not be differentiable at x_s, leading to numerical difficulties and divergence of the Newton–Raphson algorithm when the seeking solution x_i crosses the point x_s. As a result, the so-called safe version of Newton–Raphson, described in Figure 1, resulting from a combination of the Newton–Raphson and the bisection methods is proposed here after.

If there is a root in the interval [x₀,  x₁] satisfying f(x₀) ⁡ f(x₁) < 0, the safe version of Newton–Raphson method regards the midpoint of [x₀,  x₁] as the first guess of the root and the Newton–Raphson iteration starts. After each iteration, the interval is updated by replacing one of the two boundaries with the new solution. If the iteration is out of the interval, it is disregarded and replaced with a bisection step to re-position the initial guess. The Newton–Raphson algorithm requires to calculate the derivative of the yield function f(Γ) with respect to the Γ parameter when solving equation (18), so that: (24) with the following definition for the derivative of the yield function: (25)

A common method to calculate the derivatives of the yield function σ^y with respect to , and T is to use an analytical method, which is to calculate the analytical expression for every partial derivative based on the hardening flow law of the material. However, for most yield functions, it may be difficult to obtain derivatives using the analytical method. Therefore, in such cases, we propose here after a numerical solution as an alternative. This numerical solution is obtained by adding a small increment to the equivalent plastic strain, plastic strain rate and temperature respectively in order to calculate the three partial derivatives , and in equation (25): (26) (27) (28)

Of course, accurate results depends on a correct choice for the three increments , and ΔT. In all following numerical tests, the three increments have been arbitrary fixed to the same value Δx = 10⁻⁶.

3 Hardening flow law

3.1 The Johnson–Cook hardening flow law

In this paper, the Johnson–Cook hardening flow law [4,16] has been selected to be implemented in Abaqus/Explicit through VUMAT and VUHARD subroutines. Although the proposed approach can be applied to implement other elastoplastic constitutive laws, the Johnson–Cook law is the best choice for validating the proposed approach, because it has already been implemented in Abaqus/Explicit code, which means the benchmark tests can be conducted between the Abaqus Built-In constitutive law and our FORTRAN implementations.

The Johnson–Cook hardening flow law is probably the most widely used flow law for the simulation of high strain rate deformation processes taking into account plastic strain, plastic strain rate and temperature effects. Since a lot of efforts have been made in the past to identify the constitutive flow law parameters for many materials, it is implemented in numerous finite element codes such as Abaqus [1]. The general formulation is given by the following equation: (29) where is the reference strain rate, T₀ and T_m are the reference temperature and the melting temperature of the material respectively and A, B, C, n and m are the five constitutive flow law parameters. The multiplicative formulation of this flow law allows for the following form: (30)where, according to [4,13,17], the dependence on the plastic strain rate is only taken into account if : (31)and the dependence on temperature is defined so that, if T < T₀ there is no temperature dependence of the yield stress and if T ≥ T_m the material is assumed to behave like liquid: (32)

Those two conditions defined by equations (31) and (32) lead to some discontinuities in the hardening relation , its derivative and the yield function itself resulting in numerical difficulties in the iterative solving procedure. The yield function is therefore not differentiable at and T₀. Nevertheless, the robust Newton–Raphson procedure proposed in Section 2.3.2 has been found sufficiently robust to overcome those problems. Analytical forms for the derivatives of the Johnson–Cook flow law σ^y with respect to , and T are given by the three following equations: (33) (34) (35)

The last problem usually encountered is that, because of the term in equation (33), h tends to infinity when the first plastic increment occurs (i.e. when ). Therefore, we will have to apply a special treatment in the hardening parameter computation on the first plastic step. This will be presented further in Section 4.3.

3.2 VUHARD implementation in Abaqus/Explicit

The VUHARD subroutine is a straightforward approach to implement a new constitutive flow law in Abaqus/Explicit by just implementing a FORTRAN subroutine to compute the yield stress of the material and its derivatives with respect to , and T. The main part of the Built-In constitutive law is used for time integration of the stress, for a given time increment, and the provided user subroutine is used to compute the hardening flow law. Very few details are given about this implementation in the Abaqus documentation, but some useful informations are available in Jansen van Rensburg et al. [18].

The numerical implementation of the VUHARD subroutine for the Johnson–Cook flow has been done through a FORTRAN program defining the hardening flow law according to equation (29) and the three analytical derivatives of σ^y with respect to , and T defined by equations (33)–(35). This part of the code is exactly the same as the one that will be used in the VUMAT subroutine presented in the next section of this paper.

4 VUMAT implementation in Abaqus/Explicit

All details concerning the implementation of the radial return mapping algorithm using a Fortran VUMAT subroutine for the Johnson–Cook flow law are detailed in this section. The detailed flowchart algorithm is illustrated in Figure 2. The first block of the proposed algorithm “Start of VUMAT” is used to get the material properties defined as user material constants. Abaqus/Explicit provides the user subroutine VUMAT the quantities defined here after:

• the strain increment for the current time step;

• the stress tensor σ₀ and the temperature T₀ at the beginning of the current increment;

• the time increment Δt corresponding to the current time step;

• a table of solution dependent state variables (SDVs) used to store important data such as , , Γ and transfer them from one increment to the other.

The VUMAT subroutine must compute and return the value of the stress σ₁ and the SDVs variables at the end of the increment for each integration point. The internal and dissipated energies have also to be evaluated in order to compute the temperatures in the model.

4.1 Abaqus package subroutine

The package part in the Abaqus software is a mandatory step used to compute the starting values (compute the value of the time increment Δt for example by taking into account the material behaviour) for a reference time t = 0. For this we need to compute the elastic stress due to a virtual strain increment provided by Abaqus using the following expression: (36)

The computed stress σ₁ will be discarded after the end of the package step.

4.2 Elastic predictor

The aim of this subsection is to calculate the deviatoric part of the predicted elastic stress and the von Mises equivalent stress and compare it with the yield stress at the beginning of the increment in order to test if the current step is fully elastic or partly plastic. Therefore, we need firstly to decompose the initial stress σ₀ into its hydrostatic pressure (p₀) and its deviatoric () parts: (37)

Then, we compute the new pressure (p₁) and the deviatoric part of the trial stress tensor () using: (38) and compare the yield stress at the beginning of the increment to the von Mises equivalent trial stress (see Eq. (15) for this later). Depending on this test, we will then correct or not the final stress.

• if , the plastic corrector is zero, so the plastic correction steps can be skipped. We assume that the predicted stress is the final one , the plastic corrector Γ = 0, the yield stress remains unchanged and we can go directly to the final computations in Section 4.4;

• if , the step is at least partly plastic and the plastic corrector described in Section 4.3 has to be computed in order to draw back the predicted stress onto the yield surface of the material at the end of the current increment.

4.3 Plastic corrector

The safe version of the Newton–Raphson algorithm is used to recover the stress in accordance with the elastoplastic constitutive behavior law. This Newton–Raphson method is used to compute the Γ parameter defining the correction due to the increase of the strain. In order to enhance the computations efficiency, we initialize the value of the Γ parameter to its value at the end of the previous increment (Γ = Γ₀). If the current increment is the first plastic increment (i.e. if ), as we cannot compute the value of h(Γ) when the plastic strain is zero, we initialize Γ to an initial value different of zero (Γ = 10⁻⁸). The interval for the bisection part is initialized so that . The predicted equivalent plastic strain , plastic strain rate and temperature T₁ at the end of the increment are computed thanks to the system of equation (17). The yield stress , the yield function f(Γ) and its derivative f^′(Γ) are then computed. Then, we test the convergence of the Newton–Raphson algorithm by computing the increment ΔΓ of the Γ parameter: (39) and comparing it to the Newton–Raphson precision ϵ_NR defined by the user:

• if ΔΓ > ϵ_NR we need to iterate to compute the correction of the Γ value using the following relation:

(40) and re-evaluate the values of f(Γ) and f ´ (Γ);

• if ΔΓ ≤ ϵ_NR we have obtained the final value of the Γ parameter.

The final deviatoric part of the stress tensor is computed from the predicted value and the correction term Γ using: (41) with , within the framework of the radial return algorithm, defined by: (42)

4.4 Final computations

The main work of the final computations is to update the state variables (SDVs) and the energies. The final stress tensor at the end of the increment σ₁ is computed from the hydrostatic pressure p₁ and the deviatoric part of the stress tensor thanks to equation (20). The equivalent plastic strain , equivalent plastic strain rate and final temperature T₁ at the end of the increment are stored in their respective SDVs variables for subsequent re-use in the next increment. We also have to compute the new specific internal energy e₁ from: (43) and the dissipated inelastic energy from: (44)

At this point, the VUMAT subroutine comes to an end. At this point of the flowchart, the final temperature is not computed, since the software uses a subsequent thermal step to evaluate the temperature raise due to the dissipated inelastic energy and the conduction part. The internal energy is modified during this thermal step and this seems to be taken into account during this extra step (out of the VUMAT subroutine) since the tests on 1 element with imposed displacement of all nodes gives the exact same results as the Build-in routine.

5 Validation of the proposed implementations

In this section, the performance of the Johnson–Cook law programmed in the VUMAT subroutine is compared with the performance of the Abaqus/Explicit Built-In Johnson–Cook law and of the VUHARD Johnson–Cook implementation, in order to validate the proposed implementation approach. The benchmark tests consist of two different one element tests (tensile test and shear test), the necking of a circular bar and the well known Taylor impact test. The same material, a 42CrMo4 steel, has been selected for all those tests, and material properties are reported in Table 1. The proposed benchmarks are using the “dynamic temperature-displacement, Explicit” procedure of Abaqus explicit FEM code. The inelastic heat fraction parameters has been set to its default value of η = 0.9. This thermomechanical coupling option allows heat to be generated by plastic dissipation or viscoelastic dissipation. No other thermal boundary conditions are applied a globally adiabatic situation.

All benchmarks tests have been solved using Abaqus/Explicit v.6.14 on a Dell Precision T7500 computer running Ubuntu 16.04 64bits with 12 Gb of Ram and two 4 core E5620 Intel Xeon Processors. All computations have been done using the double precision option of Abaqus, with one CPU and the VUMAT and VUHARD subroutines have been compiled using the Intel Fortran 64 v.14 compiler. The following models have been tested for each of the four presented benchmarks:

• A-N-R model: VUMAT model with Newton–Raphson procedure and an analytical computation of the derivatives using equations (33)–(35);

• N-N-R model: VUMAT model with Newton–Raphson procedure and a numerical computation of the derivatives using equations (26)–(28);

• Direct model: VUMAT model with a direct evaluation of Γ using equation (19) as presented in other papers [2];

• VUHARD model: Only the Johnson–Cook constitutive flow law defined by equation (29) and its three analytical derivatives (33)–(35) have been implemented using a FORTRAN subroutine;

• Built-In model: native implementation of the Built-In Johnson–Cook constitutive law, in order to compare the results with a reference solution.

Tests have also been done using the Bisection method, but we have made the choice to omit them in this section to reduce the number of results because those results are very closed to the A-N-R model but computing time is 10–30% longer than the A-N-R and the N-N-R models. At this point, it has to be noted that, Abaqus/Explicit does not use the same objective stress rate when using the Built-In implementation and the VUMAT subroutine. As reported in documentation [1] and confirmed by our own tests based on a one element hyperelastic shear test, a Jaumann rate is used for both the Built-In formulation and the VUHARD subroutine while a Green–Naghdi rate is used for the VUMAT subroutine. Therefore, a straightforward comparison of Built-In and VUMAT results is not possible when large rotations occur.

5.1 One element benchmark tests

One element test, also called single element test, is a very simple and practical method to investigate the accuracy and sensitivity of the behavior of an element to the external loading. In this subsection, the deformation of a 4-node bilinear displacement and temperature, reduced integration with hourglass control element CPE4RT will be simulated using the proposed VUMAT subroutine and the Abaqus Built-In model. As only one under-integrated element is used, we only have one integration point located at the center of the element. In this kind of test, all nodes of the element are constrained with a prescribed displacement. As the geometry change is exactly the same in each of the two following benchmarks, it will be easy to compare the results in terms of plastic deformations, stresses and temperatures. The original size of the element is 10 mm × 10 mm.

5.1.1 One element tensile test

In this first benchmark, the two left nodes of the element are encastred and a prescribed horizontal displacement d = 10 mm is applied on the two right nodes of the same element as illustrated in Figure 3. As we are using an explicit integration scheme, the total simulation time is set to t = 0.01 s. A perfect match between all five results has been found for the one element tensile test, with a final value of the plastic strain and a final temperature T = 164.09 °C. Differences in computational time are not appreciable in this test since the final result is obtained in less than 1 s.

5.1.2 One element shear test

This second benchmark is similar to the previous one, but now, the two bottom nodes of the element are encastred and a prescribed horizontal displacement d = 10 mm is applied on the two top nodes of the same element as illustrated in Figure 3. Again, a perfect match between all five results has been found, with a final value of the plastic strain and a final temperature T = 192.22 °C.

5.2 Necking of a circular bar

The necking of a circular bar test is useful to evaluate the performance of VUMAT subroutine for materials in presence of plasticity and large deformation [12,8]. Because of the symmetric structure, an axisymmetric quarter model of the specimen is established. Dimensions of the specimen are reported in Figure 4. The loading is realized through an imposed total displacement of 7 mm along the axis on the left side of the specimen while the radial displacement of the same edge is supposed to remain zero. On the opposite side the axial displacement is restrained while the radial displacement is free. The mesh consists of 400 CAX4RT elements with a refined zone of 200 elements on the right side on 1/3 of the total height. Again, the total simulation time is set to t = 0.01 s because of the explicit integration scheme adopted.

Figure 5 shows the equivalent plastic strain contourplot of the deformed bar for two models: the Built-In model (left side) and the N-N-R model (right side). The maximum equivalent plastic strain is located into the center of the bar (so we have chosen to record the time-history evolutions for the red element in Fig. 4) and all the models give quite the same values as reported in Table 2 for , and T. Figure 6 shows the evolution of the von Mises stress with the displacement applied at the end of the specimen for the different models. As reported in this figure, the Built-In model, the VUHARD model and the two versions of the Newton–Raphson model give almost the same results. A slight difference can be seen in the Direct model, where there is a little over-estimation of the von Mises stress in the middle of the simulation.

Table 2 reports some results at the end of the computation concerning the total number of increments and the total computing time for all models. As reported in this table, the results of the two versions of the Newton–Raphson model are identical, this tends to prove that there is no difference in using numerical or analytical derivatives of the flow law. The T/I column show the ratio of the total computing time over the total number of increments, i.e. the computing time/increment. As reported in Table 2, the native procedure corresponds to the fastest algorithm in terms of computing time. This is easy to explain because this internal procedure is natively optimized within the Abaqus code. The simple fact of transferring all the data to a VUMAT or VUHARD subroutine increases the computing time independently of the optimization of the implemented stress evaluation algorithm. This translates in particular into the fact that the increase in CPU time due to the use of the VUHARD subroutine is about 18%, while the increase for the different versions of the VUMAT routines is about 12–15% in terms of T/I ratios.

It is noticeable in this case that the total number of increments needed to perform the whole simulation is lower for the Newton–Raphson models than for the three other models. This difference is illustrated in Figure 7 where the evolution of the time increment Δt during the computation is reported. The smoother variation of the time increment, and greatest values, are obtained with both versions of Newton–Raphson, leading to the minimal number of increments to complete the simulation. Using the Built-In model and the VUHARD model, a reduction of the stable time increment from Δt = 6.5 × 10⁻⁸ s to Δt = 5.2 × 10⁻⁸ s is noticed after a displacement of 2.52 mm, while using the Direct method, a large reduction of the stable increment to Δt = 3.4 × 10⁻⁸ s with some residual oscillations is encountered after a displacement 2.03 mm. The time increment using both versions of the Newton–Raphson model is constant up-to a displacement of 3.36 mm and reduces smoothly after this point because of the striction of the specimen. During the last third of the simulation time increments of the N-R models are almost equal to the one of the Built-In model and the VUHARD model.

From those results we can conclude that the integration of the constitutive equation has an influence on the stable time increment Δt. We can also note, from the results reported in Table 2 that, using the VUMAT Newton–Raphson does not increase a lot the total computational time (around 10% more time in this case) as the total number of increments has been reduced by a factor of 3.8%. We can also note that the VUMAT subroutine is faster than the VUHARD subroutine. Of course, computational cost of the VUMAT model can be reduced by optimizing the FORTRAN routine (removing the numerous number of tests inside the subroutine, optimizing the computations and the mathematical expressions, etc.). It has also to be noted that the requested precision for the Newton–Raphson subroutine ϵ_NR = 10⁻⁸ has also an influence on the computational time, reducing it reduces also the computational time.

5.3 Taylor impact test

5.3.1 Axisymmetric Taylor impact test

Finally, the performance of the proposed VUMAT subroutine is validated under high deformation rate with the simulation of the Taylor impact test [19]. In the Taylor impact test, a cylindrical specimen is launched to impact a rigid target with a prescribed initial velocity. The numerical model, reported in Figure 8 is established as axisymmetric. The height is 32.4 mm and the radius is 3.2 mm. The axial displacement is restrained on the right side of the specimen while the radial displacement is free (to figure a perfect contact without friction of the projectile onto the target). A predefined velocity of V_c = 287 m/s is imposed on the specimen. The mesh consists of 250 CAX4RT elements (5 × 50 elements). The total simulation time for the Taylor impact test is t = 8.0 × 10⁻⁵ s.

Figure 9 shows the equivalent plastic strain contourplot of the deformed rod for two models: the Built-In model (left side) and the N-N-R model (right side). The distributions of the equivalent plastic strain are almost the same for both models. The maximum equivalent plastic strain is located into the center element of the model (the red element in Fig. 8) and the models give quite the same value as reported in Table 3 for , T and the final dimensions of the specimen L_f (final length) and D_f (final diameter of the impacting face). Figure 10 shows the evolution of the equivalent plastic strain with time for the different models for the element at the center of the impacting face (the red element in Fig. 8). As reported in this figure and according to the results presented in Table 3, all models give almost the same results.

It can be noticed from this later that the Direct model (i.e. the classic way to program a VUMAT subroutine assuming that the explicit time integration scheme allows to use a straightforward evaluation of the plastic strain) does not provide results in accordance with the other models concerning the internal fields and the final geometry of the specimen. This tends to prove that if one wants to perform an inverse characterization of material constitutive parameters for dynamic applications based on a Taylor impact test and a FORTRAN VUMAT subroutine for implementing an exotic constitutive law, this can be source of errors in the identification of the constitutive equation parameters.

5.3.2 3D Taylor impact test

In order to have relatively longer computational times while remaining within the framework of the numerical simulation of the Taylor impact test, the choice was made for a 3D modeling. Therefore, a 3D quarter model of the Taylor cylindrical specimen, with the same dimensions as the Taylor 2D model, as been set-up as reported in Figure 11. The new mesh consists of 37 422 C3D8RT elements which is a quite large simulation benchmark.

The results concerning the CPU times (increments, time and T/I ratio), the final values of deformations and temperatures as well as the geometrical results are shown in Table 4. We can note a very good correlation of the results concerning geometric dimensions, whereas the results in terms of deformation and temperature (very localized results on a very constrained element) differ due to a better taking into account of the constitutive law by the VUMAT subroutine (this has been confirmed by refining the mesh). Similar to the results established in Section 5.2 concerning the necking of a circular bar, the T/I ratio is again increased by about 15% for the different versions of the VUMAT subroutines, whereas the overall increase in CPU time is only about 5% due to the reduction of the total number of increments required for simulation by about . Concerning the VUHARD approach, the increase in the T/I ratio is around 10% with approximately the same number of required increments as the Built-In method, leading to a total CPU increase of about 10% also. This proves again the efficiency of the proposed VUMAT algorithm.

5.4 VUMAT implementation of the alternative constitutive laws

In this section, the Johnson–Cook constitutive law and three other constitutive laws are implemented in the N-N-R VUMAT and the Direct VUMAT subroutines to simulate Taylor compression test, in order to validate the application of the proposed algorithm to alternative elastoplastic constitutive laws which follow J₂ plasticity and isotropic hardening and have the general form of . The main advantage of using our N-N-R approach is that it does not require to compute explicitly the derivatives of the hardening law, they are numerically computed using equations (26)–(28). The other constitutive laws are the TANH constitutive law, modified TANH constitutive law and Bäker constitutive law.

• Calamaz et al. [20] proposed the so called TANH constitutive law by adding a term modeling the strain softening to the Johnson–Cook constitutive law, given by:

(45) with (46)in which σ_JC represents the original Johnson–Cook constitutive law. The constant can modulate the strain corresponding to the peak stress, p and q are the additional constitutive law parameters, and T_rec is the onset temperature for the strain softening phenomenon.

• Hor et al. [21] proposed later a constitutive law by modifying the TANH constitutive law, given by:

(47) where (48)with (49)and A, B, C₁, C₂, m₁, m₂, n and p are the constitutive law parameters. As with the TANH model, the constant can modulate the strain corresponding to the peak stress and T_rec is the onset temperature for the strain softening phenomenon.

• the last constitutive law we have implemented using the VUMAT subroutine is the one proposed by Bäker [22,21] and given by:

(50) with (51)where A, n₀ and m are the temperature-dependent material parameters, and C and are constants.

The 2D model used for the simulation is the same as the one introduced in Section 5.3.1. The material used in this case is the 42CrMo4 steel with ferrito-perlitic (referred to as 42CrMo4-FP). The parameters of the four constitutive laws for 42CrMo4-FP steel were proposed by Hor et al. [21]. The predefined velocity is set to V_c = 200 m/s. The equivalent plastic strain contour-plots of the deformed rod for the four models are shown in Figure 12.

The plastic deformation is concentrated at the bottom, and maximum equivalent plastic strain is located in the center element of the specimen. More details concerning the total number of increments, total computing time, maximum equivalent plastic strain , final length L_f, final radius R_f of the bottom and maximum temperature T are reported in Table 5. The results of the Johnson–Cook and TANH models are almost identical. Compared with the previous two models, the modified TANH model achieves larger equivalent plastic strain at the bottom, and correspondingly it also achieves larger final radius and higher temperature at the bottom. However, as we can see, this model has less plastic deformation in the other part. No matter in terms of the maximum equivalent plastic strain and temperature or in terms of the geometric responses, the Bäker model has the least plastic deformation. Concerning the Direct approach, it can be seen that this one usually overestimates the plastic deformation and the temperature for the tree first models. In terms of global results, the three first models give almost the same results, while the Bäker model gives quite different results with regards to the previous ones.

5.5 Discussion on the benchmarks

In all the proposed benchmark tests, it is noteworthy that the results of the VUMAT subroutines and the Abaqus Built-In model are very closed but some slight differences can be pointed. Those differences can mainly be explained by the following remarks:

• the Johnson–Cook constitutive law implemented through the VUMAT subroutine does not have exactly the same expression as the Abaqus/Explicit Built-In model because the dependence on the deformation rate has been modified in our implementation as discussed before in Section 3. This can be seen in the Bar Necking benchmark where the VUHARD and the Built-In models does not provide exactly the same results;

• the Built-In model and the VUHARD model are integrated through an explicit central-difference time integration rule, while the radial return method, which belongs to an implicit integration algorithm, is employed in the VUMAT subroutines;

• as the root of the function is not an exact solution, but an approximation, the choice of the precision tolerance has an effect on the final results;

• the Jaumann stress rate is used for the Built-In model and the VUHARD subroutine while the VUMAT routine employs the Green–Naghdi stress rate. This can cause differences in the results in large deformations models when finite rotation of a material point is accompanied by finite shear [1].

6 Conclusions

In this paper, an approach has been proposed for the implementation of elastoplastic constitutive laws in Abaqus/Explicit through VUMAT subroutine. The proposed implementation is robust and efficient, and can be applied to simulate the behavior of various materials defined by any elastoplastic constitutive model following J₂ plasticity and isotropic hardening. In the proposed algorithm, an implicit time integration scheme based on the radial return method is employed, which is unconditionally stable and can ensure that the von Mises yield criterion is always satisfied during computation. Validated through numerical benchmarks, a robust and efficient root-finding method, the so-called safe version of Newton–Raphson method, has been implemented to compute the plastic corrector term.

This algorithm can now be re-used by only changing the expression of the flow law in order to implement alternative constitutive law, such as the revised form of the Johnson–Cook law proposed by Rule et al. [16], or any other one of the same family. Applications of this kind of developments mainly concerns inverse identification of constitutive law parameters based on a dynamic impact tests, or programming of alternative constitutive laws in Abaqus/Explicit. The proposed implementation has been found more precise with regards to other kind of VUMAT or VUHARD routines based on a direct evaluation of the plastic strain rate in order to provide accurate results for the identification procedure.

Acknowledgements

This work is partly supported by the scholarship from China Scholarship Council (CSC) under Grant CSC N0. 201406290010.

References

All Tables

All Figures