On mixed isotropic-distortional hardening

On mixed isotropic-distortional hardening

International Journal of Mechanical Sciences 92 (2015) 259–268 Contents lists available at ScienceDirect International Journal of Mechanical Science...

2MB Sizes 6 Downloads 34 Views

International Journal of Mechanical Sciences 92 (2015) 259–268

Contents lists available at ScienceDirect

International Journal of Mechanical Sciences journal homepage: www.elsevier.com/locate/ijmecsci

On mixed isotropic-distortional hardening Rikard Rentmeester a,n, Larsgunnar Nilsson b a b

SAAB AB Aeronautics, Bröderna Ugglas gata, SE-58254 Linköping, Sweden Division of Solid Mechanics, Linköping University, SE-58183 Linköping, Sweden

art ic l e i nf o

a b s t r a c t

Article history: Received 3 June 2014 Accepted 18 September 2014 Available online 28 September 2014

Mixed isotropic-distortional hardening allows for individual stress–plastic strain relations in different straining directions. Such hardening can be obtained by allowing the parameters in the effective stress function depend on anisotropy functions of the equivalent plastic strain. A methodology to calibrate these anisotropy functions is proposed in this work, and is demonstrated on an austenitic strainless steel. A high exponent eight parameter effective stress function for plane stress states is utilised. The anisotropy functions are calibrated by the use of experimental data from uniaxial tensile test data in three material directions and a balanced biaxial test. The plastic anisotropy is evaluated at a finite number of plastic strains, and it is assumed to vary piecewise linearly with respect to the equivalent plastic strain. At each level of plastic strain, the anisotropy is correctly represented, even if rather large increments in plastic strain are used in the calibration. It was found that there are at least two sets of anisotropy functions which satisfy the conditions in the calibration procedure. The resulting uniaxial stress–strain relations from the two sets of anisotropy functions in four additional straining directions, not included in the calibration set, were compared to the corresponding experimental data. From this validation, one of the anisotropy function sets could be discarded, whereas the other one gave a good prediction of the stress–strain relations in all the four additional directions. & 2015 Elsevier Ltd. All rights reserved.

Keywords: Plastic anisotropy Distortional hardening

1. Introduction The demand for lighter and stronger products has increased the use of high strength steels. In general, high strength steels, such as martensitic steels, have a high initial yield strength but suffer from a less ductility, which results in less formability. Thus, the choice of material in a sheet metal application is often a compromise between strength and formability. In order to fully utilise the formability of the material, there is a need for accurate formability predictions. One major material failure mode which limits the formability is plastic instability, as this causes strain localisation and eventually a subsequent material fracture. In order to accurately predict this instability, the plastic hardening and anisotropy of the material are important. Plastic anisotropy in metals has traditionally been described by anisotropic effective stress functions, giving anisotropic yield surfaces. Recent development has resulted in sophisticated multiparameter functions [1,2], which are able to account both for anisotropy in the yield stresses and in the plastic flow. The response under proportional loading is generally governed by the shape of


Corresponding author. Tel.: þ 46 13 186639. E-mail address: [email protected] (R. Rentmeester).

http://dx.doi.org/10.1016/j.ijmecsci.2014.09.013 0020-7403/& 2015 Elsevier Ltd. All rights reserved.

the yield surface and the plastic hardening function. However, it has been shown that both a translation and a distortion of the yield surface take place during plastic deformation, see [3]. Thus, kinematic and distortional hardening can be mixed in order to describe the material response under non-proportional loadings [3,4]. A similar effect can be achieved with a mixed isotropic-kinematic hardening since kinematic hardening is a special case of anisotropic hardening [5]. The yield stress changes differently at other stress states than the state defined by the current loading. The introduction of a back stress tensor, such as is used in the kinematic hardening approach, removes the homogeneity of the effective stress function. A new approach was proposed in [6] for describing the Bauschinger effect, while keeping the homogeneity of the effective stress function by introducing an additional fluctuating tensor. By their approach, the resulting yield surface contains the origin and is thus unable to account for reversed plastic yielding while unloading along the loading path. The calibration procedure of the anisotropy parameters is often based on a least square fit to a set of target experimental values. The validity of the resulting yield condition for stress states not included in the fit has usually been investigated by comparing predicted results to corresponding experimental data. For example, the variations of the uniaxial yield stress and R-value with angle to the rolling direction, RD, have been investigated in [7].


R. Rentmeester, L. Nilsson / International Journal of Mechanical Sciences 92 (2015) 259–268

