Multi-level deep drawing simulations of AA3104 aluminium alloy using crystal plasticity finite element modelling and phenomenological yield function

Abstract

This paper proposed a hierarchical multi-level model to study the crystallographic texture induced mechanical anisotropy of AA3104-H19 aluminium sheet from mesoscale to continuum scale. In the mesoscale, full-field crystal plasticity finite element method (CPFEM) was used to provide both in-plane and out-of-plane yield stresses and plastic potential points in various deformation modes. In the continuum scale, these materials sampling points were used to determine the parameters of two phenomenological yield functions (Yld2000-2d in plane stress space and Yld2004-18p in 3D stress space) using associated flow rule (AFR) and non-associated flow rule (non-AFR). The results indicate that higher accuracy obtained by Yld2000-2d and Yld2004-18p yield functions associated with non-AFR in comparison with AFR. These phenomenological models were successfully implemented into finite element (FE) code using an explicit integration scheme to simulate sheet metal forming. It is found that the 3D Yld2004-18p model involved with both in-plane and out-of-plane anisotropies is superior to 2D Yld2000-2d model which only accounts for in-plane anisotropy.

Introduction

Aluminium alloy sheets generally exhibit a considerable mechanical anisotropy due to the strong crystallographic texture induced by heavy cold rolling reductions. This property has a significant influence on the sheet metal forming operations, for example, the final shape and spring-back level of the formed part, reaction force and wear of tooling system. Therefore, the investigation of mechanical anisotropy is important for sheet metal forming.

In conventional phenomenological approach, the yield surface plays a pivotal role in describing the plastic mechanical anisotropy. Mathematically, the shape of yield surface is determined by a convex function using a set of parameters which obtained by fitting the yield stress and plastic potential points in various stress states provided by well-designed physical tests. ¹ The phenomenological approach shows strong robustness for the simulation of sheet metal forming, and high computational efficiency on plane-stress problems. However, it is hard to characterise the out-of-plane anisotropy of sheet metals via physical tests. Thus, early yield function is only used for plane-stress problems, for example, the widely used Yld2000-2d plane stress yield function.^2,3 To obtain the mechanical anisotropy in 3D stress case, Barlat et al. proposed a convex yield function named Yld2004-3D ⁴ in which the out-of-plane mechanical anisotropy is determined by an analytical crystal plasticity method so-called TBH. ⁵ This yield function has been successfully implemented into FE code and shown an accurate earing prediction of strong textured AA2090-T3 in deep drawing process. ⁶ However, the assumption of the analytical TBH crystal plasticity (CP) model cannot satisfy the stress equilibrium and compatibility simultaneously. ⁷ The TBH model showed great deviation in the prediction of directional yield stresses and plastic potentials against experimental tests. ⁸

In classic flow theory of plasticity, the associated flow rule (AFR) is generally employed in the yield function to describe the yielding and plastic potential for strain increment simultaneously. However, the plastic potentials are difficult to associate with the yield surface when describing a highly anisotropic material with an identical function.^9–11 Therefore, the non-associated flow rule (non-AFR) can be applied to overcome this limitation which eliminates the constraint of the correlation between yield surface and plastic potential.^12–14 In this flow theory, two independent functions are used to determine the yielding and plastic potential, respectively. Various combinations of phenomenological yield functions associated with non-AFR were studied by previous works. For instance, the yld2000-2d ² and Hill48 ¹⁵ yield functions using non-AFR were superior to the conventional AFR with respect to the fitting of mechanical anisotropy for AA2090-T3, AA5042, AA5754-O and AA5182 sheets.^10,11,16,17 The non-AFR also showed an adequate capability in description of anisotropic and kinematic hardening behaviour of various strong textured sheet metals.^11,16–20 However, these previous works were based on the phenomenological method using plane stress assumption which were relied on the physical tensile experiments and only determines the in-plane mechanical anisotropy.

Recently, the full-field CPFEM was developed to simulate the polycrystalline aggregates under various loading conditions which not only accounted for the crystallographic texture but also incorporated with additional morphological information, for example, the grain size, boundary and orientation. Therefore, the local grain interaction and intra-grain inhomogeneities are simulated in the crystal plasticity modelling.^7,21 Compared with analytical CP model, the main advantages of full-field CPFEM were reflected in their robust solution to solve polycrystal micromechanics problems and convergence under complicated boundary conditions. Moreover, CPFEM considered inter-grain and intra-grain micromechanical interactions in the physics of crystal mechanics. ²² However, the main drawback of CPFEM was associated with high computational costs when fine meshes are assigned to construct a high resolution representative volume element (RVE).^21,23,24

Since the full-field CPFEM showed a good capability in the description of micromechanics of polycrystal under external loads, this advanced method can successfully reproduce the crystallographic texture induced mechanical anisotropy of as-rolled sheet metal observed by experimental uniaxial tensile tests in various material directions.^7,8,25–29 Furthermore, CPFEM can be used to predict the out-of-plane anisotropy due to no constraint on deformation mode. This property was difficult to be obtained by physical tests due to the thin wall structure of sheet metal. These advantages enabled a large number of in-plane and out-of-plane yield stresses and plastic potentials points to be predictable.

To comprehensively study the plastic behaviours of strong textured aluminium alloy sheet, a hierarchically structured multi-scale constitutive model (from mesoscopic crystallographic scale to macro continuum scale) was developed to accurately describe the in-plane and out-of-plane yielding and plastic potentials under multi-axial loading conditions in this study. In this method, an RVE of the strong textured AA3104-H19 aluminium alloy was constructed to predict the mechanical anisotropy of polycrystalline aggregates under various in-plane and out-of-plane loading conditions using full-field CPFEM. These data were employed to determine the parameters of the phenomenological yield and plastic potential functions using non-AFR. This hierarchically structured model (please see Figure 1) has been successfully implemented into FE code and can be readily adapted for forming simulations of different aluminium sheets.

Figure 1.

The schematic configuration of the multiscale and multiphysics CP modelling framework used for sheet metal forming simulations.

Crystal plasticity modelling

Constitutive formulations

The crystal kinematics of isothermal finite deformation describes the process $\mathbf{x}\in {\mathcal{B}}_0\to \mathbf{y}\in \mathcal{B}$ which maps material points $\mathbf{x}$ in a reference configuration ${\mathcal{B}}_0$ , to points $\mathbf{y}$ in the current configuration $\mathcal{B}$ , by externally applied forces and displacements over a period of time. To compute this deformation, the positions of the neighbouring material points are related to an arbitrary origin in the reference configuration by the vector $\mathrm{d}\mathbf{x}$ . As a result of deformation, this vector is mapped into its state in the current configuration, $\mathrm{d}\mathbf{y}=\mathrm{d}\mathbf{x}+\mathrm{d}\mathbf{u}$ , where $\mathrm{d}\mathbf{u}$ is the differential total displacement vector. These vectors are related by the total deformation gradient, ²¹

$$\mathrm{d}\mathbf{y}=\left(\frac{\partial \mathbf{y}}{\partial \mathbf{x}}\right)\mathrm{d}\mathbf{x}=(\mathbf{I}+\frac{\partial \mathbf{u}}{\partial \mathbf{x}})\mathrm{d}\mathbf{x}=\mathbf{F}\mathrm{d}\mathbf{x}$$

(1)

where $\mathbf{I}$ is the second rank identity tensor. The nomenclature used in this study is summarized in Table 1. The deformation resulted gradient defined as $\mathbf{F}=\mathrm{d}\mathbf{y}/\mathrm{d}\mathbf{x}$ then can be written as an elastic ${\mathbf{F}}_e$ and a plastic ${\mathbf{F}}_p$ component based on the local multiplicative decomposition ³⁰ :

$$\mathbf{F}={\mathbf{F}}_e{\mathbf{F}}_p={\mathbf{V}}_e{\mathbf{R}}_e{\mathbf{F}}_p$$

(2)

where the elastic deformation ${\mathbf{F}}^e={\mathbf{V}}^e{\mathbf{R}}^e$ describes the elastic stretching and the rotation of the lattice due to external applied forces and displacements, while ${\mathbf{F}}_p$ is an irreversible permanent deformation that describes the motion of dislocations (plastic slip) on crystallographic planes (slip planes) leaving the crystal lattice unchanged. The decomposition (equation (2)) introduces two intermediate configurations between the undeformed ${\mathcal{B}}_0$ and the current configurations $\mathcal{B}$ which are defined as $\overline{\mathcal{B}}$ and $\tilde{\mathcal{B}}$ (please see Figure 2).

Table 1.

The table of nomenclature used in this study.

Variable	Definition	Units
$\mathbf{I}$	Second rank identity tensor
$\mathbf{F}$	Deformation gradient
${\mathbf{F}}^e$	Elastic deformation gradient
${\mathbf{F}}_p$	Plastic deformation gradient
$\mathbf{L}$	Spatial velocity gradient
${\mathbf{L}}_e$	Elastic velocity gradient
${\mathbf{L}}_p$	Plastic velocity gradient
${\stackrel{·}{\gamma }}^{\alpha }$	Shear strain rate	s−1
${\stackrel{·}{\gamma }}^0$	Reference shear strain rate	s−1
$m$	rate sensitivity exponent
$g^{\alpha }$	Slip resistance	MPa
${\tau }^{\alpha }$	Applied shear stress	MPa
${\tau }_0$	Initial hardness	MPa
$h_{\alpha \beta }$	Instantaneous strain hardening matrix	MPa
$q$	Latent hardening parameter
$h_0$	Initial hardening modulus	MPa
$a$	Estimated constant
$g_{\infty }$	Saturated slip resistance	MPa
$\mathbf{\sigma }$	Cauchy stress tensor	MPa
$\mathbf{S}$	Piola–Kirchhoff stress tensor	MPa
${\mathbf{E}}_e$	Elastic Green strain tensor
$\mathbb{C}$	Elasticity matrix	GPa
$\overline{\sigma }$	Equivalent stress	MPa
$w$	Yield function exponent
$\phi$	Unconditionally convex function
$X^{(k)}$	Yld2000-2d principal values
${\tilde{\mathbf{x}}}^{(k)}$	Yld2000-2d transformed stress tensor
$\mathbf{L}$	Yld2000-2d linear transformations
${\tilde{S}}^{(k)}$	Yld2004-3d principal values
${\tilde{\mathbf{s}}}^{(k)}$	Yld2004-3d transformed stress tensor
${\mathbf{C}}^{(k)}$	Yld2004-3d linear transformations
$\mathbf{T}$	Constant linear transformation
$\mathrm{\Delta }{\epsilon }^p$	Plastic strain increment
$\mathrm{\Delta }{\overline{\epsilon }}^p$	Equivalent plastic strain increment
$\stackrel{·}{\lambda }$	Plastic strain-rate multiplier
$\mathrm{\Delta }{\sigma }_{\mathrm{ij}}^b$	Elastic predictor stress	MPa
$\mathrm{\Delta }{\epsilon }_{\mathrm{ij}}$	Strain tensor increment
$\nu$	Poisson’s ratio
${\mathbf{C}}_{\mathrm{ijkl}}$	Elastic stiffness tensor	GPa
$\mathrm{\Delta }Y$	Hardening modules increment	GPa
$\mathrm{\Delta }\lambda$	Consistent modules	GPa

Figure 2.

The schematics of lattice deformation for two intermediate configurations between the undeformed and the current configurations governed by elastic stretching, rotation and the motion of dislocation.

The crystal kinematics of finite deformations requires an constitutive equation for the time rate of change of deformation gradient $\mathbf{F}$ . The velocity of each material point of a body in motion forms a vector field measured in the current configuration. Using the kinematic decomposition (equation (2)), we can write the spatial velocity gradient $\mathbf{L}$ in the current configuration $\mathcal{B}$ given as:

$$\mathbf{L}=\stackrel{·}{\mathbf{F}}{\mathbf{F}}^{-1}={\mathbf{L}}_e+{\mathbf{F}}_e{\mathbf{L}}_p({\mathbf{F}}_e{)}^{-1}$$

(3)

where ${\mathbf{L}}_p$ is the plastic velocity gradient in $\mathcal{B}$ . The ${\mathbf{L}}_p$ only describes the motion of dislocations on crystallographic planes where dislocation slip is the only deformation process. Therefore, ${\mathbf{L}}_p$ can be defined as sum of the shear rates on all slip systems:

$${\mathbf{L}}_p={\stackrel{·}{\mathbf{F}}}_p{\mathbf{F}}_p^{-1}=\sum _{\alpha }^N{\stackrel{·}{\gamma }}^{\alpha }{\mathbf{s}}_0^{\alpha }\otimes {\mathbf{m}}_0^{\alpha }$$

(4)

where ${\stackrel{·}{\gamma }}^{\alpha }$ is the slip rate of the slip system $\alpha$ and $N$ is the number of slip systems. The unit vectors ( ${\mathbf{s}}_0^{\alpha }$ , ${\mathbf{m}}_0^{\alpha }$ ) define the orientation of the $\alpha$ -slip system in ${\mathbf{\beta }}_0$ through the dyadic product ${\mathbf{S}}_0^{\alpha }={\mathbf{s}}_0^{\alpha }\otimes {\mathbf{m}}_0^{\alpha }$ , the Schmid tensor. Kinetics of a slip system (flow rule) assume that the resolved shear stress on the $\alpha$ -slip system can be decomposed as ${\tau }^{\alpha }$ which is the applied shear stress needed to overcome the barriers to dislocation motion. The kinetics of slip is based on the well-known power law^31,32 and given as:

$${\stackrel{·}{\gamma }}^{\alpha }={\stackrel{·}{\gamma }}^0{\left|\frac{{\tau }_a}{g^{\alpha }}\right|}^{1/m}\mathrm{sign}({\tau }_a)$$

(5)

where ${\stackrel{·}{\gamma }}^0$ is a reference shear strain rate and exponent $m$ determines the rate sensitivity of slip. $g^{\alpha }$ is the slip resistance which characterises the strain hardened state of crystal in the current configuration. ³³ The evolution of $g^{\alpha }$ is governed by:

$${\stackrel{·}{g}}^{\alpha }=\sum _{\beta }^N{h}_{\alpha \beta }{\stackrel{·}{\gamma }}^{\beta }$$

(6)

where $g^{\alpha }(0)={\tau }_0$ is the initial hardness for slip system of each grain, and $h_{\alpha \beta }$ is the instantaneous strain hardening matrix which given by:

$$h_{\alpha \beta }=h_0[q+(1-q){\delta }_{\alpha \beta }]{|1-\frac{g^{\beta }}{g_{\infty }}|}^a\mathrm{sgn}(1-\frac{g^{\beta }}{g_{\infty }})$$

(7)

where $q$ is a latent hardening parameter. For a coplanar slip system $q=1$ , and for a non-coplanar slip system $q=1.4$ , ³² where $h_0$ , $a$ and $g_{\infty }$ are the hardening parameters which determine the hardening behaviour of all the slip systems.

To implement the crystal plasticity constitutive model into the numerical solver, the applied shear stress ${\tau }^{\alpha }$ is related to the Cauchy stress tensor $\mathbf{\sigma }$ and the second Piola–Kirchhoff stress tensor $\mathbf{S}$ ²¹ and is expressed as:

$${\tau }^{\alpha }=\mathbf{\sigma }:\left({\mathbf{F}}_e{\mathbf{S}}_0^{\alpha }{\mathbf{F}}_e^{-1}\right)=\left({\mathbf{F}}_e^T{\mathbf{F}}_e\mathbf{S}\right):{\mathbf{S}}_0^{\alpha }$$

(8)

where $\mathbf{S}$ can be calculated using small elastic strain assumption:

$$\mathbf{S}=\mathbb{C}:{\mathbf{E}}_e$$

(9)

where ${\mathbf{E}}_e$ is the elastic Green strain tensor given by:

$${\mathbf{E}}_e=({\mathbf{F}}_e^T{\mathbf{F}}_e-\mathbf{I})/2$$

(10)

$\mathbb{C}$ is the elasticity matrix in the sample coordinate system. With respect to the crystal axis, this matrix is orthogonal symmetrical for the case of cubic crystal and is fully specified by three elastic modules $C_{11},C_{12}$ and $C_{44}$ .

Preprosessing

To build the full-field CPFEM model, a cube of dimension 1 mm³ was considered as the RVE, and a periodical Voronoi tessellation of 500 randomly placed seeds was used to generate 500 random grains using the open-source polycrystal generation and meshing software Neper. A total of 500 crystallographic orientations generated in accordance with orientation distribution functions (ODFs) which is calculated from the electron backscatter diffraction (EBSD) measurement. These grain orientations were randomly assigned to the individual grain cell of RVE. Zhang et al. ⁷ show the low root mean squared deviations (RMSD) between the reconstructed ODFs and experimentally determined ODFs, indicating that 500 crystallographic orientations are enough to represent the texture of materials.

