An improved yield criterion characterizing the anisotropic and tension-compression asymmetric behavior of magnesium alloy

Zhigang Li,Haifeng Yang,Jianguang Liu,Fu Liu

a School of Mechanical,Electronic and Control Engineering,Beijing Jiaotong University,Beijing,100044,PR China

b Automotive Research Institute of Beijing Automotive Co.,Ltd.,Beijing 101300,China

cBeijing Key Laboratory of Civil Aircraft Structures and Composite Materials,COMAC Beijing Aircraft Technology Research Institute,Beijing,102211,PR China

d Shanghai Aircraft Design and Research Institute,Shanghai 201210,PR China

Abstract A novel yield criterion based on CPB06 considering anisotropic and tension-compression asymmetric behaviors of magnesium alloys was derived and proposed (called M_CPB06).This yield criterion can simultaneously predict the yield stresses and the Lankford ratios at different angles (if any) under uniaxial tension,compression,equal-biaxial and equal-compression conditions.Then,in order to further describe the anisotropic strain-hardening characteristics of magnesium alloy,the proposed M_CPB06 criterion was further evolved to the M_CPB06ev model by expressing the parameters of the M_CPB06 model as functions of the plastic strain.As the model was developed,the stresses and Lankford ratios of AZ31B and ZK61M magnesium alloys at different angles under tensile,compressive and through-thickness compressive conditions were used to calibrate the M_CPB06/M_CPB06ev and the existing CPB06ex2 model.Calibration results reveal that compared with the CPB06ex2 yield criterion with equal quantity of coefficients the M_CPB06 criterion exhibits certain advancement,and meanwhile the M_CPB06ev model can relatively accurately predict the change of the yield locus with increase of the plastic strain.Finally,the M_CPB06ev model was developed through UMAT in LS-DYNA.Finite element simulations using the subroutine were conducted on the specimens of different angles to the rolling direction under tension and compression.Simulation results were highly consistent with the experimental results,demonstrating a good reliability and accuracy of the developed subroutine.

Keywords: Magnesium alloy;Anisotropy;Tension-compression asymmetry;Improved yield criterion;UMAT.

1.Introduction

Owing to the increasing energy demand and environmental protection pressure,lightweighting of industrial product is of great significance In areas such as aeronautics,astronautics,and automobiles etc.,utilization of lightweight alloys is a key option to realize structural lightweighting[30].In virtue of the low density and high specific strength,magnesium alloys are regarded as lightweight materials with great potential for wide application [16].However,magnesium alloys exhibit a unique hexagonal close packed (HCP) crystal structure,which is different from face-centered cubic (FCC) or body-centered cubic(BCC) crystal structures of conventional metals.As a result,magnesium alloy exhibits both texture-induced anisotropy and tension-compression asymmetry [12].

Finite element simulation has become an important method for predicting the yielding and strain hardening of structures.The accuracy of the numerical simulation is generally dependent on the used material constitutive yield model.Over the past decades,various yield criteria have been proposed and improved,including conventional isotropic yield criteria (e.g.,von Mises,Tresca) and typical anisotropic yield criteria (e.g.,Hill1948,Barlat1989,YLD2000-2D) [1].However,these criteria are only applicable for FCC and BCC crystal materials such as steel and aluminum alloys,but not for asymmetric mechanical behavior of HCP materials (e.g.magnesium alloys) under tension and compression.In order to characterize the tension-compression asymmetry besides anisotropy,Hosford [7] introduced hydrostatic pressure term into the Hill48 anisotropic yield criterion and thereafter this method has been widely used by Liu et al.[14],Stoughton and Yoon [25],Kim et al.[11],and Lou et al.[15] who proposed a common model.Yoon et al.[29] proposed a yield criterion based on the stress invariants and extended the criterion to anisotropy by using two different linear transformations.Hu et al.[8]and Tang et al.[26] proposed new tension-compression asymmetric yield criteria and extended them to anisotropy using linear transformation methods.