The yield surface is usually calibrated to experimental data at the onset of plastic deformation. If isotropic or mixed isotropickinematic hardening is used, the plastic hardening under proportional loading becomes similar regardless of stress state. Even though the initial plastic anisotropy may be well described by the yield surface, this approach is insufficient to describe different hardening rates in different directions. However, such anisotropic hardening has been considered as important for formability predictions, see [8,9]. Plunkett et al. [10] described the varying anisotropy by letting the anisotropy parameters of the effective stress function depend on the equivalent plastic strain. Such varying anisotropy parameters governing the plastic anisotropy are henceforth denoted as anisotropy functions. Using this approach, not only the initial plastic anisotropy, but also the subsequent anisotropic plastic hardening, is represented. Aretz [9] used a similar approach with the eight parameter effective stress function YLD2003 [11], in which the anisotropy functions varied with plastic work, and showed that mixed isotropic-distortional hardening is important for predicting the anisotropy in forming processes. The hardening approach presented in this paper gives a similar response as ordinary isotropic hardening in the case of non-proportional hardening, however, a combination with kinematic hardening is straightforward. Mixed isotropic-distortional hardening is further examined herein. The eight parameter effective stress function YLD2003 has been used. YLD2003 is defined by eight anisotropy parameters and requires at least eight experimental data points for calibration. The anisotropy functions were calibrated to the evolution of four yield stress ratios and four plastic strain ratios. In this work, the anisotropy functions were assumed to be piecewise linear, i.e. the anisotropy functions were calibrated at a finite number of equivalent plastic strain levels. If the functions are calibrated at a small number of plastic strains, this will affect the result. It was shown that the anisotropy functions depend not only on the anisotropy but also on the increment in equivalent plastic strain between the vertices. The anisotropy functions are calibrated by a novel calibration method, which takes the piecewise linear variation of the anisotropy functions into account. With this method, only a small number of calibrations have to be conducted in order to define the mixed isotropic-distortional hardening. The method is demonstrated on an austenitic stainless steel, HyTens 1000, within the EN 1.4310 standard. Moreover, it was found that for a set of eight experimental target values describing the plastic anisotropy, two solutions to the anisotropy functions exist, both of which describe the target values exactly. Therefore, to find the unique solution, uniaxial tensile tests in four additional material directions and shear tests in two material directions were carried out. In Section 2, the experimental procedures and data used in the calibration are presented. In the subsequent sections, the constitutive model along with the calibration procedure is described. It was found that two alternative sets of anisotropy functions fulfill the target values in the calibration procedure of the anisotropy functions. In Section 6, the predicted uniaxial stress–strain relations in four additional directions, which not are included in the calibration procedure, are compared to the corresponding experiments in order to test the validity of the two solutions. Furthermore, the force–displacement relations in two shear tests obtained from the two solutions were compared to the corresponding experimental results in order to further distinguish between the two solutions.

2. Experimental work 2.1. Material The material under investigation is an austenitic stainless steel of the brand HyTensX, see [13], within the EN 1.4310/AISI 301 standard. The material is cold rolled, and the resulting steel grade is denoted as HyTens 1000 or 1.4310 C1000. Its chemical composition is presented in Table 1. The thickness of the steel is 1.5 mm. HyTensX is meta-stable, i.e. phase transformation takes place during deformation. The authors have previously investigated strain ageing in this material, see [14]. It was found that the material is sensitive to both dynamic and static strain ageing. The phase transformation model proposed by [15], combined with a linear mixture rule of the yield stresses in the two phases, was used in order to account for martensitic transformation. Recently, the martensitic transformation in the same material (HyTensX) was investigated in [16] using an enhanced macroscopic model, in which the effect of the triaxiality on phase transformation was included. 2.2. Tensile tests Uniaxial tensile tests were conducted in the directions ϕ ¼0, 151, 301, 451, 601, 751 and 901 with respect to the RD . The geometry of the test specimen is shown in Fig. 1. A mechanical extenso-meter with initial length l0 ¼ 25 mm was attached to the specimen in order to measure the longitudinal strain, εl ¼ lnðl=l0 Þ. A second extenso-meter was attached in order to measure the width change, w  w0 , in order to evaluate the transversal strain εw ¼ lnðw=w0 Þ. The deformation was controlled by applying a constant crosshead velocity v¼3 mm/min, resulting in an approximate overall strain rate of ε_ l  10  3 s  1 . Furthermore, the force, F, was continuously measured in order to evaluate the stress, σ ¼ F expðεl Þ=A0 , where A0 is the undeformed cross section area. It can be noticed that the local strain rates were significantly higher due to serrated yielding, something which is not considered in the present work. Three repetitions were conducted in each direction. The specimen to specimen variation of the stress–strain relations was found to be small. In fact, it may be noticed that the stairs in the stress–

Fig. 1. Geometry of the tensile specimen.

Table 1 Chemical composition of the investigated steel grade, HyTensX (%). C
















Fig. 2. True stress–longitudinal strain relations from tensile tests. Three tests were conducted in each material direction.

R. Rentmeester, L. Nilsson / International Journal of Mechanical Sciences 92 (2015) 259–268

strain relations, due to serrated flow, obtained from different specimens are close to each other in the same ϕ-direction. This indicates that the hardening rates in different samples are very close, since the serrated yielding is an instability phenomenon which is related to the plastic hardening rate. The stress–strain relations for ϕ ¼0, 451 and 901 are presented in Fig. 2. Experimental data from these specific directions are included in the subsequent calibration, and the stress–strain relations are valid up to the onset of diffuse necking. The εw to εl relations are presented in Fig. 3. These relations are truncated prior to diffuse necking due to insufficient range of the extensometer used for the width measurement. The corresponding results for all other directions, i.e. ϕ ¼151, 301, 601 and 751, are used for validation, i.e. these results are not included in the calibration and are therefore presented in a subsequent section together with the corresponding numerical predictions.

2.3. Bulge test A balanced biaxial bulge test in which the pressure was applied by a punch made from silicon was performed. The in-plane strains and the radius of the bulge were measured using Digital Image Correlation. From the in-plane strains, the hydrostatic pressure in the silicon and the bulge radius the biaxial stress–strain relation could be found, as presented in Fig. 4(a). It should be noted that the stress levels are not to be trusted at small strains, due to the error in measuring the sheet radius at small curvatures and the bending of the sheet at initial loading. The latter effect results in inhomogeneous strain and stress distributions through the thickness. The bending effect decreases with the punch motion, resulting in a more uniform stress and strain distribution through the thickness. The friction between the silicon bulge and the sheet is ignored, since it is assumed to have a minor influence on the result. Furthermore, a relationship between the plastic strain in the RD, εRD , and the plastic strain in the transversal


direction, TD, εTD , was obtained, as shown in Fig. 4(b). Aretz and Keller [12] demonstrated the difficulties concerning the evaluation of the relation between the strain components, especially for an anisotropic sheet. The balanced biaxial stress state requires very low friction and axisymmetry, i.e. the same radii at bend in all directions. It was concluded that this is not the case for anisotropic materials. This fact affects the stress–strain relation as well as the relation between the strain components. However, in the present study, the material displays low plastic anisotropy, and thus, the method is considered as sufficiently good. 2.4. Shear tests Shear tests were conducted in the ϕ-directions 451 and 901. The geometry of the shear test specimen is shown in Fig. 5. The test specimen was mounted with clevis pins through holes in the specimen, which were lubricated to minimise friction and to permit possible rotations. The deformation was applied by a constant crosshead velocity v¼ 0.03 mm/min, which resulted in a strain rate of ε_ p  10  4 s  1 . The displacement was measured on the clamps in which the clevis pins were mounted. Three repetitions in each direction were performed.

3. Constitutive model Since the elastic part of the deformation is considered to be small compared to the plastic part, a hypo-elastic relation is assumed, and an additive decomposition of the rate of deformation tensor, D, was introduced: D ¼ De þ Dp e

ð1Þ p

where D and D are the elastic and plastic parts, respectively. The plastic rate of deformation is given by an associated flow rule: ∂f Dp ¼ λ_ ð2Þ ∂σ where λ_ is the plastic multiplier and f is the flow function, i.e. f ¼ σ ; ε p  σ y ðε p Þ


and where the reference yield stress, σ y , is a function of the equivalent plastic strain ε p . The yield condition f¼ 0 defines a yield surface in the stress space at which σ ; ε p ¼ σ y ðε p Þ, and thus it constitutes a surface with constant equivalent plastic strain. 3.1. Plastic anisotropy The effective stress function YLD2003 [11] is defined as Fig. 3. Relations between longitudinal and transversal strains from tensile tests. Three tests were conducted in each material direction.

2σ ðσ ; ε p Þ ¼ jσ 01 ja þjσ 02 ja þ jσ ″1  σ ″2 ja a

Fig. 4. Results from the balanced biaxial test: (a) Stress–plastic strain relationship. (b) Relation between εpRD and εpTD.



R. Rentmeester, L. Nilsson / International Journal of Mechanical Sciences 92 (2015) 259–268

Fig. 5. (a) Geometry of the shear test specimen. Dimensions in mm. (b) Finite element model of the shear test.


σ 01 σ 02


σ ″1 σ ″2

sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  2 A8 ðε p Þσ 11 þ A1 ðε p Þσ 22 A2 ðε p Þσ 11  A3 ðε p Þσ 22 7 ¼ þ A4 ðε p Þ2 σ 12 σ 21 2 2


sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  2 σ 11 þ σ 22 A5 ðε p Þσ 11  A6 ðε p Þσ 22 7 ¼ þ A7 ðε p Þ2 σ 12 σ 21 2 2


In Aretz's original work, Ai were constant parameters governing the plastic anisotropy, whereas here they are functions of equivalent plastic strain, thus Ai ¼ Ai ðε p Þ. The effective stress function σ is a first order homogeneous function of the stress tensor, i.e.


∂σ ¼σ ∂σ


From the definition of ε p and by using Eqs. (2) and (6), the rate of _ p , is obtained as specific plastic work, w   ∂σ _ p ¼ σ ε_ p ¼ σ : Dp ¼ σ : λ_ ð7Þ ¼ σ λ_ w ∂σ giving the relation ε_ p ¼ λ_ . 4. Model calibration This section gives a description of the procedures used to calibrate the yield criterion and its evolution during plastic deformation. To begin with, the analytical function in Eq. (8) was fitted to the stress–strain relation from the tensile tests such that the influence from the serrated plastic flow was avoided. In Section 4.2, the ratios between the uniaxial stress–strain relations, i.e. the yield stress ratios, and the plastic strain ratios are evaluated. The anisotropy functions are presented in the last two subsections. 4.1. Uniaxial hardening Analytical functions were fitted to the experimental stress–strain data. In a general case this is not necessary, although it is recommended in cases where the stress–strain relations are serrated. Such serration should be excluded from the stress–strain relation in the model, since it is a structural phenomenon associated with the

Table 2 Plastic hardening parameters. ϕ ð1Þ




σ 0;ϕ (MPa) Q 1;ϕ (MPa) C 1;ϕ (–) Q 2;ϕ (MPa) C 2;ϕ (–)

753 412 3.19 150 624 0.0135

775 344 9.65 142 528 0.0350

758 536 7.54 149 436 0.0453

909 8204 1.53 0.0889

941 3113 1.11 0.115

903 1985 0.814 0.0798

593 2196 0.597

634 1852 0.532

905 2003 0.821

εðt1Þ (–) ϕ G1;ϕ (MPa) K 1;ϕ (MPa) n1 (–) (–) εðt2Þ ϕ G2;ϕ (MPa) K 2;ϕ (MPa) n2 (–)

uniaxial tensile stress state rather than a physical material property. Any type of function able to represent the characteristics of the plastic hardening may be adopted. In this work, an extended Voce hardening function was combined with two powerlaw functions, i.e.

σ y;ϕ ¼

8 2  C i;ϕ εpϕ > > Þ; > > σ 0;ϕ þ ∑ Q i;ϕ ð1  e > < i¼1 G1;ϕ þ K 1;ϕ ðεpϕ Þn1 ; > > > > > : G2;ϕ þ K 2;ϕ ðεp Þn2 ; ϕ

εpϕ r εðt1Þ ϕ p ðt2Þ εðt1Þ ϕ r εϕ r εϕ


p εðt2Þ ϕ r εϕ

where σ 0;ϕ , Q i;ϕ , C i;ϕ , Gi;ϕ , K i;ϕ , ni and ε p r εðtiÞ ϕ , i¼1,2, are material parameters which are individual for each ϕ-direction. C1 continuity was required at both the transition strains. To afford such a continuity, two constraints were applied at both transition strains in the parameter evaluation procedure. As a consequence, the number of independent parameters decreases by four in each direction. The parameter values were found from a least square fit procedure, and are presented in Table 2. The corresponding stress–plastic strain relations are shown in Fig. 6.

R. Rentmeester, L. Nilsson / International Journal of Mechanical Sciences 92 (2015) 259–268


of constant plastic work, as well as a surface of constant plastic strain. σ ref ðε p Þ may be chosen arbitrarily. Here, it is set to be σ ref ¼ σ 90 , and thus, y90 ¼ 1 and εp90 ¼ ε p in this uniaxial stress state. In other directions, the yield stress ratios are evaluated at a finite number of equivalent plastic strain increments, ε pi ; i ¼ 1; …nε , and are assumed to vary linearly within each increment, i.e. yϕ ðε p Þ ¼ yϕ;i  1 þ

ε p  ε pi 1 ðy  y Þ; ε pi  ε pi 1 ϕ;i ϕ;i  1

ε pi 1 r ε p r ε pi


where the initial values are evaluated from yϕ jε p ¼ 0 ¼ Fig. 6. Stress–strain hardening relations from the tensile tests and the corresponding fitted functions.

dε p

Conventionally, anisotropic yield criteria have been calibrated with uniaxial tensile tests in various ϕ-directions. From these tests, the yield stress ratio, y, and plastic strain ratio, k, are defined as  dεpwϕ  σ ϕ ðε p Þ p p ð9Þ yϕ ðε Þ ¼ ; kϕ ðε Þ ¼  ; σ ref ðε p Þ dεp  p lϕ ε

where σ ref ðε Þ denotes the reference yield stress as a function of the equivalent plastic strain, ε p . εlϕ denotes the longitudinal logarithmic strain in the ϕ-direction, and εwϕ denotes the corresponding transversal strain. The superscripted p denotes the plastic part of the strain. If plastic incompressibility is assumed, kϕ is related to the Lankford parameter Rϕ by


If Eq. (13) is inserted into Eq. (12a), one has dεplϕ

4.2. Plastic anisotropy

σ 0;ϕ σ 0;ref


1 ¼ yϕ y


ε p  ε pi 1 ϕ;i  1 þ ε p  ε p ðyϕ;i  yϕ;i  1 Þ i i1

Rϕ ¼




p lϕ


1 þ dεp


kϕ 1 þ kϕ


where the subscript ND denotes the normal direction. Typically, tensile tests in the ϕ ¼0, 451, and 901 directions are used in the calibration procedure. Furthermore, a balanced biaxial stress state, σ TD ¼ σ RD , has often been used in order to further evaluate the anisotropy, e.g. see [7]. Accordingly, the yield stress ratio, yb, and plastic strain ratio, kb, of the biaxial stress state are defined as  dεpTD  σ b ðε p Þ p p ð11Þ ; k b ðε Þ ¼ p  yb ðε Þ ¼ σ ref ðε p Þ dεRD  p ε

4.3. Evaluation of the yield stress ratios Relations between the uniaxial stress and the equivalent plastic strain, i.e. the σ ϕ to ε p relations, are found from the definition of equivalent plastic strain, i.e. σ : Dp ¼ σ ε_ p . In the special cases of uniaxial tension in the ϕ-direction and a balanced biaxial stress state, the definition yields, see [5],

σ ϕ ε_ plϕ ¼ σ ε_ p ¼

ε pi 1 r ε p r ε pi

εplϕ ðε p Þ ¼ εplϕ ðε pi 1 Þ þ

ε pi  ε pi 1

yϕ;i  yϕ;i  1


ðyϕ;i  yϕ;i  1 Þε p þ yϕ;i  1 ε pi  yϕ;i ε pi 1 yϕ;i  1 ðε

p i 


p i  1Þ

ε pi 1 r ε p r ε pi In particular


p lϕ;i



p ð pÞ ¼ lϕ;i i

εϕ;i  1 þ p

ε pi  ε pi 1

yϕ;i  yϕ;i  1



yϕ;i  1

yb dt

yϕ;i σ ref ðε pi Þ ¼ σ ϕ ðεpϕ;i ðε pi ÞÞ


