Temperature-dependent constitutive modeling of a magnesium alloy ZEK100 sheet using crystal plasticity models combined with in situ high-energy X-ray diffraction experiment

Hyuk Jong Bong,Xiaohua Hu,Xin Sun,Yang Ren

aMaterials Deformation Department,Korea Institute of Materials Science,Changwon,Gyeongnam 51508,South Korea

b Energy and Transportation Science Division,Oak Ridge National Laboratory,Oak Ridge,TN 37831,USA

cX-ray Science Division,Argonne National Laboratory,IL 60439,USA

Abstract A multiscale crystal plasticity model accounting for temperature-dependent mechanical behaviors without introducing a larger number of unknown parameters was developed.The model was implemented in elastic-plastic self-consistent(EPSC)and crystal plasticity finite element(CPFE)frameworks for grain-scale simulations.A computationally efficient EPSC model was first employed to estimate the critical resolved shear stress and hardening parameters of the slip and twin systems available in a hexagonal close-packed magnesium alloy,ZEK100.The constitutive parameters were thereafter refined using the CPFE.The crystal plasticity frameworks incorporated with the temperature-dependent constitutive model were used to predict stress–strain curves in macroscale and lattice strains in microscale at different testing temperatures up to 200 °C.In particular,the predictions by the crystal plasticity models were compared with the measured lattice strain data at the elevated temperatures by in situ high-energy X-ray diffraction,for the first time.The comparison in the multiscale improved the fidelity of the developed temperature-dependent constitutive model and validated the assumption with regard to the temperature dependency of available slip and twin systems in the magnesium alloy.Finally,this work provides a time-efficient and precise modeling scheme for magnesium alloys at elevated temperatures.

Keywords:High-energy X-ray diffraction;Crystal plasticity finite element;Elastic-plastic self-consistent model;Twin;Temperature.?Corresponding author.

1.Introduction

The increasing demand for lightweight metals in vehicle body applications has garnered attention in recent years for magnesium alloys with a high strength-to-density ratio.Although magnesium alloys possess good tensile and fatigue properties,the formability of magnesium alloys at room temperature(RT)is poorer than that of other structural metals[1–3].This drawback hinders the practical application of magnesium alloys.

The major reason for the poor formability of magnesium alloys at RT originates from a lack of independent deformation modes to accommodate plastic deformation.The generally accepted plastic deformation mechanisms of magnesium alloys are basal,prismatic,and pyramidal slip systems,along with tensile twinning.The slip resistance of the basal slip is much lower than that of other slip and twin systems at RT;however,the basal slip alone cannot provide five independent slip systems,which are the minimum number of systems required for arbitrary plastic deformation according to the von Mises criterion[4].Therefore,the activation of other slip and twin systems is necessary to secure the ductility or formability of magnesium alloys and facilitate the application of magnesium alloys.It is well-knownthat the basal and tensile twin systems in magnesium are nearly temperature-independent while the non-basal slip systems in magnesium,prismatic and pyramidal,are significantly temperature-dependent.Therefore,one of the strategies to overcome the limited formability of the magnesium alloys can be temperature-assisted forming utilizing the diminished slip resistance of the non-basal slip systems[5,6].The significant activation of the non-basal slip systems at elevated temperatures offers other independent deformation modes beyond the basal slip and twin.Consequently,the activation of the non-basal slip systems leads to the enhanced formability.In addition to temperature-assisted forming,significant interest has been directed to the addition of rare-earth elements to magnesium alloys to reduce the unfavorable mechanical characteristics from the alloying perspective[7–11].The addition of rare-earth elements led to a weakened basal texture,which is commonly observed in AZ31 alloy and detrimental to formability.Therefore,the rare-earth elements added magnesium alloys exhibit enhanced formability compared to the AZ31 alloy[9].Particularly,ZEK100 alloy,which is a typical example of rare-earth element added magnesium alloys,exhibits enhanced formability at considerably lower temperatures than those required for AZ31[6,12].

For successful process design of the temperature assistedforming or alloy design for the initial texture modification,a fundamental understanding of the role of the deformation modes on the macroscopic mechanical behaviors is vital.To date,numerical models describing the complex mechanical responses of the magnesium alloys have been developed.The models include continuum-based phenomenological models[12–16],and grain-scale crystal plasticity models[1,2,17–25].The continuum-based models are practical in terms of computing efficiency and can be readily employed in actual metal forming simulations.On the contrary,the crystal plasticity models are physically driven and suitable for understanding the complex mechanical responses of the magnesium alloys in fundamental level.However,it demands much expensive computing cost than the continuum-based models,and the efforts should be put into solving the issue.

Numerous grain-scale crystal plasticity models for various structural metals[1,2,18–23,26–30]including hexagonal close-packed(HCP)metals with the heterogeneous deformation modes[1,2,17–25].Nevertheless,few theoretical attempts have been made to model magnesium alloys at elevated temperatures[2,18,31–34].Jain and Agnew[31]investigated the temperature-dependent mechanical behavior of AZ31B in the temperature range from RT to 250 °C.The strain hardening parameters for available slip/twin systems at each testing temperature were identified individually,and their variation as a function of temperature was discussed.A similar approach was employed in other studies[32,33].These works are sufficient to reveal that the non-basal slip systems are temperature-dependent,whereas the basal slip and tensile twin are nearly temperature-independent.However,this approach requires multiple processes to identify the hardening parameters of slip/twin systems at discrete temperatures.Systematic temperature-dependent constitutive models within crystal plasticity frameworks have also been proposed[2,18].In these approaches,the hardening parameters are explicitly expressed in temperature-dependent forms.The models were successful in predicting the macroscopic stress–strain data in the temperature range-200–150 °C[18]and RT to 250 °C[2].However,the complexity of a large number of constitutive parameters hinders practical application of the proposed models.Singh et al.[35]proposed a phenomenological polycrystal model incorporating temperature-dependent strain hardening for non-basal slip systems by introducing a scaling factor.The model reproduced the stress–strain behaviors of an HCPstructured zirconium alloy from-150 to 350 °C.The model is the simplest among the aforementioned models and can facilitate practical use.However,the predictability of the model in microscale has not yet been assessed.

