Mechanical Characterization of the Elastoplastic Response of a C11000-H2 Copper Sheet

This work presents an elastoplastic characterization of a rolled C11000-H2 99.90% pure copper sheet considering the orthotropic non-associated Hill-48 criterion together with a modified Voce hardening law. One of the main features of this material is the necking formation at small strains levels causing the early development of non-homogeneous stress and strain patterns in the tested samples. Due to this fact, a robust inverse calibration approach, based on an experimental–analytical–numerical iterative predictor–corrector methodology, is proposed to obtain the constitutive material parameters. This fitting procedure, which uses tensile test measurements where the strains are obtained via digital image correlation (DIC), consists of three steps aimed at, respectively, determining (a) the parameters of the hardening model, (b) a first prediction of the Hill-48 parameters based on the Lankford coefficients and, (c) corrected parameters of the yield and flow potential functions that minimize the experimental–numerical error of the material response. Finally, this study shows that the mechanical characterization carried out in this context is capable of adequately predicting the behavior of the material in the bulge test.


Introduction
In the sheet metal forming industry, the creation of new parts is common for various devices such as, e.g., containers, pots, radiator components, cooling fins, and conductors. In general, these parts have intricate geometries and/or low dimensional tolerances and, consequently, the design of the forming process is fundamental for manufacturing a good quality component. Sheet forming processes are complex problems due mainly to the nonlinearity produced by friction, the development of large displacements and strains and, in addition, to the difficulty of modeling the anisotropic elastoplastic material behavior. Owing to this, the analysis of these processes is nowadays routinely tackled via numerical simulation using the finite element method [1][2][3][4]. In this context, an accurate elastoplastic modeling of the forming material is required in order to achieve an adequate description of its response during the deformation stages [2].
For this last purpose, a large number of models have been proposed to predict the evolution of hardening and plastic deformation during sheet forming operations. Hardening has been usually Lankford coefficients; and (c) corrected parameters of the yield and flow potential functions that minimize the experimental-numerical error of the material response given in terms of the engineering stress-strain and width-axial true strain curves for the RD, DD, and TD samples. In order to get accurate measurements of the strain field evolution during the tests, the DIC technique is used. Finally, the model parameters obtained with the presented mechanical characterization are assessed in the analysis of the material response in the bulge test by comparing numerical predictions with experimental data of the pressure-dome height curve and final thickness profiles along RD and TD.

Material
The material used in this work is a 0.5 mm thick rolled C11000-H2 99.90% pure copper sheet. As this material is cold rolled, it is assumed to have an orthotropic behavior because the material grains have a preferred orientation after the rolling [35]. In this case, the commonly defined directions are the rolling (RD), diagonal (DD), transverse (TD), and normal (ND), as shown in Figure 1a.

Tensile Test
The tensile test was performed following the ASTM standards [36]. The dimensions of the samples are depicted in Figure 1b, where a gradual reduction of the width (from 12.6 mm to 12.5 mm) that fits the standards has been considered to induce the formation of the necking in the center of the extensometer length. Three samples for each direction, i.e., RD, DD, and TD, were tested. The speed at which the tests were performed was 2 mm/min. As usual, the evolution of the axial forces and extensometer displacements were recorded during the test. In addition, the speckled pattern printed in the samples allowed to track the evolution of the displacement field through images obtained by DIC. The final setup of the test is shown in Figure 2. The following parameters are obtained from the tensile test measurements: Young modulus E, yield strength σ y , ultimate tensile strength σ uts , Lankford coefficients R, strain at the onset of necking ε neck , and fracture strain ε f . In this context, the engineering stress is defined as σ eng = P/A 0 , where P is the instantaneous load during the test and A 0 is the initial transversal area computed as w 0 t 0 , where w 0 is the initial width (at the center of the sample) and t 0 is the initial thickness. Moreover, the axial engineering strain is defined as ε eng = (L − L 0 )/L 0 , where L and L 0 are the current and initial extensometer lengths, respectively.
Assuming volume invariance and homogeneous stress and strain patterns before the necking formation, the true stress (which is equal to the equivalent stressσ for this specific condition), true (logarithmic) strain, and equivalent plastic strain are, respectively, defined as σ = σ eng (1 + ε eng ), According to the works in [37,38], the equivalent stress in terms of the equivalent plastic strain curve should be used to characterize the material hardening response.
To characterize the plastic anisotropy, use is made of the Lankford coefficients which are defined as R θ =˙ε ẇ ε t , where θ is the angle with respect to the rolling direction and the true strain ratesε w andε t are, respectively, related to the width and thickness evolutions of the sample. Once again, assuming an isochoric and homogeneous plastic flow,ε t = −(ε l +ε w ), whereε l is the true strain rate along the axial direction of the sample. The true strain ratesε w andε l are calculated on the basis of the displacement field obtained by DIC at some specific points located within the extensometer length.

