Finite difference based stress integration algorithm for crystal plasticity finite element method

Abstract

In this study, we present a Finite Difference Method (FDM)-based stress integration algorithm for Crystal Plasticity Finite Element Method (CPFEM). It addresses the complexity of computing the first derivative of resolved shear stress in the Euler backward stress integration algorithm with Newton-Raphson method. The proposed FDM-based model was verified by evaluating its accuracy, convergence and computational efficiency through single-element simulations. The developed FDM-based model can be easily applied to various constitutive models for CPFEM, overcoming the problem of deriving complex derivative regardless of constitutive models. Additionally, the proposed FDM-based model was validated with the reduced texture approach using AA 2090-T3. Specific parameters including crystallographic orientations were calibrated and the plastic anisotropy was successfully described. In addition, the earing profiles were compared using various stress integration methods. As a result, the proposed FDM-based model can be used as an alternative to the Euler backward method using analytic derivatives with the compatible accuracy, convergence, computational efficiency along with easy implementation within the CPFEM framework.

Introduction

Crystal plasticity has been widely used to predict mechanical plastic properties. Many studies have been conducted especially to deal with anisotropy of metal since polycrystalline metal generally shows an anisotropic behavior due to a crystallographic texture [1, 2]. Five crystal plasticity models were compared for the anisotropy of AA1050, and the predicted anisotropic properties were used to characterize the continuum-level advanced yield functions [3]. The anisotropic properties of AA3103 under the cold-rolled (H18) and fully annealed (O) conditions were also studied using five crystal plasticity models, and the effect of thermo-mechanical processing conditions was investigated [4]. The Yld2004-18p anisotropic yield function was calibrated by the predictions from CPFEM and full-constraint Taylor model using the rolling and recrystallization textures [5]. In addition, the mechanical tests were replaced by CPFEM, and the mechanical properties from hot band and various cold rolled reduction conditions were obtained [6]. This multi-scale modeling was applied for the drawing processes, and the texture effect was investigated in terms of cold-rolled reductions. More applications of CPFEM were well summarized by Roters et al. [7].

A constitutive model with a stress integration method is required to implement crystal plasticity theory into finite element analysis [8,9,10,11,12,13,14,15]. In order to make it possible to update the material states, such as slip resistance and crystal lattice reorientation, various stress integration methods have been developed and can be classed as the Euler forward, semi-implicit, Euler backward (fully-implicit) methods. In the Euler forward method, the slip rates of active slip systems and elastic deformation gradients are predetermined from the last converged state, so they are all known quantities. In the semi-implicit method, only elastic deformation gradients are predetermined, and the slip rates are unknown in the current configuration [13, 16,17,18]. Therefore, the Newton-Raphson method is normally employed to satisfy the consistency conditions at the end of the time increment. In the Euler backward method, both the slip rates and elastic deformation gradients are unknown quantities in the current configuration. The Newton-Raphson method demands the Jacobian matrix which requires the derivative of the resolved shear stress at each iteration. Given that the exact derivative plays an important role in the convergence of CPFEM with a quadratic convergence rate, it is necessary to use the exact derivative. However, the stress integration methods normally include a tensor exponential function making the volume constant in plastic region. In this case, the calculation of tensor exponential function and its derivative is a challenging issue. Although the exact computation of exponential function and its derivative guarantees an accurate result with a superior convergence rate, many researchers have approximated the exponential function and its derivative to avoid the complicated calculation of them [19, 20]. The exact calculation of exponential function in single crystal plasticity was implemented by Miehe [21], and the exact calculation of its derivative was derived by de Souza Neto [22]. However, the analytical derivative of the resolved shear stress and its implementation are somewhat complicated. In addition, the derivative of the resolved shear stress highly depends on the constitutive model. Therefore, the stress integration method using the Finite Difference Method (FDM) was developed to skip the complex calculation of derivative of the resolved shear stress. Three finite difference methods were investigated in this work; the central, backward, and forward difference methods.

Crystal plasticity is a powerful tool for describing the mechanical behavior of metal. However, numerous grains exist in macroscale engineering applications, so it is not practical to employ the intensive information of crystallographic orientations. Therefore, it is desirable to implement the reasonable number of crystallographic orientations for computational efficiency while maintaining accuracy. In order to effectively reduce the computational cost, many researchers have introduced specific methods. Among them, the texture component method was employed for CPFEM [23,24,25,26]. In this method, a small number of predefined texture components obtained from mathematical schemes are considered instead of all crystallographic orientations. On the other hand, Rousselier et al. [27,28,29] developed the reduced texture methodology, where some sets of crystallographic orientations are calibrated until the satisfying description of anisotropic behavior is achieved. Moreover, the reduced texture methodology was employed to describe the multi-axial deformation of anisotropic metal with strength-differential effect [30] and ductile fracture for sheet metal [31]. Kim et al. [32] compared the results from the random mapping (RM) and misorientation mapping (MM) using five crystallographic orientations to investigate the earing profiles of baked hardening steel. In addition, the methodology for compacting crystallographic orientations to a smaller set of representative orientations was developed to decrease the computational cost [33], and generalized for any crystal structure [34]. In order to deal with a macroscale engineering application, the proposed FDM-based model was also employed with a small number of crystallographic orientations.

In the present paper, the semi-implicit method and Euler backward method using analytical derivative were summarized in a comprehensive way. The stress integration method using the FDM in the framework of Euler backward method was also newly developed. Several single element simulations with the various stress integration methods were compared in terms of accuracy, convergence, and computational efficiency. Furthermore, the time-efficient reduced texture approach was adopted to validate the proposed FDM-based model. The plastic anisotropy of AA2090-T3 was described using the reduced texture approach. Then, the circular cup drawing simulation was carried out, and the results from the FDM and analytical derivative were compared. It is shown that the proposed FDM-based model could be widely used for CPFEM.

Crystal plasticity model

Kinematics

The crystal plasticity model can predict the material deformation by the crystallographic slip and reorientation of crystal lattice. The basic kinematics of single crystal plasticity theory was proposed by Peirce et al. [11, 35]. It could be assumed that plastic deformation occurs due to the dislocation slip which is normally produced on the close-packed plane (slip plane) along the close-packed direction (slip direction). The plastic deformation generates irreversible shape change without any change of crystallographic orientation, while the elastic deformation leads to reversible crystallographic lattice deformation with the change of crystallographic orientation. Although the elastic and plastic deformations occur simultaneously, the assumption of intermediate configuration makes it possible to decompose the total deformation gradient ($\mathbf{F}$) into the elastic deformation gradient (${\mathbf{F}}^{e}$) and the plastic deformation gradient (${\mathbf{F}}^{p}$) as described in Eq. (1).

$$\mathbf{F}={\mathbf{F}}^{e}\cdot {\mathbf{F}}^{p}.$$

(1)

In the decomposition scheme, the elastic deformation gradient indicates the combination of rigid body motion and elastic deformation of the crystallographic lattice. The plastic deformation gradient means the summation of shear strains in each slip system. Figure 1 shows the decomposition of deformation gradient in the reference, intermediate, and spatial configurations. The intermediate configuration ($\overline{\Omega }$) is mapped from the reference configuration (${\Omega }_{0}$) by the plastic deformation gradient. The spatial configuration ($\Omega$) is mapped from the intermediate configuration by the elastic deformation gradient [36]. Similar to the decomposition of deformation gradient, the velocity gradient ($\mathbf{L}$) can also be divided into the elastic and plastic parts, i.e.,

Fig. 1

Multiplicative decomposition of deformation gradient