From the perspective of micromechanical responses of magnesium alloys at elevated temperatures,the majority of the aforementioned models have been tested in terms of the texture evolution during deformation or they have not been evaluated at such small scales.However,more comprehensive microstructural characteristics can be represented by lattice strain,which has been widely used for various metals[22,23,29,36,37].This can offer a much stronger linkage between the microscopic and macroscopic behaviors of metals.The evolution of lattice strains during deformation is normally measured usingin situneutron diffraction or high-energy Xray diffraction(HEXRD).A recent publication by the authors[29,38]stated that thein situHEXRD enables continuous diffraction measurement during loading without interrupting the test,whereas neutron diffraction tests are encumbered by the stress relaxation effect owing to its long collection time.Moreover,the volumetric measurement of the HEXRD with its typical energy,which is orders of magnitude higher than conventional X-ray diffraction with Cu–Kα[38,39],provides much more accurate statistical measurements than the surface measurements by conventional X-rays.In situHEXRD has been utilized to investigate the micromechanical behaviors of various materials,such as microscopic lattice strain[29,38–43].Furthermore,several of the lattice strain comparisons betweenin situHEXRD measurements and crystal plasticity model predictions are available[29,38].Recently,in situHEXRD techniques equipped with heating devices have been used to unveil the temperature-dependent micromechanical behaviors of structural alloys[40–42].However,sophisticated techniques have not been linked with constitutive models such as the crystal plasticity model,and a comparison of the HEXRD-measured and model-predicted lattice strains at elevated temperatures has not been reported in other studies.

The current work focuses on the development of a timeefficient temperature-dependent modeling scheme and the close correlation of the microscopic and macroscopic behaviors of HCP-structured magnesium alloys at elevated temperatures.The macroscopic stress–strain curves and microscopic lattice strains of a twin roll-cast ZEK100 alloy were measured by thein situHEXRD experiment during uniaxial tensions at three discrete temperatures,RT,100,and 200 °C.In modeling perspective,the temperature-dependent hardening modelproposed by Singh et al.[35]was implemented in elasticplastic self-consistent(EPSC)and crystal plasticity finite element(CPFE)models.The EPSC and CPFE models were then further evaluated at the macroscale by comparing stress–strain behaviors as well as at a smaller scale by comparing lattice strains.The CPFE enforces grain-to-grain equilibrium and compatibility in a finite element(FE)sense[44],thus explicitly treating the interactions more realistically.However,it is CPU-intensive for use with a large number of polycrystalline materials.To overcome this problem,the procedure outlined in previous works for a multi-phase material(QP steel)[38]and for the ZEK100 alloy in which multiple slip and twin systems with dissimilar critical resolved shear stress(CRSS)are activated,[29]was followed,that is,the CRSSs and hardening parameters of available slip/twin systems at RT were first identified using a computationally efficient EPSC model using the trial-and-error approach.The identified parameters were thereafter informed of the CPFE model and refined.This approach could save computing time and cost significantly with reasonable accuracy.The scaling factor,which was proposed by Singh et al.[35],representing the temperature dependency of non-basal slip systems,was subsequently identified by the EPSC simulations.The stress–strain curves and lattice strains at the three testing temperatures were predicted by the EPSC and CPFE using the identified scaling factor.Finally,the developed temperature-dependent crystal plasticity models were validated on multiple scales.This coupledin situHEXRD and crystal plasticity modeling approach enabled a close tie between the microscopic and macroscopic behaviors of the HCP structured ZEK100 at elevated temperatures.

The remainder of this paper is organized as follows:In Section 2,the details of thein situHEXRD experiments at various temperatures are presented.In Section 3,the theoretical background of the temperature-dependent EPSC and CPFE models for HCP metals are provided.The EPSC/CPFE predicted and HEXRD measured results are compared in both micro-and macro-perspectives and discussed in Section 4.

2.Experimental details and results

A twin roll-cast ZEK100 magnesium alloy with a thickness of 1.5 mm was investigated in this study.The chemical composition of ZEK100 is summarized in Table 1.

Table 1Chemical composition of ZEK100 in wt.%.

Table 2Experimental matrix used for parameter identification and for model validation.

2.1.In situ HEXRD tensile test at elevated temperature

Thein situHEXRD tensile test setup performed at beamline 11-1D-C of the Advanced Photon Source(APS)at Argonne National Laboratory is depicted in Fig.1(a).The distance from the sample to the area detector D was set to 1.816 m,calibrated via diffraction of the standard CeO2sample.Sub-sized tensile samples with a gauge length of 18 mm made of the 1.5 mm thick ZEK100 alloy,depicted in Fig.1(b),were used.A monochromatic synchrotron Xray beam with a wavelength of 0.01173 nm radiated on the sample mounted on a custom-built tensile frame with a loading capacity of 13 kN.The high-energy X-ray beam penetrates through the entire thickness of the ZEK100 alloy sample;therefore,the entire transmission geometry was used for the investigation.For the high-temperature testing at 100 and 200 °C,the sample temperature was increased by a spot-type(in the order of a few millimeters)heat gun in the experiment.

The incident beam diffracts during the penetration of polycrystalline materials,in accordance with the following Bragg’s law:

wherehkildenotes the Miller–Bravais indices,which will be used throughout this paper,of the lattice plane.Termsdhkilandθhkildenote the lattice spacing and diffraction angle of the(hkil)planes,respectively.