As stated above, the cubic RVE with the crystallographic texture was used to perform the virtual uniaxial tensile tests and to obtain the uniaxial yield stresses and plastic potentials in various directions for the identification of various yield functions. In addition, the virtual plan-strain, pure shear and equibiaxial tensile tests were performed to obtain the anisotropic properties in biaxial stress states. To reproduce the experimental stress state in virtual tests, periodic boundary conditions applied on the faces of the RVE.

To perform the virtual tensile tests in various in-plane and out-of-plane directions with respect to the RD, applying periodic boundary conditions in loading direction was inconvenient to calculate the r-values due to the non-rectangular shapes of the deformed RVE. Therefore, the boundary conditions of the RVE were remain constant in space, while rotating the Euler angle of 500 grains is a equivalent method as rotating the loading direction for the cubic-orthorhombic FCC.^7,8

Phenomenological plasticity constitutive model

Yield function

The crystal plasticity method gives the essential understanding of plasticity in metal materials. However, this method requires large computational costs for the large-scale engineering problem. The phenomenological plasticity method is widely used to describe the yielding and plastic deformation of materials due to its simple expression and user-friendly implementation in FE modelling. An advanced phenomenological plane stress anisotropic yield function so-called Yld2000-2d proposed by Barlat et al. ² gained considerable popularity mainly because of its accurate prediction of yield stresses and r-values for various strong textured aluminium alloys. This yield function is based on multiple linear transformations of two unconditionally convex functions $\phi \prime$ and $\phi ″$ of deviatoric stress tensor. ²

$$f\left(\mathbf{\sigma }\right)={\left[\frac{1}{2}(\phi \prime +\phi ″)\right]}^{1/w}=\overline{\sigma }$$

(11)

where

$$\phi \prime ={({X\prime }_1-{X\prime }_2)}^w,\phi ″={(2{X″}_2+{X″}_1)}^w+{(2{X″}_1+{X″}_2)}^w$$

(12)

where the exponent $w$ is six for BCC and eight for FCC materials, respectively. The $X^{(k)}$ are the two principal values of the transformed stress tensor ${\tilde{\mathbf{x}}}^{(k)}$ which are obtained by employing linear transformations on Cauchy stress deviators in plane stress state, $\mathbf{\sigma }=[{\sigma }_{11}{\sigma }_{22}{\sigma }_{12}{]}^T$ :

$$\tilde{\mathbf{x}}\prime =\mathbf{L}\prime :\mathbf{\sigma }$$

(13)

$$\tilde{\mathbf{x}}″=\mathbf{L}″:\mathbf{\sigma }$$

(14)

where the linear transformations were defined as:

$$\left[\begin{array}{c}{L\prime }_{11}\hfill \\ \multicolumn{1}{c}{{L\prime }_{12}}\\ \multicolumn{1}{c}{{L\prime }_{21}}\\ \multicolumn{1}{c}{{L\prime }_{22}}\\ \multicolumn{1}{c}{{L\prime }_{66}}\end{array}\right]=\left[\begin{array}{ccc}2/3\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{-1/3}& 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& -1/3\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 2/3\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 1\hfill \end{array}\right]\left[\begin{array}{c}{a}_1\hfill \\ \multicolumn{1}{c}{a_2}\\ \multicolumn{1}{c}{a_7}\end{array}\right],$$

(15)

$$\left[\begin{array}{c}{L\prime }_{11}\hfill \\ \multicolumn{1}{c}{{L\prime }_{12}}\\ \multicolumn{1}{c}{{L\prime }_{21}}\\ \multicolumn{1}{c}{{L\prime }_{22}}\\ \multicolumn{1}{c}{{L\prime }_{66}}\end{array}\right]=\frac{1}{9}\left[\begin{array}{ccccc}2\hfill & -2\hfill & 8\hfill & -2\hfill & 0\hfill \\ \multicolumn{1}{c}{1}& -4\hfill & -4\hfill & 4\hfill & 0\hfill \\ \multicolumn{1}{c}{4}& -4\hfill & -4\hfill & 1\hfill & 0\hfill \\ \multicolumn{1}{c}{-2}& 8\hfill & 2\hfill & -2\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 0\hfill & 9\hfill \end{array}\right]\left[\begin{array}{c}{a}_3\hfill \\ \multicolumn{1}{c}{a_4}\\ \multicolumn{1}{c}{a_5}\\ \multicolumn{1}{c}{a_6}\\ \multicolumn{1}{c}{a_8}\end{array}\right]$$

(16)

The Yld2000-2d has shown the adequate flexibility for the calibration of in-plane mechanical anisotropy and superior accuracy of prediction in FE simulation compared with classical quadric yield functions.^34,35 In addition, Safaei et al. ¹⁰ presented the description of anisotropic plasticity of Yld2000-2d associated with non-AFR is more accurate than AFR. However, this yield function cannot account for the out-of-plane mechanical anisotropy. To fully describe the anisotropy of materials in 3D stress space, an anisotropic yield function so-called Yld2004-18p was developed by Barlat et al. ⁴ based on the established concept of multiple linear transformations of the stress deviator. This yield function was found to have good flexibility and high accuracy in describing mechanical anisotropy of strongly textured aluminium alloy sheets in 3D stress space.^4,6,36 The formulation of Yld2004-18p yield function is expressed as, ⁴

$$f\left(\mathbf{\sigma }\right)={\left(\frac{1}{4}\sum _{i=1}^3\sum _{j=1}^3{|{{\tilde{S}}_i}^{(1)}-{{\tilde{S}}_j}^{(2)}|}^w\right)}^{1/w}=\overline{\sigma }$$

(17)

The ${\tilde{S}}^{(k)}$ are the three principal values of the transformed stress tensors ${\tilde{\mathbf{s}}}^{(k)}$ which are obtained by employing linear transformations on Cauchy stress deviators in 3D stress space, $\mathbf{\sigma }=[{\sigma }_{11}{\sigma }_{22}{\sigma }_{33}{\sigma }_{12}{\sigma }_{13}{\sigma }_{23}{]}^T$ :

$${\tilde{\mathbf{s}}}^{(k)}={\mathbf{C}}^{(k)}:\mathbf{T}:\mathbf{\sigma }={\mathbf{L}}^{(k)}:\mathbf{\sigma }$$

(18)

where ${\mathbf{C}}^{(k)}$ are the associated linear transformations on the stress deviators which are used to describe the material anisotropy and $\mathbf{T}$ is the constant linear transformation for stress deviators. The multiple linear transformations ${\mathbf{L}}^{(k)}={\mathbf{C}}^{(k)}:\mathbf{T}$ are given by:

$$\begin{array}{c}{\mathbf{L}}^{(k)}=\left[\begin{array}{cccccc}0\hfill & -c_{12}^{(k)}\hfill & -c_{13}^{(k)}\hfill & 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{-c_{21}^{(k)}}& 0\hfill & -c_{23}^{(k)}\hfill & 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{-c_{31}^{(k)}}& -c_{32}^{(k)}\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & c_{44}^{(k)}\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 0\hfill & c_{55}^{(k)}\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 0\hfill & 0\hfill & c_{66}^{(k)}\hfill \end{array}\right]\hfill \\ ·\frac{1}{3}\left[\begin{array}{cccccc}2\hfill & -1\hfill & -1\hfill & 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{-1}& 2\hfill & -1\hfill & 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{-1}& -1\hfill & 2\hfill & 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 3\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 0\hfill & 3\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 0\hfill & 0\hfill & 3\hfill \end{array}\right]\hfill \end{array}$$

(19)

Since the full-field CP model can provide a larger number of yield stress points and plastic potentials data in 3D stress space, the Yld2004-3D yield function associated with AFR is difficult to accurately calibrate these materials sampling points simultaneously. Nevertheless, this limitation can be solved by the non-AFR in which the Yld2004-3D function is separately used as yield and plastic potential functions.

Non-associated flow rule framework

In classical plasticity theory, the AFR (please see Figure 3(a)) is employed to determine the flow direction of plastic strain increment $\mathrm{\Delta }{\epsilon }^p$ and equivalent plastic strain increment $\mathrm{\Delta }{\overline{\epsilon }}^p$ based on the first differential of the yield function $f_y$ ,