$$\mathbf{L}=\dot{\mathbf{F}}\cdot {\mathbf{F}}^{-1}=\left({\dot{\mathbf{F}}}^{e}\cdot {\mathbf{F}}^{p}+{\mathbf{F}}^{e}\cdot {\dot{\mathbf{F}}}^{p}\right)\cdot \left\{{\left({\mathbf{F}}^{p}\right)}^{-1}\cdot {\left({\mathbf{F}}^{e}\right)}^{-1}\right\}={\dot{\mathbf{F}}}^{e}\cdot {\left({\mathbf{F}}^{e}\right)}^{-1}+{\mathbf{F}}^{e}\cdot {\dot{\mathbf{F}}}^{p}\cdot {\left({\mathbf{F}}^{p}\right)}^{-1}\cdot {\left({\mathbf{F}}^{e}\right)}^{-1}={\mathbf{L}}^{e}+{\mathbf{L}}^{p}.$$

(2)

Since the plastic deformation is considered to be similar to the simple shear, the plastic velocity gradient (${\mathbf{L}}^{p}$) is the summation of shear strains over all the slip systems, so that the plastic velocity gradient can be given as:

$${\mathbf{L}}^{p}=\sum_{\alpha =1}^{NSYS}{\dot{\gamma }}^{(\alpha )}{\mathbf{b}}^{(\alpha )}\otimes {\mathbf{n}}^{(\alpha )},$$

(3)

where ${\mathbf{b}}^{(\alpha )}$ is the slip direction, ${\mathbf{n}}^{(\alpha )}$ is the normal direction of slip plane, $NSYS$ is the number of slip systems, and ($\alpha$) indicates the ($\alpha$)th slip system. The slip direction and normal direction of slip plane in the spatial configuration can be achieved from those in the reference configuration as follows:

$${\mathbf{b}}^{(\alpha )}={\mathbf{F}}^{e}\cdot {\mathbf{b}}_{0}^{(\alpha )},$$

(4)

$${\mathbf{n}}^{(\alpha )}={\mathbf{n}}_{0}^{(\alpha )}\cdot {\left({\mathbf{F}}^{e}\right)}^{-1}.$$

(5)

From Eq. (3), the symmetric and skew-symmetric parts of plastic velocity gradient can be defined as:

$${\mathbf{D}}^{p}=\frac{1}{2}\left({\mathbf{L}}^{p}+{\left({\mathbf{L}}^{p}\right)}^{T}\right)=\sum_{\alpha =1}^{NSYS}{\dot{\gamma }}^{(\alpha )}\frac{1}{2}\left({\mathbf{b}}^{(\alpha )}\otimes {\mathbf{n}}^{(\alpha )}+{\mathbf{n}}^{(\alpha )}\otimes {\mathbf{b}}^{(\alpha )}\right)=\sum_{\alpha =1}^{NSYS}{\dot{\gamma }}^{(\alpha )}{\mathbf{P}}^{(\alpha )},$$

(6)

$${\mathbf{W}}^{p}=\frac{1}{2}\left({\mathbf{L}}^{p}-{\left({\mathbf{L}}^{p}\right)}^{T}\right)=\sum_{\alpha =1}^{NSYS}{\dot{\gamma }}^{(\alpha )}\frac{1}{2}\left({\mathbf{b}}^{(\alpha )}\otimes {\mathbf{n}}^{(\alpha )}-{\mathbf{n}}^{(\alpha )}\otimes {\mathbf{b}}^{(\alpha )}\right)=\sum_{\alpha =1}^{NSYS}{\dot{\gamma }}^{(\alpha )}{{\varvec{\Omega}}}^{(\alpha )},$$

(7)

where ${\mathbf{P}}^{(\alpha )}$ and ${{\varvec{\Omega}}}^{(\alpha )}$ are introduced for convenience.

Constitutive model

It is assumed that the stress arises from the elastic deformation of crystallographic lattice. Therefore, the pull-back and push-forward by the elastic deformation gradient are used to ensure the objectivity of formulation. In order to formulate the constitutive model in the intermediate configuration ($\overline{\Omega }$), Cauchy strain tensor (${\overline{\mathbf{C}} }^{e}$) and Green strain tensor (${\overline{\mathbf{E}} }^{e}$) are defined as:

$${\overline{\mathbf{C}} }^{e}={\left({\mathbf{F}}^{e}\right)}^{T}\cdot {\mathbf{F}}^{e},$$

(8)

$${\overline{\mathbf{E}} }^{e}=\frac{1}{2}\left({\overline{\mathbf{C}} }^{e}-\mathbf{I}\right).$$

(9)

Second Piola-Kirchhoff stress tensor $\overline{\mathbf{S} }$ can be obtained by using Green strain tensor:

$$\overline{\mathbf{S} }={\mathcal{L}}^{{\varvec{e}}}:{\overline{\mathbf{E}} }^{{\varvec{e}}},$$

(10)

where ${\mathcal{L}}^{{\varvec{e}}}$ is the elasticity tensor.

Kirchhoff stress tensor (${\varvec{\uptau}}$) can also be given by the push-forward of second Piola-Kirchhoff stress:

$${\varvec{\uptau}}={\mathbf{F}}^{e}\cdot \overline{\mathbf{S} }\cdot {({F}^{e})}^{T}.$$

(11)

In this frame of the work, it is convenient to formulate the constitutive model in terms of Cauchy stress (${\varvec{\upsigma}}$) because User MATerial interface (UMAT) accepts Cauchy stress instead of Kirchhoff stress. However, it is still useful to take into account Kirchhoff stress since Cauchy stress can be approximated to Kirchhoff stress by ignoring the elastic volume change:

$${\varvec{\uptau}}=J{\varvec{\upsigma}}=\left(\mathrm{det }{\mathbf{F}}^{e}\right) \, {\varvec{\upsigma}}=\left(\frac{dV}{d{V}_{0}}\right){\varvec{\upsigma}}\approx {\varvec{\upsigma}}.$$

(12)

The objectivity should be considered when the constitutive model is employed. Jaumann rate is usually adopted as the objective stress rate in the constitutive models and defined as:

(13)

From Eq. (11), the rate of Kirchhoff stress can be written as:

$$\begin{array}{l}\dot{{\varvec{\uptau}}}={\mathbf{F}}^{e}\cdot \dot{\overline{\mathbf{S}} }\cdot {\left({\mathbf{F}}^{e}\right)}^{T}+{\dot{\mathbf{F}}}^{e}\cdot \left\{{\left({\mathbf{F}}^{e}\right)}^{-1}\cdot {\mathbf{F}}^{e}\right\}\cdot \overline{\mathbf{S} }\cdot {\left({\mathbf{F}}^{e}\right)}^{T}+{\mathbf{F}}^{e}\cdot \overline{\mathbf{S} }\cdot \left\{{\left({\mathbf{F}}^{e}\right)}^{T}\cdot {({\mathbf{F}}^{e})}^{-T}\right\}\cdot {\left({\dot{\mathbf{F}}}^{e}\right)}^{T}\\ \;\;\;={\mathbf{F}}^{e}\cdot \left[{\mathcal{L}}^{{\varvec{e}}}:\left\{{\left({\mathbf{F}}^{e}\right)}^{T}\cdot {\mathbf{D}}^{e}\cdot {\mathbf{F}}^{e}\right\}\right]\cdot {\left({\mathbf{F}}^{e}\right)}^{T}+{\mathbf{L}}^{e}\cdot {\varvec{\uptau}}+{\varvec{\uptau}}\cdot {\left({\mathbf{L}}^{e}\right)}^{T}.\end{array}$$

(14)

By substituting Eq. (14) into Eq. (13), Jaumann rate of Kirchhoff stress can be expressed as:

(15)

By writing the elastic part of deformation rate as the total and plastic parts, Jaumann rate of Kirchhoff stress becomes

(16)

Stress integration method