Eq. (18) is solved for yϕ;i with εpϕ;i according to Eq. (17) for each increment of equivalent plastic strain. The resulting yield stress ratios, with an equivalent plastic strain increment 0.01, are shown in Fig. 7. Other increment lengths would result in slightly different yϕ functions. However, with the proposed method, the stress representation in the ϕ-direction would still have been correct at the vertices of the yϕ -function despite a slight discrepancy in between, see Section 4.3.1. A similar procedure may be used for the biaxial (or any other) stress state. However, due to the uncertainty of the biaxial yield stress at small plastic strains and motivated by the similarity in test results between σ b ðjεpND jÞ and σ ref ðε p Þ, the biaxial yield stress ratio was set constant, i.e. yb ¼ 1, and thus dε p ¼ jdεpND j. 4.3.1. Sensitivity to increment size Different values of the equivalent plastic strain increment, Δε py ¼ ε pi  ε pi 1 , used to evaluate yϕ in Eq. (18), will result in slightly different relations between the uniaxial yield stress ratios and the equivalent plastic strain. The linear variation of the yield stress ratios will affect the stress–strain relations. This possible linear variation is taken into account in the methodology



ð12Þ where plastic incompressibility has been assumed, i.e. dε dεpRD þ dεpND ¼ 0. The relation between the plastic work and the equivalent plastic strain implies that the yield surface constitutes a surface


From Eq. (9), the following equation is found




σ ϕ _ p dε p ε ) ¼ dεplϕ ðaÞ 



σ dε dε ) σ b ε_ pRD þ ε_ pTD ¼ σ ε_ p ¼ b ¼ dεpTD þ dεpRD ¼  dεpND ðbÞ 


The result after integration is




p TD þ

Fig. 7. Yield stress ratios as a function of the equivalent plastic strain.


R. Rentmeester, L. Nilsson / International Journal of Mechanical Sciences 92 (2015) 259–268

Fig. 8. The uniaxial yield stress ratio in the RD for different values of Δε py . Fig. 11. Anisotropy functions according to Solution 1.

Fig. 9. The uniaxial yield stress in the RD for different values of Δε py .

Fig. 12. The evolution of the yield locus according to Solution 1.

Fig. 10. Biaxial plastic strain ratio.

presented here. The obtained yield stress ratio, y00, using different values of Δε py is shown in Fig. 8. The corresponding stress–plastic strain relations are shown in Fig. 9. It is noted that a large plastic strain increment, e.g. Δε py ¼ 0:2, is not able to predict the target stress–strain relation between the evaluated plastic strain levels. However, the predicted stress–plastic strain relation coincides with the target relation at the equivalent plastic strain levels at which the yield stress ratio has been evaluated, i.e. ε p ¼ 0, 0.2, 0.4, 0.6, 0.8 and 1.0. Thus, the choice of Δε py is a matter of plastic strain levels at which the predicted stress will coincide with the target stress. Furthermore, the smaller the equivalent plastic strain increments that are used in the evaluation of the yield stress ratios, the closer the predicted stress–plastic strain relation becomes to the target curve. If instead the yield stress ratio is evaluated at a finite number of levels of plastic work, this linear variation is not taken into account. If large increments are used, the accuracy of the predicted stress–strain relation will be lost.

Fig. 13. Two alternative yield surfaces fulfilling the same set of target values y00 ; y45 ; y90 ; yb ; k00 ; k45 ; k90 and kb.

evaluate kϕ , see [5]. Following reference [14], the biaxial plastic strain ratio, kb, is a function of the plastic strain as shown in Fig. 10. When the εpϕ to kϕ and εpϕ to ε p relations are known, the evaluation of the kϕ to ε p relation is straightforward.

4.4. Plastic strain ratios 4.5. Yield surface exponent The uniaxial plastic strain ratio is assumed to be constant for each material direction ϕ, i.e. kϕ ðε p Þ ¼ kϕ , since the serrated plastic flow in the experiment is ignored, see Fig. 3. A linear function, εpw;ϕ ¼ kϕ εpl;ϕ þ mϕ , is fitted to the experimental data in order to

For the yield surface exponent, a, in Eqs. (4) and (5), the general recommendations are a¼ 6 and 8 for BCC and FCC crystal structures, respectively [17]. No such general recommendation has been

R. Rentmeester, L. Nilsson / International Journal of Mechanical Sciences 92 (2015) 259–268


Table 3 Two sets of anisotropy parameters fulfilling the same set of target values y00 ; y45 ; y90 ; yb ; k00 ; k45 ; k90 and kb.

Solution 1 Solution 2










0.8913 0.0624

1.0960 0.6690

1.1348 1.1034

1.0802 1.1747

1.0474 0.5170

0.9856 1.0892

1.0018 0.9696

1.1061 1.6782

8 8

found for multiphase materials [18]. Larsson and Nilsson [14] proposed a linear mixture where a ¼8 was used for the austinitic fraction, and a ¼6 was assumed for the martensitic fraction. However, a ¼8 (constant) is assumed in the present work due to the mainly austenitic initial microstructure. 4.6. Evaluation of yield surface functions Given a set of target values: y00 ðε p Þ, y45 ðε p Þ, y90 ðε p Þ, yb ðε p Þ, k00 ðε p Þ, k45 ðε p Þ, k90 ðε p Þ and kb ðε p Þ, the functions A1 ðε p Þ; …; A8 ðε p Þ are found from the following optimisation problem: find

A1 ðε p Þ; …; A8 ðε p Þ


eA ¼


∑ ϕ