$$\mathrm{\Delta }{\epsilon }^p=\stackrel{·}{\lambda }·\frac{\partial f_y}{\partial \mathbf{\sigma }}$$

(20)

$$\mathrm{\Delta }{\overline{\epsilon }}^p=\frac{\mathbf{\sigma }:\mathrm{\Delta }{\epsilon }^p}{f_y}=\frac{\mathbf{\sigma }:\frac{\partial f_y}{\partial \mathbf{\sigma }}·\stackrel{·}{\lambda }}{f_y}=\frac{f_y·\stackrel{·}{\lambda }}{f_y}=\stackrel{·}{\lambda }$$

(21)

Therefore, the yielding and plastic potential are determined by the same yield surface. However, for strongly textured aluminium alloy, using the same yield surface is difficult to simultaneously calibrate the directional yield stresses and r-values.^10,13,20 Alternatively, these issues have been successfully resolved using non-AFR in which the flow stresses and plastic flow directions were separately determined by yield function and plastic potential function. In non-AFR (please see Figure 3(b)), the plastic strain increments are determined by the first differential of the plastic potential function $f_p$ :

$$\mathrm{\Delta }{\epsilon }^p=\stackrel{·}{\lambda }·\frac{\partial f_p}{\partial \mathbf{\sigma }}$$

(22)

$$\mathrm{\Delta }{\overline{\epsilon }}^p=\frac{\mathbf{\sigma }:\mathrm{\Delta }{\epsilon }^p}{f_y}=\frac{\mathbf{\sigma }:\frac{\partial f_p}{\partial \mathbf{\sigma }}·\stackrel{·}{\lambda }}{f_y}=\left(\frac{f_p}{f_y}\right)\stackrel{·}{\lambda }$$

(23)

Figure 3.

The schematics of stress integration return mapping for bringing back elastic predictor on objective yield surface following: (a) AFR framework (return mapping normal to yield surface) and (b) non-AFR framework (return mapping normal to plastic potential surface).

In addition, the selection of combination of yield and plastic potential functions under non-AFR determines the accuracy of plastic constitutive model. There were various combinations have been developed including quadratic Hill series functions ³⁷ and non-quadratic convex functions.^10,38

The next increment corrects error (NICE) explicit integration scheme

The proposed anisotropic yield and plastic potential functions associated with non-AFR can be implemented in the commercial FE software Abaqus/Explicit via a user-defined material subroutine (VUMAT). Generally, a computationally efficient classical forward-Euler (CFE) scheme is applied as the standard explicit integration method. However, the satisfaction of consistency condition, $d\mathcal{F}=0$ is in-accuracy. In contrary, the classical backward-Euler (CBE) scheme is employed as a semi-explicit method in which a large number of iterations are computed to fully satisfy the consistency condition at the end of each increment, while demands the high computational cost. To overcome the limitations of CFE and CBE, Halilovic et al. ³⁹ introduced an improved explicit integration method so-called a next increment corrects error (NICE) scheme, which can be employed to obtain the accurate stress integration, while saves the computational cost. In NICE scheme, the consistency condition, $\mathrm{d}\mathcal{F}=0$ , is replaced by:

$$\mathcal{F}+\mathrm{d}\mathcal{F}=0$$

(24)

The application of NICE has shown that (i) compared with CFE scheme, the relative errors of the final stress, plastic strain and porosity are significantly reduced; (ii) when small strain increment is adopted to achieve the fulfilment of consistency condition, the significant reduction of the CPU time consumption is achieved against CBE scheme. Therefore, the NICE technique is employed in the elastoplastic constitutive model associated with non-AFR via VUMAT subroutine.

At the beginning of each time increment in VUMAT, the total strain increment $\mathrm{\Delta }{\epsilon }_{\mathrm{ij}}$ are calculated from the computed increments of the displacement field. If the assumption of plane stress is accepted, the strain increment $\mathrm{\Delta }{\epsilon }_{33}$ in thickness direction needed to be determined in VUMAT. For elastic deformation mode, the $\mathrm{\Delta }{\epsilon }_{33}$ is given by,

$$\mathrm{\Delta }{\epsilon }_{33}=\frac{-\nu }{1-\nu }(\mathrm{\Delta }{\epsilon }_{11}+\mathrm{\Delta }{\epsilon }_{22})$$

(25)

while the fully elastic predictor stress $\mathrm{\Delta }{\sigma }_{\mathrm{ij}}^b$ (linear orthotropic) is determined:

$$\mathrm{\Delta }{\sigma }_{\mathrm{ij}}^b={\mathbf{C}}_{\mathrm{ijkl}}\mathrm{\Delta }{\epsilon }_{\mathrm{ij}}$$

(26)

where ${\mathbf{C}}_{\mathrm{ijkl}}$ is the elastic stiffness tensor. Here, the yield criterion is defined as:

$${\mathcal{F}}_{n+1}=f_y({\sigma }_{\mathrm{ij}}^n+\mathrm{\Delta }{\sigma }_{\mathrm{ij}}^b)-{\sigma }_0({\overline{\epsilon }}_n^p)$$

(27)

when yield criteria ${\mathcal{F}}_{n+1}\le 0$ in which stress state point is within the limit of yield surface in stress space, the deformation mode is supposed to be fully elastic, and hence, the stress tensor is updated using equation (26). In contrary, the plastic deformation mode is active when ${\mathcal{F}}_{n+1}>0$ . In this case, a corrected stress tensor increment brings the elastic predictor back to the yield surface. A truncated Taylor expansion is applied to the consistency condition (equation (24)) to obtain the expression of corrected stress tensor increment: ³⁹

$$\mathcal{F}+\frac{\partial \mathcal{F}}{\partial {\sigma }_{\mathrm{ij}}}\mathrm{\Delta }{\sigma }_{\mathrm{ij}}+\frac{\partial \mathcal{F}}{\partial \mathrm{Y}}\mathrm{\Delta }\mathrm{Y}=0$$

(28)

where the corrected $\mathrm{\Delta }{\sigma }_{\mathrm{ij}}$ and the hardening modules increment $\mathrm{\Delta }\mathrm{Y}$ is computed by:

$$\mathrm{\Delta }{\sigma }_{\mathrm{ij}}={\mathbf{C}}_{\mathrm{ijkl}}:\mathrm{\Delta }{\epsilon }_{\mathrm{ij}}^e={\mathbf{C}}_{\mathrm{ijkl}}:(\mathrm{\Delta }{\epsilon }_{\mathrm{ij}}-\mathrm{\Delta }{\epsilon }_{\mathrm{ij}}^p)$$

(29)

$$\mathrm{\Delta }\mathrm{Y}=\frac{\partial \mathrm{Y}}{\partial \overline{{\epsilon }^p}}·\mathrm{\Delta }{\overline{\epsilon }}^p$$

(30)

Here, the plastic tensor increment $\mathrm{\Delta }{\epsilon }_{\mathrm{ij}}^p$ and equivalent strain increment $\mathrm{\Delta }{\overline{\epsilon }}^p$ associated with the non-AFR are shown in equations (22) and (23). Noted that, only one independent variable $\mathrm{\Delta }\lambda$ so-called plastic multiplier is undetermined. Substituting the equations (22), (23), (29) and (30) into the expansion of consistency condition equation (28), the $\mathrm{\Delta }\lambda$ is given by,

$$\mathrm{\Delta }\lambda =\frac{\mathcal{F}+{\left(\frac{\partial f_y}{\partial {\sigma }_{\mathrm{ij}}}\right)}^{\mathrm{T}}:{\mathbf{C}}_{\mathrm{ijkl}}:\mathrm{\Delta }{\epsilon }_{\mathrm{ij}}}{{\left(\frac{\partial f_y}{\partial {\sigma }_{\mathrm{ij}}}\right)}^{\mathrm{T}}:{\mathbf{C}}_{\mathrm{ijkl}}:\left(\frac{\partial f_p}{\partial {\sigma }_{\mathrm{ij}}}\right)-\frac{\partial f_y}{\partial \mathrm{Y}}·\frac{\partial \mathrm{Y}}{\partial \overline{{\epsilon }^p}}·\frac{{\sigma }_{\mathrm{ij}}:\frac{\partial f_p}{\partial {\sigma }_{\mathrm{ij}}}}{f_y}}$$

(31)