To measure the initial texture as an input for the crystal plasticity models,that is,EPSC and CPFE,the sample was rotated near the loading direction(LD)from 0 to 90 °with an interval of 4 °.The initial texture was measured at RT.From the HEXRD measurements,a two-dimensional(2D)Debye ring image can be obtained.A typical example of the image is shown in Fig.1(c),along with the corresponding indices of the lattice planes.The obtained 2D Debye ring images were thereafter analyzed using the fit2D software[45].From the set of measured pole figures of the{100},{0002},and{110}lattice planes,the orientation distribution function(ODF)of the undeformed ZEK100 was estimated using the crystallographic texture analysis package MTEX[46].The coefficients of the ODF were identified by minimizing the error between the experimentally collected and MTEXcalculated pole figures using least squares fitting.In the current analysis,the de la Valle–Poussin kernel and a 5 ° halfwidth were used.The calculated basal{0002}and prismatic{100}pole figures with intensity ranging 0–3.5 and 0–3,respectively,are shown in Fig.2.Unlike typical AZ31 alloys,ZEK100 exhibits a weakened basal texture.Most of the grains are aligned such that their c-axes(normal of basal plane)are(nearly)perpendicular to the rolling direction(RD).In addition,the c-axes are spread along the plane of transverse direction(TD)–normal direction(ND).Two major strong basal poles were aligned at approximately±25 ° from the ND.The ZEK100 alloys with slight variation of the alloying contents reported in previous works exhibited similar textural characteristics[1,9,47].Note that intensity maxima of the basal poles of the investigated ZEK100 alloy are～3.5,which is approximately 2,3 times lower than those of the commercial AZ31 alloy exhibiting strong basal texture[12,47–50].

After the diffraction patterns of the undeformed sample were measured at various angles at RT,a tensile test was performed.The tests were conducted at constant grip speed of 30 μm/s,approximately equivalent to the nominal strainrate of 10-4/s.The tensile tests at the other temperatures,i.e.,100 and 200 °C,by the application of a heat gun were also performed under the same testing conditions.Continuous HEXRD measurements were viable with sufficiently low imposed strain rates without interruption(or stopping)of the tensile test.To investigate the anisotropic behavior of ZEK100,tests were conducted along two directions,RD and TD1It should be noteworthy that tensile tests along TD were conducted only at RT and 200 °C,while the test at 100 °C was not conducted..During the test,the tensile load was recorded from a load cell attached to the tensile frame,and the displacement of the tensile frame was also collected.The engineering stress and strain can be calculated using the load and displacement data.Notably,the stiffness of the tensile frame should be insufficient.The procedure of stiffness correction introduced by Hu et al.[36]was employed as follows:

Fig.1.(a)Schematic diagram of in situ HEXRD tensile test setup,where Shkil denotes the diffraction vector[38](d:magnitude of a vector from the center to the circumference of the Debye ring;D:Distance from sample to the detector),(b)sub-sized uniaxial tensile sample geometry,and(c)a Debye ring image recorded before the loading of the ZEK100 alloy.

Fig.2.Calculated basal and prismatic plane pole figures of undeformed ZEK100.

where superscript‘u’represents‘uncorrected’,ethe engineering strain,Sthe engineering stress,andEthe Young’s modulus(45 GPa for magnesium).Therefore,eudenotes the uncorrected engineering strain,whereasEudenotes the initial slope of the uncorrected stress–strain curve.It is worthy to note that the corrected stress–strain curve calculated using the grip displacement and load data for a high strength steel sheet,1 mm thick QP980 steel sheet,correlated very well with the results obtained from strain measurement by digital image correlation(DIC)[38,51].The stiffness-corrected strain and stress were thereafter converted to the true stress and strain.

2.2.Lattice strain calculation

When the sample is under tensile deformation,the lattice of grains gradually deforms such that it elastically stretches along the LD and shrinks along the transverse loading direction(TLD).Lattice strains are the elastic strains of the lattice plane families of grains with similar crystallographic orientations.The extent of‘similarity’in terms of crystallo-graphic orientation was determined by a range of scattering vector(Shkilin Fig.1),here,set as±2 °.In other words,to calculate the lattice strains along the LD,the diffraction peak intensities from the measured Debye ring images were analyzed using fit2D software to integrate along the arc withη=88–92 °(here,the angleηwas defined as the angle from the TLD along the counterclockwise direction,as depicted in Fig.1(c)),denoted by I88–92.I88–92was thereafter calculated as a function of 2θhkilat discrete macroscopic strains,as shown in Fig.3.The I88–92vs.2θhkilplot was fitted by OriginPro 2017[52]using a Gaussian distribution function to determine the diffraction angle of the intensity peak.Using the determined diffraction angle at the intensity peak at each macroscopic strain,that is,2θεhkil,the lattice strains as a function of macroscopic strain can be calculated using the following equation:

Fig.3.I88–92 vs.2θ plots of an undeformed sample(ε=0)and a deformed sample with a tensile strain of 0.1(ε=0.1)along the RD at RT for the{100}lattice plane peak[29].

where superscriptεdenotes the(macroscopic)strain along LD.

The lattice strains along the TLD could also be calculated following the same procedure described above with a tolerance range ofη=-2,2°.Finally,the lattice strains along the two measurement directions,LD and TLD,and under the two tensile directions,RD and TD,could be obtained.

3.Modeling

The EPSC is a standalone code and is much more computationally efficient than the CPFE;EPSC was 2,3 orders faster than the CPFE for solving problems in the current work.However,the CPFE explicitly treats grain-to-grain and phase-to-phase interactions more realistically.Furthermore,the CPFE is capable of solving complicated internal and/or external boundary problems with the advantage of the finite element method.A strategy to overcome the computational inefficiency of the CPFE was presented in author’s previous works for a QP steel[38]and a magnesium alloy[29].In the approach,the constitutive parameters were first estimated by the EPSC model.The estimated parameters were fine-tuned later via the CPFE.The same approach was employed in the current study.

The following two sections,that is,Sections 3.1 and 3.2,provide the essential perspective of the EPSC and CPFE models incorporating the deformation twinning kinematics.The interested reader is directed to the detailed paper of the EPSC model[53]and CPFE model[29,54–57].

3.1.Elastic-plastic self-consistent(EPSC)model

Plasticity models for polycrystals provide a coupling mechanism of the single-crystal constitutive behavior with that of the overall aggregate.In the EPSC model,an extended Voce hardening model was used to express the strain-induced CRSS evolution of theα-th slip orβ-th twin systems2Throughout the manuscript,the slip and twin systems are distinguished by superscripts α and β.,τc,as follows:

wherehijdenotes a hardening coefficient represented as follows:

whereqdenotes the ratio of the‘latent’hardening rate to the‘self’hardening rate,andδijdenotes the Kronecker delta.As‘self’and‘latent’hardenings could not be distinguished in the diffraction data,a value of 1.0 was used for all components ofhijof all slip/twin systems considered in the present study.Coupling the single-crystal constitutive behavior to the response of the aggregate is accomplished using a selfconsistent algorithm based on the Eshelby ellipsoidal inclusion formalism[58].Each crystallographic orientation(or grain)is treated as an ellipsoidal elasto-plastic inclusion with instantaneous moduliMgembedded in a homogeneous effective medium(HEM)with overall moduliM,which reproduces the overall response of the polycrystal.To reproduce the overall texture,the set of discretized Euler angles was weighted,which specifies the anisotropic properties of the polycrystalline.The interaction equation linking the microscopic strain rate tensor at the grain scaleand the macroscopic strain rate tensor of the polycrystallineis derived as follows:

withAgexpressed as follows:

whereM?denotes a constraint tensor for a matrix containing an ellipsoidal inclusion with the same orientation and shape as the grain.At each time step,the constitutive equation for the single crystal and self-consistent criteria are solved simultaneously.Thus,the grain-level mechanical responses were consistent with the boundary condition imposed on the HEM.

Based on experimental evidence from neutron diffraction studies[59,60]and other deformation twinning modeling[2,18,20–23,37,61],deformation twinning was assumed to be pseudo slip,which relies on the Schmid law.

The twin volume fractionfβin a grain evolves as follows:

whereγTWdenotes the twinning shear determined by the aspect ratio of the crystal lattice,[62]:

Sinceχ=1.624 for magnesium,the twinning shears can be calculated as 0.129(for tensile twin),and 0.138(for compressive twin).Only the positive resolved shear stress is allowed to account for the polar nature of the twin[61,63].The compressive twin is known to be rarely active in magnesium alloys because of its high CRSS[62],and it is not considered further.

An enhanced predominant twinning reorientation model(ePTR)proposed by Bong et al.[29]was employed to represent the grain re-orientation induced by deformation twinning.An abrupt grain re-orientation is introduced when the accumulated twin volume fraction exceeds the threshold valuefcr,which ranges from 0 to 1,at an integration point.According to[29],fcris expressed as follows:

where,P1,P2,andP3are material constants.Note thatP3is the newly introduced exponent term in the ePTR model.TermsFEandFRdenote the volume fraction of fully twinned grains and the volume fraction of the twinned region over all grains,respectively.A single twin system with the highest shear strain increment was allowed to activate.In the ePTR model,P1andP2control the twinning rate in the early and later stages of deformation,respectively.The exponent termP3offers flexibility of the twinning in the middle of the deformation,and the experimentally observed twinning rate was much better reproduced with the newly introduced term for the ZEK100 alloy at RT[29].

TermsFEandFRcan be mathematically expressed as follows[29,64]:

whereNGdenotes the total number of grains assumed in the simulations,NGtwinnedis the total number of fully twinned grains whose twin volume fraction exceedsfcraccording to the ePTR model,andwdenotes the volume fraction of the corresponding grain over all existing grains.

Once the twin volume fraction on all the twin systems(fβ)reachesfcrthe entire grain is completely reoriented by a transformation matrix,Q,between the lattice orientation in the matrix and that in the twinned region as follows:

whereIdenotes the identity tensor,andnrepresents the twin plane normal vector of the twin system,which has the highest contribution to the total twin volume fraction.

3.2.Crystal plasticity finite element(CPFE)model

A classical isothermal crystal plasticity constitutive equation proposed by[54,55]was employed.The total deformation gradient is decomposed into Feassociated with elastic distortion and rigid body rotation,whereas Fpis associated with the deformation gradient representing plastic deformation caused by dislocation motions[65]as follows:

Anisotropic elastic-viscoplastic constitutive equation for the single crystal is represented as follows:

where S denotes the second Piola–Kirchhoff stress,Cethe fourth-order anisotropic elasticity stiffness matrix,E=(1/2)(FeTFe-I)the elastic Lagrangian strain tensor,andσthe Cauchy stress tensor.

The rate of Fpin Eq.(14)is expressed as follows:

The velocity gradient Lpis associated with the dislocation glide on the active slip planeαor on the active twin systemβwith its normaland direction.Using the approach introduced by Kalidindi[63],the velocity gradientLpcan be represented by the sum of shear strain rates,˙γ(αorβ),over the total number of available slip and twin systems as follows:

where the first and second terms in the right-hand side are associated with the slip and twin in the untwinned region,respectively.It is worth mentioning that this does not necessarily mean that twinned grain,in which the twin volume fraction exceeds a critical value,fcr,does not deform by slip.There can be untwinned region in the fully-twinned grain(fβ≥fcr)as thefcrhas the value between 0 and 1 according to the ePTR scheme.Therefore,the slip can be active in the fully-twinned grain with new crystallographic orientation.

The shear rates on theαth slip systems,˙γα,andβth twin system,,are expressed as follows[51,53]:

The extended Voce hardening model and its evolution adopted in the EPSC model,Eqs.(4)–(6)),were employed in the CPFE model to represent the strain-induced evolution of CRSS in single crystal.For the CPFE modeling,same twinning kinematics through Eqs.(9)–(13)was employed.The Cauchy stress tensor in the grain can be calculated as the volumetric average of stresses in the matrix and twinned regions as follows:

whereσMandσβdenote the Cauchy stress tensors in the matrix and twinned regions,respectively.

For a fair comparison between the EPSC and CPFE models,the extended Voce hardening law used for the EPSC model in Eqs.(4)and(5)was also adopted in the CPFE model.The‘self’and‘latent’hardenings are also treated following Eq.(6).The model described above is implemented in ABAQUS/standard through a user-defined material subroutine(UMAT)following the procedure outlined in previous studies[54].

