Issue 
Mechanics & Industry
Volume 19, Number 5, 2018



Article Number  501  
Number of page(s)  11  
DOI  https://doi.org/10.1051/meca/2018031  
Published online  03 December 2018 
Regular Article
Implementation issues of Yld20002d model under larger biaxial yield stress
School of Mechanical Engineering, Tongji University,
No. 4800 Caoan Road,
Shanghai
201804, PR China
^{*} email: sun1979@sina.com
Received:
12
January
2018
Accepted:
4
June
2018
In the field 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 efficiency and high accuracy. During recent years, applications of the Yld20002d 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 (σ_{b}/σ_{0} > 1.2). To solve the new finding 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 σ_{b}/σ_{0} = 1.4. The simulation of the Erichsen test shows the applicability of the established Yld20002d model combined with a line search algorithm in the Newton–Raphson iteration process.
Key words: Sheet metal forming / nonquadratic yield function / finite element method / line search
© AFM, EDP Sciences 2018
1 Introduction
The sheet forming process is one of the most essential technologies in automobile and aviation industries [1]. The traditional trialmanufacture design method for the forming process is usually expensive and timeconsuming. 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 trialmanufacture 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 threedimensional yield models.
The frequentlyused plane stress yield models for sheet forming simulations are von Mises, Hill48 [4], Barlat89 [5], and Yld20002d [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 AA5754H22 sheet at different temperatures using the Barlat89 model [9]. In recent years, the applications of the advanced yield model Yld20002d have become popular. The increasing utilization of the Yld20002d may be due to its excellent capacity for modeling the yield phenomena of aluminum and steel sheets. MyoungGyu 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 Yld20002d model. Many other published papers also revealed the advantage of the Yld20002d over the von Mises, Hill48, and Barlat89 [13–17]. In the field of warm forming, many researchers have studied the application feasibility of Yld20002d [18–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 nonlinear 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 Yld20002d, 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 Yld20002d using the traditional return mapping algorithm. Recently, Scherzinger [25] presented a line search algorithm in the Newton–Raphson iteration process for the implementation of Yld200418p and Hosford yield models. This robust method was previously proposed by PérezFoguet and Armero [26] for application to classical yield models. However, as described by Barlat et al. [27], neither Yld200418p nor Yld200413p could degenerate to Yld20002d. The algorithmic issues for Yld20002d should be investigated and understood separately.
This paper has highlighted the implementation issues of Yld20002d under greater biaxial yield stress. The Newton–Raphson convergence problem is likely to arise with the implementation of the Yld20002d model when the curvature of the Yld20002d 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.
2 Materials and methods
2.1 Yld20002d model
The Yld20002d 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 Yld20002d can be proven, and the numerical implementation of this model is convenient.
The yield condition of Yld20002d is $$\varphi ={\left{X\prime}_{1}{X\prime}_{2}\right}^{M}+{\left{X\mathrm{\u2033}}_{1}+2{X\mathrm{\u2033}}_{2}\right}^{M}+{\left2{X\mathrm{\u2033}}_{1}+{X\mathrm{\u2033}}_{2}\right}^{M}2{\overline{\mathit{\sigma}}}^{M}=0,$$(1)
where the model exponent M is a material coefficient, $\overline{\sigma}$is the uniaxial yield stress (effective stress) in the rolling direction and ${{X}^{\prime}}_{1}$, ${{X}^{\prime}}_{2}$, ${\mathit{X}\mathrm{\u2033}}_{1}$ and ${\mathit{X}\mathrm{\u2033}}_{2}$ are principal values of $X\prime $ and $X\mathrm{\u2033}$ $${X}_{1}=\frac{1}{2}({X}_{xx}+{X}_{yy}+\sqrt{{({X}_{xx}{X}_{yy})}^{2}+4{X}_{xy}^{2}})\text{,}$$(2a) $${X}_{2}=\frac{1}{2}({X}_{xx}+{X}_{yy}\sqrt{{({X}_{xx}{X}_{yy})}^{2}+4{X}_{xy}^{2}})\text{,}$$(2b)
where the subscripts x and y indicate the rolling and transverse directions of the sheet, and the components of $X\prime $ and $X\mathrm{\u2033}$ are $$\left[\begin{array}{c}{X\prime}_{xx}\\ {X\prime}_{yy}\\ {X\prime}_{xy}\end{array}\right]=\mathit{L}\prime \mathit{\sigma}=\left[\begin{array}{ccc}2{a}_{1}/3& {a}_{1}/3& 0\\ {a}_{2}/3& 2{a}_{2}/3& 0\\ 0& 0& {a}_{7}\end{array}\right]\mathit{\sigma},$$(3a) $$\left[\begin{array}{c}{X\mathrm{\u2033}}_{xx}\\ {X\mathrm{\u2033}}_{yy}\\ {X\mathrm{\u2033}}_{xy}\end{array}\right]=\mathit{L}\mathrm{\u2033}\mathit{\sigma}=\left[\begin{array}{ccc}(2{a}_{3}+2{a}_{4}+8{a}_{5}2{a}_{6})/9& ({a}_{3}4{a}_{4}4{a}_{5}+4{a}_{6})/9& 0\\ (4{a}_{3}4{a}_{4}4{a}_{5}+{a}_{6})/9& (2{a}_{3}+8{a}_{4}+2{a}_{5}2{a}_{6})/9& 0\\ 0& 0& {a}_{8}\end{array}\right]\mathit{\sigma},$$(3b)
where σ is the real stress tensor, $L\prime $ and $L\mathrm{\u2033}$ are linear transformation matrices, and a_{1} ∼ a_{8} are the coefficients of Yld20002d.
For the purpose of numerical implementation, it is more convenient to describe (1) equivalently as follows $$\phi (\sigma ,{\displaystyle \overline{\sigma}})=g(\sigma ){\displaystyle \overline{\sigma}}\text{,}$$(4)
where $$g(\sigma )={\left[f(\sigma )/2\right]}^{1/M}\text{,}$$(5)
and $$f(\sigma )={\phi}^{\prime}+{\phi}^{\u2033}={\left{{X}^{\prime}}_{1}{{X}^{\prime}}_{2}\right}^{M}+{\left{{X}^{\u2033}}_{1}+2{{X}^{\u2033}}_{2}\right}^{M}+{\left2{{X}^{\u2033}}_{1}+{{X}^{\u2033}}_{2}\right}^{M}\text{.}$$(6)
In this paper, the associative flow rule is applied. Thus, the plastic strain rate is given by $${{\displaystyle \dot{\mathit{\epsilon}}}}^{p}={\displaystyle \dot{\gamma}}N\text{,}$$(7)
where $\dot{\gamma}$ is the plastic multiplier and N is the plastic flow vector defined by $$N=\partial \phi /\partial \sigma =\partial g/\partial \sigma =\frac{1}{2M}{\left(\frac{f}{2}\right)}^{\frac{1M}{M}}\frac{\partial f}{\partial \sigma}\text{,}$$(8)
in which $$\frac{\partial f}{\partial \sigma}=\frac{\partial f}{\partial {{X}^{\prime}}_{1}}\frac{\partial {{X}^{\prime}}_{1}}{\partial \sigma}+\frac{\partial f}{\partial {{X}^{\prime}}_{2}}\frac{\partial {{X}^{\prime}}_{2}}{\partial \sigma}+\frac{\partial f}{\partial {{X}^{\u2033}}_{1}}\frac{\partial {{X}^{\u2033}}_{1}}{\partial \sigma}+\frac{\partial f}{\partial {{X}^{\u2033}}_{2}}\frac{\partial {{X}^{\u2033}}_{2}}{\partial \sigma}\text{.}$$(9)
As M = 6 and M = 8 are selected in most applications, $$\frac{\partial f}{\partial {{X}^{\prime}}_{1}}=M{({{X}^{\prime}}_{1}{{X}^{\prime}}_{2})}^{M1}\text{,}$$(10a) $$\frac{\partial f}{\partial {{X}^{\prime}}_{2}}=M{({{X}^{\prime}}_{1}{{X}^{\prime}}_{2})}^{M1}\text{,}$$(10b) $$\frac{\partial f}{\partial {{X}^{\u2033}}_{1}}=2M{(2{{X}^{\u2033}}_{1}+{{X}^{\u2033}}_{2})}^{M1}+M{({{X}^{\u2033}}_{1}+2{{X}^{\u2033}}_{2})}^{M1}\text{,}$$(10c) $$\frac{\partial f}{\partial {{X}^{\u2033}}_{2}}=M{(2{{X}^{\u2033}}_{1}+{{X}^{\u2033}}_{2})}^{M1}+2M{({{X}^{\u2033}}_{1}+2{{X}^{\u2033}}_{2})}^{M1}\text{.}$$(10d)
In (9) $$\frac{\partial {{X}^{\prime}}_{1}}{\partial \sigma}=\frac{\partial {{X}^{\prime}}_{1}}{\partial {{X}^{\prime}}_{xx}}\frac{\partial {{X}^{\prime}}_{xx}}{\partial \sigma}+\frac{\partial {{X}^{\prime}}_{1}}{\partial {{X}^{\prime}}_{yy}}\frac{\partial {{X}^{\prime}}_{yy}}{\partial \sigma}+\frac{\partial {{X}^{\prime}}_{1}}{\partial {{X}^{\prime}}_{xy}}\frac{\partial {{X}^{\prime}}_{xy}}{\partial \sigma}\text{,}$$(11a) $$\frac{\partial {{X}^{\prime}}_{2}}{\partial \sigma}=\frac{\partial {{X}^{\prime}}_{2}}{\partial {{X}^{\prime}}_{xx}}\frac{\partial {{X}^{\prime}}_{xx}}{\partial \sigma}+\frac{\partial {{X}^{\prime}}_{2}}{\partial {{X}^{\prime}}_{yy}}\frac{\partial {{X}^{\prime}}_{yy}}{\partial \sigma}+\frac{\partial {{X}^{\prime}}_{2}}{\partial {{X}^{\prime}}_{xy}}\frac{\partial {{X}^{\prime}}_{xy}}{\partial \sigma}\text{,}$$(11b) $$\frac{\partial {{X}^{\u2033}}_{1}}{\partial \sigma}=\frac{\partial {{X}^{\u2033}}_{1}}{\partial {{X}^{\u2033}}_{xx}}\frac{\partial {{X}^{\u2033}}_{xx}}{\partial \sigma}+\frac{\partial {{X}^{\u2033}}_{1}}{\partial {{X}^{\u2033}}_{yy}}\frac{\partial {{X}^{\u2033}}_{yy}}{\partial \sigma}+\frac{\partial {{X}^{\u2033}}_{1}}{\partial {{X}^{\u2033}}_{xy}}\frac{\partial {{X}^{\u2033}}_{xy}}{\partial \sigma}\text{,}$$(11c) $$\frac{\partial {{X}^{\u2033}}_{2}}{\partial \sigma}=\frac{\partial {{X}^{\u2033}}_{2}}{\partial {{X}^{\u2033}}_{xx}}\frac{\partial {{X}^{\u2033}}_{xx}}{\partial \sigma}+\frac{\partial {{X}^{\u2033}}_{2}}{\partial {{X}^{\u2033}}_{yy}}\frac{\partial {{X}^{\u2033}}_{yy}}{\partial \sigma}+\frac{\partial {{X}^{\u2033}}_{2}}{\partial {{X}^{\u2033}}_{xy}}\frac{\partial {{X}^{\u2033}}_{xy}}{\partial \sigma}\text{.}$$(11d)
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 $$\partial N/\partial \sigma =\partial \left[\frac{1}{2M}{\left(\frac{f}{2}\right)}^{\frac{1M}{M}}\frac{\partial f}{\partial \sigma}\right]/\partial \sigma =\frac{1}{2M}{\left(\frac{f}{2}\right)}^{\frac{1M}{M}}\frac{{\partial}^{2}f}{\partial {\sigma}^{2}}+\frac{1M}{4{M}^{2}}{\left(\frac{f}{2}\right)}^{\frac{12M}{M}}\frac{\partial f}{\partial \sigma}\otimes \frac{\partial f}{\partial \sigma}\text{,}$$(12)
in which $$\frac{{\partial}^{2}f}{\partial {\sigma}^{2}}=\frac{\partial}{\partial \sigma}\left(\frac{\partial f}{\partial {{X}^{\prime}}_{1}}\frac{\partial {{X}^{\prime}}_{1}}{\partial \sigma}+\frac{\partial f}{\partial {{X}^{\prime}}_{2}}\frac{\partial {{X}^{\prime}}_{2}}{\partial \sigma}+\frac{\partial f}{\partial {{X}^{\u2033}}_{1}}\frac{\partial {{X}^{\u2033}}_{1}}{\partial \sigma}+\frac{\partial f}{\partial {{X}^{\u2033}}_{2}}\frac{\partial {{X}^{\u2033}}_{2}}{\partial \sigma}\right)\text{,}$$(13)
and $$\frac{\partial}{\partial \sigma}\left(\frac{\partial f}{\partial {{X}^{\prime}}_{1}}\frac{\partial {{X}^{\prime}}_{1}}{\partial \sigma}\right)=\frac{\partial}{\partial \sigma}\left(\frac{\partial f}{\partial {{X}^{\prime}}_{1}}\right)\otimes \frac{\partial {{X}^{\prime}}_{1}}{\partial \sigma}+\frac{\partial f}{\partial {{X}^{\prime}}_{1}}\frac{\partial}{\partial \sigma}\left(\frac{\partial {{X}^{\prime}}_{1}}{\partial \sigma}\right)\text{,}$$(14a) $$\frac{\partial}{\partial \sigma}\left(\frac{\partial f}{\partial {{X}^{\prime}}_{2}}\frac{\partial {{X}^{\prime}}_{2}}{\partial \sigma}\right)=\frac{\partial}{\partial \sigma}\left(\frac{\partial f}{\partial {{X}^{\prime}}_{2}}\right)\otimes \frac{\partial {{X}^{\prime}}_{2}}{\partial \sigma}+\frac{\partial f}{\partial {{X}^{\prime}}_{2}}\frac{\partial}{\partial \sigma}\left(\frac{\partial {{X}^{\prime}}_{2}}{\partial \sigma}\right)\text{,}$$(14b) $$\frac{\partial}{\partial \sigma}\left(\frac{\partial f}{\partial {{X}^{\u2033}}_{1}}\frac{\partial {{X}^{\u2033}}_{1}}{\partial \sigma}\right)=\frac{\partial}{\partial \sigma}\left(\frac{\partial f}{\partial {{X}^{\u2033}}_{1}}\right)\otimes \frac{\partial {{X}^{\u2033}}_{1}}{\partial \sigma}+\frac{\partial f}{\partial {{X}^{\u2033}}_{1}}\frac{\partial}{\partial \sigma}\left(\frac{\partial {{X}^{\u2033}}_{1}}{\partial \sigma}\right)\text{,}$$(14c) $$\frac{\partial}{\partial \sigma}\left(\frac{\partial f}{\partial {{X}^{\u2033}}_{2}}\frac{\partial {{X}^{\u2033}}_{2}}{\partial \sigma}\right)=\frac{\partial}{\partial \sigma}\left(\frac{\partial f}{\partial {{X}^{\u2033}}_{2}}\right)\otimes \frac{\partial {{X}^{\u2033}}_{2}}{\partial \sigma}+\frac{\partial f}{\partial {{X}^{\u2033}}_{2}}\frac{\partial}{\partial \sigma}\left(\frac{\partial {{X}^{\u2033}}_{2}}{\partial \sigma}\right)\text{,}$$(14d)
The secondorder partial derivatives on the right side of (14a) can be unfolded similar to (13), and the values can be finally attained.
2.2 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 $$\overline{\sigma}}=548.57{(0.003+{{\displaystyle \overline{\u03f5}}}^{p})}^{0.27}\text{.$$(15)
The coefficients of Hill48 was obtained and further identified. The input data σ_{0}, σ_{45}, σ_{90}, σ_{b} and r_{0}, r_{45}, r_{90}, r_{b} of Yld20002d 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 Yld20002d, including σ_{0}, σ_{45}, σ_{90}, σ_{b} and r_{0}, r_{45}, r_{90}, r_{b}. This could be acceptable in view of the computation. The coefficients of Yld20002d are given in Table 1. A comparison of the two yield surfaces is shown in Figure 1. The curvature of the Yld20002d yield surface in the biaxial region is higher than that of the Hill48 yield surface. This is attributable to the features of the Yld20002d model, which include the material property of the biaxial yield stress. In contrast to the Yld20002d surface, the Hill48 surface still maintains moderate curvature in the biaxial region.
The Yld20002d 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 Yld20002d. However, a larger biaxial yield stress may also lead to a higher curvature of Yld20002d. 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 nonlinear 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 Yld20002d in the event of a larger biaxial yield stress.
Coefficients of Yld20002d for DC06.
Fig. 1 Comparison of the two yield surfaces Yld20002d and Hill48. X axis: Normalized σ_{11}, Y axis: Normalized σ_{22}, Continuous line: Yield locus of Yld20002, Dashed line: Yield locus of Hill48. 
2.3 Numerical integration for Yld20002d
The implicit numerical integration for the Yld20002d 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érezFoguet and Armero [26] and Scherzinger [25].
At the initialization of an increment step [t_{n}, t_{n+1}], the values of elastic strain ${\mathit{\epsilon}}_{n}^{e}$, strain increment $\mathrm{\Delta}\mathit{\epsilon}$ and plastic strain ${{\displaystyle \overline{\u03f5}}}_{n}^{p}$ are given. The successful implementation of a yield model should find the unique solution ${\mathit{\epsilon}}_{n+1}^{e}$ and ${{\displaystyle \overline{\u03f5}}}_{n+1}^{p}$. There are three steps for the solution of ${\mathit{\epsilon}}_{n+1}^{e},{{\displaystyle \overline{\epsilon}}}_{n+1}^{p}$, namely, the elastic trial step, the return mapping step and the line search step. In the first step, it is assumed that $${\u03f5}_{n+1}^{etr}={\u03f5}_{n}^{e}+\mathrm{\Delta}\u03f5\text{,}$$(16) $${{\displaystyle \overline{\u03f5}}}_{n+1}^{ptr}={{\displaystyle \overline{\u03f5}}}_{n}^{p}\text{,}$$(17)
where the superscript tr is the abbreviation of trial, i.e., this step is assumed to be purely elastic. Accordingly $${\sigma}_{n+1}^{tr}=\mathbb{C}:{\epsilon}_{n+1}^{e\text{\hspace{0.17em}}tr}\text{,}$$(18) $${{\displaystyle \overline{\sigma}}}_{n+1}^{tr}={\displaystyle \overline{\sigma}}({{\displaystyle \overline{\epsilon}}}_{n+1}^{p\text{\hspace{0.17em}}tr})\text{,}$$(19)
where $\mathbb{C}$ in (18) is the plane stress elasticity matrix. If $${\phi}^{tr}=\phi ({\sigma}_{n+1}^{tr},{{\displaystyle \overline{\sigma}}}_{n+1}^{tr})\le 0\text{,}$$(20)
then this step is truly elastic, and the solution ${\u03f5}_{n+1}^{e},{{\displaystyle \overline{\u03f5}}}_{n+1}^{p}$ can be determined as $${(\cdot )}_{n+1}={(\cdot )}_{n+1}^{tr}\text{.}$$(21)
The return mapping step is followed if (20) is not satisfied. Then, the return mapping (nonlinear) equations can be established as $${\epsilon}_{n+1}^{e}{\epsilon}_{n+1}^{etr}+\mathrm{\Delta}\gamma {N}_{n+1}=0\text{,}$$(22) $$\left[g({\sigma}_{n+1}){\displaystyle \overline{\sigma}}({{\displaystyle \overline{\u03f5}}}_{n+1}^{p})\right]/(2\mu )=0\text{,}$$(23)
where Δγ is the difference form of the plastic multiplier, and the shear modulus μ is taken into account to normalize the value of φ_{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., ${{\displaystyle \overline{\u03f5}}}_{n+1}^{p}{{\displaystyle \overline{\u03f5}}}_{n}^{p}=\mathrm{\Delta}\gamma $. Then, (22) and (23) can be rewritten as $${\mathbb{C}}^{1}:({\sigma}_{n+1}{\sigma}_{n+1}^{tr})+\mathrm{\Delta}\gamma {N}_{n+1}=0\text{,}$$(24) $$\left[g({\sigma}_{n+1}){\displaystyle \overline{\sigma}}({{\displaystyle \overline{\u03f5}}}_{n}^{p}+\mathrm{\Delta}\gamma )\right]/(2\mu )=0\text{.}$$(25)
The unknown values are σ_{n+1} and Δγ. The Newton–Raphson iteration is used to solve the nonlinear equations above. The root update scheme after each Newton–Raphson iteration process is as follows $${\sigma}_{n+1}^{k+1}={\sigma}_{n}^{k}+\mathrm{\Delta}\sigma \text{,}$$(26a) $$\mathrm{\Delta}{\gamma}^{k+1}=\mathrm{\Delta}{\gamma}^{k}+\mathrm{\Delta}(\mathrm{\Delta}\gamma )\text{,}$$(26b)
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 $${f}^{k}=\frac{1}{2}\left[{R}^{k}:{R}^{k}+{({\mathrm{\Phi}}^{k})}^{2}\right]\text{,}$$(27) in which $${R}^{k}={\u2102}^{1}:({\sigma}_{n+1}^{k}{\sigma}_{n+1}^{tr})+\mathrm{\Delta}{\gamma}^{k}{N}_{n+1}^{k}\text{,}$$(28) $${\mathrm{\Phi}}^{k}=\left[g({\sigma}_{n+1}^{k}){\displaystyle \overline{\sigma}}({{\displaystyle \overline{\u03f5}}}_{n}^{p}+\mathrm{\Delta}{\gamma}^{k})\right]/(2\mu )\text{.}$$(29)
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] $${f}^{k+1}<(12\beta ){f}^{k}\text{,}$$(30)
where β = 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 $${\sigma}_{n+1}^{k+1}={\sigma}_{n+1}^{k}+{\alpha}^{k}\mathrm{\Delta}\sigma \text{,}$$(31a) $$\mathrm{\Delta}{\gamma}^{k+1}=\mathrm{\Delta}{\gamma}^{k}+{\alpha}^{k}\mathrm{\Delta}(\mathrm{\Delta}\gamma )\text{,}$$(31b)
where α^{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 α^{k} that could sufficiently decrease f^{k}. The relationship between f^{k} and α^{k} is expressed as $${f}^{k}({\alpha}^{k})=\frac{1}{2}\left[{R}^{k}({\alpha}^{k}):{R}^{k}({\alpha}^{k})+{({\mathrm{\Phi}}^{k}({\alpha}^{k}))}^{2}\right]\text{,}$$(32)
The derivative of (32) with respect to α^{k} can be obtained as $${f}^{k}{\displaystyle \dot{(}}{\alpha}^{k})=2{f}^{k}({\alpha}^{k})\text{,}$$(33)
It is also found from (33) that the Newton–Raphson iteration process always finds a smaller f^{k} because ${f}^{k}{\displaystyle \dot{(}}{\alpha}^{k})<0$ and α^{k} = 1 regardless of the average decrease rate of f^{k}. Consequently, the line search algorithm should find a proper α^{k} between 0 and 1 to improve the average decrease rate of f^{k}. As α^{k} ∈ (0, 1], the quadratic approximation of (32) in this interval can be written as $${f}^{k}({\alpha}^{k})\approx [{f}^{k}(1){f}^{k}(0){f}^{k}{\displaystyle \dot{(}}0)]{({\alpha}^{k})}^{2}+{f}^{k}{\displaystyle \dot{(}}0){\alpha}^{k}+{f}^{k}(0)\text{.}$$(34)
Combined with (33), f^{k}(α^{k}) turns into a minimum when $${\alpha}^{k}=\frac{{f}^{k}(0)}{{f}^{k}(0)+{f}^{k}(1)}\text{.}$$(35)
Then, f^{k+1} = f^{k}(α^{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, α^{k} is updated as [25] $${\alpha}_{j+1}^{k}=\frac{{f}^{k}(0)}{{f}^{k}(0)+{f}^{k}({\alpha}_{j}^{k})}\text{,}$$(36)
where j(j ≥ 1) denotes the iteration number of the updated α^{k}. In addition, α^{k} should not be too small $${\alpha}_{j+1}^{k}=\mathrm{max}\left\{\eta {\alpha}_{j}^{k},{\alpha}_{j+1}^{k}\right\}\text{,}$$(37)
and η = 0.1 is advised for (37). Another condition is introduced to estimate whether the new α^{k} is feasible [26] $${f}_{j+1}^{k}({\alpha}_{j+1}^{k})<(12\beta {\alpha}_{j}^{k}){f}_{j}^{k}({\alpha}_{j}^{k})\text{.}$$(38)
At this point, the return mapping step and the line search step are both completed. The unknown values σ_{n+1} and Δγ can be obtained. Then, the required values ${\u03f5}_{n+1}^{e}$ and ${{\displaystyle \overline{\u03f5}}}_{n+1}^{p}$ can be given as $${\u03f5}_{n+1}^{e}={\u2102}^{1}:{\sigma}_{n+1}\text{,}$$(39a) $${{\displaystyle \overline{\u03f5}}}_{n+1}^{p}={{\displaystyle \overline{\u03f5}}}_{n}^{p}+\mathrm{\Delta}\gamma \text{.}$$(39b)
The robust implementation of Yld20002d 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.
Box 1 Numerical integration for Yld20002d. 
3 Results and discussion
3.1 Influence of the biaxial yield stress on the convergence
It is clear that the curvature of the Yld20002d 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 Yld20002d, as described in Section 2.2. The relative value of σ_{b} is often changeable for different steel sheets, whereas the relative values of σ_{0}, σ_{45}, σ_{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 σ_{b} is given different values to obtain various Yld20002d surfaces. Seven levels of σ_{b}/σ_{0} and the corresponding coefficients of Yld20002d are given in Table 2. All of the yield surfaces are expressed in Figure 3. As σ_{b}/σ_{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 σ_{b}/σ_{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 σ_{b}/σ_{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 σ_{b}/σ_{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 σ_{b}/σ_{0}. In contrast, the line search algorithm converges successfully even though σ_{b}/σ_{0} increases to 1.4. Additionally, the increasing speed of the iteration times is much slower as σ_{b}/σ_{0} becomes larger. Only 21 iterations are performed by this algorithm when σ_{b}/σ_{0} = 1.4.
Fig. 2 Mesh and load of the biaxial tension FE model. (a) Three elements of the 1/4 biaxial tension FE model. A global coordinate system and a local coordinate system, (b) Boundary and loading conditions of the 1/4 biaxial tension FE model. A global coordinate system and a local coordinate system. 
Coefficients of Yld20002d under different relative values of σ_{b}.
Fig. 3 Yld20002d yield surface under different relative values of σ_{b}. X axis: Normalized σ_{11}, Y axis: Normalized σ_{22}, Seven Yld20002d yield loci, from inside to outside corresponding to: σ_{b}/σ_{0} = 1.10, 1.15, 1.20, 1.25, 1.30, 1.35, 1.40. 
Fig. 4 The iteration times vs. σ_{b}/σ_{0} of different algorithm. X axis: Normalized σ_{b}, Y axis: Iteration times, Solid square line: The traditional Newton–Raphson algorithm, Solid circle line: Line search algorithm coupled with Newton–Raphson iteration. 
3.2 Practical application of Yld20002d 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 σ_{b}/σ_{0} is 1.33. The convergence problem for the implementation of Yld20002d 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, Yld20002d coefficient and the intrinsic FE computation error could account for the inconsistent part of Figure 7.
Fig. 5 The geometry and FE model of the Erichsen test. (a) Basic dimensions of the Erichsen test tools (two blank holders and a punch), (b) FE model FE model of the Erichsen test. A global coordinate system, a local coordinate system and axial symmetry line of the punch. 
Fig. 6 The computation result for equivalent plastic strain. The plastic strain distribution of the deformed blank. 
Fig. 7 Comparison of ϵ_{11} between measurement and simulation. (a) The distribution of measured ϵ_{11}, (b) The distribution of computational ϵ_{11}. 
4 Conclusions
The algorithm problem for the implementation of Yld20002d 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 Yld20002d 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 σ_{b}/σ_{0} 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 σ_{b}/σ_{0} increases to 1.4. In terms of practical applications, the established Yld20002d 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 σ_{b}/σ_{0} is 1.33. The results presented in this paper may enrich the computational properties of Yld20002d and could be applied by simulation engineers in the automotive stamping department to accurately predict sheet forming results.
Nomenclature
$a,b,c,f,\mathit{R},\mathrm{\Phi}$: Replacements of mathematical expressions
a_{1} ∼ a_{8}: Coefficients of Yld20002d
$g,\mathit{\varphi}\prime ,\mathit{\varphi}\mathrm{\u2033}$: Functions of σ
$\mathit{L}\prime ,\mathit{L}\mathrm{\u2033}$: Linear transformations matrices
r_{0}, r_{45}, r_{90}, r_{b}: Anisotropy coefficients along the rolling, diagonal, transverse and biaxial directions
$\mathit{X}\prime ,\text{\hspace{0.17em}}\mathit{X}\mathrm{\u2033}$: Linearly transformed stress tensor
${{X}^{\prime}}_{xx},\text{\hspace{0.17em}}{{X}^{\prime}}_{yy},\text{\hspace{0.17em}}{{X}^{\prime}}_{xy}$: Components of $\mathit{X}\prime $
${{X}^{\prime}}_{xx},\text{\hspace{0.17em}}{{X}^{\prime}}_{yy},\text{\hspace{0.17em}}{{X}^{\prime}}_{xy}$: Components of $\mathit{X}\mathrm{\u2033}$
$\mathit{X}{\prime}_{1},\text{\hspace{0.17em}}\mathit{X}\prime {}_{2}\text{\hspace{0.17em}and}\mathit{X}{\mathrm{\u2033}}_{1},\text{\hspace{0.17em}}\mathit{X}{\mathrm{\u2033}}_{2}$: Principles values of $\mathit{X}\prime $ and $\mathit{X}\mathrm{\u2033}$
$\beta ,\eta ,{\mathit{\in}}_{\mathrm{t}\mathrm{o}\mathrm{l}}$: Constants
$\dot{\gamma}},\mathrm{\Delta}\gamma $: Plastic multiplier and its difference form
${{\displaystyle \overline{\u03f5}}}^{p}$: Equivalent plastic strain
${{\displaystyle \dot{\mathit{\epsilon}}}}^{p}$: Plastic strain rate
$\mathrm{\Delta}\mathit{\epsilon}$: Strain increment
ϵ_{11}: A component of strain $\mathit{\epsilon}$
$\overline{\sigma}$: Uniaxial yield stress (effective stress), Pa
σ_{0}, σ_{45,} σ_{90,} σ_{b}: Yield stresses along the rolling, diagonal, transverse and biaxial directions, Pa
σ_{11}, σ_{22}: Components of σ, Pa
$\mathbb{C}$: Plane stress elasticity matrix
Acknowledgment
This work is supported by National Science and Technology Support Program of China (No. 2015BAF06B05), and the grant from the Science and Technology Commission of Shanghai Municipality, China (grant number 17DZ1204602).
References
 N.H. Kacem, N. Haddar, R. Elleuch, Failure analysis of an automotive shock absorber cup during manufacturing process, Mechanics & Industry 17 (2016) 604 [Google Scholar]
 L.B. Said, J. Mars, M. Wali, F. Dammak, Effects of the tool path strategies on incremental sheet metal forming process, Mechanics & Industry 17 (2016) 411 [Google Scholar]
 E.d.S. Neto, D. Peric, D. Owen, Computational methods for plasticity: theory and applications, Wiley, Chichester, 2008 [CrossRef] [Google Scholar]
 R. Hill, A theory of the yielding and plastic flow of anisotropic metals, Proc. R. Soc. Lond. A 193 (1948) 281–297 [CrossRef] [Google Scholar]
 F. Barlat, J. Lian, Plastic behavior and stretchability of sheet metals. Part I: a yield function for orthotropic sheets under plane stress conditions, Int. J. Plast. 5 (1989) 51–66 [CrossRef] [Google Scholar]
 F. Barlat, J.C. Brem, J.W. Yoon, K. Chung, R.E. Dick, D.J. Lege, F. Pourboghrat, S.H. Choi, E. Chu, Plane stress yield function for aluminum alloy sheets–part 1: theory, Int. J. Plast. 19 (2003) 1297–1319 [CrossRef] [Google Scholar]
 Y. Wang, S. Coppieters, P. Lava, D. Debruyne, Anisotropic yield surface identification of sheet metal through stereo finite element model updating, J. Strain Anal. Eng. Des. 51 (2016) 598–611 [Google Scholar]
 N. Hedayati, R. Madoliat, R. Hashemi, Strain measurement and determining coefficient of plastic anisotropy using digital image correlation (DIC), Mechanics & Industry 18 (2017) 311 [Google Scholar]
 S.S. Panicker, S. Kumar Panda, Improvement in material flow during nonisothermal warm deep drawing of nonheat treatable aluminum alloy sheets, J. Manuf. Sci. Eng. 139 (2016) 031013 [Google Scholar]
 M.G. Lee, D. Kim, C. Kim, M.L. Wenner, R.H. Wagoner, K. Chung, Springback evaluation of automotive sheets based on isotropickinematic hardening laws and nonquadratic anisotropic yield functions, Int. J. Plast. 21 (2005) 883–914 [Google Scholar]
 M.G. Lee, D. Kim, C. Kim, M.L. Wenner, K. Chung, Springback evaluation of automotive sheets based on isotropickinematic hardening laws and nonquadratic anisotropic yield functions, part III: applications, Int. J. Plast. 21 (2005) 915–953 [CrossRef] [Google Scholar]
 S. Kim, J. Lee, F. Barlat, M.G. Lee, Formability prediction of advanced high strength steels using constitutive models characterized by uniaxial and biaxial experiments, J. Mater. Process. Technol. 213 (2013) 1929–1942 [CrossRef] [Google Scholar]
 T. Kuwabara, T. Mori, M. Asano, T. Hakoyama, F. Barlat, Material modeling of 6016O and 6016T4 aluminum alloy sheets and application to hole expansion forming simulation, Int. J. Plast. 93 (2017) 164–186 [CrossRef] [Google Scholar]
 H. Tian, B. Brownell, M. Baral, Y.P. Korkolis, Earing in cupdrawing of anisotropic Al6022T4 sheets, Int. J. Mater. Form. 10 (2016) 329–343 [CrossRef] [Google Scholar]
 H.J. Choi, K.J. Lee, Y. Choi, G. Bae, D.C. Ahn, M.G. Lee, Effect of evolutionary anisotropy on earing prediction in cylindrical cup drawing, JOM 69 (2017) 915–921 [CrossRef] [Google Scholar]
 K. Charoensuk, S. Panich, V. Uthaisangsuk, Damage initiation and fracture loci for advanced high strength steel sheets taking into account anisotropic behaviour, J. Mater. Process. Technol. 248 (2017) 218–235 [CrossRef] [Google Scholar]
 M. Rossi, F. Pierron, M. Štamborská, Application of the virtual fields method to large strain anisotropic plasticity, Int. J. Solids Struct. 97–98 (2016) 322–335 [Google Scholar]
 R. Bagheriasl, K. Ghavam, M.J. Worswick, Formability improvement with independent die and punch temperature control, Int. J. Mater. Form. 7 (2014) 139–154 [CrossRef] [Google Scholar]
 K. Ghavam, R. Bagheriasl, M.J. Worswick, Analysis of nonisothermal deep drawing of aluminum alloy sheet with induced anisotropy and rate sensitivity at elevated temperatures, J. Manuf. Sci. Eng. 136 (2013) 011006 [Google Scholar]
 N. Abedrabbo, F. Pourboghrat, J. Carsley, Forming of AA5182O and AA5754O at elevated temperatures using coupled thermomechanical finite element models, Int. J. Plast. 23 (2007) 841–875 [CrossRef] [Google Scholar]
 Z. Cai, M. Wan, Z. Liu, X. Wu, B. Ma, C. Cheng, Thermalmechanical behaviors of dualphase steel sheet under warmforming conditions, Int. J. Mech. Sci. 126 (2017) 79–94 [CrossRef] [Google Scholar]
 M. G. Lee, C. Kim, E. J. Pavlina, F. Barlat, Advances in sheet forming—materials modeling, numerical simulation, and press technologies, J. Manuf. Sci. Eng. 133 (2011) 061001 [Google Scholar]
 R.D. Bors, P.H. Feenstra, Studies on anisotropic plasticity with reference to the Hill criterion, Int. J. Numer. Methods Eng. 29 (1990) 315–336 [Google Scholar]
 M. Dutko, D. Peric, D.R.J. Owen, Universal anisotropic yield criterion based on superquadric functional representation: Part 1. Algorithmic issues andaccuracy analysis, Comput. Methods Appl. Mech. Eng. 109 (1993) 73–93 [Google Scholar]
 W.M. Scherzinger, A return mapping algorithm for isotropic and anisotropic plasticity models using a line search method, Comput. Methods Appl. Mech. Eng. 317 (2017) 526–553 [Google Scholar]
 A. PérezFoguet, F. Armero, On the formulation of closestpoint projection algorithms in elastoplasticitypart II: globally convergent schemes, Int. J. Numer. Methods Eng. 53 (2002) 331–374 [Google Scholar]
 F. Barlat, H. Aretz, J.W. Yoon, M.E. Karabin, J.C. Brem, R.E. Dick, Linear transfomationbased anisotropic yield functions, Int. J. Plast. 21 (2005) 1009–1039 [Google Scholar]
 F. Barlat, Y. Maeda, K. Chung, M. Yanagawa, J.C. Brem, Y. Hayashida, D.J. Lege, K. Matsui, S.J. Murtha, S. Hattori, R.C. Becker, S. Makosey, Yield function development for aluminum alloy sheets, J. Mech. Phys. Solids 45 (1997) 1727–1763 [Google Scholar]
 P. Peters, N. Manopulo, C. Lange, P. Hora, A strain rate dependent anisotropic hardening model and its validation through deep drawing experiments, Int. J. Mater. Form. 7 (2013) 447–457 [CrossRef] [Google Scholar]
 J.W. Yoon, F. Barlat, R.E. Dick, K. Chung, T.J. Kang, Plane stress yield function for aluminum alloy sheets–part II: FE formulation and its implementation, Int. J. Plast. 20 (2004) 495–522 [CrossRef] [Google Scholar]
 M. Safaei, J.W. Yoon, W. De Waele, Study on the definition of equivalent plastic strain under nonassociated flow rule for finite element formulation, Int. J. Plast. 58 (2014) 219–238 [CrossRef] [Google Scholar]
 W.H. Press, S.A. Teukolsky, W.T. Vetterling, B.P. Flannery, Numerical recipes, Cambridge Unidersity, New York, 1986 [Google Scholar]
 M.C. Butuc, F. Barlat, J.J. Gracio, Study on plastic flow localization prediction using a physicallybased hardening model, Comput. Mater. Sci. 50 (2011) 2688–2697 [Google Scholar]
Cite this article as: Y. Zhang, Q. Zhang, Y. Sun, Implementation issues of Yld20002d model under larger biaxial yield stress, Mechanics & Industry 19, 501 (2018)
All Tables
All Figures
Fig. 1 Comparison of the two yield surfaces Yld20002d and Hill48. X axis: Normalized σ_{11}, Y axis: Normalized σ_{22}, Continuous line: Yield locus of Yld20002, Dashed line: Yield locus of Hill48. 

In the text 
Box 1 Numerical integration for Yld20002d. 

In the text 
Fig. 2 Mesh and load of the biaxial tension FE model. (a) Three elements of the 1/4 biaxial tension FE model. A global coordinate system and a local coordinate system, (b) Boundary and loading conditions of the 1/4 biaxial tension FE model. A global coordinate system and a local coordinate system. 

In the text 
Fig. 3 Yld20002d yield surface under different relative values of σ_{b}. X axis: Normalized σ_{11}, Y axis: Normalized σ_{22}, Seven Yld20002d yield loci, from inside to outside corresponding to: σ_{b}/σ_{0} = 1.10, 1.15, 1.20, 1.25, 1.30, 1.35, 1.40. 

In the text 
Fig. 4 The iteration times vs. σ_{b}/σ_{0} of different algorithm. X axis: Normalized σ_{b}, Y axis: Iteration times, Solid square line: The traditional Newton–Raphson algorithm, Solid circle line: Line search algorithm coupled with Newton–Raphson iteration. 

In the text 
Fig. 5 The geometry and FE model of the Erichsen test. (a) Basic dimensions of the Erichsen test tools (two blank holders and a punch), (b) FE model FE model of the Erichsen test. A global coordinate system, a local coordinate system and axial symmetry line of the punch. 

In the text 
Fig. 6 The computation result for equivalent plastic strain. The plastic strain distribution of the deformed blank. 

In the text 
Fig. 7 Comparison of ϵ_{11} between measurement and simulation. (a) The distribution of measured ϵ_{11}, (b) The distribution of computational ϵ_{11}. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
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.
Initial download of the metrics may take a while.