Characterization of the Elastoplastic Response of Low Zn-Cu-Ti Alloy Sheets Using the CPB-06 Criterion

Unlike other HCP metals such as titanium and magnesium, the behavior of zinc alloys has only been modeled in the literature. For the low Zn-Cu-Ti alloy sheet studied in this work, the anisotropy is clearly seen on the stress-strain curves and Lankford coefficients. These features impose a rigorous characterization and an adequate selection of the constitutive model to obtain an accurate representation of the material behavior in metal forming simulations. To describe the elastoplastic behavior of the alloy, this paper focuses on the material characterization through the application of the advanced Cazacu-Plunket-Barlat 2006 (CPB-06 for short) yield function combined with the well-known Hollomon hardening law. To this end, a two-stage methodology is proposed. Firstly, the material characterization is performed via tensile test measurements on sheet samples cut along the rolling, diagonal and transverse directions in order to fit the parameters involved in the associate CPB-06/Hollomon constitutive model. Secondly, these material parameters are assessed and validated in the simulation of the bulge test using different dies. The results obtained with the CPB-06/Hollomon model show a good agreement with the experimental data reported in the literature. Therefore, it is concluded that this model represents a consistent approach to estimate the behavior of Zn-Cu-Ti sheets under different forming conditions.


Introduction
Zinc is commonly used as a corrosion-resistant coating. However, it is also produced as thin sheets, mainly used in architecture and construction as roofing material, rain gutters and decorative products. In addition to its corrosion resistance property, zinc shows high malleability, ductility and a high quality and durable surface finish. Despite these wide uses, there is a lack of studies and information with respect to zinc sheet formability, in which the high c/a ratio may lead to a marked and evolving anisotropy in the plane of the blank as a consequence of the texture modification [1][2][3][4].
Zinc has a Hexagonal Close Packed structure (HCP), for which the rolling process leads to a strong texture, and slight local changes in the material induced by the manufacturing process (non-homogeneous cooling rates, local microsegregation of alloys, etc.) often generate significant

Material
The material used in this work is the low Zn-Cu-Ti alloy commercially known as Zn-20. The RD, DD and TD tensile samples were gathered from cold-rolled sheets of 1.0 mm thickness tested at a strain rate of 0.007 s −1 [1]. The experimental true stress-strain tensile curves obtained and reported in [1] are presented in Figure 1. strain rate of 0.007 s −1 [1]. The experimental true stress-strain tensile curves obtained and reported in [1] are presented in Figure 1. The mechanical properties, i.e., yield strengths and Lankford coefficients, of the Zn-20 alloy sheet are presented in Table 1. The Young modulus and Poisson ratio with respective values of 127.7 GPa and 0.23 were taken from [30]. These data are used in the fitting procedure to be presented in Section 2.3.

CPB-06/Hollomon Elastoplastic Model
The constitutive model used in this work is defined in the context of the associated flow rule and rate-independent plasticity with the standard elastoplastic strain decomposition [31]. It was assumed that RD and TD are aligned with the x and y axes in the material reference system; thus, the z axis defines the out-of-plane component. The CPB-06 yield criterion adopted to describe the material response is written as [19,20]: The mechanical properties, i.e., yield strengths and Lankford coefficients, of the Zn-20 alloy sheet are presented in Table 1. The Young modulus and Poisson ratio with respective values of 127.7 GPa and 0.23 were taken from [30]. These data are used in the fitting procedure to be presented in Section 2.3.