Two different stress integration methods were reviewed and extensively compared in “Semi-implicit method” and “Fully-implicit method (Euler backward method)” section: the semi-implicit method [13, 16] and fully-implicit method [9]. For the semi-implicit method, the Newton-Raphson method is utilized to find only the slip rates of slip systems, while the elastic deformation gradients come from the last convergence configuration. On the other hand, the Newton-Raphson method is used to find both the slip rates of slip systems and elastic deformation gradients in the fully-implicit method. In the fully-implicit method, the derivative of the resolved shear stress is also required to employ the Newton-Raphson method. The calculation of derivative is complicated to implement, and could be changed according to constitutive models. Therefore, the stress integration using the FDM for CPFEM was proposed to approximate the exact derivative of the resolved shear stress in “Fully-implicit method (Euler backward method) based on finite difference method” section, which allows for ease of implementation, regardless of the complexity of its derivative.

Semi-implicit method

In the semi-implicit method, the kinematics is assumed to evolve linearly over a time increment by using the forward gradient projection. The resolved shear stress on the slip system (α) can be determined as:

$${{\varvec{\uptau}}}^{(\alpha )}={\varvec{\uptau}}:{\mathbf{P}}^{(\alpha )}={\mathbf{n}}^{(\alpha )}\cdot {\varvec{\uptau}}\cdot {\mathbf{b}}^{(\alpha )}={\mathbf{b}}_{0}^{(\alpha )}\cdot{\varvec{\Sigma}}\cdot {\mathbf{n}}_{0}^{(\alpha )},$$

(17)

with

$${\varvec{\Sigma}}={\overline{\mathbf{C}} }^{e}\cdot \overline{\mathbf{S} }.$$

(18)

The derivative of Eq. (17) with respect to time yields

$${\dot{\tau }}^{(\alpha )}={\mathbf{R}}^{(\alpha )}:\mathbf{D}-{\mathbf{R}}^{(\alpha )}:\sum_{\beta =1}^{NSYS}{\dot{\gamma }}^{(\beta )}{\mathbf{P}}^{(\beta )}.$$

(19)

By using Eq. (19), the resolved shear stress can be updated at the end of the time increment as follows:

$${\tau }_{n+1}^{(\alpha )}={\tau }_{n}^{(\alpha )}+\Delta t\cdot {\dot{\tau }}^{(\alpha )}=\left[{\tau }_{n}^{(\alpha )}+\Delta t\cdot {\mathbf{R}}^{(\alpha )}:\mathbf{D}\right]-\Delta t\cdot {\mathbf{R}}^{(\alpha )}:\sum_{\beta =1}^{NSYS}{\dot{\gamma }}_{n+1}^{(\beta )}{\mathbf{P}}^{(\beta )}.$$

(20)

The term in the square bracket can be regarded as a trial state because ${\tau }_{n}^{(\alpha )}$ and ${\mathbf{R}}^{(\alpha )}$ are the known values and $\mathbf{D}$ is the input value. In other words, the trial state is directly updated to the next stress state when all the slip systems are inactive or have very small values of the slip rates. Besides, the updated resolved shear stresses ensure the consistency conditions for all slip systems at the end of the time increment:

$${\phi }_{n+1}^{(\alpha )}={\tau }_{n+1}^{(\alpha )}-{\left.{\tau }_{rss}^{(\alpha )}\right|}_{n+1}=\left[{\tau }_{n}^{(\alpha )}+\Delta t\cdot {{\varvec{R}}}^{(\alpha )}:{\varvec{D}}\right]-\Delta t\cdot {{\varvec{R}}}^{(\alpha )}:\sum_{\beta =1}^{NSYS}{\dot{\gamma }}_{n+1}^{(\beta )}{\mathbf{P}}^{(\beta )}-{\tau }_{rss}^{(\alpha )}\left({\dot{\gamma }}_{n+1}^{(\alpha )},{g}_{n+1}^{(\alpha )}\right).$$

(21)

Therefore, the unknown values (slip rates) are obtained after the iterative process based on the Newton-Raphson method for highly non-linear equations. Furthermore, Jacobian matrix is required when the Newton-Raphson method is employed, thus the derivative with respect to the slip rate is needed:

$${\phi }_{k}^{(\alpha )}+\sum_{\beta =1}^{NSYS}\frac{d{\phi }_{k}^{(\alpha )}}{d{\dot{\gamma }}^{(\beta )}}\cdot \delta {\dot{\gamma }}_{k}^{(\beta )}=0,$$

(22)

with

$$\frac{d{\phi }_{k}^{(\alpha )}}{d{\dot{\gamma }}^{(\beta )}}=-\left[\Delta t\cdot {\mathbf{R}}^{(\alpha )}:{\mathbf{P}}^{(\beta )}+\frac{\partial {\tau }_{rss}^{(\alpha )}}{\partial {\dot{\gamma }}^{(\beta )}}+\right.\sum_{\delta =1}^{NSYS}\left.\frac{\partial {\tau }_{rss}^{(\alpha )}}{\partial {g}^{(\delta )}}\frac{\partial {g}^{(\delta )}}{\partial {\dot{\gamma }}^{(\beta )}}\right].$$

(23)

In Eq. (22), the subscript k is the iteration number at each time increment and $\delta {\dot{\gamma }}_{k}^{(\beta )}$ is the correction added to the current estimate of slip rate to obtain the next estimate. The iterations are repeated until the residual is smaller than a numerical tolerance. Then, the grain level stress at (n + 1) time increment is updated as:

(24)

where $\mathbf{R}$ is the rotation tensor of material. Later, the grain level stresses are volume-averaged to achieve the integration point level stress in this study.

Fully-implicit method (Euler backward method)

In the fully-implicit method, both the slip rates and elastic deformation gradients are unknown quantities. Therefore, an iterative method is performed to ensure that the quantities satisfy the consistency conditions at the end of time increment. Therefore, the derivative of the resolved shear stress is required for the Newton-Raphson method. In this section, the derivative is fully summarized. Additionally, the elastic deformation gradients are directly updated to reduce the storage requirements because the elastic deformation gradients can evaluate the crystal reorientations.

Tensor exponential function

The exact computation of tensor exponential functions is recommended for high accuracy of the Euler backward method, and it can be defined by the series representation [21]:

$${\text{exp}}[\mathbf{X}]=\sum_{n=0}^{\infty }\frac{1}{n!}{\mathbf{X}}^{n}.$$

(25)

In addition, the tensor exponential function makes the volume constant during the plastic deformation. The determinant of tensor exponential function becomes [37]

$${\text{det}}[{\text{exp}}[\mathbf{X}]]={\text{exp}}[{\text{tr}}[\mathbf{X}]].$$

(26)

In Eq. (26), ${\text{det}}[{\text{exp}}[\mathbf{X}]]=1$ when ${\text{tr}}\left[\mathbf{X}\right]=0$. Given these relationships with Eq. (31), $det\left[{\mathbf{F}}_{n+1}^{p}\right]=1$ since ${\text{tr}}\left[\sum\limits_{\alpha =0}^{NSYS}{\dot{\gamma }}_{n+1}^{(\alpha )}\cdot {\mathbf{b}}_{0}^{(\alpha )}\otimes {\mathbf{n}}_{0}^{(\alpha )}\right]=0$. For the numerical computation of series expansion, a convergence criterion is required, and the convergence criterion is defined as follows:

$${\text{exp}}[\mathbf{X}]=\sum_{n=0}^{{n}_{{\text{max}}}}\frac{1}{n!}{\mathbf{X}}^{n},$$

(27)

with

$$\frac{1}{{n}_{{\text{max}}}!}\Vert {\mathbf{X}}^{{n}_{{\text{max}}}}\Vert <{\varepsilon }_{tol},$$

(28)

where ${\varepsilon }_{tol}={10}^{-10}$ was adopted in this paper. The exact derivative of tensor exponential function can also be expressed using the series representation [38]:

$${\left[\frac{\partial {\text{exp}}[\mathbf{X}]}{\partial \mathbf{X}}\right]}_{ijkl}=\sum_{n=1}^{{n}_{{\text{max}}}}\frac{1}{n!}\sum_{m=1}^{n}{[{\mathbf{X}}^{m-1}]}_{ik}{[{\mathbf{X}}^{n-m}]}_{lj},$$

(29)

