Stefan C Soare
email: stefsoare@yahoo.com

effects, yield surface influence, and material heterogeneity

Keywords: Marciniak-Kuczynski Model, Forming Limit Diagram, Sheet Metal

The paper presents a study of the Marciniak and Kuczynski (MK for short) model for the prediction of limit strains of orthotropic sheet metal under in-plane proportional biaxial stretching. In two particular cases analytical results can be obtained if the groove of the MK model is oriented along one of the in-plane symmetry axes. The first case is the plane strain loading mode. Necessary and sufficient conditions are derived for the MK-predicted plane strain limit strain to match exactly the experimentally measured limit strain. An example of material, the AA5182-O aluminum alloy, that does not satisfy these conditions is discussed. It is shown then that if a power-law strain rate sensitivity is included in the hardening law then the MK-model can match exactly any target plane strain limit strain. The second case is the non-hardening case for positive strain ratios. This case allows for an insight into the way the MK-predicted limit strains depend upon the yield function. Based on the theory developed for the plane strain case, material heterogeneity as a possible cause for unstable plastic flow is further discussed. It is shown that such heterogeneities can be modeled by perturbing the rate of deformation with an eigenstrain. This allows for an extension of the MK-model to sheets of uniform thickness.

Introduction

As is well known, the constitutive equations of classical plasticity, convex yield surface coupled with the normality flow rule, lead to the paradoxical result that one can stretch indefinitely a metal sheet, without breaking it. In more precise terms, if one assumes a homogeneous continuum under homogeneous boundary conditions, the deformation field predicted by the mentioned classical theory is homogeneous. On the other hand, these equations are hyperbolic in the rigidplastic approximation, [START_REF] Hill | The mathematical theory of plasticity[END_REF], and hence they allow, at least in principle, for the propagation of velocity discontinuities as required by neck formation during biaxial stretching of metal sheets. However, since the characteristics are constant for proportional loading, this theory does not allow for a progressive development of instability: in Hill's theory, necking is instantaneous and takes place along zero extension directions. Hence the phenomenology of necking developed by [START_REF] Hill | On discontinuous plastic states, with special reference to localized necking in thin sheets[END_REF][START_REF] Hill | On the mechanics of localized necking[END_REF] describes the bifurcation point along the loading history.

A more successful, from a practical point of view, model for sheet failure was proposed later by [START_REF] Marciniak | Limit strains in the process of stretch-forming sheet metal[END_REF] (although the initial focus was on positive strain ratios only, the model was later extended by the Brown school of the 1970's to negative strain ratios by allowing oblique grooves, see next). While retaining the homogeneous continuum assumption, and the classical rigid-plastic constitutive framework, the cited work introduces a geometrical heterogeneity in order to simulate neck formation: the sheet is assumed to have a thickness imperfection in the form of a groove across its width, Fig. 1, groove that under external in-plane loads will develop ultimately into a neck (a direction of zero extension, see next). Remarkable, the neck now develops gradually: the groove is nothing but a new characteristic direction etched at the outset onto the material, direction along which discontinuities can propagate. It will be shown here that Hill's failure criterion is nothing but the MK's failure criterion for the material in the groove (neck). This explains why Hill's criterion always overestimates the limit strains for negative strain ratios. Most current calculations of forming limit diagrams (FLD, for short) are based on the MK model. Two features of this model explain its popularity: it is capable of predicting the entire FLD (that is, for positive and negative strain ratios along a proportional strain path), and the initial depth of the groove can serve as a parameter of the model to adjust its predictions with measured data. With few exceptions, e.g., [START_REF] Sowerby | Failure in sheet metal in biaxial tension[END_REF], Hutchinson and Neale (1978 a,b), much of our current understanding about the MK-model can be safely regarded as empirical (based on numerical calculations). This approach cannot explain, for example, why for some materials the parameter of the MK-model cannot be adjusted for an exact prediction of the plane strain limit strain, see for example [START_REF] Wu | Evaluation of anisotropic yield functions for aluminum sheets[END_REF]. One objective of this work is to report some new analytical formulas regarding the plane strain strain case and the non-hardening case of the MK-model. A critique often brought to the MK model, e.g., [START_REF] Storen | Localized necking in thin sheets[END_REF], Barata Da Rocha et al (1984), is that the initial depth of the groove is usually too big to be correlated with a thickness imperfection observed in actual metallic sheets. In other words, the groove of the MK model is far from being an explanation or a real model for the initiation of unstable plastic flow. One can say that something else triggered the instability and then the MK model takes over the description from a later stage when a groove of a certain depth has already developed. The second main objective of this report is to show that the MK-model can indeed be extended to sheets of uniform thickness.

The MK-model

Before proceeding with the announced program, let us first describe briefly the MK-model and the constitutive framework we shall adopt here. Detailed presentations, in a more general constitutive framework, can be found in [START_REF] Kuroda | Forming limit diagrams for anisotropic metal sheets with different anisotropic yield criteria[END_REF], [START_REF] Wu | Evaluation of anisotropic yield functions for aluminum sheets[END_REF], Aretz (2007a). In what follows, the sheet is assumed orthotropic and the orthogonal coordinate system xyz is always oriented along the symmetry axes: x along the rolling direction, y along the transverse direction and z is the normal to the plane of the sheet direction. The material of the sheet is assumed homogeneous, incompressible and rigid-plastic, the yield surface being also flow potential (normality rule). The yield surface is assumed first order positive homogeneous. A dot over a variable/function will always denote material derivative with respect to time. Two spatial areas are distinguished: one is the groove, which will be denoted from now on as zone B, the other is the bulk of the sheet. The values taken by a variable in zone B will be distinguished by a right superscript B attached to them. The sheet is loaded in its plane by a Figure 1: Geometry of the MK-model: the shaded area represents the initial groove, the zone where the thickness of the sheet is assumed to be smaller than the thickness of the rest of the sheet. system of forces collinear with the symmetry axes of the sheet, Fig. 1. Since the boundary conditions are assumed (spatially) uniform and the sheet is homogeneous, we can safely assume the stress-strain fields in the bulk of the sheet as being uniform. The loading process is viewed as a monotonic succession of equilibrium states through infinitesimal increments. The strain increment tensor will be denoted d . The loading conditions are such that the strain path is proportional:

d y d x =: ρ = constant (1)
From incompressibility it follows that the strain increment tensor is

d = d x               1 0 0 0 ρ 0 0 0 -(1 + ρ)               (2)
whereas from the orthotropic symmetry of the sheet it follows that the shear stress components are zero and hence, with the plane stress assumption, the stress state in the bulk of the sheet is

σ = σ x               1 0 0 0 t 0 0 0 0               (3)
Above we have denoted t := σ y /σ x the constant stress ratio in the bulk of the sheet. Given a strain ratio ρ, t is determined using the flow rule, the first order homogeneity of the yield function, and eq. ( 1):

∂ f ∂σ y (1, t) = ρ ∂ f ∂σ x (1, t) (4)
Conversely, for a given stress ratio t, the above relation becomes an equation determining explicitly the strain ratio ρ. Due to the symmetry of the sheet and of the loads, only stress ratios t ∈ [0, 1] will be considered, that is, the range between uniaxial and balanced-biaxial stressing where σ x is the dominant stress component. In this text, f = f (σ) = f (σ x , σ y , σ xy ) will always denote the yield surface. Above, the third argument of the yield function, the shear component σ xy , has been omitted. The equivalent stress σ of a certain stress state σ is defined as the value of the yield function on that stress state: σ := f (σ). The equivalent strain increment, which in our case is also equal to the equivalent plastic strain increment, is defined as the plastic multiplier of the corresponding stress state:

d = dλ, where d = dλ ∂ f ∂σ (σ) (5)
Using the first order homogeneity of the yield function, the flow rule and the above definitions, we have then the following work-equivalence relationship:

σ • d = f (σ)d (6)
With eqs. ( 2) and (3) one obtains a differential equation for ( x ) which integrated with the initial condition

| x =0 = 0 leads to f (1, t)d = (1 + tρ)d x =⇒ = 1 + tρ f (1, t) x (7)
In what follows we shall assume the hardening properties of the material isotropic and hence there exists a universal loading curve σ = H( ) (strain rate sensitivity will be also included later). For the monotonic loading paths considered in this work this assumption is in general satisfactory. Then, from the yielding condition the stress is determined as

f (σ x , σ y ) = H( ) =⇒ σ x = H( ) f (1, t) , σ y = tσ x (8)
Formulas ( 7) and (8) determine the stress in the bulk of the sheet for a given loading path and strain level. Next, we turn our attention to zone B of the sheet. This is characterized relative to the bulk of the sheet by the following imperfection ratio:

f h := h B h ⇐⇒ ln f h = ln(h B /h B 0 ) + ln(h B 0 /h 0 ) -ln(h/h 0 ) ⇐⇒ f h = f 0 exp B z -z (9)
with h B denoting the current (uniform) thickness of zone B, and h the current (uniform) thickness of the sheet outside zone B, while a subscript pertains to the initial values of the two thicknesses, and f 0 is the initial imperfection ratio. According to the phenomenology of necking, the width of zone B can be considered small, of magnitude comparable with the thickness of the sheet, throughout the loading process. Then we shall assume the stress-strain state in zone B as being spatially uniform. This consideration also allows us to take the spatial rigid motion of zone B identical with that of a material line in the bulk of the sheet of the same inclination with zone B. With reference to Fig. 1, with t = (-sin ψ, cos ψ) and n = (cos ψ, sin ψ) we denote the tangential and normal to the groove directions, respectively. Then the orientation of zone B is described by the formula:

tan ψ = exp (1 -ρ) x tan ψ 0 (10)
Since a material element along the boundary between the two zones must experience a unique state of strain, the following compatibility equation must hold

d B t = d t ⇐⇒ d B ∂ f ∂σ σ B : t • t = d ∂ f ∂σ (σ) : t • t (11)
At the boundary between the two zones one must also have force equilibrium. Since in both areas of the sheet the stress state is uniform, this condition amounts to:

h B σ B : n = h (σ : n) ⇐⇒        σ B n = σ n / f h σ B nt = σ nt / f h (12)
With further use of the yielding condition in zone B, eqs. (9-12) determine the stress-strain state in zone B (the groove) at any moment during the loading process. The resulting system of equations is highly nonlinear and hence, in the general case, it is solved numerically. Also numerical is the failure criterion: strain increments are performed upon the bulk of the sheet until the following inequality holds true:

d B /d > N f (13)
with N f a big enough number. The failure strains are the (logarithmic) strains accumulated in the bulk of the sheet up to the failure moment. To the author's knowledge, there is no proof in the literature that the above ratio will always grow unbounded.

The Plane Strain Case

This case is usually employed to determine the MK-parameter, the initial thickness imperfection ratio f 0 , so that the MK-FLD prediction best fits the experimentally determined plane strain limit strain. It is possible in this case to integrate the MK evolution equations into an explicit system of equations determining the failure strains, if one assumes the groove along one of the in-plane symmetry axes. Most materials attain their ρ = 0 MK-limit-strain when zone B is oriented along one of the symmetry axes. For the few exceptions known to the author the deviation of the MK-groove from a symmetry axis is small, and hence the case studied here is relevant even for these materials. Since in this study the major strain (or stress component) is along the x-axis, we assume the groove parallel to the y-axis. According to eq.( 10), it will remain so during the entire loading process. Furthermore, the numerical failure criterion (13) is restated as:

d B d = +∞ (14)
We have then the following result:

Theorem 1. Let ρ = 0, t PS denote the corresponding stress ratio, k := f (1, t PS ), and zone B oriented along the y-axis. For a material described with a hardening law in the form σ = H( ), with the function H strictly concave, the MK-limit strain is x = k F with the equivalent strain at failure F determined by the following system of equations:

H B F -kH B F = 0 (15) f 0 H B F exp -k B F -H ( F ) exp (-k F ) = 0 (16)
Proof. Let us first note that the stress ratio in zone B, t B = σ B y /σ B x , is constant and equal to the one in the bulk of the sheet. Indeed, this follows from the flow rule and compatibility eq. ( 11):

d y = 0 = d B y =⇒ ∂ f ∂σ y (1, t) = 0 = ∂ f ∂σ y (1, t B ) =⇒ t = t B = t PS (17)
Above we have assumed that the plane strain point is unique on the yield surface (in the positive octant of the (σ x , σ y ) plane between the σ y = 0 and σ x = σ y axes). Furthermore, using Euler's identity for first order homogeneous functions we get

σ B • ∂ f ∂σ (σ B ) = f (σ B ) =⇒ ∂ f ∂σ x (1, t) = f (1, t) = k (18)
Then from incompressibility and flow rule we have B z = -k B and z = -k , and hence eq. ( 9) becomes in this case

f h = f 0 exp -k B - (19)
Next, from the yielding conditions in the two zones of the sheet and equilibrium eq. ( 12) we obtain

f σ x f h , σ x f h t = H B ⇐⇒ σ x f (1, t) = f h H B ⇐⇒ H ( ) = f h H B ( 20
)
Substituting above formula (19) leads to eq. ( 16). This relationship allows us to solve for B once is known, as long as solutions exits. We can therefore speak of a one-to-one relationship B = B ( ). Then differentiating in ( 16) with respect to leads to:

1 f 0 H ( ) -kH( ) exp(k ) = H ( B ) -kH( B ) exp(k B ) d B d (21)
Now, the following sequence of inequalities holds true:

k < k f 0 ≤ 1 f 0 exp k B - -1 B - (22)
and we also have, using the concavity of H in eq. ( 16):

1 f 0 exp k B - = H B H ( ) ≤ 1 + H ( ) H ( ) B - (23)
The two inequalities ( 22) and ( 23) imply that for all for which eq. ( 16) has solutions there holds the inequality:

H ( ) -kH( ) > 0 (24)
Since B (0) = 0 and the functions involved in eq. ( 16) are continuous, the same inequality will hold true for all B in some maximal interval I B := [0, B F ), with B F > 0. This interval is bounded, for inequality (23) holds true only for a bounded set of values of B . Hence we also have B < +∞. Using the above in eq. ( 21) it follows that

d B d ( ) > 0, (∀) B ∈ [0, F ) ( 25
)
where F is the equivalent strain in the bulk of the sheet corresponding to B F . Hence the function B ( ) is monotonically strictly increasing. According to our definition of the interval I B , for B F there must hold eq. ( 15). Finally, since F < B F , the left-hand member in eq. ( 21) is strictly positive for F , and then

d B d ( F ) = +∞ ( 26
)
Thus F is the equivalent failure strain in the bulk of the sheet.

As a first remark, in the above proof we have also shown that eq. ( 14) is indeed a failure criterion: the nature of the failure point on the graph B = B ( ) is cuspidal.

A second remark is in order here: the failure strain in zone B does not depend on the initial imperfection ratio f 0 . It is an intrinsic property of the material. We have then the following result:

Corollary 1.1 In the conditions of Theorem 1, given a (target) value T x , there exists f 0 < 1 such that for the predicted MK-strain x we have x = T

x , if and only if the failure strain in zone B satisfies B

x > T x . Proof. For a given B F satisfying eq. ( 15), the function -→ φ( ) = H( ) exp(-k ) is strictly increasing on the interval [0, B F ] and strictly decreasing on [ B F , +∞). For an exact match of the target strain we must have

f 0 = φ T x /k /φ B F < 1, that is, T x < k B x = B x .
Most materials encountered in practice satisfy the condition of the above Corollary. As a first example, let us consider the useful power law for which the calculation is explicit.

Corollary 1.2 In the conditions of Theorem 1, for a hardening law in the form H ( ) = K ( 0 + ) n , there exists an initial imperfection ratio f 0 < 1 such that MK-predicted plane strain limit strain matches the target strain T x if and only if the following inequality holds true:

T x < n -k 0 (27)
Proof. The solution of eq. ( 15) is in this case

B F = n/k -0 .
We remark that the influence of the yield function, through the parameter k, is very small, since 0 is usually close to zero. For example, for the AISI 304 steel described in [START_REF] Campos | Theoretical and experimental determination of the forming limit diagram for the AISI 304 stainless steel[END_REF] the hardening law was given in the form required by Corollary 1.2, with n = 0.47 and 0 = 0.01 (a complete strain rate dependent form is discussed later). The experimentally measured plane strain limit strain, the target strain, was T x = 0.32; after identifying the yield surface from directional and biaxial data, the parameter k is identified as k = 0.9; the MK failure strain in zone B is in this case B x = nk 0 = 0.4610 and hence satisfies eq. ( 27). The initial imperfection ratio can then be determined as (with φ defined in Corollary 1.1) f 0 = φ T x /k /φ B F = 0.9737 As a second example, we consider the aluminum alloy AA3104-H19 described in [START_REF] Wu | Evaluation of anisotropic yield functions for aluminum sheets[END_REF]. For this alloy the cited work employs a power hardening law with n = 0.07 and 0 = 0. The target strain is in this case T x = 0.042 < B x = n. The MK model can then match exactly its prediction with the target strain T

x for an initial imperfection

f 0 = T x /n n exp n -T x = 0.
9923, precisely the value identified in the cited work through a trial and error procedure (0.992). Noteworthy, the estimation B x = n has also been obtained in Hutchinson and Neale (1978a) based on a combined perturbation MK-analysis. However, not every material satisfies the condition of Corollary 1.1. The aluminum alloy AA5182-O is described in [START_REF] Wu | Evaluation of anisotropic yield functions for aluminum sheets[END_REF] by the following hardening law in incremental form:

dσ = (1 -σ/K 1 ) n d (28)
After integration it takes the form:

σ = H( ) = K 1 1 -(1 -σ 0 /K 1 ) 1-n + K 2 (n -1) /K 1 1/(1-n) ( 29
)
with σ 0 = H(0). The function H is identified as usual by matching it to the results of a uniaxial traction test along the rolling direction, and then σ 0 is the initial yield stress in this direction. To compute the failure strain in zone B we solve eq. ( 15) which in this case reduces to

g(a) := a n + Kk(a -1) = 0 (30)
where we have denoted

K := K 1 /K 2 , and a 1-n := (1 -σ 0 /K 1 ) 1-n + K 2 (n -1) /K 1 .
The function g above is strictly increasing and convex. Additionally, g(0) = -K < 0 and g(1) = 1 > 0, and hence the above equation has a unique solution a ∈ (0, 1). It can be found numerically with the Newton algorithm. For σ 0 = 120.0 MPa, K 1 = 385.0 MPa, K 2 = 4100.0 MPa, n = 1.2 and k = 0.9383, one obtains a = 0.1189, value which substituted into

B x = k B = a 1-n -(1 -σ 0 /K 1 ) 1-n Kk (n -1) (31)
leads to the result: B x = 0.1998. The target strain reported in Wu et al. ( 2003) is T x = 0.22. This explains why in the cited work it was not possible in the case of the AA5182-O alloy to fit the MK-prediction to the measured data: according to Corollary 1.1, there exits no f 0 < 1 to allow this exact match. Although no explicit formulas are available for the pure tension case, it can be shown, by using the classical trial and error procedure, that a similar phenomenon takes place if one tries to fit the MK-predicted limit strain in pure tension (stress ratio t = 0) to the experimentally measured value: the later is greater than the MK-prediction for any f 0 < 1. Hence differences in hardening properties between pure tension and plane strain cannot offer a complete explanation either. Recent calculations, [START_REF] Signorelli | Predictions of forming limit diagrams using a rate-dependent polycrystal self-consistent plasticity model[END_REF], suggest that consideration of anisotropic hardening, by incorporating texture effects into the constitutive model, may lead to a better MK-prediction of the FLD of the AA5182-O alloy. While texture effects have been shown to influence considerably the MK-prediction of the FLD, e.g., [START_REF] Wu | On crystal plasticity FLD analysis[END_REF], these effects alone cannot explain entirely the mechanism by which a sheet metal fails since hardening is a material property whose basic features are decided at crystal level. In particular, it is well known that the hardening behavior of AA5182-O is strongly influenced by a particular form of interaction between dislocations and solute atoms, [START_REF] Picu | Strain rate sensitivity of the commercial aluminum alloy AA5182-O[END_REF], phenomenon reflected at macroscale by the serrated aspect of the hardening curve. This serrated aspect is featured by the experimentally measured hardening curves reported in [START_REF] Wu | Evaluation of anisotropic yield functions for aluminum sheets[END_REF] and [START_REF] Picu | Strain rate sensitivity of the commercial aluminum alloy AA5182-O[END_REF]. It seems then appropriate to include here a discussion on strain rate sensitivity for the plane strain case of the MK model. While closed form equations for the failure strains seem no longer possible, even for the simplest case considered here, it will be shown that including strain rate dependence in the hardening law of AA5182-O solves the above problem of the predictability of the plane strain limit strain (although not in a coherent manner). The discussion here parallels that in Hutchinson and Neale (1978b), where similar conclusions are drawn. We shall account for strain rate dependence by extending the hardening law to:

σ = H( ) ˙ m (32)
By a similar argument that lead to eq. ( 16) one obtains: with m > 0. We require H(0) 0. We have now an ordinary differential equation describing uniquely the evolution of B = B ( ), given the additional initial condition B (0) = 0 (for increased precision one could also consider a strictly positive initial value for B (0), since zone B is thinner and hence yields earlier than the rest of the sheet; however, this correction has negligible effects upon the overall results). With this new constitutive model, there appears now a remarkable difference in comparison with the rate independent case: the failure criterion ( 14) is fulfilled for an infinite equivalent strain in zone B, or equivalently, the equivalent strain at failure F in the zone outside the groove is a vertical asymptote for the solution B = B ( ). This can be shown by employing the equation for the inverse function = B , with the initial condition (0) = 0. Then, by repeated use of the Cauchy-Lipschitz theorem of existence and uniqueness and the concavity of H, the solution = B can be extended to the maximal interval [0, +∞).

d B d =         G( ) f 0 G B         1/m , with G(x) := H(x) exp(kx) (33)

Since =

B is strictly increasing and (+∞) = 0, it has an upper horizontal asymptote F < +∞. Next, since [0, F ) is the maximal interval on which the solution B = B ( ) can be extended, there is a one to one relationship between the parameter f 0 and F . Hence given a target strain T x it uniquely determines the parameter f 0 that matches it through the MK-prediction. Since for f 0 = 1 the solution of eq. ( 33) with B (0) = 0 is B ( ) = , defined over the interval [0, +∞), by the continuity of the solution of (33) with respect to the parameter f 0 , for a sequence ( f (k) 0 ) k -→ 1 there exists a corresponding sequence of maximal intervals [0, (k) F ) such that (k) F -→ +∞. We conclude that in the rate-dependent case for every target strain T x > 0 there exists an initial imperfection f 0 such that the MK-prediction of the plane strain limit strain is T

x . In the conditions imposed by the shape of the hardening curve (strictly concave) and f 0 < 1, the behavior of the solution of eq. ( 33) determined by B (0) = 0 is quite typical: strictly increasing convex function with explosive growth near the right-end F . We illustrate it for two materials. The first is the AISI304 steel for which the complete description in [START_REF] Campos | Theoretical and experimental determination of the forming limit diagram for the AISI 304 stainless steel[END_REF] was in the form of eq. ( 32) with H the power law with material parameters detailed earlier and m = 0.012. Then, with a trial end error procedure (and with a numerical Runge-Kutta o.d.e. integrator) we find that for f 0 = 0.950 the MK-prediction equals the target strain T x = 0.32, close to the value of f 0 = 0.955 employed in the cited work. The graph of the B = B ( ) function is shown in Fig. 2. The second application is for AA5182-O. Following data reported in Fig. 2 of [START_REF] Picu | Strain rate sensitivity of the commercial aluminum alloy AA5182-O[END_REF], the strain rate exponent for this material is approximately m = -0.0075. While a negative m is perfectly legit in eq.( 32), the differential equation ( 33) is no longer applicable in this case since it becomes identical with the equation for the inverse function = B , for which the (0) = 0 solution can be extended ad infinitum. Furthermore, the constitutive formula (32) may be satisfactory for constant strain rate loading histories but it becomes inadequate for loading histories with variable strain rates: the strain rate in zone B is strictly increasing and then eq.( 32) predicts a decrease of the yield stress in zone B, for m < 0, leading further to the loss of force equilibrium between zone B and the rest of the sheet. This remark applies equally well for positive m-exponents also, since in this case eq.( 32) may predict a too stiff stress-strain response.

In absence of relevant experimental data we refrain from further considerations. Then, just for the purpose of illustration of eq.( 33), we shall proceed against experimental evidence and take m = 0.0075 [START_REF] Signorelli | Predictions of forming limit diagrams using a rate-dependent polycrystal self-consistent plasticity model[END_REF] used m = 0.02). Readjusting the material parameters to fit the same hardening curve as in [START_REF] Wu | Evaluation of anisotropic yield functions for aluminum sheets[END_REF] but in the strain rate dependent form of eq. ( 32) with H given by eq.( 29) and ˙ = 0.001/s we get: σ 0 = 126.4 MPa, K 1 = 398.5 MPa, K 2 = 4100 MPa, n = 1.1. Illustrated in Fig. 2 is the B = B ( ) evolution which for f 0 = 0.998 fits exactly the target strain T x = 0.22. However, besides the constitutive inconsistency signaled here, there is another troubling consequence of eq. ( 33): the thickness strain in zone B is infinite at failure. The sheet fails when the thickness in zone B has become zero, and hence eq. ( 14) is in this case equivalent to a fracture criterion.

About the relationship between the Hill-Swift and MK failure criteria and a simple approximation for negative strain ratios

Eq. ( 15) is nothing but Hill, or Swift criterion when ρ = 0, [START_REF] Hill | On discontinuous plastic states, with special reference to localized necking in thin sheets[END_REF], [START_REF] Swift | Plastic instability under plane stress[END_REF]. Thus in this case Hill's failure criterion is equivalent with the MK failure criterion. It is then natural to inquire further about the relationship between these criteria for the whole range of strain ratios. First, let us make a simple remark about the structure of the stress tensor in zone B. According to the equilibrium equation ( 12) its components σ B n and σ B nt are proportional with the corresponding components of the stress tensor σ in the bulk of the sheet. To discern completely the contribution of the stress state in the bulk of the sheet to that in zone B we write the σ B t component in the form σ B t = σ A t + q / f h , where q is the supplementary stress component in the t-direction. Then with respect to the xy coordinates we have:

σ B xy = 1 f h        [σ] xy + q         n 2 2 -n 1 n 2 -n 1 n 2 n 2 1                = 1 f h (σ + qT ) ( 34
)
with an obvious definition for the matrix T . We state next a generalization of Theorem 1 that will be needed in what follows. Its proof follows similar ideas and therefore we defer it to an Appendix.

