Stress-Invariants-Based Anisotropic Yield Functions and Its Application to Sheet Metal Plasticity

: The yield criterion, or so-called yield function, plays an important role in the study of the plastic working of a sheet because it governs the plastic deformation properties of the sheet during the plastic-forming process. In this paper, we propose a novel anisotropic yield function useful for describing the plastic behavior of various anisotropic sheets. The proposed yield function includes the anisotropic version of the second stress invariant J 2 and the third stress invariant J 3 . The proposed yield function can explain the anisotropic plastic behavior of various sheets by introducing the parameters α and β and also exhibits both symmetrical and asymmetrical yield surfaces. The parameters included in the proposed model were determined with an optimization algorithm from uniaxial and biaxial experimental data under a proportional loading path. In this study, the validity of the proposed anisotropic yield function was veriﬁed by comparing the yield surface shape, normalized uniaxial yield stress value, and Lankford anisotropic coefﬁcient R-value derived from the experimental results. Applications of the proposed anisotropic yield functions to an aluminum sheet showed symmetrical yielding behavior and, to pure titanium sheets, showed asymmetric yielding behavior; thus, it was shown that the yield curve and yield behavior of various types of sheet materials can be predicted reasonably by using the proposed new yield anisotropic function.


Introduction
Sheet metal forming is being applied to the manufacturing of various products such as interior and exterior panels of automobiles, aircraft fuselages, and beverage cans.As a mass production process, sheet material processing technology is one of the technologies for manufacturing products with complex shapes from thin sheets by using the plastic deformation properties of the sheet material.
In this process, it is necessary to optimize the manufacturing process to secure product quality, increase productivity, shorten the product development period, and reduce the manufacturing cost.To this end, some trial error methods based on the experience of skilled technicians have been traditionally used in manufacturing sites, but virtual simulation technology based on finite element analysis (FEA), which is a computational simulation numerical analysis technology, has recently been widely adopted [1][2][3].
To increase the accuracy of the sheet metal forming process simulation using this FEA, it is necessary to select a suitable finite element type, to discretize a sufficiently small finite element mesh, and to introduce suitable boundary conditions and an accurate plastic deformation model of the material.In particular, in order to understand the plastic deformation behavior of a material well, an accurate understanding of work hardening, the plastic yield criterion, and the flow rule is basically required.As part of this, in order to highly predict the plastic deformation behavior of sheet materials through FEA, various types of work hardening laws [4][5][6], plastic flow laws [7,8], and yield functions have been proposed by many researchers.In addition, combinations of various components have been attempted to improve the accuracy of FEA [9][10][11][12].
Among the plastic deformation models of materials, the anisotropic yield criterion serves to well describe the anisotropic yield behavior of cold-rolled anisotropic sheets, and it is one of the important factors in FEA analyses because it determines the elastic and plastic states of the sheets subjected to various stress ratios.
Various anisotropic yield conditions have been proposed by many scholars in previous studies.For isotropic materials, since the yield function must be invariant for all orthogonal transformations, it is possible to utilize the second invariant J 2 and the third invariant J 3 for the deviating stress as first proposed by Prager in the yield condition [13]: These J 2 and J 3 are expressed in terms of the deviation stress tensor as follows: where σ k and s k (k = 1, . . .3) represent the main values of the Cauchy stress tensor and the deviation stress tensor, respectively.Tresca [14] first proposed an isotropic yield function by considering that the yielding occurs when the maximum shear stress reaches a constant value.Von Mises [15] proposed the following yield function, which is widely used for isotropic materials: where k is the yield stress in pure shear, calculated as k = σ T / √ 3 , and σ T is the yield stress in uniaxial tension.Drucker [16] proposed a yield function that includes J 3 as well as the stress invariant J 2 as bellow: A notable property of Drucker's yield criterion is that the yield surface lies between the Tresca and Von Mises yield surfaces.Drucker's yield function is known to predict the plastic behavior of aluminum sheets particularly well.
Various types of anisotropic yield functions have been proposed to explain the intrinsic plastic anisotropy properties of cold-rolled steel sheet materials (Hill-48 (1948) [17], Hill-79 (1979) [18], Hosford (1979) [19], Yld2000-2d (2001) [20], Cazacu-Barlat (2001) [21], Yld2004-18p (2005) [22], and Banabic (2010 and 2020) [23,24]).These yield functions can be classified into two groups: one is the second-order yield functions and the other is the non-second-order yield functions.Von Mises and Hill-48 (1948) are examples of quadratic yield functions, and the other examples mentioned above are nonquadratic yield functions.The Cazacu-Barlat (2001) yield function, one of the nonquadratic yield functions, was developed by extending the isotropic Drucker yield criterion.They replaced the isotropic versions of J 2 and J 3 with anisotropic versions of J 0 2 and J 0 3 , respectively: where c is a constant.The Cazacu-Barlat yield function showed reasonable predictions for both the yield stress and Lankford anisotropic coefficient R-value in each direction for the aluminum sheet.However, since this yield function uses only one function for the two terms J 2 and J 3 , it is not suitable for explaining the anisotropy of the sheet material well.[25] On the other hand, the Yld2004-18p yield function can accurately predict the anisotropic plastic deformation behavior, but it requires many experiments to identify the In addition, the validity of the nonassociated flow rule (non-AFR), which introduced the yield function and other plastic potential functions, was also discussed.Meanwhile, Karafillis and Boyce [27] first introduced the concept of anisotropic linear transformation to construct anisotropic yield functions.
Cazacu (2018) [28] proposed the following yield function and equivalent stress for textured metallic materials with isotropic invariants of J 0 2 and J 0 3 : 1/8 (7) In the case of HCP (hexagonal close-packed) metals such as titanium and magnesium alloys, the plastic deformation caused by the twin mechanism ultimately showed a distinct difference in yield strength under tensile and compressive loading conditions.This phenomenon is called the strength difference effect (SD effect) [29][30][31][32][33][34][35][36].The asymmetric plastic behavior of the yield stress in tension and compression can be easily confirmed from the shape of the yield surface.In order to confirm the deformation behavior of these materials, experiments must be performed under several deformation modes such as uniaxial and biaxial tension as well as uniaxial compression.[37] Cazacu and Barlat (2004) [38] proposed the following anisotropic yield function based on the stress invariants J 2 and J 3 .In this yield function, the SD effect exhibited in HCP materials was explained by expressing the second and third invariants of the stress deviation as odd functions of the stress: In order to express both anisotropy and tension/compression asymmetry for pressureinsensitive metals, Cazacu et al. (2006) [39] proposed another yield function (CPB'06) using a quadratic linear transformation for stress deviation.Additionally, Gao (2011) [40] proposed a yield function combining the effects of all three of the first invariant I 1 , the second invariant J 2 , and the third invariant J 3 as follows: Meanwhile, Khan et al., (2012) [41] proposed an anisotropic yield condition including a normalized third invariant to explain the SD effect based on the experimental investigation of Ti-6Al-4V alloys in tension and compression tests.Yoon (2014) [42] proposed the following pressure-dependent yield function based on three invariants: (10) In this study, we propose a new isotropic yield function combining the second and third anisotropic invariants (i.e., J 0 2 and J 0 3 ) with two parameters α and β and explain the effects of these parameters in this new yield function.In addition, in order to confirm the validity of the proposed yield condition, the yield curve shape was predicted for AA6016-T4, AA2090-T3, and pure titanium sheets.In addition, after taking a specimen every 15 • from the rolling direction (RD) to the transverse direction (TD), uniaxial tension was performed to compare the anisotropy (normalized yield stress and R-value of Lankford coefficient of anisotropy) with experimental data.

Proposed Model 2.1. Proposed Anisotropic Yield Function
In this study, we propose the following yield function for isotropic sheet material (named the Kim-Van yield function, simply KV '12), which includes J 2 and J 3 as follows and does not depend on the hydrostatic pressure magnitude: where J 2 and J 3 are the invariants, k is the yield stress at pure shear, and α and β j are the material parameters to define the shape of the yield function.Here, when β 1 = 0 at j = 1 and β 2 = 0, the yield function of the Equation ( 11) reduces to This proposed yield function represents an asymmetric yield function (named KV'12A) because the third term of the above equation is an odd function of the stress and thus captures the SD effect.This asymmetric yield function is suitable for HCP materials such as titanium and magnesium.
Meanwhile, when j = 2 and when β 1 = 0 and β 2 = 0, the yield function of the Equation ( 11) reduces to The proposed yield function represents a symmetric yield function (named KV'12S) because all terms of the above equation are an even function of the stress and thus cannot represent the SD effect.This symmetrical yield function is suitable for BCC (Body-Centered Cubic) and FCC (Face-Centered Cubic) materials such as steel and aluminum.
This equation can be simplified with the following equation Mathematically, α and β 2 in Equation ( 13) have a mutually dependent relationship so that α and β 2 may be expressed as α = α 2 1 and β 2 = 2α 1 , respectively.However, the α and β 2 in Equation ( 13) and the α and β 1 in Equation (12) are assigned independently of each other in this study so that the shape of various yield surfaces could be described.
Here, the J 0 2 and J 0 3 invariants are generalized to consider anisotropy for orthogonal anisotropic sheet materials and are expressed as follows [14,19]: In the above equation, a i (i = 1-6) and b j (j = 1-11) are anisotropic coefficients.
In Equation ( 17), the effective stress (σ) is expressed as here, V is a constant for correcting the deviation of effective stress from uniaxial tension along the rolling direction and can be expressed as follows: From the associated flow rule, which regards the yield function f in Equation ( 17) as the plastic potential g, the plastic strain increment is the partial derivative of the plastic potential stress and is expressed as follows: The advantage of the proposed yield function compared with other types of anisotropic yield functions previously proposed using the principal deviation stress values (s 1 , s 2 , s 3 ) is that ∂f/∂σ ij for determining the strain increment can be calculated directly.
The normalized yield stress at σ xx = σ yy = σ b , under equi-biaxial tension in the plane stress state (σ 3 = 0) is expressed as follows: On the other hand, the Lankford coefficient value according to the material orientation, R θ , indicating the anisotropy characteristics of the anisotropic sheet material is expressed as follows: The Lankford modulus value in equi-biaxial tension, R b , is calculated in the following form: Additionally, the normalized uniaxial yield stress σ θ /σ 0 (see Equation ( 17)) according to the angle θ between the rolling direction and the loading direction is σ xx = σ θ • cos 2 θ, σ yy = σ θ • sin 2 θ.It can be expressed using the relationship θ and σ xy = σ θ • sin θ cos θ:

The Convexity Condition of the Yield Function
It is assumed that the yield surface is always convex because the internal plastic work must be dissipated during the plastic deformation of the material.This is called the convexity condition of the yield surface.In order to confirm the convexity of the yield condition for the proposed yield criterion in this study, the stress range in which the yield surface is convex was checked.The convexity of the proposed yield function is satisfied by confirming the convexities of f σ xx , σ yy , f σ xx , σ xy , and f σ yy , σ xy .Mathematically, to ensure the convexity of each f σ i , σ j function, it is necessary to check whether the corresponding Hessian H is positive semidefinite [20,43,44]: The positive semidefinite property of a Hessian matrix H can be proved by showing that the matrix eigenvalues are positive.Based on this analysis, it was confirmed that the yield function proposed in this study satisfies the convexity condition.

Effect of α and β on Yield Surface under Plane Stress
In order to confirm the role of the material parameters α and β of the proposed yield function in the isotropic form, a i = b j = 1, various combinations of these parameters were adopted and the shape of the yield surface in the space of principal stresses under plane stress condition was observed.Figure 1 shows the predicted yield locus of Equation ( 17) for various combinations of α and β parameters to confirm the effect of each parameter.If both parameters α and β are 0, the proposed yield function returns to the Mises yield function.
If α and β are both negative, the yield surface is inside the Mises yield surface, whereas if α and β are both positive, the yield surface is outside the Mises yield surface.Nevertheless, in both cases, the yield trajectory always passes through three points: (1, 0), (0, 1), and the equi-biaxial load (1, 1).When α is changed to α = −12 and α = 12, the shape of the yield surface does not change significantly (Figure 1a).In contrast, when β is changed to β = −3 and β = 3, the shape of the yield surface changes significantly (Figure 1b).That is, it can be seen that the parameter β affects the curvature of the yield surface more than α in the yield surface proposed in this study.

Effect of α and β on Yield Surface under Plane Stress
In order to confirm the role of the material parameters α and β of the proposed yield function in the isotropic form,  =  = 1, various combinations of these parameters were adopted and the shape of the yield surface in the space of principal stresses under plane stress condition was observed.Figure 1 shows the predicted yield locus of Equation ( 17) for various combinations of α and β parameters to confirm the effect of each parameter.If both parameters α and β are 0, the proposed yield function returns to the Mises yield function.If α and β are both negative, the yield surface is inside the Mises yield surface, whereas if α and β are both positive, the yield surface is outside the Mises yield surface.Nevertheless, in both cases, the yield trajectory always passes through three points: (1, 0), (0, 1), and the equi-biaxial load (1, 1).When α is changed to α = −12 and α = 12, the shape of the yield surface does not change significantly (Figure 1a).In contrast, when β is changed to β = −3 and β = 3, the shape of the yield surface changes significantly (Figure 1b).That is, it can be seen that the parameter β affects the curvature of the yield surface more than α in the yield surface proposed in this study.In addition, various combinations of the two parameters were reviewed to understand the flexibility of the proposed yield function. Figure 2 shows yield surfaces of different shapes according to parameters α and β.It can be seen that when α and β are negative numbers (α = −2 and β = −2), the curvature of the yield surface becomes smaller in the biaxial tensile region.In addition, various combinations of the two parameters were reviewed to understand the flexibility of the proposed yield function. Figure 2 shows yield surfaces of different shapes according to parameters α and β.It can be seen that when α and β are negative numbers (α = −2 and β = −2), the curvature of the yield surface becomes smaller in the biaxial tensile region.
For the case of Equation ( 16), the asymmetric type of yield function, Figure 3 shows the effect of material parameter β with α = 0.The shape of the yield surface captures well the SD effect similar to Cazacu and Barlat (2004) and Khan et al. (2012) in that the yield stress in compression is larger than that in tension.Additionally note the asymmetric yield shape in that the proposed yield function reproduces the higher strength in equi-biaxial tension than in equi-biaxial compression observed in a titanium sheet of HCP crystal structure.
Therefore, it can be accepted that the yield function proposed in this study can express various types of yield surfaces by properly combining the two parameters α and β.For the case of Equation ( 16), the asymmetric type of yield function, Figure 3 shows the effect of material parameter β with α = 0.The shape of the yield surface captures well the SD effect similar to Cazacu and Barlat (2004) and Khan et al. (2012) in that the yield stress in compression is larger than that in tension.Additionally note the asymmetric yield shape in that the proposed yield function reproduces the higher strength in equi-biaxial tension than in equi-biaxial compression observed in a titanium sheet of HCP crystal structure.Therefore, it can be accepted that the yield function proposed in this study can express various types of yield surfaces by properly combining the two parameters α and β.

AA6016-T4 Sheet
In conventional methods used to determine parameters of the yield function, the uniaxial tension, pure shear, and hydraulic bulge test and biaxial tensile experiment with cruciform have been widely used [45].In order to identify the parameters of the yield function developed by Cazacu-Barlat et al. (2001) for the AA6016-T4 sheet, the uniaxial tensile and hydraulic bulge test data were used [21,46,47].
In order to identify the parameters of the proposed yield function, the following error function was used to minimize the error value: (26) Here, n and m represent the number of yield stress σ θ and Lankford coefficient value R θ in the direction obtained in the uniaxial tensile test, respectively.σ b and R b are the yield stress and Lankford coefficient of equi-biaxial tension, respectively.t denotes the amount of biaxial tension and compression tested at various load ratios, and the superscripts denotes whether each value is an experimental value (data) or a calculated value (th).
Additionally, η i , γ i , δ, ν, and χ k are weight factors for each term on the right side contributing to the error function.These weights are values arbitrarily determined by the user according to which factors are given priority to identify a i and b j because the experimental results in this study are greater than the number of equations [48][49][50].Habraken et al. [49] set all of these weight factors to 1.0, and Nomura and Kuwabara [50] used the values that best fit the yield stress and anisotropy coefficient in the range of 1000~0.1.In this study, an appropriate value was used according to the problem so as not to fall into the local minimum value.(l) data k and (l ) th k are the distances between the origin and the experimental and predicted values in the principal stress space, respectively.That is, l = σ 2 1 + σ 2 2 .Using KV'12S and the Cazacu-Barlat yield function (2001), the yield surface, normalized yield stress, and Lankford anisotropy coefficients for the AA6016-T4 sheet were predicted.Figure 4a,c,d are shown together with the experimental data.The coefficients used in the KV'12S yield function for in-plane stress σ 3 = 0 were identified using 16 experimental data of uniaxial tensions, i.e., the yield stresses and the Lankford coefficients at θ = 0 • , 15 • , 30 • , 45 • , 60 • , 75 • , and 90 • , and the equiyield tension stress, i.e., equi-biaxial yield stress (σ b /σ 0 = 1.0 [21]) and the Lankford coefficient (R b = 1.05 [46]).These data were also used to minimize the error function, Equation (26).These coefficients are shown in Table 1.were predicted.Figure 4a,c,d are shown together with the experimental data.The coefficients used in the KV'12S yield function for in-plane stress  = 0 were identified using 16 experimental data of uniaxial tensions, i.e., the yield stresses and the Lankford coefficients at θ = 0°, 15°, 30°, 45°, 60°, 75°, and 90°, and the equiyield tension stress, i.e., equi-biaxial yield stress ( / = 1.0 [21]) and the Lankford coefficient ( = 1.05 [46]).These data were also used to minimize the error function, Equation (26).These coefficients are shown in Table 1.Similarly, in the Cazacu-Barlat yield function, uniaxial tensile yield stress and Lankford modulus values tested in the directions θ = 0°, 30°, 45°, 75°, and 90° and biaxial tensile yield stress  / = 1.0 were used to identify the coefficients of the Cazacu-Barlat yield function.Figures 4a,c,d show the predicted yield surface, normalized yield stress, and Lankford coefficient values in the KV'12S model using the Cazacu-Barlat yield Similarly, in the Cazacu-Barlat yield function, uniaxial tensile yield stress and Lankford modulus values tested in the directions θ = 0 • , 30 • , 45 • , 75 • , and 90 • and biaxial tensile yield stress σ b /σ 0 = 1.0 were used to identify the coefficients of the Cazacu-Barlat yield function. Figure 4a,c,d show the predicted yield surface, normalized yield stress, and Lankford coefficient values in the KV'12S model using the Cazacu-Barlat yield function, respectively.Figure 4b shows the shape of the yield surface considering the effect of shear stress on the KV'12S yield function.
From the results in Figure 4, the value of the Lankford coefficient according to the yield surface, normalized yield stress, and angle predicted by the KV'12S model had a low error value (i.e., F = 0.0006 for KV'12S and F = 0.0007 for the Cazacu-Barlat model), and it can be seen that the experimental results are well represented.However, the predicted uniaxial yield stress σ θ /σ 0 according to the angle was partially inconsistent with the corresponding experimental data at θ = 0 • , 45 • , and 90 • .This is thought to be because overfitting occurred when 16 equations were used to fit the experimental data.
In order to prevent such overfitting at θ = 0 • , 45 • , and 90 • , the anisotropy coefficients in the Kim-Van yield function were identified using a total of 12 experimental data, i.e., uniaxial tensile yield stress and Lankford values at θ = 0 • , 15 • , 45 • , 75 • , and 90 • , and equi-biaxial tension, which were σ b /σ 0 = 1.0 and R b =1.05.The anisotropic coefficients identified are summarized in Table 2 and the predictions are depicted in Figure 5. overfitting occurred when 16 equations were used to fit the experimental data.
In order to prevent such overfitting at θ = 0°, 45°, and 90°, the anisotropy coefficients in the Kim-Van yield function were identified using a total of 12 experimental data, i.e., uniaxial tensile yield stress and Lankford values at θ = 0°, 15°, 45°, 75°, and 90°, and equibiaxial tension, which were  / = 1.0 and  =1.05.The anisotropic coefficients identified are summarized in Table 2 and the predictions are depicted in Figure 5.In Figure 5, it can be seen that the predicted yield locus, surface, normalized yield stress, and Lankford coefficient values from the KV'12S yield function agreed well with the experiment at the angles taken to identify the parameters.However, the error value in this case (F = 0.0013) was larger than that of the Cazacu-Barlat model (F = 0.0007).The reason is incorrect predictions in the other data θ = 30 • and 60 • .Therefore, it can be concluded that it is necessary to take appropriate experimental data for each material to improve the prediction of material behavior.

AA2090-Ts3 Sheet
In previous studies, in the case of the AA2090-T3 [51,52] sheet, the Yld2004-18p function was used to successfully describe the plastic anisotropy properties of strongly anisotropic materials.However, the Yld2004-18p yield function requires many experiments to identify parameters.
Here, to confirm the applicability of the KV'12S yield function to strongly anisotropic material, the predicted yield surface, normalized uniaxial yield stress, and Lankford coefficient for the AA2090-T3 sheet were compared with experimental results.
To identify the parameters of the KV'12S yield function, uniaxial tensile yield stress and Lankford modulus values according to θ = 0 • , 15 • , 45 • , and 90 • ; biaxial tensile yield stress σ b /σ 0 = 1.035; and biaxial yield stress values at three different load ratios of biaxial tension (i.e., σ x : σ y = 1 : 2, 4 : 3, 2 : 1) were used.On the other hand, the coefficient of the Cazacu 2018 yield function (refer to Equation ( 7)) was identified by minimizing the error function (22) using the uniaxial tensile yield stress and Lankford coefficient value for θ = 0 • , 30 • , 45 • , 75 • , and 90 • and biaxial yield stress σ b /σ 0 = 1.035.Table 3 shows the coefficients of the two yield functions. Figure 6a,c,d show the predicted yield surface, normalized yield stress, and Lankford coefficient values from the KV'12S and Cazacu 2018 yield functions together with experimental data.Figure 6a represents the yield locus with different shear stress σ xy /σ 0 and shows clearly that the yield surfaces predicted by Equation ( 12) did not exhibit the same shape for different shear stresses; additionally, it shows couplings between shear and normal components of stress, as pointed in many crystal plasticity studies [53].

AA2090-Ts3 Sheet
In previous studies, in the case of the AA2090-T3 [51,52] sheet, the Yld2004-18p function was used to successfully describe the plastic anisotropy properties of strongly anisotropic materials.However, the Yld2004-18p yield function requires many experiments to identify parameters.
Here, to confirm the applicability of the KV'12S yield function to strongly anisotropic material, the predicted yield surface, normalized uniaxial yield stress, and Lankford coefficient for the AA2090-T3 sheet were compared with experimental results.
To identify the parameters of the KV'12S yield function, uniaxial tensile yield stress and Lankford modulus values according to θ = 0°, 15°, 45°, and 90°; biaxial tensile yield stress  / = 1.035; and biaxial yield stress values at three different load ratios of biaxial tension (i.e.,  :  = 1: 2, 4: 3, 2: 1) were used.On the other hand, the coefficient of the Cazacu 2018 yield function (refer to Equation ( 7)) was identified by minimizing the error function (22) using the uniaxial tensile yield stress and Lankford coefficient value for θ = 0°, 30°, 45°, 75°, and 90° and biaxial yield stress  / = 1.035 .Table 3 shows the coefficients of the two yield functions.6a represents the yield locus with different shear stress  / and shows clearly that the yield surfaces predicted by Equation ( 12) did not exhibit the same shape for different shear stresses; additionally, it shows couplings between shear and normal components of stress, as pointed in many crystal plasticity studies.[53]   Moreover, the shape of the yield surface was quite sharp in equi-biaxial tension, and strong anisotropic behavior is shown in the Lankford coefficient value (Figure 6d).The prediction of the KV'12S yield function agreed well with the experiment.Error values were calculated to evaluate the predictability of material behavior for two yield functions.The error value of 0.0837 was calculated from the KV'12S yield function, while the Cazacu 2018 yield function could not be calculated.The reason is that the convexity of this yield function was not satisfied in the case of the Cazacu 2018 yield function as shown in the yield surface (Figure 6a).
To compare the difference in convexity between the two yield functions, a Hessian matrix (H) was observed.The first eigenvalue of a matrix H is always greater than or equal to the second eigenvalue.Therefore, only the second eigenvalue ( ) was considered in this study.
In the case of the KV'12S yield function, it can be seen that the second eigenvalue of matrix H was always greater than or equal to 0 in the range of  = [−200, 200] MPa and  = [−200, 200] MPa (Figure 6a).On the other hand, the Cazacu 2018 yield function did not satisfy the convexity under some conditions (Figure 7b).In particular, in the range  Moreover, the shape of the yield surface was quite sharp in equi-biaxial tension, and strong anisotropic behavior is shown in the Lankford coefficient value (Figure 6d).The prediction of the KV'12S yield function agreed well with the experiment.Error values were calculated to evaluate the predictability of material behavior for two yield functions.The error value of 0.0837 was calculated from the KV'12S yield function, while the Cazacu 2018 yield function could not be calculated.The reason is that the convexity of this yield function was not satisfied in the case of the Cazacu 2018 yield function as shown in the yield surface (Figure 6a).
To compare the difference in convexity between the two yield functions, a Hessian matrix (H) was observed.The first eigenvalue of a matrix H is always greater than or equal to the second eigenvalue.Therefore, only the second eigenvalue (λ 2 ) was considered in this study.
Figure 7c,d show the shape of the second invariant in the principal stress plane.The KV'12S yield function is a smooth shape without vertices as shown in Figure 7c.However, the Cazacu 2018 yield function as shown in Figure 7d has many protrusions.This phenomenon indicates that the yield curve predicted using the Cazacu 2018 model did not satisfy the convexity condition for the strongly anisotropic material.From this result, it can be concluded that the combination of the two parameters α and β in the KV'12S yield function proposed in this study can control the shape of the yield surface that satisfies the convexity of the strongly anisotropic material.In other words, the KV'12S yield function can control the various shapes of the yield surface in a wider range.
second eigenvalue.Therefore, only the second eigenvalue ( ) was considered in this study.
In the case of the KV'12S yield function, it can be seen that the second eigenvalue of matrix H was always greater than or equal to 0 in the range of  = [−200, 200] MPa and  = [−200, 200] MPa (Figure 6a).On the other hand, the Cazacu 2018 yield function did not satisfy the convexity under some conditions (Figure 7b).In particular, in the range   7d has many protrusions.This phenomenon indicates that the yield curve predicted using the Cazacu 2018 model did not satisfy the convexity condition for the strongly anisotropic material.From this result, it can be concluded that the combination of the two parameters α and β in the KV'12S yield function proposed in this study can control the shape of the yield surface that satisfies the convexity of the strongly anisotropic material.In other words, the KV'12S yield function can control the various shapes of the yield surface in a wider range.

Pure Titanium Sheet
A commercially available pure titanium (CP-Ti) sheet having an HCP crystal structure exhibited an asymmetric shape in the yield curve.To explain the asymmetric plastic yield behavior of CP-Ti materials, Yoon et al. [42] proposed a new yield function having a coupling effect of all three stress invariants as shown in Equation (10).To describe the SD effect of the yield stress difference appearing in tension and compression tests, we adopted the KV'12A yield function, and the anisotropic parameters were determined from the experimental results of uniaxial tension and compression tests and equi-biaxial tension and compression tests.
Figure 8 shows experimental data from [42] and the Yoon 2014 yield function together to confirm the applicability of the KV'12A asymmetric yield function of Equation ( 12) for HCP material.

Pure Titanium Sheet
A commercially available pure titanium (CP-Ti) sheet having an HCP crystal structure exhibited an asymmetric shape in the yield curve.To explain the asymmetric plastic yield behavior of CP-Ti materials, Yoon et al. [42] proposed a new yield function having a coupling effect of all three stress invariants as shown in Equation (10).To describe the SD effect of the yield stress difference appearing in tension and compression tests, we adopted the KV'12A yield function, and the anisotropic parameters were determined from the experimental results of uniaxial tension and compression tests and equi-biaxial tension and compression tests.
Figure 8 shows experimental data from [42] and the Yoon 2014 yield function together to confirm the applicability of the KV'12A asymmetric yield function of Equation ( 12) for HCP material.
(2) KV'12S can also well represent various types of yield curves including strongly anisotropic materials.It was shown that the proposed yield function can implement various types of yield surface shapes by combining the two parameters α and β. (3) To verify the flexibility of the proposed model, the convexity of the yield surface was checked and compared with the Cazacu 2018 yield function.The convexity of KV'12S was satisfied in most stress ranges, but the Cazacu 2018 yield function for the material covered in this study did not satisfy the convexity in some stress ranges.(4) As a result of applying the KV'12S yield function proposed in this study to the AA2090-T4 material, which exhibited very strong anisotropy, the results deviated from some experimental data.Attention should be paid to the selection of experimental data used for the parameter optimization of the yield function proposed in this study.In this case, twelve experimental data (e.g., ten experimental data from the uniaxial tensile test and two experimental data from the equi-biaxial tension test) are recommended.(5) The KV'21A yield function was used to predict the yield surface of a commercial pure titanium sheet that exhibited an asymmetrical yield surface shape, and it was compared with the Yoon 2014 yield function and experimental data.Both yield functions were in good agreement with the experiment overall.
In a future study, the anisotropic yield function proposed in this study will be applied to the forming limit curve prediction.In addition, the effectiveness of the proposed yield function is evaluated by applying it to the finite element analysis of the sheet material processing process such as the cup drawing process and the incremental sheet forming process.

Figure 1 .
Figure 1.Yield locus according to proposed isotropic yield function for different values of parameters α and β: (a) effect of parameter α, and (b) effect of parameter β.

Figure 1 .
Figure 1.Yield locus according to proposed isotropic yield function for different values of parameters α and β: (a) effect of parameter α, and (b) effect of parameter β.

Figure 6 .
Figure 6.Predicted anisotropy according to the Kim-Van and the Cazacu 2018 yield functions for AA2090-T3 sheet with 12 equations for fitting (a) yield locus, (b) 3D yield surface, (c) normalized uniaxial yield stresses, and (d) R-values.

Figure 6 .
Figure 6.Predicted anisotropy according to the Kim-Van and the Cazacu 2018 yield functions for AA2090-T3 sheet with 12 equations for fitting (a) yield locus, (b) 3D yield surface, (c) normalized uniaxial yield stresses, and (d) R-values.

Figure
Figure7c,d show the shape of the second invariant in the principal stress plane.The KV'12S yield function is a smooth shape without vertices as shown in Figure7c.However, the Cazacu 2018 yield function as shown in Figure7dhas many protrusions.This phenomenon indicates that the yield curve predicted using the Cazacu 2018 model did not satisfy the convexity condition for the strongly anisotropic material.From this result, it can be concluded that the combination of the two parameters α and β in the KV'12S yield function proposed in this study can control the shape of the yield surface that satisfies the convexity of the strongly anisotropic material.In other words, the KV'12S yield function can control the various shapes of the yield surface in a wider range.
[26] form similar to the Cazacu-Barlat yield function (2001), Lou and Yoon (2017)[26]modified the isotropic Drucker yield function by expressing s = L σ and the second and third invariants as linearly transformed stress tensors:

Table 3 .
Anisotropic coefficients of the Kim-Van and the Cazacu 2018 yield functions for AA2090-T3 sheet.