where ${n}_{{\text{max}}}$ comes from Eq. (28).

Stress integration algorithm

By substituting Eq. (3) into Eq. (2), the derivative of plastic deformation gradient can be defined as:

$${\dot{\mathbf{F}}}^{p}=\left[\sum_{\alpha =1}^{NSYS}{\dot{\gamma }}^{(\alpha )}{\mathbf{b}}_{0}^{(\alpha )}\otimes {\mathbf{n}}_{0}^{(\alpha )}\right]\cdot {\mathbf{F}}^{p}.$$

(30)

With the backward exponential map integrator [37], the plastic deformation gradient at the current configuration (n + 1) can be defined as:

$${\mathbf{F}}_{n+1}^{p}={\text{exp}}\left[\sum_{\alpha =1}^{NSYS}\left(\Delta t\cdot {\dot{\gamma }}_{n+1}^{(\alpha )}\cdot {\mathbf{b}}_{0}^{(\alpha )}\otimes {\mathbf{n}}_{0}^{(\alpha )}\right)\right]\cdot {\mathbf{F}}_{n}^{p}.$$

(31)

Substituting Eq. (31) into Eq. (1) gives the updated elastic deformation gradient at the current configuration:

$${\mathbf{F}}_{n+1}^{e}={\mathbf{F}}_{n+1}\cdot {\left({\mathbf{F}}_{n+1}^{p}\right)}^{-1}={\mathbf{F}}_{n+1}\cdot {({\mathbf{F}}_{n})}^{-1}\cdot {\mathbf{F}}_{n}^{e}\cdot {\text{exp}}\left[\sum_{\alpha =1}^{NSYS}\left(-\Delta t\cdot {\dot{\gamma }}_{n+1}^{(\alpha )}\cdot {\mathbf{b}}_{0}^{(\alpha )}\otimes {\mathbf{n}}_{0}^{(\alpha )}\right)\right],$$

(32)

where ${\mathbf{F}}_{n+1}$ is the input value, and ${\mathbf{F}}_{n}^{e}$, ${\mathbf{F}}_{n}$ are the known quantities. Detailed stress integration algorithm in the fully-implicit method is summarized in Appendix 1.

Fully-implicit method (Euler backward method) based on finite difference method

Finite Difference Method (FDM)

It is complicated to compute the derivative of the resolved shear stress as derived in “Fully-implicit method (Euler backward method)” section. Also, other constitutive models are sometimes introduced, resulting in new calculation of the derivative. In order to avoid this issue, the analytical derivative of the resolved shear stress was directly approximated by using the numerical derivative: Finite Difference Method (FDM). The FDM-based model makes the implementation of conventional and new constitutive models easier. The numerical approach was also implemented for the return mapping method used in the continuum-level constitutive models for sheet metal forming [38, 39].

In this study, three finite difference methods were investigated. The central difference method was adopted since it has a higher-order truncation error. The backward and forward difference methods were also employed due to computational efficiency. For the central difference method, the first derivative of the resolved shear stress can be derived by using the Taylor series:

$$\begin{array}{l}{\tau }^{(\alpha )}\left({\dot{\gamma }}^{(1)},{\dot{\gamma }}^{(2)},\cdots ,{\dot{\gamma }}^{(\beta )}+\Delta {\dot{\gamma }}^{(\beta )},\cdots ,{\dot{\gamma }}^{(NSYS)}\right)={\tau }^{(\alpha )}\left({\dot{\gamma }}^{(1)},{\dot{\gamma }}^{(2)},\cdots ,{\dot{\gamma }}^{(\beta )},\cdots ,{\dot{\gamma }}^{(NSYS)}\right)\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;+\frac{\delta {\tau }^{(\alpha )}}{\delta {\dot{\gamma }}^{(\beta )}}\left({\dot{\gamma }}^{(1)},{\dot{\gamma }}^{(2)},\cdots ,{\dot{\gamma }}^{(\beta )},\cdots ,{\dot{\gamma }}^{(NSYS)}\right)\cdot\Delta {\dot{\gamma }}^{(\beta )}\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;+\frac{1}{2!}\frac{{\delta }^{2}{\tau }^{(\alpha )}}{\delta {\dot{\gamma }}^{(\beta )2}}\left({\dot{\gamma }}^{(1)},{\dot{\gamma }}^{(2)},\cdots ,{\dot{\gamma }}^{(\beta )},\cdots ,{\dot{\gamma }}^{(NSYS)}\right)\cdot {\left(\Delta {\dot{\gamma }}^{\left(\beta \right)}\right)}^{2},\end{array}$$

(33)

$$\begin{array}{l}{\tau }^{(\alpha )}\left({\dot{\gamma }}^{(1)},{\dot{\gamma }}^{(2)},\cdots ,{\dot{\gamma }}^{(\beta )}-\Delta {\dot{\gamma }}^{(\beta )},\cdots ,{\dot{\gamma }}^{(NSYS)}\right)={\tau }^{(\alpha )}\left({\dot{\gamma }}^{(1)},{\dot{\gamma }}^{(2)},\cdots ,{\dot{\gamma }}^{(\beta )},\cdots ,{\dot{\gamma }}^{(NSYS)}\right)\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;-\frac{\delta {\tau }^{(\alpha )}}{\delta {\dot{\gamma }}^{(\beta )}}\left({\dot{\gamma }}^{(1)},{\dot{\gamma }}^{(2)},\cdots ,{\dot{\gamma }}^{(\beta )},\cdots ,{\dot{\gamma }}^{(NSYS)}\right)\cdot\Delta {\dot{\gamma }}^{(\beta )}\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;+\frac{1}{2!}\frac{{\delta }^{2}{\tau }^{(\alpha )}}{\delta {\dot{\gamma }}^{(\beta )2}}\left({\dot{\gamma }}^{(1)},{\dot{\gamma }}^{(2)},\cdots ,{\dot{\gamma }}^{(\beta )},\cdots ,{\dot{\gamma }}^{(NSYS)}\right)\cdot {\left(\Delta {\dot{\gamma }}^{\left(\beta \right)}\right)}^{2}\cdots .\end{array}$$

(34)

By subtracting Eq. (34) from Eq. (33), the first derivative of the resolved shear stress can be achieved:

$$\frac{\partial {\tau }^{(\alpha )}}{\partial {\dot{\gamma }}^{(\beta )}}\approx \frac{{\tau }^{(\alpha )}\left({\dot{\gamma }}^{(1)},{\dot{\gamma }}^{(2)},\cdots ,{\dot{\gamma }}^{(\beta )}+\Delta {\dot{\gamma }}^{(\beta )},\cdots ,{\dot{\gamma }}^{(NSYS)}\right)-{\tau }^{(\alpha )}\left({\dot{\gamma }}^{(1)},{\dot{\gamma }}^{(2)},\cdots ,{\dot{\gamma }}^{(\beta )}-\Delta {\dot{\gamma }}^{(\beta )},\cdots ,{\dot{\gamma }}^{(NSYS)}\right)}{2\cdot\Delta {\dot{\gamma }}^{(\beta )}}$$

(35)

For the backward and forward difference methods, the first derivatives of resolved shear stress can also be expressed in Eqs. (36) and (37), respectively.

$$\frac{\partial {\tau }^{(\alpha )}}{\partial {\dot{\gamma }}^{(\beta )}}\approx \frac{{\tau }^{(\alpha )}({\dot{\gamma }}^{(1)},{\dot{\gamma }}^{(2)},\cdots ,{\dot{\gamma }}^{(\beta )},\cdots ,{\dot{\gamma }}^{(NSYS)})-{\tau }^{(\alpha )}({\dot{\gamma }}^{(1)},{\dot{\gamma }}^{(2)},\cdots ,{\dot{\gamma }}^{(\beta )}-\Delta {\dot{\gamma }}^{(\beta )},\cdots ,{\dot{\gamma }}^{(NSYS)})}{\Delta {\dot{\gamma }}^{(\beta )}}$$