ðyϕ ðε p Þ  r~ ϕ ðε p ÞÞ2 þ


∑ ϕ

ðkϕ ðε p Þ  k~ ϕ ðε p ÞÞ2

þ ðyb ðε p Þ  r~ b ðε p ÞÞ2 þ ðkb ðε p Þ  k~ b ðε p ÞÞ2


where the superscript ð~Þ denotes a quantity evaluated from the yield criterion with a set of function values A1 ðε p Þ…A8 ðε p Þ. For a given set of target values, i.e. y00 ; y45 ; y90 ; yb ; k00 ; k45 ; k90 and kb, the resulting error eA becomes small ðeA  10  12 Þ. The standard Matlab [19] function fmincon was used in order to find the optimal function values Ai ðε p Þ with an initial guess corresponding to isotropy, i.e. A01 ¼ ⋯ ¼ A08 ¼ 1, where the superscript 0 indicates the initial guess at ε p ¼ 0. The final function values were used as initial guess in the next increment. Like the yield stress ratios, the anisotropy functions are assumed to vary linearly within these equivalent plastic strain increments. Observe that a linear variation in the anisotropy functions within an increment may result in a slight non-linearity of the resulting yield stress ratios. However, this possible non-linearity has been found to be negligible, since small equivalent plastic strain increments, i.e. Δε pA ¼ 0:001, have been used. Note also that different plastic strain increments can be used in the calibration of the yield stress ratios and the calibration of the anisotropy functions, but Δε pA r Δε py is recommended. It has been observed that more than one such set of anisotropy function values may exist, see Section 4.6.1. The first solution of anisotropy functions, denoted as Solution 1, is shown in Fig. 11, and the corresponding evolution of the yield locus is shown in Fig. 12. 4.6.1. Uniqueness of yield surface parameters It was found that the function values obtained depend on the initial guess, A0i . Two solutions at ε p ¼ 0:05 are shown in Fig. 13. The corresponding anisotropy function values are given in Table 3. Whether more than two or exactly two possible solutions exist has not been further investigated. Since both sets of anisotropy functions given in Table 3 represent the experimental target values included in the fitting procedure (19), neither of the two sets of parameter values can be directly rejected. Furthermore, σ is convex for m Z1, for all values of Ai. Thus, additional experiments are needed in order to choose the correct parameter set. Such additional experimental data may be obtained from a uniaxial tensile test in an additional direction. Here, such additional experimental data have been excluded from the parameter fitting but will be used in order to validate the evolution of the yield surface. Other experiments may be used for choosing the correct parameter values, e.g. shear tests and plane strain tests. However, such tests often tend to

Fig. 14. Predicted and experimental stress–strain relations in the ϕ¼ 151, 301, 601 and 751 directions: (a) Solution 1 and (b) Solution 2.

result in a complex stress and strain distribution and the evaluation of these experimental results often requires inverse analyses [5].

5. Numerical implementation The material model was implemented as a user defined material subroutine in LS-DYNA [20]. This section describes the numerical implementation. Rotations of tensors from the global to the corotated element coordinate system, and back, are performed before and after the stress update, respectively. Henceforth, the corotational superscript ð^Þ will be excluded and all subsequent constitutive relations will be related to this configuration. The material model was implemented for use with a shell element with a through thickness stress component, in which the normal stress components, i.e. σ11 and σ22 in Eq. (5), were substituted with σ n11 ¼ σ 11  σ 33 and σ n22 ¼ σ 22  σ 33 , respectively. The stress update algorithm is based on the tangent cutting plane, TCP, algorithm, see [21]. A brief description of the stress update algorithm and the inclusion of distortional hardening can also be found in [9]. The stress tensor is evaluated by assuming an elastic strain increment:

σ trial n þ 1 ¼ σ n þ C : dε



R. Rentmeester, L. Nilsson / International Journal of Mechanical Sciences 92 (2015) 259–268

Fig. 15. Predicted and experimental k-values using (a) the proposed effective stress functions and (b) the alternative functions.