Aforementioned criteria describing the tensioncompression asymmetry were proposed by introducing a hydrostatic pressure term.In addition,Cazacu and Barlat[4] modified the Drucker yield criterion by changing the exponent terms of the equation to odd numbers to fulfill the asymmetry of the yield locus.In succession,Cazacu and Barlat [4] introduced more coefficient to invariants of stress deviator to change the criterion from isotropy to anisotropy;and Nixon et al.[17] also changed the model to anisotropy through performing linear transformation on the second and third invariants of stress deviator.Cazacu et al.[5] further proposed a typical anisotropic and tension-compression asymmetric yield criterion (CPB06),which was applied to magnesium and titanium alloys.

In order to improve the prediction accuracy of the CPB06 criterion,Plunkett et al.[20] put forward a Cazacu-Plunkett-Barlat (CPB)-type yield function using linear transformation of the stress deviator (e.g.CPB06ex2).Raemy et al.[23] proposed an anisotropic and tension-compression asymmetric yield criterion based on Fourier series.Khan et al.[9] proposed a yield criterion related to the strain rate and temperature,and meanwhile this criterion can describe both the anisotropy and tension-compression asymmetry.Recently,a new criterion is proposed to describe the anisotropy and tension-compression asymmetry by increasing two additional terms referred to as scaling and asymmetry functions [18].

In addition to accurate representation of the initial yield of the material,anisotropic hardening is also an important issue requiring attention [19].The conventional isotropic hardening model cannot describe this phenomenon.Therefore,Kim et al.[10] described the hardening of HCP materials by combinations of CPB06 yield criterion with two different hardening models,and the prediction results were verified against experimental data under three-point bending condition.Plunkett et al.[[21],[22]] expressed the parameters that describe the anisotropy in CPB06 criterion as an equation of plastic strain to describe the variation of the hardening curves under different strain paths.The results revealed that the modified yield criterion could predict the hardening characteristics of HCP materials.Recently,the anisotropic hardening behavior of AZ31B [[27],[28]],TA-6 Vtitanium alloys [6],highly pure titanium alloys [17] and ZM21 magnesium alloy [24] were characterized using a method similar to the one proposed by Plunkett et al.[[21],[22]].

In summary,developing a constitutive model that can well describe the yielding and hardening characteristics of magnesium alloys has been an important topic.As a typical criterion,the aforementioned associated CPB06 criterion can characterize the behaviors of anisotropy and tension-compression asymmetry.However,studies revealed that the description on the plastic behavior of AZ31B by this criterion was still unsatisfying.Therefore,it is necessary to develop a new yield criterion to adequately describe the anisotropic and tensioncompression asymmetric yielding and strain-hardening characteristics of AZ31B magnesium alloy.

2.Development of an improved CPB06 model

In this section,the typical CPB06 yield criterion was modified to a new model to improve the prediction accuracy on the yield stress and the Lankford ratio of AZ31B magnesium alloy under different conditions.The CPB06 is an anisotropic yield criterion that can describe the tensioncompression asymmetry [5],with the equation as:

wherekis the coefficient describing tension-compression asymmetry andais the order of the homogeneous equation.-1≤k≤1 anda≥1 should be satisfied in order to guarantee the convexity of the yield locus.is the effective stress;βis a constant.Σ1,Σ2,Σ3are the three principal values of the transformation tensor (Σ) by linear transformation of the stress deviator.Σcan be obtained according to Eq.(2).

whereCis a fourth-order tensor of matrix.Relative to the orthotropic axes (x,y,z),the matrixCcan be expressed as Eq.(3).Regarding to the rolling sheet,the axis ofx,y,zrefers to the rolling direction,width direction and thickness direction of the sheet,respectively.

In plane stress condition,components ofΣcan be expressed as:

Fig.1.Normalized yield stresses and Lankford ratios of AZ31B magnesium alloy at different strains under different conditions.

Fig.2.Calibration results of the hardening curves by Ludwik equation.

According to Eq.(1),the yield criterion can be regarded as an expression of effective stress as:

In this study,the linear transformation of stress tensor was employed to obtain a new model.The prediction accuracy of the yield stress and the Lankford ratio can be improved by increasing the quantity of model coefficients Meanwhile,this method did not alter the convexity of the yield surface.Based on above description,a new expression of effective stress can be obtained using different linear transformation matrix (C′)(Eq.(7)) and coefficient (k′) in the CPB06 as Eq.(8).

As different transformation matrix (C′) was used,CijinΣ1,Σ2,Σ3was substituted byto obtain expressions of