(36)

$$\frac{\partial {\tau }^{(\alpha )}}{\partial {\dot{\gamma }}^{(\beta )}}\approx \frac{{\tau }^{(\alpha )}\left({\dot{\gamma }}^{(1)},{\dot{\gamma }}^{(2)},\cdots ,{\dot{\gamma }}^{(\beta )}+\Delta {\dot{\gamma }}^{(\beta )},\cdots ,{\dot{\gamma }}^{(NSYS)}\right)-{\tau }^{(\alpha )}\left({\dot{\gamma }}^{(1)},{\dot{\gamma }}^{(2)},\cdots ,{\dot{\gamma }}^{(\beta )},\cdots ,{\dot{\gamma }}^{(NSYS)}\right)}{\Delta {\dot{\gamma }}^{(\beta )}}$$

(37)

Equations (35)–(36) are the key equations to be implemented. Two additional resolved shear stresses are computed for the central difference method, while one additional resolved shear stress is required for the backward and forward difference methods. The Jacobian matrix consists of $12\times 12=144$ terms for FCC materials leading to 144 approximations per one Newton-Raphson iteration.

Step size for finite difference method

The step size ($\Delta {\dot{\gamma }}^{(\beta )}$) for the FDM is required to compute the derivative of the resolved shear stress, and it is critical for the FDM. In this research, a constant step size for all slip rates was used due to simplicity, and expressed in Eq. (38). Here, the step size for the FDM is different from the time increment size ($\Delta t$) in the finite element simulation.

$$\Delta {\dot{\gamma }}^{(\beta )}=\delta .$$

(38)

The reason why the step size was investigated is that a very small step size generally does not ensure the accuracy of the FDM. The accuracy of the FDM depends on the summation of the truncation and round-off errors. The truncation error results from the Taylor series expansion in Eqs. (33) and (34). For the central difference method in Eq. (35), the truncation error has the order of ${({\dot{\gamma }}^{(\beta )})}^{2}$. Therefore, the truncation error is reduced when the step size decreases. The forward and backward difference methods have the same tendency although the order of truncation error is ${\dot{\gamma }}^{(\beta )}$. On the other hand, the round-off error results from the limited storing space for numbers in computer, and the round-off error becomes larger when the step size decreases. Therefore, the truncation error from the series expansion and the round-off error from the numerical problem are inversely proportional to each other. Additionally, in order to check whether the step size depends on the input value, two cases were investigated.

In Case 1, ${\mathbf{F}}_{n}$ was an identity and ${\mathbf{F}}_{n+1}$ was set to cause 10% uniaxial tensile state. For the crystal orientation, a single crystal with a cube texture, $\left({\varphi }_{1},\phi ,{\varphi }_{2}\right)=\left({0}^{\circ },{0}^{\circ },{0}^{\circ }\right)$, was used. Also, ${\dot{\gamma }}^{\left(\alpha \right)}=({10}^{-10},0.5,-0.5,{10}^{-10},-\mathrm{0.5,0.5},-\mathrm{0.5,0.5},$ ${10}^{-10},0.5,{10}^{-10},-0.5)$ was employed for the slip rates because the 2,6,8 and 10th slip systems are activated in the positive direction, and the 3,5,7 and 12th slip systems are activated in the negative direction. The twelve slip systems for FCC materials are shown in Fig. 2. For the elasticity tensor constants, $\left({C}_{1111},{C}_{1122},{C}_{1212}\right)=\left(108000,62000,28300\right)$ [MPa] was used [40]. In Case 2, only ${\mathbf{F}}_{n+1}$ was changed to cause 0.1% uniaxial tension to investigate the size effect of the input value.

Fig. 2

Twelve slip systems for FCC materials

Figure 3 shows the average relative errors of the FDM (144 terms) with respect to the analytical derivative according to the step size and FDM type; the central, backward, and forward difference methods. The results show two characteristics. First, the error decreases to a specific step size value but increases as the step size further decreases due to the truncation and round-off errors. Second, a specific step size for the minimum error does not exist. Therefore, the constant step size approach with $\delta ={10}^{-5}$ was adopted for the central difference method, and $\delta ={10}^{-7}$ was employed for the backward and forward difference methods considering the minimum errors in Fig. 3

Fig. 3

Averaged relative errors of resolved shear stress derivatives by three difference methods: (a) Central difference method; (b) Backward difference method; (c) Forward difference method

Comparison of stress integration methods

Numerical examples in this paper are based on the face-centered cubic (FCC) lattice which has twelve slip systems. For the polycrystalline aggregates, the Taylor assumption was employed. The main objective in this section is the comparison of various stress integration methods including the proposed FDM-based model. Therefore, the semi-implicit, Euler backward, and proposed FDM-based methods were compared in terms of accuracy, convergence, and computational efficiency. All simulations were performed by using a commercial software Abaqus2016/Standard with UMAT. When UMAT for Abaqus is developed and the local orientation is employed, the input form of deformation gradient should be precisely figured out. The input forms for the deformation gradient depending on the local orientation were well explained by Nolan et al. [41].

Accuracy

Uniaxial tensile and simple shear simulations were performed in this section. The results from various stress integration methods were compared in terms of accuracy. For the crystal orientation, a Copper texture component, $\left({\varphi }_{1},\phi ,{\varphi }_{2}\right)=\left({90}^{\circ },{35}^{\circ },{45}^{\circ }\right)$, was employed. There are two types of hardening models: phenomenological [11, 13, 42] and physical-based models [43,44,45,46,47]. In this paper, a phenomenological model was used. The resolved shear stress can be defined as:

$${\tau }_{rss}^{(\alpha )}={g}^{(\alpha )}{\left(\frac{\left|{\dot{\gamma }}^{(\alpha )}\right|}{{\dot{\gamma }}_{0}}\right)}^{m}sign({\dot{\gamma }}^{(\alpha )}),$$

(39)

where ${g}^{\left(\alpha \right)}$ is the slip resistance, ${\dot{\gamma }}^{\left(\alpha \right)}$ is the slip rate on the slip system $\left(\alpha \right)$, ${\dot{\gamma }}_{0}$ is the reference shear rate, and $m$ is the strain-rate sensitivity exponent. The slip resistance can be evolved by using the rate form:

$${\dot{g}}^{(\alpha )}=\sum_{\beta =1}^{NSYS}{h}^{\alpha \beta }\left|{\dot{\gamma }}^{(\beta )}\right|,$$

(40)

where ${h}^{\alpha \beta }$ represents the hardening matrix. The hardening matrix can be defined as:

$${h}^{\alpha\beta}={q}^{\alpha\beta }{h}^{\beta }\left(\mathrm{no\;sum\;on}\;\beta\right),$$

(41)

with

$${q}^{\alpha \beta }=\left[\begin{array}{cccc}\mathbf{A}& q\mathbf{A}& q\mathbf{A}& q\mathbf{A}\\ q\mathbf{A}& \mathbf{A}& q\mathbf{A}& q\mathbf{A}\\ q\mathbf{A}& q\mathbf{A}& \mathbf{A}& q\mathbf{A}\\ q\mathbf{A}& q\mathbf{A}& q\mathbf{A}& \mathbf{A}\end{array}\right],$$

(42)

where $\mathbf{A}$ is the $3\times 3$ unity matrix, and $1.0\le q\le 1.4$ [48]. In this section, the isotropic hardening model was assumed by setting $q=1$. Finally, the slip hardening rate becomes as follows [49, 50]:

$${h}^{\beta }={h}_{s}+({h}_{0}-{h}_{s}){{\text{sech}}}^{2}\left\{\left(\frac{{h}_{0}-{h}_{s}}{{g}_{s}-{g}_{0}}\right)\gamma \right\},\text{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}}\gamma ={\int }_{0}^{t}\sum_{\alpha =1}^{NSYS}\left|{\dot{\gamma }}^{(\alpha )}\right|dt,$$