3.3.Temperature-dependent modeling

The aforementioned formulations complete the constitutive models,i.e.,EPSC and CPFE models,under isothermal conditions.This section provides modeling scheme relevant to temperature-dependency employed in both EPSC and CPFE models.To systematically model the temperature-dependent mechanical behavior of prismatic and pyramidal slip systems,the instantaneous CRSS is formulated as a function of temperature as follows[35]:

where the Voce hardening parameters of the prismatic and pyramidal slip systems are expressed as a function of temperature,that is,τ0(T),τ1(T),θ0(T),andθ1(T).

These parameters are scaled by a scaling factor,k(T),as follows:

wherek(T)denotes a newly introduced scaling factor to model the temperature dependency of the non-basal slip systems,whereasT0denotes the reference temperature set to RT(or 25 °C).

The scaling factork(T)ranges from 0 to 1,and it is expressed as follows[14]:

wherek1andk2denote the material constants to be identified.The specific form of Eq.(22)is proposed based on reported values[31,66–68].In the corresponding references,experimentally measured and viscoplastic self-consistent(VPSC)simulated CRSSs of non-basal slip systems at various temperatures were reported.On the experimental side,Yoshinaga and Horiuchi[68]and Wonsiewics and Backofen[67]measured the CRSSs of prismatic slip,and Obara et al.[66]measured the CRSSs of pyramidal＜c+a＞slip of magnesium single crystals.On the prediction side,Jain and Agnew[31]predicted the CRSSs of non-basal slip systems using VPSC simulations.The CRSSs,either measured or predicted,were normalized by the maximum CRSS at the lowest testing temperature.Note that,the normalized CRSS corresponds to the scaling factork(T)in Eq.(21).They were thereafter fitted using Eq.(22).The results are shown in Fig.4.The best-fit curvereproduces well the thermal softening behavior of non-basal slip systems;therefore,Eq.(22)was adopted to represent the temperature dependency of non-basal slip systems.

Fig.4.Scaling parameter k(T)of non-basal slip systems,i.e.,prismatic and pyramidal＜c+a＞,at various temperatures and their fit using Eq.(22).

It is also noteworthy that the basal slip and tensile twin exhibit athermal behavior in the temperature range of interest[2,31,32,69];thus,the two systems were essentially modeled to be temperature-independent.

The elastic stiffness constants were assumed as C11=58,C12=25,C13=20.8,C33=61.2,C44=16.6 GPa[70]owing to their weak temperature dependency.The strain rate sensitivity of magnesium alloys is known to be markedly dependent on temperature.The experimentally measured data of a ZEK100 alloy at different temperatures ranging from RT to 150 °C[1]was employed to identify the strain rate sensitivity exponentmformulated according to the following equation[16]:

wherem1andm2denote material constants,and they are determined as 0.02 and 2.147×10-4,respectively,from the reported data.

3.4.Simulation setup

Following previous modeling incorporating the crystal plasticity models for magnesium alloys[20–23],four distinct available slip and twin systems were considered in the EPSC and CPFE models for the ZEK100 alloys:(a)basal slip{0001}＜110＞,(b)prismatic slip{101}＜110＞,(c)pyramidal＜c+a＞slip{112}＜113＞,and(d)tensile twins{102}＜011＞.Using the calculated ODF described in Section 2.1,the discrete Euler angles in the Bunge definition were extracted.A total of 24,000 discrete Euler angles were used to represent the polycrystal for both the EPSC and CPFE models.

Fig.5.FE model for the tension along RD with 24,000 grains.

The EPSC model was simulated using incremental stress boundary conditions.In each time step,the stress and strain of each grain and the average response of the HEM were calculated when the self-consistency condition was met.The measured stress–strain curves and lattice strains by thein situHEXRD tensile experiments were used as objective values to identify the CRSS and hardening parameters of the available slip and twin systems,along with ePTR parameters at RT using the trial-and-error approach.

The CPFE simulations,incorporated with the commercial software ABAQUS/standard 6.12,were performed with the developed user subroutine.An FE model for the tension along the RD shown in Fig.5 was developed with a similar length ratio along theX,Y,andZdirections in the actual tensile sample,that is,12:2:1.The three-dimensional array is discretized by 120×20×10(X×Y×Z)FEs,that is,24,000 in total,using C3D8R(eight-node solid element with reduced integration)type mesh.Each FE represents a grain with a specific crystallographic orientation;thus,the FE model represents a polycrystal with 24,000 grains3The use of the reduced 3D element,C3D8R,did not cause hourglass problem over all CPFE simulations without locally concentrated deformation in solving the uniaxial tension problems..A prescribed displacementboundary condition with a macroscopic strain rate of approximately 10-4/s was imposed on the nodes on the positive X plane,as depicted in Fig.5 whereas the nodes on the X-plane on the other side were fixed in all directions.No constraints along theY-andZ-directions were applied.Therefore,the FE model represents an arbitrary location of the gauge area of the sub-sized actual tensile sample.The foregoing EPSC and CPFE simulation procedures were repeated at three different testing temperatures using the determined constitutive parameters at the corresponding temperature.For the tension simulation along the TD,theX-andY-directions were set parallel to the TD and RD,respectively.

In both the EPSC and CPFE models,the macroscopic stress–strain response of the polycrystalline was calculated using the corresponding values of individual grains in the averaging sense as follows:

whereNGdenotes the total number of grains in the polycrystalline aggregate,andσnandεndenote the stress and strain tensors of thenthgrain,respectively.

The lattice strains were calculated using the elastic strain tensor of each grain from the EPSC and CPFE predictions.For the lattice strain along the LD,the principal elastic strain components along the LD of grains whose lattice plane normal is almost parallel to the LD(within a tolerance ofη=88–92 °)were averaged.The lattice strain along the TLD was also calculated in the same manner with a tolerance ofη=-2–2 °.The calculated average lattice strains were compared with insituHEXRD data.

4.Results and discussion