Theorem 1.b. Assuming a strain rate independent hardening law, then for every stress ratio t ∈ [0, 1] and every initial orientation of zone B the equivalent strain at failure in zone B, defined by eq. ( 14), is finite: B F < +∞.

Returning to the investigation of this section, we have: Theorem 2. For a strain rate independent hardening law, in zone B the following implication holds true

d B d = +∞ =⇒ dσ B n d B n = σ B n , dσ B nt d B n = σ B nt , dσ B t d B n -σ B t = 1 f h dq d B n ( 35
)
Proof. Let us consider zone B oriented at an arbitrary angle as in Fig. 1. Differentiating the first of the equilibrium equations in ( 12) with respect to B , the flow "time" within the groove, and since, according to Theorem 1.b, the failure strains in zone B are finite, the following equality must hold at the failure moment:

dσ B n d d d B = d d B f h σ B n =⇒ d d B f h σ B n B = B F = 0 ( 36
)
Differentiating eq. ( 9):

d f h d B = f h d d B B z -z (37)
From incompressibility and compatibility, eq. ( 11), we have:

d B z -d z = -d B n + d B t + (d n + d t ) = -d B n + d n (38)
and hence

d d B B z -z B = B F = - d B n d B B = B F ( 39
)
where above we have used the flow rule and the MK failure criterion:

d A n d = n 2 1 + ρn 2 2 ∂ f ∂σ x (1, t) < +∞ =⇒ d n d B B = B F = 0 (40)
Returning then to eq.(36), after performing the limit B -→ B F (to legitimate the use of the chain rule for derivatives), we have:

f h - d B n d B σ B n + dσ B n d B B = B F = 0 ⇐⇒ f h -σ B n + dσ B n d B n B = B F = 0 (41)
Simplifying above with f h (which is nonzero since the failure strain B is finite) gets the first implication in (35). The second is obtained similarly by using the second equilibrium equation in ( 12) while the third by using eq.( 34).

The first of the implications in eq. ( 35) becomes a Swift type criterion if the groove is oriented along the transverse axis: dσ B x /d B x = σ B x . However, this failure criterion holds for the stressstrain state in zone B. Applying it requires knowledge of the entire stress-strain history, as in the MK model, and as remarked in the early work of [START_REF] Swift | Plastic instability under plane stress[END_REF]. Instead, the current practice features Swift criterion in combination with the approximation of the stress-strain state in the neck with the one outside it. To recover a Hill type of failure criterion we need the following: Corollary 2.1 For a strain rate independent hardening law, the following implication holds true

d B d = +∞ =⇒ dσ B n d B z = -σ B n ( 42
)
Proof. From Theorem 1.b we know that the failure moment is cuspidal. Hence it can be approached only with a decreasing to zero sequence of increments d . Since every corresponding d B increment is always finite, from the compatibility equation ( 11) it follows that at failure there holds: d B t = 0. Then from incompressibility we also have

d B z = -d B n .
The conclusion follows from Theorem 2. Equation ( 42) is Hill's failure criterion written for the stress-strain state in zone B. If one makes the approximation σ B = σ, then with σ n = (1 + t tan 2 ψ)σ x /n 2 1 , and since the neck angle is stationary at the failure instant, one obtains the more familiar form dσ x /d z = -σ x . In what follows we shall weaken this assumption to obtain a better approximation for the MK-predicted limit strains in the range of negative strain ratios. We assume:

σ B = 1 f h σ (43)
As a consequence, we have t B = t: the stress ratio in zone B is equal to the one in the bulk of sheet. For negative strain ratios this assumption is close to the actual stress path predicted by the MK model. For example, for the previously described AISI 304 sheet in pure tension, t = 0, the MK model predicts at failure t B = (σ y + qn 2 1 )/(σ x + qn 2 2 ) ≈ 0.08, while for most of the loading process it satisfies: t B < 0.04. One also neglects the shear component σ B xy . For example, for the same pure tension experiment with AISI 304, the MK model predicts σ B xy /σ A x < 0.06 for the whole loading process. Since the stress ratio is the same and constant in both areas of the sheet we have

d z = -(d x + d y ) = -(1 + ρ)d x = -(1 + ρ) f (1, t) 1 + tρ (44)
with a similar formula for d B z . Hence assumption (43) has also the consequence that the thinning rate in both areas is the same. The initial conditions are different: for zone B we have f 0 < 1, determined from the plane strain case discussed earlier, while outside zone B we have f 0 = 1. Then with a reasoning similar to the one followed in the proof of Theorem 1, from (43) the equations for the failure strains are now

H B F -kH B F = 0, with: k := 1 + ρ 1 + tρ f (1, t) (45) T exp(k F ) -H( F ) = 0, with: T := f 0 H B F exp -k B F (46) x = f (1, t) 1 + tρ F , y = ρ x (47)
As illustrated in the next section, the above equations improve considerably the Hill approximation to the failure strains for negative strain ratios, coinciding with the MK prediction for ρ = 0. Therefore we shall refer to the above formulas as the HMK approximation. Unfortunately, such simple formulas are no longer applicable for positive strain ratios, since in this case the stress ratio in zone B varies considerably during the loading process and hence assumption (43) is no longer valid. However, as shown next, for positive strain ratios there is another approach available, if simple approximations are desirable, based on another explicit formula.

