On the applicability of multi-surface, two-surface and non-linear kinematic hardening models in multiaxial fatigue

In this work, a comparison between NLK and Mróz-Garud’s multi-surface formulations is presented. A unified common notation is introduced to describe the involved equations, showing that the Mróz-Garud model can be regarded as a particular case of the NLK formulation. It is also shown that the classic two-surface model, which is an unconventional simplified plasticity model based on the translation of only two surfaces, can also be represented using this formulation. Such common notation allows a direct quantitative comparison among multi-surface, two-surface, and NLK hardening models.


INTRODUCTION
he Bauschinger effect, commonly called kinematic hardening, can be modeled in stress spaces by allowing the yield surface to translate with no change in its size or shape. So, in the deviatoric stress space, kinematic hardening maintains the radius S of the yield surface fixed while its center is translated, changing the associated generalized plastic modulus P that defines the slope between stress and plastic strain increments in the Prandtl-Reuss plastic flow rule, also known as the normality rule. There are several models to calculate the current value of P as the yield surface translates, as well as the direction of such translation, to obtain the associated plastic strain increments. Most of these hardening models can be divided into three classes: Mróz multi-surface (or multi-linear), non-linear, and two-surface kinematic models. Mróz [1] defined in 1967 the first multi-surface kinematic hardening model to approximately describe the behavior of elastoplastic solids through a family of nested surfaces in the stress space, the innermost being the yield surface associated with the material yield strength. It assumes that P is piecewise constant, resulting in a multi-linear description of the stressstrain curve, i.e. the non-linear shape of the elastoplastic stress-strain relation is approximated by several linear segments. The Mróz model can induce a few numerical problems, which can result in hardening surfaces improperly intersecting in more than one point under finite strain increments. Garud [2] proposed a geometrical correction that avoids intersection problems even for coarse integration increments. In spite of its limitations, several multiaxial fatigue works use the Mróz-Garud model to predict the stress-strain behavior under combined loading, especially due to its ability to store plastic memory effects under variable amplitude loading. This T popularity is understandable, since the multi-linear stress-strain curves generated by the Mróz-Garud model provide good results for balanced loadings. However, such multi-linear models cannot predict any uniaxial ratcheting or mean stress relaxation caused by unbalanced loadings, since their idealized hysteresis loops always close because they are unrealistically assumed as perfectly symmetric. In addition, under several non-proportional loading conditions, these models predict multiaxial ratcheting with a constant rate that never decays, severely overestimating the ratcheting effect measured in practice [3]. As a result, multi-surface kinematic hardening models should only be confidently applied to balanced proportional loading histories. To correctly predict the stress-strain history associated with unbalanced loadings, it is necessary to introduce non-linearity in the hardening surface translation equations, the main characteristic of the non-linear kinematic (NLK) models. NLK models are more general than Mróz-Garud because they use non-linear equations to describe the surface translation direction and the value of P, leading thus to a more precise description of the non-linear stress-strain curves. Armstrong and Frederick's original formulation [4] was improved by Chaboche [5], which indirectly introduced the Mróz nestedsurface idea to NLK models, however in a non-linear instead of multi-linear formulation. Therefore, both NLK and Mróz-Garud formulations have several common features, as discussed further in this work. A third class of kinematic hardening models involves the so-called two-surface models, which use a rather simplified formulation that combines elements of both non-linear and Mróz multi-surface kinematic models. In this work, instead of defining the nested hardening surfaces in the 6D stress or 6D deviatoric stress spaces, a 5D reduced order deviatoric stress space E 5s is adopted, using the Mises yield function to describe each surface. This 5D space has two advantages over the usual 6D formulations: it is a non-redundant representation of the deviatoric stresses, which decreases the computational cost of stress-strain calculations; and the radius S of the yield surface is equal to the yield strength without the need to include the scaling factor 2 3 required in 6D formulations. Besides, even though all kinematic hardening equations are presented here in the 5D space, their conversion to 6D versions is trivial.

