Characterization and modeling of biaxial plastic anisotropy in metallic sheets Zhenkai Mu a, Jiale Liu a, Wei Wang a, Xuerui Dai a, Shibo Ma a,*, Yong Hou b,* a Hebei Key Laboratory of Material Near-Net Forming Technology, School of Materials Science and Engineering, Hebei University of Science and Technology, Shijiazhuang 050018, PR China b Institute of Forming Technology and Lightweight Components (IUL), TU Dortmund University, Dortmund 44227, Germany A R T I C L E I N F O Keywords: Yield criterion Plastic anisotropy Associated flow rule Cruciform biaxial tension Principal stress direction Anisotropic hardening A B S T R A C T The anisotropic behavior of cold-rolled sheet metals has been extensively studied, typically characterized by uniaxial loading tests in different directions to determine yield strengths and plastic flow (Lankford r-values). However, biaxial principal stress states often focus solely on yield loci in the RD/TD (rolling/transverse) di- rections, with limited studies on other loading directions. This study analyzes the relationship between the principal stress yield loci and the material principal anisotropic stress directions based on the Hill48 yield function. Biaxial tensile tests are conducted on DP590 high-strength steel using cruciform specimens cut in various directions different from the sheet RD/TD directions. The evolution of the anisotropic yield locus and plastic flow of DP590 under biaxial tensile directions beyond the traditional RD/TD anisotropic directions is revealed. Furthermore, a modified Hill48 model that considers any principal stress direction of loading is pro- posed, accurately describing the continuous evolution of equi-biaxial tension stress, near-plane strain stress, and plastic flow direction for different principal stress directions of loading. The issues of failing to predict plastic flow under biaxial loadings in the original Hill48 model, and the degraded representation of existing constitutive models when accounting for out-of-plane anisotropy are addressed by decoupling the relationship between anisotropic parameters and stress state. In addition, the proposed calibration strategy of anisotropic parameters, using experimental data from various principal stress directions of loading, effectively captures the anisotropic evolution caused by changes in principal stress directions. Comprehensive evaluation and verification of the plasticity model should consider yield locus for any principal stress directions of loading, and also the normal and shear planes. 1. Introduction With the requirements of energy conservation and emission reduc- tion, the demand for weight reduction in the automotive industry has increased significantly. Advanced high-strength steels are widely used to cope with the increasingly fierce competitive environment [1,2]. How- ever, the sheet metals produced by rolling process have oriented texture, which results in obvious anisotropic behavior. Under the complex stress states and strain paths, there also has significant distortional hardening effect [3–6], tension-compression asymmetry [7–11], anisotropic creep deformation [12], et al., which significantly affect the formability of the sheet metals [13–15].The numerical simulation technology provides an important support for the development of sheet metal stamping process [16–19]. Due to the simulation results strongly depend on the constitutive models, the improve of the prediction accuracy under complex plastic deformation need the anisotropic constitutive model with higher accuracy. The anisotropic behavior received attention from the famous scholar Rodney Hill as early as 1948 [20], and the first anisotropic yield model was established based on the Mises isotropic yield function. After that, the anisotropic behavior of sheet metal was gradually taken into account in the stamping. A series of improved Hill yield criteria have appeared successively [21–24]. With the development of new materials, aluminum alloy with face-centered cubic (FCC) crystal structures and Mg alloy with hexagonal close-packed (HCP) crystal structures show anisotropic properties different from traditional steel materials. There- fore, Hosford [25] established a non-quadratic yield function that can adjust the curvature of the yield locus by introducing exponential * Corresponding author. E-mail addresses: mashibo1980@163.com (S. Ma), Yong.Hou@iul.tu-dortmund.de (Y. Hou). Contents lists available at ScienceDirect International Journal of Mechanical Sciences journal homepage: www.elsevier.com/locate/ijmecsci https://doi.org/10.1016/j.ijmecsci.2024.109640 Received 23 June 2024; Received in revised form 29 July 2024; Accepted 9 August 2024 International Journal of Mechanical Sciences 282 (2024) 109640 Available online 11 August 2024 0020-7403/© 2024 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY license ( http://creativecommons.org/licenses/by/4.0/ ). mailto:mashibo1980@163.com mailto:Yong.Hou@iul.tu-dortmund.de www.sciencedirect.com/science/journal/00207403 https://www.elsevier.com/locate/ijmecsci https://doi.org/10.1016/j.ijmecsci.2024.109640 https://doi.org/10.1016/j.ijmecsci.2024.109640 https://doi.org/10.1016/j.ijmecsci.2024.109640 http://crossmark.crossref.org/dialog/?doi=10.1016/j.ijmecsci.2024.109640&domain=pdf http://creativecommons.org/licenses/by/4.0/ parameter. Based on this, the Hosford series of yield criteria were developed. Among them, a landmark method was proposed by Barlat et al. [26] to construct an anisotropic yield model based on the linear transformation of the Cauchy stress tensor. Based on this method, the scholars constructed various advanced yield models with more aniso- tropic parameters, which aimed to make up for the deficiency that the Hill48 model can only predict the stresses or r-values in 0◦, 45◦, and 90◦ to the rolling direction [27–33]. Based on the Barlat89 yield models, Banabic et al. [34,35] proposed BBC families anisotropic yield models. After that, Comsa and Banabic [36] established the BBC2008 model for strongly anisotropic materials, it was expressed in the form of a finite series that can be expanded to retain more or fewer terms. These func- tions can modify the strength under shear and plane strain tension by changing the value of the exponent. Although the linear transformation method is very effective, the solution of its derivations is more compli- cated. Soare et al. [37–40] developed a new method to establish anisotropic yield functions based on homogeneous polynomials. The yield models can be polynomial functions of order 4, 6 and 8, respec- tively, which showed satisfactory results in describing the anisotropic behavior of BCC materials. In addition, the Druker series yield models based on stress invariants have also received extensive attention and research. Cazacu and Barlat [41] established a yield model based on Drucker yield function. In this model, the J2 and J3 invariants are replaced by J02 and J03 with aniso- tropic parameters. Cazacu and Barlat [42] proposed an odd yield func- tion with stress invariants to capture the SD effect. Then, they further proposed a new isotropic yield criterion with homogeneous function and extend it to anisotropy [43]. Yoshida et al. [44] proposed a 3D yield function of sixth-order polynomial type to describe the anisotropic behavior accurately. Lou et al. [45] coupled the pressure effect to the Yld2000–2d function to model the anisotropy and SD effect of sheet metals. Yoon et al. [46] modeled the anisotropy and SD effect by introduced the first invariant. Hu et al. [47] proposed a normalized stress invariant-based yield criterion for anisotropic-asymmetric mate- rials. Lou and Yoon [32] extended the Drucker yield function with non-associate flow rule (non-AFR) and the sum of n-components of the anisotropic Drucker function. Cazacu [48] extended the Drucker func- tion to eight order to improve the description of FCC metals. Lou et al. [49] coupled the Drucker yield function with the pressure-sensitivity to describe strength from shear to plane strain tension under proportional loadings. The advantage of these yield functions is that the yield func- tions can be utilized under both plane stress condition and the spatial loading cases conveniently, and the derivatives can be easily computed analytically. In the early construction of the constitutive model, the associated flow criterion was usually used. To describe the anisotropic properties more comprehensively, the anisotropic parameters of the yield models are continuously increased, which makes the order and form of the yield functions extremely complicated. Under this condition, it is often necessary to solve the anisotropic parameters by means of optimization algorithms, which is unfriendly to the prediction of subsequent yield evolution. The non-associated flow criterion improves the construction strategy of the constitutive model [50–52]. Not only the form of constitutive function is simplified, but also the prediction accuracy is improved. In addition, the development of analytical yield models has been paid more attention. Hu and Yoon [53] proposed the analytical method for asymmetric yield function Yoon2014 by considering aniso- tropic hardening under non-associated flow rule. In recent years, the method of coupling different yield functions to establish anisotropic models is widely used, and many models with excellent prediction accuracy have been derived. Lee et al. [54] pro- posed the coupling of quadratic and non-quadratic yield functions (CQN) to describe anisotropic hardening and control the shape of the yield surface. Park et al. [55] proposed a model with the multiplication of two additional terms referred to as scaling and asymmetry functions, which had noticeable flexibility and general applicability in describing the anisotropic/asymmetric yielding behavior. Hu et al. [56] coupled Poly4 with Hosford models to improve the prediction accuracy of anisotropy. The framework by Lee et al. [54] was modified by Hou et al. [57,58], the new yield criterion increased the flexibility and range in controlling the shape of the yield surface of the entire model. Subse- quently, the second and third stress invariants of the deviator stress tensor were also used to develop the coupled yield function [59,60]. Based on the multiplicative coupling method, additive coupling models have been proposed and show more excellent performance. Hou et al. [61] proposed a new anisotropic-asymmetric criterion by adding the asymmetric Poly4 function and the isotropic stress invariant-based yield function. Then, they further proposed a generalized versatile plasticity model framework to describe the tension and compression asymmetry [62]. In addition, Lou and Yoon [63] introduced other novel approach, employing a Lode-dependent function to extend anisotropic yield functions into asymmetric ones, focusing on the Hill48 function and Yld91 function [64]. Zhou et al. [65] extended the symmetric yield criterion SY2009 to an asymmetric one using shape functions, accurately predicted tensile and compressive r-values along different loading di- rections. Mu et al. [66] improved the prediction accuracy and flexibility of Hill48 model by introducing stress correction coefficient and linear interpolation method. Zhang and Lou [67] proposed a new yield func- tion for anisotropic-differential hardening by introducing a stress state-sensitive (SSS) function. Then, they further proposed a five stress states-sensitive (Five-SSS) yield function by considering the non-linear effect of stress triaxiality η on plasticity, which was used to model the springback of V-bending [68]. The accurate description of anisotropic evolution will directly determine the forming accuracy of sheet metal, especially the yield surface evolution under biaxial loadings. Therefore, the prediction of anisotropic evolution is an important aspect in evaluating the ability of constitutive models, including the directional hardening curves, r- values, and the subsequent yield loci under proportional and non- proportional loading path [69–72]. In addition, the distortional hard- ening of yield surface under non-linear loading conditions, variation of strain rate and temperature have also been paid attention and studied [73]. For this, numerous phenomenon-based constitutive models were proposed to describe hardening behaviors under non-proportional loadings by shifting the yield surfaces, i.e., kinematic hardening, by distorting the yield surfaces or by the combination of them. A homo- geneous anisotropic hardening (HAH) model was developed by Barlat et al. [74] to replace the kinematic hardening model. On this basis, Lee et al. [75] combined the HAH model with the CQN model to accurately predict the yield evolution behavior and Bauschinger effect. Shi et al. [76] applied thermodynamically consistent directional distortional hardening model for AZ31. Zhang et al. [77] introduced a generalized François’s distorted stress term to any non-quadratic yield functions. Raj et al. [78] calculated the evolution of plastic work contour for different stress path. He et al. [79] developed a phase transformation-informed non-associated flow model incorporated with an asymmetric yield cri- terion, distortional anisotropic hardening, and phase transformation-related hardening law to capture complex anisotropic evolutions for QP steel. Lee and Rubin [80] developed a method to capture strength-differential effect (SDE) based on the distortion of yield locus, which used a SDE parameter to replace the first stress invariant. For modeling yield surface evolution, the constitutive parameters are usually interpolated by linear or non-linear equations [81–83]. In order to overcome the complexity and time-consuming of the optimization method in engineering applications, the analytical parameter compu- tations of yield functions were developed [54,60,84]. Then, Hu et al. [85] further gave an analytical solution strategy for Poly6 yield func- tion. Du et al. [86] proposed a new parameter calibration strategy for analytical Poly6 yield criterion, and gave the systematic evaluation strategy for material models. Recently, the recognition of identification of anisotropic mechanical properties and establishment of constitutive Z. Mu et al. International Journal of Mechanical Sciences 282 (2024) 109640 2 model combined with artificial intelligence algorithm have been paid more attention. Jeong et al. [87] calibrated the parameters of Poly6 yield criterion using indentation plastometry based on a neural network (NN) system. Liu et al. [88] calibrated the phenomenological constitu- tive model and hardening model based on the established virtual labo- ratory. Wessel et al. [89] introduced a new machine learning-based sampling approach to perform virtual experiments with respect to the full stress state to identify the parameters of anisotropy yield functions. However, most of the current studies focus on the principal plane yield loci [60,63,86,90–96]. The materials experience complex loading histories in stamping, and the stress state is constantly changing. The yield models with anisotropic parameters calibrated only by the prin- cipal plane yield loci are insufficient to describe the yield behavior under other stress states. In other words, the principal stress component obtained from the stress components on the yield surface through the coordinate transformation will not completely coincide with the yield locus in the RD/TD direction. In order to accurately describe the whole yield surface, it is necessary to calibrate the anisotropic parameters ac- cording to the yield stress data under different principal stress directions of loading. However, there are few studies on the difference of principal plane yield loci in different principal stress directions of loading and their evolution with plastic deformation. In addition, existing constitu- tive models have fixed anisotropic parameters, which can only be cali- brated with specific stress or deformation data. It cannot use the stress or deformation data in other biaxial loading directions for parameters calibration, which makes it impossible to predict the diversity of yield loci related to the principal stress directions of loading, but can only describe a few specific yield loci. In this work, the evolution of anisotropic behavior in sheet metals under various biaxial tensile directions is examined, and a novel method for describing these variations is proposed. In Section 2, the relationship between principal stress yield locus and principal stress direction of loading is analyzed using the Hill48 yield function. In Section 3, biaxial tensile tests are conducted on DP590 high-strength steel using cruciform specimens cut in various directions. The evolution of principal plane yield locus and plastic flow directions with changes in the PSD is investigated. In Section 4, a modified Hill48 model that considers the principal stress direction of loading is established. The anisotropic pa- rameters of this model can be calibrated using stress or deformation data under any biaxial loading directions. The reliability and prediction ac- curacy of the modeling strategy are verified and discussed by comparing it with recently proposed advanced yield models. 2. Effect of principal stress direction of loading on yield loci For the yield function established with the stress components in the reference coordinate system, the relationship between the principal stress and the anisotropic principal axis should be discussed in solving the yield locus. When the principal stress directions of loading coincides with the sheet material’s principal anisotropic directions (RD/TD axis), the (x,y)-axis stress components in the reference coordinate system are equivalent to the principal stress components. In this case, the principal plane yield locus can be obtained by the yield locus function. While the principal stress directions of loading do not coincide with the material principal anisotropic direction, it is necessary to transform the stress tensor to the principal axis coordinate system to obtain the principal stress yield locus function. Obviously, the principal stress yield locus function must contain the parameter of the principal stress direction of loading. Therefore, in this section, the diversity of principal stress yield locus related to anisotropic parameters and loading directions will be clarified. 2.1. Yield function expressed with principal stresses Under the material anisotropic principal (x,y)-axis coordinate sys- tem, the Hill48 yield function for plane stress states can be expressed as follows 2fy ( σij ) = Fσ2 y + Gσ2 x + H ( σx − σy )2 + 2Nτ2xy = 1. (1) When the angle between principal stress direction of loading (σ1) and the rolling direction (x-direction) is α, according to the stress Mohr circle, the relationship between the principal stresses and the aniso- tropic stress components in the (x,y)-axis coordinate system can be expressed as ⎧ ⎪⎪⎪⎪⎪⎪⎨ ⎪⎪⎪⎪⎪⎪⎩ σx = σ1 + σ2 2 + σ1 − σ2 2 cos2α σy = σ1 + σ2 2 − σ1 − σ2 2 cos2α τxy = σ1 − σ2 2 sin2α . (2) By substituting Eq. (2) into the yield function, the yield locus func- tion expressed by the principal stress can be obtained as F [(σ1 + σ2 2 )2 − 1 2 ( σ2 1 − σ2 2 ) cos2α + (σ1 − σ2 2 )2 cos22α ] + G [(σ1 + σ2 2 )2 + 1 2 ( σ2 1 − σ2 2 ) cos2α + (σ1 − σ2 2 )2 cos22α ] + H [ (σ1 − σ2) 2cos22α ] + 2N [ (σ1 − σ2) 2sin22α / 4 ] = 1. (3) The Eq. (3) can be further simplified as k1σ2 1 + k2σ2 2 + k3σ1σ2 = 1, (4) where ⎧ ⎨ ⎩ k1 = Fsin4α + Gcos4α + 2(N − 2H)sin2αcos2α + H k2 = Gsin4α + Fcos4α + 2(N − 2H)sin2αcos2α + H k3 = 2 [ (F + G+ 4H − 2N)sin2αcos2α − H ] . (5) The Eq. (4) is a typical quadratic function. When k1 > 0, k2 > 0, the shape of yield locus is ellipse. If the angle between ellipse major axis and σ1 is β, the yield locus in the principal stress space can be shown as Fig. 1. 2.2. Diversity of yield locus in the principal plane If the coordinates are transformed to the directions of the major and minor axes of the ellipse, the locus becomes a standard elliptic function. Therefore, according to Eq. (4) and the properties of the elliptic func- tion, there is Fig. 1. Schematic diagram of yield locus in principal stress coordinate system. The red circle is the intersection point of the yield locus and the coordinate axis, representing the uniaxial loading stress states. Z. Mu et al. International Journal of Mechanical Sciences 282 (2024) 109640 3 tan2β = k3 k1 − k2 = (F + G+ 4H − 2N)sin22α − 4H 2(G − F)cos2α . (6) For isotropic sheet metal, N = 3F = 3G = 3H. Set the yield stress in the RD as σs, according to Eq. (1), thus, σs = 1/ ̅̅̅̅̅̅̅̅̅̅̅̅ F + H √ = 1 / ̅̅̅̅̅̅ 2F √ . Then, k1 = k2 = 1/σ2 s , k3 = − 1/σ2 s . There is ( σ1 σs )2 − ( σ1 σs )( σ2 σs ) + ( σ2 σs )2 = 1. (7) Under this condition, the angle between the major axis of ellipse and the principal stress axis is π/4, as shown in Fig. 2(a). The yield locus is only related to the yield stress, and the subsequent yield locus is an expanding ellipse. For sheet metal with in-plane isotropy and thickness anisotropy, the necessary and sufficient conditions to meet axisymmetric anisotropy isN = F + 2H = G + 2H. Then, it can be calculated as k1 = k2 = F + H, k3 = − 2H. By substituting it into Eq. (4), there is σ2 1 − 2H F + H σ1σ2 + σ2 2 = 1 F + H . (8) Setting the yield stress as σs and the anisotropy coefficient as r, it can be obtained as σs = 1/ ̅̅̅̅̅̅̅̅̅̅̅̅ F + H √ , r = H/F. According to the Eq. (8), there is ( σ1 σs )2 − 2r 1+ r ( σ1 σs )( σ2 σs ) + ( σ2 σs )2 = 1. (9) Under this condition, the yield locus is also elliptic. The angle be- tween the major axis of ellipse and the σ1 axis is also π/4, as shown in Fig. 2(b). Different from the isotropic plates, the yield stresses under symmetric biaxial loading (equi-biaxial tension, pure shear, equi-biaxial Fig. 2. Schematic diagram of yield locus in principal stress coordinate system for (a) isotropic plate; (b) sheet metal with in-plane isotropy and thickness anisotropy. The red circles represent the uniaxial loading stress states, and the green circles represent the symmetric biaxial loading stress states. Fig. 3. Schematic diagram of yield locus in different principal stress directions of loading (The angle between the principal loading direction and the anisotropic RD/ TD directions): (a) α = 0◦; (b) α = 45◦; (c) α = 90◦. Z. Mu et al. International Journal of Mechanical Sciences 282 (2024) 109640 4 compression) are no longer equal to the uniaxial tensile yield stress, but is directly determined by the thickness anisotropic coefficient r-value. When r > 1, the major axis of ellipse is lengthened and the minor axis is shortened. When r < 1, the major axis of ellipse is shortened and the minor axis is lengthened. The subsequent yield loci are still an expanding ellipse. For orthogonal anisotropic sheet, according to the Eqs. (4) and (5), the yield locus depends not only on the anisotropic parameters F, G, H and N, but also on the principal stress direction of loading. The yield locus has different forms for different principal stress directions of loading. Thus, several typical cases are discussed as below. (1) α = 0◦ (principal loading directions coincident with the aniso- tropic RD/TD directions), thus, k1=G+H; k2= F+H; tan2β = − 2H/ (G − F). Under this condition, the yield locus function can be expressed as (G+H)σ2 1 − 2Hσ1σ2 + (F+H)σ2 2 = 1. (10) If G > F, thus π/4 < β < π/2. If G < F, thus 0 < β < π/4. With the increase of plastic deformation, F, G and H monotonically decrease, the subsequent yield locus continues to expand, as shown in Fig. 3(a). (2) α = 45◦, thus, k1= k2= (F+ G+ 2N)/4; k3= (F+ G − 2N)/2; β = π/4. Under this condition, the yield locus function can be expressed as (F + G+ 2N) 4 ( σ2 1 + σ2 2 ) + 1 2 (F+G − 2N)σ1σ2 = 1. (11) The yield locus geometry is similar to that of isotropic plates and in- plane isotropic thick anisotropic plates, and the subsequent yield locus is an expanding isocline ellipse, as shown in Fig. 3(b). (3) α = 90◦, thus, k1 = F + H; k2 = G + H; tan2β = − 2H/(G − F). Under this condition, the yield locus function can be expressed as (F+H)σ2 1 − 2Hσ1σ2 + (G+H)σ2 2 = 1. (12) If G > F, thus 0 < β < π/4; if G < F, thus π/4 < β < π/2, as shown in Fig. 3(c). 3. Mechanical experiments for different principal stress directions of loading In this section, a DP590 steel plate with a thickness of 1.0 mm was investigated to verify the diversity of yield loci. Experiments were conducted for uniaxial tensile (UT) and biaxial tensile (BT) loading under various principal stress directions (PSD). Strain distributions were measured using digital image correlation (DIC). The yield loci and plastic flow were obtained and compared along the RD/TD, 30◦/120◦, and 45◦/135◦ directions of biaxial loadings. The results were analyzed to assess the variation in anisotropic behavior. 3.1. Uniaxial tension Uniaxial tensile tests along 0◦, 15◦, 30◦, 45◦, 60◦, 75◦ and 90◦ to the rolling direction (RD) were carried out by the Inspekt 100 kN electronic universal material testing machine, in accordance with ISO 6892–1 standards at room temperature. The tensile rate was 5 mm/min corre- sponding to the quasi-static loading. The curves of uniaxial tensile true stress with equivalent plastic strain (EPS) in different directions and the curves of longitudinal-transverse plastic strain were obtained by the plastic work equivalence principle, as shown in Fig. 4. Then, the experimental flow curves is smoothed and fitted using the hardening equation shown in Eq. (13) to characterize the hardening behavior accurately. σα = K+M1 ( 1 − e− εP/b1 ) +M2 ( 1 − e− εP/b2 ) , (13) where K, M1,M2, b1 and b2 are the fitted coefficients; εP is the equivalent plastic strain. The fitted result in seven loading directions for DP590 are shown in Table 1. The Lankford coefficient (r-value) characterizing plastic deformation is expressed as Fig. 4. (a) Uniaxial tensile true stress to equivalent plastic strain curves and (b) the transverse to longitudinal plastic strain curves at 0◦, 15◦, 30◦, 45◦, 60◦, 75◦ and 90◦ with respect to RD. Table 1 Fitted coefficients in the hardening equation for DP590 obtained from uniaxial tensile tests under seven loading directions with respect to RD. α (◦) DP590 K (MPa) M1 (MPa) b1 M2 (MPa) b2 0 324.03 345.40 0.2068 222.78 0.0236 15 325.59 348.69 0.2036 216.45 0.0234 30 332.22 368.11 0.2287 219.64 0.0245 45 337.54 225.95 0.0259 399.91 0.2808 60 341.01 241.91 0.0280 590.97 0.5269 75 340.98 473.29 0.3880 238.85 0.0275 90 343.92 238.08 0.0276 438.56 0.3520 Z. Mu et al. International Journal of Mechanical Sciences 282 (2024) 109640 5 r = εw εz = − εw εw + εl . (14) Based on the plastic strain in Fig. 4(b) and the plastic work equiva- lence principle, the discrete r-value data related to the equivalent plastic strain can be obtained. Furthermore, Eq. (15) is used to describe the continuous change of r-value with EPS. rα = Aα + Bα ∗ exp ( − εP /Cα ) , (15) where Aα, Bα and Cα are the fitted coefficients. The coefficients for fitting r-value of DP590 in seven loading directions are shown in Table 2. 3.2. Biaxial tension In order to verify the above analysis about the evolution behavior of the yield locus at different principal stress directions of loading, cruci- form specimens were cut at 0◦, 30◦ and 45◦ from the rolling direction, respectively, as shown in Fig. 5. Following ISO 16842, The biaxial ten- sile tests with seven load ratios were carried out on the cruciform specimens in the three directions respectively. To ensure continuous proportional loading during testing, the biaxial tensile tests were per- formed under force control mode, the average strain rate was 3.0e- 6~1.6e-5 s− 1 corresponding to the quasi-static loading. The load ratios of two loading directions were F1: F2= 4: 0, 4: 1, 4: 2, 4: 3, 4: 4, 3: 4, 2: 4, 1: 4, 0: 4. The biaxial tensile tests with load ratios of 4:0 and 0:4 were replaced by the uniaxial tensile tests along the di- rections of the specimen. When the principal stress direction of loading is α = 0◦, F1: F2 = 4: 0and F1: F2 = 0: 4are uniaxial tension along the principal anisotropic directions RD and TD axis respectively. When the principal stress direction of loading is α = 30◦, F1: F2= 4: 0and F1: F2= 0: 4are uniaxial tension in the 30◦ and 120◦ to the principal anisotropic directions RD axis. When the principal stress direction of loading is α = 45◦, F1: F2= 4: 0and F1: F2= 0: 4are uniaxial tension in the 45◦ and 135◦ to the principal anisotropic directions RD axis. Fig. A1 (in the Appendix A) shows the biaxial tensile true stress-strain curves in different prin- cipal stress directions of loading, which were used to calculate the yield locus. 3.2.1. Yield locus Based on the plastic work equivalence principle, the stress-strain relationship between the biaxial loading stress state and the uniaxial loading stress state is satisfied as w = ∫ε P 1 0 σ1dεP 1 + ∫ε P 2 0 σ2dεP 2 = ∫ε P α 0 σαdεP α. (16) According to Eq. (16), the yield loci of three principal stress di- rections of loading at different EPS levels are obtained as shown in Fig. 6. With the increase of plastic deformation, the plastic work con- tours of these three directions expands outward. Due to the distortional hardening, the shape of the yield locus is also constantly changing, that is, the yield locus in the subsequent yield deviates from the initial yield locus. In addition, due to the planar anisotropy in the rolled sheet, the upper and lower contours of the yield loci bounded by EBT are asym- metrical. In Fig. 6(d), the plastic work contours of the three principal stress directions of loading under different EPS are difference, which is consistent with the previous analysis. These differences are not only reflected in the initial yield, but also change with the increase of plastic strain, especially in the region of near-plane strain stresses state. Therefore, in addition to the evolution of yield loci caused by hardening, the variation of yield loci caused by the change of principal stress orientation should also be considered in the subsequent yield. In this way, the prediction accuracy of yield model can be further improved, and more accurate finite element simulation results can be obtained. To analyze the difference of yield loci in different stress principal axes more clearly, the relative relationship between the biaxial tension Table 2 Coefficients describing evolution of r-values for DP590 obtained from uniaxial tensile tests under seven loading directions with respect to RD. α (◦) Aα Bα Cα ra 0 0.7629 − 0.1561 0.0229 0.75 15 0.7724 − 0.2939 0.0043 0.75 30 0.843 − 0.3486 0.0083 0.82 45 0.8136 − 0.1342 0.0274 0.81 60 0.9762 − 0.1714 0.0469 0.96 75 0.9702 − 0.1291 0.0236 0.96 90 0.9576 − 0.0897 0.0181 0.95 a The average r-values calculated by linear regression in plastic strain range of 2 %− 12 %. Fig. 5. (a) Biaxial tensile testing machine, (b) dimensions of cruciform specimen and (c) cruciform specimens cut along 0◦, 30◦ and 45◦ principal stress directions of loading. The biaxial tensile tests were performed under force control mode. Strain distributions were measured using DIC. Z. Mu et al. International Journal of Mechanical Sciences 282 (2024) 109640 6 Fig. 6. Yield loci of DP590 along different principal stress directions of loading: (a) RD/TD; (b) 30◦/120◦; (c) 45◦/135◦; (d) Comparison of the yield loci obtained from three principal stress directions at different EPS. Fig. 7. Comparison between the biaxial tension stresses of yield locus obtained along (a) 30◦/120◦ and (b) 45◦/135◦ sampling direction and that in the RD/TD direction at different EPS. Z. Mu et al. International Journal of Mechanical Sciences 282 (2024) 109640 7 stresses of yield locus obtained along 30◦/120◦ and 45◦/135◦ sampling direction and that in the RD/TD direction are calculated by Eq. (17): δα PWC− βi = ̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅( σi 1− α − σi 1− 0∘ )2 + ( σi 2− α − σi 2− 0∘ )2 √ ̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅( σi 1− 0∘ )2 + ( σi 2− 0∘ )2 √ , (17) where “σi 1− α” and “σi 1− 0∘ ” indicate the first principal stress of biaxial tension along the α and 0◦ respectively; “i” indicates the different stress state. βi = arctan(γi), and “γi” indicates the ith principal stress ratio. Note that the “σi 1− α” and “σi 1− 0∘ ” have the same stress ratio (γi), that is γi = σi 2− α σi 1− α = σi 2− 0∘ σi 1− 0∘ . (18) Then δα PWC− βi is actually the ratio of the “distance” between the stress data points (with the same γi) to the “distance” between the stress data point and the origin in the yield locus diagram. By substituting Eq. (18) into Eq. (17), δα PWC− βi at a specific EPS level εP can be simplified as δα PWC− βi ( εP) = σi 1− P(ε P) σi 1− E(εP) − 1. (19) The comparison between the plastic work contours obtained by BT along the 30◦/120◦ direction and that in the RD/TD direction is shown in Fig. 7. The relative distance of the points on the plastic work contour of different principal stress directions decreases first and follow increases, and then decreases again change with φ. For the plastic work contours under 45◦/135◦ BT, due to the rotation of its principal axis is further increased, the variation law of the plastic work contours relative to RD/TD are more obvious. 3.2.2. Plastic flow In order to further analyze the plastic deformation behavior under different principal stress loading directions, biaxial tensile plastic strains under different loading stress states were obtained based on DIC, as shown in Fig. 8. Among them, the strain relationship curves of σ1: σ2= 0: 4 and σ1: σ2 = 4: 0 were derived from uniaxial tension under different principal stress directions of loading. According to the above plastic strain relationship data, the plastic flow directions under different biaxial tension stress states were calcu- lated, as shown in Fig. A2 (in the Appendix A). The results show that the change of β is obvious in the initial yield stage. With the increase of plastic deformation, the flow directions under different stress states tends to be saturated, that is, evolution does not develop. In addition, there are also obvious differences in the plastic flow directions under different principal stress directions of loading, as shown in Fig. 9. This difference can be easily inferred from the uniaxial tension r-values. The plastic flow directions corresponding to φ = 0◦ and φ = 90◦ were calculated according to the r-values under UT stress state. For aniso- tropic metal sheets, uniaxial tensile r-values change constantly with the angle to the RD. Therefore, for biaxial loadings along different PSD, the Fig. 8. The curves of true plastic strain with seven load ratios under different principal stress directions of biaxial loadings: (a) RD/TD; (b) 30◦/120◦; (c) 45◦/135◦. Z. Mu et al. International Journal of Mechanical Sciences 282 (2024) 109640 8 corresponding plastic flow directions of φ = 0◦ and φ = 90◦ must be different. Then, for the NPS and EBT stress states, it is also inevitable that the plastic flow directions under different PSD show significant difference. The results show that there is a strong correlation between the yield stresses (plastic flow directions) and the PSD under NPS and EBT. Therefore, to accurately describe the yield and deformation behavior, it is necessary to fully consider the stress state related to the PSD. Only a few scholars have paid attention to yield loci in the 45◦/135◦ direction, which is still insufficient according to the above experimental results. In addition, the anisotropic parameters of existing yield functions are usually calibrated only by the experimental data under biaxial loadings along the RD/TD direction, it cannot be possible to consider the test data under other principal stress directions of loading. Therefore, the existing models will not be able to describe the differences in yield surface caused by the variation of the principal stress directions of loading. 4. A new yield model coupled with the principal stress loading direction This section begins by outlining the process of modifying the Hill48 model to incorporate the principal stress loading direction in the anisotropic principal (x,y)-axis coordinate system. The calibration strategy for its anisotropic parameters under different stress states is then detailed. Following this, the modified Hill48 model is generalized to the principal stresses of loading coordinate system. The feasibility of the new model and the parameter calculation method for predicting the yield surface under a full stress state is discussed by comparing the yield loci and plastic flow directions predicted under different principal stress directions of loading with those predicted by advanced anisotropic yield models. 4.1. Modified Hill48 model in anisotropic principal axis coordinate system In plane stress state, the modified Hill48 yield function considering principal stress direction of loading and subsequent yielding can be expressed as [ G ( εP, α ) +H ( εP,α )] σ2 x − 2H ( εP, α ) σxσy + [ F ( εP,α ) +H ( εP, α )] σ2 y + 2N ( εP,α ) τ2xy = 1, (20) where F, G, H and N are the anisotropic parameters, which are related to the equivalent plastic strain and the principal stress direction of loading. Fig. 9. Variation of the direction of plastic strain rate with loading directions under several EPS levels and principal stress directions: (a) RD/TD; (b) 30◦/120◦; (c) 45◦/135◦. (d) Comparison of the directions of plastic strain rate obtained from three principal stress directions of loading at εP=0.005. Z. Mu et al. International Journal of Mechanical Sciences 282 (2024) 109640 9 ⎧ ⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎨ ⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎩ H ( εP,α ) = χ2h ( εP, α ) G ( εP,α ) = χ2 [σ0(εP)] 2 − H F ( εP,α ) = χ2 [σ90(εP)] 2 − H N ( εP,α ) = 2χ2 [σ45(εP)] 2 − F + G 2 , (21) where α is the angle between the principal stress axis and the rolling direction, it can be easily obtained from the stress Mohr’s circle ac- cording to the stress components; σ0(εp), σ45(εp) and σ90(εp) are the experimental yield stress with different equivalent plastic strains ob- tained by uniaxial tension along 0◦, 45◦ and 90◦ to RD, respectively. h(εp, α) is the intermediate parameter. χ is the stress correction coeffi- cient. They have different expression forms for different stress states. 4.1.1. Uniaxial loading stresses state According to the modified Hill48 yield function Eq. (20), the aniso- tropic parameters under the uniaxial stress state can be expressed as ⎧ ⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎨ ⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎩ Hu ( εP,α ) = χ2 u ( εP,α ) hu ( εP,α ) Gu ( εP, α ) = χ2 u ( εP,α ) 1 [σ0(εP)] 2 − Hu Fu ( εP, α ) = χ2 u ( εP, α ) 1 [σ90(εP)] 2 − Hu Nu ( εP,α ) = χ2 u ( εP,α ) 2 [σ45(εP)] 2 − Fu + Gu 2 , (22) where χu is the correction factor of uniaxial tensile stress, which can be solved by the experimental uniaxial tensile stresses as where σexp α (εp) is the experimental UT stress related to the EPS along α direction of loading. hu(εp, α) can be calculated from the uniaxial tensile r-values rexpα (εp) as 4.1.2. Biaxial loading stress state In the biaxial stress states, in addition to accurately predicting the yield loci, the plastic flow directions should also be described with sufficient accuracy, which will directly determine the flow behavior of the material. However, it is difficult to predict both yield loci and plastic flow based on the AFR. There are only a few higher-order yield func- tions, such as Yld2000–2d, BBC05 and Ploy6 et al., can better describe the plastic flow direction under EBT, due to their anisotropic parameters are calibrated by the ratio of principal strain of EBT. The establishment of non-AFR provides a good method to solve this problem, but its pa- rameters calibration is often complicated. At present, only the plastic flow of EBT and NPS under BT along RD/TD direction (normal plane) can been well captured. It is difficult to predict the plastic flow in other principal stress directions of loading. In order to describe the yield loci and the plastic flow directions simultaneously, the anisotropic parameters under the biaxial loadings including the biaxial stress correction factor and the plastic flow direc- tion can be expressed as ⎧ ⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎨ ⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎩ Hb ( εP,α ) = χ2 b ( εP, α ) hb ( εP, α ) Gb ( εP, α ) = χ2 b ( εP,α ) 1 [σ0(εP)] 2 − Hb Fb ( εP, α ) = χ2 b ( εP,α ) 1 [σ90(εP)] 2 − Hb Nb ( εP,α ) = χ2 b ( εP, α ) 2 [σ45(εP)] 2 − Fb + Gb 2 , (25) where χb is stress correction coefficients of the biaxial loadings, which can be calculated by χb = aγ2 + bγ + c, (26) where a, b and c can be determined by UT stress correction coefficient χu, EBT stress correction coefficient χeb, and NPS stress correction coeffi- cient χps. ⎧ ⎪⎪⎪⎪⎪⎪⎪⎪⎪⎨ ⎪⎪⎪⎪⎪⎪⎪⎪⎪⎩ a = χps − γpsχeb + ( γps − 1 ) χu γ2ps − γps b = − χps + γ2psχeb + ( 1 − γ2ps ) χu γ2ps − γps c = χu , (27) where hb is the intermediate variable for solving the plastic strain ratio χu ( εP, α ) = 1 σexp α (εP) ̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅ 1 [σ90(εP)] 2 sin2α + 1 [σ0(εP)] 2 cos2α + ( 4 [σ45(εP)] 2 − 2 [σ0(εP)] 2 − 2 [σ90(εP)] 2 ) sin2αcos2α √ √ √ √ , (23) hu ( εP, α ) = 1 1+ rexpα (εP) {( sin2α [σ90(εP)] 2 + cos2α [σ0(εP)] 2 ) rexpα ( εP) − ( 4 [σ45(εP)] 2 − 2 [σ0(εP)] 2 − 2 [σ90(εP)] 2 ) sin2αcos2α } . (24) Z. Mu et al. International Journal of Mechanical Sciences 282 (2024) 109640 10 of biaxial loading stress states (the derivation process is shown in Ap- pendix B) hb ( εP,α ) = ξ12(εP,α) / σ2 0 − γ / σ2 90 γξ12(εP,α) − 1 , (28) where ξ12 is the ratio of principal strain increment under biaxial load- ings. It can be calculated by the ratio of principal strain increment under UT (ξu), NPS (ξps), and the EBT (ξeb) as ξ12 = Aγ2 + Bγ + C, (29) with ⎧ ⎪⎪⎪⎪⎪⎪⎪⎪⎪⎨ ⎪⎪⎪⎪⎪⎪⎪⎪⎪⎩ A = ξps − γpsξeb + ( γps − 1 ) ξu γ2ps − γps B = − ξps + γ2psξeb + ( 1 − γ2ps ) ξu γ2ps − γps C = ξu , (30) where γ = |σ2|/|σ1|, σ1 and σ2 are the two principal stresses in plane stress state. σ2 is the smaller one of the two principal stresses, σ1 is the bigger one. γ = 0 represents the uniaxial loading stress state, γ = 1 represents the symmetric biaxial loading stress states. 0 < γ < 1 repre- sents the asymmetric biaxial stress state between the uniaxial and symmetric biaxial stress. γps is the NPS stress state. When the biaxial tension is along the RD/TD direction, the axis of the first principal stress coincides with the RD direction, α = 0◦ Therefore, according to Eqs. (20) and (21), the correction coefficients for sym- metric biaxial stress χeb and plane strain stress correction coefficient χps can be expressed as Fig. 10. Comparison of yield loci predicted by different models along (a) RD/TD, (b) 30◦/120◦ and (c) 45◦/135◦ principal stress directions of biaxial loadings at εP=0.002, εP=0.01, εP=0.02, and εP=0.03. The blue circles are the experimental data. The different curves are the results predicted by the new model and advanced yield models. Z. Mu et al. International Journal of Mechanical Sciences 282 (2024) 109640 11 where σexp eb(T− T)(ε p, 0), σexp eb(C− C)(ε p, 0) are the experimental values of equi- biaxial tension and equi-biaxial compassion yield stress along RD/TD direction at different equivalent plastic strains. σexp eb(T− C)(ε p, 0), σexp eb(C− T)(ε p, 0) are the experimental values of simple shear stress in the fourth and second quadrants at different equivalent plastic strains. σexp nps0− 1(εp, 0), and σexp nps90− 1(εp, 0) are the near-plane strain 1-stresses of biaxial tension along RD/TD direction. σexp nps0− 2(ε p, 0), and σexp nps90− 2(ε p, 0) are the near-plane strain 2-stresses of biaxial tension along RD/TD di- rection. The stress ratio under the near-plane strain stress state can be expressed γ1 = σexp nps0− 2/σexp nps0− 1 and γ2 = σexp nps90− 2/σexp nps90− 1. When the biaxial loading tests are carried out along the RD/TD di- rection (α = 0◦), two principal stresses under the plane stress state can be obtained according to Eqs. (20) and (21) ⎧ ⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎨ ⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎩ σ1= 1 χb(εP,0) ̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅ 1 [σ0(εP,0)]2 − 2hb ( εP,0 ) γ+ 1 [σ90(εP,0)]2 γ2 √ ,σ2=γσ1 |σ2|≤|σ1| σ1=γσ2,σ2= 1 χb(εP,0) ̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅ 1 [σ0(εP,0)]2 γ2 − 2hb ( εP,0 ) γ+ 1 [σ90(εP,0)]2 √ |σ2|≥|σ1| . (33) 4.2. The modified Hill48 model generalized to the principal stresses coordinate system When the principal stress direction of loading α ∕= 0◦, the yield loci in the principal stress coordinate system can be obtained by coordinate transformation. Then, the principal stress yield locus function including the direction angle of the principal stress axis can be expressed as k1σ2 1 + k2σ2 2 + k3σ1σ2 = 1, (34) where ⎧ ⎨ ⎩ k1 = Fsin4α + Gcos4α + 2(N − 2H)sin2αcos2α + H k2 = Gsin4α + Fcos4α + 2(N − 2H)sin2αcos2α + H k3 = 2 [ (F + G+ 4H − 2N)sin2αcos2α − H ] . (35) Then, the two principal stresses under BT in the α loading direction can be expressed as ⎧ ⎪⎪⎪⎪⎨ ⎪⎪⎪⎪⎩ σ1 = 1 ̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅ k1 + k2γ2 + k3γ √ , σ2 = γσ1 |σ2| ≤ |σ1| σ1 = γσ2, σ2 = 1 ̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅ k1γ2 + k2 + k3γ √ |σ2| ≥ |σ1| . (36) Under this condition, the stress correction coefficients of symmetric biaxial loadings for the arbitrary principal stress directions of loading can be expressed as Fig. 11. 3D yield surfaces predicted by the new modified Hill48 model at εp = 0.002. NPS0, NPS30 and NPS45 represent the near-plane strain stress along 0◦, 30◦, 45◦ to RD, respectively; The blue curve represents the normal plane, the orange curve represents the 30◦ plane, the burgundy curve represents the di- agonal plane, and the purple dashed curve represents the shear plane. χeb ( εP,α=0 ) = ⎧ ⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎨ ⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎩ 1 σexp eb(T− T)(ε P,0) ̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅ 1 / [σ0(εP)] 2 + 1 / [σ90(εP)] 2 − 2hb(εP, reb,0) √ (σ1 > 0andσ2 > 0) 1 σexp eb(C− C)(ε P, 0) ̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅ 1 / [σ0(εP)] 2 + 1 / [σ90(εP)] 2 − 2hb(εP, reb, 0) √ (σ1 < 0andσ2 < 0) 1 σexp eb(T− C)(ε P,0) ̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅ 1 / [σ0(εP)] 2 + 1 / [σ90(εP)] 2 + 2hb(εP, reb,0) √ (σ1 > 0andσ2 < 0) 1 σexp eb(C− T)(ε P,0) ̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅ 1 / [σ0(εP)] 2 + 1 / [σ90(εP)] 2 + 2hb(εP, reb,0) √ (σ1 < 0andσ2 > 0) , (31) χps ( εP, α=0 ) = { 1 σexp nps0− x(εP, 0) ̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅ 1 [σ0(εP,0)]2 − 2γ1hb ( εP, rps,0 ) + 1 [σ90(εP, 0)]2 γ21 √ NPS ∈ 0∘ 1 σexp nps90− x(εP, 0) ̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅ 1 [σ0(εP, 0)]2 − 2hb ( εP, rps, 0 ) γ2 + 1 [σ90(εP,0)]2 γ22 √ NPS ∈ 90∘ , (32) Z. Mu et al. International Journal of Mechanical Sciences 282 (2024) 109640 12 χα eb = ⎧ ⎪⎪⎪⎪⎨ ⎪⎪⎪⎪⎩ 1 σexp eb(T− T) ̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅ K1 + K2 + K3 √ (σ1 > 0 and σ2 > 0) 1 σexp eb(C− C) ̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅ K1 + K2 + K3 √ (σ1 < 0 and σ2 < 0) . (37) The stress correction coefficients for the near-plane strain state can be calculated as ⎧ ⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎨ ⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎩ χα ps = 1 σexp∼α nps− x ̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅ K1 + ( γα ps )2 K2 + γα psK3 √ χ90+α ps = 1 σexp∼(90+α) nps− x ̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅ K1 + K2 ( γ90+α ps )2 + K3 γ90+α ps √ √ √ √ , (38) where ⎧ ⎨ ⎩ K1 = fsin4α + gcos4α + 2(n − 2h)sin2αcos2α + h K2 = gsin4α + fcos4α + 2(n − 2h)sin2αcos2α + h K3 = 2 [ (f + g + 4h − 2n)sin2αcos2α − h ] , (39) { h = h(α) f = 1 / σ2 90 − h g = 1 / σ2 0 − h n = 2 / σ2 45 − (f + g) / 2 , (40) where σexp∼α eb(T− T), σexp∼α eb(C− C) are the experimental values of equi-biaxial ten- sion and equi-biaxial compassion yield stress along α principal stress direction of loading. The stress ratio under the near-plane strain stress state can be expressed γα ps = σexp∼α nps− 2/σexp∼α nps0− 1 and γ90+α ps = σexp∼(90+α) nps− 2 /σexp∼(90+α) nps0− 1 . σexp∼α nps0− 1, and σexp∼(90+α) nps− 1 are the near-plane strain 1- Fig. 12. Variation of the slope of biaxial tensile plastic strain curves with loading directions for different levels of EPS under (a) PSD: RD/TD, (b) PSD: 30◦/120◦ and (c) PSD: 45◦/135◦ ξ12 = ε2/ε1 at 0◦ < φ < 45◦; ξ12 = ε1/ε2 at 45◦ < φ < 90◦ Z. Mu et al. International Journal of Mechanical Sciences 282 (2024) 109640 13 stresses of biaxial tension along α and 90 + α direction respectively. σexp∼α nps0− 2, and σexp∼(90+α) nps− 2 are the near-plane strain 2-stresses of biaxial tension along α and 90 + α direction respectively. According to the above establishment process, compared with the original Hill48 model, the main improvements of the modified Hill48 model are as follows: (1) The decoupling between anisotropic parame- ters and stress state is realized by introducing the principal stress di- rection of loading, and the parameters calibration can be carried out based on the stress state data of uniaxial loadings and biaxial loadings respectively. The material parameters associated with in-plane and out- of-plane anisotropy are independent. (2) By performing quadratic interpolation between UT, NPS, and EBT, the flexibility of the model in predicting NPS stress has been significantly improved. This method is also friendly to the materials with asymmetric yield locus shape, including tension-compression asymmetry and distortion hardening asymmetry. (3) The original Hill48 model is unable to calibrate aniso- tropic parameters using the r-value of biaxial tension, resulting in its inability to accurately capture the direction of plastic flow such as EBT stress state. By further introducing the biaxial tensile plastic strain ratio in the modified model, the limitation of the original Hill48 model in the prediction of biaxial tensile plastic flow direction is effectively broken, and it can adapt to different principal stress loading directions. 4.3. Calibration and verification 4.3.1. Yield loci In this section, the yield loci along three principal stress directions of loading are investigated to validate the proposed general framework. For comparison, some advanced yield models Yld2000–2D, Poly4- &Hosford [56], Poly6 [85], and Hou2022-CQI [59] models are also considered. As shown in Fig. 10, the modified Hill48 model significantly im- proves the prediction accuracy of biaxial loading stresses under EBT and NPS, and can capture the evolution of yield loci. In addition, the anisotropic parameters can be calibrated using the experimental data of different principal stress directions of loading, resulting in the modified Hill48 model can accurately capture the change of yield loci with the principal stress directions of loading. However, because other yield models only use the biaxial loading data in the RD/TD direction, they can describe the principal plane yield loci in the RD/TD direction well. Unfortunately, there are significant deviations in the prediction of EBT Fig. 13. Comparison of the experimental and calculated directions of plastic strain rate by the newly modified Hill48 model along (a) RD/TD, (b) 30◦/120◦, and (c) 45◦/135◦ directions. (d) Comparison of the directions of plastic strain rate obtained from three principal stress directions of loading at εP=0.02. Experimental data are indicated by symbols (labeled as Exp.). The various types of curves represent the predicted results of the new model at different EPS (labeled as Pre.). Z. Mu et al. International Journal of Mechanical Sciences 282 (2024) 109640 14 and NPS stresses of the yield loci along the 30◦/120◦ and 45◦/135◦ loading directions. The Yld2000–2D model is based on linear transformations of the deviatoric stress tensor, while the Poly6 model utilizes homogeneous polynomials. The Poly4&Hosford and Hou2022-CQI models are devel- oped on the coupling multiplication of two additional terms. In these models, the material parameters associated with in-plane and out-of- plane anisotropy are interdependent rather than independent when considering the full stress state. Consequently, they fail to account for out-of-plane anisotropy, resulting in a degraded representation of this aspect. The yield surface in three-dimensional space predicted by the modified Hill48 model is shown in Fig. 11. Compared with other yield models, the modified model can take the NPS0, NPS30 and NPS45 as constraint points, which further improves the prediction accuracy of yield surface. Therefore, calibrating and verifying yield models solely based on the principal stress yield loci in the RD/TD direction is insuf- ficient. To improve the prediction of yield surface anisotropy under the full stress state, it is essential to develop a method that allows for independent parameter identification. This approach would ensure a more accurate and comprehensive representation of anisotropic behavior in sheet metals. 4.3.2. Description of plastic flow under different principal stress directions of loading According to the biaxial tension plastic strain given in Fig. 8 under different stress states, the slope of the plastic strain curves ξ12 can be calculated. Furthermore, the relationship curves between ξ12 and prin- cipal stress ratios at different EPS levels as shown in Fig. 12 can be ob- tained. The change of ξ12 with principal stress ratio approximates the shape of a quadratic curve. The variation of ξ12 with the EPS is relatively small under specific principal stress directions of loading, while there are significant differences in ξ12 under different principal stress di- rections of loading. This is the direct reason for the significant differ- ences in plastic flows under different principal stress directions of loading, as shown in Fig. 9(d). According to the new yield model and its parameters calibration method, the coefficients in Eq. (30) about ξ12 can be calculated by the Fig. 14. Comparison of the predicted directions of plastic strain rate by various yield models at εp = 0.002(a) RD/TD, (b) 30◦/120◦, and (c) 45◦/135◦ directions. Experimental data are indicated by black circles (labeled as Exp.). The various types of curves represent the predicted results of different yield models. The new model significantly improves the prediction accuracy of plastic flow direction under EBT and NPS stress states for different principal stress directions of loading. Z. Mu et al. International Journal of Mechanical Sciences 282 (2024) 109640 15 ratio of principal strain increment under UT (ξu), NPS (ξps), and the EBT (ξeb). Then, the anisotropic parameters under biaxial loading stress states can be obtained by substituting Eq. (30) into Eqs. (25), (28) and (29). Further, according to the calculation method of the biaxial tensile plastic flow, the comparison of βbetween the prediction and the experimental results are shown in Fig. 13. The limitation of original Hill48 model being unable to consider the plastic flow of biaxial loadings was solved by introducing the ratio of principal strain increment. In addition, the modified Hill48 model can not only accurately describe the plastic flow direction of the EBT stress, but also capture that of the NPS stress state. The plastic strain flow under different principal stress directions of loading is also well predicted because the newly modified model can cover the variation of principal stress directions. Compared with the advanced yield models developed in recent years, the new modified model also shows excellent predictive ability, espe- cially in near-plane strain stress state, as shown in Fig. 14. This advan- tage will certainly further improve the accuracy for describing the plastic flow behavior of sheet metal stamping. 5. Conclusions To reveal the evolution of the anisotropic yield locus and plastic flow in sheet metals under biaxial tensile directions beyond the traditional RD/TD anisotropic directions, the principal stress yield function of the Hill48 model was reconstructed by incorporating principal stress di- rections of loading in relation to RD/TD. Theoretical analysis demon- strated that the diversity of principal stress yield loci depends not only on anisotropic parameters but also on the principal stress direction of loading. Based on this, biaxial tensile tests of DP590 sheet metal were con- ducted using cruciform specimens cut along different principal stress directions of loading. Clear differences were observed in EBT stresses, NPS stresses, and plastic strain flow directions under biaxial tension along the RD/TD, 30◦/120◦, and 45◦/135◦ directions. To predict the evolution and rotation of the yielding ellipse for these various biaxial loading directions, a new modified Hill48 model under AFR was pro- posed. This model addressed the significant shortcoming of traditional yield models, namely the inability to calibrate anisotropic parameters with stress of deformation data from different principal stress directions of loading. Comparisons with several advanced yield models, including Yld2000–2D, Poly4&Hosford, Poly6, and Hou2022-CQI, showed that EBT and NPS stresses under RD/TD, 30◦/120◦, and 45◦/135◦ directions were accurately captured by introducing the general principal stress direction of loading and utilizing a quadratic interpolation method. This approach also solved the issue of degraded representation in constitutive models when accounting for out-of-plane anisotropy. For the plastic flow direction, the modified Hill48 model overcame the limitation of the original Hill48 model, which was unable to account for plastic flow under biaxial loadings, by introducing the ratio of principal strain increment. The prediction accuracy of plastic flow direction under EBT and NPS stress states in directions other than the RD/TD loading di- rections was significantly improved. Evaluating yield models for predicting anisotropic behavior should not be limited to the normal plane and shear plane yield loci. It is essential to verify the yield loci under other principal stress axes. The calibration strategy that considers the principal stress direction of loading and the discretization of anisotropic parameters offers a novel approach for developing future anisotropic constitutive models ac- counting for full stress states anisotropy. CRediT authorship contribution statement Zhenkai Mu: Writing – original draft, Methodology, Investigation, Conceptualization. Jiale Liu:Writing – review & editing, Methodology, Investigation. Wei Wang: Writing – review & editing, Validation, Formal analysis. Xuerui Dai: Software, Investigation. Shibo Ma: Su- pervision, Project administration, Funding acquisition. Yong Hou: Writing – review & editing, Validation, Project administration, Funding acquisition. Declaration of competing interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. Data availability Data will be made available on request. Acknowledgments The authors would like to acknowledge the financial support for this research provided through National Natural Science Foundation of China (no. 52205353), Natural Science Foundation of Hebei Province (no. E2021208025), and Hebei Province Innovation Capability Improvement Plan (no. 225A2201D). Yong Hou acknowledges the support of the research fellowship from Alexander von Humboldt Foundation. Appendix A. Experimental stress-strain curves and plastic flow direction of biaxial tensile under different principal stress directions of loading The true stress-strain curves along two biaxial tensile direction for three principal stress directions of loading were calculated as shown in Fig. A1. Based on the biaxial tensile strain data measured using DIC, the variation curves of the plastic flow direction with the EPS were calculated as shown in Fig. A2. Z. Mu et al. International Journal of Mechanical Sciences 282 (2024) 109640 16 Fig. A1. Biaxial tensile true stress-strain curves with seven load ratios in different principal stress axis directions of loading: (a-b) RD/TD; (c-d) 30◦/120◦; (e-f) 45◦/ 135◦. 4:0 represent the uniaxial tension along the loading directions. Z. Mu et al. International Journal of Mechanical Sciences 282 (2024) 109640 17 Fig. A2. Variation of the direction of plastic strain rate with EPS under different principal stress directions of loading: (a) RD/TD; (b) 30◦/120◦; (c) 45◦/135◦. Different colored curves represent different stress states. The direction of plastic flow fluctuates greatly near the initial yield and becomes saturated with the increase of deformation. Appendix B. The prediction of plastic flow under biaxial loadings The prediction of plastic flow under biaxial loadings is β = arctan(dε2 / dε1), (B1) where, the ratio of principal strain increment is defined as ξ12 = dε2/dε1. (B2) According to the yield function and the plastic flow rule, there is ξ12 = 2(F + H)σ2 − 2Hσ1 2(G+ H)σ1 − 2Hσ2 = σ2 / σ2 90 − hbσ1 σ1 / σ2 0 − hbσ2 . (B3) Then, the hb can be obtained as hb(εp, α) = σ1ξ12 / σ2 0 − σ2 / σ2 90 ξ12σ2 − σ1 . (B4) Let γ = σ2/σ1, the Eq. (B4) can be simplified as hb(εp, α) = ξ12 / σ2 0 − γ / σ2 90 γξ12 − 1 . (B5) Appendix C. Convexity analysis of the modified Hill48 model (1) Convexity analysis based on EBT stress correction Z. Mu et al. International Journal of Mechanical Sciences 282 (2024) 109640 18 According to the modified Hill48 model, the experimental EBT and NPS stresses will directly determine the shape of the yield surface. In order to ensure the convexity at the EBT stress position of yield locus, it is necessary to determine the maximum range of the NPS stress when the EBT stress is given. According to the anisotropic parameters of the modified model, there is a corresponding relationship between the EBT stress and the r-value. That is, combining Eqs. (34) and (35), the equal biaxial tensile stress can be expressed as σb = 1 ̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅ k1 + k2 + k3 √ = 1 ̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅ 1 / σ2 0 + 1 / σ2 90 − 2h̃ √ , (C1) where h̃ ( εP) = 1 1+ r̃ {( sin2α [σ90(εP)] 2 + cos2α [σ0(εP)] 2 ) r̃ − ( 4 [σ45(εP)] 2 − 2 [σ0(εP)] 2 − 2 [σ90(εP)] 2 ) sin2αcos2α } . (C2) By combining Eqs. (C1) and (C2), the r̃ corresponding to experimental EBT stress can be obtained. It is used to determine the anisotropic pa- rameters that can accurately predict EBT stress. r̃ = ( 1 / σ2 0 + 1 / σ2 90 − 1 / σ2 b ) + 2 ( 4 / σ2 45 − 2 / σ2 90 − 1 / σ2 0 ) sin2αcos2α 2 ( sin2α / σ2 90 + cos2α / σ2 0 ) − ( 1 / σ2 0 + 1 / σ2 90 − 1 / σ2 b ) . (C3) Furthermore, according to Eq. (31), it can be concluded that χb(α) = 1. Combining Eq. (25), and when α = 0, there is ⎧ ⎪⎪⎪⎨ ⎪⎪⎪⎩ H̃b ( εP) = 1 [σ0(εP)] 2 r̃ 1+ r̃ G̃b ( εP) = 1 [σ0(εP)] 2 − H̃b F̃b ( εP) = 1 [σ90(εP)] 2 − H̃b Ñb ( εP) = 2 [σ45(εP)] 2 − G̃b + F̃b 2 . (C4) Then, the principal stress yield loci equation and its anisotropy parameters can be expressed as k1σ2 1 + k2σ2 2 + k3σ1σ2 = 1, (C5) where ⎧ ⎨ ⎩ k1 = F̃Bsin4α + G̃Bcos4α + 2(ÑB − 2H̃B)sin2αcos2α + H̃B k2 = G̃Bsin4α + F̃Bcos4α + 2(ÑB − 2H̃B)sin2αcos2α + H̃B k3 = 2 [ (F̃B + G̃B + 4H̃B − 2ÑB)sin2αcos2α − H̃B ] . (C6) Then, the yield locus with accurately EBT stress can be obtained with the above anisotropic parameters. Under this condition, the convexity of the yield locus can be analyzed based on the Hessian matrix Hessian = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ∂2f ∂σ2 1 ∂2f ∂σ1∂σ2 ∂2f ∂σ1∂σ12 ∂2f ∂σ2∂σ1 ∂2f ∂σ2 2 ∂2f ∂σ2∂σ12 ∂2f ∂σ12∂σ1 ∂2f ∂σ12∂σ2 ∂2f ∂σ2 12 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ = ⎡ ⎣ 2k1 k3 0 k3 2k2 0 0 0 0 ⎤ ⎦. (C7) The three leading principal minors should be positive ⎧ ⎪⎪⎨ ⎪⎪⎩ H1 = 2k1 H2 = 4k1k2 − k23 H3 = 0 . (C8) (i) The analysis of H1 According to Eq. (C6), there is H1 = 2k1 = 1 σ2 90 sin2α ( 1 − 2cos2α ) + 1 σ2 0 cos2α ( 1 − 2sin2α ) + 4 σ2 45 sin2αcos2α. (C9) The surfaces ofH1 under different anisotropy conditions are obtained by giving different normalized stress ratios, as shown in Fig. C1. When α = 0◦, H1>0. Although the relationship between H1 and normalized stress ratios (σ45/σ0 and σ90/σ0) varies with the principal stress directions of loading, H1> 0 is satisfied in the given range. Z. Mu et al. International Journal of Mechanical Sciences 282 (2024) 109640 19 Fig. C1. . Values of the principal minors H1calculated by Eq. (C9) under various stress states and principal stress directions of loading. H1 > 0 is satisfied under a wide stress state based on the r̃ and anisotropic parameters (Eq. (C4)). (ii) The analysis of H2 According to Eq. (C6), there is H2 = 1 σ2 0σ2 90 A+ ( a σ2 0 + a σ2 90 + 1 σ2 0σ2 45 + 1 σ2 90σ2 45 − 1 4σ4 0 − 1 4σ4 90 − 2a σ2 45 ) B − C2, (C10) whereA = 1 − 6sin2αcos2α,B = 4sin2αcos2α, C = r̃/ [ σ2 0(1 + r̃) ] . The ranges of r̃, σ90/σ0, and σ45/σ0 meet the convexity condition are determined as shown in Fig. C2. For the DP590 material studied in this paper, the relationship of r̃-value, σ90/σ0 and σ45/σ0 located in the region where convexity is satisfied. Fig. C2. . Values of the principal minors H2calculated by Eq. (C10) under various stress states and principal stress directions of loading. (a) α = 0◦(b) α = 15◦(c) α = 30◦(d) α = 45◦. The red curves represent H2 = 0; H2 < 0when situated below the designated red curve. Z. Mu et al. International Journal of Mechanical Sciences 282 (2024) 109640 20 (2) Convexity analysis based on NPS stress correction According to the quadratic function interpolation correction method established in this paper, the size of the plane strain stress will directly determine the shape of the yield surface. According to the yield locus function, the plane strain stress corresponding to the r̃ can be expressed as σ̃α nps− x = 1 ̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅ k1 + ( γα ps )2 k2 + γα psk3 √ . (C11) Then, the maximum ratio of plane strain stress to uniaxial tensile stress can be obtained as PTRmax = σ̃α nps− x σ0 . (C12) The PTRmax is the upper limit of plane strain stress that can meet the convexity condition. The variation law of PTRmax with EBT stress is shown in the Fig. C3. That is, beyond this value, concave will occur at the EBT stress point. Therefore, when the ratio of experimental plane strain stress to uniaxial tensile stress is greater than PTRmax, the maximum correction about the plane strain stress should take the corresponding PTRmax. In the Fig. C3, PTRmin is the lower limit of the plane strain stress, which is obtained according to the geometry-inspired numerical convex analysis (GINCA) approach proposed by Lou et al. [97]. Under the certain EBT stress, the PTR value is continuously reduced to determine the minimum PTR, which its convex can be satisfied in the whole yield locus. Then, through changing the value of σb/σ0 (σb/σ0=1, σb/σ0>1 and σb/σ0<1), then the relationship between PTRmin and σb/σ0 under a wider scope is obtained. For the yield loci on the different principal stress directions of loading for DP590, their near-plane strain stresses PTR under different EPS lie between the upper and lower limit of convexity requirements. Fig. C3. The limit values of PTR to meet the convexity of yield locus. The red curve represents the lower limit for convex, and the blue curve represents the upper limit for convex. The PTR of near-plane strain stresses are indicated by symbols under different EPS and PSD. References [1] Li Q, Zhang H, Chen F, Xu D, Sui D, Cui Z. Study on the plastic anisotropy of advanced high strength steel sheet: experiments and microstructure-based crystal plasticity modeling. Int J Mech Sci 2020;176:105569. https://doi.org/10.1016/j. ijmecsci.2020.105569. [2] Zhang W, Xu J. Advanced lightweight materials for Automobiles: a review. Mater Des 2022;221:110994. https://doi.org/10.1016/j.matdes.2022.110994. [3] Yang C, Shi B, Peng Y, Pan F. Loading path dependent distortional hardening of Mg alloys: experimental investigation and constitutive modeling on cruciform specimens. Int J Mech Sci 2019;160:282–97. https://doi.org/10.1016/j. ijmecsci.2019.06.046. [4] He J, Han G, Guo C. Non-associated anisotropic plasticity of metal sheets based on the distortional concept. Thin-Wall Struct 2021;161. https://doi.org/10.1016/j. tws.2021.107523. [5] Zheng L, Cheng C, Wan M, Meng B. Experimental characterization and theoretical modeling of size-dependent distortional hardening behavior of ultrathin metal sheets under multi-axial loading. Eur J Mech A/Solid 2022;92. https://doi.org/ 10.1016/j.euromechsol.2021.104461. [6] Han G, He J, Li S, Lin Z. Simple shear methodology for local structure–property relationships of sheet metals: state-of-the-art and open issues. Prog Mater Sci 2024; 143:1–57. https://doi.org/10.1016/j.pmatsci.2024.101266. [7] Habib SA, Khan AS, Gnäupel-Herold T, Lloyd JT, Schoenfeld SE. Anisotropy, tension-compression asymmetry and texture evolution of a rare-earth-containing magnesium alloy sheet, ZEK100, at different strain rates and temperatures: experiments and modeling. Int J Plast 2017;95:163–90. https://doi.org/10.1016/j. ijplas.2017.04.006. [8] Pai N, Samajdar I, Patra A. Microstructural and mechanistic insights into the Tension–Compression asymmetry of rapidly solidified Fe–Cr alloys: a phase field and strain gradient plasticity study. J Mech Phys Solid 2024;189:105695. https:// doi.org/10.1016/j.jmps.2024.105695. [9] Li Z, Yang H, Liu J, Liu F. An improved yield criterion characterizing the anisotropic and tension-compression asymmetric behavior of magnesium alloy. J Magnes Alloy 2021. https://doi.org/10.1016/j.jma.2021.05.005. [10] Park J, Hou Y, Min J, Hou Z, Han HN, He B, et al. Understanding plasticity in multiphase quenching & partitioning steels: insights from crystal plasticity with stress state-dependent martensitic transformation. Int J Plast 2024:104075. [11] Wang Z, Jiang B, Wu S, Liu W. Anisotropic tension-compression asymmetry in SLM 316L stainless steel. Int J Mech Sci 2023;246:108139. https://doi.org/10.1016/j. ijmecsci.2023.108139. [12] Wang X, Shi Z, Lin J. Experimental study and modelling of anisotropic behaviour of aluminium-lithium alloys in creep age forming. Int J Mech Sci 2023;260:108659. https://doi.org/10.1016/j.ijmecsci.2023.108659. [13] He J, Feng Y. Determining localized necking in polycrystalline sheet metals using the bifurcation phenomenon in strain evolution. Crystals 2023;13. https://doi.org/ 10.3390/cryst13020272. [14] Banabic D, Kami A, Comsa DS, Eyckens P. Developments of the Marciniak- Kuczynski model for sheet metal formability: a review. J Mater Process Technol 2021;287:116446. https://doi.org/10.1016/j.jmatprotec.2019.116446. [15] Pham QT, Islam MS, Sigvant M, Caro LP, Lee MG, Kim YS. Improvement of modified maximum force criterion for forming limit diagram prediction of sheet metal. Int J Solids Struct 2023;273:112264. https://doi.org/10.1016/j. ijsolstr.2023.112264. [16] Esmaeilpour R, Kim H, Park T, Pourboghrat F, Xu Z, Mohammed B, et al. Calibration of Barlat Yld2004-18P yield function using CPFEM and 3D RVE for the simulation of single point incremental forming (SPIF) of 7075-O aluminum sheet. Int J Mech Sci 2018;145:24–41. https://doi.org/10.1016/j.ijmecsci.2018.05.015. [17] Yoon SY, Lee SY, Barlat F. Numerical integration algorithm of updated homogeneous anisotropic hardening model through finite element framework. Z. Mu et al. International Journal of Mechanical Sciences 282 (2024) 109640 21 https://doi.org/10.1016/j.ijmecsci.2020.105569 https://doi.org/10.1016/j.ijmecsci.2020.105569 https://doi.org/10.1016/j.matdes.2022.110994 https://doi.org/10.1016/j.ijmecsci.2019.06.046 https://doi.org/10.1016/j.ijmecsci.2019.06.046 https://doi.org/10.1016/j.tws.2021.107523 https://doi.org/10.1016/j.tws.2021.107523 https://doi.org/10.1016/j.euromechsol.2021.104461 https://doi.org/10.1016/j.euromechsol.2021.104461 https://doi.org/10.1016/j.pmatsci.2024.101266 https://doi.org/10.1016/j.ijplas.2017.04.006 https://doi.org/10.1016/j.ijplas.2017.04.006 https://doi.org/10.1016/j.jmps.2024.105695 https://doi.org/10.1016/j.jmps.2024.105695 https://doi.org/10.1016/j.jma.2021.05.005 http://refhub.elsevier.com/S0020-7403(24)00681-7/sbref0010 http://refhub.elsevier.com/S0020-7403(24)00681-7/sbref0010 http://refhub.elsevier.com/S0020-7403(24)00681-7/sbref0010 https://doi.org/10.1016/j.ijmecsci.2023.108139 https://doi.org/10.1016/j.ijmecsci.2023.108139 https://doi.org/10.1016/j.ijmecsci.2023.108659 https://doi.org/10.3390/cryst13020272 https://doi.org/10.3390/cryst13020272 https://doi.org/10.1016/j.jmatprotec.2019.116446 https://doi.org/10.1016/j.ijsolstr.2023.112264 https://doi.org/10.1016/j.ijsolstr.2023.112264 https://doi.org/10.1016/j.ijmecsci.2018.05.015 Comput Method Appl Mech Eng 2020;372:113449. https://doi.org/10.1016/j. cma.2020.113449. [18] Yoon SY, Barlat F, Lee SY, Kim JH, Wi MS, Kim DJ. Finite element implementation of hydrostatic pressure-sensitive plasticity and its application to distortional hardening model and sheet metal forming simulations. J Mater Process Technol 2022;302. https://doi.org/10.1016/j.jmatprotec.2022.117494. [19] Li Z, Kang Q, Sui X, Xu X, Zhan L, Wang G. Analysis of the mechanism of orientations evolution during hot rolling and mechanical properties of TiBw/TA15 composites based on crystal plasticity finite element model. J Mater Sci Technol 2023;167:137–51. https://doi.org/10.1016/j.jmst.2023.05.023. [20] Hill R. A theory of the yielding and plastic flow of anisotropic metals. Proc R Soc Lond A 1948;313:509–29. [21] Hill R. Theoretical plasticity of textured aggregates. Math Proc Cambrid Philos Soc 1979;85:179–91. https://doi.org/10.1017/S0305004100055596. [22] Hill R. Constitutive modelling of orthotropic plasticity in sheet metals. J Mech Phys Solid 1990;38:405–17. https://doi.org/10.1016/0022-5096(90)90006-P. [23] Hu W. Characterization behavior and corresponding yield criterion of anisotropic sheet metals. Mater Sci Eng A 2003;345:139–44. https://doi.org/10.1016/S0921- 5093(02)00453-7. [24] Cardoso RPR, Adetoro OB. A generalisation of the Hill’s quadratic yield function for planar plastic anisotropy to consider loading direction. Int J Mech Sci 2017; 128–129:253–68. https://doi.org/10.1016/j.ijmecsci.2017.04.024. [25] Hosford. A generalized isotropic yield criterion. J Appl Mech Trans ASME 1972;39: 1172. https://doi.org/10.1115/1.3422872. [26] Barlat F, Lege DJ, Brem JC. A six-component yield function for anisotropic materials. Int J Plast 1991;7:693–712. https://doi.org/10.1016/0749-6419(91) 90052-Z. [27] Barlat F, Maeda Y, Chung K, Yanagawa M, Brem JC, Hayashida Y, et al. Yield function development for aluminum alloy sheets. J Mech Phys Solid 1997;45: 1727–63. https://doi.org/10.1016/S0022-5096(97)00034-3. [28] Barlat F, Brem JC, Yoon JW, Chung K, Dick RE, Lege DJ, et al. Plane stress yield function for aluminum alloy sheets - part 1: theory. Int J Plast 2003;19:1297–319. https://doi.org/10.1016/S0749-6419(02)00019-0. [29] Aretz H. A non-quadratic plane stress yield function for orthotropic sheet metals. J Mater Process Technol 2005;168:1–9. https://doi.org/10.1016/j. jmatprotec.2004.10.008. [30] Barlat F, Aretz H, Yoon JW, Karabin ME, Brem JC, Dick RE. Linear transfomation- based anisotropic yield functions. Int J Plast 2005;21:1009–39. https://doi.org/ 10.1016/j.ijplas.2004.06.004. [31] Aretz H, Barlat F. New convex yield functions for orthotropic metal plasticity. Int J Non Linear Mech 2013;51:97–111. https://doi.org/10.1016/j. ijnonlinmec.2012.12.007. [32] Lou Y, Yoon JW. Anisotropic yield function based on stress invariants for BCC and FCC metals and its extension to ductile fracture criterion. Int J Plast 2018;101: 125–55. https://doi.org/10.1016/j.ijplas.2017.10.012. [33] Lou Y, Zhang S, Yoon JW. A reduced Yld2004 function for modeling of anisotropic plastic deformation of metals under triaxial loading. Int J Mech Sci 2019;161–162: 105027. https://doi.org/10.1016/j.ijmecsci.2019.105027. [34] Banabic D, Kuwabara T, Balan T, Comsa DS, Julean D. Non-quadratic yield criterion for orthotropic sheet metals under plane-stress conditions. Int J Mech Sci 2003;45:797–811. https://doi.org/10.1016/S0020-7403(03)00139-5. [35] Banabic D, Aretz H, Comsa DS, Paraianu L. An improved analytical description of orthotropy in metallic sheets. Int J Plast 2005;21:493–512. https://doi.org/ 10.1016/j.ijplas.2004.04.003. [36] Comsa D, Banabic D. Plane-stress yield criterion for highly-anisotropic sheet metals. Numisheet 2008;0:43–8. 2008. [37] Soare S, Yoon JW, Cazacu O. On the use of homogeneous polynomials to develop anisotropic yield functions with applications to sheet forming. Int J Plast 2008;24: 915–44. https://doi.org/10.1016/j.ijplas.2007.07.016. [38] Soare S, Barlat F. Convex polynomial yield functions. J Mech Phys Solids 2010;58: 1804–18. https://doi.org/10.1016/j.jmps.2010.08.005. [39] Soare SC. A parameter identification scheme for the orthotropic Poly6 yield function satisfying the convexity condition. Eur J Mech A/Solids 2022;92:104467. https://doi.org/10.1016/j.euromechsol.2021.104467. [40] Soare SC. Bezier5YS and SHYqp: a general framework for generating data and for modeling symmetric and asymmetric orthotropic yield surfaces. Eur J Mech A/ Solids 2023;97. https://doi.org/10.1016/j.euromechsol.2022.104781. [41] Cazacu O, Barlat F. Generalization of Drucker’s yield criterion to orthotropy. Math Mech Solid 2001;6:613–30. https://doi.org/10.1177/108128650100600603. [42] Cazacu O, Barlat F. A criterion for description of anisotropy and yield differential effects in pressure-insensitive metals. Int J Plast 2004;20:2027–45. https://doi. org/10.1016/j.ijplas.2003.11.021. [43] Cazacu O, Plunkett B, Barlat F. Orthotropic yield criterion for hexagonal closed packed metals. Int J Plast 2006;22:1171–94. https://doi.org/10.1016/j. ijplas.2005.06.001. [44] Yoshida F, Hamasaki H, Uemori T. A user-friendly 3D yield function to describe anisotropy of steel sheets. Int J Plast 2013;45:119–39. https://doi.org/10.1016/j. ijplas.2013.01.010. [45] Lou Y, Huh H, Yoon JW. Consideration of strength differential effect in sheet metals with symmetric yield functions. Int J Mech Sci 2013;66:214–23. https:// doi.org/10.1016/j.ijmecsci.2012.11.010. [46] Yoon JW, Lou Y, Yoon J, Glazoff MV. Asymmetric yield function based on the stress invariants for pressure sensitive metals. Int J Plast 2014;56:184–202. https://doi. org/10.1016/j.ijplas.2013.11.008. [47] Hu Q, Li X, Han X, Li H, Chen J. A normalized stress invariant-based yield criterion: modeling and validation. Int J Plast 2017;99:248–73. https://doi.org/10.1016/j. ijplas.2017.09.010. [48] Cazacu O. New yield criteria for isotropic and textured metallic materials. Int J Solids Struct 2018;139–140:200–10. https://doi.org/10.1016/j. ijsolstr.2018.01.036. [49] Lou Y, Zhang S, Yoon JW. Strength modeling of sheet metals from shear to plane strain tension. Int J Plast 2020;134:102813. https://doi.org/10.1016/j. ijplas.2020.102813. [50] Stoughton TB. A non-associated flow rule for sheet metal forming. Int J Plast 2002; 18:687–714. https://doi.org/10.1016/S0749-6419(01)00053-5. [51] Kim JJ, Pham QT, Kim YS. Thinning prediction of hole-expansion test for DP980 sheet based on a non-associated flow rule. Int J Mech Sci 2021;191:106067. https://doi.org/10.1016/j.ijmecsci.2020.106067. [52] Mu Z, Zhao J, Meng Q, Sun H, Yu G. Anisotropic hardening and evolution of r- values for sheet metal based on evolving non-associated Hill48 model. Thin-Wall Struct 2022;171:108791. https://doi.org/10.1016/j.tws.2021.108791. [53] Hu Q, Yoon JW. Analytical description of an asymmetric yield function (Yoon2014) by considering anisotropic hardening under non-associated flow rule. Int J Plast 2021;140:102978. https://doi.org/10.1016/j.ijplas.2021.102978. [54] Lee EH, Stoughton TB, Yoon JW. A yield criterion through coupling of quadratic and non-quadratic functions for anisotropic hardening with non-associated flow rule. Int J Plast 2017;99:120–43. https://doi.org/10.1016/j.ijplas.2017.08.007. [55] Park N, Stoughton TB, Yoon JW. A criterion for general description of anisotropic hardening considering strength differential effect with non-associated flow rule. Int J Plast 2019;121:76–100. https://doi.org/10.1016/j.ijplas.2019.04.015. [56] Hu Q, Yoon JW, Manopulo N, Hora P. A coupled yield criterion for anisotropic hardening with analytical description under associated flow rule: modeling and validation. Int J Plast 2021;136:102882. https://doi.org/10.1016/j. ijplas.2020.102882. [57] Hou Y, Min J, Stoughton TB, Lin J, Carsley JE, Carlson BE. A non-quadratic pressure-sensitive constitutive model under non-associated flow rule with anisotropic hardening: modeling and validation. Int J Plast 2020;135:102808. https://doi.org/10.1016/j.ijplas.2020.102808. [58] Hou Y, Min J, Lee MG. Non-associated and Non-quadratic Characteristics in Plastic Anisotropy of Automotive Lightweight Sheet Metals. Automot Innov 2023;6: 364–78. https://doi.org/10.1007/s42154-023-00232-5. [59] Hou Y, Min J, Lin J, Lee MG. Modeling stress anisotropy, strength differential, and anisotropic hardening by coupling quadratic and stress-invariant-based yield functions under non-associated flow rule. Mech Mater 2022;174:104458. https:// doi.org/10.1016/j.mechmat.2022.104458. [60] Chen Z, Wang Y, Lou Y. User-friendly anisotropic hardening function with non- associated flow rule under the proportional loadings for BCC and FCC metals. Mech Mater 2022;165:104190. https://doi.org/10.1016/j.mechmat.2021.104190. [61] Hou Y, Min J, El-Aty AA, Han HN, Lee MG. A new anisotropic-asymmetric yield criterion covering wider stress states in sheet metal forming. Int J Plast 2023;166: 103653. https://doi.org/10.1016/j.ijplas.2023.103653. [62] Hou Y, Du K, Min J, Lee H, Lou Y, Park N, et al. A generalized, computationally versatile plasticity model framework - Part I : theory and verification focusing on tension ‒ compression asymmetry. Int J Plast 2023;171:103818. https://doi.org/ 10.1016/j.ijplas.2023.103818. [63] Lou Y, Yoon JW. Lode-dependent anisotropic-asymmetric yield function for isotropic and anisotropic hardening of pressure-insensitive materials. Part I: quadratic function under non-associated flow rule. Int J Plast 2023;166:103647. https://doi.org/10.1016/j.ijplas.2023.103647. [64] Wang S, Shang H, Han M, Zhou C, Chen Q, Lou Y. Lode-dependent Yld91 function for anisotropic-asymmetric hardening modeling of metals under non-associated flow rule. J Mater Process Technol 2024;325:118298. https://doi.org/10.1016/j. jmatprotec.2024.118298. [65] Zhou Y, Hu Q, Chen J. A concise analytical framework for describing asymmetric yield behavior based on the concept of shape functions. Int J Plast 2023;164. https://doi.org/10.1016/j.ijplas.2023.103593. [66] Mu Z, Zhao J, Meng Q, Zhang Y, Yu G. Limitation analysis of the Hill48 yield model and establishment of its modified model for planar plastic anisotropy. J Mater Process Technol 2022;299:117380. https://doi.org/10.1016/j. jmatprotec.2021.117380. [67] Zhang C, Lou Y. Characterization and modelling of evolving plasticity behaviour up to fracture for FCC and BCC metals. J Mater Process Technol 2023;317. https://doi. org/10.1016/j.jmatprotec.2023.117997. [68] Zhang C, Lou Y. Influences of the evolving plastic behavior of sheet metal on V- bending and springback analysis considering different stress states. Int J Plast 2024;173:103889. https://doi.org/10.1016/j.ijplas.2024.103889. [69] Hou Y, Min J, Guo N, Lin J, Carsley JE, Stoughton TB, et al. Investigation of evolving yield surfaces of dual-phase steels. J Mater Process Technol 2021;287: 116314. https://doi.org/10.1016/j.jmatprotec.2019.116314. [70] Hou Y, Min J, Guo N, Shen Y, Lin J. Evolving asymmetric yield surfaces of quenching and partitioning steels: characterization and modeling. J Mater Process Technol 2021;290:116979. https://doi.org/10.1016/j.jmatprotec.2020.116979. [71] Hou Y, Lee M-G, Lin J, Min J. Experimental characterization and modeling of complex anisotropic hardening in quenching and partitioning (Q&P) steel subject to biaxial non-proportional loadings. Int J Plast 2022;156:103347. https://doi.org/ 10.1016/j.ijplas.2022.103347. [72] Mu Z, Zhao J, Yu G, Huang X, Meng Q, Zhai R. Hardening model of anisotropic sheet metal during the diffuse instability necking stage of uniaxial tension. Thin- Walled Struct 2021;159:107198. https://doi.org/10.1016/j.tws.2020.107198. Z. Mu et al. International Journal of Mechanical Sciences 282 (2024) 109640 22 https://doi.org/10.1016/j.cma.2020.113449 https://doi.org/10.1016/j.cma.2020.113449 https://doi.org/10.1016/j.jmatprotec.2022.117494 https://doi.org/10.1016/j.jmst.2023.05.023 http://refhub.elsevier.com/S0020-7403(24)00681-7/sbref0020 http://refhub.elsevier.com/S0020-7403(24)00681-7/sbref0020 https://doi.org/10.1017/S0305004100055596 https://doi.org/10.1016/0022-5096(90)90006-P https://doi.org/10.1016/S0921-5093(02)00453-7 https://doi.org/10.1016/S0921-5093(02)00453-7 https://doi.org/10.1016/j.ijmecsci.2017.04.024 https://doi.org/10.1115/1.3422872 https://doi.org/10.1016/0749-6419(91)90052-Z https://doi.org/10.1016/0749-6419(91)90052-Z https://doi.org/10.1016/S0022-5096(97)00034-3 https://doi.org/10.1016/S0749-6419(02)00019-0 https://doi.org/10.1016/j.jmatprotec.2004.10.008 https://doi.org/10.1016/j.jmatprotec.2004.10.008 https://doi.org/10.1016/j.ijplas.2004.06.004 https://doi.org/10.1016/j.ijplas.2004.06.004 https://doi.org/10.1016/j.ijnonlinmec.2012.12.007 https://doi.org/10.1016/j.ijnonlinmec.2012.12.007 https://doi.org/10.1016/j.ijplas.2017.10.012 https://doi.org/10.1016/j.ijmecsci.2019.105027 https://doi.org/10.1016/S0020-7403(03)00139-5 https://doi.org/10.1016/j.ijplas.2004.04.003 https://doi.org/10.1016/j.ijplas.2004.04.003 http://refhub.elsevier.com/S0020-7403(24)00681-7/sbref0036 http://refhub.elsevier.com/S0020-7403(24)00681-7/sbref0036 https://doi.org/10.1016/j.ijplas.2007.07.016 https://doi.org/10.1016/j.jmps.2010.08.005 https://doi.org/10.1016/j.euromechsol.2021.104467 https://doi.org/10.1016/j.euromechsol.2022.104781 https://doi.org/10.1177/108128650100600603 https://doi.org/10.1016/j.ijplas.2003.11.021 https://doi.org/10.1016/j.ijplas.2003.11.021 https://doi.org/10.1016/j.ijplas.2005.06.001 https://doi.org/10.1016/j.ijplas.2005.06.001 https://doi.org/10.1016/j.ijplas.2013.01.010 https://doi.org/10.1016/j.ijplas.2013.01.010 https://doi.org/10.1016/j.ijmecsci.2012.11.010 https://doi.org/10.1016/j.ijmecsci.2012.11.010 https://doi.org/10.1016/j.ijplas.2013.11.008 https://doi.org/10.1016/j.ijplas.2013.11.008 https://doi.org/10.1016/j.ijplas.2017.09.010 https://doi.org/10.1016/j.ijplas.2017.09.010 https://doi.org/10.1016/j.ijsolstr.2018.01.036 https://doi.org/10.1016/j.ijsolstr.2018.01.036 https://doi.org/10.1016/j.ijplas.2020.102813 https://doi.org/10.1016/j.ijplas.2020.102813 https://doi.org/10.1016/S0749-6419(01)00053-5 https://doi.org/10.1016/j.ijmecsci.2020.106067 https://doi.org/10.1016/j.tws.2021.108791 https://doi.org/10.1016/j.ijplas.2021.102978 https://doi.org/10.1016/j.ijplas.2017.08.007 https://doi.org/10.1016/j.ijplas.2019.04.015 https://doi.org/10.1016/j.ijplas.2020.102882 https://doi.org/10.1016/j.ijplas.2020.102882 https://doi.org/10.1016/j.ijplas.2020.102808 https://doi.org/10.1007/s42154-023-00232-5 https://doi.org/10.1016/j.mechmat.2022.104458 https://doi.org/10.1016/j.mechmat.2022.104458 https://doi.org/10.1016/j.mechmat.2021.104190 https://doi.org/10.1016/j.ijplas.2023.103653 https://doi.org/10.1016/j.ijplas.2023.103818 https://doi.org/10.1016/j.ijplas.2023.103818 https://doi.org/10.1016/j.ijplas.2023.103647 https://doi.org/10.1016/j.jmatprotec.2024.118298 https://doi.org/10.1016/j.jmatprotec.2024.118298 https://doi.org/10.1016/j.ijplas.2023.103593 https://doi.org/10.1016/j.jmatprotec.2021.117380 https://doi.org/10.1016/j.jmatprotec.2021.117380 https://doi.org/10.1016/j.jmatprotec.2023.117997 https://doi.org/10.1016/j.jmatprotec.2023.117997 https://doi.org/10.1016/j.ijplas.2024.103889 https://doi.org/10.1016/j.jmatprotec.2019.116314 https://doi.org/10.1016/j.jmatprotec.2020.116979 https://doi.org/10.1016/j.ijplas.2022.103347 https://doi.org/10.1016/j.ijplas.2022.103347 https://doi.org/10.1016/j.tws.2020.107198 [73] Iftikhar CMA, Brahme A, Inal K, Khan AS. An evolution of subsequent yield loci under proportional and non-proportional loading path of ‘as-received’ extruded AZ31 magnesium alloy: experiments and CPFEM modeling. Int J Plast 2022;151: 103216. https://doi.org/10.1016/j.ijplas.2022.103216. [74] Barlat F, Gracio JJ, Lee MG, Rauch EF, Vincze G. An alternative to kinematic hardening in classical plasticity. Int J Plast 2011;27:1309–27. https://doi.org/ 10.1016/j.ijplas.2011.03.003. [75] Lee EH, Choi H, Stoughton TB, Yoon JW. Combined anisotropic and distortion hardening to describe directional response with Bauschinger effect. Int J Plast 2019;122:73–88. https://doi.org/10.1016/j.ijplas.2019.07.007. [76] Shi B, Peng Y, Yang C, Pan F, Cheng R, Peng Q. Loading path dependent distortional hardening of Mg alloys: experimental investigation and constitutive modeling. Int J Plast 2017;90:76–95. https://doi.org/10.1016/j. ijplas.2016.12.006. [77] Zhang W, Zhuang X, Zhang Y, Zhao Z. An enhanced François distortional yield model: theoretical framework and experimental validation. Int J Plast 2020;127: 102643. https://doi.org/10.1016/j.ijplas.2019.102643. [78] Raj A, Verma RK, Singh PK, Shamshoddin S, Biswas P, Narasimhan K. Experimental and numerical investigation of differential hardening of cold rolled steel sheet under non-proportional loading using biaxial tensile test. Int J Plast 2022;154: 103297. https://doi.org/10.1016/j.ijplas.2022.103297. [79] He J, Han G, Feng Y. Phase transformation and plastic behavior of QP steel sheets: transformation kinetics-informed modeling and forming limit prediction. Thin- Wall Struct 2022;173. https://doi.org/10.1016/j.tws.2022.108977. [80] Lee EH, Rubin MB. Eulerian constitutive equations for the coupled influences of anisotropic yielding, the Bauschinger effect and the strength-differential effect for plane stress. Int J Solid Struct 2022;241:111475. https://doi.org/10.1016/j. ijsolstr.2022.111475. [81] Plunkett B, Lebensohn RA, Cazacu O, Barlat F. Anisotropic yield function of hexagonal materials taking into account texture development and anisotropic hardening. Acta Mater 2006;54:4159–69. https://doi.org/10.1016/j. actamat.2006.05.009. [82] Yoshida F, Hamasaki H, Uemori T. Modeling of anisotropic hardening of sheet metals including description of the Bauschinger effect. Int J Plast 2015;75:170–88. https://doi.org/10.1016/j.ijplas.2015.02.004. [83] Lou Y, Zhang C, Zhang S, Yoon JW. A general yield function with differential and anisotropic hardening for strength modelling under various stress states with non- associated flow rule. Int J Plast 2022;158:103414. https://doi.org/10.1016/j. ijplas.2022.103414. [84] Stoughton TB, Yoon JW. Anisotropic hardening and non-associated flow in proportional loading of sheet metals. Int J Plast 2009;25:1777–817. https://doi. org/10.1016/j.ijplas.2009.02.003. [85] Hu Q, Yoon JW, Stoughton TB. Analytical determination of anisotropic parameters for Poly6 yield function. Int J Mech Sci 2021;201:106467. https://doi.org/ 10.1016/j.ijmecsci.2021.106467. [86] Du K, Huang S, Li X, Wang H, Zheng W, Yuan X. Evolution of yield behavior for AA6016-T4 and DP490 - Towards a systematic evaluation strategy for material models. Int J Plast 2022;154:103302. https://doi.org/10.1016/j. ijplas.2022.103302. [87] Jeong K, Lee K, Kwon D, Lee MG, Han HN. Parameter determination of anisotropic yield function using neural network-based indentation plastometry. Int J Mech Sci 2024;263:108776. https://doi.org/10.1016/j.ijmecsci.2023.108776. [88] Liu W, Li X, Liu M, Cui H, Huang J, Pang Y, et al. Virtual laboratory enabled constitutive modelling of dual phase steels. Int J Plast 2024;175. https://doi.org/ 10.1016/j.ijplas.2024.103930. [89] Wessel A, Morand L, Butz A, Helm D, Volk W. Machine learning-based sampling of virtual experiments within the full stress state. Int J Mech Sci 2024;275:109307. https://doi.org/10.1016/j.ijmecsci.2024.109307. [90] Cheng C, Wan M, Meng B, Zhao R, Han WP. Size effect on the yield behavior of metal foil under multiaxial stress states: experimental investigation and modelling. Int J Mech Sci 2019;151:760–71. https://doi.org/10.1016/j.ijmecsci.2018.12.031. [91] Zhang K, He Z, Zheng K, Yuan S. Experimental verification of anisotropic constitutive models under tension-tension and tension-compression stress states. Int J Mech Sci 2020;178:105618. https://doi.org/10.1016/j. ijmecsci.2020.105618. [92] Zhu H, He Z, Lin Y, Zheng K, Fan X, Yuan S. The development of a novel forming limit diagram under nonlinear loading paths in tube hydroforming. Int J Mech Sci 2020;172. https://doi.org/10.1016/j.ijmecsci.2019.105392. [93] He Z, Zhang K, Zhu H, Lin Y, Fu MW, Yuan S. An anisotropic constitutive model for forming of aluminum tubes under both biaxial tension and pure shear stress states. Int J Plast 2022;152:103259. https://doi.org/10.1016/j.ijplas.2022.103259. [94] Hu Q, Yoon JW, Chen J. Analytically described polynomial yield criterion by considering both plane strain and pure shear states. Int J Plast 2023;162:103514. https://doi.org/10.1016/j.ijplas.2022.103514. [95] Du K, Huang S, Hou Y, Wang H, Wang Y, Zheng W, et al. Characterization of the asymmetric evolving yield and flow of 6016-T4 aluminum alloy and DP490 steel. J Mater Sci Technol 2023;133:209–29. https://doi.org/10.1016/j. jmst.2022.05.040. [96] Hu Q, Maier L, Nishiwaki T, Hartmann C, Volk W, Yoon JW. User friendly FE Formulation for anisotropic distortional hardening model based on non-associated flow plasticity and its application to springback prediction. Thin-Wall Struct 2024; 202. https://doi.org/10.1016/j.tws.2024.112142. [97] Lou Y, Zhang C, Wu P, Yoon JW. New geometry-inspired numerical convex analysis method for yield functions under isotropic and anisotropic hardenings. Int J Solid Struct 2024;286–287:112582. Z. Mu et al. International Journal of Mechanical Sciences 282 (2024) 109640 23 https://doi.org/10.1016/j.ijplas.2022.103216 https://doi.org/10.1016/j.ijplas.2011.03.003 https://doi.org/10.1016/j.ijplas.2011.03.003 https://doi.org/10.1016/j.ijplas.2019.07.007 https://doi.org/10.1016/j.ijplas.2016.12.006 https://doi.org/10.1016/j.ijplas.2016.12.006 https://doi.org/10.1016/j.ijplas.2019.102643 https://doi.org/10.1016/j.ijplas.2022.103297 https://doi.org/10.1016/j.tws.2022.108977 https://doi.org/10.1016/j.ijsolstr.2022.111475 https://doi.org/10.1016/j.ijsolstr.2022.111475 https://doi.org/10.1016/j.actamat.2006.05.009 https://doi.org/10.1016/j.actamat.2006.05.009 https://doi.org/10.1016/j.ijplas.2015.02.004 https://doi.org/10.1016/j.ijplas.2022.103414 https://doi.org/10.1016/j.ijplas.2022.103414 https://doi.org/10.1016/j.ijplas.2009.02.003 https://doi.org/10.1016/j.ijplas.2009.02.003 https://doi.org/10.1016/j.ijmecsci.2021.106467 https://doi.org/10.1016/j.ijmecsci.2021.106467 https://doi.org/10.1016/j.ijplas.2022.103302 https://doi.org/10.1016/j.ijplas.2022.103302 https://doi.org/10.1016/j.ijmecsci.2023.108776 https://doi.org/10.1016/j.ijplas.2024.103930 https://doi.org/10.1016/j.ijplas.2024.103930 https://doi.org/10.1016/j.ijmecsci.2024.109307 https://doi.org/10.1016/j.ijmecsci.2018.12.031 https://doi.org/10.1016/j.ijmecsci.2020.105618 https://doi.org/10.1016/j.ijmecsci.2020.105618 https://doi.org/10.1016/j.ijmecsci.2019.105392 https://doi.org/10.1016/j.ijplas.2022.103259 https://doi.org/10.1016/j.ijplas.2022.103514 https://doi.org/10.1016/j.jmst.2022.05.040 https://doi.org/10.1016/j.jmst.2022.05.040 https://doi.org/10.1016/j.tws.2024.112142 http://refhub.elsevier.com/S0020-7403(24)00681-7/sbref0097 http://refhub.elsevier.com/S0020-7403(24)00681-7/sbref0097 http://refhub.elsevier.com/S0020-7403(24)00681-7/sbref0097 Characterization and modeling of biaxial plastic anisotropy in metallic sheets 1 Introduction 2 Effect of principal stress direction of loading on yield loci 2.1 Yield function expressed with principal stresses 2.2 Diversity of yield locus in the principal plane 3 Mechanical experiments for different principal stress directions of loading 3.1 Uniaxial tension 3.2 Biaxial tension 3.2.1 Yield locus 3.2.2 Plastic flow 4 A new yield model coupled with the principal stress loading direction 4.1 Modified Hill48 model in anisotropic principal axis coordinate system 4.1.1 Uniaxial loading stresses state 4.1.2 Biaxial loading stress state 4.2 The modified Hill48 model generalized to the principal stresses coordinate system 4.3 Calibration and verification 4.3.1 Yield loci 4.3.2 Description of plastic flow under different principal stress directions of loading 5 Conclusions CRediT authorship contribution statement Declaration of competing interest Data availability Acknowledgments Appendix A Experimental stress-strain curves and plastic flow direction of biaxial tensile under different principal stress ... Appendix B The prediction of plastic flow under biaxial loadings Appendix C Convexity analysis of the modified Hill48 model References