Numerical research on springback law of tension-compression yield asymmetry sheet in multi-point forming

Abstract

Springback of SUS430/AA1050/TA1 (SAT) laminated sheet in multi-point forming (MPF) is difficult to predict due to its complex microstructure. In this work, the polycrystalline texture and crystal structure of the three materials were measured firstly by electron backscatter diffraction (EBSD). Crystal plasticity finite element method (CPFEM) considering twinning and slip was used to predict the anisotropy and tension-compression yield asymmetry (TCYA) of the materials. The data obtained by virtual simulation are validated and used to calibrate the phenomenological yield criteria. In the springback prediction model, the three materials were modeled using the layered modeling method. The anisotropic yield criteria were implemented into finite element software ABAQUS via a user-defined material subroutine VUMAT to fully consider the anisotropy of the SAT laminated sheet. The comparison between the forming experiment and simulation results shows that the model can predict the springback accurately. Finally, the effect of plastic deformation, strain path, and TA1 layer thickness on springback were discussed.

Introduction

SAT laminated sheet which combines the high strength, wear resistance, corrosion resistance of stainless steel and titanium, the small specific gravity of aluminum, as well as easy processing and high thermal conductivity, is widely applied in shipbuilding industry. ¹ However, it is difficult to predict the springback because of the complex composition of SAT laminated sheet, which makes it difficult to predict the springback, which makes the design of forming tools tougher and costlier. One of the methods to effectively deal with the problems in the design of tools and process optimization for these advanced materials is to apply accurate calculation simulations, which requires the correct description of the materials’ deformation characteristics. ²

Strong in-plane anisotropy, which has a non-negligible effect on springback during the stamping forming process, is found in stainless steel, aluminum, and titanium. In particular, the deformation of titanium with TCYA is more complicated due to its hexagonal close-packed (HCP) structure. ³ Li et al.^4,5 concluded that for FCC and BCC materials, it is sufficient to describe the anisotropy of aluminum and stainless steel based on the linear transformation method and the generalized stress invariance method. Nevertheless, for HCP materials, twinning is required for the lack of slip systems to accommodate plastic deformation, which causes the most significant anisotropy and asymmetry such as Mg-alloy^6–9 and Ti-alloy.^10–13 In earlier research, TCYA was mainly studied with the two-surface plasticity model^14,15 and it is more widely used to incorporate the third stress invariant during the last decade. Lee et al. ¹⁶ introduced a uniaxial analytical method to study the elastic-plastic pure bending and bending-under-tension of AZ31B magnesium alloy sheets. The comparison between the experiments and the analytical method illustrates that the proposed model has the capability to predict the law of springback in Mg-alloy sheets. Nixon et al. ¹⁷ presented a comprehensive tentative and theoretical research on the deformation behavior of $\mathrm{\alpha }$ -titanium under quasi-static conditions. The results of simulations and that of experiment are identical. Barros et al. ¹⁸ described the main features of the implementation of the Cazacu-Plunkett-Barlat2006 (CPB06) yield criterion in the implicit in-house finite element solver DD3IMP. It is reflected in the results that the CPB06 can accurately describe both the r-values and the normalized yield stresses for tensile and compression cases. Li and Fang ¹⁹ conducted a research on the 3D plastic anisotropy and TCYA of extruded Mg-alloy and explored the neutral layer shifts law in three-point bending. The results show that the plastic responses can be described accurately by Yoon2014 in all compression and tension states. Mehrabi et al. ²⁰ developed a novel theoretical approach in view of Cazacu-Barlat2004 (CB2004) yield function under the uniaxial and plane strain cases. Large quantity improvements are displayed in the results with the consideration of the plane strain condition and nonlinear unloading. The above work show that the asymmetric yield criterion based on the third stress invariant is more effective for springback prediction.

Except for experimental purposes (usually involving uniaxial tension and compression in different directions),^21–23 CPFEM is also used as a good example to study the inhomogeneous mechanical response of the materials. Gurao et al. ²⁴ used a viscoplastic self-consistent (VPSC) model with a secant method to corroborate the evolution of microstructure and texture under homeothermy compression of pure-Ti with several different initial values. Knezevic et al. ²⁵ adopted a polycrystal self-consistent model embedded in finite elements to simulate the deformation of textured α-titanium under quasi-static conditions. Lim et al. ²⁶ develop an atomistically informed CPFEM model for $\mathrm{\alpha }$ -Fe, which presents the temperature and strain rate dependent yield surfaces in Fe that diverge from conventional CPFE models. Rodríguez-Galán et al. ²⁷ studied the behavior of nanostructured pure titanium using CPFEM. Becker and Lloyd ²⁸ proposed a plasticity approach for HCP structure. The simplified CP method maintains the characteristic of crystallographic that control macroscale plastic deformation of HCP metals and promote the computational efficiency. Muhammad et al. ²⁹ calculated the wrap bending using CPFEM with a 3D multi-scale modeling framework. The calculated results for texture evolution were consistent with the corresponding experimental results. Kondori et al. ³⁰ presented a reduced-order model of quasi-crystal plasticity whereby the two activation surfaces stand for the slip and twinning-dominated flow. The application of the model for a Mg-alloy plate demonstrates the practicality of the proposed method. Yang et al. ³¹ developed a theoretical approach to provide a basis for producing a whole spectrum of knowledge on the TCYA in plastic deformation of the high-performance materials. The results show that the model is efficient for producing the above findings.