CPB-06/Hollomon Elastoplastic Model
The constitutive model used in this work is defined in the context of the associated flow rule and rate-independent plasticity with the standard elastoplastic strain decomposition [31]. It was assumed that RD and TD are aligned with the x and y axes in the material reference system; thus, the z axis defines the out-of-plane component. The CPB-06 yield criterion adopted to describe the material response is written as [19,20]: where σ is the equivalent stress, σ is the Cauchy stress tensor, Y is the isotropic hardening stress and ε p is the equivalent plastic strain. The equivalent stress is given by: such that f (χ) , for χ = Σ or χ = γ, is defined as: where a is the degree of homogeneity, Σ i are the principal components of the transformed stress tensor, γ i are the modified anisotropic coefficients and k is the asymmetry parameter (related, as already mentioned, to the SD effect). The reported expression for the transformed stress tensor Σ is given by: where the components of tensor L are the anisotropic coefficients and σ m is the deviatoric part of the Cauchy stress tensor expressed in the material reference system. The modified anisotropic coefficients γ i are: Moreover, the hardening behavior is described by means of the Hollomon power law written for RD as: where K is the strength coefficient, n is the hardening exponent (note that in this context, unlike other approaches [1,7,14,16,17], these two coefficients are only defined for RD), and with ε 0 = σ yp RD K 1 n , σ yp RD being the yield strength for RD (see Table 1). In addition, the rate of the equivalent plastic strain is . ε p = σ m : . ε p σ , such that ε p is the plastic strain tensor whose rate obeys the classical (objective, i.e., frame-indifferent) associated flow rule λ is the plastic consistency parameter. This model was implemented in an in-house finite element code with a radial-return scheme based on the Newton-Raphson iterative method [31]. The proposed model is used to describe different strain path-dependent behaviors in a complete set of the bulge test. The computed numerical results show good agreement with the experiments, as will be discussed in Section 4. Additionally, the present work improves previous referenced studies on Zn-Cu-Ti sheet formability, by allowing the fitting of all directions with the use of a unique set of parameters for both the yield function and Lankford coefficients.

Fitting Procedure via the Tensile Test
The fitting procedure is based on the analytical expression for the stress and strain behavior on the unidirectional tensile test adopting the plane stress assumption. Thus, only σ xx , σ xy and σ yy are different from zero. The steps involved in the methodology, summarized by the flow diagram in Figure 2, are described below.

Data Preparation
The experimental data was considered until the Ultimate Tensile Stress (UTS) in the axial true stress-true strain ( ° − ° ) curves of the samples (i.e., = 0°, 45° and 90° for RD, DD and TD, respectively), for which a homogeneous state is assumed [32]. For simplicity, the same number of experimental ( ° − ° ) values m were considered for the curves of the three samples. To obtain the plastic component of the axial strain , ° for a stress beyond the yield strength, a simple decomposition was used: where is the Young's modulus.

Hardening Fitting
The hardening parameters (K and n) were obtained through the minimization of the following objective function: where ° is the numerical axial stress for the RD sample computed with the expression of ° given in Equation (16).