Based on the method proposed by Bron and Besson [3],the effective stress in the yield criterion can be define as Eq.(9).The convexity of the yield criterion can be kept using this kind of accumulation function.

whereσkdenotes a function with convexity;αkis a positive coefficient and=1;bis a constant.

In Eq.(9),letb=1 andn=2.Substituting the Eqs.(6) and (8) into Eq.(9) to obtain a new expression of yield criterion,which was define as M_CPB06.Because two identical forms of functions of effective stress were used,the effect by coefficienαkwas negligible.To simplify the equation,αkis assumed to be 1.The final expression of M_CPB06 is as Eq.(10).

whereλis a non-independent parameter that can be represented by other parameters via uniaxial tensile tests.This equation contains totally 21 coefficients including 18 coefficient inCandC′,the coefficientkandk’that describe the asymmetry of tension-compression and a coefficiena.

In the state of plane stress,onlyσxx,σyy,andσxyare not 0,therefore,four coefficient ofC44,C55andC′44,C′55would not work and 17 coefficient were left to be determined in the equation.

As specimens were obtained at different angles with the rolling direction,the coordinate system of the specimen should be transformed to the original sheet coordinate system as the specimen was described by an anisotropic yield criterion.σθrefers to the uniaxial yield stress with an angle ofθto the rolling direction.Under the condition of uniaxial loading,the component of stress tensor can be expressed as[1]:

Substituting Eq.(11) into Eq.(10) and taking the yield stress (σ0) along the rolling direction as the reference stress for normalization,the expressions of the uniaxial tensile and compressive yield stresses (and) of M_CPB06 yield criterion at different angles can be obtained as Eq.(12).

The yield stress under equal-biaxial tension and equalbiaxial compression can be computed as:

According to the definition of the Lankford ratio and the associate flow criterion,the Lankford ratio of specimens with different angles to the rolling direction under tensile and compressive conditions can be calculated by:

Under equal-biaxial tension,the Lankford ratio can be calculated by:

In Eqs.(18) and (19),the partial derivative of yield stress to stress component can be calculated based on the chain derivation rule [20] as:

3.Calibration of the improved CPB06 model (M_CPB06)

3.1.Experimental data of magnesium alloys

To calibrate M_CPB06 model,and evaluate the advantage of the improved model on the applicability of magnesium alloys,the experimental data of AZ31B and ZK61M alloys under uniaxial tension,uniaxial compression,and throughthickness compression tests from our previous study[13]were used.Considering the strain of ZK61M under compression at an angle of 90° is very short (only about 4%?5%) and only limited data can be used to calibrate M_CPB06,therefore,in this section,only the data of AZ31B at different strains and angles was presented and the data of ZK61M at limited strain was shown in theDiscussionsection to further examine the advantage of the proposed model.

The normalized yield stresses and the Lankford ratios at plastic strains of 0,1,3,5 and 7% at different angles relative to the sheet rolling direction under tension,compression,and through-thickness compression were extracted,and the results are shown in Fig.1.The plastic strain of AZ31B at 90° to the rolling direction is only about 8%,thus the plastic strain up to 7% was used to calibrate the model.To reflect the capability of the yield criterion in characterizing the yield stress under equal-biaxial compression,an equivalent method reported by Yoon et al.[29] was used to calculate the yield stress under equal-biaxial compression,according to the formula of/4,in whichare the yield stresses at an angle of 0,45° and 90° to the rolling direction,respectively,under uniaxial compression.

In order to describe the hardening characteristics of AZ31B,the hardening curves at different angles to the rolling direction were fitte using the Ludwik equation(see Eq.(24)).Nin the hardening equation can reflect the variation trend of the curve with strain.As observed in Fig.2,the hardening equation fit the experimental data fairly well.Additionally,the fittin results demonstrate thatn,which reflect the hardening characteristics,is dependent on the angle and the difference among different angles was relatively large,suggesting that AZ31B does not follow the isotropic hardening law.

whereA,B,andnare the coefficient of the Ludwik equation.

3.2.Calibration on the M_CPB06 yield model