The experimental data used for fitting in constitutive parameter identification and for validating the developed model are summarized in Table 2.The temperature-dependent constitutive behaviors of ZEK100 were characterized as follows:(1)The constitutive parameters of available slip and twin systems at the reference temperature,that is,RT,were first identified using the EPSC model by fitting to both experimental stress–strain curves and lattice strains;(2)The EPSC-identified constitutive parameters were refined for the CPFE model;(3)The temperature-dependent scaling factor,introduced in Eq.(22)was determined using the measured stress–strain data at three testing temperatures:RT,100,and 200 °C;(4)The EPSC and CPFE simulations at each temperature assuming isothermal condition were followed using the identified constitutive parameters,and the prediction results were compared with the corresponding experimental data.In particular,the predictive capability of the developed model was assessed by comparing with the measured data which have not been used in the fitting process,i.e.,stress–strain curve during TD tension at 200 °C,and lattice strain evolutions at elevated temperatures,100,and 200 °C.

Section 4 is organized as follows:The details of the aforementioned procedures(1)and(2)are provided in the Section 4.1.In Section 4.2,the details in the process of the identification of the temperature-dependent scaling factor are presented.In Section 4.3,developed model is validated by comparing with independent experimental data from the fitting.

4.1.Parameter identification at reference temperature

The material parameters for the extended Voce hardening model(τ0(T0),τ1(T0),θ0(T0),andθ1(T0))were identified by the trial-and-error approach using the computationally efficient EPSC model.The four Voce hardening-related parameters assigned for each deformation mode constitute the fitting parameters;16 parameters in total,were fitted to the macroscopic stress–strain responses and the lattice strains measured by thein situHEXRD tensile experiments up to measured uniform elongation(i.e.,8.4 and 20.5% for the tension along the RD and TD,respectively).In addition to the 16 hardening parameters,the ePTR model parameters(P1,P2,andP3)related to the tensile twinning kinematics in Eq.(11)were simultaneously identified using the basal plane{0002}pole intensity data measured byin situHEXRD during the tensionalong the TD at RT.Note that the intensity variation in the{0002}pole is associated with the twinning rate,as discussed in[20],and the determined constitutive parameters for the extended Voce hardening and ePTR models at RT are listed in Table 3.

Table 3Best-fit constitutive parameters of slip/twin systems at RT for ZEK100 via EPSC model.

Table 4Best-fit constitutive parameters of slip/twin systems at RT for ZEK100 via CPFE model.

Fig.6.Comparisons between measurements from in situ HEXRD tensile tests and EPSC/CPFE predictions at RT:(a)macroscopic stress–strain curves during tension along the RD and TD at RT,(b)lattice strain as a function of macroscopic strain at RT during tension along the RD,(c)lattice strain as a function of macroscopic strain at RT during tension along the TD,and(d)normalized intensity of grains,whose basal plane normal is parallel to the LD within±2°scatter during tension along the TD.

The same set of constitutive parameters determined for the EPSC model listed in Table 3 were thereafter employed in the CPFE simulations as the first estimate,and they were further refined using the CPFE model.As discussed in previous research works by the authors[23],CPFE generally predicts faster strain hardening in two regimes compared with the EPSC:(a)regime of early deformation where elastic-plastic transition occurs,and(b)large strain regime.These two regimes are relevant to the initial and linear hardening parameters,θ0andθ1in Eq.(4).To mitigate the faster strain hardening in the CPFE,two parameters were simply scaled in sucha way that=0.5×and=0.7×.The ePTR parameters were also slightly refined.The refined parameters via the CPFE model are listed in Table 4,where the CPFE prediction refers to the prediction results using the refined constitutive parameters).

Fig.7.Comparison of stress–strain behaviors at three different testing temperatures:measurement from the tensile tests vs.predictions by the EPSC simulations.

The predicted stress–strain curves by the EPSC using the constitutive parameters at RT in Table 3 and by the CPFE using the refined constitutive parameters at RT in Table 4 are compared with the measured tensile test results in Fig.6(a).The measured stress–strain curves revealed anisotropic stress–strain behaviors.Because twinning governs the plastic deformation under tension along the TD,the yield stress under such loading conditions is approximately 80 MPa lower than that under tension along the RD.However,the non-basal slip systems,which are harder than the basal slip or twin,are more likely to be active in the twinned region as the deformation proceeds.Therefore,strain hardening is faster during the tension along the TD than during the tension along the RD.The CPFE and EPSC predictions are almost identical,and they are consistent with the measured data.In particular,both models accurately predict the anisotropic stress–strain behaviors along two different LDs.

The EPSC-and CPFE-predicted lattice strains of various lattice planes as a function of macroscopic strain during the tension along the RD and TD are shown in Fig.6(b)and(c),respectively,along with thein situHEXRD measured data.The EPSC and CPFE models reproduce the lattice strain evolution during tension along the RD,as shown in Fig.6(b).The lattice strain evolutions during the tension along the TD shown in Fig.6(c)also show that the EPSC and CPFE models reproduce the measured data well with marginal errors.Experimentally,the lattice strain of the{0002}plane,whose normal is parallel to the TLD or RD,abruptly becomes zero at a(macro)strain of～0.13.During tensile loading along the TD,the grains whose c-axes are nearly parallel to the TD are gradually twinned.Finally,no more grains are left within the angular tolerance(±2 °)at a certain instance;consequently,the lattice strain becomes zero.Although the EPSC and CPFE predict early drops in the lattice strain at～8%(macro)strain,two models capture this microscopic behavior.

Fig.8.Variation in k(T)as a function of temperature at RT,100 and 200 °C under the tension along the RD and its fit using Eq.(22).