when the plane stress assumption is accepted, the strain increment in thickness direction $\mathrm{\Delta }{\epsilon }_{33}$ is solved using the condition $\mathrm{\Delta }{\sigma }_{33}=0$ and is computed as:

$$\mathrm{\Delta }{\epsilon }_{33}=\frac{{\mathrm{\Theta }}_{11}+{\mathrm{\Theta }}_{22}}{C_{3333}-{(C_{3311}\frac{\partial f_p}{\partial {\sigma }_{11}}+C_{3322}\frac{\partial f_p}{\partial {\sigma }_{22}})}^2\mathrm{\Gamma }}$$

(32)

where the variables ${\mathrm{\Theta }}_{\mathrm{ij}}$ and $\mathrm{\Gamma }$ are given by,

$${\mathrm{\Theta }}_{\mathrm{ij}}=C_{33\mathrm{ij}}\mathrm{\Delta }{\epsilon }_{\mathrm{ij}}-\frac{\partial f_p}{\partial {\sigma }_{\mathrm{ij}}}(f_y+\left(\frac{\partial f_y}{\partial {\sigma }_{\mathrm{ij}}}\right)·{\mathbf{C}}_{33\mathrm{ij}}·\mathrm{\Delta }{\epsilon }_{\mathrm{ij}})\mathrm{\Gamma }$$

(33)

Additionally, to clarify the effectiveness of non-AFR in prediction of plasticity, the AFR based NICE explicit stress integration also has been developed by replacing the equations (22) and (23) with equations (20) and (21).

The hierarchical framework of the multi-level model

The CPFEM constitutive model was implemented into the general commercial FE solvers, ABAQUS/Standard. The anisotropic properties of textured aluminium alloy were explored by the virtual material test by loading the RVE in many deformation modes, for example, uniaxial tension, plane-strain tension, equibiaxial tension and pure shear. After each simulation, post-processing was implemented to read the stress tensor and strain rate for each loading case at the desired level of plastic work and respectively used as yield stress and plastic potential points for identification of yield functions. To find the optimal parameters of the selected yield functions, a fitting procedure was employed using the Levenberg-Marquardt nonlinear least square method (NLSM) in which the number of yield stress and plastic potential points is equal to or greater than the number of parameters. It was solved by minimising the residual value of $\mathcal{R}(a)$ and $\mathcal{R}({\mathbf{C}}^{(k)})$ for Yld2000 and Yld2004, respectively:

$$\mathcal{R}({\mathbf{C}}^{(k)}),\mathcal{R}(a)=\frac{1}{{\mathrm{n}}_{\mathrm{Y}}+{\mathrm{n}}_{\mathrm{r}}}\left\{\begin{array}{c}{\sum }_{i=1}^{{\mathrm{n}}_{\mathrm{Y}}}{(\frac{{\sigma }_f^{\mathrm{Y}}}{{\sigma }_{\mathrm{CP}}^{\mathrm{Y}}}-1)}^2+\hfill \\ \multicolumn{1}{c}{{\sum }_{i=1}^{{\mathrm{n}}_{\mathrm{r}}}{(\frac{r_f}{r_{\mathrm{CP}}}-1)}^2}\end{array}\right\}\Rightarrow \mathrm{Min}$$

(34)

where ${\mathrm{n}}_{\mathrm{Y}}$ and ${\mathrm{n}}_{\mathrm{r}}$ are the number of yield stress and plastic potential points obtained from the virtual material tests, respectively. ${\sigma }_f^{\mathrm{Y}}$ and $r_f$ are the yield function calculated yield stress and plastic potential, respectively. Their counterparts predicted by the CP model are respectively defined as ${\sigma }_{\mathrm{CP}}^{\mathrm{Y}}$ and $r_{\mathrm{CP}}$ . As a result, the calibrated yield functions were implemented into FEM solvers by using AFR and Non-AFR stress integration schemes to perform the sheet metal forming simulations.

Materials of multi-level model

Material characterisation

Here, the full-filed CPFEM was applied to investigate the texture induced anisotropic plasticity of strong textured AA3104 aluminium alloy sheet (an Al-alloy with about 1% manganese and 1% magnesium^40,41). This sheet was a typical canbody stock and was widely used in metal bodymaker industry due to its good formability and strength. To meet the strength (≥275 MPa) and low gauge (≤0.4 mm) requirements of aluminium can, canbody stock AA3104 is fabricated through a thermo-mechanical rolling process route from the direct chill, casting, homogenised, ingot that is hot rolled to 2.0 mm with a heavily strain hardened, hard rolled H19 temper (87–89% cold rolling reduction) to final gauge about 0.27–0.40 mm. Therefore, a pronounced rolling texture is acquired during thermo-mechanical processing, where the preferred crystallographic orientations are assembled along the so-called $\beta$ fibre which starts from the copper orientation $\{112\}〈111〉$ through the s texture $\{123\}〈634〉$ to the brass texture $\{011\}〈211〉$ in 3D Euler angle space.

To measure the crystallographic texture, a specimen with dimensions of 5 mm × 3 mm × 0.4 mm was prepared for crystallographic texture and orientation analyses using a FEI Quanta 3D field emission gun scanning electron microscope (FEG-SEM) equipped with an EDAX-TSL Hikari EBSD system. During the measurement, the EBSD mapping area was assigned to be 1 × 1 mm with a step size of 1 µm. The sliced orientation density maps were obtained from ODFs through the reduced Euler space $0°\le \{{\varphi }_1,\Phi ,{\varphi }_2\}\le 90°$ in Bunge’s convention due to the cubic-orthorhombic symmetry of aluminium sheet. The sliced orientation density maps with respect to ${\varphi }_2$ are given in Figure 4. A Gaussian spread of 12° was used to quantify the volume fractions of the different texture components using texture analysis software MTEX ⁴² and given by:

1. 18.51% $\{011\}〈211〉$ Brass texture

2. 22.51% $\{123\}〈634〉$ S texture

3. 4.91% $\{112\}〈111〉$ Copper texture

4. 11.26% $\{001\}〈100〉$ Cube texture

Figure 4.

The orientation density map of 500 reconstructed grains from the orientation density function (ODF) of cold-rolled AA3104-H19 aluminium sheet as sections through the reduced Euler space for cubic–orthorhombic symmetry (constant ${\varphi }_2$ sections, 0° to 85° in 5° steps, Euler angles in Bunge’s convention).

Therefore, the volume fractions 45.93% of cold rolled texture, in which the crystal orientation density is high along the $\beta$ fibre, is acquired by the cold rolling process from recrystallised hot rolled band.

Here, a total of 500 crystallographic orientations generated by ODFs were randomly assigned into these grain cells of the RVE. To obtain a high resolution for each grain cell, the RVE discretised with 2,032,594 tetrahedron elements with about 4065 elements per grain in average and was shown in Figure 5.

Figure 5.

A RVE generated by periodical Voronoi tessellation of 500 randomly placed seeds at resolution of 4065 tetrahedron elements each cell and each Voronoi cell (shown in various colour) represents one grain orientation for AA3104-H19 with a total of 500 crystallographic orientations generated by ODFs.

Material constants of CPFEM

The constitutive material constants of CPFEM were fitted through adjusting the predicted stress-strain curve to the experimental tests by applied three periodical boundary conditions on the surfaces of RVE, for example, uniaxial, plane-strain tension and pure shear in RD. Such experimental tests were performed on the well-designed AA3104-H19 specimens cut from the same sheet, that is, uniaxial tensile experiments were conducted following the ASTM-E8-E8M13 standard ⁴³ and the specimen was tested using an Instron 1342 servo-hydraulic test system at the loading speed of 0.005 mm/s. Artificial speckle patterns were sprayed on the specimen surfaces and the loading processes were recorded for strain measurement using digital image correlation (DIC) method (please see Figure 6(a)). The CCD camera frame rate was configured at 20fps, hence, the uniaxial stress versus Misses equivalent strain curve was measured. The results show the uniaxial yield stress of AA3104-H19 in RD is about 262 MPa. The calculation of uniaxial r-value using the plastic incompressibility condition is given by:

$$r=\frac{{\stackrel{·}{\epsilon }}_w}{{\stackrel{·}{\epsilon }}_t}=-\frac{{\stackrel{·}{\epsilon }}_w}{{\stackrel{·}{\epsilon }}_l+{\stackrel{·}{\epsilon }}_w}$$

(35)

