Issue 
yabo亚博
Volume19, Number3，2018 


Article Number  308  
Number of page(s)  15  
DOI  https://doi.org/10.1051/meca/2018021  
Published online  12 September 2018 
Regular Article
An efficient and robust VUMAT implementation of elastoplastic constitutive laws in Abaqus/Explicit finite element code
Laboratoire Génie de Production, Université de Toulouse, INPENIT, Toulouse, France
^{*}email:Olivier.Pantale@enit.fr
Received:19 May 2017
Accepted:2018年4月6日
This paper describes the development of an efficient and robust numerical algorithm for the implementation of elastoplastic constitutive laws in the commercial nonlinear finite element software Abaqus/Explicit through a VUMAT FORTRAN subroutine. In the present paper, while the Abaqus/Explicit uses an explicit time integration scheme, the implicit radial return mapping algorithm is used to compute the plastic strain, the plastic strain rate and the temperature at the end of each increment instead of the widely used forward Euler approach. This more complex process allows us to obtain more precise results with only a slight increase of the total computational time. Corrector term of the radial return scheme is obtained through the implementation of a safe and robust Newton–Raphson algorithm able to converge even when the piecewise defined hardening curve is not derivable everywhere. The complete method of how to implement a userdefined elastoplastic material model using the radial return mapping integration scheme is presented in details with the application to the widely used Johnson–Cook constitutive law. Five benchmark tests including one element tests, necking of a circular bar and 2D and 3D Taylor impact tests show the efficiency and robustness of the proposed algorithm and confirm the improved efficiency in terms of precision, stability and solution CPU time. Finally, three alternative constitutive laws (the TANH, modified TANH and Bäker laws) are presented, implemented through our VUMAT routine and tested.
关键词：Radial return algorithms / elastoplasticity / finite element / straindependent / thermomechanical
© 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 nonlinear 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 socalled consistent Jacobian matrix, but only the evaluation of the stress tensorσ_{1}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而不需要任何本地迭代算法[2，3]。
In this paper, a different approach has been chosen and consist in evaluating, in the VUMAT subroutine, the plastic strain incrementusing the implicit radial return algorithm presented inSection 2。In order to illustrate the proposed approach and to be able to validate it, the Johnson–Cook constitutive flow law [4]中提出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 inSection 4。数值效率和所提出的实现的精确度最终调查中，Section 5，through some dynamic benchmarks tests and have shown the efficiency and robustness of the proposed integration scheme.Section 5also present some alternative constitutive laws results for the Taylor impact test.
2 Time integration algorithm ofJ_{2}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。相对于在时间的时间导数的位置和一些t，the integration scheme calculates the same quantities at a later timet + Δt。Through the iteration of the procedure, the time evolution of those quantities can be traced.
2.1J_{2}plasticity model
The elastic stress–strain equation is usually written in the following incremental form:(1)where是一个客观压力为了考虑吗account the objectivity in large deformations, : is the double dot product,是由下式给出的线性各向同性第四阶弹性张量：(2)同the second order identity tensor,the fourth order identity tensor, ⊗ the Dyadic product,Gthe shear modulus andK体积模量挂钩的杨氏模量Eand Poisson’s ratiov通过以下：(3)and是the socalled elastic part of the rate of deformation tensorgenerally assumed [5] for hypoelastic material to conform to the following additive decomposition equation:(4)where是the plastic part of the rate of deformation tensor。According to the literature,J_{2}可塑性，也叫冯米塞斯屈服准则，在1904年，冯·米塞斯提出的胡贝尔在1913年[6]。It is widely applied in mechanical engineering and metal forming. It assumes the existence of a scalar yield functionf计算公式如下：（5）where是the von Mises equivalent stress, or effective stress, defined from the deviatoric stressas:（6）andσ^{y}是the socalled current yield stress of the material, depending in our kind of applications on the the equivalent plastic strain，the equivalent plastic strain rateand the temperatureT计算公式如下：（7）（8）whereη是the Taylor–Quinney [7]潜系数限定的塑料的工作的量转换成热能，C_{p}是the specific heat coefficient andρ是the density of the material. The socalled yield functionfdefines the yield surface (whenf = 0), the elastic domain (whenf < 0) and the noaccess domain (whenf > 0). Assuming an associative plastic flow rule, the plastic strain ratecan be expressed with the following equation:（9）whereγ是a scalar representing the flow intensity and是二阶张量（单位法向流动应力在试验弹性应力方面[专门确定8]) representing the flow direction given by:（10）
From equation（7）one can easily obtain:（11）
此处上述的塑性模型必须被集成在时间相对于增量客观算法。Cauchy应力张量的时间导数的不同的表情在方程(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 BuiltIn constitutive models and a Green–Naghdi stress rate in case of VUMAT user subroutines [1]。这可能会导致一些困难，当一个人想VUMAT结果与比较内置机型，我们将看得更远。
2.2 Radial return mapping algorithm
对的时间积分一种广泛使用的方法J_{2}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]或Ponthot [12]。We assume that, thanks to the finite element formulation, the strain increment between timet_{0}(the last equilibrated configuration) andt_{1}(the new equilibrated configuration witht_{1} = t_{0} + Δt) isand the deviatoric stress tensor at timet_{0}是。应力张的审判偏部分和最终的静水压力p_{1}are calculated from the strain incrementassuming that the whole step is fully elastic in a first time, so that:（12）and（13）
Conforming to equation（5），at timet_{1}，屈服函数f因此：（14）同:（15）andthe yield stress of the last equilibrated configuration at the beginning of the current increment (t = t_{0}）。In order to know if the predicted trial stress是admissible or not, a test has to be performed on the sign of the yield functionf导致两个选项后：
如果f ≤ 0, the whole step is fully elastic, which means that the predicted stress is admissible and we conclude that;
如果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 parameterf(Γ) = 0 at the end of the current increment (so that timet = t_{1}）：（16）同（17）
Combination of equations（16）and（17）leads to the following final form of the consistency parameterf(Γ):(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 timestep, 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)同，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 nonlinear terms of the constitutive flow law becomes important. In our approach, we have therefore chosen to solve equation(18)使用类似于在Zaera等人提出的一个的方法。[13] the so called root finding methods used to solve equation（11）将呈现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而更新的内部变量计算由于方程（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 Rootfinding of the nonlinear equation
Several rootfinding 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]。不幸的是，在保证收敛的方法被称为是较慢的，而那些表现出较快的收敛速度也能迅速冲至无穷远没有警告，如果不采取措施，避免这种行为。Among the large list of available methods (Chord, Bisection, RegulaFalsi, 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 rootfinding method presented in this paper is the socalled bisection method. This method is one that cannot fail, but with a slow rate of convergence. It’s an iterative rootfinding 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ε_{二}是二分法的容错。假设根被发现最初处于区间[括号x_{0}， x_{1}] satisfyingf(x_{0})f(x_{1}) < 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 functionf(x) itself and not its derivativef^{′}(x) as in other methods such as the Newton–Raphson. Therefore, this method can become competitive when the computation of the derivative of the functionf(x) becomes long to compute; i.e. when the expression off(x) is complex.
2.3.2牛顿迭代方法的安全版本
The Newton–Raphson method [14，15)是最著名的方法找到根由于e 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 derivativef^{′}(x）功能f(x), so it’s only usable in problems wheref^{′}(x) can be readily computed, or numerically evaluated as we will see further. From Taylor series expansion off(x) aboutx，we obtain the following expression:（22）
The Newton–Raphson is an iterative process starting with a first guessx_{0}该函数的一个根f(x) and the process is repeated until the convergence criterion given by:（23）where是牛顿迭代法的误差容限，到达。
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, whenf(x）是由两个不同的表达式上的给定点的左侧和右侧限定的分段函数x_{s}，f(x）将不会在微x_{s}，导致数值困难和牛顿迭代算法的分歧时，寻求解决方案x_{i}crosses the pointx_{s}。As a result, the socalled safe version of Newton–Raphson, described inFigure 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_{0}， x_{1}] satisfyingf(x_{0}) f(x_{1}) < 0, the safe version of Newton–Raphson method regards the midpoint of [x_{0}， x_{1})作为第一个猜的根和NewtonRaphson 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 reposition the initial guess. The Newton–Raphson algorithm requires to calculate the derivativeof the yield functionf求解方程时（Γ）相对于所述Γ参数(18)，so that:（24）用下面的定义为产率函数的导数：（25）
一种常见的方法来计算产率函数的导数σ^{y}同respect to，andT是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，andin 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^{−6}。
Fig. 1 流程安全NR算法。 
3 Hardening flow law
3.1 The Johnson–Cook hardening flow law
在本文中，约翰逊  库克硬化流动规律[4，16)已经被选择实现有限元分析/交货plicit 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 BuiltIn 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是given by the following equation:（29）where是the reference strain rate,T_{0}andT_{m}are the reference temperature and the melting temperature of the material respectively andA，B，C，nandmare 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是only taken into account if:(31)和对温度的依赖关系是defined so that, ifT < T_{0}there is no temperature dependence of the yield stress and ifT ≥ T_{m}the material is assumed to behave like liquid:(32)
Those two conditions defined by equations(31)and(32)导致一些不连续的硬化关系，its derivative且产函数本身导致迭代求解过程中的数值困难。因此，产率函数不是在可微andT_{0}。Nevertheless, the robust Newton–Raphson procedure proposed inSection 2.3.2has been found sufficiently robust to overcome those problems. Analytical forms for the derivatives of the Johnson–Cook flow lawσ^{y}同respect to，andTare given by the three following equations:（33）（34）（35）
The last problem usually encountered is that, because of the termin equation（33），htends 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 inSection 4.3。
3.2中的Abaqus VUHARD执行/显式
该VUHARD子程序是刚刚实施FORTRAN子程序来计算材料的屈服应力，以实现ABAQUS /新组成的流动法律明确一个简单的方法及其相对于衍生物，andT。内置的本构关系用于应力时间积分，对于给定的时间增量，以及所提供的用户子程序用于计算硬化流动规律的主要部分。很少给出细节有关的Abaqus文档在此实现，但一些有用的信息在詹森说伦斯堡等人都可用。[18]。
The numerical implementation of the VUHARD subroutine for the Johnson–Cook flow has been done through a FORTRAN program defining the hardening flow lawaccording to equation（29）and the three analytical derivatives ofσ^{y}同respect to，andTdefined 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 inFigure 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对于the current time step;
应力张量σ_{0}and the temperatureT_{0}at the beginning of the current increment;
the time increment Δt对应于当前时间步长;
溶液依赖状态变量（SDVs）的表用于存储诸如重要数据，，Γ从一个增量转移他们到其他。
The VUMAT subroutine must compute and return the value of the stressσ_{1}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.
Fig. 2 Flowchart of VUMAT implementation. 
4.1 Abaqus package subroutine
在软件的Abaqus的包的部分是用于计算的初始值（计算时间增量值的强制步骤Δt对于example by taking into account the material behaviour) for a reference timet = 0. For this we need to compute the elastic stress due to a virtual strain increment使用下面的表达式通过的Abaqus提供：（36）
The computed stressσ_{1}封装步骤结束之后将被丢弃。
4.2 Elastic predictor
本款的目的是计算预测弹性应力的偏部分和冯米塞斯等效应力and compare it with the yield stress at the beginning of the incrementin order to test if the current step is fully elastic or partly plastic. Therefore, we need firstly to decompose the initial stressσ_{0}into its hydrostatic pressure (p_{0}) and its deviatoric () parts:（37）
Then, we compute the new pressure (p_{1}) and the deviatoric part of the trial stress tensor () using:（38）而在增量开始时比较屈服应力to the von Mises equivalent trial stress(see Eq.（15）这个更新版本）。根据本次测试中，我们将那么正确与否最终压力。
如果，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 unchangedand we can go directly to the final computations inSection 4.4;
如果的步骤至少部分地塑料和塑料校正器中所描述Section 4.3has 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 (Γ = Γ_{0}）。如果当前的增量是第一塑料增量（即，如果), as we cannot compute the value ofh(Γ) when the plastic strain is zero, we initialize Γ to an initial value different of zero (Γ = 10^{−8}）。The interval for the bisection part is initialized so that。预测的等效塑性应变，plastic strain rate和温度T_{1}at the end of the increment are computed thanks to the system of equation（17）。屈服应力，屈服函数f（Γ）及其衍生物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:
如果ΔΓ > ε_{NR}we need to iterate to compute the correction of the Γ value using the following relation:
如果ΔΓ ≤ ε_{NR}we have obtained the final value of the Γ parameter.
应力张量的最后部分偏是computed from the predicted valueand the correction term Γ using:(41)同中，径向返回算法的框架内，通过定义：(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σ_{1}从流体静压计算p_{1}and the deviatoric part of the stress tensor由于公式(20)。The equivalent plastic strain，等效塑性应变率and final temperatureT_{1}at the end of the increment are stored in their respective SDVs variables for subsequent reuse in the next increment. We also have to compute the new specific internal energye_{1}from:(43)而消散能量无弹性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 Buildin routine.
所提出的实施方式的验证5
在本节中,约翰逊公司的性能ok law programmed in the VUMAT subroutine is compared with the performance of the Abaqus/Explicit BuiltIn 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 inTable 1。The proposed benchmarks are using the “dynamic temperaturedisplacement, 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:
ANRmodel: VUMAT model with Newton–Raphson procedure and an analytical computation of the derivatives using equations (33)–(35);
NNRmodel: VUMAT model with Newton–Raphson procedure and a numerical computation of the derivatives using equations (26)–(28);
Direct模式：与Γ用方程的直接评价VUMAT模型(19)as presented in other papers [2];
VUHARDmodel: Only the Johnson–Cook constitutive flow law defined by equation（29）及其三种分析衍生物（33）  （35）已使用FORTRAN子程序执行;
内建的型号：本机实现内置约翰逊  库克本构关系，为了比较的参考解决方案的结果。
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 ANR model but computing time is 10–30% longer than the ANR and the NNR models. At this point, it has to be noted that, Abaqus/Explicit does not use the same objective stress rate when using the BuiltIn implementation and the VUMAT subroutine. As reported in documentation [1]和由我们自己的基于一个元件超弹性剪切试验测试证实，一个焦曼速率被用于内置制剂和VUHARD子程序而绿Naghdi速率被用于VUMAT子程序两者。因此，当大的转动发生的内置和VUMAT结果的直接的比较是不可能的。
Material properties of the 42CrMo4 steel.
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 4node bilinear displacement and temperature, reduced integration with hourglass control element CPE4RT will be simulated using the proposed VUMAT subroutine and the Abaqus BuiltIn model. As only one underintegrated 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
在该第一基准，该元件的两个左节点encastred和规定的水平位移d= 10毫米两个正确的节点的应用same element as illustrated inFigure 3。由于我们使用一个明确的集成方案，总的模拟时间设置为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 strainand a final temperatureT = 164.09 °C. Differences in computational time are not appreciable in this test since the final result is obtained in less than 1 s.
Fig. 3 数值模型的一个元件的拉伸和剪切试验。 
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 displacementd = 10 mm is applied on the two top nodes of the same element as illustrated inFigure 3。Again, a perfect match between all five results has been found, with a final value of the plastic strainand a final temperatureT = 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 inFigure 4。The loading is realized through an imposed total displacement of 7 mm along theaxis 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 tot = 0.01 sbecause of the explicit integration scheme adopted.
Figure 5shows the equivalent plastic strain contourplot of the deformed bar for two models: the BuiltIn model (left side) and the NNR model (right side). The maximum equivalent plastic strain是located into the center of the bar (so we have chosen to record the timehistory evolutions for the red element inFig. 4) and all the models give quite the same values as reported inTable 2对于，andT。Figure 6shows the evolution of the von Mises stress与位移施加在样品的不同模型的端部。如该图报，内置的模型，该模型VUHARD和牛顿迭代模型的两个版本给几乎相同的结果。在微小的差别可以直接模型，其中有von Mises应力在模拟的中间多一点，估计可以看出。
Table 2reports 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. TheT/Icolumn show the ratio of the total computing time over the total number of increments, i.e. the computing time/increment. As reported inTable 2中，天然过程对应于最快的算法在计算时间方面。这是很容易解释，因为该内部过程在Abaqus的代码内本地优化。所有的数据传送到一个或VUMAT子程序VUHARD的简单事实独立地增加了实现的应力评估算法的优化的计算时间。这意味着特别是到一个事实，即由于使用VUHARD子程序的增加CPU时间大约是18％，而增加了不同版本的VUMAT例程的约1215％来讲T/Iratios.
正是在这种情况下显着，为执行整个仿真所需增量的总数是针对牛顿迭代模型比其他三个车型更低。这种差异是在所示Figure 7where the evolution of the time increment Δtduring 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 BuiltIn model and the VUHARD model, a reduction of the stable time increment from Δt = 6.5 × 10^{−8} s to Δt = 5.2 × 10^{−8} 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^{−8} 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 upto 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 NR models are almost equal to the one of the BuiltIn model and the VUHARD model.
从这些结果中我们可以得出结论，本构方程的集成对稳定时间增量的影响Δt。我们还可以注意到，从结果中报告Table 2that, 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^{−8}has also an influence on the computational time, reducing it reduces also the computational time.
Fig. 4 Numerical model for the necking of a circular bar. 
Fig. 5 Equivalent plastic straincontourplot for the necking of a circular bar. 
Comparison of results for the necking of a circular bar benchmark.
Fig. 6 Von Mises stressvs. displacement for the necking of a circular bar. 
Fig. 7 Time increment Δtvs. displacement for the necking of a circular bar. 
5.3泰勒冲击试验
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]。在泰勒冲击试验中，一个圆柱形试样推出与一个规定的初始速度撞击的刚性目标。该数值模型，报道Figure 8被确定为对称。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 ofV_{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 ist = 8.0 × 10^{−5} s.
Figure 9shows the equivalent plastic strain contourplot of the deformed rod for two models: the BuiltIn model (left side) and the NNR model (right side). The distributions of the equivalent plastic strain are almost the same for both models. The maximum equivalent plastic strain位于到模型的（中心元件的红色元件中Fig. 8) and the models give quite the same value as reported inTable 3对于，Tand the final dimensions of the specimenL_{f}（最终长度）和D_{f}（撞击面的最终直径）。Figure 10shows the evolution of the equivalent plastic strain随着时间对于不同的型号用于在冲击面的中心处的元件（红色元件中Fig. 8）。As reported in this figure and according to the results presented inTable 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.
Fig. 8 Model of the Taylor impact test. 
Fig. 9 Equivalent plastic straincontourplot for the Taylor impact test. 
Comparison of results for the Taylor impact test.
Fig. 10 Equivalent plastic strain对时间的泰勒冲击试验。 
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 setup as reported inFigure 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 andT/Iratio), the final values of deformations and temperatures as well as the geometrical results are shown inTable 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 in5.2节关于圆形杆的缩口，该T/Iratio 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 theT/Iratio is around 10% with approximately the same number of required increments as the BuiltIn method, leading to a total CPU increase of about 10% also. This proves again the efficiency of the proposed VUMAT algorithm.
Fig. 11 三维泰勒冲击试验模型。 
Comparison of results for the 3D Taylor impact test.
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 NNR 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 followJ_{2}plasticity and isotropic hardening and have the general form of。使用我们的NNR方法的主要优点是，它不需要计算明确地硬化法的衍生物，它们被使用等式数值地计算（26）  （28）。其他的本构关系是TANH组织法，修改TANH本构关系和贝克本构关系。
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:
Hor et al. [21]通过修改TANH本构关系，给出建议后，一个本构关系：
the last constitutive law we have implemented using the VUMAT subroutine is the one proposed by Bäker [22，21] and given by:
The 2D model used for the simulation is the same as the one introduced inSection 5.3.1。在这种情况下所使用的材料是42CrMo4钢ferrito  珠光体（简称42CrMo4FP）。四个本构关系42CrMo4FP钢的参数进行了贺等人提出。[21]。The predefined velocity is set toV_{c} = 200 m/s. The equivalent plastic strain contourplots of the deformed rod for the four models are shown inFigure 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 lengthL_{f}，final radiusR_{f}of the bottom and maximum temperatureTare reported inTable 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.
Fig. 12 Equivalent plastic strain contourplot of the Taylor specimen for the four constitutive models using the NNR VUMAT. 
Results for the Taylor compression test.
5.5论基准
In all the proposed benchmark tests, it is noteworthy that the results of the VUMAT subroutines and the Abaqus BuiltIn 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 BuiltIn model because the dependence on the deformation rate has been modified in our implementation as discussed before inSection 3。This can be seen in the Bar Necking benchmark where the VUHARD and the BuiltIn models does not provide exactly the same results;
内置的模型和VUHARD模型是通过显式的中心差分时间积分规则地结合，而所述径向返回方法，属于一个隐含的积分算法，采用在VUMAT子程序;
作为函数的根是不是一个确切的解决方案，但是近似，精度公差的选择对最终结果产生影响;
the Jaumann stress rate is used for the BuiltIn 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
本文提出了一种方法为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 followingJ_{2}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 rootfinding method, the socalled safe version of Newton–Raphson method, has been implemented to compute the plastic corrector term.
This algorithm can now be reused 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.
参考
 Abaqus v. 6.14 documentation, Dassault Systemes Simulia Corporation[Google Scholar]
 C. Gao, Fe realization of a thermoviscoplastic constitutive model using VUMAT in Abaqus/Explicit program, in: Computational Mechanics, Springer, 2007, pp. 301–301[Google Scholar]
 Abaqus v. 6.14 user subroutines reference guide, Dassault Systemes Simulia Corporation[Google Scholar]
 G.R。约翰逊，W.H.厨师，对经受大的应变，高应变速率和高温，在金属构模型和数据：第七届国际研讨会论文集上弹道，卷。21，海牙，荷兰，1983年，第541547[Google Scholar]
 S. NematNasser, On finite deformation elastoplasticity, Int. J. Solids Struct. 18 (1982) 857–872[Google Scholar]
 M.H. Yu, Generalized plasticity, Springer Science and Business Media, 2006[Google Scholar]
 G.I. Taylor, H. Quinney, The latent energy remaining in a metal after cold working, Proc. R. Soc. Lond. Ser. A, 143 (1934) 307–326 (Containing Papers of a Mathematical and Physical Character)[交叉引用][Google Scholar]
 J.C. Simo, T.J.R. Hughes, Computational inelasticity, Springer, 1998[Google Scholar]
 M.L.威尔金斯，弹性塑性流动，技术的计算。代表，在：B.桤木（编辑），的方法在计算物理卷。3，科学出版社，1963年，第211263[Google Scholar]
 g . Maenchen麻袋,张量的代码,:“肾上腺脑白质退化症”er, S. Fernback, M. Rotenberg (Eds.), Methods of Computational Physics, Vol. 3, Academic Press, New York, 1964, pp. 181–210[Google Scholar]
 F.邓恩，N. Petrinic，介绍计算可塑性，牛津大学出版社，2005年，纽约[Google Scholar]
 J.P. Ponthot, Unified stress update algorithms for the numerical simulation of large deformation elastoplastic and elastoviscoplastic processes, Int. J. Plast. 18 (2002) 91–126[交叉引用][Google Scholar]
 R. Zaera，J.费尔南德斯Sáez研究，隐式一致算法在绝热条件和有限变形，  传热诠释本构方程的整合。J.固体结构。43（2006）1594年至1612年[Google Scholar]
 F.S.阿克顿，数值方法的工作，MAA，1990年[Google Scholar]
 W.H. Press, Numerical recipes, in: The Art of Scientific Computing, 3rd Edition, Cambridge University Press, 2007[Google Scholar]
 W.K.规则，S.琼斯，修改后的形式约翰逊  库克强度模型，诠释。J.影响工程。21（1998）609624[Google Scholar]
 L. Schwer，可选的应变率形式的约翰逊  库克本构模型和参数小量0的作用，在：第六届欧洲LSDYNA用户大会，Anwenderforum，弗兰肯塔尔，德国，2007年的论文集[Google Scholar]
 G.詹森面包车伦斯堡，S.角，教程基于状态变量可塑性：一个ABAQUS UHARD子程序[Google Scholar]
 G.I. Taylor, The testing of materials at high rates of loading, J. Inst. Civ. Eng. 26 (1946) 486–519[交叉引用][Google Scholar]
 M. Calamaz, D. Coupard, F. Girot, Numerical simulation of titanium alloy dry machining with a strain softening constitutive law, Mach. Sci. Technol. 14 (2010) 244–257[交叉引用][Google Scholar]
 A. Hor, F. Morel, J.L. Lebrun, G. Germain, Modelling, identification and application of phenomenological constitutive laws over a large strain rate and temperature range, Mech. Mater. 64 (2013) 91–110[Google Scholar]
 M. BAKER，高速切割力的有限元模拟，J.母校。PROC。TECHNOL。176（2006）117126[交叉引用][Google Scholar]
引用本文为：L.明，O.Pantalé的一种高效率和强大的执行VUMAT在Abaqus的弹塑性本构关系/显式有限元软件，机械和工业yabo亚博19，308（2018）
All Tables
All Figures
当前使用指标显示文章点击的累积计数（全文文章的观点，包括HTML视图，PDF和ePub下载，根据现有数据）上Vision4Press平台和摘要视图。
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
指标初始下载可能需要一段时间。