An Extended Drucker Yield Criterion to Consider Tension–Compression Asymmetry and Anisotropy on Metallic Materials: Modeling and Verification

: Pressure sensitive asymmetric Drucker yield criterion is developed to deal with pressure dependent sheet metals for instance steels and aluminum alloys. The sensitivity to pressure is conserved by introducing three-dimensional anisotropic parameters in the first stress invariant; while the third deviatoric stress invariant is remained in odd function form to consider the strength differential effect (SDE). To describe the flow stress directionalities of metallic materials, the Drucker yield function is extended using two transformation matrix consisting of anisotropic parameters. The proposed Drucker yield criterion is utilized to predict the anisotropic yield and plastic deformation of aluminum alloys with weak SDE: AA 2090-T3 with face-centered cubic (FCC) crystal systems and AA 2008-T4 with body-centered cubic (BCC) crystal systems as well as metals with strong SDE: Zirconium clock-rolled plate with hexagonal close packing (HCP) crystal systems. The comparison between the predicted anisotropic behavior and experimental results reveals that the extended anisotropic Drucker yield criterion can precisely model the anisotropy for FCC, BCC and HCP metals. The proposed function is implemented into ABAQUS VUMAT subroutines to describe the four-point bending test which is used to consider the effect of various yield functions and material orientations on deformation behavior. The obtained contours of the cross-section, strain components distribution and also the shift of neutral layer indicate that the extended Drucker yield function can well predict the final geometric configuration of the deformed Zirconium beam


Introduction
With the increasing demand for lightweight, fuel consumption and crashworthiness of automobiles, more and more lightweight materials such as aluminum alloys are increasingly demanded in the car industry to produce body parts in view of their excellent mechanical properties. Besides, texture patterns produced during the rolling process lead to vastly different properties in the transverse and rolling directions. In order to simulate the sheet metal forming process as precise as possible, a precisely calibrated plastic constitutive is needed to describe the mechanical properties as well as deformation behavior under various stress states. Therefore, it is required to carefully study the anisotropic yielding evolution of these sheets under various forming scenarios. In 1948, Hill [1] put forward his first phenomenological yield function which was quadratic. In 1979, Hill [2] proposed a generalized anisotropic yield function which was non-quadratic and applied to non-steel materials. Similar to Hill's 48, Hosford [3] proposed generalization for materials exhibiting orthotropic symmetry using Hershey's isotropic criterion [4]. For the sake of involving shear components in the yield functions, Barlat et al. [5] modified Hosford's criterion by introducing shear stresses in the formulation. Barlat et al. [6] and Karafillis et al. [7] introduced the linear transformations on the stress component to consider the effect of material deformation directionality. Consequently, Barlat proposed a well-known plane stress yield criteria [8] expressed in form of principal values of linear transformation tensors and an anisotropic one for spacial 3D stress conditions [9] using two transformations of stress deviator. Banabic et al. [10] independently developed BBC family criteria since 2000 and the latest version, BBC2005 criterion, was, in fact, the same as that of Yld2000-2d. Bron and Besson et al. [11] utilized two stress tensor transformations to substitute in expression the Karafillis and Boyce yield function [7], introducing twelve anisotropy coefficients to describe anisotropic materials. Another alternative method to extend the existing isotropic yield function was done by Cazacu et al. [12] who generalized the second and third invariants of deviatoric stress tensor of Drucker's 3D yield criterion. Therefore, the isotropic yield criterion was transformed into anisotropic yield criterion by using anisotropic generalizations instead of stress deviator invariants. Gao et al. [13] proposed a yield function considering three stress invariants. Yoshida et al. [14] developed a high order polynomial based spacial yield criterion according to the work of Cazacu and Barlat [12] and assessed the performance of the criterion by comparing the predicted stress directionalities and r-values with experimental observations. Lou et al. [15] proposed their anisotropic yield function by adding up n-components of Drucker's function and validated it by simulation and experiment.
Lou et al. [16] articulated that during in-plane deformation, the average tensile yield stress of magnesium alloy AZ31B was twice of the compressive yield stress approximately which was caused due to strength differential effect (SDE). Metals, especially hexagonal close-packed (HCP) crystal structures, undergo plastic deformation through not only slip and but also twinning. SDE is induced mainly by directional twinning: one-directional shear leading to twinning formation, but not vice versa. Cazacu and Barlat [17] put forward anisotropic yield criterion expressed by the second and third invariants of the stress deviator and later extended it using the generalized invariants to describe SDE as well as the anisotropy of magnesium and its alloys. Nixon et al. [18] used the strategy of stress tensor transformation to extend the isotropic yield function of Cazacu et al. [18] to orthotropic function and it was recommended to use in case of limited data [19]. Cazacu et al. [20] developed another yield function which was isotropic and insensitive to pressure. He also used the tactics of stress deviator transformation to capture tension-compression asymmetry and anisotropy simultaneously. Plunkett et al. [21] incorporated another linear transformation to improve the accuracy of anisotropy representation.
In addition to the SDE caused by directional twinning, the plastic yield of single crystal and polycrystalline of steel and aluminum alloys was significantly affected by superimposed hydrostatic pressure which was observed by Spitzig and his colleagues [22,23]. To predict the pressure-sensitive deformation behavior, Stoughton and his colleague [24] developed a yield criterion under nonassociated flow rule (non-AFR) to model the SDE. Lou et al. [25] extended the Yld2000-2d symmetric yield criterion by incorporating the first invariant and validated it to aluminum alloys and HCP metals. Yoon et al. [26] presented a yield criterion that was asymmetric based on three stress tensor invariants to model SDE of pressure-sensitive metals. Due to simplicity purpose, the pressure sensitivity parameters of the first stress invariants were simplified to be uniform.
In cases such as incremental sheet forming, hydroforming, and earing during cup drawing where hydrostatic pressure plays a nonnegligible role, it is vitally important to take into account the stress along the thickness direction when modeling the constitutive model of deformation. In the paper, the first stress invariant is incorporated into the proposed Drucker yield function, considering the influence of external hydraulic pressure. In order to reflect the difference of pressure directionality during loading, the first stress invariant is updated with three-directional coefficients. The yield function also relies on the third invariant of deviatoric stress tensor, which accounts for SDE. In order to meet the demand of postulated incompressibility during plastic deformation of voidfree materials, proposed three-dimensional yield criterion which is sensitive to hydrostatic pressure is used with non-AFR. Consequently, an extended Drucker yield function considering three-directional sensitivities to the first stress invariant and regarding SDE due to the second deviatoric stress invariant is proposed. The yield function under AFR and non-ARF is given using different transformation matrix. Constitutive equations along different directions are briefly introduced. The downhill simplex algorithm is utilized to solve high non-linear equations in order to obtain the anisotropic parameters. BBC metal AA 2008-T4 and FCC metal AA 2090-T3 are used to predict anisotropic yield and directional plastic flow for cold-rolled metals with weak SDE. The proposed Drucker yield criterion is also used to represent the asymmetric yield and plastic anisotropy of the Zirconium clock-rolled plate with strong SDE. Four-point bending test is simulated on Zirconium clock-rolled plate to show its reliability of reflecting the SD effect and predicting the shift of the neutral plane.

General Description of Drucker Yield Function
The stress state denoted by σ or ij σ for anisotropic material can be used to determine three stress invariants shown below: where Mises σ indicates the stress of von Mises.
Drucker [27] gave a yield criterion in Equation (5) to consider the influence of the second and the third deviatoric stress invariants on plastic yield locus: Three stress invariants were considered in the yield function proposed by Gao et al. [13] in Equation (7): It was reported in [26] that the Gao's yield function was not quite reasonable and there was no evidence from the experimental results.
In order to employ the dependence of yield on the first stress invariant based on experimental findings of Spitzig and his colleagues [22,23] and to capture the SDE, the yield function which gave a linear dependence of the first stress invariant and asymmetry of the third stress invariant is in Equation (8).
with b and c are material parameters to adjust the influence of hydrostatic pressure and deviatoric stress tensor, respectively, and need to be calibrated with experimental data along different orientations, a is a material parameter defined by the uni-axial stress-strain curve and in Equation (9): The convexity of the proposed yield criterion in Equation (8) and the influence of material parameter b which was sensitive to hypostatic pressure on yield surface of the presented yield criterion was briefly discussed in Reference [26]. It was concluded that if parameter c satisfied 3 3 4,3 3 4 c   ∈ −   , the convexity of the yield function was guaranteed.
To precisely capture the influence of anisotropic pressure, in the paper the sensitivity parameters, such as x b , y b and z b , which represent the effect of pressure along the rolling direction (RD), transverse direction (TD) and normal direction (ND), are introduced into the Drucker yield function in Equation (10): Similarly, the material parameter a can be determined based on uni-axial stress-strain curve in Equation (11):

