Implementation issues of Yld2000-2d model under larger biaxial yield stress

. In the ﬁ eld of sheet forming simulation, yield models serve as one of the most crucial factors for accurate computational results, and plane stress yield models have the capacity for both high ef ﬁ ciency and high accuracy. During recent years, applications of the Yld2000-2d model to sheet forming simulation of steel and aluminum have become increasingly popular due to its outstanding ability in describing these materials ’ yield phenomena. For the computational implementation of this model, the Newton – Raphson iteration can correctly obtain the solutions of return mapping equations in most cases. However, it has been found in this work that the traditional iteration process may fall into a convergence problem when the yield stress is prominent ( s b / s 0 > 1.2). To solve the new ﬁ nding problem, a line search algorithm is added to the Newton – Raphson iteration process. Biaxial tension simulation results show that the line search algorithm could converge successfully even when s b / s 0 = 1.4. The simulation of the Erichsen test shows the applicability of the established Yld2000-2d model combined with a line search algorithm in the Newton – Raphson iteration process.


Introduction
The sheet forming process is one of the most essential technologies in automobile and aviation industries [1]. The traditional trial-manufacture design method for the forming process is usually expensive and time-consuming. Fortunately, it has been possible to carry out the complicated calculation of the forming process in recent decades, because the performance of computers has been developing quickly since the second half of the 20th century. The calculation generally refers to the finite element (FE) model solution of the specific problem. Applications of accurate FE models may replace the traditional trial-manufacture method and decrease the design costs [2].
As the development of the FE method itself is satisfactory for plastic deformation problems nowadays, the primary computational error for the sheet forming process may stem from the selected constitutive model which is assumed and included in the FE model. A particular constitutive model of one material should be able to accurately describe the mechanical behavior of this material [3]. At present, three factors, including the yield model, the plastic flow rule, and the hardening law in the constitutive model, are still research hotspots. In this paper, the focus lies in the field of the yield model. It is appropriate to set the yield model as the plane stress state for the sheet forming process since the stress along the thickness direction of the sheet is negligible. The plane stress yield models can also decrease the computing time significantly compared with three-dimensional yield models.
The frequently-used plane stress yield models for sheet forming simulations are von Mises, Hill48 [4], Barlat89 [5], and Yld2000-2d [6]. Yueqi Wang et al. [7] applied the inverse method to identify the parameters of Hill48 for a DC06 sheet. Hedayati et al. [8] also conducted a similar investigation. Sudhy et al. studied the formability of an AA5754-H22 sheet at different temperatures using the Barlat89 model [9]. In recent years, the applications of the advanced yield model Yld2000-2d have become popular. The increasing utilization of the Yld2000-2d may be due to its excellent capacity for modeling the yield phenomena of aluminum and steel sheets. Myoung-Gyu et al. [10,11] used this model to investigate an automotive sheet spring back effect. Seoknyeon et al. [12] improved the formability prediction for advanced high strength steels by utilizing the Yld2000-2d model. Many other published papers also revealed the advantage of the Yld2000-2d over the von Mises, Hill48, and Barlat89 [13][14][15][16][17]. In the field of warm forming, many researchers have studied the application feasibility of Yld2000-2d [18][19][20][21].
In addition to the selection of a proper yield model, the implementation of the yield model into FE codes successfully and accurately is also quite crucial. Among all of the solution methods for the established FE model, an implicit approach with the superior accuracy is highlighted [22]. Consequently, in this work, a fully implicit algorithm is considered for the numerical integration of constitutive equations [3]. In the implicit algorithm, the Newton-Raphson iteration is usually applied to solve non-linear equation(s) with the constraint of the yield model. The von Mises model could converge unconditionally. Borst et al. [23] confirmed that the traditional Newton-Raphson (return mapping) algorithm could also ensure the convergence of Hill48, although the local curvature of this yield surface became stronger. Martin et al. [24] presented a line search procedure in the return mapping algorithm to get out of the convergence bowl caused by the larger exponent of Barlat89. For the Yld2000-2d, almost no related investigation has paid attention to the convergence issue of this model. It seems that there is no algorithmic problem for the implementation of Yld2000-2d using the traditional return mapping algorithm. Recently, Scherzinger [25] presented a line search algorithm in the Newton-Raphson iteration process for the implementation of Yld2004-18p and Hosford yield models. This robust method was previously proposed by Pérez-Foguet and Armero [26] for application to classical yield models. However, as described by Barlat et al. [27], neither Yld2004-18p nor Yld2004-13p could degenerate to Yld2000-2d. The algorithmic issues for Yld2000-2d should be investigated and understood separately. This paper has highlighted the implementation issues of Yld2000-2d under greater biaxial yield stress. The Newton-Raphson convergence problem is likely to arise with the implementation of the Yld2000-2d model when the curvature of the Yld2000-2d yield surface in the biaxial region is higher, i.e., when the value of the biaxial yield stress is sizeable. A line search algorithm is included in the Newton-Raphson iteration process to solve the new finding problem, and the influence of the biaxial yield stress on the convergence is analyzed. The results could be used by simulation engineers to predict the sheet forming results accurately.

