Comprehensive Modelling of the Hysteresis Loops and Strain–Energy Density for Low-Cycle Fatigue-Life Predictions of the AZ31 Magnesium Alloy

Magnesium is one of the lightest metals for structural components. It has been used for producing various lightweight cast components, but the application of magnesium sheet plates is less widespread. There are two reasons for this: (i) its poor formability at ambient temperatures; and (ii) insufficient data on its durability, especially for dynamic loading. In this article, an innovative approach to predicting the fatigue life of the AZ31 magnesium alloy is presented. It is based on an energy approach that links the strain–energy density with the fatigue life. The core of the presented methodology is a comprehensive new model for tensile and compressive loading paths, which makes it possible to calculate the strain–energy density of closed hysteresis loops. The model is universal for arbitrary strain amplitudes. The material parameters are determined from several low-cycle fatigue tests. The presented approach was validated with examples of variable strain histories.


Introduction
To achieve a lightweight design that will fulfil the ever-growing ecological requirements and reduce fuel consumption, especially in the automotive and aerospace industries, new materials are being sought that demonstrate a good strength-to-weight ratio [1][2][3][4][5][6]. Magnesium alloys appear as a promising solution for these requirements as magnesium is one of the lightest metals, approximately 1.6 times lighter than aluminium alloys and about 4.5 times lighter than steel, while it has high strength, good machinability and good recyclability [4,5,[7][8][9]. Various magnesium alloys exist. They are mostly suitable for casting purposes, e.g., AZ63 or AZ91 [10], although magnesium alloys such as AZ31 are available as sheet plates of several thicknesses that are suitable for forging purposes [4,8]. The use of the latter as the casting is not as widespread, mostly for two reasons. First, as AZ31 is poorly formable at ambient temperatures, it always requires forming at raised temperatures, typically over 230 • C, which increases the production costs. Second, there is a lack of experimental data on dynamic loading, which makes the durability predictions of structural components manufactured from AZ31 either very challenging or not reliable [5,6,10].
Three deformation mechanisms interchange during the dynamic loading of AZ31, i.e., sliding, twinning and untwining [6], although twinning and untwining could be interpreted as a single mechanism. Sliding is also the dominant deformation mechanism in conventional metallic structural components such as steel and aluminium alloys, whereas twinning and untwining are more typical of magnesium and titanium alloys [4,6,[11][12][13][14]. They are a consequence of the hexagonal close-packed crystal structure, where the twinning mechanism requires lower activation energy than pyramidal and prismatic sliding [4,6,[13][14][15]. Therefore, the crystal structure will include an amount of twinned crystals before or with the activation of the sliding mechanism, especially in the compressive direction [10]. If the crystal structure has already been twinned under compressive loading and then the load is inverted into the tensile direction, the untwining mechanism occurs prior to the sliding mechanism in tension as the activation energy for the untwining is lower than the activation energy for the sliding. In contrast, sliding requires the smallest activation energy in the body-centred and face-centred cubic crystal structures typical for steel and aluminium alloys, regardless of the loading direction [6]. On the macroscopic scale, a symmetrical stress-strain response is obtained for the body-centred and the face-centred cubic crystal structures if the material is loaded with an alternating strain [16]. The interchange of the three mechanisms in the hexagonal close-packed crystal structure causes an asymmetrical shape of the stress-strain curve, typically also observed for AZ31 [6]. If such a response is to be simulated, the simulation model has to enable distinguishing between the mechanisms which occur either during the initial loading or during the tensile/compressive loading paths of the hysteresis loops.
The dynamic loading of metallic structures can lead to fatigue failure. Under severe loads, the material will fail due to low-cycle fatigue. The strain-life relation is usually used to describe the behaviour of the materials for the low-cycle fatigue. The strain-life relation provides the influence of the strain amplitude to the fatigue life of the material for the given dynamic ratio (e.g., typically R = −1). In addition to the strain amplitude, the stress amplitude and the mean stress drastically influence the fatigue life of the material, too. The prediction models for low-cycle fatigue life typically assume that the strain is the known variable, whereas the stress response has to be modelled. Especially for variable load histories, neglecting the influence of the stress amplitude or the mean stress, erroneous fatigue-life predictions might be observed. Damage parameters have therefore been introduced which consider the influence of the strain amplitude, the stress amplitude and the mean stress (e.g., Smith-Watson-Topper and Bergmann damage parameters [17]). Damage parameters however are based on the dissipated energy within a hysteresis loop and are therefore directly connected to the strain-energy density calculation. By obtaining the hysteresis loops in a given variable loading history, it is possible to calculate the damage contribution of every hysteresis loop to the total accumulated damage. Due to the rather complicated shape of the closed hysteresis loops for the AZ31 magnesium alloy, we present in this paper how an energy approach to calculating the fatigue life can be applied. In the literature, different energy approaches were applied to calculate the fatigue damage. They differ according to the definition of the critical damage. Some approaches link the fatigue-life damage to the static tensile curve (e.g., see Letcher et al. [18] and Ozaltun et al. [19]). Another approach is to link the energy to the Coffin-Manson durability curve (e.g., see Park and Nelson [20], Kabir and Yeo [21], and Jahed et al. [22]). We decided to follow a different approach in which the fatigue life is calculated according to the energy curve of the material, because this approach was already applied in the past to calculate the fatigue life for both the constant and variable strain amplitude loading of magnesium alloys (e.g., see Jahed et al. [23], Park et al. [24], Wang et al. [25] or Dallmeier et al. [26,27]).
As we focused on the fatigue-life prediction under a variable loading history, it is important to be able to properly calculate the energy density of a hysteresis loop considering its position in the stress-strain space. Sigmoidal functions are used to describe the S-shaped stress-strain response of magnesium alloys [27][28][29][30]. Some researchers have decomposed the closed hysteresis loop of magnesium alloys into its elastic and plastic parts (e.g., Lee et al. [28] and Muhammad et al. [29]). They put a special emphasis on the modelling of the plastic part of the hysteresis loop. Since their experimental data ranged up to a strain amplitude of 6%, they observed the twinning and sliding plasticity mechanisms in both the tensile and the compressive loading paths. Even though the models of Lee et al. [28] and Muhammad et al. [29] are somewhat different, their basic idea is to have three different models for the hardening law(s). Their models can be summarised as follows (see also Figure 1): 1. A monotonic tensile curve (0-1) is a sum of the constant monotonic tensile yield stress σ 0,T and the exponential curve with the generic form A · 1 − exp −B · ε p . 2. A compressive reverse loading path (1-2) is a combination of the constant monotonic compressive yield stress σ 0,RC , the exponential curve A · 1 − exp −B · ε p for modelling the twinning process and the sigmoidal function for modelling the sliding plasticity during the compression. 3. A tensile reverse loading path (2-3) is a combination of the constant monotonic tensile yield stress σ 0,RT , the exponential curve A · 1 − exp −B · ε p for modelling the untwining process and the sigmoidal function for modelling the sliding plasticity during the compression. The parameters for this model are different from the parameters of the second model. For the modelling of the closed hysteresis loop, only the last two models are important, because our objective is to calculate the fatigue-life damage due to cyclic loading. The main problem linked to the models of Lee et al. [28] and Muhammad et al. [29] is that the elastic and plastic parts of the hysteresis loop are separated, since the Masing rule cannot be applied for the magnesium alloy. Furthermore, we noticed that the initial plastic part of the tensile and the compressive loading path is not well modelled with the sum of the initial yield strength and the exponential function: It turned out that the plastic part of the Ramberg-Osgood equation [31] fits the experimental data much better. In addition, it can be applied to the elastic-plastic material curve relatively easily, which is well suited to our problem. Consequently, we used the full Ramberg-Osgood equation in the following form (see Dowling [32]) as a base function for the compressive and tensile loading paths: where E is the material's Young's modulus, K is a function of the yield strength and n is a hardening exponent. By adding a sigmoidal function, the transition from exhaustion of the twinning into to the sliding plastification mode during the tensile loading path can be described, as also noted by Dallmeier et al. [27]. However, the main advantage of the new model is its comprehensiveness. Namely, the cyclic behaviour of AZ31 can be correctly described parametrically and phenomenologically, regardless of the size and the position of the hysteresis loop in the stress-strain space, i.e., a single function with a single set of material parameters, without the need for interpolation between the hysteresis loops of different sizes, is required for the stress-strain modelling. The energy-density calculation is then possible using this new model for the parametric description of the compressive and tensile loading paths. For a variable loading history, a rule of closed hysteresis loops must also be considered [33,34]. Furthermore, special attention has been given to the large nested cycles that do not originate at the envelope of the maximum strain, including a variety of even smaller nested cycles that do not finish at the envelope.
Such cycles indicate the modelling capabilities of a material model. The paper is structured as follows. First, the theoretical background for the calculation of the strain-energy density of closed hysteresis loops is given in Section 2. The experimental results are given in Section 3. The parameters of the hysteresis loops for AZ31 are determined in Section 4. In the same section, a comparison between the obtained and predicted fatigue-life results for a variable loading history is given. Section 5 summarises the research outcomes.

Strain-Energy Density of a Closed Hysteresis Loop
Our model for the calculation of the strain-energy density of the closed hysteresis loops for the magnesium alloy is a combination of two different stress-strain functions: a model of the compressive loading path and a model of the tensile loading path. Since the maximum strain amplitudes were smaller than 2%, in our case, the sliding-plasticity mechanism was not activated during the compression phase of the loading cycle. For this reason, only the Ramberg-Osgood equation was applied for the compressive loading path of the stabilised magnesium hysteresis loops: where stabilisation represents a half of the fatigue lifetime. On the other hand, the sliding-plasticity mechanism was always activated after the untwining process during the tensile phase of the loading cycle. This is why it was decided on the basis of the literature (Dallmeier et al. [27], Lee et al. [28], and Muhammad et al. [29]) and our own research that a composite model for a tensile loading path of the stabilised Mg-hysteresis loop is used. It is composed of two terms: the Ramberg-Osgood relationship, which describes the untwining mode during the tensile reloading, and the sigmoidal relationship, which describes the transition from the exhaustion of the twinning into the sliding plastification mode: The term RO −1 T (ε) in Equation (4) represents an inverse function of the Ramberg-Osgood equation with the material parameters that are characteristic for the tensile loading path of the hysteresis loop: Since Equation (5) cannot be inverted, the stress σ as a function of strain ε should be calculated numerically. A Newton-Raphson scheme was applied for this purpose, because the function is general and is valid for the complete range of strains for the tensile loading paths. The two parameters B and F depend on the maximum strain ∆ε max,cyc , related to a starting point (SP) of the tensile loading path from the most outside compressive loading path (see the left-hand diagram in Figure 2): F(∆ε max,cyc ) = f 1 · (∆ε max,cyc ), ∆ε max,cyc < f 2 ; The parameter B describes the increasing and decreasing of the logistic function along the stress σ, while the parameter F defines a position of the sigmoidal function along the strain ε. The parameter f 2 represents a saturation strain for the untwining process within the tensile loading path. In this manner, the complete model for the stabilised, closed Mg hysteresis loop is composed of only ten material constants: E, K C , n C , K T , n T , b 1 , b 2 , D, f 1 and f 2 . The parameters can be estimated from the measured closed hysteresis loops at different strain amplitudes with numerical optimisation algorithms. Knowing the material constants, the plastic strain-energy density of a closed hysteresis loop is calculated as where σ is the modelled stress and ε is the given strain in this loop.

Calculating Fatigue-Life with an Energy Approach
According to the selected approach, the energy fatigue-life curve is determined first from the constant strain-amplitude fatigue-life experiments. If the effect of the mean stress is not significant, i.e., if the fatigue loading cycles have approximately equal dynamic coefficients R, the energy fatigue-life curve is defined as follows: where N f is the number of loading cycles to failure, ∆W p is the plastic strain-energy density of the stabilised hysteresis loops, and m p and C p are material constants. If the loading-cycle mean stress needs to be considered, the energy fatigue-life curve is related to the total strain-energy density ∆W t : The total strain-energy density is the sum of the plastic strain-energy density and the tensile elastic strain energy, as presented in the right-hand diagram in Figure 2 (see also Park et al. [24]): where σ max,cyc is the loading-cycle maximum tensile stress, E is the Young's modulus, and m t and Ct are material constants. To calculate the fatigue life for a variable strain time series, the closed hysteresis loops are first extracted from the time series (see Section 2.3 for details) with the rainflow counting method according to the ASTM E1049 standard [35] or Amzallag et al. [36]. Then, the fatigue life is calculated with the linear damage-accumulation rule (Palmgren-Miner) on the basis of the energy fatigue-life curve from Equation (9) or Equation (10): where Dmg is the fatigue-life damage and n i is the number of loading cycles at the loading level k. In our case, we repeated a block of variable-strain time series until the specimens were broken. In that case, the fatigue life is represented as the number of repeated variable-strain blocks until failure. If the critical fatigue-life damage is defined as Dmg c , the number of variable-strain block repetitions to failure rep is equal to: In Equation (15), k represents the number of closed hysteresis lops in one variable-strain block.

Extracting Closed Hysteresis Loops of Loading Cycles
The closed hysteresis loops for the variable ε(t) time series are modelled on the basis of the rainflow counting algorithm according to Amzallag et al. [36]. To avoid the formation of the residuum after the counting, the ε(t) time series is transformed so that it starts and ends with the maximum strain ε max (see Figure 3). With the exception of the largest and the outermost loading cycle, which is defined with Points 1, 4 and 7 in Figure 3, all other loading cycles that form a closed hysteresis loop in the σ-ε diagram must fulfil the condition: where ε i−2 to ε i+1 are four consecutive points in the strain time series. ε i corresponds to the peak of the standing loading cycle with ε i > ε i−1 or to the valley of the hanging loading cycle with ε i < ε i−1 (see Figure 3). The standing loading cycles are modelled so that first the tensile loading path σ T;i−1,i (ε) is modelled between the strain points ε i−1 and ε i with a reference to the preceding compressive loading path that is spanned between the strain points ε i−2 and ε i−1 . Then, the reversed compressive loading path σ C;i−1,i (ε) between the strain points ε i and ε i−1 is modelled so that the standing loading cycle is completely closed. Similarly, the hanging loading cycles are modelled so that first the compressive loading path σ C;i−1,i (ε) is modelled between the strain points ε i−1 and ε i with a reference to the preceding tensile loading path that is spanned between the strain points ε i−2 and ε i−1 . Then, the reversed tensile loading path σ T;i−1,i (ε) between the strain points ε i and ε i−1 is modelled so that the hanging loading cycle is completely closed. Figure 3 schematically explains this modelling.
If four consecutive strains ε i−2 to ε i+1 in the strain time series are found that contain the closed hysteresis loop according to Equation (16), the loading cycle is completed, its parameters are saved and the strain energy is calculated as presented in Section 2.2. Then, the two points related to the strains ε i and ε i−1 are removed from the strain time series ε(t) and the current strain point index is returned to the value of i = 2. By following this procedure, the closed hysteresis loops are extracted from the strain time series until only three strain points remain. These three points form the outer-most hysteresis loop that is spanned over the complete strain range from ε min to ε max (see Points 1, 4 and 7 in Figure 3). If the outermost hysteresis increases, then also the minimum and maximum strains ε min and ε max change.
The tensile and compressive loading paths that are based on Equations (3) and (4) are modelled as described below. To simplify the procedure, the individual elastic-plastic loading paths are shifted fully into the first quadrant of the σ-ε space, as presented in the right-hand diagram of Figure 1 for the plastic loading paths.

Modelling the Tensile Loading Paths
The tensile loading path is modelled according to the following procedure: • If the tensile loading path starts from the value of ε min , it is the outermost tensile loading path σ T;i−1,i (ε) with the parameter ∆ε max,cyc in Equation (4) being equal to ∆ε max,cyc = |ε max − ε min | (see the σ-ε curve segment 4-7 in Figure 3). • If the tensile loading path originates from the outermost compressive loading path, but not from ε min (as is the case for the segment 2-3 in Figure 3), its tensile loading path σ T;i−1,i (ε) is modelled with Equation (4) and the parameter ε max,cyc is calculated as follows: According to the case in Figure 3, ∆ε max,cyc for the first loading cycle is ∆ε max,cyc, If the tensile loading path originates from the inner compressive loading path, the corresponding parameter ∆ε max,cyc is determined first by extending the line with a slope of the elastic modulus E from the point (ε i−1 , σ i−1 ) to the outermost compressive loading path (see the left diagram in Figure 4): where ε * connects the reversal point with the intersection of the tangent to the tensile loading path and the outermost compressive loading path, as shown in Figure 5. The tensile loading path σ T;i−1,i (ε) of the loading cycle (e.g., the segment 4-5 in Figure 4) is then modelled as a linear combination of the tensile loading path σ T (ε) from Equation (4) with the corresponding parameter ∆ε max,cyc from Equation (18) and the compressive loading path σ C;i−2,i−1 (ε) of the previous strain range (the segment 3-4 in Figure 4) with both of the curves modelled in the first quadrant of the σ-ε space: The model for the compressive loading path σ C;i−2,i−1 (ε) is described below. The mixing weight w T;i is defined as follows (see also Figure 4): • For the inner loading cycles, part of which is also the tensile loading path, the rule of closed hysteresis loops must hold, as defined by Jayakumar [34]. In our case, this means that a strain shift ε shift;i must be determined for the inner tensile loading path σ T;i−1,i (ε) to ensure that its model, which starts at the point ε i−1 , reaches exactly the point ε i−2 of the preceding compressive loading path σ C;i−2,i−1 (ε). To calculate the strain shift ε shift;i , the tensile loading path is modelled in the first quadrant of the σ-ε space (see the right-hand diagram in Figure 4 (In general, the quantity ε shift;i is not equal to the quantity ε * in Figure 4)). • Given that the time series of the previous stresses σ i−1 to σ 1 corresponding to the strain values ε i−1 to ε 1 are known, the highest stress σ i in the current tensile loading path that corresponds to the strain ε i is calculated as follows: After the tensile loading path σ T;i−1,i (ε) of the loading cycle and the value of σ i are known, either the next compressive loading path σ C;i,i+1 (ε) is determined if it does not close the loading cycle, or the compressive loading path σ C;i,i−1 (ε) is determined that closes the loading cycle.

Modelling the Tensile Loading Paths
The compressive loading path is modelled according to the following procedure: • If the compressive loading path starts from the value of ε max , it is the outermost compressive loading path σ C;i−1,i (ε), which is modelled by Equation (3) (e.g., see the σ-ε curve segment 1-2 in Figure 3). • If the compressive loading path does not originate from the strain ε max , regardless of its basic tensile loading path, its model depends on the value of the reference strain If the reference strain ε ref is smaller than the value ∆ε max,cyc − ε max + ε slpT,min from the preceding tensile loading path (e.g., the σ-ε curve segment 2-3 in Figure 5), the model of the compressive loading path is equal to the model of the preceding tensile loading path, ε slpT,min;i−2,i−1 represents the strain that corresponds to the minimum slope of the preceding tensile loading path σ T;i−2,i−1 (ε). If the reference strain ε ref is larger than the value ∆ε max,cyc − ε max + ε slpT,min from the preceding tensile loading path (e.g., the segment 4-5 in Figure 3), the compressive loading path σ C;i−1,i (ε) is a linear combination of the outermost compressive loading path from Equation (3) (the segment 1-2 in Figure 3 or Figure 5) and the Ramberg-Osgood term from Equation (4): The mixing weight w C;i is defined as follows (see also Figures 3 and 5): • As with the tensile loading path model, the strain shift ε shift;i must also be determined for the compressive loading path σ C;i−1,i (ε) in order for its extension to reach exactly the point (ε i−2 , σ i−2 ) of the preceding tensile loading path σ T;i−2,i−1 (ε) (see Figure 5). • Given that the time series of the previous stresses σ i−1 to σ 1 corresponding to the strain values ε i−1 to ε 1 are known, the lowest stress σ i in the current compressive loading path that corresponds to the strain ε i is calculated as follows: In Equation (25), the compressive loading path must be modelled in the first quadrant of the σ-ε space.
After the compressive loading path σ C;i−1,i (ε) of the loading cycle and the value of σ i are known, either the next tensile loading path for a particular loading cycle σ T;i,i+1 (ε) is determined, if it does not close the loading cycle, or the tensile loading path σ T;i,i−1 (ε) that closes the loading cycle is determined.

Strain-Energy Density Calculation
When all the closed hysteresis loops have been modelled, their plastic strain-energy densities can be calculated according to Equation (8). Using the linear damage-accumulation rule (Palmgren-Miner), the fatigue life can be calculated from Equation (13).

Experimental Data
To build a model for predicting the low-cycle fatigue life for variable-strain time series, the following experiments involving the AZ31 magnesium alloy were carried out: • Static tensile tests were performed up to the rupture of the specimen with an on-line measurement of the engineering stress σ and strain ε during the experiment (four specimens). One specimen was used for metallographic experiments in the zone of maximum strain. • Strain-controlled fully reversal (R = −1) low-cycle fatigue experiments with constant strain amplitudes ε a were performed to determine the shape of the closed hysteresis loops at different loading levels and the energy fatigue-life curve: ε a = 0.25% (four specimens), ε a = 0.5% (three specimens), ε a = 0.75% (one specimen), ε a = 1.00% (three specimens), ε a = 1.25% (one specimen) and ε a = 1.5% (one specimen). In each experiment, the number of loading cycles to a break N was acquired. One specimen was also used for the metallographic experiments. • Strain-controlled, fully reversal (R = −1), step-wise, low-cycle fatigue experiments with increasing strain amplitudes ε a were performed until failure. At each strain amplitude level, 20 loading cycles were completed to stabilise the hysteresis loops. The initial strain amplitude was ε a = 0.1% and the second strain amplitude level was ε a = 0.25%. After that loading block, the strain-amplitude level was increased by 0.25% in each of the next loading blocks. The step-wise test was terminated after the rupture of the specimen. The purpose of this test was to determine the shape of the closed hysteresis loops and the hardening rule (isotropic or kinematic). Altogether, five step-wise experiments were performed. • Strain-controlled, variable, low-cycle fatigue experiments for two different strain time histories were performed, as presented in Figure 6. Altogether, two variable-load, low-cycle fatigue experiments were performed between the limit strains of ε min = −1.5% and ε max = 1.5% to validate the proposed model for predicting the low-cycle fatigue life of the magnesium alloy AZ31. All specimens were cut from a 2-mm-thick sheet plate of the AZ31 alloy in the "as-received" condition without any additional processing. The dimensions of the sheet plate were 200 mm × 200 mm. The specimens were fine cut with a water jet in the direction of rolling and in the direction perpendicular to the rolling. The specimen surface was not additionally mechanically treated. Their shape and dimensions are presented in Figure 7. The basic specimen shape, as defined in the ASTM E606 [35], was modified with its dimensions selected in such a manner that: (i) the number of specimens was maximised; (ii) the specimens can be attached to the testing machine; and (iii) an extensometer can be mounted onto the specimens. The tensile and fatigue experiments were carried out on a 100-kN MTS hydraulic testing machine. To measure the strains, an MTS 834.11F-24 extensometer was attached to the middle part of the specimens. The force was measured with a 100-kN tensile-compressive load cell that is integrated into the testing machine. An in-house-developed guiding device was used to prevent any buckling of the specimens [37,38]. The experimental arrangement is presented in Figure 7. The testing frequency for the low-cycle fatigue experiments was varied according to the applied strain-amplitude levels ε a . For the smallest strain amplitudes, the testing frequency was 1 Hz, and, for the highest strain amplitudes, it was 0.1 Hz. The relative experimental errors of the hydraulic test machine were 0.2% and 0.3% for the force and the strain measurements, respectively. The maximum obtained stress was therefore 250 ± 0.5 MPa and the maximum obtained strain was 0.015 ± 0.000045. It was assumed that the experimental error negligibly affected the measurements.
The analysis of the microstructure was performed to corroborate the difference between the deformation mechanisms in AZ31 during the static and dynamic loading, which was the initial motivation for the development of the material model. For the metallographic analysis, the samples were cut in the transversal direction and mounted in Bakelite. The samples were first plane-ground with SiC paper and later fine-ground with a composite disc. For that, a diamond suspension and lubricant were used. Later on, the samples were polished. For polishing, a diamond suspension with a grain size of 1.0 µm was used. The samples were chemically etched with a solution of 99.5% ethanol, distilled water, acetic acid and picric acid.

Experimental Results and Their Modelling
After the static tensile and low-cycle fatigue tests were performed, a portion of each specimen, which was cut out from the most deformed part of the specimen, was crystallographically checked. The pictures of the deformed AZ31 specimens made with the optical microscope and the SEM are presented in Figure 8. The microstructure after the tensile test (Figure 8, top left) is almost free of twins, since it is known that in magnesium alloys with a c/a ratio of less than 1.732 the twinning occurs under compressive loading perpendicular to the c-axis. The microstructure after the low-cycle fatigue is shown in Figure 8 (top right). Twinning hardly occurs when the sample is under tensile stress. However, when the compressive stress is high enough to activate the twinning, the fraction of twins increases with a higher strain. In tension, most twins disappear; only a small portion of twins remain in the alloy. The final microstructure with a large portion of twins shows that the samples were in the compression state. The fracture surfaces, after both tensile and low-cycle fatigue tests, from the SEM show the characteristics of a trans-granular ductile fracture with micro-void coalescence (Figure 8, bottom left and bottom middle). Additionally, fatigue striations can also be seen in Figure 8 (bottom right).
When performing the low-cycle constant-amplitude tests and low-cycle step-wise tests, a considerable scatter in the number of load repetitions until failure was observed. An example of a hysteresis-loop scatter between different tests at an ε a = 0.75% strain amplitude is presented in the left-hand diagram of Figure 9. A similar scatter was also observed for the other strain-amplitude levels. There was also a considerable scatter in the measured fatigue lives for the specimens with the same loading level. The range of the fatigue-life scatter can be seen in Figure 10, in which the number of loading cycles to failure is presented for different plastic and total strain-energy densities. We can see in this figure that the scatter is almost half of the order of magnitude along the abscissa axis. This phenomenon was also observed during the step-wise tests that were repeated for five different specimens. Namely, only twice was the fatigue failure at the same loading level, i.e., at ε a = 1.25%. For the other three specimens, the loading levels at rupture were ε a = 0.75%, 1.0% and 1.75%. One cause for such fatigue-life scatter is definitely the specimen orientation. For this reason, some specimens were cut from the AZ31 plate at an angle of 90 degrees relative to the rolling direction. These results imply that it is difficult to predict the exact fatigue life for a variable-loaded AZ31 specimen. However, despite the scatter of the fatigue lives, it turned out that the Young's modulus E from the static and the low-cycle tests at small values of the strain was always between E = 41 GPa and 45 GPa, with an average value of 43.5 GPa. The latter value was used in further data processing. Figure 10. Energy fatigue-life curves: plastic strain-energy density (circles) and total strain-energy density (squares).
In the right-hand diagram of Figure 9, a family of stabilised hysteresis loops from one step-wise experiment is presented. It can be concluded from this diagram that the plastic hardening mechanism is changed if the strain amplitude levels are increased. For small strain amplitudes up to the value of ε a = 1.0%, the prevailing hardening mechanism is isotropic hardening. Between the strain amplitudes of ε a = 1.0% and ε a = 1.25%, there is a transition from isotropic hardening to kinematic hardening, which prevails at strain amplitudes larger than 1.25%. This conclusion is supported by the right-hand diagram in Figure 11, which shows compressive loading paths from the low-cycle experiments in the first quadrant of the σ-ε space. The changing hardening mechanisms have the following implications for predicting the fatigue life for variable loads: 1.
The stabilised hysteresis loops only (up to the outermost hysteresis-loop envelope) should be considered for predicting the number of loading cycles to failure.

2.
After a certain strain-amplitude level is reached during the cyclic loading, it is its (the outermost) hysteresis loop that governs the stress-strain behaviour, also for the smaller strain-amplitude loading cycles.

3.
The parameters of the Ramberg-Osgood relationship for the outermost compressive loading path that is modelled with Equation (3) should be determined either just for the highest strain amplitude loading level or for all the compressive loading paths from the domain of the kinematic hardening, i.e., ε a ≥ 1.25% in our case.

4.
To build a general model for the tensile loading paths from Equation (4) by considering all the experimental data, the tensile loading paths as presented in the left-hand diagram of Figure 11 must first be extended to the outer-most compressive loading path. Only after this transformation can the parameters of Equation (4) be determined. Finally, two variable-strain time-series experiments were performed. The stabilised σ-ε response for the complete period of the simpler variable test is presented at the top of Figure 12. This response corresponds to the ε(t) time series from the left-hand diagram of Figure 6, which is derived from one of the time series from [27]. The stabilised σ-ε response for the complete period of the second variable signal in Figure 6 is presented at the bottom of Figure 12. We can see from Figure 12 that the tensile and compressive loading paths form an envelope for all the other hysteresis loops that are present in the two variable signals. We included the outermost tensile and compressive loading paths from the two variable tests into Figure 11 to compare them with the constant strain-amplitude experiments. It can be seen from both diagrams in Figure 11 that the loading paths from the variable test agree well with all the other low-cycle fatigue data. It took 29 repetitions of the period for the first variable load until the fatigue failure, but only 13 repetitions until failure for the second variable load. In both cases, the σ-ε response for the complete loading period was stabilised after three loading cycles. For this reason, the stabilised σ-ε response at the half repetitions of the variable signals was considered for the data analysis, i.e., the 15th loading period for the first variable load and the 7th loading period for the second variable load. The model for the compressive loading paths according to Equation (3) was built first. The parameters n C and K C were estimated with a real-valued genetic algorithm (RVGA) that was applied before for modelling different fatigue-life curves and their scatter. For details on the theoretical background and the selected parameters that govern the optimisation procedure, see the works of Klemenc and Fajdiga [39] and Klemenc [40]. The objective function was the sum of the squared distances between the measured and the modelled data points: The training database consisted of k C = 704 data points (ε j , σ j ) that belong to all the measured compressive loading paths of the constant amplitude tests, the step-wise tests and the compressive envelopes of the two variable σ-ε responses for ε a ≥ 1.25%. Ten RVGA runs with variable initial conditions and 10,000 iterations were made to reach the final estimates of the n C and K C parameters for the value of the Young's modulus E = 43,500 MPa (see Table 1). The model of the compressive loading path is presented in the right-hand diagram of Figure 11. It can be seen in this figure that the model represents the data well and is a good trade-off between the compressive loading paths from different kinds of low-cycle fatigue experiments. 28.395 D [-] 523.29 f 1 [-] 0.95959 f 2 [-] (0.38102) After the model for the compressive loading paths was built, the tensile loading paths from the left-hand diagram in Figure 11 were first scaled so that their highest points coincide with the modelled compressive loading path. Then, the scaled tensile-loading paths from all the experiments (the constant amplitude tests, the step-wise tests and the compressive envelopes of the two variable σ-ε responses) were joined into a training database with k T = 3632 data points (ε j , σ j ). To each data point, a weight was added to improve the parameter-estimation procedure, since the density of the data points is much smaller at the end of the loading paths. The data-point weight w DP;j depended on the data point's stress magnitude: With such a modification, the contribution of the data points at different stress levels is approximately equalised, because the point density along the stress axis varies significantly.
Finally, the parameters K T , n T , b 1 , b 2 , D, f 1 and f 2 of the tensile loading paths were estimated with the same version of the RVGA as before. The objective function was a weighted sum of the squared distances between the measured and the modelled data points, but for the tensile loading path in this case the following was applied: (29) Ten RVGA runs with variable initial conditions and 10,000 iterations were made to reach the final estimates of the seven parametersK T , n T , b 1 , b 2 , D, f 1 and f 2 (see the last column in Table 1). Again, the model represents the data well and is a good trade-off between the tensile loading paths from different types of low-cycle fatigue experiments. The tensile loading paths are not presented in the article in order to avoid data overload in the figures. The value of the parameter f 2 in Table 1 is put into parentheses because it is less reliable than the other parameters. The reason for this is that there was no well-expressed untwining saturation in the tensile loading paths below the maximum strain amplitude of ε a = 2.0%.

Prediction of Fatigue Life Using the Strain-Energy Density
From the 13 fatigue-life experiments with a constant amplitude strain, the energy fatigue-life curves were estimated. As described in Section 2.2, the plastic strain-energy density ∆W p and total strain-energy density ∆W t were calculated from the stabilised hysteresis loops at the half fatigue life. In Figure 10, the (∆W p/t , N f ) points are presented together with the corresponding energy fatigue-life curves. The material-dependent parameters of the two energy fatigue-life curves are as follows: • plastic strain-energy density ∆W p : C p = 537.52, m p = 1.0705; and • total strain-energy density curve ∆W t : C t = 153.80, m t = 0.7627.
These energy curves were applied to predict the fatigue-life damage for the two variable-strain time series. To predict the fatigue life, the closed hysteresis loops were first extracted from the ε(t) history, which was followed by modelling their tensile and compressive loading paths, as described in Section 2.3. There are k cyc = 10 closed hysteresis loops in the first variable signal and k cyc = 21 hysteresis loops in the second variable signal. The modelled hysteresis loops for both the variable signals are presented in the right-hand diagrams of Figure 12. We can conclude from the results in Figure 12 that the agreement between the measured and the modelled σ-ε responses is very good if we consider the amount of scatter related to the fatigue-life experiments. The only significant deviation between the measured and modelled responses is in the lower-right-hand corner of the two diagrams, where the compressive loading path changes its slope. Even in this area, the deviation between the measured and modelled responses is less than 10% along the stress axis. The quality of the agreement is confirmed by the data in Tables 2 and 3. If we compare the maximum stresses σ max of the closed hysteresis loops, we can see that the deviation between the measured and modelled data is always less than 10% for the first variable strain history, and that only four out of 21 hysteresis loops have the deviation of the maximum stress larger than 10%. However, these are the smallest loading cycles that do not influence the fatigue life significantly.  After the hysteresis loops were extracted from the measured and modelled σ-ε responses, the strain-energy densities ∆W p and ∆W t were calculated for each hysteresis loop. In Table 2, the calculated plastic and total strain-energy densities are listed for the measured and modelled hysteresis loops of the first variable-strain history. In Table 3, the same data are listed for the second variable-strain history. The comparison between the measured and the modelled values are in good agreement, apart from a few smallest cycles which negligibly contribute to the total damage.
The energy fatigue-life curves in Figure 10 were applied to calculate the damage Dmg i on the basis of the data in Tables 2 and 3. The fatigue-life damage that is caused by one repetition of each variable-strain history is equal to: If the critical fatigue damage Dmg c is equal to 1.0, the predicted number of variable-signal repetitions until fatigue failure is (see also Section 2.2): During the low-cycle fatigue testing, the first (simpler) variable-strain history was repeated 29 times until the failure of the specimen. The predicted fatigue life on the basis of the plastic strain-energy density ∆W p was 60 repetitions until failure for the measured and modelled data. On the other hand, the predicted fatigue life on the basis of the total strain-energy density ∆W t was 74 repetitions for the measured data and 73 repetitions for the modelled data. The second variable-strain history was repeated 13 times until failure. The predicted fatigue life on the basis of the plastic strain-energy density was 15 (measured data) and 16 (modelled data) loading cycles. If the total strain-energy density was considered, the predictions were 18 (measured data) and 20 (modelled data) loading cycles. We can see that the predictions were always non-conservative. Regardless of the method, the fatigue-life predictions from the measured and modelled data differed by no more than two repetitions of the variable-strain history. This is another confirmation that the presented model of the AZ31 hysteresis loops is very good. From the presented results, it can be concluded that in the case of the low-cycle fatigue regime the AZ31 magnesium alloy does not show a significant mean-stress effect, because for both variable-strain histories consistently better fatigue-life predictions were made on the basis of the plastic strain-energy density ∆W p . Since the scatter of the energy fatigue-life curve along the number of loading cycles to failure N is considerable (see Figure 10), the prediction is in its scatter range, even for the first variable signal, despite the fact that the predicted fatigue life was twice the measured fatigue life. This means that the energy-based approach combined with the presented model of the hysteresis loops is appropriate for predicting the fatigue life of the AZ31 magnesium alloy.
The prediction model has been applied to AZ31 magnesium alloy during this study. However, the methodology is applicable also to other materials that exhibit a similar asymmetric stress-strain response. The material parameters of the model have to be determined from the LCF tests for the material under investigation and then the stress-strain response and the fatigue damage can be predicted.

Conclusions
In the article, a method for predicting the fatigue life of magnesium alloys is presented, based on the energy approach. With the introduction of an analytical model for describing the tensile and compressive loading paths of the closed hysteresis loops, it is possible to predict the fatigue life for an arbitrary strain-loading history on the basis of the plastic or total strain-energy density. The research outcomes that are presented in the article can be summarised as follows:

•
A new model for tensile and compressive elastic-plastic hysteresis loops is presented, which enables the modelling of an arbitrary hysteresis loop in the low-cycle fatigue domain with only 10 parameters. • Static, dynamic and metallographic experiments were carried out to determine the elastic-plastic response of the AZ31 alloy and to estimate the influence of the loading direction on the microstructure of thin AZ31 sheets. • The presented hysteresis-loop model was validated against experimental fatigue-life data that emerge from different tests: constant strain-amplitude testing, step-wise variation of the strain-amplitude levels and the arbitrary variation of strain amplitudes in a variable signal. • An energy fatigue-life curve was determined on the basis of the experimental data that were applied for predicting the fatigue life in the low-cycle fatigue domain. • It was proven that the fatigue lives that are predicted on the basis of the modelled hysteresis loops are equivalent to the fatigue lives that are predicted on the basis of the measured stress-strain response. • It was shown that better fatigue-life predictions were obtained on the basis of the plastic strain-energy densities than the total strain-energy densities.

Acknowledgments:
The authors would like to thank Tomaž Bešter for his technical support in performing the low-cycle fatigue experiments.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

Abbreviations
The following abbreviations are used in this manuscript: