Development of the subroutine library ‘UMMDp’ for anisotropic yield functions commonly applicable to commercial FEM codes

Numerous types of yield functions have been proposed to describe the shape of a realistic yield surface. Major commercial finite element codes include few anisotropic functions. Alternatively, the codes allow users to implement material models through user- subroutines. We develop the Unified Material Model Driver for Plasticity (UMMDp) subroutine library, which enables users to implement an arbitrary yield function easily. In this paper, the framework of the UMMDp is presented and its availabilities is shown through examples of sheet metal forming analyses.


Introduction
Yield functions generally represent the properties of plasticity in a numerical simulation. The yield functions under the associated flow rule describe not only relationship between the stress components in the plastic state, but also the direction of the plastic strain increment. Therefore, employment of appropriate anisotropic yield function, which can quantitatively express the realistic plasticity, is required for accurate prediction. Numerous types of yield functions have been proposed for this purpose so far in the past research [1]. But, general-purpose FE codes include few classical anisotropic yield functions. Alternatively, most FE codes allow users to implement their own material models as user-subroutines. However, it is not easy for the standard user to write subroutines by themselves, because detailed knowledge of continuum mechanics and numerical procedures is required to write subroutines appropriately.
The present paper describes the development of an easy-to-use user-subroutine library for anisotropic yield functions, which is widely applicable to major commercial FE codes.

Elasto-plastic constitutive model in FE codes
The classical FE codes for elasto-plastic deformation use a hypoelastic material model, and it is used for stress integration and calculation of the consistent tangent matrix. Although the above-mentioned roles of the material model are identical for all FE codes, detailed applications of subroutines in those codes are only discussed in closed communities such as user's conferences held by the software vendor. In order to promote industrial use of anisotropic yield functions, a universal easy-to-use subroutines package is needed to implement various yield functions for those codes.

Stress integration and consistent tangent matrix
In the backward Euler scheme, all of the constitutive equations are satisfied in the increment n+1. The state in which the stress is located on the isotropically hardened yield surface is described as 0 ) ( ) ( (1) where e σ is an arbitrary yield function, is the equivalent plastic strain, and ) ( p Y σ is the hardening function driven by p. By assuming the associated flow rule, the plastic strain increment } { p ε ∆ can be described as follows in the Voigt notation: The direction of the plastic strain increment coincides with the outward normal to the yield surface of increment n+1. Here, {N} is the 1st-order partial differential of e σ with respect to the stress. If the elastic deformation is small enough, the additive decomposition of the total strain increment can be applied as . Thus, the updated stress can be written respectively. These equations are linearized to obtain corrections } { σ δ∆ and p ∆ δ by Newton-Raphson iteration, as follows:   (1) and (3) are respectively.
From the above, the updated stress and the consistent tangent matrix can be calculated only if the yield function e σ and its 1st-and 2nd-order partial differentials, {N} and ] , respectively, are given.

Framework of the UMMDp subroutine library
As mentioned above, the roles of the material user-subroutines are identical for all FE codes, and the numerical procedures in the user-subroutines are independent of the type of yield functions. Based on these two features, the material model can be unified as a common numerical procedure, externalizing the variety of FE codes and the variety of yield functions. Figure 1 shows the framework of the Unified Material Model Driver for Plasticity (UMMDp). The core that executes the common numerical procedure is called from particular user-subroutines of FE codes. Each user-subroutine acts as a 'plug-in', which connects the UMMDp core to the FE codes. In the right branch, various yield functions are modularized as subroutines, in which e σ and its differentials are calculated as the arguments. The various yield functions were implemented in the UMMDp: Hill's 1948 and 1990, Gotoh's bi-quadratic (2D), Hu's bi-quadratic (3D), Yoshida's 6thorder polynomial, Bralat's Yld89, Yld2000-2d, and Yld2004, Banabic's BBC2005 and BBC2008, Cazacu's CPB2006, Karafills-Boyce, and Vegter's spline. Furthermore, as shown in the left branch, various hardening rules, including kinematic hardening, are modularized for future expansion.

Working example of the UMMDp
As an example, FEA for a hole-expansion test is performed using the UMMDp. In this test, the material anisotropy is emphasized in the thickness distribution around the hole. Figure 2 shows the experimental set-up of the test [2]. Five yield functions Yld2000-2d [3], Yld2004 [4], CPB2006 [5], Yoshida's 6th-order polynomial [6], and BBC2005 [7] are used in the simulations. The coefficients of these functions are determined based on the results of biaxial tensile tests using cruciform specimens [8]. Figure 3 shows the yield loci together with experimental plots of the normalized plastic work   contour in which the axes of anisotropy coincide with that of the principal stress. The yield loci of all functions are in agreement with the experimental results. Figure 4 shows the results of FEA performed using four commercial FE codes (Abaqus, ANSYS, LS-DYNA, and Marc) in conjunction with the same UMMDp. Figure 5 shows the circumferential distributions of the thickness strain along the edge of the hole. The results obtained using each FE code with Yld2004 were very similar. Also figure 4 shows the simulation results of various yield functions using Marc. The results differed depending on the type of yield function used in analyses.

Conclusion
A universal user-subroutine library for various anisotropic yield functions was developed. This library is available for major FE codes and enables users to easily implement their yield functions. Verification of the library is performed using an isotropic hardening model. The library is currently being expanded to include various kinematic hardening rules.