In this study,the optimization was conducted using MATLAB.The combination of genetic algorithm and single point algorithm was used.First,rapid global optimization was achieved using genetic algorithm to obtain initial optimized results.Then,the results obtained by the genetic algorithm were used as inputs of single point algorithm (fminsearch function in MATLAB) to obtain the final values of the coefficients This method can combine the advantages of both algorithms.First,genetic algorithm is an intelligent method,which requires no manual assignment of initial values,and it can fulfill rapid global optimization,while single point algorithm can achieve further optimization of problems with reasonable initial values.Tests revealed that the optimized result by this combination method was superior to that of either algorithm alone.The optimized results of the coefficient in the M_CPB06 model at plastic strains of 0,1,3,5 and 7%are listed in Table 1.

Table 1The calibrated coefficient values of the M_CPB06 model.

Table 2The calibrated coefficient values of the M_CPB06ev model.

Fig.3 showed the prediction results of the yield stress and the Lankford ratio of AZ31B obtained based on the calibrated coefficient and the comparison with experimental data.As observed,at all the plastic strains of 0,1,3,5 and 7%,the M_CPB06 yield model exhibited favorable consistency in its prediction of yield stresses and Lankford ratios at different angles (if any) with the experimental results under uniaxial tension,uniaxial compression,equal-biaxial tension and compression conditions.

4.Anisotropic hardening of AZ31B

4.1.Evolution of the yield criterion on describing the anisotropic strain-hardening behavior of AZ31B

For the materials with isotropic hardening behavior,the calibration on the model can be achieved using initial yielding data or the yielding data at any plastic strain.However,as aforementioned,the hardening curves of AZ31B at different angles varied with the increase of the plastic strain.The Lankford ratio followed the same trend.Therefore,the proposed yield criterion (M_CPB06) was further evolved to try to describe the variation of the flow stress and the Lankford ratio with plastic strain at different angles.To achieve this point,Tari et al.[27] employed the coefficient ofCijandkof the Cazacu-Plunkett-Barlat (CPB)-type yield model to represent by functions of plastic strain.In this study,Cijwas denoted by a second-order exponential function,different from the function used in Tari et al.[27].This change can keep the advantage ofCijin Tari et al.[27] (continuous,differentiable and stable) and meanwhile it can improve the flexibility of the original expression.The expression of the coefficienkwas kept the same with that in Tari et al.[27].

Fig.3.Prediction of the M_CPB06 model on the yield stresses and Lankford ratios at different plastic strains.

Fig.3.Continued

According to the calibration in Section 3.2,it is found thatain the yield criterion was not a constant parameter at different strains.Therefore,awas regarded as a function of plastic strain.The same form of equation askwas adopted fora.kandain Eq.(26) increased monotonically with plastic strain,thus the values ofkandaat all strains can be easily confine in appropriate ranges during optimization and meanwhile the yield locus can be retained convexity.whereprefers to effective plastic strain andaij,bij,cij,dij,e,f,g,j,h,n,mare material parameters.

Substituting Eq.(26) into the yield criterion (Eq.(11)),then an evolutional model accounting for the hardening with plastic strain was obtained,denoted as M_CPB06ev.

4.2.Calibration on the evolutional yield criterion(M_CPB06ev)

As the evolutional yield criterion (M_CPB06ev) was obtained,the new coefficient were calibrated using the optimization method similar with the one in Section 3.2.The optimization target is minimizing the error equation as Eq.(27).Different from Eqs.(25)and(27)is a summation form,which means optimizing the coefficient at five strains of 0,1,3,5 and 7% simultaneously.In other words,the M_CPB06ev model needs to predict a large amount of experimental data of yield stress and the Lankford ratio than that of the M_CPB06 model.Considering the prediction ability of M_CPB06 on yield stress at specific strain under equal-biaxial compression has been verified and the data under this condition was not obtained experimentally.Therefore,in optimizing the coeffi cients of M_CPB06ev model,the data was removed to guarantee the prediction accuracy of M_CPB06ev on other experimental data points.

wherekindicates the plastic strain of 0,1,3,5 and 7%;andirepresents the angle of 0°,22.5°,45°,67.5°,and 90°.