(43)

where ${h}_{s}$,${h}_{0}$ are the asymptotic and initial hardening rates, and ${g}_{s}$,${g}_{0}$ are the saturated and initial values of the slip resistance. The material constants for CPFEM are shown in Table 1. Boundary conditions for the uniaxial tensile and simple shear simulations are illustrated in Fig. 4. Also, 1 [mm]$\times$ 1 [mm]$\times$ 1 [mm] solid element with reduced integration (C3D8R element in Abaqus) was employed.

Table 1 Material constants for comparison used in crystal plasticity finite element analysis

Fig. 4

Boundary conditions for single element simulations: a Uniaxial tension; b Simple shear

Various time increment sizes were imposed for the simulations where the total time increment was 1 s in order to evaluate the accuracy of various stress integration methods. Figures 5 and 6 show the stress-strain curves from the uniaxial tensile and simple shear simulations, respectively. In each simulation, the time increment was set to constant. In other words, 0.001 s constant time for ‘1000 steps’, 0.1 s constant time for ‘10 steps’, and 0.25 s constant time for ‘4 steps’ were imposed. The uniaxial tensile simulation results from the FDM-based model show good agreement with those from the Euler backward method using the analytical derivative and those from the semi-implicit method even for the large time increment size. On the other hand, a large error occurs when the time increment size becomes larger in the simple shear simulation results from the semi-implicit method. Since the quantities of kinematics are not from the current configuration but from the last convergence state, the semi-implicit method is converged with a large error or diverged when the large change in the elastic deformation gradient exists. For the Euler backward method using the analytical derivative and FDM-based model, they showed good agreement with each other regardless of time increment sizes. Therefore, the Euler backward method using the analytical derivative and FDM-based model are beneficial when a large time increment and rotation are imposed. In addition, the average relative errors of first derivatives of resolved shear stresses were evaluated between the FDM and analytical derivative. The errors were compared at the first Newton-Raphson iteration at the half-time step (0.5 s) of the simulation where total four steps were imposed. The estimated errors are listed in Table 2, and the errors are negligible.

Fig. 5

Comparisons of stress-strain curves from uniaxial tensile simulations according to the number of steps and various stress integration methods: a various methods with thousand increments; b various steps under semi-implicit method; c various steps under Euler backward method using analytical derivative; d various steps under the proposed method using central difference method; e various steps under the proposed method using backward difference method; b various steps under the proposed method using forward difference method

Fig. 6

Comparisons of stress-strain curves from simple shear simulations according to the number of steps and various stress integration methods: a various methods with thousand increments; b various steps under semi-implicit method; c various steps under Euler backward method using analytical derivative; d various steps under the proposed method using central difference method; e various steps under the proposed method using backward difference method; b various steps under the proposed method using forward difference method

Table 2 Relative errors [%] of three finite difference methods with respect to analytical derivatives

Convergence

The uniaxial tensile simulation with the same boundary condition in Fig. 4a was performed to evaluate the convergence of various stress integration methods. One hundred random orientations were used, and the (1 1 1) pole figure is shown in Fig. 7a. The material constants for CPFEM were the same with Table 1 except for the strain-rate sensitivity exponent ($m$) which was set to 0.01 and 0.001 in this section. The results of $m=0.01$ from various stress integration methods are shown in Fig. 8. The time increment was 0.001 s for ‘1000 steps’. For the ‘auto’ scheme, the time increment was bounded between ${10}^{-6}$ s and $1$ s with an initial time increment of ${10}^{-3}$ s. With small time increments, every stress integration method predicts the same result. However, the number of steps and predicted results become different when the time increment increases. For the semi-implicit method, 27 steps were required to be converged, and significant errors also occurred. It is because the elastic deformation gradients are from the last convergence state, and the kinematics is assumed to evolve linearly over a time increment. It resulted in the deviated solution when the time increment size was large. For the Euler backward method using the analytical derivative and FDM-based model, only 25 steps were used for convergence, and the accuracy was also maintained. The convergence issue was more pronounced at lower strain-rate sensitivity ($m=0.001$). Although the semi-implicit method failed to converge, the Euler backward method using the analytical derivative and FDM-based model predicted the same results and converged with 91 steps. In conclusion, the proposed FDM-based model has the same accuracy and convergence with the Euler backward method using analytical derivative.

Fig. 7

(1 1 1) pole figures of random orientations: a 100 grains; b 300 grains

Fig. 8

Comparisons of stress-strain curves from uniaxial tensile simulations according to the number of steps and various stress integration methods: a various methods with thousand increments; b under semi-implicit method; c under Euler backward method using analytical derivative; d under the proposed method using central difference method; e under the proposed method using backward difference method; b under the proposed method using forward difference method

Computational efficiency

The uniaxial tensile simulation was performed, as shown in Fig. 4a to investigate computational efficiency. In order to evaluate the computational time according to the derivatives more clearly, three hundred random orientations were employed, and the (1 1 1) pole figure for the random grains is shown in Fig. 7b. The material constants for CPFEM were the same with Table 1 except for the strain-rate sensitivity exponent ($m$) which was set to 0.01. All the simulations took a constant time increment, and thousand steps were imposed in total. Since the simulation time depends on the terms for computing the tensor exponential functions and their derivatives, as shown in Eqs. (27) and (29), the computation time according to the terms was compared. The wall-clock simulation times from all the cases for terms and stress integration methods were averaged after performing simulations three times. Also, the averaged simulation time was normalized by that from the Euler backward method using analytical derivative with ${n}_{{\text{max}}}=2$ in Eqs. (27) and (29). As shown in Fig. 9, the simulation time from the FDM-based model is faster than that from the Euler backward method using the analytical derivative in most cases. Besides, the Euler backward method using analytical derivative highly depends on how many terms are used for computing the tensor exponential functions and their derivatives, while the number of terms has a negligible effect on the FDM-based model. Therefore, the FDM approach is more efficient than the analytical derivatives when the number of terms for computing the series expansion is large.

Fig. 9

Normalized wall-clock simulation time depending on ways of derivatives and terms for computing tensor exponential functions

Validation with the reduced texture approach

In order to validate the proposed FDM-based model, the reduced texture approach was implemented due to the benefit of computational efficiency in large-scale CPFEM. In the reduced texture approach, the specific parameters including crystallographic orientations are calibrated to characterize the plastic anisotropy of matal. The main objective of the reduced texture approach is not a complete representation of the real texture behavior but is to use the crystal plasticity model in large-scale engineering applications with a reasonable computational cost by minimizing the number of crystallographic orientations. The reduced texture methodology was proposed with the self-consistent model [27,28,29]. In this chapter, the Euler backward method using analytical derivative and the proposed FDM-based model were applied.

Experimental data

For experimental data, a strongly anisotropic aluminum alloy (AA 2090-T3) data was used. The r-values and yield stress ratios (${s}_{\theta }={\sigma }_{\theta }/{\sigma }_{0}$) where ${\sigma }_{\theta }$ indicates the yield stress along the $\theta$ degrees from the rolling direction were obtained from the tensile test carried out along the rolling (RD), diagonal (DD), and transverse (TD) directions in the previous research, and experimental results are specified in Table 3 [51, 52]. The reference data were generated by considering the experimental results; the r-values, yield stress ratios, and stress-strain hardening curve along the RD [27, 28]. Therefore, the stress-strain hardening curves along the $\theta$ degrees from the RD are written in Eq. (44) with the assumption of isotropic hardening.

Table 3 Experimental results of AA 2090-T3

$${\tilde{\sigma }}_{\theta }=646{\left(0.025+{\overline{\varepsilon }}^{p}\right)}^{0.227}\cdot {s}_{\theta }\left[{\text{MPa}}\right].$$

(44)

The transverse strain curves were also built with the assumption that the strain value has a very small effect on the r-values. All of the data were reconstructed up to ${\varepsilon }_{11}=0.28$ to consider the prediction of the large strain value in the circular cup drawing simulation. Figure 10 shows the reconstructed experimental data.