MRÓZ MULTI-SURFACE MODELS
n Fig. 1, the first (and innermost) circle is the monotonic yield surface, with radius r 1  S Y . In addition, M 1 hardening surfaces with radii r 1 < r 2 < … < r M  1 are defined, along with an outermost failure surface whose radius r M  1 is equal to the true rupture stress U of the material. Their centers are located at points ci s  I Except for the failure surface, all other hardening surfaces can translate as the material strain-hardens, as shown in Fig.  1(right). The centers of the hardening surfaces (which are circles in the 2D example from Fig. 1) move as the material plastically deforms and hardens, because they are successively pushed by the inner surfaces. The radii r i of the various hardening surfaces are equal to the stress levels associated with the plastic strains  xpli that delineate the multi-linear representation of the stress-strain curve, fitted to properly describe the stress-strain  behavior of the material. The difference between the radii of each pair of consecutive surfaces is defined as r i  r i  1  r i . In principle, all hardening surfaces radii ri may change during plastic deformation as a result of isotropic and NP hardening effects.
) trying to move outwards; and (ii) the hardening surfaces cannot cross through one another, therefore they gradually become mutually tangent to one another at the current stress point s  as the material plastically deforms. In the Mróz multi-surface formulation, the outermost surface that is moving at any instant is called the active surface, denoted here as the surface with index iA. Any changes in the stress state that happen inside the yield surface are assumed elastic, not resulting in any surface translation as long as , therefore no surface would be active in this case and thus i A  0. Each surface is associated with a generalized plastic modulus P i (i  1, 2, …, M  1), which altogether define a field of hardening moduli. The value of P is then chosen as the Pi from the active surface i iA. Fig. 3 illustrates two consecutive surfaces i and i  1 (i  1), with radii r i and r i  1 in the E 5s 5D deviatoric stress space. Assuming that the current active surface is i  i A , then according to the Mróz model all outer hardening surfaces do not translate, therefore the increments of the respective backstresses are 1 2 ... 0 The Mróz multi-surface formulation assumes that, during plastic straining, all inner surfaces 1, 2, …, i A 1 must translate altogether with the active surface i  iA, therefore their centers do not move relatively to each other, resulting in Thus, translation rules in the Mróz multi-surface formulation only need to be applied to the evolution of the backstress For Mises materials, the translation direction of Prager-Ziegler's kinematic rule is then But Prager-Ziegler's translation rule is too simplistic to model kinematic hardening. An improved rule was adopted by Mróz [1], who assumed that the translation i d    of the active surface occurs in a direction i v  defined by the segment that joins the current stress state with the corresponding "image stress point" that has the same normal unit vector n  at the next hardening surface i 1, see Fig. 3. The Mróz translation direction is then Note that the Mróz translation direction is the combination of Prager-Ziegler's hardening term with a "dynamic recovery" direction i     . This last term induces the center of the active surface to translate back towards the center of the next hardening surface, attempting to dynamically recover from the hardening state described by the surface backstress i    that separates their centers. This dynamic recovery term is able to consider a fading memory of the plastic strain path, which is necessary to model mean stress relaxation and ratcheting effects. However, the Mróz rule can induce a few numerical problems, which can result in surfaces intersecting in more than one point under finite load increments. Garud proposed an empirical correction that avoids intersection problems even for coarse integration increments, as detailed in [2]. A major concern of the Mróz multi-surface formulation is that the directions of the calculated stress or strain paths may significantly vary depending on the number of surfaces used. Even worse, better predictions are not necessarily obtained from using a larger number of surfaces. As a result, the number of hardening surfaces that results in the best calculation accuracy is a finite number that would also need to be calibrated, an undesirable feature. Moreover, the multi-linear formulation adopted by Mróz and Garud models is not able to correctly predict ratcheting and mean stress relaxation effects. A better approach is to replace such multi-linear models with a non-linear kinematic hardening formulation, described next.