Yld2000-2d model
The Yld2000-2d yield model proposed by Barlat et al. could be used for the sheet forming simulation of different materials [6]. Compared with Yld96 [28], the convexity of Yld2000-2d can be proven, and the numerical implementation of this model is convenient.
The yield condition of Yld2000-2d is where the model exponent M is a material coefficient, s is the uniaxial yield stress (effective stress) in the rolling direction and X 0 1 , X 0 2 , X 00 1 and X 00 2 are principal values of X 0 and X 00 where the subscripts x and y indicate the rolling and transverse directions of the sheet, and the components of X 0 and X 00 are where s is the real stress tensor, L 0 and L 00 are linear transformation matrices, and a 1 ∼ a 8 are the coefficients of Yld2000-2d. For the purpose of numerical implementation, it is more convenient to describe (1) equivalently as follows 'ðs; sÞ ¼ gðsÞ À s; ð4Þ where and In this paper, the associative flow rule is applied. Thus, the plastic strain rate is given by where _ g is the plastic multiplier and N is the plastic flow vector defined by in which As M = 6 and M = 8 are selected in most applications, ∂f ∂X 00 In (9) ∂X 0 The values of the factors in (11a) can be easily obtained according to (2a) and (3a).
When using the Newton-Raphson iteration in the implicit algorithm, the derivative of the flow vector should be acquired in which and ∂ ∂s The second-order partial derivatives on the right side of (14a) can be unfolded similar to (13), and the values can be finally attained.

Material
A 0.7 mm steel sheet DC06 is applied throughout this work. The hardening and anisotropy characteristics of this material were identified by standard uniaxial tensile tests [7]. Swift's model was used to describe the strain hardening phenomena s ¼ 548:57ð0:003 þ e p Þ 0:27 : ð15Þ The coefficients of Hill48 was obtained and further identified. The input data s 0 , s 45 , s 90 , s b and r 0 , r 45 , r 90 , r b of Yld2000-2d can be obtained according to the Hill48 model. The corresponding advanced yield model can thus be determined. Experimental evidence may be lacking for this operation. However, most features of Hill48 are reserved by Yld2000-2d, including s 0 , s 45 , s 90 , s b and r 0 , r 45 , r 90 , r b . This could be acceptable in view of the computation. The coefficients of Yld2000-2d are given in Table 1. A comparison of the two yield surfaces is shown in Figure 1. The curvature of the Yld2000-2d yield surface in the biaxial region is higher than that of the Hill48 yield surface. This is attributable to the features of the Yld2000-2d model, which include the material property of the biaxial yield stress. In contrast to the Yld2000-2d surface, the Hill48 surface still maintains moderate curvature in the biaxial region.
The Yld2000-2d exponent M is mainly related to the crystal structure of the material. M = 6 is applied throughout this paper. There are nearly no applications of larger values of M, which may result in higher curvature of the yield surface corners. An analysis of this case may be of no significance for Yld2000-2d. However, a larger biaxial yield stress may also lead to a higher curvature of Yld2000-2d. Larger values of the biaxial yield stress could be easily found in steel sheets such as DC05 [29]. In this work, the DC06 biaxial yield stress is 33% higher than the uniaxial yield stress in the rolling direction. The Newton-Raphson iteration is commonly applied to solve the non-linear equations. Barlat89 has difficulty converging when the model exponent is larger, which produces high curvature in the biaxial region [24]. A convergence problem may also be encountered for the implementation of the Yld2000-2d in the event of a larger biaxial yield stress.