Table 2 listed the coefficient obtained by optimization.Figs.4 and 5 illustrated the hardening curves and the Lankford ratios at different angles predicted by M_CPB06ev model and the experimental results.As observed,the evolutional model can well predict the variation of the stress with strain under tension and compression at all five angles (Fig.4).For the prediction on the Lankford ratio curves,in general,the results predicted by the evolutional model were fairly agreement with the experiments,except certain deviation was observed at the angle of 67.5° under tension.This is because,as determining the model coefficients the corresponding between the experimental data and the model prediction at five plastic strains were optimized simultaneously,resulting in a large challenge in balancing the model predictions.The prediction results for the equal-biaxial state were not shown because the yield stress has already clearly shown in the yield locus as Fig.6,which showed the yield loci predicted by the M_CPB06ev model at different plastic strains.As observed,the M_CPB06ev model can well predict the yield stresses under different conditions and characterize the shape change of the yield locus with plastic strain.

Fig.4.Prediction of the M_CPB06ev model on the strain hardening curves at different angels.

Fig.5.Prediction of the M_CPB06ev model on the Lankford ratio curves at different angels.

Fig.6.The yield loci at different plastic strains predicted by the M_CPB06ev model (σxy=0).

5.Development of the UMAT for M_CPB06ev model and simulation

5.1.Development of the UMAT

The UMAT for the proposed M_CPB06ev criterion was developed and it was coded into LS-DYNA.At each step,the main program of the LS-DYNA is responsible for inputting the element strain increment,material constants and the historical variables of the last step into the subroutine.The subroutine outputs the element stress component and the historical variables of this step,and returns to the main program.The flow chart of the UMAT is shown in Fig.7.First,at elastic stage,the trial stress was calculated according to the stiffness matrix and then,it was checked whether it yielded.If not,the subroutine would update the current stress and historical variable values and wait for the next call;otherwise,the subroutine would enter the plastic stage and calculate the plastic flow direction based on the used flow rule as well as the plastic parameter increment through numerical iteration.Finally,the convergence was checked.If the result was convergence,the subroutine would update the stress and return to the yield surface;otherwise,it will direct update the historical variables and wait for the next call.

Fig.7.The flow chart of the UMAT.

In this study,the stress integration method used the return algorithm,which mainly consists of elastic prediction step and plastic adjustment step that made the trial stress beyond yield surface return to yield surface.At each time step,the developed yield criterion of M_CPB06ev was used to determine whether the material has yielded.If Eq.(29) is satisfied the stress and relevant historical variables are updated.If not,the yielding has reached and the trial stress should return to the yield surface.

At the plastic adjustment step,the associate flow criterion was employed as Eq.(30).

Eq.(30) can be kept at a very small value (less than 10-6)by solving the increment of the plastic parameter,which can be calculated according to Eq.(31) [27].

5.2.Results of simulation

In order to verify the reliability of the developed subroutine,simulations were conducted on the specimens under tensile and compressive conditions.Simulation results of the stress-strain curves and the Lankford ratios at different angles were compared with the prediction results by the theoretical model as shown in Figs.8 and 9.As observed,the simulation results were highly consistent with the theoretical results in both stress-strain curves and the Lankford ratios,demonstrating a good reliability of the subroutine.

Fig.8.Comparison of hardening curves predicted by the theoretical model (M_CPB06ex) and the simulation using subroutine.

Fig.9.Comparison of the Lankford ratio curves predicted by the theoretical model (M_CPB06ex) and the simulation using subroutine.

6.Discussion

To verify the advantage of the proposed M_CPB06 ctierion,the current M_CPB06 model was compared with the classical CPB06ex2 model,both of which has 17 coefficients The CPB06ex2 model is also an improved model based on CPB06 [20].In calibrating the coefficient of the CPB06ex2 model,the same optimization method as M_CPB06 was used.The errors of the combined yield stress and Lankford ratio(Eq.(25)) between the experiments and the prediction results by the CPB06ex2 and M_CPB06 models at different strains are shown in Fig.10.It is seen that the error between the experiment and that predicted by the current M_CPB06 model was significantly lower than that by CPB06ex2.In order to facilitate showing the comparison between CPB06ex2 and M_CPB06,the results of the yield stresses and the Lankford ratios at different angles at plastic strains of 0,1,3,5 and 7% predicted by these two models were compared with the experimental data.To simplify,only the results at the plastic strain of 0 and 7% were presented as shown in Fig.11.As observed,both M_CPB06 and CPB06ex2 can simultaneously characterize the plastic behavior of AZ31B under tension and compression,however,the M_CPB06 model exhibited certain advantage in the prediction of yield stresses and the Lankford ratios at different angles over CPB06ex2.