The basal pole intensity variation,which gives insight into the rate of twinning,were also predicted and compared with those measured over the deformation history in Fig.6(d).Again,the basal plane peak intensity variation provides useful information on the deformation twinning rate.Under the tensile stress along the TD,the grains whose c-axis is within the angular tolerance(η=88–92 °)with the tensile LD deforms primarily by tensile twinning,and they are eventually reoriented.Therefore,the intensity gradually decreased and finally became zero.In previous works by the authors[23],efforts were devoted to accurately reproduce thein situHEXRDmeasured intensity change using the crystal plasticity models incorporating the PTR model,which was proposed to offer flexibility by introducing the exponent termP3.The comparison of the basal pole intensity in Fig.6(d)shows that the EPSC and CPFE incorporated with the ePTR model reproduce the measured data well.In terms of prediction accuracy,the CPFE considering the grain-to-grain interaction and spatial distribution of the individual grain can slightly better predict the measured data compared to the EPSC.It is worth note that the stepwise variation of the EPSC/CPFE predicted{0002}intensity is anticipated to the less number of grains assumed in the simulations than the experiments.Among 24,000 grains assumed in the simulations,only 20 grains are oriented such that their c-axis is within the tolerance range ofη=88–92 °.In contrast,more than 80 grains are oriented in such condition in thein situHEXRD experiment.

4.2.Temperature dependent constitutive parameter identification

Once the constitutive parameters for the reference temperature are identified,the temperature-dependent scaling factor in Eq.(21)along with its explicit form in Eq.(22)are straightforward.To model the thermally activated prismatic and pyra-midal slip systems,the scaling factors at elevated temperatures(100 and 200 °C)under tension along the RD were identified using a computationally efficient EPSC model.As mentioned previously,the basal slip and tensile twins were assumed to be temperature-independent.The scaling factor was set to 1 for the EPSC modeling at the reference temperature(RT),whereas the scaling factors were adjusted for the other testing temperatures to best estimate the measured stress–strain behaviors of the ZEK100 polycrystal.The fitting yields the scaling factor of 0.78 and 0.59 for the testing temperature of 100 °C and 200 °C,respectively.The fitted stress–strain curves are shown in Fig.7.Although,discrepancies can be found in terms of the yield stresses,the stress–strains behaviors are well reproduced by the EPSC simulations incorporating the simple temperature dependent scaling factor.

Fig.9.Comparison of stress–strain curves at various temperature:EPSC predictions vs.experiments(a)during tension along RD,and(b)during tension along TD,and CPFE predictions vs.experiments(c)during tension along RD,and(d)during tension along TD.

Using the scaling factors at the three testing temperatures,the material constantsk1andk2in Eq.(22)can be obtained.The two constants were determined by fitting the scaling factors as a function of temperature using a least-squares method,as shown in Fig.8.The fitting yieldsk1=246.36 andk2=1.06.

The proposed approach can save substantial computing time and cost for temperature-dependent modeling.Notably,the other previous approaches require individual fitting to each temperature using the trial-and-error approach[31,32]or complex fitting procedure to identify a large number of model parameters[18,33].Therefore,if it is successful in capturing the temperature dependence,the model proposed in this study can be meaningful,particularly in terms of practical applications.Hereafter,the proposed modeling approach will be comprehensively examined in two main streams:(a)macroscopic stress–strain behaviors and(b)lattice strains in the microscale at elevated temperatures along the two LDs.Note that the samek1andk2determined by the EPSC model will be adopted in the CPFE model without further refinement.

Fig.10.EPSC/CPFE-predicted and in situ HEXRD-measured lattice strains at 100 °C under the tension along the RD.

4.3.Validation

4.3.1.Macroscopic stress–strain curves

The stress–strain behaviors predicted by the EPSC and CPFE models along the two LDs at the three testing temperatures utilizing the identifiedk(T)are compared with the experimental data in Fig.9.Experimental data were plotted to measure the uniform elongation at each testing condition.Notably,the uniform elongation diminished as the temperature increased,regardless of the LD.The root cause of this phenomenon is the non-uniform temperature distribution generated by the spot-type heat gun used to increase the sample temperature during the experiment.Consequently,when the heat gun was used in the experiment,samples exhibited premature necking initiated from where the heat gun impinged,and thus,the commonly observed improved ductility of the magnesium alloys at elevated temperatures[71–73]was not observed in thein situHEXRD tension experiments.Nevertheless,the prediction results exhibit reasonably good consistency when compared with the predictions and experimental data.The yield stress drops at elevated temperatures in the RD tension and relatively smaller yield stress drops in the TD tension at elevated temperatures in the early state of the deformation are well captured.The hardening rate decreased slightly during tension along the RD and significantly decreased during tension along the TD at higher temperatures.Regardless of the LD,the EPSC/CPFE models captured the reduced hardening rates well.The results validate the fidelity of the simple temperature-dependency modeling scheme using Eqs.(20)–(22).

4.3.2.Lattice strain

The variations in the lattice strains for several major lattice planes were measured byin situHEXRD during deformation at elevated temperatures,that is,100 and 200 °C,and the two LDs were compared with the EPSC and CPFE predictions.The comparison of the lattice strains up to uniform elongations in the HEXRD experiments during tension along the RD and TD are shown in Figs.10 and 11,respectively.Overall,the results reveal that the EPSC and CPFE reasonably well captured the general trends of the lattice strain variation at the elevated temperatures observed from the HEXRD experiments with marginal discrepancies.

Fig.11.EPSC/CPFE-predicted and in situ HEXRD-measured lattice strains at 200 °C under the tension along(a)the RD(b)the TD.

The EPSC and CPFE underpredict the lattice strain evolution of the{0002}plane whose normal is(nearly)parallel to the TLD during the tension along the TD at 200 °C,as shown in Fig.11(b).This error might be attributed to 1the assumption of temperature independence with regard to tensile twinning.However,as reported by Jain and Agnew[31],the CRSSs of tensile twinning can be slightly temperaturedependent.Although the slight temperature dependency of the tensile twin has not been studied deeply to date,better prediction can be achieved if the constraints,that is,temperature independence of the tensile twinning,are relaxed.However,further simulations with refined parameters were not performed to maintain the simplicity of the developed models.

Fig.12.Relative activities and twin volume fraction calculated from EPSC predictions under the tension along(a)the RD,and(b)the TD,and from CPFE predictions under the tension along(c)the RD and(d)the TD.