Numerical integration for Yld2000-2d
The implicit numerical integration for the Yld2000-2d model into FE codes was first conducted by Yoon et al. [30] without considering the convergence problem mentioned above. In this section, a line search algorithm is added to the fully implicit return mapping algorithm. The line search algorithm was previously focused and applied by Pérez-Foguet and Armero [26] and Scherzinger [25].
At the initialization of an increment step [t n , t n+1 ], the values of elastic strain e e n , strain increment De and plastic strain e p n are given. The successful implementation of a yield model should find the unique solution e e nþ1 and e p nþ1 . There are three steps for the solution of e e nþ1 ; e p nþ1 , namely, the elastic trial step, the return mapping step and the line search step. In the first step, it is assumed that e e tr nþ1 ¼ e e n þ De; ð16Þ where the superscript tr is the abbreviation of trial, i.e., this step is assumed to be purely elastic. Accordingly where C in (18) is the plane stress elasticity matrix. If then this step is truly elastic, and the solution e e nþ1 ; e p nþ1 can be determined as The return mapping step is followed if (20) is not satisfied. Then, the return mapping (non-linear) equations can be established as e e nþ1 À e e tr nþ1 þ DgN nþ1 ¼ 0; ð22Þ where Dg is the difference form of the plastic multiplier, and the shear modulus m is taken into account to normalize the value of f n+1 [25]. It has been confirmed by Mohsen et al. [31] that the equivalent plastic strain increment is equal to the plastic multiplier if the associative flow rule is applied, i.e., e p nþ1 À e p n ¼ Dg. Then, (22) and (23) can be rewritten as gðs nþ1 Þ À sðe p n þ DgÞ The unknown values are s n+1 and Dg. The Newton-Raphson iteration is used to solve the non-linear equations above. The root update scheme after each Newton-Raphson iteration process is as follows where the superscript k means the number of the iteration process.
However, not all of the iterations can be converged (easily), as is presented in Section 2.2. Hence, a line search step is added to the Newton-Raphson iteration if this situation occurs. Each Newton-Raphson iteration process will be focused and examined, and the objective is to let the left sides of (24) and (25) approach zero as the iteration process continues. To express this idea mathematically, the following description is given in which The Newton-Raphson iteration without the line search step could also let the right side of (27) approach zero if there is no convergence problem. However, if the curvature of the yield surface corners is higher, the solved root may wander off as the iteration continues. In this situation, f still directs along a descent direction due to the Newton-Raphson iteration [32], whereas the decrease in f is limited, and its value may be far from zero, i.e., the local optimum may be encountered for the optimization of f using the traditional algorithm. To examine whether the Newton-Raphson iteration is effective, the following expression can be used [25,32] where b = 10 À4 is recommended. If (30) is satisfied, indicating that the decrease in f in this Newton-Raphson iteration process is sufficient, the iteration will go on. Otherwise, the line search algorithm is applied to revise the root update scheme after the Newton-Raphson iteration process as where a k ∈ (0, 1] denotes the modified step size. The problem that lets the left sides of (24) and (25) approach zero could then be transformed by finding a proper a k that could sufficiently decrease f k . The relationship between f k and a k is expressed as The derivative of (32) with respect to a k can be obtained as It is also found from (33) that the Newton-Raphson iteration process always finds a smaller f k because _ f k ða k Þ < 0 and a k = 1 regardless of the average decrease rate of f k . Consequently, the line search algorithm should find a proper a k between 0 and 1 to improve the average decrease rate of f k . As a k ∈ (0,1], the quadratic approximation of (32) in this interval can be written as Combined with (33), f k (a k ) turns into a minimum when : ð35Þ Then, f k+1 = f k (a k ), and (30) would be checked again. At this time, if the decrease in f is sufficient, the Newton-Raphson iteration process will go on.
If not, a k is updated as [25] a k jþ1 ¼ where j(j ≥ 1) denotes the iteration number of the updated a k . In addition, a k should not be too small and h = 0.1 is advised for (37). Another condition is introduced to estimate whether the new a k is feasible [26] At this point, the return mapping step and the line search step are both completed. The unknown values s n+1 and Dg can be obtained. Then, the required values e e nþ1 and e p nþ1 can be given as The robust implementation of Yld2000-2d is summarized in Box 1. The limit times of the Newton-Raphson iteration Newton is set as 50, and tol = 10 À6 is defined in this paper. 3 Results and discussion