Fig.10.The error between the experiment and the prediction by CPB06ex2 and M_CPB06 at different strains.

Fig.12 showed the yield loci predicted by the M_CPB06 and CPB06ex2 models and the experimental data.The yield stress points under equal-biaxial tension and equal-biaxial compression were also included.As observed,In general,the M_CPB06 exhibited advantage over CPB06ex2 in the prediction of yield stress under most conditions,especially under the condition of equal-biaxial compression.

In order to further examine the advantage of proposed M_CPB06 on the applicability of other magnesium alloys,the experimental data of ZK61M at the plastic strain of 2%from the study of Li et al.[13] was used to calibrate the models of M_CPB06 and CPB06ex2.The yield stresses and the Lankford ratios at different angles predicted by these two models are shown in Fig.13.As observed,M_CPB06 exhibited advantage over CPB06ex2 in predictions of both yield stress and the Lankford ratio.Fig.14 showed the yield loci predicted by these two models and the comparison with the experimental data.It is seen that the yield locus predicted by the M_CPB06 is slightly superior to that of the CPB06ex2.

According to above comparisons of M_CPB06 and CPB06ex2 on the prediction of plastic behavior of AZ31B and ZK61M magnesium alloys,in summary,the proposed M_CPB06 exhibits advantage in the prediction accuracy over CPB06ex2,both of which have the same number of parameters.

Fig.11.The experimental normalized yield stress and Lankford ratio and the prediction by the M_CPB06 and CPB06ex2 models.

Fig.12.Comparison of the yield loci of AZ31B predicted by the M_CPB06 and CPB06ex2 models.

Fig.13.The experimental normalized yield stress and Lankford ratio at the plastic strain of 2% and the prediction by the M_CPB06 and CPB06ex2 models.

Fig.14.Comparison of the yield loci of the ZK61M magnesium alloy predicted by the M_CPB06 and CPB06ex2 models.

In this study,some limitations still exist which need further study.(1) In calibrating and validating the proposed M_CPB06,due to the limitation by the small maximal strain value (approximate 8%?9%) at 90° under uniaxial tension,the stresses and Lankford ratios corresponding to relative small strains (0%,1%,3%,5%,and 7%) under different conditions were selected to calibrate the M_CPB06 model.For other magnesium alloys with higher ductility,this improved model should be also applied.once the stresses and Lankford ratios corresponding to large wide range of strains can be obtained,they will be used to calibrate the model in future study to further reflect the performance of the proposed model on predicting the variation of the yield locus with increase of the plastic strain.(2) The UMAT subroutine was only verified with the simple specimens.Next step,more tests on some components,such as cup drawing,three-bending tests etc.,should be conducted and further validate the accuracy of the subroutine under more complex conditions.

7.Conclusions

(1) Based on CPB06 criterion,a novel anisotropic and tension-compression asymmetric yield criterion called M_CPB06 criterion was proposed.Calibration results revealed that the M_CPB06 yield criterion can well predict the yield stresses and the Lankford ratios of AZ31B and ZK61M under different conditions.Additionally,M_CPB06 yield criterion exhibits advantage over the CPB06ex2 model with equal quantity of coefficients.

(2) The M_CPB06 yield criterion was further evolved to M_CPB06ev criterion to describe the anisotropic strainhardening characteristics of magnesium alloy,which can describe the variation law of the yield stress and the Lankford ratio with plastic strain.Calibration results showed that the M_CPB06ev criterion can adequately describe the shape change of the yield locus with plastic strain.

(3) The UMAT of the proposed M_CPB06/M_CPB06ev model was developed and coded into LS-DYNA.The developed subroutine was proven to demonstrate a good reliability and accuracy.

Acknowledgments

This study was supported by Beijing Natural Science Foundation (No.L201010),the United Fund of Ministry of Education for Equipment Pre-Research (Grant No.6141A02033121),and National Natural Science Foundation of China (Grant No.51975041).The authors would like to thank Yuanli Bai and Yueqian Jia from University of Central Florida for their kindness help.