The non-hardening case (positive strain ratios)

The complicated nature of the general equations of the MK model makes almost impossible their explicit integration. Explicit solutions, even in particular cases, are of great help in discerning the role of each constitutive ingredient and the degree to which it influences the MK-predicted limit strains. The non-hardening case is amenable to explicit integration and its solution brings further support for the importance of accurate modeling of the biaxial yield curve when predicting limit strains with the MK model, [START_REF] Bassani | On the prediction of necking in anisotropic sheets[END_REF], [START_REF] Barlat | Crystallographic texture, anisotropic yield surfaces and forming limit diagrams[END_REF], [START_REF] Lian | Plastic behavior and stretchability of sheet metals. Part II: effect of yield surface shape on sheet forming limit[END_REF]. This case was first analyzed in qualitative terms by [START_REF] Sowerby | Failure in sheet metal in biaxial tension[END_REF]. Here we make more precise statements. We have:

Theorem 3. Let us assume zone B oriented along the transverse symmetry axis, and positive strain ratios. Then for a perfectly rigid-plastic material, and a given f 0 < 1, there exists ρ * > 0 such that for strain ratios ρ ∈ [0, ρ * ] the material fails at the beginning of the loading process.

For ρ > ρ * the limit strains are:

x = t B 0 t PS f y (1, τ) 2 dτ f (1, τ) ρ f x (1, τ) -f y (1, τ) , y = ρ x ( 48
)
where f x and f y denote the partial derivatives ∂ f /∂σ x and ∂ f /∂σ y , respectively, t PS is the stress ratio corresponding to ρ = 0, and t B 0 > t PS is the initial stress ratio in zone B, determined from the condition f (1, t B 0 ) = f 0 f (1, t ρ ), with t ρ the stress ratio in the bulk of the sheet. Proof. If the constant yielding limit is denoted H 0 , yielding of the entire sheet requires that

σ x f h f (1, t B ) = H 0 = σ x f (1, t ρ ) =⇒ f 0 exp( B z -z ) = f (1, t B ) f (1, t ρ ) (49)
The above equality cannot be satisfied for any strain ratio ρ ≥ 0. For example, for ρ = 0 the left hand member is strictly less than one, while the right hand member is equal to one. There exists then a maximal interval [0, ρ * ] of strain ratios for which the sheet necks as soon as the loading begins. ρ * corresponds to the stress ratio t * defined by

f 0 = f (1, t PS ) f (1, t * ) and then ρ * = f y f x (t * ) (50)
As before, we assume that yielding is simultaneous in both areas of the sheet and hence neglect the strains accumulated in zone B while the bulk of the sheet was rigid. Next, we assume ρ > ρ * . The compatibility equation enforces the constraint:

d B y = d y ⇐⇒ f y (1, t B )d B = f y (1, t)d (51)
Adding the incompressibility condition we get

B z -z = - B 0 f x (1, t B ) - 1 ρ f y (1, t B ) de ( 52
)
where B is the current accumulated plastic strain in zone B and t B = t B (e) is the evolution of the stress ratio in the groove for e ∈ [0, B ]. Substituting the above formula into eq. ( 49) we obtain the following integral equation for the function [0, B ] e -→ t B (e) for any accumulated plastic strain B :

B 0 f x (1, t B ) - 1 ρ f y (1, t B ) de = ln f 0 f (1, t ρ ) f (1, t B ( B )) ( 53
)
Differentiating above with respect to B yields the following differential equation with separated variables:

dt B d B = f (1, t B ) 1 ρ - f x f y (1, t B ) ⇐⇒ d B dt B = ρ f y (1, t B ) f (1, t B ) f y (1, t B ) -ρ f x (1, t B ) ( 54
)
which now defines the function B = B (τ), for τ ∈ [t B , t B 0 ], with t B 0 the initial stress ratio in zone B. It is computed from the yielding condition at the beginning of the loading process:

f (1, t B 0 ) = f 0 f (1, t ρ ).
With the initial condition satisfying t B 0 < t ρ , it follows that the solution t B ( B ) of the first equation above is strictly decreasing (an information already contained in the compatibility equation, as remarked in the analysis of [START_REF] Sowerby | Failure in sheet metal in biaxial tension[END_REF]). Its maximal range is [t PS , t B 0 ], where t PS is the stress ratio in zone B at the failure moment. We denote

B F := B (t PS ).
We are now in a position that allows us calculate the limit strains in the bulk of the sheet:

d B z -d z = -d B x + d x =⇒ x = B x = B F + ( B z -z ) = B F (55)
With the change of variable ( 54) we have:

B x = B F = B F 0 f x (1, t B )d B = t B 0 t PS ρ f x (1, τ) f y (1, τ)dτ f (1, τ) ρ f x (1, τ) -f y (1, τ) (56)
( B z -z ) = B F = - B F 0 f x (1, t B ) - 1 ρ f y (1, t B ) d B = - t B 0 t PS f y (1, τ) f (1, τ) dτ ( 57
)
Substituting the last two formulas into eq. ( 55) leads to the result in (48).

To test formula (48) against actual experimental data we incorporate hardening by making the following approximation of the limit strains when ρ ≥ 0:

x = PS x , y = ρ x , for ρ ∈ [0, ρ * ] ( 58
) x = PS x + t B 0 t PS f y (1, τ) 2 dτ f (1, τ) ρ f x (1, τ) -f y (1, τ) , y = ρ x , for ρ > ρ * (59)
where PS x is the plane strain limit strain as predicted with the HMK formulas (45-47). The same HMK approximation is used for ρ < 0.

Besides the hardening law, the above formulas require only the description of the biaxial curve in the first quadrant (σ x , σ y ). We shall use an in-plane isotropic sixth order polynomial (since only half of that quadrant is actually needed):

f (σ x , σ y ) 6 = a 1 σ 6 x + σ 6 y + a 2 σ 5 x σ y + σ x σ 5 y + a 3 σ 4 x σ 2 y + σ 2 x σ 4 y + a 4 σ 3 x σ 3 y ( 60
)
Figure 3: Influence of the shape parameter P upon the biaxial yield curve of AISI 304.

A detailed description of the identification procedure for the four parameters a i can be found in [START_REF] Soare | A four parameter in-plane isotropic yield function[END_REF]. Four data points are required: the yielding stress and the r-value along the rolling direction, σ 0 and r 0 , the balanced biaxial yield stress, σ b , and another yielding stress along a constant stress ratio path. For the purpose of illustration we replace the fourth data point with a shape parameter P and illustrate the effect of its variation upon the biaxial curve in the case of the AISI 304 steel, Fig. 3, where σ 0 = σ b , and r 0 = 0.92, data from [START_REF] Campos | Theoretical and experimental determination of the forming limit diagram for the AISI 304 stainless steel[END_REF]. The shape parameter P is similar to the parameter P = σ PS /σ b introduced by [START_REF] Barlat | Crystallographic texture, anisotropic yield surfaces and forming limit diagrams[END_REF], and it can be used to control the shape of the biaxial yield curve. The results of the HMK and non-hardening approximations for the MK-predictions in the case of AISI 304 are featured in Fig. 4. For this material the plateaus defined by eq. ( 58) are quite large: for P = 0.928 we have ρ * = 0.36, for P = 0.900 we have ρ * = 0.28, and for P = 0.875 we have ρ * = 0.3. This is caused by the large initial imperfection (i.e., f 0 small) required by this steel. For aluminum alloys ρ * is much closer to zero, since f 0 is much closer to 1.0. Besides this restriction in applicability, formula (59) captures very well the sensitivity of the MK model to the shape of the yield surface. For comparison, the actual MK-predictions are shown in Fig. 5. The AISI 304 steel turns out to be a tough test for these approximations. For smaller initial imperfections and shorter hardening times they are much closer to the MK-predictions.