The present work is aimed to establish a model to predict the springback of SAT laminated sheet. First, the layered method was used to model materials. And the in-plane anisotropy of stainless steel and aluminum was described using Yld2004-18p. ³² For titanium, the anisotropic description of the TCYA is performed using the asymmetric yield criterion CB2004 based on the third invariant. CPFEM based on the polycrystalline texture and crystal structure of the three materials is used to predict the anisotropic mechanical properties. ³³ The tensile and compressive yield stresses in various directions were obtained for parameter identification of Yld2004-18p and CB2004. Finally, MPF experiments and simulations of SAT laminated sheet were carried out and compared so that the results of springback prediction were verified.

Constitutive equation

Crystal plasticity

As in previous studies, ³⁴ details of the crystallographic constitutive model for FCC and BCC materials were presented. However, for the HCP crystal structure, the plastic part of the velocity gradient is comprised of contributions from twinning, basal slip, <c+a> pyramidal slip and <a> prismatic slip, which are respectively denoted as:

$${\mathbf{L}}^{\mathbf{p}}={\mathbf{L}}_{\mathbf{tw}}^{\mathbf{p}}+{\mathbf{L}}_{\mathbf{bas}}^{\mathbf{p}}+{\mathbf{L}}_{\mathbf{py}}^{\mathbf{p}}+{\mathbf{L}}_{\mathbf{pr}}^{\mathbf{p}}$$

(1)

where ${\mathbf{L}}_{\mathbf{tw}}^{\mathbf{p}}$ , ${\mathbf{L}}_{\mathbf{bas}}^{\mathbf{p}}$ , ${\mathbf{L}}_{\mathbf{py}}^{\mathbf{p}}$ , and ${\mathbf{L}}_{\mathbf{pr}}^{\mathbf{p}}$ are the plastic part of the velocity gradient for twinning, basal slip, pyramidal slip, and prismatic slip. The velocity gradient associated with these approximations is given as:

$$\begin{array}{c}\multicolumn{1}{c}{{\mathbf{L}}^{\mathbf{p}}=\sum _{\mathit{\alpha }=1}^6{\stackrel{·}{\gamma }}_{\mathrm{tw}}^{\alpha }({\mathit{s}}_{\mathit{tw}}^{\mathit{\alpha }}\otimes {\mathit{m}}_{\mathit{tw}}^{\mathit{\alpha }})+{\stackrel{·}{\gamma }}_{\mathrm{bas}}({\mathit{s}}_{\mathit{bas}}^{\mathit{\alpha }}\otimes {\mathit{m}}_{\mathit{bas}}^{\mathit{\alpha }})}\\ \multicolumn{1}{c}{+\frac{3}{2}\stackrel{.}{\overline{\epsilon }}\frac{\mathbf{\sigma }\prime }{{\sigma }_{\mathrm{vm}}}}\end{array}$$

(2)

where $\stackrel{.}{\overline{\epsilon }}$ is an equivalent plastic strain rate for the pyramidal plus prismatic slip contribution and ${\sigma }_{\mathrm{vm}}$ is the von Mises equivalent stress. ${\mathbf{s}}_{\mathbf{tw}}^{\mathbf{\alpha }}$ and ${\mathit{m}}_{\mathit{tw}}^{\mathit{\alpha }}$ are the Schmidt vector of twinning system. ${\mathit{s}}_{\mathit{bas}}^{\mathit{\alpha }}$ and ${\mathit{m}}_{\mathit{bas}}^{\mathit{\alpha }}$ are the Schmidt vector of slip system. $\mathbf{\sigma }\prime$ is the deviatoric part of the stress which is given by:

$$\mathbf{\sigma }\prime =\mathbf{\sigma }-\frac{1}{3}\mathbf{I}(\mathbf{I}~\mathbf{\sigma })$$

(3)

The stress is based on an unspecified elastic strain energy function Ø, which depends on the elastic part of the deformation gradient, ${\mathbf{F}}^{\mathbf{e}}$ . The Cauchy stress in terms of deformation from the intermediate configuration is:

$$\mathbf{\sigma }\prime =\frac{1}{det\left({\mathit{F}}^{\mathit{e}}\right)}{\mathbf{F}}^{\mathbf{e}}·\frac{\partial Ø}{\partial {\mathbf{E}}^{\mathbf{e}}}$$

(4)

According to prior works, basal <a> slip and twinning, which were assumed to be rate-independent, were modeled respectively using Voce hardening:

$$k_{\mathrm{ss}}^{\mathrm{bas}}=k_{s0}^{\mathrm{bas}}{\left[\frac{\left|{\stackrel{·}{\gamma }}^{\mathrm{bas}}\right|}{{\stackrel{·}{\gamma }}_{S0}^{\mathrm{bas}}}\right]}^{m\prime }$$

(5)

$$k_{\mathrm{ss}}^{\mathrm{tw}}=k_{s0}^{\mathrm{tw}}{\left[\frac{\left|{\stackrel{·}{\gamma }}^{\mathrm{tw}}\right|}{{\stackrel{·}{\gamma }}_{S0}^{\mathrm{tw}}}\right]}^{m\prime }$$

(6)

where $k_{\mathrm{ss}}^{\mathrm{bas}}$ and $k_{\mathrm{ss}}^{\mathrm{tw}}$ are respectively the flow strengths for basal slip and twinning, $k_{s0}^{\mathrm{bas}}$ and $k_{s0}^{\mathrm{tw}}$ are the initial flow strengths. ${\stackrel{·}{\gamma }}^{\mathrm{bas}}$ and ${\stackrel{·}{\gamma }}^{\mathrm{tw}}$ are the shear strain rate for basal slip and twinning. ${\stackrel{·}{\gamma }}_{S0}^{\mathrm{bas}}$ and ${\stackrel{·}{\gamma }}_{S0}^{\mathrm{tw}}$ are the initial shear strain rate. $m\prime$ is a temperature-related correction factor.

Yield function

In this research, the CB2004 yield function is applied to capture the tension-compression asymmetry, which is given by:

$$f={\left(J_2^{\frac{3}{2}}-c{J}_3\right)}^{\frac{1}{3}}={\tau }_y$$

(7)

where, ${\tau }_y$ is the yield stress in pure shear, c is a material parameter limited to a given numerical range $\mathrm{c}ϵ\left[\begin{array}{cc}-3.375\hfill & 2.25\hfill \end{array}\right]$ , which is expressed as:

$$\mathrm{c}=\frac{3\sqrt{3}({\sigma }_t^3-{\sigma }_c^3)}{2({\sigma }_t^3+{\sigma }_c^3)}$$

(8)

$J_2$ and $J_3$ are respectively the second and third invariants of the stress deviator tensor, which are given by:

$$\begin{array}{c}\multicolumn{1}{c}{J_2=\frac{a_1}{6}{({\sigma }_{\mathrm{xx}}-{\sigma }_{\mathrm{yy}})}^2+\frac{a_2}{6}{({\sigma }_{\mathrm{xx}}-{\sigma }_{\mathrm{zz}})}^2}\\ \multicolumn{1}{c}{+\frac{a_1}{6}{({\sigma }_{\mathrm{yy}}-{\sigma }_{\mathrm{zz}})}^2+a_4{\sigma }_{\mathrm{xy}}^2+a_5{\sigma }_{\mathrm{xz}}^2+a_6{\sigma }_{\mathrm{yz}}^2}\end{array}$$

(9)

$$\begin{array}{c}\multicolumn{1}{c}{J_3=\frac{1}{27}(b_1+b_2){\sigma }_{\mathrm{xx}}^3+\frac{1}{27}(b_3+b_4){\sigma }_{\mathrm{yy}}^3}\\ \multicolumn{1}{c}{+\frac{1}{27}[2(b_1+b_4)-b_2-b_3]{\sigma }_{\mathrm{zz}}^3+2b_{11}{\sigma }_{\mathrm{xy}}{\sigma }_{\mathrm{xz}}{\sigma }_{\mathrm{yz}}}\\ \multicolumn{1}{c}{+\frac{2}{9}(b_1+b_2){\sigma }_{\mathrm{xx}}{\sigma }_{\mathrm{yy}}{\sigma }_{\mathrm{zz}}-\frac{1}{9}(b_1{\sigma }_{\mathrm{yy}}+b_2{\sigma }_{\mathrm{zz}})}\\ \multicolumn{1}{c}{{\sigma }_{\mathrm{xx}}^2-\frac{1}{9}(b_3{\sigma }_{\mathrm{zz}}+b_4{\sigma }_{\mathrm{xx}}){\sigma }_{\mathrm{yy}}^2}\\ \multicolumn{1}{c}{-\frac{1}{9}[(b_1-b_2+b_4){\sigma }_{\mathrm{xx}}+(b_1+b_3+b_4){\sigma }_{\mathrm{yy}}]}\\ \multicolumn{1}{c}{{\sigma }_{\mathrm{zz}}^2-\frac{{\sigma }_{\mathrm{yz}}^2}{3}[(b_6+b_7){\sigma }_{\mathrm{xx}}-b_6{\sigma }_{\mathrm{yy}}-b_7{\sigma }_{\mathrm{zz}}]}\\ \multicolumn{1}{c}{-\frac{{\sigma }_{\mathrm{xz}}^2}{3}[2b_9{\sigma }_{\mathrm{yy}}-b_8{\sigma }_{\mathrm{zz}}-(2b_9-b_8){\sigma }_{\mathrm{xx}}]}\\ \multicolumn{1}{c}{-\frac{{\sigma }_{\mathrm{xy}}^2}{3}[2b_{10}{\sigma }_{\mathrm{zz}}-b_5{\sigma }_{\mathrm{yy}}-(2b_{10}-b_5){\sigma }_{\mathrm{xx}}]}\end{array}$$

(10)

where x, y, and z represent RD, TD, and ND respectively. All the coefficients $a_k$ (k = 1,…6) and $b_k$ (k = 1,…11) reduce to unity (c = 0) for isotropic conditions. In particular, the section of the yield locus with ${\sigma }_3=0$ is:

$$\begin{array}{c}{[\frac{1}{6}(a_1+a_3){\sigma }_1^2-\frac{a_1}{3}{\sigma }_1{\sigma }_2+\frac{1}{6}(a_1+a_2){\sigma }_2^2]}^{3/2}\hfill \\ -c[\frac{1}{27}(b_1+b_2){\sigma }_1^3+\frac{1}{27}(b_3+b_4){\sigma }_2^3\hfill \\ \multicolumn{1}{c}{-\frac{1}{9}(b_1{\sigma }_1+b_4{\sigma }_2){\sigma }_1{\sigma }_2]={\tau }_y^3}\end{array}$$

(11)

In the case of a sheet, the anisotropy coefficients of the CB2004 yield criterion are fitted based on the experimental yield stress and the single r-value available using the Levenberg-Marquardt (LM) algorithm. The error function was defined as:

$$E\left(a_k,b_k,c\right)={\displaystyle \sum }_p{w}_p{\left({\sigma }_p^{pr}-{\sigma }_p^{ex}\right)}^2-{\displaystyle \sum }_q{w}_q{\left(R_q^{pr}-R_q^{ex}\right)}^2$$

(12)

where $\mathrm{p}$ and $\mathrm{q}$ is the number of available experimental flow stress and r-values. The yield stress in tension ( ${\sigma }_R^T$ ) and compression ( ${\sigma }_R^C$ ) along the RD follow that:

$${\sigma }_R^T={\tau }_y{\{{\left[\frac{1}{6}(a_1+a_3)\right]}^{3/2}-\frac{c}{27}(b_1+b_2)\}}^{-1/3}$$

(13)

$${\sigma }_R^C={\tau }_y{\{{\left[\frac{1}{6}(a_1+a_3)\right]}^{3/2}+\frac{c}{27}(b_1+b_2)\}}^{-1/3}$$

(14)

The yield stress in tension ( ${\sigma }_T^T$ ) and compression ( ${\sigma }_T^C$ ) along the TD follow that:

$${\sigma }_T^T={\tau }_y{\{{\left[\frac{1}{6}(a_1+a_2)\right]}^{3/2}-\frac{c}{27}(b_3+b_4)\}}^{-1/3}$$

(15)

$${\sigma }_T^C={\tau }_y{\{{\left[\frac{1}{6}(a_1+a_2)\right]}^{3/2}+\frac{c}{27}(b_3+b_4)\}}^{-1/3}$$

(16)

Yielding under equibiaxial tension occurs when ${\sigma }_1={\sigma }_2={\sigma }_b^T$ :

$${\sigma }_b^T={\tau }_y{\{{\left[\frac{1}{6}(a_3+a_2)\right]}^{3/2}-\frac{c}{27}(-2b_1+b_2+b_3-2b_4)\}}^{-1/3}$$

(17)

while yielding under equibiaxial compression occurs when ${\sigma }_1={\sigma }_2={\sigma }_b^C$ :

$${\sigma }_b^C={\tau }_y{\{{\left[\frac{1}{6}(a_3+a_2)\right]}^{3/2}+\frac{c}{27}(-2b_1+b_2+b_3-2b_4)\}}^{-1/3}$$

(18)

The plastic anisotropy parameter $R_{\theta }$ in equation (12) is calculated as:

$$R_{\theta }=-\frac{\frac{\partial f}{\partial {\sigma }_{\mathrm{xx}}}{\sin }^2\theta +\frac{\partial f}{\partial {\sigma }_{\mathrm{yy}}}{\cos }^2\theta -\frac{\partial f}{\partial {\sigma }_{\mathrm{xy}}}\sin 2\theta }{(\frac{\partial f}{\partial {\sigma }_{\mathrm{xx}}}+\frac{\partial f}{\partial {\sigma }_{\mathrm{yy}}})}$$

(19)

Hereinafter the CB2004 anisotropic criterion will be applied to describe the anisotropy and asymmetry of the TA1. The Yld2004-18p anisotropic criterion will be applied to SUS430 and AA1050 sheet which present the obvious anisotropy and the neglected asymmetry.

Determination of phenomenological constitutive model using CPFEM

Texture based virtual mechanical tests

The crystallographic texture of AA1050, SUS430, and TA1 were measured by EBSD to construct the representative volume elements (RVEs). The IPF of TA1 obtained by EBSD and RVEs are shown in Figure 1. ³ The pole figure of TA1 is revealed in Figure 2. The IPF and ODF of AA1050 and SUS430 are respectively displayed in Figures 3 to 6. As identified in Section 2, the dominant mechanisms for TA1 include five gliding systems {0001} <11–20> (Bas<A>), {−1100} <11–20> (Pr<A>), {10–11} <11–20> (Py<A>), {10–11} <11–23> (Py<C+A>), and {11–22}<−1–123> (Py2<C+A>), two twinning modes tensile twinning {10–12} <10–1–1> (Ttw) and compressive twinning {11–22} <11–2–3> (Ctw). Moreover, it shows in Figures 7 and 8 the stress-strain curves for the three materials obtained by CPFEM and tensile test. The parameters used in the simulation are listed in Tables 1 and 2. The CPFEM which was consistent with the tensile test can be applied to determine the parameters of CB2004 and Yld2004-18p.

Figure 1.

Comparison of the IPF between (a) EBSD and (b) RVE of TA1.

Figure 2.

Pole figure of TA1.

Figure 3.

Comparison of the IPF between (a) EBSD and (b) RVE of AA1050.

Figure 4.

ODF of AA1050.

Figure 5.

Comparison of the IPF between (a) EBSD and (b) RVE of SUS430.

Figure 6.

ODF of SUS430.

Figure 7.

Comparison between CPFEM results and experimental data corresponding to the uniaxial loading in RD for SUS430 and AA1050 and their RVEs.

Figure 8.

Comparison between CPFEM results and experimental data corresponding to the uniaxial tension (T) and compression (C) in RD, TD, and ID for TA1 and its RVE.

Table 1.

Elastic stiffness parameters utilized in CPFEM (unit: GPa).

	$C_{11}$	$C_{12}$	$C_{13}$	$C_{44}$	$C_{55}$	$C_{66}$
SUS430	335.2	176.3	–	43.8	–	–
AA1050	108.2	61.3	–	28.5	–	–
TA1	162.4	92.0	69.0	180.7	117.0	35.2

Table 2.

Voce hardening parameters for CPFEM of the materials.

Slip/twin mode	${\stackrel{·}{\gamma }}_0$ (s−1)	$k_{s0}$ (MPa)	$k_{\mathrm{ss}}$ (MPa)	$h_0$ (MPa)	${\stackrel{·}{\gamma }}_{S0}$ (s−1)	$m^{\prime }$	m	${\gamma }_{\mathrm{twin}}$	$V_{\mathrm{th}1}$	$V_{\mathrm{th}2}$
SUS Slip&Tw	1	88.2	132.1	110.3	5.0e10	0.0005	0.05	0.577	0.25	0.85
AA Slip&Tw	1	55.2	89.7	63.2	5.0e10	0.0004	0.05	0.707	0.1	0.5
Bas <A>	0.001	200	250	19.6	5.0e10	0.0001	0.05	–	–	–
Pr <A>	0.001	80	100	37.5	5.0e10	0.0001	0.05	–	–	–
Py <A>	1	750	937.5	320	5.0e10	0.005	0.05	–	–	–
Py <C+A>	1	670	850	320	5.0e10	0.005	0.05	–	–	–
Py2 <C+A>	0.001	690	845	18.8	5.0e10	0.0001	0.05	–	–	–
Ttw	1	200	320	245	5.0e10	0.005	0.05	0.174	0.2	0.5
Ctw	1	380	530	245	5.0e10	0.005	0.05	0.219	0.2	0.5

Determination of the yield surfaces with CPFEM results

The yield stress and r-values in different loading cases predicted by CPFEM for the three materials are shown in Tables 3 and 4. The parameters for Yld2004-18p and CB2004 yield functions were summarized in Tables 5 and 6, respectively.