Influence of the biaxial yield stress on the convergence
It is clear that the curvature of the Yld2000-2d yield surface corners is stronger under larger biaxial yield stress. There may be a convergence problem using the traditional algorithm if the biaxial tension stress state is prominent. In this section, the influence of the biaxial yield stress which is one of the material properties, on convergence is emphasized. It is tested with a simple biaxial tension numerical experiment. A 1/4 biaxial tension FE model is established, as is shown in Figure 2. There are only three square elements in this model to reduce the computation cost. The generalpurpose shell S4R with one integration point and hourglass control in Abaqus is selected. The length of each element side is 0.5 (standard units). Symmetric constraints are applied to the bottom and left sides of the FE model. To ensure the occurrence of plastic deformation, a 0.01 displacement is imposed on the top and right sides of the FE model at the same time. A local coordinate system is built to assign the material orientation, which is used in the User Material (UMAT).
The Hill48 model identified by Yueqi Wang et al. [7] is converted to Yld2000-2d, as described in Section 2.2. The relative value of s b is often changeable for different steel sheets, whereas the relative values of s 0 , s 45 , s 90 and r 0 , r 45 , r 90 , r b are mostly stable in realistic applications [7,29,33]. In this section the biaxial yield stress s b is given different values to obtain various Yld2000-2d surfaces. Seven levels of s b /s 0 and the corresponding coefficients of Yld2000-2d are given in Table 2. All of the yield surfaces are expressed in Figure 3. As s b /s 0 increases, the curvature of the yield surface becomes larger in the biaxial region.
The algorithm iteration times including the line search or not are expressed as a function of s b /s 0 in Figure 4. The iteration times are recorded once the first equilibrium iteration (the most difficult convergence iteration for both methods) is accomplished for the center element in Abaqus/Standard. Both iteration methods, i.e., the Newton-Raphson iteration and the Newton-Raphson iteration with the line search algorithm, successfully converge if s b /s 0 < 1.2. However, the iteration times of the traditional algorithm increase more quickly than the iteration times with the line search algorithm for this condition. When s b /s 0 is larger than 1.2, the Newton-Raphson method enters into its limit iteration times quickly without convergence. The convergence issue proposed in Section 2.2 can be explained here quantitatively due to the value of 1.33 for s b /s 0 . In contrast, the line search algorithm converges successfully even though s b /s 0 increases to 1.4. Additionally, the increasing speed of the iteration times is much slower as s b /s 0 becomes larger. Only 21 iterations are performed by this algorithm when s b /s 0 = 1.4.

Practical application of Yld2000-2d under larger biaxial yield stress
In this section, a complex numerical Erichsen test is applied to verify the practical performance of the robust algorithm. The biaxial tension stress state is prominent in this test. The FE model of the Erichsen test is built based on the work done by Yueqi Wang et al. [7]. The basic dimensions of the tools are shown in Figure 5a. A sufficient force is applied to the flange region to keep the sheet fixed during the bulging process. The punch moves 5.4 mm once it contacts the blank. The FE model is as shown in Figure 5b. 7416 S4R elements are used for the blank. The two blank holders are modeled as discrete rigid bodies, and the punch is modeled as analytical rigid body. Since the top surface deformation of the blank is concerned, the blank nodes are shifted from the default middle surface to the top side. The friction between the two holders and the blank is 0.7, and between the punch and the blank is 0.02 [7]. A local coordinate is also built to assign the material orientation. As is shown in Table 1, the value of s b /s 0 is 1.33. The convergence problem for the implementation of Yld2000-2d may occur when applying the traditional Newton-Raphson method. Therefore, the Newton-Raphson iteration with the line search algorithm is applied instead. The computation process is carried out by Abaqus/Standard. All of the iteration processes with the line search algorithm converge successfully for all of the blank elements. As expected, the convergence problem frequently occurs when using the traditional Newton-Raphson algorithm. The distribution of the equivalent plastic strain of the deformed blank is as shown in Figure 6. The accuracy of the computation result can be verified by the strain distribution comparison between the measurement and simulation, as is shown in Figure 7. The measurement data were obtained by digital image correlation technology. The measurement, friction coefficient, Yld2000-2d coefficient and the intrinsic FE computation error could account for the inconsistent part of Figure 7.

Conclusions
The algorithm problem for the implementation of Yld2000-2d model is highlighted in this paper. It is found that the traditional Newton-Raphson iteration process may fall into convergence problem with the implementation of the Yld2000-2d model when the biaxial yield stress is prominent. To solve this new finding problem, a line search algorithm is added to the Newton-Raphson method. The effect of the biaxial yield stress on the convergence is analysed through a simple biaxial tension test. When the value of s b /s 0  (a) measurement [7] (b) simulation is higher than 1.2, the traditional Newton-Raphson algorithm is mostly unable to achieve its solution, whereas the line search algorithm converges successfully even though s b /s 0 increases to 1.4. In terms of practical applications, the established Yld2000-2d model associated with the line search algorithm in the Newton-Raphson iteration process is successfully applied to the simulation of an Erichsen test in which the value of s b /s 0 is 1.33. The results presented in this paper may enrich the computational properties of Yld2000-2d and could be applied by simulation engineers in the automotive stamping department to accurately predict sheet forming results. Anisotropy coefficients along the rolling, diagonal, transverse and biaxial directions X 0 ; X 00 Linearly transformed stress tensor X 0 xx ; X 0 yy ; X 0 xy Components of X 0 X 0 xx ; X 0 yy ; X 0 xy Components of X 00 X 0 1 ; X 0 2 and X 00 1 ; X 00 2 Principles values of X 0 and X 00 Yld2000-2d yield function C Plane stress elasticity matrix