Stress update algorithm based on finite difference method and its application to homogenous anisotropic hardening (HAH) model

This paper deals with stress update algorithm based on finite difference method. The proposed algorithm is based on Euler backward method with multi-step return mapping approach. Central difference method was utilized for the approximation of the first and second derivatives of yield function. With the proposed algorithm, it is possible to perform elastic-plastic finite element simulation without any analytical derivatives of yield function. General yield functions including Hill’s (1948) and Yld2000-2d, and also the HAH distortion hardening model were implemented to an implicit finite element code (ABAQUS/STANDARD). For the verification purpose, various finite element simulations were performed. With the anisotropic functions, single element loading-unloading and cup-drawing simulation were carried out. The results obtained from the proposed algorithm were compared with the results from analytical derivatives and reference data. The availability of the proposed algorithm for distortional plasticity such as HAH model was evaluated by single element loading-reloading simulations: tension-compression and tension-orthogonal tension. The effectiveness of the proposed algorithm compared to the classical Euler backward method was identified from the simulation results.


Introduction
Stress update algorithm is an indispensable step to integrate the advanced constitutive models with rate or incremental form of equation. Among many algorithms, Euler backward method based on multi-step return mapping method [1] called hereafter the Euler backward method (analytical derivatives) has gained popularity because of its accuracy and unconditional stability. For example, by using the algorithm, Yoon et al. [2,3] successfully implemented Yld2000-2d [4] and Yld2004-18p [5] anisotropic functions under isotropic hardening assumption and they accurately predicted severe anisotropy of aluminium alloy. Recently, Lee et al. [6] derived the Euler backward method for the homogenous anisotropic hardening (HAH) model [7]. Throughout the stress update procedure using the Euler backward method, the calculation and implementation of first and second derivatives are the most time consuming and tedious components. In this work, finite difference method was applied to the Euler backward method to skip the calculation and implementation of any derivatives of yield function. Single element loading-unloading and cup-drawing simulations were performed with Hill's (1948) [8] and Yld2000-2d [4] anisotropic yield functions based on the proposed numerical algorithm and the Euler backward method (analytical derivatives). The accuracy and efficiency of the proposed algorithm were evaluated by comparing the simulation results with the two numerical and analytical methods. Finally, the strength of the proposed algorithm for the HAH model which has quite complicated forms of analytical derivatives was estimated with a single element loading-reloading simulation.

Stress update algorithm based on finite difference method (FDM)
In the Euler backward method (analytical derivatives), increment of the plastic multiplier is obtained at each sub step (k) and iteration (i) during stress update procedure as follows: and H is the instantaneous slope of hardening curve. Based on the updated at each iteration (i) and sub step (k), other variables are updated as: It is note that additional state variables (X) for the HAH model consisting of equivalent plastic strain could be updated. The stress update procedure will be end when the potential residual defined as Eq. (3) is satisfied within the numerical tolerance. Consistent tangent modulus is also defined to ensure the quadratic convergence rate of Newton`s iteration as follows: As shown in Eqs. (1) and (8), both first and second derivatives of yield function are required during the stress integration procedure. Central difference method was utilized for the approximation of the derivatives. With the method, both derivatives are approximated under the plane stress assumption as follows :   11  11  22  12  11  11  22  12 11 11  2  11  11  22  12  11  22  12  11  11  22  12  22  11  11   (  , , ) As indicated in Eqs. (9) and (10), three sets of stress components are required for the approximation and the way to define three sets during stress integration procedure strongly influences the convergence of FE simulation. The method to define three sets is explained with a simple example as below.
Let`s assume that the total number of sub steps is only one (N=1). Then, the stress ( 11 n  ) is updated by satisfying the consistency condition. For the definition during stress integration procedure, the stress lower limit should be defined based on the step size of stress as follows: where k indicates the corresponding sub step. Since the total number of sub steps is one in this demonstration, only one stress limit is defined, which is the same as the last converged stress ( n  ). For the first iteration at the step (i=1, k=1), the step size (


) is newly defined for three sets as follows: (13) Then, three sets for the first iteration are obtained as: By the three sets, both derivatives are approximated for the first iteration. If the updated stress ( 1.new  ) at the first iteration does not satisfy the convergence condition, the sets for the next iteration (i=2) are defined based on the step size for the next iteration as follows: Since the definition of the three sets is based on the step size of stress components, both denominator and numerator of the derivatives can be zero resulting in divergence of FE simulation. Therefore, with the numerical tolerance of step size of stress components, the derivatives within the tolerance step size should be zero as following example:

Single element loading-unloading simulation
The r-values of AA6022-T4E32 were predicted by a single element loading-unloading simulation. Hill's (1948) and

Cup-drawing simulation
Cup-drawing simulation was performed with two algorithms to predict the earing profile for AA2090-T3. The same geometry, tool dimension, and material properties were utilized as reported by Yoon et al. [3]. Figure 3 shows the predicted results based on Hill's (1948) and Yld2000-2d anisotropic functions, respectively. In the graph, black and red color indicate the simulation results based on each anisotropic yield functions. In addition, solid lines and symbols are the simulation results obtained from the Euler backward method (analytical derivatives) and the proposed numerical algorithm. As shown in the figure, the proposed algorithm shows the same level of accuracy as the Euler backward method (analytical derivatives). As for time efficiency, the proposed algorithm takes a longer simulation time as 1.2 times as the Euler backward method (analytical derivatives) for both constitutive models.

Single element loading-reloading simulation
Single element loading-reloading simulations were performed with two loading conditions: tensioncompression and tension-orthogonal tension. The original HAH model and 'Generic material' proposed by Barlat et al. [7] were utilized for the simulations. Figures 4(a) and 4(b) show the predicted true stresseffective plastic strain curves for each loading condition. The predicted curves well describe Bauschinger effect and material behavior after loading path change for both loading conditions. In addition, the distortion of yield surface during the simulations were well described as shown in figure 5.

Conclusions
Stress integration algorithm based on finite difference method was proposed. The accuracy and time efficiency was identified by comparing the simulation results obtained from the proposed algorithm with those obtained from the Euler backward method (analytical derivatives). The applicability of the proposed algorithm for distortional plasticity was successfully evaluated by applying the algorithm to the original HAH model. Considering accuracy and convenience for implementation, the proposed algorithm could be an efficient way to implement advanced constitutive models with moderate delay of simulation time.