Fig. 10

Reconstructed experimental data: a stress-strain curves; b transverse strain curves

Calibration of parameters

Calibration of parameters was conducted with Isight2016 and Abaqus2016/Standard. As mentioned in “Experimental data” section, the experimental results prepared for the calibration are defined as:

• True stress-strain curves(RD, DD, TD; three curves)

• Transverse strain curves(RD, DD, TD; three curves)

For the material characterization, the single element tensile simulations were performed along the RD, DD, and TD. 1 [mm]$\times$ 1 [mm]$\times$ 1 [mm] solid element with reduced integration (C3D8R element in Abaqus) was elongated to produce the strain (${\varepsilon }_{11}$) value of 0.28. Total time increment was set to be 2,800 s to impose ${\dot{\varepsilon }}_{11}={10}^{-4}$ ${s}^{-1}$. The fixed material constants are shown in Table 4. Among the parameters to be calibrated, there are specific parameters that imply the information of orientations. In the same way as Rousselier et al. [27], two texture representatives were adopted; $\left({\varphi }_{1}\left(1\right)\phi \left(1\right){\varphi }_{2}\left(1\right)\right)$ and $\left({\varphi }_{1}\left(2\right)\phi \left(2\right){\varphi }_{2}\left(2\right)\right)$. Each texture representative consists of four crystallographic orientations defined as $\left({\varphi }_{1}\phi {\varphi }_{2}\right)$,$\left(-{\varphi }_{1}\phi -{\varphi }_{2}\right)$,$\left(-{\varphi }_{1}-\phi -{\varphi }_{2}\right)$, and $\left({\varphi }_{1}-\phi {\varphi }_{2}\right)$ for the orthotropic symmetry. Total eight crystallographic orientations are employed at each integration point. In addition, other parameters related to the slip resistance and latent hardening exist, thus all of the parameters to be optimized are as follows:

• Latent hardening parameter ($q$; one parameter),

• Slip resistance parameters (${h}_{0}$,${h}_{s}$,${g}_{0}$, and ${g}_{s}$; four parameters),

• Orientation parameters

  $$\left(\left({\varphi }_{1}\left(1\right)\phi \left(1\right){\varphi }_{2}\left(1\right)\right), \left({\varphi }_{1}\left(2\right)\phi \left(2\right){\varphi }_{2}\left(2\right)\right),\mathrm{\;and\;}f\left(1\right);\mathrm{seven\;parameters}\right),$$

where $f\left(1\right)$ represents the volume fraction of grains for texture representative (1), and $q$ is limited between 1 and 1.4. Even though it is challenging to identify the parameter of latent hardening from the only monotonic tensile experiment where strain path is not changed, it was adopted for the parameter as a fitting scheme.

Table 4 Material constants for characterizing AA 2090-T3

Table 5 shows the calibrated parameters, and the corresponding simulation results from the integration methods using analytical derivative and FDM are shown in Fig. 11. The integration methods using analytical derivative and FDM predicted the same results. In addition, the predicted stress-strain and transverse strain curves show good agreement with the experimental data. With eight crystallographic orientations, the FDM-based model successfully predicted the anisotropic behavior of AA 2090-T3.

Table 5 Results of calibrated parameters for AA 2090-T3

Fig. 11

Experimental data and simulation results from: a stress-strain curves under Euler backward method using analytical derivative; b transverse strain curves under Euler backward method using analytical derivative; c stress-strain curves under the proposed FDM-based model; d transverse strain curves under the proposed FDM-based model

Cup drawing simulation

A circular cup drawing simulation was performed to estimate the applicability of FDM-based model for large-scale engineering problems. Figure 12a shows the schematic view, tool dimensions, and process variables for the circular cup drawing analysis [52]. In the circular cup drawing analysis, a quarter part of the blank was analyzed due to the assumption that the metal has orthotropic symmetry, and the mesh of the blank is shown in Fig. 12b. For the elements, 100 and 2400 elements were utilized for C3D6 and C3D8R, respectively. Total time was 5000 s, and the time increment was allowed between 0.001 s and 30 s.

Fig. 12

a Schematic view, tool dimensions, and process variables for circular cup drawing; b Quarter of blank mesh for circular cup drawing simulation

The cup drawing simulations were conducted with the Euler backward method using analytical derivative and FDM; the central, backward, and forward difference methods. Therefore, four cases were performed and compared in terms of accuracy and computational efficiency. Figure 13a shows the earing profiles from the experimental results and predictions from the Yld2000-2d [38] and the reduced texture approach with the integration method using the analytical derivative. Figure 13b shows the earing profiles predicted by the FDM-based model. The numbers of simulation steps and normalized times from various methods are specified in Table 6. Three comparisons are following. First, all of the predicted cup heights from the reduced texture approach are the same regardless of the ways of derivatives; the analytical derivative, central, backward, and forward difference methods. Second, the predicted earing profiles show good agreement with the experimental result. The cup height from the reduced texture approach was also compared with that from the Yld2000-2d [38] since both of them employed the yield stress ratios and r-values along the RD, DD, and TD. For the Yld2000-2d, a small ear in the 15 degrees from the RD was not predicted. On the other hand, the reduced texture approach could predict six ears including the small ear in the 15 degrees, in spite of the fact that the experimental r-value along the 75 degrees was not utilized. As shown in Fig. 14, r-value at 75 degrees affects the small ear in 15 degrees, as explained by Yoon et al. [52]. It was also shown that r-value has a large effect on the cup height in other research [53, 54]. In Fig. 14, the r-values were obtained every 15 degrees in the tensile strain range between 0.03 and 0.2, and the spline interpolation was used to draw the line. Also, the FDM-based model using the central difference method was employed. Third, the numbers of simulation steps and normalized times (wall-clock time) of all the cases were compared, which were averaged after performing three times per case. As shown in Table 6, the numbers of simulation steps are the same regardless of the ways of derivatives. In the view of computational efficiency, the FDM-based model had a shorter simulation time than the integration method using analytical derivative. Considering computational efficiency, the backward and forward difference methods could be a better choice since their accuracy is also quite high.

Fig. 13

Cup height from: a experiment, Yld2000-2d, and reduced texture approach with Euler backward method using analytical derivatives; b reduced texture approach with the proposed FDM-based model

Table 6 Results of the number of steps and normalized wall-clock simulation time

Fig. 14

Simulation result of r-value predicted by the FDM-based model using central difference method

Summary and conclusion

The Euler backward stress integration algorithm using the FDM was proposed for CPFEM. The FDM successfully approximated the analytical derivative of the resolved shear stress. Three methods were investigated for the FDM: the central, backward, and forward difference methods. Three finite difference methods were adopted because the central difference method showed high accuracy because of the higher-order truncation error. On the other hand, the backward and forward difference methods showed high computational efficiency.

The proposed FDM-based model was verified through the single element simulations in terms of accuracy, convergence, and computational efficiency. The accuracy and convergence of the FDM-based model were almost the same with the Euler backward method using analytical derivative. In the view of computational efficiency, simulation time depends on how many terms are used for computing the tensor exponential functions and their derivatives. The results indicate that the larger terms are employed, the more efficient the FDM-based model is. Moreover, the FDM-based model can be easily applied to any form of constitutive model for CPFEM (such as Eq. (10)), so that the restriction to derive the complex derivative could be solved regardless of the constitutive model. Taking into account the single element simulation results, the proposed FDM-based model could be the alternative to the Euler backward method using analytical derivative in the view of accuracy, convergence, computational efficiency, and easy implementation.

For further validation of the FDM-based model, the reduced texture approach was adopted to implement the crystal plasticity model into a large-scale metal forming example. In this work, the anisotropic behavior of AA 2090-T3 was characterized. Circular cup drawing simulations were also performed, and the earing profiles were compared. The same earing profiles were predicted regardless of the ways of derivatives. However, the FDM-based model showed a shorter simulation time than the Euler backward method using analytical derivative. The proposed FDM-based model was successfully implemented in CPFEM and validated through the prediction of plastic anisotropy and earing profile.