CPB-06 Fitting
The objective function proposed in [33] is also used here to obtain, through its minimization, the parameters involved in the CPB-06 model. A symmetric material response, i.e., k = 0, is assumed, since there is no experimental evidence of twining for this alloy for low strain rates. Moreover, L11 = 1 was chosen [20,[23][24][25][26][27][28][29], while L55 and L66 were also set to 1 due to the unavailability in this study of experimental results associated with the out-of-plane stress components. In summary, the CPB-06

Data Preparation
The experimental data was considered until the Ultimate Tensile Stress (UTS) in the axial true stress-true strain (σ exp θ • − ε exp θ • ) curves of the θ samples (i.e., θ = 0 • , 45 • and 90 • for RD, DD and TD, respectively), for which a homogeneous state is assumed [32]. For simplicity, the same number of experimental (σ exp θ • − ε exp θ • ) values m were considered for the curves of the three samples. To obtain the plastic component of the axial strain ε exp p,θ • for a stress beyond the yield strength, a simple decomposition was used: where E is the Young's modulus.

Hardening Fitting
The hardening parameters (K and n) were obtained through the minimization of the following objective function: where σ num 0 • is the numerical axial stress for the RD sample computed with the expression of σ num θ • given in Equation (16).

CPB-06 Fitting
The objective function proposed in [33] is also used here to obtain, through its minimization, the parameters involved in the CPB-06 model. A symmetric material response, i.e., k = 0, is assumed, since there is no experimental evidence of twining for this alloy for low strain rates. Moreover, L 11 = 1 was chosen [20,[23][24][25][26][27][28][29], while L 55 and L 66 were also set to 1 due to the unavailability in this study of experimental results associated with the out-of-plane stress components. In summary, the CPB-06 parameters to be obtained are six L coefficients and exponent a. The objective function is written as: where R exp θ • and R num θ • are the experimental and numerical Lankford coefficients of a θ sample, respectively, and W is a weighting factor. For simplicity, the weights W were set to 1 in this work for the five terms of Equation (11).
The fitting routine, to minimize Expressions (10) and (11), is based on the non-linear Levenberg-Marquardt algorithm.
Numerical Stress σ num θ • : The numerical stress term is obtained from the general form of the equivalent stress given by Equation (2), where f (Σ) can be written for a uniaxial tensile loading in the form of: where the expressions for ϕ 1 , ϕ 2 and ϕ 3 are: From (12), it can be seen that for θ • = 0 (RD), f (ϕ) becomes f (γ) and σ (σ) = σ 0 • is fulfilled. In addition, σ num θ • can be written as: where includes the set of L coefficients to be fitted.
Numerical Lankford Coefficients R num θ • : Considering the inherent plastic incompressibility of the model, the numerical Lankford coefficients are written as: where the superscript r denotes the tensile test reference system such that the sample is loaded in the x direction. Although the uniaxial test is very important, it is also relevant to assess the proposed model under loading conditions that are more representative of real applications. For this reason, the bulge test will be used to study the mechanical response of the material under biaxial loading.

Numerical Simulations of the Bulge Test
According to the bulge tests carried out in [1], three different dies with the following minor to major axis ratios β were used: 1.00 (equibiaxial), 0.50 and 0.33. For the β = 0.5 and β = 0.33 dies, samples with the major axis aligned with RD and TD were considered (the largest of the three dies was 120 mm). Therefore, five simulations were performed in order to replicate the experimental strain paths reported in [1].
The complete domain was meshed with three sub-sets: the sheet sample, the die and the sheet contact interface. The die was assumed to be rigid. For the sheet, 10800 trilinear 8-noded hexahedral elements with B-bar integration to avoid numerical locking [31] were used (considering 6 elements along the thickness), while the die and contact interface were discretized with bilinear 4-noded quadrilateral elements, 2160 for the die and 3600 for the interface. The geometrical models and finite element meshes of the bulge test for the different analyzed dies are plotted in Figure 3. along the thickness), while the die and contact interface were discretized with bilinear 4-noded quadrilateral elements, 2160 for the die and 3600 for the interface. The geometrical models and finite element meshes of the bulge test for the different analyzed dies are plotted in Figure 3. As in the experiment, an internal pressure was prescribed on one side of the sheet, with displacements restrained at the edges of the sheet. Coulomb friction is considered with a friction coefficient value of 0.3 between the sheet and the die for all simulations.

Fitting Procedure
The obtained Hollomon and CPB-06 fitted coefficients are respectively presented in Tables 2 and  3.   As in the experiment, an internal pressure was prescribed on one side of the sheet, with displacements restrained at the edges of the sheet. Coulomb friction is considered with a friction coefficient value of 0.3 between the sheet and the die for all simulations.

Fitting Procedure
The obtained Hollomon and CPB-06 fitted coefficients are respectively presented in Tables 2 and 3. The adjusted true stress-strain curves, based on the CPB-06/Hollomon model, are displayed in Figure 4. The numerical Lankford coefficients obtained with the parameters reported in Tables 2 and 3 are summarized in Table 4. The error of the fitting procedure in the true stress-strain curves and Lankford coefficients for the three test directions can be assessed through the Root Mean Square Error (RMSE) given by the expressions: The obtained RMSEs for the true stress-strain curves and Lankford coefficients for the three test directions are shown in Table 5. Table 5. RMSE of the fitting procedure in the true stress-strain curves and Lankford coefficients. The numerical Lankford coefficients obtained with the parameters reported in Tables 2 and 3 are summarized in Table 4. The error of the fitting procedure in the true stress-strain curves and Lankford coefficients for the three test directions can be assessed through the Root Mean Square Error (RMSE) given by the expressions: The obtained RMSEs for the true stress-strain curves and Lankford coefficients for the three test directions are shown in Table 5.  Figure 5 shows the plane stress yield envelope in the σ xx and σ yy plane (with σ xy = 0) at the initiation of yielding for the von Mises, Hill-48 and CPB-06 criteria. The Hill-48 function is computed based on the R values shown in Table 1.

Bulge Test
The experimental and numerical strain paths on the major and minor strains diagram for the different dies and sample orientations are plotted in Figure 6 (the results from the tensile tests are also included for completeness). The experimental measurements correspond to those reported in [1]. The numerical results were gathered from the central element of the external side of the sheet.

Bulge Test
The experimental and numerical strain paths on the major and minor strains diagram for the different dies and sample orientations are plotted in Figure 6 (the results from the tensile tests are also included for completeness). The experimental measurements correspond to those reported in [1]. The numerical results were gathered from the central element of the external side of the sheet.

Discussion
The complexity of the anisotropy shown by zinc alloys requires the use of more elaborate elastoplastic constitutive laws. Thanks to the recent advances of the material science community, we can dispose right now of a rather large amount of the various constitutive laws that can be used to study anisotropic materials like zinc alloys. These tools vary in complexity, and of course in precision. As one might expect, more complex constitutive relationships often convey a better precision, but also, in a general way, it can be said that the more complex the material, the larger the number of material parameters that should be identified [34][35][36]. Needless to say, that large number of material parameters relies on complex and expensive experimental campaigns, which often do not meet the requirements of competitive industries. Additionally, the identification process of these material parameters is carried out through the use of inverse analysis tools, typically leading to ill-posed problems. The challenge consists in obtaining a balance between complexity and precision.
As already mentioned, previous studies on zinc alloy formability carried out material characterization by means of the Hill-48 yield function and the Hollomon or Swift hardening laws, separately fitted for each sample direction [1,5,7,[14][15][16]. Although the Hill-48 criterion is a simple anisotropic plastic model that requires a low number of material parameters, using independent models for the different loading directions drastically increases the number of material parameters. Additionally, the implementation of such an approach in some numerical codes could be cumbersome.
The use of the Hollomon hardening law simplifies the implementation of the constitutive model and the fitting process, showing good agreement between the experimental and numerical stressstrain curves, with less than 4 MPa of RMSE for the RD (see Table 5). It is important to point out that the change in mechanical response in directions other than RD is only driven by the yielding criterion CPB-06. In addition, as an interesting result, using the Hollomon hardening law identified from the RD data combined with the CPB-06 flow rule makes it possible to improve the fit of the hardening curve in the DD and TD. Thus, the RMSE in the DD and the TD is reduced by 50% in comparison to the one in the RD (see Table 5). This improvement can also be seen in a qualitative way when comparing the numerical predictions and experimental data in the stress-strain curve plotted in

Discussion
The complexity of the anisotropy shown by zinc alloys requires the use of more elaborate elastoplastic constitutive laws. Thanks to the recent advances of the material science community, we can dispose right now of a rather large amount of the various constitutive laws that can be used to study anisotropic materials like zinc alloys. These tools vary in complexity, and of course in precision. As one might expect, more complex constitutive relationships often convey a better precision, but also, in a general way, it can be said that the more complex the material, the larger the number of material parameters that should be identified [34][35][36]. Needless to say, that large number of material parameters relies on complex and expensive experimental campaigns, which often do not meet the requirements of competitive industries. Additionally, the identification process of these material parameters is carried out through the use of inverse analysis tools, typically leading to ill-posed problems. The challenge consists in obtaining a balance between complexity and precision.
As already mentioned, previous studies on zinc alloy formability carried out material characterization by means of the Hill-48 yield function and the Hollomon or Swift hardening laws, separately fitted for each sample direction [1,5,7,[14][15][16]. Although the Hill-48 criterion is a simple anisotropic plastic model that requires a low number of material parameters, using independent models for the different loading directions drastically increases the number of material parameters. Additionally, the implementation of such an approach in some numerical codes could be cumbersome.
The use of the Hollomon hardening law simplifies the implementation of the constitutive model and the fitting process, showing good agreement between the experimental and numerical stress-strain curves, with less than 4 MPa of RMSE for the RD (see Table 5). It is important to point out that the change in mechanical response in directions other than RD is only driven by the yielding criterion CPB-06. In addition, as an interesting result, using the Hollomon hardening law identified from the RD data combined with the CPB-06 flow rule makes it possible to improve the fit of the hardening curve in the DD and TD. Thus, the RMSE in the DD and the TD is reduced by 50% in comparison to the one in the RD (see Table 5). This improvement can also be seen in a qualitative way when comparing the numerical predictions and experimental data in the stress-strain curve plotted in Figure 4. The good agreement shown in Figure 4 up to the Ultimate Tensile Stress (UTS) is an encouraging result, since it means that the damaging process of the material could eventually be captured by coupling the presented approach with some coupled non-local damage models [37][38][39].
An important feature of this approach is related to the change of the yielding locus induced by the CPB-06 yielding criterion. Figure 5 shows a comparison between the yielding surface in plane bi-axial stress (no shear) of different classic flow rules. The key features of the proposed approach are obviously the anisotropic nature of the yielding criterion and also the Strength Differential (SD) effect. In preliminary fitting runs, where the k parameter was set to 0 (neglecting SD effect), the error based on stress and Lankford increased in all three directions. especially for TD, where the fitted stress-strain curve was over-estimated for the entire range and the Lankford value decreased from 0.60 to 0.51 in TD. Assuming an asymmetric behavior (presence of SD effect) with a fitting of the k parameter, it is possible, at the same time, to match the stress-strain curves without compromising the estimation of the Lankford coefficients for all three directions (RD, DD and TD).
It is also worth mentioning that the strain paths predicted by the model in the case of uniaxial loading present good agreement with the experimental data. Figure 6 shows the experimental and the numerical predictions of the strain path corresponding to the RD, the DD and the TD uniaxial tests. A more quantitative way of looking at the quality of the prediction in terms of transversal strain is to look at the different Lankford coefficients (R). It can be seen in Table 4 that the predicted R values are in adequate agreement with the corresponding experimental measurements given in Table 1. The relative error in each of the three Lankford coefficients is lower than 1%, which stems from the way they have been defined; thus, the RMSE is close to 0 in all three directions. The definition of the Lankford coefficients as functions of the L, k and a result in values that are almost the same as the ones determined experimentally. The previous results prove that the proposed approach is able to successfully predict the anisotropic mechanical response of the studied zinc alloy over different uniaxial loading directions. Furthermore, these classic and simple experimental tests provide all the information required in order to calibrate the model. However real-life applications involve mechanical loadings that are much more complex. For instance, biaxial loading conditions are common in many material forming industrial processes. The bulge test simulation (Section 3.2) is a very interesting application involving biaxial loading of the material sheet. The strain paths (experimental and numerical predictions) corresponding to the different elliptical dies used in the bulge test are plotted in Figure 6. On top of the aforementioned uniaxial strain paths, the numerical prediction of the equibiaxial loading condition also presents an excellent agreement with experimental measurements (see blue data series in Figure 6). Concerning the bulge tests with elliptical dies, they can be divided into two sets of experiments by using the material direction (RD or TD) oriented with the long axis of the ellipse. For the sake of simplicity, these two sets of bulge tests will be referred to as RD and TD, respectively.
The bulge RD tests present a slight deviation to the right of the experimental cloud point for the 0.5 die, but highly displaced to the left for the 0.33. In the case of the TD, the numerically obtained curves deviated slightly to the left for the 0.33 die but, contrary to the RD situation, were highly displaced for the 0.5 die. The slope of the different strain paths denotes the behavior described above. In particular, a satisfactory experimental validation of the numerical model was obtained for the tensile test, bulge equibiaxial and bulge for paths β = 0.50-RD and β = 0.33-TD, where only the cases β = 0.50-TD and β = 0.33-RD exhibit small differences.

Conclusions
The CPB-06/Hollomon associate constitutive model, in addition with the proposed fitting procedure, proves to be a valid and robust way to describe the elastoplastic anisotropic behavior of the Zn-20 alloy. In this context, a unique set of anisotropic coefficients was able to reproduce the experimental tensile stress-strain curves and Lankford coefficients. The strain paths in the bulge test using different dies were properly validated for the equibiaxial, β = 0.50-RD and β = 0.33-TD cases while only approximate results were obtained for the β = 0.50-TD and β = 0.33-RD cases. These results, together with the good approach of the stress-strain curves, reinforce the use of an associated flow rule to reproduce the anisotropy behavior of Zn-20 sheets. This improvement is closely related with the use of a specific yield function that considers the SD effect, as does the CPB-06. The use of an associated flow rule simplifies the implementation of the constitutive models, gives mathematical and physical consistency to the solution and reduces the complexity of the fitting process because a reduction in the number of coefficients to be defined.
Finally, the present work sets new steps to improve the predictability of more general forming conditions including combined hardening laws and damage criteria. Funding: This work has been supported by the funds provided by the "Comisión Nacional de Investigación Científica y Tecnológica" CONICYT-Chile, through the ECOS-SUD Chile project C16E02 and the Fondecyt Project 1180591.