Extension to Anisotropy Yield Function
Cazacu [12] extended the Drucker yield criterion to consider material anisotropy on the basis of representation theory and tress tensor transformations. Yoshida [14] developed the yield criterion as a summation of limited components of the Cazacu and Barlat [12] which were determined by the stress tensor linear transformation. Based on the on-going works, the yield function proposed will be extended by the linear transformation of the Cauchy stress tensors to describe metal anisotropy during plastic deformation in Equations (12)- (14): x xx y yy z zz where coefficients ′ C and ′′ C are two virtue matrix, T is a matrix used to transform Cauchy stress coordinate to deviatoric stress coordinate. The material parameters a and c in Equation (10)  Lou et al. [15] assumed the two linear transformed stress tensors to be the same, i.e., ′ ′′ = s s , leading to six reduced anisotropic parameters excluding parameter c which is regarded as 2 for FCC metals and 1.226 for BCC metals. In order to improve the simulation accuracy, they suggested a 7-parameter Drucker yield criterion to be described in the form of summation [14] which in turn inducing a pronounced increase of anisotropic parameters.
If non-AFR is used for modeling the anisotropy of metals, Equation (12) has the following form and the plastic potential is determined by r-values along different orientations in Equation (18): Similarly, two different fourth-rank linear transformation tensors are used to transform the Cauchy stress tensor σ to two stress tensors ′  s and ′′  s defined in Equations (21)- (23): L as 12 anisotropic parameters to be determined.
Among these twelve anisotropic parameters, eight parameters listed as ( )  are used to determine the anisotropy of in-plane plastic deformation. These parameters are calibrated by r-values along different directions and the uni-axial tensile yield stresses (also uni-axial compressive yield stresses if available). Because of the difficulty to measure the mechanical properties along the thickness direction of the thin sheet, the rest parameters are simplified as 4

Constitutive Equations along Different Directions
The paremeter Tθ σ is denoted as an un-axial tensile yield stress measured along a direction of θ in relation to the RD. Each stress component along various directions can be described below: Substituting Equation (24) into Equation (12), the tensile yield stress of Tθ σ will be provided in Equations (25)-(29): /3  3 2  2  2  2  2  2  2  11  21  12  12  22   1   cos  sin  3 Thus, the un-axial tensile yield stress along 0°, 45° and 90° directions are expressed in Equations (30)-(32), respectively: In the case of in-plane biaxial tensile, the stress components are reduced to Equation (33): and the biaxial tensile yield stress of Tb σ is expressed in Equation (34): When the specimen is compressed, the uni-axial compressive yield stress along RD is denoted as 0 C σ , and the stress components are expressed: Substituting Equation (35) into Equation (12), the uni-axial compressive yield stress along RD is expressed in Equation (36): Similarly, the uni-axial compressive yield stress along TD and DD are given in Equations (37)-(38), respectively: The in-plane equi-axial compressive yield stress named as Cb σ is obtained in Equation (39):