Appendices

Appendix 1. Detailed stress integration algorithm in fully-implicit method

In order to use the Newton-Raphson method, the resolved shear stress derivative is required for the Jacobian matrix. Differentiating Eq. (17) with respect to the slip rate gives:

$$\frac{\partial {\tau }_{n+1}^{(\alpha )}}{\partial {\dot{\gamma }}_{n+1}^{(\beta )}}={\mathbf{b}}_{0}^{(\alpha )}\cdot \frac{\partial {{\varvec{\Sigma}}}_{n+1}}{\partial {\dot{\gamma }}_{n+1}^{(\beta )}}\cdot {\mathbf{n}}_{0}^{(\alpha )}.$$

(45)

This equation is followed by the derivative of Eq. (18):

$$\frac{\partial {{\varvec{\Sigma}}}_{n+1}}{\partial {\dot{\gamma }}_{n+1}^{(\beta )}}=\frac{\partial {\overline{\mathbf{C}} }_{n+1}^{e}}{\partial {\dot{\gamma }}_{n+1}^{(\beta )}}\cdot {\overline{\mathbf{S}} }_{n+1}+{\overline{\mathbf{C}} }_{n+1}^{e}\cdot \frac{\partial {\overline{\mathbf{S}} }_{n+1}}{\partial {\dot{\gamma }}_{n+1}^{(\beta )}}.$$

(46)

Calculations of the derivatives of Eqs. (8)–(9) yield

$$\frac{{\partial \overline{\mathbf{S}} }_{n+1}}{\partial {\dot{\gamma }}_{n+1}^{(\beta )}}={\mathcal{L}}^{{\varvec{e}}}:\frac{\partial {\overline{\mathbf{E}} }_{n+1}^{e}}{\partial {\dot{\gamma }}_{n+1}^{(\beta )}},$$

(47)

$$\frac{\partial {\overline{\mathbf{E}} }_{n+1}^{e}}{\partial {\dot{\gamma }}_{n+1}^{(\beta )}}=\frac{1}{2}\frac{\partial {\overline{\mathbf{C}} }_{n+1}^{e}}{\partial {\dot{\gamma }}_{n+1}^{(\beta )}},$$

(48)

$$\frac{\partial {\overline{\mathbf{C}} }_{n+1}^{e}}{\partial {\dot{\gamma }}_{n+1}^{(\beta )}}=\frac{\partial \left\{{\left({\mathbf{F}}_{n+1}^{e}\right)}^{T}\right\}}{\partial {\dot{\gamma }}_{n+1}^{(\beta )}}\cdot {\mathbf{F}}_{n+1}^{e}+{\left({\mathbf{F}}_{n+1}^{e}\right)}^{T}\cdot \frac{\partial {\mathbf{F}}_{n+1}^{e}}{\partial {\dot{\gamma }}_{n+1}^{(\beta )}}.$$

(49)

Differentiation of Eq. (32) gives

$$\frac{\partial {\mathbf{F}}_{n+1}^{e}}{\partial {\dot{\gamma }}_{n+1}^{(\beta )}}={\mathbf{F}}_{n+1}\cdot {({\mathbf{F}}_{n})}^{-1}\cdot {\mathbf{F}}_{n}^{e}\cdot \frac{\partial {\text{exp}}\left[\sum\limits_{\alpha =1}^{NSYS}\left(-\Delta t\cdot {\dot{\gamma }}_{n+1}^{(\alpha )}\cdot {\mathbf{b}}_{0}^{(\alpha )}\otimes {\mathbf{n}}_{0}^{(\alpha )}\right)\right]}{\partial \left[\sum\limits_{\alpha =1}^{NSYS}\left(-\Delta t\cdot {\dot{\gamma }}_{n+1}^{(\alpha )}\cdot {\mathbf{b}}_{0}^{(\alpha )}\otimes {\mathbf{n}}_{0}^{(\alpha )}\right)\right]}:\left(-\Delta t\cdot {\mathbf{b}}_{0}^{(\beta )}\otimes {\mathbf{n}}_{0}^{(\beta )}\right),$$

(50)

where the derivative of tensor exponential function is specified in Eq. (29).

The unknown quantities are obtained by using the Newton-Raphson method in a similar way to the semi-implicit method. Kirchhoff stress at the current configuration can be updated by using Eqs. (8)–(9) and (10). The grain level stress at the current configuration is updated as:

$${{\varvec{\uptau}}}_{n+1}={\mathbf{F}}_{n+1}^{e}\cdot {\overline{\mathbf{S}} }_{n+1}\cdot {\left({\mathbf{F}}_{n+1}^{e}\right)}^{T}.$$

(51)

In the Euler backward method, the stress state can be computed directly, as shown in Eq. (51), so there is no need for the space to store the stress state. On the other hand, the stress state is updated by adding Jaumann rate from the previous stress state in the semi-implicit method, thus the storage space for the stress state is required. Therefore, the Euler backward method requires less storage space than the semi-implicit method. Further calculation for the tangent modulus can be found in Appendix 2.

Appendix 2. Tangent modulus

The tangent modulus is required to ensure the quadratic convergence rate in CPFEM. The derivative of the consistency condition takes the following form:

$${\dot{\phi }}^{(\alpha )}={\dot{\tau }}^{(\alpha )}-{\dot{\tau }}_{rss}^{(\alpha )}.$$

(52)

The critical resolved shear stress (${\tau }_{rss}^{(\alpha )}$) is the function of the slip rate (${\dot{\gamma }}^{(\alpha )}$) and the slip resistance (${g}^{(\alpha )}$). Therefore, the derivative of critical resolved shear stress becomes

$${\dot{\tau }}_{rss}^{(\alpha )}=\sum_{\beta =1}^{NSYS}\frac{\partial {\tau }_{rss}^{(\alpha )}}{\partial {\dot{\gamma }}^{(\beta )}}\frac{{\dot{\gamma }}^{(\beta )}}{\Delta t}+\sum_{\beta =1}^{NSYS}\sum_{\delta =1}^{NSYS}\frac{\partial {\tau }_{rss}^{(\alpha )}}{\partial {g}^{(\delta )}}\frac{\partial {g}^{(\delta )}}{\partial {\dot{\gamma }}^{(\beta )}}\frac{{\dot{\gamma }}^{(\beta )}}{\Delta t}.$$

(53)

Substituting Eqs. (19) and (53) into Eq. (52) gives

$${\mathbf{R}}^{(\alpha )}:\left(\mathbf{D}\Delta t\right)=\sum_{\beta =1}^{NSYS}{\chi }^{(\alpha \beta )}\cdot {\dot{\gamma }}^{(\beta )},$$

(54)

with

$${\chi }^{(\alpha \beta )}=\frac{\partial {\tau }_{rss}^{(\alpha )}}{\partial {\dot{\gamma }}^{(\beta )}}+\sum_{\delta =1}^{NSYS}\frac{\partial {\tau }_{rss}^{(\alpha )}}{\partial {g}^{(\delta )}}\frac{\partial {g}^{(\delta )}}{\partial {\dot{\gamma }}^{(\beta )}}+\Delta t\cdot {\mathbf{R}}^{(\alpha )}:{\mathbf{P}}^{(\beta )}.$$

(55)

Taking the inverse matrix of Eq. (55) ($\left[\psi \right]={\left[\chi \right]}^{-1}$) on both sides of Eq. (54) gives the expression of the slip rate in terms of deformation rate as follows:

$${\dot{\gamma }}^{(\alpha )}=\sum_{\beta =1}^{NSYS}{\psi }^{(\alpha \beta )}{\mathbf{R}}^{(\beta )}:(\mathbf{D}\Delta t).$$

(56)

Finally, the tangent modulus can be defined from Eq. (16):

(57)