where ${\stackrel{·}{\epsilon }}_t$ , ${\stackrel{·}{\epsilon }}_w$ and ${\stackrel{·}{\epsilon }}_l$ are the true strain rates in the thickness, width and longitudinal directions, respectively. Here, the ${\stackrel{·}{\epsilon }}_w$ and ${\stackrel{·}{\epsilon }}_l$ are determined by the true strain increment measured using DIC between each CCD camera frame. Moreover, the plane-strain specimens were developed to study the yielding in plane-strain stress state.^44,45 Figure 6(b) shows the equivalent strain distributions in plane strain region zone were measured by DIC. The yield stress of plane-strain tensile $\sigma \mathrm{ps}$ is calculated by:

$${\mathrm{F}}_{\mathrm{t}}={\sigma }_{\mathrm{ps}}·\mathrm{t}·{\mathrm{W}}_{\mathrm{c}}+2{\mathrm{F}}_{\mathrm{e}}$$

(36)

where, ${\mathrm{F}}_{\mathrm{t}}$ and ${\mathrm{F}}_{\mathrm{e}}$ are the total force and edge effects induced force, respectively. $\mathrm{t}$ and ${\mathrm{W}}_{\mathrm{c}}$ respectively denote the thickness and width of the plane strain region. Since ${\mathrm{F}}_{\mathrm{t}}$ is a linear function of ${\mathrm{W}}_{\mathrm{c}}$ , the slope ( ${\sigma }_{\mathrm{ps}}·\mathrm{t}$ ) of the function can be determined by subtracting ${\mathrm{F}}_{\mathrm{e}}$ using two different widths of plane-strain specimens, and hence, the yield stresses in plane-strain tension can be measured. In addition, pure shear tests were performed and the shear stresses were measured in accordance with ASTM-B831-05 standard. ⁴⁶ The equivalent strain distributions in the shear zone were measured by DIC as shown in Figure 6(c). As a result, the experimental relationship between the equivalent strain and flow stress in various stress states were recorded during the loading process and reported in Figure 7. Consequently, the CPFEM material constants of AA3104-H19 were determined by these experimental tensile tests and given in Table 2. The validation of CPFEM was performed by employing these boundary conditions on RVE (please see Figure 6(d) to (f)) and reading the corresponding stress-strain curves (please see Figure 7). As a result, Figure 7 indicates the hardening behaviours of AA3104-H19, which are accurately predicted by the CPFEM in various periodical boundary conditions.

Figure 6.

Mises equivalent strains measured by DIC method for: (a) uniaxial tensile experiment, (b) plane-strain tensile experiment, (c) pure shear experiment, the distribution of Mises equivalent stresses of RVE predicted by CPFEM applied with various boundary conditions for, (d) virtual uniaxial tension, (e) virtual plane-strain tension and (f) virtual pure shear.

Figure 7.

Comparison of stress–strain curves for uniaxial tension, plane-strain tension and pure shear obtained by experimental tensile tests (solid lines with various colours) and predicted by virtual CPFEM tests (dot lines with various colours).

Table 2.

CP model parameters for the cold rolled AA3104-H19 aluminium alloy.

$C_{11}$	$C_{12}$	$C_{44}$	${\stackrel{·}{\gamma }}_0$	$h_0$
206 GPa	118GPa	54 GPa	0.001 s–1	2493 MPa
${\tau }_0$	$g_{\infty }$	$m$	$a$	$q$
98 MPa	126 MPa	0.02	2.0	1.4

Forming experiment and simulation

To operate the standard cup drawing experiments of round blank, a tooling system with punch, drawing die and blank holder was designed and manufactured (please see Figure 8(a)). The configuration and dimensions of the experimental setup are summarised in Figure 8(c). The round blanks with a radius 30 mm was cut from a 0.4 mm AA3104-H19 sheet using electro-discharge machining, and then they were deep drawn using a single stroke with a speed of 0.1 mm/s. The blank holder force about 5.0 kN and lucubration were employed to avoid flange wrinkling and surface scratching, respectively.

Figure 8.

The schematics of forming experiment and simulation for: (a) deep drawing tooling system with a punch, a drawing die and a blank, (b) mesh quality of initial blank discretised with 7839 elements each layer and (c) configuration and dimension of the tooling system.

To numerically evaluate the out-of-plane anisotropy as well as plasticity flow rule, the VUMAT subroutines of CPFEM determined plane stress Yld2000-2d and 3D Yld2004-18p yield functions associated with AFR and non-AFR were separately used to simulate deep drawing of an AA3103-H19 sheet. For the mesh of blank, 3D stress C3D8R and plane stress S4R element type were employed in Yld2000-2d and Yld2004-18p models, respectively (please see Figure 8(b)). Punch, blank holder and drawing die were constrained as rigid bodies. Friction conditions between the sheet metal and tooling system are important factors that affect the prediction of sheet metal forming simulation. The isotropic Coulomb friction model with a coefficient $\mu$  = 0.025 was used to describe the contact interface between the blank and tooling system. This friction condition was estimated by adjusting the predicted punch forces to the experimental data.

Results

Validations

The experimental tensile tests in various directions were performed to validate the anisotropic properties predicted by the CPFEM. Here, the specimens used for uniaxial tension and shear tension were prepared from the same sheet at an angular interval of 15° from the RD. Plane strain specimens were prepared at an angular interval of 45° from the RD. For the CPFEM prediction, the anisotropic properties were determined by rotating the grain orientation. The r-value was calculated using the strain rates which are determined by the deformation rates of cubic RVE between $2\%$ and $10\%$ elongations. The influence of plastic work induced texture evolution on strain hardening characteristics was reported as negligible for AA3104-h19 (experimental results see reference ⁴⁷ ) due to all slip systems in the FCC lattice belong to the $\{110\}〈111〉$ family. Hence, only initial yield stresses in various directions were studied for the characterisation of strength anisotropy.

Figure 9(a) to (d) show a comparison of the anisotropic properties predicted by the virtual tests with those from mechanical experiments. It can be seen that the directional r-values and uniaxial stresses predicted by CPFEM coincided with those experimental results (please see Figure 9(a) and (b)). Additionally, the results showed r-value appear two major peaks near $\theta =55°$ , $90°$ , and a minor peak at $0°$ . AA3104-H19 sheet generally exhibits a strong rolling texture and moderate Cube texture. According to the work of Darrieulat and Piot, ⁴⁸ the r-value of material with the strong rolling texture generally appear a major peak near $\theta =45°$ and those of cube texture show two major peaks in RD and TD. Therefore, the texture components of AA3104-H19 was designed to reduce the anisotropy level through controlling the fraction of rolling and cube texture during the thermo-mechanical processing to meet the requirement of canbody maker. For the strength anisotropy, the predicted uniaxial yield stresses increased monotonically from RD to TD and showed the maximum normalised yield strength is 1.058. The comparison between the plane strain and pure shear yield stresses obtained by the virtual tests and the experimental tests were also plotted in Figure 9(c) and (d), respectively. The normalised yield stresses of virtual plane-strain tensile tests were monotonically increasing from minimum 1.088 in RD to maximum 1.203 in TD. For the directional normalised yield stresses of pure shear, the peak value of 0.567 was found near $45°$ , while minimum values round 0.532 were observed at both $0°$ and $90°$ .

Figure 9.

Comparisons of in-plane mechanical anisotropies obtained by experimental tensile tests and predicted by virtual CPFEM tests in various stress states for: (a) normalised uniaxial yield stress variation from RD to TD, (b) normalised uniaxial r-value variation from RD to TD, (c) normalised plane-strain yield stress variation from RD to TD and (d) normalised pure shear yield stress from RD to TD.

As a result, the predictions of the virtual tests operated by the full-field CPFEM showed good agreement with those experimental results. It indicates that the virtual material tests could be an alternative to the extensive and costly physical tests.

Calibration of yield and plastic potential functions

In this section, the mechanical anisotropies of the AA3104-H19 sheet were described using the advanced plan stress Yld2000-2d and 3D stress Yld2004-18p yield functions associated with both AFR and non-AFR. The calibrations of plan stress Yld2000-2d yield function were determined using the in-plane anisotropy in [RD-TD] plane predicted by the virtual tests. The calibrated parameters of Yld2000-2d are documented in Table 3. Additionally, not only the in-plane anisotropy in [RD-TD] plane but also the out-of-plane anisotropy in [RD-ND] and [TD-ND] planes were involved in the determination of parameters for 3D stress Yld2004-18p yield function. The calibrated parameters of Yld2004-18p are documented in Table 4. Here, the uniaxial yield stress at specific plastic work of 1 MPa and r-values between the plastic strain of 2%–10% were used to determine the parameters of yield functions.