Table 3.

Predicted yield stress ( ${\mathrm{\sigma }}_{\phi }/{\mathrm{\sigma }}_0$ ) in various loading cases.

	0	15	30	45	60	75	90	Biaxial	Shear
SUS430	1	1.0139	1.0482	1.0821	1.0946	1.0886	1.0833	1.1376	0.5792
AA1050	1	0.9627	0.9261	0.9438	0.9761	0.9451	0.9125	1.1084	0.5645
TA1-T	1			1.0953			1.179	1.1690	0.5519
TA1-C	0.826			1.0142			1.138	1.0090

Table 4.

Predicted r-values in various loading cases.

	0	15	30	45	60	75	90
SUS430	1.0802	1.0178	1.0131	1.1201	1.3116	1.5864	1.8005
AA1050	0.6280	0.8051	0.8905	1.1028	1.5351	0.9267	0.7779
TA1-T	1.47			1.3956			2.45
TA1-C	0.597			1.0864			2.1753

Table 5.

Yld2004-18p coefficients for AA1050 and SUS430.

	c12′	c13′	c21′	c23′	c31′	c32′	c44′
SUS430	0.9570	0.9682	0.9646	0.7891	0.8629	0.9280	0.9075
AA1050	0.9813	0.7570	1.0230	0.9820	1.1797	0.7503	0.6948
a	c12″	c13″	c21″	c23″	c31″	c32″	c44″
6	0.9360	0.9382	0.9902	0.8657	0.9998	0.7886	0.9076
8	0.8934	0.9291	0.8726	0.8889	0.4855	1.1415	1.3806

Table 6.

Identified parameters of CB2004 yield function.

$\overline{\epsilon }$	$a_1$	$a_2$	$a_3$	$b_1$	$b_2$	$b_3$	$b_4$	$c$
0.002	0.8555	0.2758	1.0306	0.3038	3.1502	0.3834	0.6882	0.3846
0.05	1.1243	0.4501	1.2884	1.2749	2.6510	1.0944	0.2817	0.3611
0.1	1.1133	0.7317	1.1455	1.4404	3.0904	4.2995	0.2313	0.0570
0.2	1.1212	0.8850	1.0042	2.2249	1.2540	1.3849	0.4140	0.1357
0.3	1.1411	0.9847	0.8806	1.5503	1.5800	0.1914	0.1618	0.2373

Figure 9(a) and (b) show the comparison of the r-values and normalized uniaxial yield stress calculated by yield functions with those obtained by CPFEM. Both Yld2004-18p and CB2004 functions are found to provide accurate evaluations of the r-values in tensile and compression deformation area. In the figures, the curves show that these models do well in capturing the anisotropic distribution of both r-values and flow stresses of the materials. It is noteworthy that the normalized stress for SUS430, TA1-T, and TA1-C are almost share the same trend with their r-values while AA1050 is otherwise, which is because that the normalized stress and r-values are mainly affected by texture. It can be seen from the ODF and pole figure that the main texture for SUS430 is the orientation with Euler angles (0,30,45) while for TA1, it is the basal texture. These textures are relatively stable during deformation and therefore exhibits regularity values versus with angle of rolling direction. Whereas for AA1050, the main texture is Cube (0,0,0), which is prone to rotation during deformation, thus leading to irregularities in the normalized stress and r-values versus with angle of rolling direction. The yield surfaces of the three materials for different plastic deformation are displayed in Figure 10. As can be seen from Figure 10(a) and (b), the yield surfaces constructed by Yld2004 has good convexity, which consequently could be used to represent the material behavior of AA1050 and SUS430. Therefore, it can be concluded that the second invariant-based Yld2004-18p function is sufficient to describe the regular yielding evolutions of AA1050 and SUS430. In Figure 10(c), CB2004 function is shown to provide an accurate prediction on yield stresses in almost all tension and compression cases. To conclude, CB2004 function is capable of successfully capturing the evolving anisotropy and asymmetry by involving the third invariant.

Figure 9.

(a) Uniaxial yield stress and (b) r-value for the three materials calculated by the yield functions and those predicted by CPFEM.

Figure 10.

Yield surfaces of (a) AA1050, (b) SUS430, and (c) TA1 for different strain level.

Flow stress

In this work, Ludwig equations are selected to illustrates the size of the yield function during deformation. The hardening models of equivalent stress ( $\overline{\mathbf{\sigma }}$ ) and equivalent plastic strain ( $\overline{\epsilon }$ ) of the three materials are presented as follows:

• SUS430: $\overline{\mathbf{\sigma }}=272.73+539.8$ · ${\overline{\epsilon }}^{0.4781}$