Along with the prediction of the stress–strain curves,the agreement between the measured and EPSC/CPFE predicted lattice strains firmly ties the microscopic and macroscopic behaviors of the HCP-structured ZEK100 at elevated temperatures.Finally,the results presented in this section validate the precision of the developed temperature-dependent crystal plasticity models at multiple scales,as well as the practicality of the models.

4.5.Micro-mechanism of deformation

The crystal plasticity frameworks offer information on the activation of deformation modes,that is,slip and twinning,under a certain deformation.The relative activity of theκth deformation mode,Aκ,is calculated from the EPSC/CPFEsimulated results as follows:

whereNmodedenotes the total number of slip or twin systems belonging to each deformation mode corresponding to basal,prismatic,and pyramidal＜c+a＞slip systems,along with tensile twins.The evolution of the twin volume fraction can be obtained from the crystal plasticity frameworks besides the relative activity.The volume fraction of the twinned region over all grains,FRin Eq.(12),was assumed as the total twin volume fraction.

The calculated relative activities and total twin volume fraction from the simulated results by the EPSC and CPFE as a function of the macroscopic strain are shown in Fig.12.The relative activity and total twin volume fraction plots at RT and 200 °C are presented,whereas those at 100 °C are omitted because the relative activities and total twin volume fractions of the slip/twin systems at 100 °C lie between the two extremes,that is,RT and 200 °C,without abnormal behavior.The relative activity and total twin volume fraction plots predicted by the EPSC during tension along the RD and TD are shown in Fig.12(a)and(b),respectively.Most of the grains in ZEK100 are oriented such that their c-axis is parallel to the TD-ND plane(see Fig.2).These grains will be compressed along their c-axis during the tension along the RD.The compressive strain is mostly accommodated by the activation of pyramidal slip,and thus,the tensile twin is rarely active in such loading condition.In contrast,during tension along the TD,tensile twinning is noticeably active to accommodate the tensile strain along the TD in the grains whose c-axes are(nearly)parallel to the TD.The facts are also revealed by larger total twin volume fraction during the tension along the TD:the total twin volume fraction atε=0.1 and at RT during the TD tension is～7 times larger than that during the RD tension.Once the grain is fully twinned according to the PTR scheme,the grain rotates～86° and its c-axis is aligned parallel to the RD.Therefore,the pyramidal slip becomes more active as the activation of the tensile twinning diminishes.

The major changes in terms of the relative activity and the total twin volume fraction at the elevated temperature,200 °C,are as follows:(a)regardless of the tensile direction,the non-basal slip systems become more active,and(b)the tensile twinning is less active,and thus,the total twin volume fraction is lower at elevated temperatures during the tension along the TD.The two major changes are anticipated to be the comparability of the CRSSs of the non-basal slip systems to that of tensile twinning at elevated temperatures.Similar trends are discernible from the plots of the calculated relative activities from the CPFE predictions,as shown in Fig.12(c)and(d).One major difference between the CPFE and EPSC is that the EPSC predicts inactive pyramidal slip up toε～0.01,whereas the CPFE does not.The inactivation most likely originates from the treatment of elastic and plastic behavior in the EPSC model,that is,the shear rate is zero if the resolved shear stressτis less than the CRSSτcin the EPSC model[74],as discussed in[29].

5.Conclusions

In the present study,a phenomenological temperaturedependent constitutive model for HCP metals introducing a scaling factor was developed and implemented in crystal plasticity codes,namely,EPSC and CPFE models.The constitutive parameters to represent the hardening of slip and twin systems at RT were first determined using the computationally efficient EPSC model by fitting to thein situHEXRDmeasured stress–strain curves and lattice strains.The CPFE simulations at RT were followed by refined hardening parameters from those identified using EPSC.Finally,temperaturedependent scaling factors were determined using the EPSC model.Two crystal plasticity models were tested for a magnesium alloy ZEK100,not only in predicting the macroscopic stress–strain behavior but also the microscopic lattice strains at three different testing temperatures,that is,RT,100 and 200 °C,under tension along two different LDs,that is,the RD and TD.The prediction results were compared with insituHEXRD measurements.Particularly,the comparison of thein situHEXRD-measured and the crystal plasticity modelpredicted lattice strains are presented for the first time.The following conclusions can be drawn from the corresponding studies:

·The temperature-dependent crystal plasticity models were able to reasonably predict both the macroscopic stress–strain behavior and the microscopic lattice strain under uniaxial tension along the RD and TD at the three tested temperatures.

·The prediction results and their comparisons with the measurements,which are independent from the fitting procedure for the constitutive parameter identification,validate the precision and practicality of the developed temperaturedependent model.

·The consistency of the lattice strain evolutions at elevated temperatures,which has not been reported in previous studies,could more closely tie the microscopic and macroscopic responses of the HCP-structured ZEK100 at elevated temperatures.

·The simulation results and data analyses revealed that the non-basal slip systems were significantly more active,whereas the tensile twins were less active at elevated temperatures.

Disclaimer

This paper was prepared as an account of work sponsored by an agency of the United States Government.Neither the United States Government nor any agency thereof,nor any of their employees,makes any warranty,express or implied,or assumes any legal liability or responsibility for the accuracy,completeness,or usefulness of any information,apparatus,product,or process disclosed,or represents that its use would not infringe privately owned rights.Reference herein to any specific commercial product,process,or service by trade name,trademark,manufacturer,or otherwise does not necessarily constitute or imply its endorsement,recommendation,or favoring by the United States Government or any agency thereof.The views and opinions of authors expressed herein do not necessarily state or reflect those of the United States Government or any agency thereof.

Declaration of Competing Interest

The authors declare that they have no conflict of interest.

Acknowledgements

H.J.Bong appreciates the supports by the Fundamental Research Program of the Korea Institute of Materials Science(KIMS,PNK7760).Oak Ridge National Laboratory is operated by UT-Battelle,LLC,for the U.S.DOE under contract DE-AC05-00OR22725.This research used resources of the Advanced Photon Source(APS),U.S.DOE Office of Science User Facility operated for the DOE Office of Science by Argonne National Laboratory under Contract No.DE-AC02-06CH11357.