Table 3.

The parameters of Yld2000-2d for AA3104-H19 aluminium alloy.

Parameters for the YLD2000 function using AFR framework
$a_1$	$a_2$	$a_3$	$a_4$	$a_5$	$a_6$	$a_7$	$a_8$
0.855021	1.005011	1.039419	0.943642	1.029384	1.032025	0.974824	1.001533
Parameters for the YLD2000 function using NAFR framework
Yield function:
$a_1$	$a_2$	$a_3$	$a_4$	$a_5$	$a_6$	$a_7$	$a_8$
1.073683	0.819553	0.737751	0.915491	1.039769	1.125420	0.775537	1.150539
Plastic potential function
$a_1$	$a_2$	$a_3$	$a_4$	$a_5$	$a_6$	$a_7$	$a_8$
1.086868	0.772367	1.567497	1.006626	1.082278	1.229948	0.817586	0.343555

Table 4.

The parameters of Yld2004-18p for AA3104-H19 aluminium alloy.

Parameters for the YLD2004 function using AFR framework
$c_{12}^{(1)}$	$c_{13}^{(1)}$	$c_{21}^{(1)}$	$c_{23}^{(1)}$	$c_{31}^{(1)}$	$c_{32}^{(1)}$	$c_{44}^{(1)}$	$c_{55}^{(1)}$	$c_{66}^{(1)}$
0.967136	0.945800	0.899031	0.954927	1.004614	0.573419	1.103712	1.135837	0.985791
$c_{12}^{(2)}$	$c_{13}^{(2)}$	$c_{21}^{(2)}$	$c_{23}^{(2)}$	$c_{31}^{(2)}$	$c_{32}^{(2)}$	$c_{44}^{(2)}$	$c_{55}^{(2)}$	$c_{66}^{(2)}$
1.022051	0.706361	0.608427	0.576953	1.156761	1.088780	0.827078	0.816403	1.073345
Parameters for the YLD2004 function using NAFR framework
Yield function
$c_{12}^{(1)}$	$c_{13}^{(1)}$	$c_{21}^{(1)}$	$c_{23}^{(1)}$	$c_{31}^{(1)}$	$c_{32}^{(1)}$	$c_{44}^{(1)}$	$c_{55}^{(1)}$	$c_{66}^{(1)}$
1.201842	0.791908	0.728014	0.645056	1.244460	0.587453	1.127995	1.466560	0.400051
$c_{12}^{(2)}$	$c_{13}^{(2)}$	$c_{21}^{(2)}$	$c_{23}^{(2)}$	$c_{31}^{(2)}$	$c_{32}^{(2)}$	$c_{44}^{(2)}$	$c_{55}^{(2)}$	$c_{66}^{(2)}$
0.890145	0.297731	0.767302	0.689291	1.320437	1.394671	0.719779	0.002018	1.271461
Plastic potential function
$c_{12}^{(1)}$	$c_{13}^{(1)}$	$c_{21}^{(1)}$	$c_{23}^{(1)}$	$c_{31}^{(1)}$	$c_{32}^{(1)}$	$c_{44}^{(1)}$	$c_{55}^{(1)}$	$c_{66}^{(1)}$
0.808495	0.732755	0.887072	1.062561	0.589963	−0.416623	1.356548	0.876497	0.960388
$c_{12}^{(2)}$	$c_{13}^{(2)}$	$c_{21}^{(2)}$	$c_{23}^{(2)}$	$c_{31}^{(2)}$	$c_{32}^{(2)}$	$c_{44}^{(2)}$	$c_{55}^{(2)}$	$c_{66}^{(2)}$
0.827461	0.473721	0.224387	1.518515	1.285679	0.882250	0.738844	1.997928	1.050068

Figure 10(a) and (b) show the yield functions calculated and CPFEM predicted uniaxial yield stresses and r-values in [RD-TD] plane, respectively. The CPFEM prediction showed that uniaxial stresses are slightly increasing from RD to TD and appear two peaks at around 30° and 75°, while r-values significantly increase from RD to TD and appear two peaks at about 50° and 90°. These results indicated that Yld2000-2d and Yld2004-18p associated with non-AFR successfully describe the variation of uniaxial stresses and r-values of AA3104-H19. In contrary, the deviation of fitting using AFR was significantly greater than non-AFR for both the Yld2000-2d and Yld2004-18p yield functions.

Figure 10.

Comparisons of in-plane and out-of-plane mechanical anisotropies predicted by virtual CPFEM tests and calculated by various yield functions for: (a) normalised uniaxial yield stresses variation in [RD-TD] plane, (b) uniaxial r-values variation in [RD-TD] plane, (c) normalised uniaxial yield stresses variation in [RD-ND] plane, (d) uniaxial r-values variation in [RD-ND] plane, (e) normalised uniaxial yield stresses variation in [ND-TD] plane and (f) uniaxial r-values variation in [ND-TD] plane.

In addition, although the accuracy of calibration using Yld2000-2d associated with non-AFR was superior to other functions in [RD-TD] planes, Yld2000-2d cannot describe the out-of-plane anisotropy. Therefore, only Yld2004-18p associated with AFR and non-AFR were used to calibrate the CPFEM predicted uniaxial yield stresses and r-values in both [RD-ND] and [TD-ND] planes, as shown in Figure 10(c) to (f). Figure 10(c) and (d) show that CPFEM predicted directional uniaxial yield stresses in [RD-ND] plane are slightly stronger at 60° and weakest at 0°, while the predicted r-values are almost constant from 0° to 50° and then significantly increase to about 1.0 in ND. For the mechanical anisotropy in [TD-ND] plane (please see Figure 10(e) and (f)), the predicted normalised uniaxial yield stress decreased from 1.061 in TD to 1.011 at 55° and then slightly increased to 1.022 in ND. Moreover, the predicted r-value is lowest at 45° and greatest in TD. The comparison of out-of-plane mechanical anisotropy between prediction of the CPFEM and the phenomenological yield functions indicated that the accuracy of fitting for out-of-plane anisotropy using Yld2004-18p yield function associated with non-AFR is superior to AFR. Figure 11(a) and (b) show the distribution of directional uniaxial yield stresses in 3D polar coordinates generated by Yld2004-18p yield functions associated with AFR and non-AFR, respectively. It can be seen that yield function associated with non-AFR provided greater flexibility to accurately capture the variation of uniaxial yield stresses aligned various 3D spatial directions (please see Figure 10(a), (c) and (e)). Essentially, it is due to non-AFR enable the yield function only account for yield stress, unlike the AFR in which the yield surface is simultaneously determined by the yield stresses and plastic potentials.

Figure 11.

The variation of uniaxial yield stresses of various materials orientations in 3D polar coordinates generated by: (a) Yld2004-18p associated with AFR and (b) Yld2004-18p associated with non-AFR.

For the mechanical anisotropy in multi-axial stress states, the yield surfaces generated by considering all phenomenological yield functions using yield stress points provided by the virtual multi-axial tensions in [RD-TD], [RD-ND] and [TD-ND], respectively. Figure 12(a) shows a comparison between the predicted yield surfaces and the 2D yield stress points in [RD-TD] plane generated by the virtual tests at a plastic work per unit volume of 1 MPa. It can be seen that Yld2004-18p yield functions associated with both AFR and non-AFR can accurately fit the 61 yield stress points provided by virtual tests in [RD-TD] plane. However, the Yld2000-2d associated with both AFR and non-AFR showed a minor difference to the CPFEM predictions being a little more ovality in the region around equibiaxial stress state. This is because the Yld2000-2d yield function is less flexible than the Yld2004-18p due to its smaller number of fitting parameters. As illustrated in Figure 12(b) and (c), 62 and 60 biaxial yield stress points in [RD-ND] and [TD]-[ND], respectively, were applied to generate various yield surfaces using the Yld2004-18p yield function associated with AFR or non-AFR. It can be seen that there is a good agreement between the yield stress points and yield surfaces determined by the Yld2004-18p yield function associated with both AFR or non-AFR.

Figure 12.