Calibration of Anisotropic Parameters
In case of Drucker yield function under associated flow rule (AFR), four tensile yield stresses ( Exp. Exp. Exp. 0 4 5 9 0 In case of non-AFR Drucker yield criterion, since the yield and the potential functions are individually expressed with stresses and r-values, the yield function parameters are determined merely depending on uni-axial tensile and compressive yield stresses at every 15 degree from the RD Exp.
Exp. 55 0 The equations listed are of high non-linearity and the anisotropic parameters in the described yield criterion are obtained by optimization strategies. The downhill simplex algorithm was widely used by authors [26,28] to apply to nonlinear optimization problems for which derivatives may not be easily obtained.

Application to AA 2090-T3 and AA2008-T4 with Weak SDE
The modified Drucker yield criterion stated is implemented to anisotropic aluminum alloys AA 2090-T3 and AA2008-T4 with weak SDE. Table 1 shows the experimental data described in detail by Yoon [28,29] and are supplied to calibrate the anisotropic material parameters in the extended Drucker yield function with the method described in Section 5. In order to clarify the simulation accuracy of the extended Drucker yield criterion, the simulation results with popular yield functions Hill's 48 and Yld2004-18p are adopted to supply anisotropic parameters. The calibration of anisotropic parameters of the Hill's 48 yield function is on the basis of yield stresses   In the paper, two types of Drucker yield functions are selected for comparative study. The one denoted as Ducker-σ-R is determined by Equation (40), while the one named as Ducker-σ is determined by Equation (41). The optimized anisotropic material parameters are listed in Table 6 for AA 2090-T3 and Table 7  and Tb σ obtained by experiment. Only the Ducker-σ can precisely describe the SDE. The pronounced difference between the yield surfaces denoted as Ducker-σ-R and Ducker-σ reveals attention must be paid when choosing objective functions to optimize material parameters.   Drucker-σ-R yield function shows a noticeable discrepancy when the angle begins from 60°, probably Hill48-σ Yld2004-18p Drucker-σ-R Drucker-σ Exp. data [28,29] due to the effect of r-values on the prediction of tensile yield stress. Due to the symmetric character of Hill's 48 and Yld2004-18p yield criteria, they are not capable of correctly predicting the uni-axial compressive yield stress. The modified Drucker yield function Drucker-σ precisely predicts the compressive yield stresses along all directions.   Hill48-σ Yld2004-18p Drucker-σ-R Drucker-σ Exp. data [28,29]  Hill48-σ Yld2004-18p Drucker-σ-R Drucker-σ Exp. data [28,29]   In case anisotropy prediction of tensile r-values, Drucker-R potential function determined by Equation (42) is used instead of Drucker-σ yield function, i.e., non-AFR will be implemented. The calibrated material parameters of Drucker-R potential function are listed in Table 8 for AA 2090-T3  and Table 9 for AA 2008-T4, respectively. Figures 7 and 8   Hill48-σ Yld2004-18p Drucker-σ-R Drucker-σ Exp. data [29,30]  Hill48-σ Yld2004-18p Drucker-σ-R Drucker-σ Exp. data [28,29]

Application to Zirconium Plate with Significant SDE
The developed Drucker yield function formulated in Section 3 is applied to predict the yield surface evolution of a clock-rolled Zirconium plate with significant SDE. A series of clock-rolling and annealing processing cycles were carried out, resulting in typical basal texture (where c -axes of Hill48-R Yld2004-18p Drucker-σ-R Drucker-R Exp. data [28,29]  Hill48-R Yld2004-18p Drucker-σ-R Drucker-R Exp. data [28,29] the crystals are predominately accordant to the plate normal direction). For sake of obtaining isotropic in-plane texture, multiple rolling passes with rotation were carried out. The visco-plastic self-consistent (VPSC) polycrystalline model [30] was used by Plunkett [31] on the basis of experimental observations (structures and textures, stress-strain relations of uniaxial tensile\compressive test and micromorphology) to model the evolution of yield surfaces. The calculation results of VPSC analytical model show that twinning induced by stretching dominates in tension at room temperature while pyramidal c a + -slip prevails along the c -axis in compression. Their experimental results have shown that the mechanical response of the plate significantly relies on the dominant orientation of the c -axis regarding the loading direction.
In Figure 9, the high-purity Zirconium indicates a comparatively severe basal texture showing almost in-plane axisymmetry. It indicates that its mechanical behavior is prone to be controlled by the direction and intensity of the applied load. To verify the predictability of the extended Drucker yield function, the strain distribution, as well as the movement of the neutral layer due to the significant difference between tensile and compressive behavior in a four-point beam bending test, were simulated. The finite element simulation was conducted based upon the work of Tomé [32] and Kaschner [33] as they have reported experimental results about the deformation anisotropy of textured polycrystalline pure Zirconium. Hill's 48, Yld2004-18p and the extended Drucker yield criteria are performed in the test simulation by ABAQUS/Explicit user subroutines. There are two upper pin dowels moving downward and two lower pin dowels holding stopped in the four-point bending fixture. The displacement of the upper pin dowels is 6 mm. The distance between the center of two upper pins is 12.7 mm and 38.1 mm of two lower pins. The schematic diagram of the four-point bending test and the geometry dimensions of the beam are shown in Figure 10. Due to geometric symmetry, only half of the bend geometry is established in the finite element simulation and symmetric boundary conditions are applied. The conditions of the four-point beam bending test, including element size and boundary conditions, were taken from [33,34]. Eight-node hexahedral elements (C3D8R in ABAQUS 6.14/CAE, Dassault Systemes Simulia Corp., Providence, RI, USA) with selective reduced integration are used to discretize the beam, shown in Figure 11. The number of elements along transverse and normal directions of the beam cross-section is six and along the longitudinal direction is 75. There are 2700 finite elements in total, which is the same as the one in Reference [34]. There are two kinds of beam orientations considered in all the yield criteria, supposing the rolling direction is in accordance with the x-axis: (i) case through-thickness bending (TTB), the normal direction of the beam is predominantly consistent with the global z-axis; and (ii) case in-plane bending (IPB), the beam normal direction is mostly coincident with the global y-axis. The initial crosssectional dimensions of the beam are 6.35 mm × 6.35 mm, with a longitudinal length of 50.8 mm. In the paper, the law of isotropic work hardening is considered, simplifying the hardening model as the defined yield surface expanding isotropically during plastic deformation. The maximum strain of the outmost layer of the bent beam is about 0.15 and as reported in Reference [35], the Zirconium plate shows the significant SDE when the plastic strain is of 15%. Thus, in the paper the anisotropic parameters corresponding to the strain of 15% are used to optimize the anisotropic parameters of the yield criteria.
According to Reference [35], the anisotropic parameters for the Zirconium plate for the strain of 15% is presented in Table 10. The data includes r-values, uniaxial tensile and compressive yield stresses along different directionalities. Besides, the equibiaxial tensile yield stress and r-value are also supplied. To describe the material plastic deformation behavior as close as possible, the above mentioned anisotropic data are used to calibrate Hill's 48, Yld2004-18p and two extended Drucker yield criteria using the functions and optimization method stated in Section 5. Tables 11-13 give the calibrated material coefficients of Yld2004-18p, Hill's 48 and two extended Drucker yield criteria for Zirconium.    Figure 12 shows the yield surfaces for Zirconium with different yield criteria. The predicted ones using extended Drucker yield function reproduce the deformation behavior quite well for all stress states. There are fewer differences between two kinds of Drucker yield functions, indicating there is no noticeable discrepancy when using AFR or non-AFR to predict the initial yield locus. Because of the symmetry of Hill's 48 and Yld2004-18p yield function, though they have the ability to well predict the stress ratio in a tensile stress state, they cannot predict the tension-compression asymmetry, leading to significant underestimation of equibiaxial compression value.      After calibration, the four-point bending test is simulated and the loading force, strain distribution and beam configuration are predicted and compared. Figure 16 presents the evolution of stroke force for the Zirconium beam with different yield functions, indicating that the punch forces of two extended Drucker evolve in almost the same way. The same tendency can be found with Hill's 48 and Yld2004-18p yield functions. In case of TTB, the load force for symmetric yield functions of Hill's 48 or Yld2004-18p is lower than the extended asymmetric Drucker yield functions, because of their underestimations of yielding in compression. In the case of IPB, the punch force evolutions obtained for Zirconium are almost the same no matter symmetric or asymmetric yield functions are implemented, indicating nearly isotropic mechanical property distribution in-plane direction. Hill48-R Yld2004-18p Drucker-σ-R Drucker-R Exp. data [35]  Hill48-σ Yld2004-18p Drucker-σ-R Drucker-σ Exp. data [35] (a) (b) Figure 16. Predicted punch loading evolution for various yield criteria for Zirconium: (a) case through-thickness bending (TTB); (b) case in-plane bending (IPB). Figure 17 shows the final configurations of the simulated results with various yield functions and superimposed to the experimentally obtained data (black dots) in Reference [33]. Regarding the predicted deformed beam central cross-section, all models predict similar wedge-shaped crosssections for the TTB case as the hard-to-deform c -axes are mostly in accordance with the z-axis.
The reason for the wedge-shaped cross-sections is due to a stress state of nearly uni-axial tension at the outmost layer along the beam longitudinal direction (the x-axis) and changing to almost uni-axial compressive stress state from the neutral plane to the innermost layer.
It is apparent that the predictions using two extended Drucker models are in coincidence with the observed data. This is due to the fact that the extended Drucker models can incorporate SDE quite well. Figure 17b is the local enlarged zone of case TTB (the dashed line rectangular at the low right corner) and it shows that the extended Drucker yield criterion under non-AFR can model the SDE better.
In case of IPB, where the c -axes are in accordance to the x-axis, the predicted contours of all yield functions well coincide with the experimental data, i.e., rectangular cross-section shown in Figure 17c, indicating correct prediction of the rigidity in the direction of hard-to-deform c -axes of all the mentioned yield functions.   Figure 18 shows the strain distributions along the direction of z-axis in the central cross-section and compared with observed experimental data obtained by a local strain measurement method determined by dot-matrix deposition and mapping [34]. Figure 18a is the case of TTB, where all the yield criteria can predict similar plastic strain distributions along the normal orientation of Zirconium plate where the hard-to-deform c -axes are predominately parallel to. But for the plastic strain distribution along x-axis, the extended Drucker yield criteria are in better coincidence with observed data, particularly with the increase of the absolute value of beam height. The unbalanced strain distributions lead to the shift of the neutral layer predicted by extended Drucker yield criteria considering SDE.
For IPB case, compared with experimental data, all the yield criteria can predict the distribution of strain components quite well. Because in IPB case, the hard-to-deform direction is predominately parallel to the y-axis, the strain distribution difference between the outmost and innermost fibers of the rectangular Zirconium beam is relatively small compared with that in the case of TTB. The maximum value of plastic strain according to Figure 18b, obtained by simulation or experiment, is about 17%, which is approaching 15% plastic strain value chosen to define the yield and hardening locus for the Zirconium plate. Beam height /mm

Conclusions
The Drucker yield criterion is extended to asymmetric yield criterion to model the anisotropic deformation behavior and SDE of metals sensitive to hydrostatic pressure. By the introduction of two fourth-order linear transformation tensors of Cauchy stress components, the extended Drucker yield function is transformed into anisotropic form to simulate the orthotropy and asymmetry of deformation behavior for cold-rolled steel sheets.
The proposed Drucker yield criterion is validated by the comparisons of the numerically predicted stress directionalities and r-values of planar anisotropy, in addition to the shape of yield surfaces, with the relevant experimental data on FCC and BCC metals. The proposed Drucker yield criterion is also used to describe the anisotropic and asymmetric yield of the clock-rolled Zirconium plate with strong SDE. The comparison shows that the extended Drucker yield criterion is feasible to model both SDE and anisotropy of cold-rolled HCP sheets with significant tension-compression asymmetry.
Four-point bending test was simulated on Zirconium clock-rolled plate, considering the plate normal direction aligned with different orientations. The comparison between simulation and experiment results indicate that the evolution of bending force is severely dependent on the accurate prediction of the tension-compression asymmetry. Regarding to the contour of the cross-section of the bent Zirconium beam, the predictions of the two extended Drucker yield functions are in better coincidence with experimental data. To clarify the differences among various yield functions and the influence of crystal orientations on deformation behavior, the strain distributions along the z-axis of the central cross-section are also compared with that of experimental data. The figures reveal that the simulation results using the extended model are closely related to the experimental data. It can capture the stiffness response of the bent beam along the hard-to-deform c -axis paralleling to the z-axis. As a result, the simulation model can reflect the effect of the SD and correctly predicts the shift of the neutral plane.
Author Contributions: B.T. designed the work; supervised; revised manuscript; Z.W. analyzed and prepared manuscript for the work; N.G. acquired, analyzed, and interpreted data for the work; Q.W. analyzed the data and improved language; P.L. contributed FE analysis tools. All authors have read and agreed to the published version of the manuscript. Beam height /mm