An MK-variant: material heterogeneity

The explanations for the diverse failure mechanisms of metals must ultimately rely on their physical structure. At crystal level: motion of dislocations, their interaction with inclusions, the penetration of dislocations through grain boundaries; at polycrystal level: the spatial orientation and rotation of the grains (the texture), and also the effect of the apparition and evolution of additional substructures like voids. One has to recognize that these complex features manifest at continuum level not only through constitutive behavior but also through heterogeneity. In a work on strain localization in single crystals [START_REF] Asaro | Strain localization in ductile single crystals[END_REF] refer to the work of [START_REF] Price | Deformation of age-hardened aluminum alloy crystals-II[END_REF] where the following interesting behavior is observed. If after stretching a single crystal strip just enough to generate step bands (which if strained further will develop into necking bands), one unloads and polishes away these bands and then reloads the strip, then new bands will form after enough straining but never at the previous locations. Continued straining will eventually result in failure of the material at these new locations. Although in [START_REF] Asaro | Strain localization in ductile single crystals[END_REF] necessary conditions for localization are sought only at constitutive level, the above experiment can be clearly interpreted as a manifestation of spatial heterogeneity of the material. At one scale step above, an interesting study showing how a non-uniform deformation regime appears naturally in textured polycrystals can be found in the paper of [START_REF] Harren | Nonuniform deformations in polycrystals and aspects of the validity of the Taylor model[END_REF].

Admitting then that the material (plastic) properties may feature spatial variations across the sheet, we ask further what is the amplitude of this variation so that it be a plausible candidate for explaining sheet failure. In other words, can we estimate/explain the limit strains observed in practice based on an assumption of material heterogeneity ? Let us consider a thin sheet of uniform thickness modeled as a rigid-plastic material with a hardening law in the form

H A ( ) = H( ) = K ( 0 + ) n .
Next, let us assume there exists a variation in the parameters of this hardening law along a narrow band of material spread across the sheet as illustrated in Fig. 6. This variation is assumed in the form

H B ( ) = K B 0 + n+δ , with B 0 = ( 0 ) n/(n+δ) (61)
in zone B, the band. A similar variation was considered in Aretz (2007b). Next, for simplicity we consider here only the plane strain strain case: the sheet is monotonically bi-axially loaded along its symmetry axes x and y such that d A y /d A x = 0. The superscript identifies the zone of the sheet where a variable/function is evaluated. Above, B 0 was chosen so that H B (0) = H A (0). Then assuming the yield surface is the same at all material points, plastic loading is simultaneous everywhere in the sheet. Furthermore, at the boundary between the two zones the compatibility of deformation requires d B y = d A y and hence, from the constitutive law, the stress ratio t := σ y /σ x is the same in both zones: t = t PS , with t PS the plane strain stress ratio. Hence the yielding condition in the two zones reads:

σ B x f (1, t) = H B ( B ) and σ A x f (1, t) = H A ( A ).
Since the two zones of the sheet work-harden at (slightly) different rates, a thickness heterogeneity will develop in zone B (assuming the material here more ductile) and as in the MK model we define the thickness ratio f h :

f h (t) = h B (t) h A (t) =⇒ ḟh = f h (˙ B z -˙ A z ) =⇒ f h = exp( B z -A z ) (62)
Note that initially h B = h = h A and hence f h (0) = f 0 = 1. Finally, equilibrium at the interface between the two zones requires that we have σ B x = σ A x / f h . From this relationship, the yielding condition in the two zones, and eq. ( 62) we can write:

H A ( A ) exp(-k A ) = H B ( B ) exp(-k B ) (63)
where we denoted k := f (1, t PS ). This relationship holds true all along the loading process. Differentiating, the failure criterion d B /d A = +∞ requires that at failure the following relation holds true:

H B ( B F ) -kH B ( B F ) = 0 (64)
with the subscript F indicating the failure moment. Equations ( 63) and (64) determine the equivalent plastic strains in the two zones at failure. The failure strain is then A x = k A F . The AISI 304 steel sheet fails in plane strain at x = 0.32. With 0 = 0.001, n = 0.47 we find that a variation in the hardening exponent of δ ≈ 0.054 leads to failure at the mentioned target strain.

The aluminum alloy AA3104-H19 sheet failed at x = 0.04. With 0 = 0 and n = 0.07, a variation of δ ≈ 0.003 leads to failure at the target strain. The reasoning can be repeated for the AA5182-O alloy, although here we will not vary the hardening exponent. With the hardening law described by formula (29), the parameters σ 0 = 120 MPa, K 1 = 385 MPa, K 2 = 4100 MPa, n = 1.2 are assumed for the bulk of the sheet, and the following variations, σ 0 = 120 MPa, K 1 = 395 MPa, K 2 = 3703 MPa, n = 1.2, for zone B. In these conditions the sheet will fail precisely at the target value T x = 0.22. In Fig. 7 are shown the two hardening curves for AA5182-O. Also shown is the set of experimental data points, digitized from [START_REF] Wu | Evaluation of anisotropic yield functions for aluminum sheets[END_REF]. The hardening curve in zone B is within the area covered by the experimental data set. Based on the above idealized analysis, we can then conclude that small variations in the plastic properties across the sheet, undetectable by macroscopic experimental techniques in the case of AA5182-O, can trigger unstable plastic flow and ultimately lead to failure at strain levels comparable to those observed in practice.

An extension of the MK-model to sheets of uniform thickness

The analysis in the previous section reveals an alternative to the thickness heterogeneity of the classical MK model: one can induce discontinuities in the thickness strain rate, and hence inhomogeneous plastic flow, by assuming a material heterogeneity. We aim next at extending the MK model so that unstable plastic flow is no longer initiated by a thickness imperfection but by a material heterogeneity. However, unlike in the previous section, and for more generality, at constitutive level we shall still work with a homogeneous material. The presence of the spatial material heterogeneity will be modeled by introducing an additional term to the global rate of deformation, an eigenstrain-rate (in the sense defined in [START_REF] Mura | Micromechanics of defects in solids[END_REF]). Let us consider a thin metallic sheet under in-plane biaxial stretching. As before, the strain path is assumed proportional (however, the following arguments are general and can be easily extended to nonlinear loading paths). At a certain point during the loading history a neck will develop across the sheet. We assume its geometry the same as in the MK model, that is, a gradual developing narrow groove in the thickness of the sheet, Fig. 1, oriented at the moment t at some angle ψ(t) with respect to the symmetry axes. The groove itself will be identified as zone B. Then, if we denote h(t) and h B (t) the current thicknesses outside zone B and in zone B, respectively, we can define the thickness imperfection factor:

f h (t) = h B (t) h(t) (65)
We assume that until the moment of groove initiation, call it t 0 , the thickness of the sheet was uniform. In particular, we have f 0 = f h (t 0 ) = 1. After groove initiation (which in this initial stage has the features of a diffuse neck, as opposed to the sharp neck characteristic of the final stages of failure) zone B is characterized by a faster thinning rate than the rest of the sheet. Differentiating above we have:

ḟh = f h ḣB h B - ḣ h = f h D B z -D z ( 66
)
where D B and D denote the rate of deformation tensors in zone B and outside it, respectively (as in the classical MK model, the strains are assumed spatially uniform in both areas of the sheet). According to our assumption, at the moment t 0 we have h B = h; however, groove initiation at the moment t 0 requires that the rate ḟh (t 0 ) be nonzero, and more precisely, it must be strictly negative. Thus we shall assume that there exists a strictly positive constant p such that:

D B z (t 0 ) -D z (t 0 ) = -p (67)
Since the material of the sheet is incompressible, it follows that there must be a strain rate jump in the in-plane strains also. Since the tangential strain t along the groove is constrained by the boundary sharing with the bulk of the sheet, it follows that the in-plane strain rate jump is restricted to the normal strain component n . Finally, assuming the jump p propagates during the entire loading history after the t 0 moment, the strain rate in the groove has the form:

D B = p (n ⊗ n -e z ⊗ e z ) + D B ( 68
)
where e z is the unit vector in the normal to the sheet direction, and D B is the usual rate of deformation tensor in zone B generated by the stressing at the boundary of zone B, and hence determined by the constitutive law and compatibility at the boundary, to be discussed next.

Depending on the phenomenology of the failure event to be modeled, a more sophisticated theory may take into account only a local effect, in time, or a more complex time-evolution of the perturbation p. However, our assumption that the constant jump p propagates during the entire loading history is consistent with the presence of a material heterogeneity. Then the initiation moment can be safely considered at the starting moment of the loading process: t 0 = 0. Then returning to eq. ( 66), we have for the imperfection ratio:

ḟh = f h -p + D B z -D z =⇒ f h (t) = exp -pt + B z -z (69)
Having established conditions that generate inhomogeneous plastic flow across the sheet, we can next proceed along the same lines as in the classical MK model to establish equations for the evolution of the stress-strain state in the two areas of the sheet. The only addition to the general theory in Section 2 is the consideration of the strain rate into the hardening law.

The description is now incremental, with emphasis on the numerical solution of the equations of the model. Thus, in the bulk of the sheet constant strain increments d x are successively added, at a constant strain rate D x = ˙ x , with the time increment dt = d x /˙ x . As for the classical MK model, accuracy demands very small strain increments: d x ≤ 10 -4 . From work equivalence, for a given strain ratio ρ the equivalent plastic strain rate in the bulk of the sheet is

˙ = 1 + t ρ ρ f (1, t ρ ) ˙ x ( 70
)
where above t ρ denotes the stress ratio in the bulk of the sheet. Next, from equilibrium and compatibility at the boundary between the two zones, from the yielding condition in the two zones, and making use of the representation (34) for the stress tensor in zone B, the following system of equations determines the accumulated plastic strain state ˙ B and the stress state σ B = (σ + qT )/ f h , in zone B, for the current strain increment in the bulk of the sheet:

f (σ + qT ) = f h H B + ˙ B dt, ˙ B (71) T : ∂ f ∂σ (σ + qT ) ˙ B = n 2 2 + ρn 2 1 ∂ f ∂σ x (1, t ρ ) ˙ (72)
The rotation of zone B follows the same law (10) as for the original MK model. The above system is solved for the unknowns q and ˙ B with the Newton-Raphson algorithm. At each increment, as starter point (initial guess) one can take the solutions at the previous iteration. For the very first increment the starter point can be (q, ˙ B ) = (0, ˙ ). To ease the computational burden, one can use in the above system the imperfection f h at the previous increment (with f h = exp(-pdt) for the very first increment). Note, however, that for very small values of the strain rate jump p, as in the example featured next, the variation of f h during one increment is no longer negligible and one has to either enhance the Newton-Raphson algorithm with a line search or to work with the full Jacobian of the above system. The loading process ends when the following (numerical) failure criterion is met:

˙ B ˙ ≥ N f (73)
with N f a big enough number (depending on the magnitude of the strain increment in the bulk of the sheet; for 10 -4 ≤ d x ≤ 10 -5 , N f = 10 allows a satisfactory closeness to the failure strains). The calculations are repeated for all possible initial orientations of zone B, ψ 0 ∈ [0 o , 90 o ] (based on information from the previous failure angles, this range can be considerably reduced), and then define the limit strain for a given strain ratio ρ as the smallest limit strain:

x (ρ) = Min { x (ψ 0 , ρ)| ψ 0 ∈ [0 o , 90 o ]} , y (ρ) = ρ x (ρ) (74)
Instead of the thickness imperfection parameter of the classical MK model, the parameter of the present model, which we shall term the pMK model, is the strain rate jump p. This parameter is determined to fit an experimentally measured limit strain. Similarly to the MK model, we take here the plane strain limit strain as data point. Assuming in this case, for simplicity, zone B along the transverse direction, with arguments similar to those in Section 3, it can be shown that eqs. (71-72) reduce to the following ordinary differential equation:

˙ B = ˙         exp(pt) G( ) G B         1/m , with G(x) := H(x) exp(kx) (75)
where above, and in what follows, we assume a hardening law in the form (32), k := f (1, t PS ), and the constant ˙ is defined by eq. ( 70) as ˙ = ˙ x /k. With a trial and error approach, the parameter p can be tuned so that the failure time t F , when the criterion (73) is met, determines a corresponding strain in the bulk of the sheet, x = ˙ x t F , equal or close to the experimentally measured plane strain limit strain T x . We first illustrate the pMK model for the AA5182-O sheet. Applying the above fitting procedure for the p parameter with ˙ x = 10 -3 /s, gets us p = 1.5 × 10 -5 /s (with an exact match of the target strain 0.22; the hardening parameters: σ 0 = 126.4 MPa, K 1 = 398.5 MPa, K 2 = 4100 MPa, n = 1.1, m = 0.0075). The pMK-predictions for the entire range of strain ratios are shown in Fig. 8. For comparison, the predictions of the classical MK model are also shown there. Since no initial imperfection f 0 < 1 can fit the target plane strain limit strain, we used f 0 = 0.9995, a value also used in [START_REF] Wu | Evaluation of anisotropic yield functions for aluminum sheets[END_REF]. For both pMK and MK models the yield surface of the AA5182-O was described with the following plane stress sixth order polynomial function

σ 6 = a 1 σ 6 x + a 2 σ 5 x σ y + ... + a 8 σ 4 x + ... σ 2 xy + a 13 σ 2 x + ... σ 4 xy + a 16 σ 6 xy (76)
with the coefficients a i listed in Table 1, obtained through optimization by using the directional and biaxial data reported in [START_REF] Wu | Evaluation of anisotropic yield functions for aluminum sheets[END_REF]. The difference between the predictions of the two models is due to the strain rate sensitivity of the hardening law used in the pMK model (while with the MK model the rate independent formula (29) was used).

The second application is for AISI 304. Its yield surface is described with the same Poly6 function, restricted to an in-plane isotropic approximation, with the polynomial coefficients listed in Table 1 (the a 1 , ..., a 7 coefficients define the P = 0.900 biaxial curve in the (σ x , σ y ) plane used in Section 5). As for AA5182-O, in the bulk of the sheet we use the low (quasistatic) strain rate ˙ x = 0.001/s. With the hardening law H( ) = K ( 0 + ) n ˙ m and hardening parameters K = 1527 MPa, 0 = 0.01, n = 0.47, m = 0.012 the target plane strain limit strain e T x = 0.32 is matched for p = 1.7 × 10 -4 . The pMK predictions for the whole range of strain ratios are shown in Fig. 9, together with the predictions of the corresponding rate-independent MK model, H( ) = K ( 0 + ) n , with the same hardening parameters as above. The differences between the two predictions suggest that the strain rate sensitivity of AISI 304, at the low strain rate considered here, may not be as high as reported in [START_REF] Campos | Theoretical and experimental determination of the forming limit diagram for the AISI 304 stainless steel[END_REF]. Indeed, taking a smaller strain-rate sensitivity, m = 0.001, and retaining the same values for the rest of the hardening parameters, the plane strain limit strain is matched for p = 9.0 × 10 -5 ; the pMK predictions in this case are closer to the experimental data, Fig. 9. However, as discussed earlier, since the strain rate in zone B is variable, and abruptly increasing near the failure point, the simple hardening law employed here can be regarded only as a rough approximation to the actual hardening behavior -21.726 34.567 -23.923 9.8762 31.809 -23.714 33.169 26.411 a 1 a 2 a 3 a 4 a 5 a 6 a 7 a 8 1. 0000 -2.8750 7.1906 -9.6311 7.1906 -2.8750 -19.881 37.774 -19.881 8.8750 27.691 -39.906 27.691 31.672 where t is the stress ratio in the bulk of the sheet that corresponds to the strain ratio ρ. From incompressibility, flow rule and compatibility we have:

B z -z = - B 0 f x (1, t B ) de + f (1, t) 1 + tρ (A.
3) with f x denoting the partial derivative with respect to σ x . Eq. (A.2) holds true during the entire loading process. Then differentiating with respect to and using eq. (A.3) we obtain that the following equality must also hold true during the loading process:

H ( ) -KH( ) f (1, t) exp(K ) = f 0 H ( B ) - 1 f (1, t B ) d d B f (1, t B ) + f x (1, t B ) H( B ) f (1, t B ) exp(-B z ) d B d (A.4)
where for simplicity we denoted K := f (1, t)/(1 + tρ).

Next we assume that ρ 0 (when ρ = 0 the result follows from Theorem 1). Then, since f y (1, t) 0, the compatibility equation can be reformulated as an ordinary differential equation for the function = ( B ): .5) with t B = t B ( B ) a function defined by, say, eq. (A.2). We assume as initial condition | B =0 = 0. Then eq. (A.2) defines the initial value of t B as f (1, t B ) = f 0 f (1, t). With f 0 < 1 it follows: f y (1, t B ) < f y (1, t). Hence d d B (0) < 1 (A.6) By continuity, the above inequality holds true for a (positive) maximal neighborhood of zero, say [0, E). Then it follows that the function f (1, t B ) is strictly decreasing on this neighborhood. Indeed, for 1 < 2 , and the corresponding B 1 < B 2 (since = ( B ) is increasing), with 1 < B 1 and 2 < B 2 , we have using eqs. (A.2) and (A.3):

d d B = f y (1, B ) f y (1, t) (A
f (1, t B 2 ) f (1, t B 1 ) = H( 1 ) H( 2 ) H( B 2 ) H( B 1 ) exp        K( 2 -1 ) - B 2 B 1 f x (1, t B )d B        (A.7)
In the present constitutive context, the function [0, 1] τ -→ f (1, τ) has positive range, is convex, and has a global unique minimum at τ = t PS ∈ (0, 1), the plane strain stress ratio. Also, f x (1, τ) has positive range and is strictly decreasing, while f y (1, τ) is strictly increasing. Then f x (1, t B ) > K, since f (1, t B ) < f (1, t) and f x (1, t) = K. Then we have the inequality: .8) with the entity between the square brackets being positive. By taking 1 close to 2 , the product of the two H-ratios is close to one. Due to the concavity of H, for the same distance 2 -1 the decrease of the exponential term is faster, and hence f (1, t B 2 ) < f (1, t B 1 ). Then, by continuity, f (1, t B (E)) < f (1, t) and then d /d B (E) < 1. The above reasoning can be repeated and hence the solution of eq. (A.5) can be again extended. The maximal interval [0, E B ) As in the case with vertical groove, due to (A.14), the left-hand term of the above equality is strictly positive only over a finite interval. Finally, the counterpart of eq. (A.5), the compatibility equation, is:

f (1, t B 2 ) f (1, t B 1 ) < H( 1 ) H( 2 ) H( B 2 ) H( B 1 ) exp -K ( B 2 -2 ) -( B 1 -1 ) (A
d d B = ∂ f ∂σ (T A + rT ) • T n 2 2 f x (1, t) + n 2 1 f y (1, t) (A.16)
The reasoning can now proceed along the same lines as for the case with vertical groove.

Figure 2 :

 2 Figure 2: B = B ( ) graphs for the AISI304 and AA5182-O alloys in the rate-dependent case. The dashed line represents the B = line.

Figure 4 :

 4 Figure 4: Approximations of the MK-predictions: HMK on the left (ρ ≤ 0) and formulas (58-59) on the right (ρ > 0).

Figure 5 :

 5 Figure 5: Predictions of the MK model in the case of AISI 304.

Figure 6 :

 6 Figure 6: In Zone B the material properties feature small deviations from those in the rest of the sheet, zone A.

Figure 7 :

 7 Figure 7: The hardening curve(s) of AA5182-O in uniaxial traction along RD.

Figure 8 :

 8 Figure 8: Predictions of the pMK and MK models in the case of AA5182-O.

Figure 9 :

 9 Figure 9: Predictions of the pMK and MK models in the case of AISI 304.

Table 1 :

 1 Poly6 coefficients, eq. (76), for the AA5182-O and AISI 304 yield surfaces.

		a 1	a 2	a 3	a 4	a 5	a 6	a 7	a 8
		1.0000 -2.5109 8.0490 -12.360 8.8226 -2.996 1.0960 9.2518
	AA5182-O	a 9	a 10	a 11	a 12	a 13	a 14	a 15	a 16

Acknowledgement. The author thanks the National Council for Research and Innovation of Romania for supporting this research (Grant 4/30.04.2008 RP4). Also, thanks are addressed to the anonymous reviewers whose comments helped improve the original version of this work.

of the material, and hence the above predictions should be judged accordingly.

Conclusions

The present analysis of the MK model considers several features of this model that were not previously addressed on a sound theoretical basis. For the rate-independent constitutive model it was shown that the strains at failure in the groove are finite (and hence also finite in the rest of the sheet) and that eq. ( 14) holds true. The nature of the failure point (when the evolution B = B ( ) is considered) is in this case cuspidal. This is not true in the case of rate-dependent materials, if rate-dependence is introduced with the simple power law (32). In this case, the nature of the failure point is asymptotic: the failure strain in the bulk of the sheet is finite while the strain in the groove is infinite. This was first proved in Hutchinson and Neale (1978b). The equations of the MK model have been integrated explicitly in the plane strain case, when the groove is along one of the symmetry axes. This allows for a simple procedure for an apriori estimation of the MK parameter, the initial thickness ratio, if the plane strain limit strain is used as calibration data point. Furthermore, the analysis of the plane strain case reveals that the Hill-Swift failure criterion holds true in the MK-groove at the failure moment. This was generalized to any groove orientation: in the MK-groove Hill's criterion holds true once eq. ( 14) is satisfied. The equations of the MK model can also be integrated in the case of a non-hardening material.

The corresponding solution provides a useful insight into the way the yield surface influences the MK-predictions of the FLD. As a byproduct of these analyses, a simple approximate analytical solution for the MK-predictions, in the case of proportional loading along the symmetry axes and isotropic hardening, termed here HMK, was then developed. The plane strain solution obtained here has then been used to make a preliminary (and also elementary) phenomenological analysis of material heterogeneity as a cause for inhomogeneous plastic flow. Based on the conclusions drawn from this analysis, the MK model was then extended to sheets of uniform thickness by introducing a local eigenstrain-rate. Comparisons between the predictions of the extended model, the pMK model, and the classical MK model for two materials, AA5182-O and AISI304 steel, point to an equivalence between the two models, when the magnitude p of the eigenstrain-rate is constant (hence modeling the presence of a material heterogeneity). However, more refined evolution laws for the p parameter may account for different failure mechanisms.

Appendix A. Proof of Theorem 1.b in Section 4

We assume first that zone B is oriented along the transverse direction. From the yielding condition in zone B and stress equilibrium at its boundary we have:

where we t B := σ B y /σ B x . Using now the yielding condition in the bulk of the sheet we get further:

onto which the solution can be extended has the property: f y (1, t B ) = 0, that is, t B = t PS . With this solution we return now to eq. (A.4). Let us denote [0, E B F ) the maximal interval on which the numerator multiplying d B /d is strictly positive. Since d B /d > 0, it follows that we also have that the numerator of the left-hand member of eq. (A.4) is strictly positive. Note however, that due to the concavity of H the interval on which H ( ) -KH( ) > 0 is finite (see the proof of Theorem 1). Let us assume

Since H is concave, we have H ( B )/H( B ) -→ 0 as B -→ ∞ and then:

The case when the groove is oblique has similar features with the case of vertical groove. We only give here the new forms of the relevant equations. The stress state in zone B has the form in eq. ( 34). To condense the notation, we shall denote in what follows: r := q/σ x , T A + rT := 1 + rn 2 2 , t + rn 2 1 , -rn 1 n 2 , where T A := (1, t). Eq. (A.2) takes now the form: .11) whereas eq. (A.3) becomes:

where we denoted .13) with the constant K defined above. With eq. ( 10) the function R can be integrated explicitly but the only information about R needed here is that it is bounded:

The counterpart of eq. (A.4) is now: .15)