The yield criterion (3) is evaluated at the trial stress according to ( r0 elastic p p ; ε Þ  σ ð ε Þ f ¼ σ ðσ trial ð21Þ y nþ1 n n 40 plastic deformation If the elastic condition is fulfilled, the stress update is finished and p p thus σ n þ 1 ¼ σ trial n þ 1 and ε n þ 1 ¼ ε n . Otherwise, plastic flow occurs and the plastic strain increment, Δε p ¼ ε pn þ 1  ε pn , has to be found. It is assumed that   ∂σ  ∂σ  p  ð22Þ pðkÞ ðkÞ σ ; ε  ∂σ n þ 1 n þ 1 ∂σ σ n þ 1 ;ε n þ 1 where an iteration index, k, has been introduced and where

pð1Þ p trial σ ð1Þ n þ 1 ¼ σ n þ 1 and ε n þ 1 ¼ ε n . The evaluation of the second order derivatives of the equivalent stress function, σ , with respect to the stress components, σ , is hereby avoided and the equation to be

solved in each iteration is reduced to a single scalar equation. Eq. (21) is linearised according to  ∂f  ðk þ 1Þ ðkÞ f n þ 1 ¼ f n þ 1 þ p  dε pðkÞ ð23Þ nþ1 ¼ 0 ∂ε ε pðkÞ nþ1

where the current plastic strain, ε


pðkÞ nþ1


p nþ


∑ dε


pðiÞ nþ1


p nþ


pðkÞ n þ 1,

is evaluated according to



and the derivative of the flow function, f, with respect to the equivalent plastic strain, ε p , is evaluated from       8 ∂σ  ∂f  ∂σ  ∂σ  ∂σ y  ∂Ai   : C :  þ ∑ p  pðkÞ ¼  p p   ∂σ σ ðkÞ ∂σ σ ðkÞ ∂ε ε pðkÞ ∂ε ε i ¼ 1 ∂Ai σ ðkÞ ∂ε ε pðkÞ nþ1






ð25Þ where the last term is the contribution from the distortional hardening. The anisotropy functions were evaluated separately and were provided to the stress update algorithm as piecewise linear relations between Ai and ε p . The consistency condition is then re-checked with the updated stress and equivalent plastic ðk þ 1Þ strain. If the flow function is within the tolerance, i.e. jf n þ 1 jr f tol , 1 the iteration procedure has finished and ðÞn þ 1 ¼ ðÞkn þ . þ 1 Otherwise a new plastic strain increment is evaluated. The tolerance f tol ¼ 10  4  σ was used. Normally, the iteration procedure converges within less than 7 iterations, but most often in less than 4.

6. Validation The uniaxial stress–strain relations in four intermediate directions and the predicted force in shear tests in two material directions were

Fig. 16. Extrapolated stress–strain relations.

compared to experimental data in order to validate the distortional hardening and distinguish between the two obtained solutions. 6.1. Uniaxial stress–strain relations The uniaxial stress–plastic strain relations in four directions

ϕ ¼151, 301, 601, and 751 with respect to the RD were predicted

and compared to the corresponding experimental data, see Fig. 14. In addition, the plastic strain ratios in the additional directions were predicted, see Fig. 15. It is concluded that the first set of anisotropy functions, i.e. Solution 1, gives a good prediction both of the stress–plastic strain relations and of the plastic strain ratios. The second set, i.e. Solution 2, gives a poor prediction of the stress–plastic strain relation in the directions ϕ ¼151 and 301. Moreover, the predicted variation of the plastic strain ratio with ϕ, as obtained with Solution 2, is poor, see Fig. 15b. 6.2. Shear test Unlike uniaxial tensile tests, the stress and strain distributions in shear tests are inhomogeneous over the cross section. Therefore, FE simulations of the shear tests were conducted and compared to the corresponding experimental data. The shear test was simulated in LS-DYNA [20]. The FE model shown in Fig. 5(b) was used. Four node shell elements with thickness stress and stretch were used, see [22]. The strain obtained in shear tests are significantly larger than the plastic strains at diffuse necking. Therefore, there is a need to extrapolate the stress–strain relations for larger strains. In this work, a third powerlaw hardening function was appended to the hardening function, see Eq. (8), i.e.

σ y;ϕ ¼ G3;ϕ þ K 3;ϕ ðεpϕ Þn3 ;

εpu;ϕ r εpϕ


R. Rentmeester, L. Nilsson / International Journal of Mechanical Sciences 92 (2015) 259–268


Fig. 17. Predicted and experimental relations between the nominal shear stress and displacement relations from the shear tests, in the directions (a) ϕ¼ 451 and (b) ϕ¼ 901.

where εpu;ϕ is approximately the limit strain of uniform elongation in the ϕ-direction. A value of 0.3 has been used in all three directions. Three constraints are required in order to evaluate the three constants in each direction, of which C1 continuity at εpϕ ¼ εpu;ϕ gives two. A third constraint was achieved by σ100, which denotes the yield stress at 100% plastic strain. This is a constant that must be found by other experiments than uniaxial tensile tests. In [14] it was found that a plausible value of σ100 is 2000 MPa for the present material. This value has been used herein in all three material directions. The obtained extrapolated stress–strain relations are shown in Fig. 16. Fig. 17 shows the predicted and experimentally found relations between the nominal shear stress and the elongation of the shear test specimens. Although the extrapolations of the uniaxial stress– strain relations to a large extent depend on the parameter σ100, it is obvious that Solution 2 over-predicts the nominal stress, while Solution 1 gives a good prediction in both material directions.

7. Conclusions This work presents a calibration methodology for mixed isotropic-distortional hardening of an anisotropic yield function. The yield stress ratios were evaluated at a finite number of plastic strain levels, with an assumed linear variation in between. Using this method, the stress–strain relations in the directions included in the calibration were represented accurately. In a similar manner, the anisotropy functions were evaluated at a number of plastic strain levels. The methodology was used together with the eight parameter effective stress YLD2003, but it can be used with any effective stress function. The specific effective stress function is able to represent three uniaxial yield stresses, the balanced biaxial yield stress, and the four corresponding plastic strain ratios, exactly. Experimental data at these stress states were used in the calibration procedure, whereas uniaxial tensile test data and shear test data in four additional material directions were used for validation. Two sets of anisotropy functions were found, both of which fulfill the experimental target values. The corresponding yield loci are different but convex. Thus, none of them could be directly rejected. The stress–strain relations and the corresponding plastic strain ratios in four additional material directions, ϕ ¼151, 301, 601 and 751, were predicted and compared to the corresponding experimental data. One of the anisotropy function sets resulted in predictions close the experimental data, whereas the second set resulted in a less accurate prediction and was therefore rejected. The same conclusion was drawn from the comparison between the numerical and the experimental results in the shear tests. It was

concluded that the yield stress in at least one additional stress state should be compared to experimental data in order to validate the obtained anisotropy parameters. Regarding the set of anisotropy functions giving the best predictions, the numerical stress–strain ratios were close to the experimental findings. However, there was a slight discrepancy between the predicted and experimentally observed plastic strain ratios. It is the authors' belief that this discrepancy may be partly caused by uncertainties in the experimental results caused by the serrated plastic flow. Despite the discrepancy between the predicted and measured plastic strain ratios, the variation of these ratios with respect to the angle ϕ seems to be plausible. It is to be noticed that it is the effective stress function, σ , itself that determines whether the prediction of the plastic strain ratios in the ϕ ¼151, 301, 601 and 751 directions is accurate or not, since these are interpolations of the mixed isotropic-distortional hardening in each test direction. Since the predicted k-values do not change significantly with plastic strain, it is also concluded that this lack of accuracy is independent of the distortional hardening, but depends on the effective stress function itself. It was also shown that the predicted stress–strain relations will coincide with the target values at all equivalent plastic strain levels at which the corresponding yield stress ratios have been evaluated. The smaller the equivalent plastic strain increments, the closer the relation between the predicted stress–plastic strain relation will become.

Acknowledgments The work presented in this paper has been conducted with funding from the SFS ProViking project Super Light Steel Structures. Outokumpu Stainless and Dr. Ramin Moshfegh are greatly acknowledged for advice and material supply. References [1] Barlat F, Aretz H, Yoon J, Karabin M, Brem J, Dick R. Linear transformationbased anisotropic yield functions. Int J Plast 2005;21:1009–39. [2] Banabic D, Barlat F, Cazacu O, Kuwabara T. Advances in anisotropy and formability. Int J Mater Form 2010;3:165–89. [3] Ortiz M, Popov E. Distortional hardening rules for metals plasticity. J Eng Mech 1983;109:1042–57. [4] Axelsson K. On constitutive modelling in metal plasticity with special emphasis on anisotropic strain hardening and finite element formulation [Thesis]. Sweden: Department of Structural Mechanics, Chalmers University of Technology; 1979. [5] Larsson R, Björklund O, Nilsson L, Simonsson K. study of high strength steels undergoing non-linear strain paths – Experiments and modelling. J Mater Process Technol 2011;211:122–32. http://dx.doi.org/10.1016/j.jmatprotec.2010.09.004.


R. Rentmeester, L. Nilsson / International Journal of Mechanical Sciences 92 (2015) 259–268

[6] Barlat F, Gracio JJ, Lee M-G, Rauch EF, Vincze G. An alternative to kinematic hardening in classical plasticity. Int J Plast 2011;27:1309–27. [7] Barlat F, Brem JC, Yoon J-W, Chung K, Dick RE, Lege DJ, Pourboghrat F, Choi S-H, Chu E. Plane stress yield function for aluminium alloy sheets—Part I: theory. Int J Plast 2003;19:1297–319. [8] Aretz H. Numerical analysis of diffuse and localized necking in orthotropic sheet metals. Int J Plast 2007;23(5):798–840. [9] Aretz H. A simple isotropic-distortional hardening model and its application in elastic-plastic analysis of localized necking in orthotropic sheet metals. Int J Plast 2008;24:1457–80. [10] Plunkett B, Lebensohn R, Cazacu O, Barlat F. Anisotropic yield function of hexagonal materials taking into account texture development and anisotropic hardening. Acta Mater 2006;54(16):4159–69. [11] Aretz H. Applications of a new plane stress yield function to orthotropic steel and aluminium sheet metals. Model Simul Mater Sci Eng 2004;12:491–509. [12] Aretz H, Keller S. On the non-balanced biaxial stress state in bulge-testing. Steel Research International 2011 Special Edition; 2011. p. 738–43. [13] Schedin E, Prentzas L, Hilding D. Finite element simulation of the trip-effect in austenitic stainless steel. In: SAE World Congress 2004, paper 2004-01-0885; 2004. [14] Larsson R, Nilsson L. On the modelling of strain ageing in a metastable austenitic stainless steel. J Mater Process Technol 2012;212:46–58.

[15] Hänsel AHC, Hora P, Reissner J. Model for the kinetics of strain-induced martensitic phase transformation at non-isothermal conditions for the simulation of sheet metal forming processes with metastable austenitic steels. In: Huétink, Baaijens (Eds.), Simulation of materials processing: theory, methods, and applications, Rotterdam; 1998. p. 373–8. [16] Lindgren L-E, Olsson M, Carlsson P. Simulation of hydroforming of steel tube made of metastable stainless steel. Int J Plast 2010;26(11):1576–90. [17] Hosford WF. The mechanics of crystals and textured polycrystals. New York: Oxford Science Publications; 1993. [18] Aretz H, Hopperstad OS, Lademo O-G. Yield function calibration for orthotropic sheet metals based on uniaxial and plane strain tensile tests. J Mater Process Technol 2007;186:221–35. [19] MATLAB, version 7.5. Natick: The Mathworks Inc. 2007. [20] Hallquist J. LS-DYNA theory manual. Livermore: Livermore Software Technology Corporation; 2006. [21] Ortiz M, Simo J. An analysis of a new class of integration algorithms for elastoplastic constitutive relations. Int J Numer Methods Eng 1986;23:353–66. [22] Hallquist J. LS-DYNA keyword user's manual. Livermore: Livermore Software Technology Corporation; 2007.