Bulge Test
The hydraulic bulge test consists in applying a pressure to a circular section of the sample whose perimeter circumference is clamped. The sample used was a 300 × 300 mm square with the dimensions shown in Figure 3a. The height of the dome and the applied pressure evolutions are measured during the test by means of a comparing dial gauge (resolution of 0.01 mm with a range of 40 mm) and an analog manometer (resolution of 5 bar with a range of 0 to 100 bar), respectively. The test setup involved an image recorder to register the motion of the dial gauges. A manual hydraulic cylinder was used to generate the pressure in the bottom chamber. The set-up is shown in Figure 3b.
The hydraulic bulge test is performed to assess the biaxial response of the sheet. As the material is anisotropic, the classical analytical equations of this test, in which the assumed isotropy leads to an equibiaxial state, are not adequate in this case. Therefore, as shown in Section 3.2, numerical simulations have been carried out in order to more accurately predict the stress and strain patterns that develop during this test.

Constitutive Modelling
To describe the mechanical response of the copper sheet, a large strain plasticity model is used that respectively assumes isotropy and orthotropy for the elastic and plastic regimes. The stress-strain law is given by [37] where σ is the Cauchy stress tensor, C is the elastic constitutive tensor, and e is the Almansi strain tensor with its plastic component e p . The adopted yield criterion is the Hill-48 function defined as where the subscripts of σ denoted by x, y, and z are, respectively, related to RD, TD, and ND shown in Figure 1. Moreover, F, G, H, L, M, and N are the Hill-48 parameters; C p is the hardening function; andσ is the equivalent (i.e., Hill-48) stress. Because of the difficulty to estimate the out-of-plane shear parameters L and M with the available instruments, it is assumed that they are equal to the in-plane shear parameter N. The expressions for each parameter as a function of the yield limits associated to the tensile (with samples cut along RD, DD and TD) and the biaxial (denoted by subscript B) tests are [39]: The hardening function used in this study is a modified Voce law [40] written for RD as where K, σ sat , and n are hardening parameters, andε p is the effective plastic strain whose rate is given below. Moreover, in order to adequately represent the strain corresponding to the onset of necking formation ε neck , Equation (7) is complemented, as proposed in [41], with the following definition of the parameter K, According to the results to be presented in Section 3.1, the following non-associated flow rule is chosen to properly capture the plastic response of the material, where L ν is the well-known Lie (frame-indifferent) derivative,λ is the plastic consistency parameter (computed according to classical concepts of the plasticity theory), and g is the flow potential chosen in this work as the Hill-48 function, but with parameters defined in terms of the Lankford coefficients (obtained from samples tested at RD, DD, and TD) as [42]: Finally, the rate ofε p is given in this context byε p = σ:L ν (e p ) σ [18,19].

Numerical Simulations
Numerical simulations of the tensile and bulge tests are, respectively, used to fit and validate the obtained results with the corresponding experimental measurements. All the computations are performed with the model described in Section 2.2 via an in-house finite element code extensively validated in many forming problems [37]. The finite element meshes are composed of linear hexahedral elements with a B-bar strain-displacement matrix to prevent the volumetric locking effect on the numerical solution when incompressible plastic flows are considered [37]. Moreover, a mesh refinement sensitive analysis was previously carried out in order to ensure mesh-independent results. For both problems, 4 elements along the thickness of the sheet were used.
The boundary conditions and symmetry planes considered in the simulation of the tensile test for samples cut along RD and TD are shown in Figure 4a while those of sample DD are depicted in Figure 4b. The same mesh is used for the RD and TD samples (obviously exchanging both directions for each case) where only one eighth of the geometry is discretized due to the existence of three symmetry planes in these cases. This last condition is not fulfilled in the DD sample for which only one half of the domain is discretized owing to the symmetry plane whose normal is ND. In the three samples, the axial displacement at the boundary is imposed up to the corresponding value of the experimental displacement at the fracture stage for each case. The boundary conditions applied in the simulation of the bulge test are shown in Figure 5. In this case, only one fourth of the domain is discretized as two symmetry planes are considered. While the outer edges of the sample are constrained, an increasing pressure is applied on the bottom surface of the sheet until a dome vertical displacement value similar to that achieved experimentally is reached.