Comparison of yield stress points predicted by virtual CPFEM tests and yield surface calculated by various yield functions for various loading rates of biaxial stress states for: (a) normalised ${\sigma }_{11}$ and ${\sigma }_{22}$ associated with [RD-TD] plane, (b) normalised ${\sigma }_{11}$ and ${\sigma }_{33}$ associated with [RD-ND] plane and (c) normalised ${\sigma }_{22}$ and ${\sigma }_{33}$ associated with [TD-ND] plane.

Figure 13(a) to (d) show the 3D yield surfaces in the stress space of [ ${\sigma }_{11}$ – ${\sigma }_{22}$ – ${\sigma }_{12}$ ] using all considered yield functions. It can be seen that the yield surfaces predicted by Yld2000-2d (AFR) and Yld2004-18p (AFR and Non-AFR) are similar, while that predicted by Yld2000-2d (NAFR) shows higher convexity around the region of pure shear tension. Compared with the experimental maximum normalised yield stress with 0.524 in pure shear tension (please see Figure 9(c)), all considered yield functions overestimate the maximum shear stress, for example, 0.558, 0.688, 0.564 and 0.589 predicted by Yld2000-2d (AFR and non-AFR) and Yld2004-18p (AFR and Non-AFR), respectively. Therefore, in contrast to the remaining yield functions, the Yld2000-2d yield function associated with non-AFR significantly overrates the effect of shear stress.

Figure 13.

The multi-axial yield surface contour maps associated with various shear stress levels plotted in the three dimensional stress space [ ${\sigma }_{11}$ – ${\sigma }_{22}$ – ${\sigma }_{12}$ ] calculated by various yield functions: (a) Yld2000-2d associated with AFR, (b) Yld2000-2d associated with non-AFR, (c) Yld2004-18p associated with AFR and (d) Yld2004-18p associated with non-AFR.

Earing profile predictions

This subsection presents the results of AA3104-H19 cup drawing simulations using both yield functions Yld2000-2d and Yld2004-18p which are implemented in the AFR and non-AFR based NICE scheme linked to the ABAQUS/Explicit via the user material subroutines VUMAT. Figure 14(a) shows the predicted earing profiles of all considered models in cylindrical cup-drawing and those results were compared with the experimentally determined cup earing profiles. All experimental cup earing profiles exhibit six major ears located at $0°$ , $50°$ , $120°$ , $180°$ , $230°$ and $310°$ which are determined by the crystallographic textures induced variation of directional yield stresses and r-values. For the numerical results, none of the AFR based models presented variations could predict more than four ears due to less flexibility, while all non-AFR based models successfully predict the number and position of ears. Nevertheless, Yld2000-2d associated with non-AFR overrates the variation of cup height, especially at around the $0°$ and $180°$ . Moreover, the earing profile obtained from the non-AFR based Yld2004-18p is in good agreement with the experimental results.

Figure 14.

The results of deep drawing simulations for AA3104-H19: (a) comparison of earing profile achieved by experimental deep drawing tests and predicted by FE simulations using various virtual CPFEM tests determined yield functions and (b) the RMSD of the earing profile predicted by various yield function using the number of cup height data from 0° to 360° at an interval of 1° in comparison with experimental results.

In order to quantify the accuracy of predictions, the RMSD of cup height is calculated using following equation,

$$\mathrm{RMSD}=\sqrt{\frac{{\sum }_{i=1}^N{(H_{\mathrm{cup}}^f(i)-H_{\mathrm{cup}}^{\exp }(i))}^2}{N}}$$

(37)

where $H_{\mathrm{cup}}^f$ and $H_{\mathrm{cup}}^{\exp }$ are the cup height at a specific position determined from simulations and experiments, respectively. Here, the RMSD of the earing profile prediction using the number of cup height data from 0° to 360° at an interval of 1° is documented in Figure 14(b). The comparison indicates that AFR based Yld2000-2d and Yld2004-18p exhibits the RMSD with 0.331 and 0.251, respectively. Additionally, when using non-AFR take the place of AFR the RMSD of Yld2000-2d and Yld2004-18p are reduced to 0.215 and 0.196, respectively. Therefore, from Figure 14(b) one may conclude that the minimum RMSD of earing profile prediction is obtained using Yld2004-18p associated with non-AFR, that is, the capability of Yld2004-NAFR in plastic deformation prediction is superior to others. The equivalent plastic strain distribution in the deformed configuration for completed stage of the cup-drawing process using Yld2000-2d and Yld2004-18p (associated with AFR and non-AFR) are shown in Figure 15(a) to (d), respectively.

Figure 15.

The configurations and equivalent plastic strain distribution of the deep drawn cups predicted by FE simulations using various virtual CPFEM tests determined yield functions: (a) Yld2000-2d associated with AFR, (b) Yld2000-2d associated with NAFR, (c) Yld2004-18p associated with AFR and (d) Yld2004-18p associated with NAFR.

Discussions

The multi-level model with hierarchical structure was used to study the crystallographic texture induced anisotropy and its effects on plastic deformation of AA3104-H19. The validation of full-field CPFEM predictions was presented and the results showed good agreement with experimental results with respect to various deformation modes for the AA3104-H19 sheet. It is interesting to compare the fitting deviation using different plastic flow rules for both Yld2000-2d and Yld2004-18p. The results obtained with the calibrated Yld2004-18p and Yld2000-2d associated with non-AFR indicated that both yield functions were capable of describing the in-plane plastic anisotropy determined by the underlying CPFEM with good accuracy. Simultaneously, CPFEM predicted the out-of-plane can be captured by Yld2004-18p using non-AFR well. For the AFR based model, neither Yld2000-2d nor Yld2004-18p can capture the sampling points of materials in 2D or 3D stress space. All considered yield functions reproduced the curvature and radius of the yield surface closely to what was determined by the CPFEM models.

It was interesting to note that the influence of out-of-plane anisotropy on plastic deformation is not negligible in sheet metal forming simulations. In other words, the comparison of earing profile prediction using Yld2000-2d (only consider in-plane anisotropy) and Yld2004-18p (consider with in-plane and out-of-plane anisotropy) in Figure 14(a) indicated that the accuracy of plastic deformation prediction is improved when the employed yield function involved with the out-of-plane anisotropy. This was because the black holder compression force is applied on the flange region of the sheet metal in the thickness direction and an out-of-plane force happens in the bending region during the deep drawing process.

Additionally, it was quite clear that the accuracy of both Yld2004-18p and Yld2000-2d associated with non-AFR is better than that of AFR. Therefore, from the results one may conclude that when the non-AFR Yld2004-18p yield function is employed the accuracy of this hierarchical multi-scale model is mainly determined by the quality of the full-field CPFEM in reproducing the experimentally measured mechanical anisotropy. In other words, when a high-resolution crystallographic texture is available for the CPFEM, the Yld2004-18p associated with non-AFR seems capable of giving a realistic description of the texture induced plastic anisotropy in 3D stress space.

Since the CP determined Yld2004-18p yield function associated with non-AFR is superior in the prediction of plastic deformation during the sheet metal forming process, this constitutive model could be readily applied to perform some design oriented parametric studies in the future investigation, for example, tooling system optimisation, initial blank shape design and texture controlling for sheet metal rolling processes.

Conclusions

This work proposed a hierarchical framework of the multi-level model, that is, a multi-scale scheme, in which phenomenological yield function and full-field CPFEM are employed at the continuum scale and mesoscale, respectively. In this method, advanced phenomenological yield functions were successfully implemented in FE models and can be used to operate the accurate metal forming simulations at the engineering scale with relatively low computational costs, while the CPFEM provided yield stress and plastic potential points to determine the parameters of advanced yield functions.

Here, the mechanical anisotropy of an AA3104-H19 aluminium alloy sheet has been studied by the experiments, full-field CPFEM as well as various phenomenological yield function in 2D and 3D stress space. The full-field CPFEM predictions exhibited a good capability in reproducing experimentally observed mechanical anisotropy of AA3104-H19. Additionally, it is noted that the prediction of 3D Yld2004-18p yield function associated with non-AFR show a good agreement with the CPFEM predicted mechanical anisotropy and experimentally measured cup earing profiles.

Since the multi-level model bridged the CP model and the phenomenological anisotropic plasticity, it would be helpful to quantify the microstructure effects of polycrystal orientations on the texture induced mechanical anisotropy of as-rolled sheet metal. In addition, the CPFEM and anisotropic elasto-plastic constitutive models were presented in a general form so that they can be readily employed for other aluminium alloys in sheet forming applications as well as texture controlling for sheet metal rolling processes.