NON-LINEAR KINEMATIC (NLK) HARDENING MODELS
he multi-linear stress-strain curves generated by the Mróz and Garud multi-surface models provide good results for balanced proportional loadings, which by definition do not induce ratcheting or mean stress relaxation. However, such piecewise linear models cannot predict any uniaxial ratcheting or mean stress relaxation effects caused by unbalanced proportional loadings. This shortcoming is due to the linearity of the Mróz and Garud surface translation rules and the resulting multi-linearity of the approximated stress-strain representation, which describes all elastoplastic hysteresis loops using multiple straight segments, instead of predicting the experimentally observed curved paths caused by those non-linear effects. Such straight segments generate unrealistic perfectly symmetric hysteresis loops that always close under constant amplitude proportional loadings. In addition, for NP loadings, the Mróz and Garud multi-surface models may predict multiaxial ratcheting with a constant rate that never decays, severely overestimating the ratcheting effect measured in practice. As a result, Mróz or Garud kinematic strain-hardening models should only be applied to balanced loading histories, severely limiting their applicability. Another serious flaw of the Mróz and Garud models becomes evident for a NP path example where the stress state s  follows the contour of the active yield surface. In this example, the stress increment ds  would always be tangent to such active surface, therefore it would induce 0 T ds n       and no plastic strain would be predicted from the normality rule. This conclusion is physically inadmissible, since it would assume zero plastic straining along a load path with radius rA larger than the yield surface radius r1. The Mróz and Garud models would predict that the yield surface can translate tangentially to the active surface without generating plastic strains. These major drawbacks are a consequence of Mróz and Garud multi-surface kinematic hardening models being of an "uncoupled formulation" type, as qualified by Bari and Hassan in [8]. Such uncoupling means that the generalized plastic modulus P  P i in this formulation is not a function of the load translation direction i v  , since it is a constant for each surface. Such "uncoupled procedure" (where P and i v  are independent) provides undesirable additional degrees of freedom to the Mróz and Garud models that allow, for instance, 90 o out-of-phase tension-torsion predictions of resulting plastic strain amplitudes that are not a monotonic function of the applied stress amplitudes, as they should be [9]. These T wrong Mróz and Garud predictions are both qualitatively and quantitatively dependent on the number of surfaces adopted in the model, without any clear convergence.
To correctly predict the stress-strain history in unbalanced loadings, P and i v  must be coupled, in addition to introducing non-linearity in the surface translation equations, generating the non-linear kinematic (NLK) models discussed next. The first non-linear kinematic (NLK) hardening model was proposed by Armstrong and Frederick in 1966 [4]. Their original single-surface model did not include any hardening surface, but their yield surface already translated according to a nonlinear rule. Chaboche et al. [5] significantly improved Armstrong-Frederick's model capabilities by indirectly introducing the concept of multiple hardening surfaces. As demonstrated by Ohno and Wang in [10], this improvement allowed the NLK formulation to use the same representation of the hardening state as the one in the Mróz multi-surface formulation, which includes one yield surface, M 1 hardening surfaces, and one failure surface, all of them nested within each other without crossing. Therefore, the improved NLK models also adopt a multi-surface formulation, using a "coupled procedure" based on non-linear translation rules. However, unless otherwise noted, the denomination "multi-surface model" is traditionally associated in the literature with the Mróz "uncoupled procedure" based on multi-linear translation rules and piecewise-constant generalized plastic moduli P i [8].
In summary, despite their significant differences, both Mróz and NLK approaches can be represented using the same multi-surface formulation. As in the Mróz models, instead of defining these surfaces in the 6D stress or deviatoric stress spaces, in this work a 5D reduced order deviatoric stress space E5s is adopted, using the Mises yield function to describe each surface.
In the multi-surface version of NLK models, the backstress   ...
Once again, the length i    of each hardening surface backstress is always between 0 (in the unhardened condition) and its saturation value ri (a maximum hardening condition when these surfaces become tangent, see Fig. 3 As usual, the hardening surfaces cannot pass through one another, remaining nested at all times, since | | . But during plastic straining in the NLK multi-surface formulation, the yield and all hardening surfaces do translate, as opposed to the Mróz formulation, where all surfaces outside the active one would not move. The NLK models do not use an "active surface" concept, since all hardening surfaces become active during plastic straining. The yield and hardening surfaces from the NLK models behave as if they were all attached to one another with non-linear spring-slider elements, causing coupled translations even before they enter in contact, i.e. even for | | i i r      . Therefore, any hardening surface translation causes all surfaces to translate, usually with different magnitudes and directions, even before they become tangent to each other. Pairs of consecutive hardening surfaces i and i  1 may eventually become mutually tangent at the saturation condition where the scalar functions  i * and m i * are defined as The calibration parameters for each surface i are the ratcheting exponent  i , the multiaxial ratcheting exponent m i , the ratcheting coefficient i, and the multiaxial ratcheting coefficient i, scalar values that are listed in Tab. 1 for several popular NLK hardening models. Note that several literature references represent the NLK hardening parameters r i , p i , and  i as r (i) , c (i) , and  (i) , but this notation is not used in this work to avoid mistaking the (i) superscripts for exponents.

Year
Chen-Jiao-Kim [20] 0   i < ∞ ∞ < m i < ∞ 1 1 The translation direction i v  of each hardening surface i, shown in Eq. (7), can be separated into three components: (i) the Prager-Ziegler term, in the normal direction n  perpendicular to the yield surface at the current stress point s  ; (ii) the dynamic recovery term, in the opposite direction of the backstress component of the considered surface, which acts as a recall term that gradually erases plastic memory with an intensity proportional to the product of the ratcheting terms  i * m i *  i  i ; and (iii) the radial return term, in the opposite direction of the normal vector n  , which mostly affects multiaxial ratcheting predictions. Fig. 4 shows the geometric interpretation of these three components.

TWO-SURFACE KINEMATIC HARDENING MODEL
he two-surface model proposed by Dafalias and Popov [21] and independently by Krieg [22] is an unconventional plasticity model based on the translation of only two moving surfaces: a single hardening surface (i  2), usually called bounding or limit surface, and an inner yield surface (i  1), shown in Fig. 5. The outer failure surface (i  3) is also present, however it does not translate.

CONCLUSIONS
n this work, it is concluded that the Mróz-Garud multi-surface formulation can lead to very poor multiaxial stressstrain predictions in the presence of significant mean stresses, severely limiting their application in multiaxial fatigue analysis. The two-surface model better accounts for unbalanced loadings, however its application is mostly recommended for monotonic plastic processes, since it does not appropriately deal with complex variable amplitude loading histories, with decreasing amplitudes and several layers of hysteresis loops within loops, which are very common under spectrum loading. On the other hand, NLK models can accurately deal with non-proportional loadings, either balanced or unbalanced. Moreover, contrary to two-surface models, they deal with complex variable amplitude spectrum loading, as long as a sufficiently high number of backstress components (i.e. number of hardening surfaces) is adopted in a refined NLK formulation, to efficiently store plastic memory effects. Therefore, this analysis shows that NLK models should be preferred over multi-surface and two-surface models in multiaxial fatigue calculations under variable amplitude loading.