Fitting Procedure
A three-step procedure based on the minimization, via the Levenberg-Marquardt algorithm [43], of the error between the experimental tensile test measurements and the corresponding analytical-numerical predictions is proposed to derive the constitutive material parameters of the constitutive model presented in Section 2.2. The error for each step is quantified by means of the normalized root mean square deviation (NRMSD), which is given by where m is the number of experimental data, y i is the experimental measurement,ŷ i is the analytical-numerical fitted value and ∆ = |y max − y min | is the interval width bounded by the maximum and minimum values of y i .
Hardening parameters derived from the RD tensile curve.
The fitting is carried out using the experimentally measured equivalent stress-equivalent plastic strain curve for the RD tensile sample and the analytical expressions given by Equation (7) together with the constraint imposed by Equation (8). Considering that σ RD and ε neck are obtained from the experimental RD tensile curve, the resulting hardening parameters to be calibrated are σ sat and n.
First prediction of the Hill-48 parameters based on the Lankford coefficients. Since the biaxial yield limit σ B is not available, the yield function parameters given by Equations (3)-(6) cannot be directly computed. Therefore, an associated plastic model is preliminary adopted as a first approximation to describe the material response (it should be mentioned that non-associated model will be used in Step 3 for the final parameter identification). In this context, the sheet orthotropy is evaluated by means of the Lankford coefficients along RD, DD, and TD assuming an isochoric strain field whose in-plane components were obtained from the DIC measurements within the interval of homogeneous strain, i.e., before the onset of necking formation. Then, Equations (11)- (14) are used to compute the initial values of the Hill-48 parameters. The model response using these parameters is assessed in the ε w -ε l (true strains along the axial and width directions of the sample, respectively) and σ eng -ε eng curves for the whole deformation range up to the rupture stage. In this work, the ε w -ε l curve is plotted at the central point of the neck in order to achieve a larger strain interval and, thus, improve the prediction capability of the calibration methodology.
Correction of the parameters of the flow and yield potential functions. In general, the model response considering the Hill-48 parameters analytically obtained in Step 2 does not accurately adjust the corresponding experimental measurements along the whole test in which a triaxial and non-homogeneous stress and strain patterns develop after the early necking formation occurring in the tensile samples. Therefore, these parameters need to be amended. To this end, the proposed correction phase is performed via numerical simulations of the tensile test for the different samples in two successive sub-steps detailed below. This correction process is repeated until the error NRMSD for each curve is minimized.
3.1. Correction of the flow potential function parameters F , G , H , and N using the experimental and computed ε w -ε l curves. Each curve allows the determination of specific parameters keeping the rest constant, i.e., the correction sequence is (a) parameters G and H with the RD curve subjected to the condition G + H = 1 (which is also fulfilled by Equations (12) and (13)), (b) parameter F with the TD curve, and (c) parameter N with the DD curve.
3.2. Correction of the yield function parameters F, G, H, and N using the experimental and computed σ eng -ε eng curves. Once again, each curve allows the determination of specific parameters keeping the rest constant, i.e., the correction sequence is (a) parameters G and H assuming G = G and H = H (it should be noted that these two parameters are not corrected with the RD curve since it has been considered as a reference data to calibrate the hardening response and, therefore, G and F do not play any role provided G + H = 1 is preserved as also stated by Equations (4) and (5), (b) parameter F with the TD curve, and (c) parameter N with the DD curve.

Tensile Test
The average experimental engineering stress-strain curves for the RD, DD, and TD tensile samples are shown in Figure 6, where marked differences in the material response can be seen for the different in-plane directions of the sheet. These data allow obtaining the elastic Young modulus E, the yield strength σ y , the ultimate tensile strength σ uts , the strain at the onset of necking formation ε neck , the strain at the fracture stage ε f , and the Lankford coefficients as shown in Table 1. According to Step 1 of Section 2.4, Figure 7 shows the fitting of the hardening function with the tensile experimental data for the RD sample. In order to guarantee a homogeneous strain pattern in the sample, this procedure is applied within a strain range up to the value of necking formationε p = 6.0%, which is equivalent to the average engineering strain of ε neck = 6.5% (see Table 1). Taking the value of σ RD as the average yield strength of Table 1, the rest of the derived hardening parameters are summarized in Table 2.   The average values of the Lankford coefficients included in Table 1 were computed at the corresponding strain level ε neck assuming, as already mentioned in Step 2 of Section 2.4, an isochoric strain field whose in-plane components were measured by DIC. These Lankford coefficients are used to obtain a first approximation of the Hill-48 parameters of the preliminary adopted associated plastic model. These parameters, summarized in Table 3, are considered in turn to simulate the tensile test from which Figures 8 and 9 show for the three sample directions the experimental and computed ε w -ε l (at the central point of the neck) and σ eng -ε eng curves, respectively. Although Figure 8 exhibits a reasonable agreement between the experimental and numerical results, this is not the case for the DD and TD results plotted in Figure 9. This fact leads, as stated in Section 2.4, to the need to recalculate the Hill-48 parameters. Table 3. Parameters of the Hill-48 yield and flow potential functions.

F/F' G/G' H/H' N/N'
Initial values for f and g  Then, Step 3 of the calibration procedure described in Section 2.4 is applied through numerical simulations of the tensile tests to obtain a more representative description of the material behavior in the plastic regime. The resulting parameters of the yield and flow potential functions are included in Table 3. As shown in Table 4, the numerical predictions obtained with these parameters minimize the experimental-numerical error for the results plotted in Figures 8 and 9. This is, in particular, more apparent for the DD and TD stress-strain curves, respectively, plotted in Figure 9b,c which, as already mentioned, exhibit the largest initial NRMSD values. Table 4. Normalized root mean square deviation (NRMSD) in the ε w -ε l and σ eng -ε eng curves.

Curve Initial NRMSD (Step 2) Final NRMSD (Step 3)
ε w -ε l for RD 0.021 0.020 ε w -ε l for DD 0.021 0.016 ε w -ε l for TD 0.050 0.007 σ eng -ε eng for RD 0.020 0.019 σ eng -ε eng for DD 0.160 0.051 σ eng -ε eng for TD 0.130 0.049 Finally, Figure 10 shows the experimentally measured via DIC and computed axial strain (ε l ) at two levels of axial engineering strain (ε eng ) for the three tensile sample directions (RD, DD, and TD). For these orientations, the neck occurs at very different levels of strain. This is the reason why different scales are defined for each orientation of the specimen. The level 4% is used to visualize strain gradients that appear at low strain levels, such as those exhibited by the DD and TD specimens, where the level 8% is chosen to visualize the same effect in the RD orientation.  Figure 11 shows the average experimental and computed pressure-dome height curve in the bulge test. The two plotted numerical predictions have been respectively obtained at the end of Steps 2 and 3 of the fitting procedure described in Section 2.4, i.e., using the Hill-48 parameters presented in Table 3. For the maximum dome height of the sample, Figure 12 depicts an experimental-numerical comparison of the thickness profiles along the RD and TD radial lines while Figure 13 shows the computed in-plane maximum and minimum principal stresses at the outer surface. These results clearly illustrate the noticeable orthotropic response exhibited by the material at the end of the test.

Discussion
An important aspect that has been taken into account in the proposed mechanical characterization methodology applied to the studied material is, as shown in Figure 6 and Table 1, the early necking formation that causes, more markedly in the DD and TD samples, non-homogeneous stress and strain patterns. Due to this, the hardening and orthotropic parameters cannot be exclusively derived by means of analytical expressions but through numerical simulations able to describe the complex triaxial state that is generated in the post-necking regime. Figure 7 shows an excellent fitting (with NRMSD < 0.01) obtained with the modified Voce law of the material response along the whole deformation range exhibited by the RD sample. This would not be the case if the Voce or Swift laws [10] had been used, for which the hardening would be under or overestimated, respectively. Moreover, it should be noted that this good agreement with the experimental data was achieved thanks to the restriction imposed on the hardening parameter K (see Equation (8)) as, if this condition had not been taken into account, the strain corresponding to the onset of necking formation would have been unrealistically predicted as 30% instead of the experimentally measured value of 6.5%.
The results plotted in Figures 8 and 9, and to a lesser extent the ultimate tensile strength (σ uts ) shown in Table 1, clearly show the orthotropic and non-associated nature of the material plastic response (this is consequently reflected in the final Hill-48 parameters shown in Table 3 that differ from those related to an isotropic and associated response, i.e., F = F = G = G = H = H = 0.5 and N = N = 1.5). They also reveal the soundness of the proposed fitting procedure in properly describing the stress and strain patterns that develop during practically the entire deformation range of the samples along the different studied directions. As the hardening model is adjusted in the RD direction, the numerical results agree well with the experimental measurements in the whole range of strain for this direction. However, some differences between the experimental and numerical results are observed for DD and TD, which are mainly manifested at the higher strain levels. Therefore, it may be possibly necessary to include more deformation paths in the characterization to assess the choice of the yield function to improve the model predictions. Moreover, it should be also noted that the different strain levels at which the neck is formed for the RD, DD, and TD samples (Table 1) is reasonably reproduced by the model.
The strain contours plotted in Figure 10 exhibit an overall good agreement between the experimental and computed results. Although the major discrepancies are related to the maximum values of the highly localized strain field, the extension of the necking zone is, however, notably well captured by the numerical predictions. This confirms the reasonable predictive capabilities achieved by the constitutive model used in this work when its parameters are properly identified.
The performance of the previous mechanical characterization is evaluated in the analysis of the bulge test in which the material is mainly subjected to biaxial loading conditions [7,44]. It is seen that the computed pressure-dome height curve and thickness profiles, respectively, shown in Figures 11  and 12 are in good agreement with the corresponding experimental data. Although no significant differences are observed for the numerical predictions computed with the material parameters derived in Steps 2 and 3 of the fitting procedure described in Section 2.4, a better experimental-numerical agreement is obtained by the results from the non-associate model. Finally, Figure 13 clearly shows the non-equibiaxial stress pattern that develops during the test due to the already mentioned noticeable orthotropy of the material, where, according to the hardening responses shown in Figure 9, the peak values of the maximum and minimum principal stress contours are, respectively, aligned with RD and TD.

Conclusions
An elastoplastic characterization of a rolled C11000-H2 99.90% pure copper sheet considering the orthotropic non-associated Hill-48 criterion together with a modified Voce hardening law has been presented. To this end, a robust inverse calibration approach, based on an experimental-analytical-numerical iterative predictor-corrector methodology, has been proposed and used with tensile test measurements of samples stretched along three in-plane directions. The main difficulty related to the early development of non-homogeneous stress and strain patterns has been successfully tackled with this approach as the calibration procedure carried out in this work was found to adequately describe the orthotropic behavior of the material during the whole deformation range up to the rupture stage. Besides, such mechanical characterization was also experimentally validated in the simulation of the bulge test where sound predictions of the pressure-dome height evolution and final thickness profiles have been achieved.
Finally, future research including more strain paths will be devoted to improve the current material characterization by assessing the effects of out-plane shear response, kinematic hardening and mechanical damage and, in addition, to explore the potential benefits of using more advanced yield criteria. Funding: This research was founded by the Agencia Nacional de Investigación y Desarrollo (ANID) of Chile.

Acknowledgments:
The authors acknowledge with thanks the support through Fondecyt Project 1180591 funded by ANID-Chile. The first author's doctoral studies are financially supported by ANID through a "Beca de Doctorado Nacional" ANID/Doctorado Nacional/2016-21161640.

Conflicts of Interest:
The authors declare no conflict of interest.
Data Availability: The raw/processed data required to reproduce these findings cannot be shared at this time as the data also forms part of an ongoing study.

Abbreviations
The following abbreviations are used in this manuscript.