• AA1050: $\overline{\mathbf{\sigma }}=29.32+85.25$ · ${\overline{\epsilon }}^{0.4445}$

• TA1: $\overline{\mathbf{\sigma }}\prime =175.76+291.2$ · ${\overline{\epsilon }}^{0.4421}$

FEM and experiment of MPF process

Verification of SAT constitutive model

To verify the accuracy of the hierarchical modeling method for describing the constitutive model of SAT laminated sheet, tensile test simulations and experiments were performed along RD, ID, and TD. In the simulations, the elastic modulus E and the poisson’s ratio ν were set, which are E = 69.3 GPa and ν = 0.33 for AA1050, E = 108.2 GPa and ν = 0.28 for TA1, while E = 215 GPa and ν = 0.34 for SUS430. The CB2004 and Yld2004 yield functions and the hardening rules were implemented in VUMAT subroutines. The stress-strain curves obtained from FEM were compared with the experiments in Figure 11. And the stress in the simulation is taken as the average stress of the elements over the parallel length. The agreement between these sets of data is so excellent that proves the accuracy of the hierarchical modeling method. In the experiment, the bonding interface is not as good as that in the simulation, which results in the amount of damage in the experiment at the initial stage of tensile test.

Figure 11.

Stress-strain curves of SAT laminated sheet in RD, ID, and TD obtained from simulation and experiment.

MPF modeling

ABAQUS/explicit was used to perform the FE simulation of the loading process. As shown in Figure 12, the model consists of multi-point die (MPD), elastic cushion, and sheet. As the effect of bonding interface on the springback prediction is ignorable, the bonding interface of the sheet was assumed perfect. The dimension of the sheet and the elastic cushion was 200 mm × 200 mm (length and width), with a thickness of 10 mm for the elastic cushion while that of the sheet was 2.5 mm, where SUS430 is 0.5 mm, AA1050 is 1.7 mm, and TA1 is 0.3 mm. The friction coefficient was set as 0.15. The element type, loading speed, and unloading process can refer to the previous work. ³⁴ The element size of the elastic cushion is 2 mm × 2 mm × 2 mm, while that of the sheet is 1 mm × 1 mm (length and width). Additionally, both TA1 and SUS430 were divided into one layer of mesh in the thickness direction, in which AA1050 was divided into three-layer grids. The total number of elements for TA1, SUS430, and AA1050 are respectively 40,000, 40,000, and 120,000.

Figure 12.

FE model of MPF.

MPF experiment

The MPF experiments were performed by the YAM-200 MPF press. The dimension of the sheet is identical to that of the FEM. The die surface is a spherical surface of R500 mm as shown in Figure 13(a). During forming, the upper die punch unit is driven by the press downwards, so that the sheet material is gradually bent to fully fit the mold surface. The elastic pad was placed between the die and the sheet. In this work, two experiments were conducted separately when the TA1 layer is put on the upper end (TA1-U) as well as the lower end (TA1-D). The formed parts are shown in Figure 13(c) and (d). After forming, the workpiece was measured by the NDI Procam 3000 3D laser scanner, which is shown in Figure 13(b).

Figure 13.

(a) MPD, (b) the 3D laser scanner, (c) the formed TA1-D, and (d) TA1-U part.

Results and discussion

Comparison and verification of TA1-D and TA1-U workpiece

In Figure 14, it shows the displacement distribution in the z-direction of TA1-D and TA1-U part during springback. And the two parts are respectively simulated for isotropic (ISO) and anisotropy (ANI). The shadow in the figure is the outline of the part before springback. The greater the springback, the greater the displacement in the z-direction. It can be seen from the figure that during the springback, the center position of the sheet moves in the negative direction of the z-axis, while the corners move oppositely, which causes the decrease in the curvature and sagitta of the sheet. The magnitude of the displacement legend in the figure can indicates the decrease in sagitta. It is obvious that the springback of TA1-U is greater than that of TA1-D, which is possibly caused by the movement of the neutral layer. Meanwhile, since the elastic strain on the tension side is easier to release than the compression side, the springback effect of the overall laminate is asymmetrical. When SUS430 is under tensile stress, a greater amount of springback can be released, which is because SUS430 has a larger elastic modulus and is more prone to elastic recovery. For both TA1-D and TA1-U part, the springback predicted by ANI model is larger than that by ISO model. For ISO model, the magnitude of the displacement legend of TA1-D part is 3.55 mm while that of the TA1-U part is 4.97 mm. The difference between the two is 1.42 mm. And for ANI model, the magnitude of the displacement legend of TA1-D part is 4.86 mm while that of the TA1-U part is 5.61 mm. The difference between the two is 0.85 mm, which is less than that of ISO model. It may be due to the TCYA of TA1, whose springback under biaxial tensile stress is greater than that of the biaxial compression for R500 mm.

Figure 14.

Displacement distribution of TA1-D and TA1-U part in the z-direction during springback.

It compares in Figure 15 the profiles of the formed TA1-D and TA1-U parts between the experiments and the FEM results after springback. The maximum error of sagitta for ISO model is less than 0.292 mm for TA1-D part and does not exceed 0.253 mm for TA1-U part. For ANI model, the maximum error of sagitta is less than 0.200 mm for TA1-D part and 0.187 mm for TA1-U part, which respectively are 68.5% and 73.9% of that of ISO model. As can be seen, for both TA1-D and TA1-U parts, the consistency of the ANI model is better than that of the ISO model. Therefore, the modeling method proposed in this work is proved to be able to effectively predict the springback of SAT laminated sheet.

Figure 15.

Profile comparison of the formed TA1-D and TA1-U part between experiments and FEM results after springback.

Effect of bending radius on springback

In general, smaller bending radius will result in larger plastic deformation. Figure 16 presents the displacement in the z-direction with varied bending radii during springback. The bending radii in the x-direction and that in the y-direction are identical. As displayed in Figure 16, the displacement versus forming radius apparently indicates the springback rise on the decreasing of plastic deformation. In addition, the springback of TA1-U is greater than that of TA1-D for all bending radii.

Figure 16.

Displacement for different bending radius.

To study the effect of strain path on springback, inconsistent bending radii in the x-direction and the y-direction are also simulated. Two cases are set respectively. One is varied in y-direction when the bending radius fixed at R500 in the x-direction (x500), while the other is varied in x-direction when the radius fixed at R500 in the y-direction (y500). Figure 17 shows the displacement in the z-direction with different bending radii during springback. It can be seen that the springback of y500 is greater than that of x500 for TA1-U, while the case is contrary for TA1-D. This difference is mainly due to the TCYA on the yield surface. Thus, it could be concluded that TCYA has the different effects on springback under various strain paths. It is significant to consider TCYA in FEM.

Figure 17.

Displacement for bending radius varied in x- and y-direction.

Effect of layer thickness on springback

The corrosion resistance of SAT laminated sheet increases with the thickness of the TA1 layer. If the total thickness of the SAT sheet and the SUS430 layer stay the same, the thickness of the TA1 layer is inversely correlated with that of the AA1050 layer. Therefore, it is necessary to study the effect of different thickness of TA1 and AA1050 layer thickness on springback. In Figure 18, it shows the displacement in the z-direction during springback with varied thickness of AA1050 and TA1 layer. It can be seen from the figure that as the thickness of AA1050 layer increases, the springback becomes smaller for both TA1-U and TA1-D part. When the thickness of the TA1 layer is greater than 0.3 mm, the springback decreases as the thickness of the TA1 layer increases for both TA1-U and TA1-D part. However, when the thickness of the TA1 layer is less than 0.3 mm, the springback of TA1-U and TA1-D show a rising trend with the increase of the thickness of the TA1 layer. The red column in Figure 18 presents the displacement in the z-direction during springback of TA1 layer varied between 0.1 and 0.5 mm. The total thickness of the SAT sheet is constant at 2.5 mm while that of the SUS430 layer is 0.5 mm. The springback becomes larger with the increasing thickness of TA1-U. Nevertheless, the springback displacement decreases when the thickness of TA1 exceeds 0.4 mm, which is because the springback restraint effect of the increase in the thickness of the TA1 layer gradually exceeds the promotion effect of the decrease of the AA1050 layer.

Figure 18.

Displacement in the z-direction during springback for varied metal thickness: (a) TA1-U part and (b) TA1-D part.

Conclusions

A new FE model based on CB2004 and Yld2004-18p yield functions has been developed to investigate the bending and the springback behaviors of SAT laminated sheet. CPFEM was verified to capture the anisotropy of the three materials. According to the current study, the following conclusions can be drawn:

1. CPFEM with consideration of twinning and slip can be used to predict the anisotropy of FCC, BCC materials and the asymmetry of deformation behavior of HCP materials. Therefore, it can be considered as an alternative to uniaxial tension, compression, and biaxial stress experiments to calibrate the orthogonal and non-orthogonal anisotropy yield criterion.

2. By the layered method based on phenomenological yield functions to model the bending process of SAT laminated sheet, the abnormal behavior of springback caused by material anisotropy can be captured. And it is found from the comparison of experimental and simulation results that the springback predicted by the proposed model is more accurately than that by the isotropy model.

3. The springback of laminate sheet is significantly influenced by the elastic modulus, the bending radius, and the thickness of the metal layer. When the metal layer with larger elastic modulus is under tensile stress, the springback is more obvious. In addition, the smaller the bending radius and the thicker the metal, the smaller the springback.