STRESS GRADIENT EFFECTS ON THE NUCLEATION AND PROPAGATION OF COHESIVE CRACKS

. The aim of the present work is to study the nucleation and propagation of cohesive cracks in two-dimensional elastic structures. The crack evolution is governed by Dugdale’s cohesive force model. Speciﬁcally, we in-vestigate the stabilizing eﬀect of the stress ﬁeld non-uniformity by introducing a length (cid:96) which characterizes the stress gradient in a neighborhood of the point where the crack nucleates. We distinguish two stages in the crack evolution: the ﬁrst one where the entire crack is submitted to cohesive forces, followed by a second one where a non-cohesive part appears. Assuming that the material characteristic length d c associated with Dugdale’s model is small in comparison with the dimension L of the body, we develop a two-scale approach and, using the methods of complex analysis, obtain the entire crack evolution with the loading in closed form. In particular, we show that the propagation is stable during the ﬁrst stage, but becomes unstable with a brutal crack length jump as soon as the non-cohesive crack part appears. We also discuss the inﬂuence of the problem parameters and study the sensitivity to imperfections.

Abstract. The aim of the present work is to study the nucleation and propagation of cohesive cracks in two-dimensional elastic structures. The crack evolution is governed by Dugdale's cohesive force model. Specifically, we investigate the stabilizing effect of the stress field non-uniformity by introducing a length which characterizes the stress gradient in a neighborhood of the point where the crack nucleates. We distinguish two stages in the crack evolution: the first one where the entire crack is submitted to cohesive forces, followed by a second one where a non-cohesive part appears. Assuming that the material characteristic length dc associated with Dugdale's model is small in comparison with the dimension L of the body, we develop a two-scale approach and, using the methods of complex analysis, obtain the entire crack evolution with the loading in closed form. In particular, we show that the propagation is stable during the first stage, but becomes unstable with a brutal crack length jump as soon as the non-cohesive crack part appears. We also discuss the influence of the problem parameters and study the sensitivity to imperfections.

Introduction.
1.1. Cohesive force models vs Griffith's model: State of the art. Griffith's theory of fracture [19] is based on the concept of critical energy release rate G c which comes from the fundamental but somewhat too restrictive assumption that the surface energy associated with a crack is proportional to the crack area (at least in a homogeneous and isotropic body), or, equivalently, that there is no interaction between the crack lips. It remains the most used approach in fracture mechanics thanks to its simplicity in terms of material behavior. However, this theory contains some major drawbacks. In particular, since Griffith's model does not contain a critical stress, (i) it allows stress singularity and (ii) cannot account for the crack nucleation in a sound body. Accordingly, cohesive-force models have been introduced in order to prohibit these types of unphysical singularities by allowing only finite stresses. Particularly, following the ideas of [14] and [4], many such models have been proposed and tested, see for instance [33,28,22,31,32,13].
In the first stage of their developments, these cohesive models were used in a restricted framework where the body contains preexisting cracks without consideration of their nucleation. In this context, various comparisons between Griffith's and Barenblatt's models were carried out to establish the precise contributions of the latter with respect to the former. In particular, under monotonic loading which does not require to introduce any irreversibility condition, it was rigorously proved, first by [26] in a restricted one-dimensional setting, then by [18] in three-dimensional, that Barenblatt's model leads to a crack propagation law which converges (in the sense of Gamma-convergence) to the Griffith's law when, for a given G c , the ratio between the material characteristic length (which is necessarily present in Barenblatt's model) and the body size goes to 0, see also [34] for a more formal proof. In essence this result means that the role of cohesive forces becomes negligible, as far as only crack propagation is concerned, once the crack length is sufficiently large. In such a case, the cohesive zone is essentially concentrated in a neighborhood of the crack tip, its size being of the order of the material characteristic length. In fact, the cohesive forces are then useful only to correct the shape of the crack opening near its tip by eliminating any stress singularity.
A first fundamental difference between the two models appears when one uses them in the case of cyclic loading. On the one hand, with Griffith's model, a crack can propagate only during the first cycle. On the other hand, after introducing a suitable irreversibility condition in Barenblatt's type model, like in the models proposed by [31,32,21], a crack can evolve from one cycle to the other and hence it becomes possible to account for the fatigue. Moreover, one can prove that the fatigue law induced by the cohesive force model converges to a Paris-like fatigue law [30] when the material characteristic length is small in comparison with the size of the body. This fundamental result was first observed in [29] by purely numerical considerations before being rigorously proved with Gamma-convergence argument by [21,5] in the restricted context of a peeling test. Its generalization to cracks in mode I or III is briefly presented in [1,2] and detailed in [3].
1.2. The issue of cracks nucleation with cohesive force models. Barenblatt's model radically differs from Griffith's model as far as the crack nucleation is concerned. Indeed, it turns out that the cohesive models, because they contain a critical stress, can also explain the process of crack nucleation in a sound body whereas Griffith's model cannot, in general. [12] was the first author who could establish such a result in a rigorous manner in an one-dimensional setting by using a variational approach and introducing a stability criterion. With such a criterion, a complete comparison between the two types of models can be carried out in the one-dimensional restricted setting, see [8]. Accordingly, assuming that the surface energy density which governs the cohesive forces is a smooth increasing concave function Φ(JuK) of the displacement jump JuK, the stability criterion requires that the stress field σ(x) at equilibrium is everywhere less than the derivative at 0 of the energy density function, i.e., σ(x) ≤ σ c := Φ (0). Therefore, σ c plays the role of the material critical stress. This result can be extended to a general three-dimensional setting via the same stability criterion. In particular, assuming that the material is isotropic and hence that the surface energy density is only a function of the normal displacement jump and of the norm of the tangential displacement jump across the crack lips, i.e., Φ(JuK · n, JuK − (JuK · n) n ), where n denotes the local unit normal vector to the crack, it is stated in [23] and proved in [9] that the crack nucleation criterion takes the form of an intrinsic curve in the Mohr stress plane which involves the directional derivatives at (0, 0) of Φ. Furthermore, when Φ admits partial derivatives at (0, 0), the nucleation criterion simply reduces to the two usual criteria of maximal shear stress and maximal tensile stress. This means that everywhere in the body the stress field must satisfy the following two inequalities where τ c and σ c denotes respectively the maximal shear stress and the maximal normal stress that the material can sustain. This result is fundamental in establishing a link between the cohesive cracks nucleation and the empirical materials strength criteria proposed by the engineers like Mohr and Caquot at the early of the twentieth century. However, this result only says that a cohesive crack will appear somewhere in the body when the stress field predicted by a pure elastic response reaches a threshold, but it says nothing about the growth process of these nucleated cracks. To treat this delicate issue, one must include in a single formulation both the nucleation and the (cohesive) cracks propagation. In essence, this is one of the main purposes of the variational approach to fracture, see [5] for an overview. In the context of cohesive force models, some partial results have already been obtained. For instance, [16,17] studied the size and shape effects of preexisting defects in the case of Dugdale's model. These authors show that the loading value at which the first cohesive crack occurs strongly depends on the preexisting defect shape. On the other hand, sufficiently small defects have practically no influence on the overall structure resistance.
Besides these first results on shape or size effects, [10,11] treated the nucleation and propagation of a cohesive crack problem at the notch tip in the same context of Dugdale's model. In such a situation, since the notch induced a stress singularity if the response is purely elastic, there exists no elastic phase in the loading process and a cohesive crack is created at the notch tip as soon as a load is applied. The crack length and its opening grow in such a manner that there is no stresses singularity. During a first stage of the loading, the crack growth is stable, but it becomes unstable when the crack opening just at the notch tip reaches the critical value δ c associated with Dugdale's model. Subsequently, a macroscopic crack is created by instability, which length is governed by an energy conservation condition. In the case where the material characteristic length is small in comparison with the overall body dimension, it is even possible, by using a two-scale technique, to obtain in closed form the formula giving the load at which a macroscopic appears at the notch tip.
1.3. The nucleation and growth of a crack at a regular point. The goal of the present work is to consider the same problem as in [10,11] but for another type of structure. In particular, we assume that the body contains neither a notch nor any corner which would induce elastic singularities. In other words the stress field associated with a pure elastic response is assumed to be smooth and bounded, but non-uniform. If we consider a symmetric structure submitted to an increasing loading and adopting Dugdale's law as the cohesive force model, a cohesive crack appears at a material point where the normal stress is maximal when the loading reaches a critical value. The question is then to study the crack propagation process and to highlight the stabilizing effects of the stress gradients. In particular, one shows that the crack growth is first progressive, which means that it depends continuously on the loading parameter, by virtue of the stress field non-uniformity. In fact, the first stage of the crack growth is controlled by the second derivatives of the stress field. In the second stage, when the loading reaches a value such that the crack opening at its center reaches the critical value δ c , then a non-cohesive zone appears in the center and the propagation becomes brutal, the crack size jumping instantaneously to a value which is fixed by the stress gradient characteristic length. This second critical loading value can be seen as the moment where a genuinely macroscopic and non-cohesive crack appears. The main feature of the paper is to obtain all the results in closed form by using the methods of complex potentials and a two-scale technique.
The paper is organized as follows. Section 2 is devoted to the problem setting and to the main assumptions, while section 3 contains its resolution and the major part of the results. In particular, we construct a solution in closed form by using a two-scale approach and the methods of complex potentials. This latter section finishes with a long discussion where one studies the influence of the parameters, specially the dependence on the material length d c and the stress gradient length , and the sensitivity of the response to the imperfections. A short comparison with Griffith's theory is also presented and, finally, the main ingredients of the resolution by using complex potentials are recalled in the appendix.
2. Problem setting and main assumptions.
2.1. The body, its elastic behavior and its loading. Throughout the paper the analysis is made in a plane strain setting. A Cartesian system (x 1 , x 2 , x 3 ) is used with its canonical orthonormal basis (e 1 , e 2 , e 3 ). The body reference configuration is the open subset Ω of R 2 in the plane (x 1 , x 2 ). The body is made of an isotropic brittle material whose elastic behavior before cracking is characterized by its Lamé coefficients λ and µ (or equivalently by its Young modulus E and its Poisson ratio ν). This cracks behavior in the material is governed by the Dugdale's model (see below for a precise statement of that model). The body is submitted to a proportional loading parameterized by the increasing parameter t > 0 called from now the time. Accordingly, if the response were purely elastic, then the displacement field u el (t) and the stress fields σ el (t) at time t would be the solutions of the following linear boundary value problem In (1), ε(u el (t)) denotes the strain field associated with the displacement field u el (t), i.e., the symmetric part of the gradient of u el (t). The body forces are neglected, ∂ N Ω represents the boundary part where the surface forces are (progressively) applied whereas ∂ D Ω represents the complementary boundary part where the displacements are prescribed. The loading is proportional in the sense that the applied forces intensity and the prescribed displacements amplitude are proportional to the parameter t.
By virtue of the linearity of the problem (1), its solution depends linearly on t and hence can be written tF O where (u el , σ el ) are solutions of the following linear elastic problem 2.2. Symmetry and smoothness assumptions. We assume that the body is symmetric with respect to the two axes x 1 = 0 and x 2 = 0. Moreover, the loading preserves this symmetry and the elastic response enjoys the following properties 1. The shear stress σ el 12 vanishes on the axes x 1 = 0 and x 2 = 0. Consequently, the stress tensor is diagonal in the basis (e 1 , e 2 ) at each point of the axes and its eigenvalues are respectively denoted σ el 1 and σ el 2 ; 2. The elastic stress field σ el (x) is a smooth function of x. The maximum of σ el nn (x) is reached at the origin O = (0, 0), in the direction n = e 2 and is positive.
These symmetry and smoothness assumptions on the elastic stress field induce some properties on the stresses distribution along the axis x 2 = 0 that will be useful in the sequel. In particular, 1. the stress vector σ el (x 1 , 0)e 2 is purely normal, say 2. the normal stress distribution Σ(x 1 ) is an even smooth function of x 1 which is maximum at x 1 = 0. Therefore, if we expand it with respect to x 1 , near x 1 = 0, up to the second order, we obtain Σ(x 1 ) = σ el 2 (0, 0) + where the normal stress at the origin σ el 2 (0, 0) is positive and its second derivative σ el 2,11 (0, 0) is negative.
the expansion of the normal stress distribution can read as This expansion up to the second order can be considered as a good approximation of Σ(x 1 ) provided that x 1 is small in comparison with .
Remark 1. If one considers that the expansion of the normal stress distribution up to the second order is the true stress distribution whatever x 1 , i.e., if Σ(x 1 ) is given by then the normal stress is maximal at 0, positive (in tension) for |x 1 | < / √ 2 and negative (in compression) for |x 1 | > / √ 2. That presence of a compression at a large distance from the origin will limit the crack propagation as we will see in Section 3.2.

2.3.
Dugdale's model of crack opening. The nucleation and the growth of cracks in the body are governed by Dugdale's cohesive force model whose main ingredients are recalled here. This model, formulated in energetic terms, is based on the fundamental assumption that the surface energy density Φ depends in a non-trivial manner on the displacement jump, unlike the Griffith's model in which Φ is assumed to be constant. So in Dugdale's model, by assuming that the crack is always in mode I, i.e., that only the normal displacement is discontinuous, the surface energy density reads : In (7), Ju n K denotes the normal displacement jump, G c is the critical energy release rate of Griffth's theory, whereas δ c is an internal length characteristic of the cohesive forces model. The ratio G c /δ c has the dimension of a stress, say σ c In terms of the cohesive forces, the normal stress σ nn giving the interaction between the crack lips is equal to σ c as long as 0 < Ju n K < δ c and vanishes as soon as Therefore, the crack lips are generally divided into two zones: the so-called cohesive zone in which the cohesive forces are equal to σ c and a so-called non-cohesive zone in which there are no cohesive forces. Indeed, if we follow [12], [8] or [5] and use a principle of energy minimization, it can be shown that the elastic response is no more a relative minimum of the total energy of the rod once the prescribed traction reaches the critical stress σ c . In a full three-dimensional context, [9] show that the elastic response is an energy local minimum only if the stresses are less than σ c everywhere in the body and thus that σ c enters in the crack initiation criterion. The direct consequence of the presence of a critical stress in the model is that a given structure can only sustain loads of limited amplitude.
Remark 3. The length δ c characterizes the crack critical opening from which no more cohesive forces exist in Dugdale's model. It is a material characteristic length. But, in the plane strain calculations, another material characteristic length appears which involves also the materials elastic properties. This length is defined by and gives the magnitude order of the cohesive zone length. In practice, since E is much greater than σ c for usual materials, d c is much greater than δ c .

2.4.
General formulation of the crack evolution problem. Owing to the symmetry and the smoothness assumptions above, we assume that a crack will nucleate at (0, 0) at the critical time t e when the maximal tensile stress associated with the elastic response reaches the critical value σ c , i.e., Then it is supposed that the crack will remain straight and along the axis x 2 = 0. Accordingly, the crack at time t, defined as the set of points where the displacement u(t) is discontinuous and denoted by S u(t) , is a subset of the part Γ of the axis x 2 = 0 included in Ω: Moreover, still by symmetry, the displacement jump across the crack lips will be assumed to be purely normal. So, the crack is in mode I and the normal displacement jump Ju 2 (t)K(x 1 ) at (x 1 , 0) ∈ S u(t) is called the crack opening at x 1 .
Under this assumption on the crack path, the problem giving the displacement field at time t and hence the crack state at that time can be formulated by using a variational approach, like in [16,5,17]. Specifically, let U ad (t) be the set of kinematically admissible displacement fields at time t, i.e., the set of smooth vector fields which satisfy the kinematic boundary conditions, which are allowed to jump only on Γ and whose normal jump is non-negative: where H 1 denotes the usual Sobolev space equipped with its natural norm · 1 . For a given t, at each u * ∈ U ad (t) is associated to the body total energy E t (u * ) as the sum of its elastic energy, its surface energy and the potential of the applied forces. Namely, E t (u * ) reads as We are now in a position to give a precise formulation of the crack evolution problem. That leads to the following Definition 1 (Variational Formulation of the crack evolution problem). At each time t ≥ 0, the displacement field u(t) is a local minimizer of the total energy E t among the set of all kinematically admissible displacement fields U ad (t). Accordingly, u(t) must be such that The stress field at time t is given by σ(t) = λtr ε(u(t))I + 2µε(u(t)), whereas the crack at time t corresponds to the jump set S u(t) .

Remark 4. (Various advantages of a variational formulation)
The variational formulation of the crack evolution problem presents several benefits in comparison with other formulations only based on equilibrium equations and constitutive conditions, namely 1. That allows to formulate the problem in a condensed form which remains valid even if one changes the behavior, the loading or the geometry; 2. It contains in a unique formulation both the equilibrium and stability concepts. Indeed, it turns out that the classical formulations based on equilibrium equations and constitutive conditions are simply first order stability conditions, as it is proved in Proposition 1; 3. The variational formulation supplies natural numerical methods to construct approximate solutions, see [25].
Remark 5 (Absence of irreversibility conditions). Let us emphasize that no irreversibility conditions have been introduced in Dugdale's model and hence in the evolution problem. The absence of an explicit irreversibility condition allows us to simplify the presentation, but can be seen as a weakness of the formulation, because that could lead to unphysical responses. We will discuss this point when we will construct a solution in the next sections. The reader interested by a complete formulation taking account of the irreversibility should refer to [21,5,3] where the fatigue modeling issue is considered.
Let us now establish the local conditions that u(t) and σ(t) must satisfy to be a solution of the variational problem stated in Definition 1.
Proposition 1 (First order stability conditions). The displacement field u(t) and its associated stress field σ(t) are solutions of the variational problem of Definition 1 only if they satisfy the following local conditions • Equilibrium equations: • Boundary conditions: • Crack path conditions: Proof. We only give a sketch of the proof which is based on classical variational arguments. However, the proof is valid for any cohesive force model and not merely for Dugdale's model. The method consists in considering kinematically admissible displacement fields of the form u * = u(t) + hv with h > 0 sufficiently small. Inserting such a u * into the stability condition gives E t (u(t)) ≤ E t (u(t) + hv). Then dividing by h and passing to the limit when h → 0 lead to Let us consider different types of fields v. 1. Let v be a smooth field such that JvK = 0 on Γ and v = 0 on ∂ D Ω. Then by classical arguments one deduces that 2. After an integration by parts of the first integral of (14) and using (15), we obtain the following inequality which must hold for any admissible v. 3. Since Jv 1 K can be chosen arbitrarily on Γ, one gets σ(t) 12 = 0 on Γ. Consequently, after inserting the latter equality in (16) and dividing Γ into S u(t) and Γ \ S u(t) , (16) becomes where σ c := Φ (0+). Finally, since Jv 2 K can be chosen arbitrarily on S u(t) whereas Jv 2 K is necessarily non-negative on Γ\S u(t) in order that Ju(t) 2 + hv 2 K ≥ 0 on Γ, one gets

TUAN HIEP PHAM, JÉRÔME LAVERNE AND JEAN-JACQUES MARIGO
The proof is complete.
Remark 6. Let us emphasize the most important results contained in Proposition 1.
• The equilibrium equations, boundary conditions and crack path conditions are only necessary conditions in order that u(t) be stable in the sense of Definition 1. In general, they are not sufficient and one must add second order stability conditions. However, in the present paper, we will not introduce these second order stability conditions and the interested reader should refer to [8,5] for more details on their use. • The crack path conditions contain not only the fact that the cohesive forces on the crack lips are given by the derivative of the surface energy density, namely σ(t) 22 = Φ (Ju(t) 2 K), but also the fact that σ c = Φ (0+) plays the role of a yield stress criterion for the crack nucleation. Indeed, one must have σ(t) 22 ≤ σ c everywhere (on Γ). Accordingly, the stress field is necessarily bounded and no singularity is allowed.
3. Resolution of the crack evolution problem.

3.1.
Reduction of the problem with the help of the symmetry assumptions. The evolution problem which is stated above in its general form can be reduced with the help of the symmetry assumptions introduced in Subsection 2.2.
Let us note however that, since the solution uniqueness is not ensured, the search for a solution respecting these symmetries constitutes additional assumptions. The first extra assumption is the following Hypothesis 1 (Centered crack path). At every time t > 0, the jump set S u(t) is either empty or an interval centered at (0, 0), i.e., there exists a(t) ≥ 0 such that S u(t) = (−a(t), a(t)) × {0}.
The second assumptions concerns the monotonicity of the crack opening.
Hypothesis 2 (Symmetry and monotonicity of Ju(t) 2 K). When the jump set S u(t) is not empty, the opening Ju(t) 2 K(x 1 ) is an even continuous function of x 1 , maximal at x 1 = 0 and decreasing to 0 when |x 1 | grows to a(t).
This second assumption limits the number of possibilities for the crack state at a given t. Specifically, we can distinguish the following three cases, see also  i. there exists no crack, i.e., a(t) = 0 and S u(t) = ∅. The response is purely elastic and the set of all crack states of this type which satisfy the first order stability conditions is called the elastic branch; ii. a crack exists but its opening at x 1 = 0 is less than the critical value δ c associated with Dugdale's model, i.e., a(t) > 0 and Ju(t) 2 K(0) ≤ δ c . That corresponds to the case where the entire crack lips are submitted to the cohesive force σ c . The set of all crack states of this type which satisfy the first order stability conditions is called the fully cohesive branch; iii. a crack exists and its opening at x 1 = 0 is greater than the critical value δ c , i.e., a(t) > 0 and Ju(t) 2 K(0) > δ c . In that case, since the opening is a monotonic function of |x 1 | decreasing to 0, there exists two symmetrical points (±b(t), 0) with 0 < b(t) < a(t) where the opening is equal to δ c . Therefore, by virtue of Dugdale's model, the crack is divided into two parts: The set of all crack states of this type which satisfy the first order stability conditions is called the partially non-cohesive branch. Of course the positions a(t) and b(t) of the tips of the cohesive zone and non-cohesive zone, when they exist, have to be determined. This is the absence of singularity which supplies the equation giving a(t), as it is shown in the following Proposition: Proposition 2 (Vanishing of the stress intensity factor K I at the tips of the cohesive crack). Since the normal stress σ(t) 22 must be bounded everywhere on Γ, no singularity can exist at the cohesive crack tips and hence the stress intensity factor K I must vanish at x 1 = ±a(t).
Proof. Let us consider the case where a crack exists, i.e., a(t) > 0. Correspondingly, the displacement field u(t) and the stress field σ(t) must satisfy in a neighborhood of the crack tips (±a(t), 0) the elasticity equations with the boundary conditions σ(t)e 2 = σ c e 2 on the lips of the crack (close to the tips). Therefore, we are in the situation of a crack in a linear elastic isotropic medium submitted to Neumann's boundary conditions. The structure of the solution is then well known, cf [6,20], and contains a priori a singular part. Specifically, by virtue of the symmetry assumptions, the crack is in mode I and the displacement u(t) in a neighborhood of the tip (a(t), 0) reads as where (r, θ) denotes the polar coordinates of x, i.e., x = r cos θe 1 + sin θe 2 . Accordingly, the normal displacement jump on the crack lips near (a(t), 0) is given by Therefore, one must have K I ≥ 0 in order that Ju(t) 2 K(x 1 ) ≥ 0. On the other hand, the normal stress field σ(t) 22 on the axis θ = 0 near the tip (a(t), 0) reads as + regular terms.
But since σ(t) 22 (x 1 ) must be not greater than σ c , one must have also K I ≤ 0 and hence finally K I = 0.
As far as the position b(t) of the non-cohesive crack tips is concerned, we simply have by construction Proposition 3 (Critical opening at the non-cohesive crack tip). The position b(t) of the non-cohesive crack tips, when they exist, must be such that the opening at these tips be equal to δ c : We are now in a position to exhibit a method for constructing a solution of the crack evolution problem. The procedure is as follows i. One solves the elastic problem and determine the fields (u el , σ el ). One deduces the elastic branch which corresponds to u(t) = tu el for 0 ≤ t ≤ t e . (Indeed, for t > t e , the elastic response tu el cannot be a solution of the crack evolution problem, because tσ el 22 > σ c somewhere on Γ by virtue of the definition of t e .) ii. One considers the case of a fully cohesive crack with a length 2a > 0 at time t > 0. For given a and t, we define the associated displacement and stress fields as the unique solution, denoted (u[t, a], σ[t, a]), of the following linear elastic problem posed on the cracked body with uniform cohesive forces on the crack lips: Note that this problem admits the same symmetries as the original elastic problem and hence its solution too. From its solution, one deduces the value of the mode I stress intensity factor K I [t, a] at the cohesive crack tips. Requiring that it vanishes, one obtains the right value a(t) of the crack tip position: Of course, it is not ensured at this stage that there exists a unique solution of (21) for a(t). If several solutions exist, then one can define several fully cohesive branches, but one can expect than only one starts from a = 0 at t = t e . One defines (u(t), σ(t)) by (u[t, a(t)], σ[t, a(t)]) and one must find for what values of t (u(t), σ(t)) satisfies all the first order stability conditions. The main remaining condition is that the opening must be positive and less than δ c everywhere on the crack lips: If this condition is satisfied only when t lies in some interval (t e , t i ), then this interval will constitute the fully cohesive branch. Finally, it will remain to check that σ(t) 22 ≤ σ c everywhere on Γ for those values of t. iii. One considers finally the case of a partially non-cohesive crack at time t whose non-cohesive length is 2b and the cohesive zones tips are at ±a. For given (a, b, t) with 0 < b < a and t > 0, we define the associated displacement and stress fields as the unique solution, denoted (u[t, a, b], σ[t, a, b]), of the following linear elastic problem posed on the cracked body with non-uniform cohesive forces on the crack lips: This problem also admits the same symmetries as the original elastic problem and hence its solution too. From the solution, one deduces the value of the mode I stress intensity factor K I [t, a, b] at the cohesive crack tips. Requiring that it vanishes, one obtains a first equation for the right values a(t) and b(t) of the crack tips position at time t. The second equation is given by the opening at ±b(t). Accordingly, the system of equations for (a(t), b(t)) reads as We are ensured neither that a solution of (24) for (a(t), b(t)) exists, nor that the solution is unique. We can expect that a solution exists only for some values of t. If several solutions exist, then one can define several partially non-cohesive branches, but one can expect than only one starts from a = a i and b = 0 at t = t i . Then, for any solution (a(t), b(t)) at time t, one defines the displacement and stress fields (u(t), σ(t)) by (u[t, a(t), b(t)], σ[t, a(t), b(t)]). It remains to check that (u(t), σ(t)) really satisfies all the first order stability conditions. Accordingly, one must verify that 3.2. Analytical calculation with a two-scale approach. In this section we construct a solution of the crack evolution problem in closed form, but under the condition that the material characteristic length d c defined in (10) is small in comparison with the characteristic length L of the body. Moreover, we assume that the stress field is genuinely non-uniform by considering that is of the same order as or much smaller than L : d c L, L.
This lengths hierarchy allows us to use a two-scale approach to construct the solution (which will be hence an only approximate solution). The construction follows the procedure described in the previous subsection.
3.2.1. Determination of the elastic branch. Once the original elastic problem (3) is solved, one obtains the normal stress distribution σ el 22 (x 1 ) along the axis Γ and therefore by (11) the time t e which corresponds to the validity limit of the elastic response.

3.2.2.
Determination of the fully cohesive branch. Let t > 0 and a such that 0 < a L. In such a case, since the crack length is small, the crack should perturb the elastic fields only in a neighborhood of the origin. Therefore, if we introduce in (20) the gaps of the solution with the elastic fields, i.e., a](x) should tend to 0 when x becomes large in comparison with a. Moreover, on the crack lips, the gap of the normal stress verifiesσ[t, a] 22 (x 1 ) = σ c − tΣ(x 1 ) where Σ(x 1 ) is given by (6). Accordingly, using (11) leads tō When a we can neglect the term o(x 2 1 ) and only consider the first two terms of the expansion. But even if a is of the same order as , one can consider that the parabolic normal stress distribution is the simplest case to study the influence of the stress gradient on the crack nucleation. These considerations allow us to write the problem giving the gaps in a neighborhood of the origin as follows Thus the construction of this simplified problem is based on an approximation and an assumption: 1. The problem is posed on the whole plane R 2 with the condition that the stresses must decrease to 0 at infinity. We use the fact that a is small in comparison with L and hence that we must recover the elastic solution far enough from the crack. The approximation simply consists in treating the problem at the scale of a and in sending to infinity the body boundary; 2. The normal stress distribution on the crack lips is assumed to be quadratic.
This can be considered either as a particular case of loading or an approximation valid when a is small in comparison with .
Owing to these approximations it becomes possible to solve the problem (26) in closed form. Indeed the simplified problem is a particular case of a family of plane elastic problems which can be solved with the methods of complex potentials developed by [27]. The main steps of the method are recalled in the appendix and we can directly use the results by identifying the normal stress distribution with The stress intensity factor K I [t, a] is given by (50) and after the calculation of the integral one gets The condition K I [t, a] = 0 gives the crack tips position in function of time: Hence the crack length is proportional to the characteristic stress gradient length, starts from 0 at t e and then increases with time. This solution is valid as long as the crack opening at x 1 = 0 remains less than δ c , and provided that the normal stress is less that σ c all along the axis. The normal stress and the opening are obtained by integration. First, from (48) one gets where a(t) is given by (27). Therefore, by virtue of (46), one obtains One deduces thatσ[t, a(t)] 22 (x 1 , 0) is a monotonically decreasing function of |x 1 |, decreasing from σ c (t/t e − 1) at the crack tips to 0 at infinity. It is then easy to check that the normal stress σ(t) 22 (x 1 , 0) is less than σ c for all x 1 .
Remark 7. The above expressions of the normal stress are based on the assumption that a(t) L and are only valid at a small scale. With these approximations we can simply conclude that the normal stress is less than σ c at small scale, i.e., in a neighborhood of the origin. It could happen that the maximal traction criterion be reached at another point (far from the origin) at time t. In such a case, another crack would nucleate at that point. But in the first stage of their growth, the cracks do not interact each other and the present procedure remains valid.
Let us now determine the crack opening. Let z ± = x 1 ± i0 be the points on the crack lips at x 1 , |x 1 | ≤ a(t). It comes from (28) that and hence that

TUAN HIEP PHAM, JÉRÔME LAVERNE AND JEAN-JACQUES MARIGO
By a straightforward integration and using the fact that JϕK(±a(t)) = 0, one gets JϕK(x 1 ). Finally, by virtue of (47), the opening reads as Let us remark that the opening is an even function of x 1 , maximal at the origin which justifies a posteriori Hypothesis 2. Using (27) and the definition (10) of d c , the opening at the origin can read as Hence the opening at x 1 = 0 is a monotonically increasing function of t for t ≥ t e . It will reach the critical value δ c at time t i given by The time t i corresponds to the end of the fully cohesive branch, after which a noncohesive zone will appear at the crack center. At that time, the cohesive crack half-length is a i = a(t i ). When d c is much smaller than , t i can be approximated by and the crack half-length at time t i is given by

Let us note that the magnitude order of a i is intermediate between d c and :
d c a i .

3.2.3.
Determination of the partially non-cohesive branch. Let us now consider, at a given time t > 0, the case of a partially non-cohesive crack whose non-cohesive length is 2b whereas the cohesive zones tips are at ±a. We assume that 0 < b < a L and hence, still, that the crack perturbs the elastic fields only in a neighborhood of the origin. Therefore, introducing in (23) the gaps of the solution with the elastic fields, i.e.,ū and using the same approximations as in the case of a fully cohesive crack, the local problem reads as follows where

COHESIVE CRACKS NUCLEATION IN NON-UNIFORM STRESS FIELDS 573
The problem (33) can be solved in closed form, because it is still a particular case of the family of plane elastic problems presented in the appendix. For a given triple (t, a, b), the solution u[t, a, b] is unique up to a rigid displacement field and the stress field σ[t, a, b] is unique. For a given t > 0, the solution belongs to the partially non-cohesive branch only if a and b are such that The stress intensity factor K I [t, a, b] is given by (50) and after the calculation of the integral one gets The vanishing condition of the stress intensity factor gives a first relation between a, b and t, namely The calculation of the opening Ju[t, a, b] 2 K(b) requires to determine the jump of the complex potential ϕ(z) on the cohesive crack lips, i.e., JϕK(x 1 ) for b < |x 1 | < a.

Using (34), (37) and (48), after a tedious calculation of the integral one eventually gets
Since JϕK(a) = 0, an integration of (38) leads to and finally, by using (51), one obtains the opening at the non-cohesive crack tip The requirement that this opening is equal to δ c gives the second relation between a, b and t. This relation reads as follows Remark 8. It is possible to prove that if a triple (t, a, b) satisfies (37) and (40), then the associated fields u[t, a, b] and σ[t, a, b] satisfy the conditions (25) and hence verify all the first order stability conditions. The proof is based on a careful study of the complex potential ϕ(z) for z = x 1 ± i0, but the calculations are too long to be reproduced here.

Representation of the three branches.
(i) The elastic branch starts at t = 0 and finishes at t = t e given by (11). All along this branch, there is no crack, a(t) = 0. Therefore, the elastic branch corresponds to the segment line [0, t e ] × {0} in the (t, a) diagram.
(ii) The fully cohesive branch starts at t = t e and finishes at t = t i given by (32). Along this branch, the length a(t) of the cohesive crack grows continuously with t from 0 to a i , a(t) being given by (27). Therefore, the fully cohesive branch corresponds to the monotonic curve represented on Figure 5 in the (t, a) diagram which starts from the point (t e , 0) and finishes at the point (t i , a i ).
(iii) For the partially non-cohesive branch, a, b and t are related by the two conditions (37) and (40). For studying these two conditions, let us set Subsequently, using (37), (40) can read as For a given α ∈ (0, 1), (41) is a cubic equation forā := a/ which depends on the parameter := d c / . It turns out that this equation admits a unique solution, a (α), whose dependence on α is non-monotonic. Indeed,ā (α) starts from a i / = 1 − t e /t i at α = 0, then is first decreasing up to a m / before being increasing and finally tends to 1 when α tends to 1, cf Figure 4. Accordingly, (37) gives t/t e as a function of α which depends also on , sayt (α) .
The functiont (α) starts from t i /t e at α = 0 and is first monotonically decreasing up to t l /t e , that minimum being reached at α l . Subsequently,t (α) grows to infinity when α grows to 1, cf Figure 4. Finally, the evolution of b with α is given by the functionb (α): b =b (α) := αā (α).
As it is shown on Figure 4, b is a monotonically increasing function of α, starting from 0 at α = 0 and tending to when α tends to 1. Accordingly, the triples (a, b, t) satisfying (37) and (40) can be seen as two parametric curves (t(α), a(α)) and (t(α), b(α)) parameterized by α ∈ (0, 1) and depending on and on the ratio d c / . In particular the curve (t(α), a(α)) represents the partially non-cohesive branch in the (t, a) diagram, cf Figure 5. Since the function a (α) andt (α) are non-monotonic and monotonically decreasing for small α, the partially non-cohesive branch contains a snap-back in a neighborhood of (t i , a i ) and a limit point (t l , a l ), both points depending on and d c . Accordingly, the branch has the shape of a loop which can be divided into two parts: the lower part between (t i , a i ) and (t l , a l ), the upper part after (t l , a l ).
• Finally the three branches can be represented in a diagram (t, a) and one obtains typically the curves plotted in Figure 5. Their dependence on and d c will be discussed in the next section.
Remark 9. The fact that a and b tends to a limit, namely , when t tends to infinity is due to the fact that the elastic response leads to a negative normal stress distribution at large distance of the origin, see Remark 1. Note however that the limit is greater than the distance / √ 2 at which a compression appears in the elastic response.   3.3. Discussion.

3.3.1.
Dependence of the curves on the characteristic lengths d c and . We assume here that the critical stress σ c is fixed and study the dependence of the Dugdale's branches on d c at fixed , or, on at fixed d c . Therefore, in any case, the loading t e at which a cohesive crack nucleates is fixed and hence the elastic branch is always the same.
• At fixed . For all d c , the fully cohesive branch is a part of the parabola a = 1 − t e /t. Only the final point (t i , a i ) depends on d c , and both t i and a i are increasing functions of d c (or ), see (31)-(32) and Figure 6. On the one hand, when d c (or ) goes to 0, then t i tends to t e and a i / tends to 0 like 1/3 . On the other hand, when d c / goes to infinity, then a i tends to and t i tends to infinity. That means that the smaller the material length d c , the weaker the stabilizing effect of the stress gradient.
In the same manner, for the partially non-cohesive branch, the smaller the material length d c , the more accentuated the snap-back and the larger the loop size. When d c tends to 0, the parameter α m of the lowest point of the loop tends to 1/e and hence the time t m tends to 2te π arccos(1/e) ≈ 0.760 t e whereas a m ≈ ed c tends to 0 like d c . The crack length a l of the limit point tends to /  • At fixed d c . For a given material, one can see the influence of the stress gradient intensity by comparing on Figure 7 the Dugdale's branches associated with different values of . Let us recall that the higher the stress gradient, the smaller the length , the case of a uniform stress field corresponding to = +∞. Accordingly, the higher the gradient, the greater the fully cohesive branch, the smaller the loop of the partially non-cohesive branch and the smaller the crack final length. For small stress gradient and hence large and small , the asymptotic behaviors when goes to 0 are the same as those presented above for fixed .

3.3.2.
The response under monotonically increasing loading. If the body is submitted to an increasing proportional loading starting from t = 0 and growing to infinity, then the response is purely elastic as long as t ≤ t e . In the interval (t e , t i ), a cohesive crack nucleates and its length continuously grows since the fully cohesive branch a(t) is monotonically increasing. At t = t i , the cohesive crack length is a i and its maximal opening (located at x 1 = 0) reaches the critical value δ c at which the cohesive forces disappear. As soon as the loading becomes greater than t i , necessarily a non-cohesive zone must appear at the crack center. But since the partially non-cohesive branch suffers a snap back, the response cannot follow its loop and the length of the crack must be discontinuous at t i . If one neglects inertial effects, then the unique possible stable configuration corresponds to a partially non-cohesive crack of total length a * i located on the upper part of the partially non-cohesive branch, cf Figure 8. Of course, the fact that one can neglect the inertial effects while the crack length suffers a jump discontinuity should be justified by a careful dynamical analysis. An alternative option would be to consider that the crack length jump is governed by an energy conservation principle. Such a study is outside the scope of the present work and the interested reader should refer to [15,7,24] for a complete analysis of dynamical crack propagation in the framework of Griffith's theory. Accordingly, if one adopts the assumption that the inertial effects are negligible, then the jump amplitude depends in particular on = d c / : the smaller , the greater the jump. Accordingly, for a given material length d c , the smaller the gradient stress length , the smaller the jump. On the one hand, for large values of (small stress gradient) and hence small values of , a i / is small and of the order of 1/3 whereas a * i is large and practically equal to . That means that the crack nucleation is brutal and the stabilizing effect of gradient stress is weak. On the other hand, for values of of the same order as d c , the jump is weak and the stabilizing effect of gradient stress is stronger, see Figure 7.
Remark 10. It would seem that the loop shape and the snap-back part of the partially non-cohesive branch do not play any role in the crack propagation under monotonic loading. Moreover one could believe that the snap-back part is a purely mathematical byproduct of our modeling where the irreversibility of the crack propagation is not taken into account. Indeed, even if one decreases the loading just after the end of the fully non-cohesive branch, it is not physically admissible to follow the partially non-cohesive branch because the total crack length a should decrease. In fact, we will show in the next paragraph that the loop can be observed and even that it plays a fundamental role in presence of imperfections.

Sensitivity to the imperfections.
Up to now, all the analysis is made in the ideal case where the body is homogeneous and does not contain any defect before the loading process. Such a situation will be called the perfect case by opposition to the preexisting defects case. In the present paper we will only consider the case where the imperfection corresponds to an initial cut along the x 2 = 0 axis, centered at 0 and of half-length a 0 < . In other words, we assume that the body contains a preexisting non-cohesive crack (−a 0 , a 0 ) × {0} whose length is a parameter. Accordingly, the elastic response is no more regular, but the stress is singular at the tips ±a 0 as soon as a loading is applied. Therefore, by virtue of Dugdale's model and Proposition 2, there exists no more an elastic branch, but a cohesive zone must nucleate ahead the tips ±a 0 as soon as t > 0 with a length a − a 0 such that the singularity vanishes at the tips ±a. Assuming that the initial crack is small in comparison with the size of the body, i.e., a 0 L, one can follow the same two-scale approach presented in Section 3.2. In particular we can use the expressions (36), at time t, for the stress intensity factor at the crack tips ±a whose non-cohesive zones are of length b, namely K I [t, a, b]. Therefore, the relation between a, b and t in order that the singularity vanishes remains given by (37). Similarly, the opening at the tips ±b, namely Ju[t, a, b] 2 K(b), is still given by (39). Equipped with those two relations, it is easy to determine the preexisting crack evolution under a monotonically increasing loading. Specifically, the evolution can be divided into the following two or three parts, according to the value of a 0 : 1. Cohesive phase: Growth of two symmetric purely cohesive zones, the noncohesive crack tips remaining at ±a 0 . For t small enough, the initial noncohesive crack does not propagate because the opening at ±a 0 remains less than δ c , but two symmetric cohesive zones grow in order to cancel the singularity at those points. The relation between a and t is given by the condition Since (44) gives a monotonically increasing relation between t and a when a ∈ [a 0 , ), the relation is invertible and hence a is an increasing function of t starting from a 0 at t = 0. That allows us to define the so-called cohesive branch associated with the initial crack length a 0 in the diagram (t, a). Moreover, for a and t satisfying (44), (39) with b = a 0 gives and hence the opening at ±a 0 is an increasing function of t starting from 0 at t = 0. By construction, it will reach the critical value δ c when the triple (a, a 0 , t) satisfies both (37) and (40). Therefore that triple is the point of the partially non-cohesive branch of the perfect case which corresponds to b = a 0 . The associated parameter α 0 is given by the equation its uniqueness being ensured by the monotonicity of the functionb (α). In other words the cohesive branch will finish when it intersects the loop of the perfect case. In conclusion, the cohesive branch starts from (0, a 0 ) and finishes at (t (α 0 )t e ,ā (α 0 ) ). During this phase, the total crack length and the opening of any point of the crack lips increase with t and hence there exists no incompatibility with an irreversibility condition. 2. Possible jump of the crack length: Brutal crack propagation if the cohesive branch intersects the lower part of the loop of the perfect case. If the final point of the cohesive branch is lower than the limit point of the loop, i.e., ifā (α 0 ) < a l , then the crack evolution will suffer a jump. Indeed, the evolution cannot follow the lower part of the loop in the increasing time direction for obvious irreversibility reasons since the crack length should decrease. Therefore, if one neglects inertial effects, the unique possibility is that the evolution be discontinuous and the point just after the jump be the point located at the same timet (α 0 )t e on the upper part of the loop. On the other hand, if the final point of the cohesive branch is at or above the limit point of the loop, i.e., ifā (α 0 ) ≥ a l , then the evolution can continuously follow that part of the curve in the increasing time direction since the crack length increases and no jump is necessary. 3. The continuous growth of a partially non-cohesive crack. Once the upper part of the loop has been reached, which can be arrived after a jump, the crack evolution simply follow that upper part of the loop in the increasing time direction and finally the crack length will tends to when t goes to infinity as in the perfect case. The system will finally forget its initial imperfection.
All these results can be seen on Figure 9 where are considered five cases of imperfection size. The first three, which correspond to a small initial crack length, lead to a jump whereas the last two, corresponding to a sufficiently large initial crack length, give rise to a continuous crack growth. Of course, the critical initial crack length above which the evolution is continuous depends both on and d c . In any case, one sees the fundamental role played by the loop of the perfect system.

3.3.4.
Comparison with Griffith's theory. To finish this discussion let us compare the evolution predicted by Dugdale's model with the one associated with Griffith's theory. To this purpose, let us consider an initial centered (non-cohesive) crack of half-length a 0 L and let us determine for which loading t 0 that crack will propagate if one uses Griffith's criterion. Since the initial crack is small, one can use the results of the two-scale approach. Since there is no cohesive forces in Griffith's theory, the stress field is singular at the crack tips ±a 0 and the stress intensity factor is given by (36) with b = a 0 . Accordingly, one gets Along the upper branch, above (t G , a G ), t 0 increases to infinity when a 0 increases to , whereas along the lower branch, below (t G , a G ), t 0 increases to infinity when a 0 decreases to 0. Accordingly, when a 0 < a G , the smaller the initial crack, the greater the loading at which it propagates. At the limit no crack can nucleate in a sound body, what is one of the main drawbacks of Griffith's theory. When the initial crack is such that 0 < a 0 < a G and if one neglects inertial effects, the crack length will jump instantaneously at t 0 to the associated point a * 0 on the upper Griffith's branch. Subsequently, the evolution will propagate continuously by following the upper Griffith's branch, see Figure 10.
Let us compare with Dugdale's law. For small values of , the upper part of the partially non-cohesive branch of Dugdale's model is close to the upper part of Griffith's branch. In particular, when tends to 0, a l tends to a G , the ratio t l /t G tends to 1 whereas both t l and t G go to 0 like √ . But the lower part of the partially non-cohesive branch and the fully cohesive branch of Dugdale's model remains different from the lower part of Griffith's branch. In particular, the loading at which a crack nucleates or a preexisting crack propagates cannot be greater than t e with Dugdale's model whereas it is not bounded but strongly dependent on the size of the preexisting crack with Griffith's model. That means that the nucleation and the first phase of the crack propagation are strongly different according to one uses Griffith's or Dugdale's model. But once the crack length is large in comparison with the Dugdale's characteristic length d c , then the cohesive zones become negligible and the two models give practically the same results.
4. Conclusion and perspectives. Let us summarize the main results obtained in this paper. First, since Dugdale's law contains a critical stress σ c , one can account for the crack nucleation in a sound body at a finite loading t e , in contrast with Griffith's law. However, only the first phase of the nucleation, at which the entire crack is submitted to cohesive forces, leads to a continuous crack length evolution with the loading. Indeed, at the loading t i when the opening reaches the critical value δ c , the evolution is necessarily discontinuous and leads to a crack length jump because of the presence of a snap-back in the equilibrium branch. Moreover, since Dugdale's model contains also a material characteristic length d c , size effects are possible. Assuming that d c is small in comparison with the size of the body, situation the most frequent in practice, the entire solutions can be obtained in closed form which renders easy the study of the size effects. In particular, one shows that, in presence of stress gradient, the response is very sensitive to the ratio between the material length d c and the stress gradient characteristic length . The smaller the stress gradient, the higher the length , the shorter the first phase of nucleation and the greater the crack length jump at t i which is of the order of . Accordingly, the loading t i at which the jump occurs can be considered as the loading at which a "macroscopic" crack nucleates in the body. The formula (32) which gives t i can be considered as universal, in the sense of not dependent on the particular problem (but specific to Dugdale's model), provided that is small. Finally, the snap-back in the partially non-cohesive branch, which is also obtained in closed form and whose first part can be also considered as universal, plays an important role in presence of imperfections. All these results which have been obtained for Dugdale's model should be extended for more general cohesive models. The two-scale approach can be followed in any case, but the difficulty will be to solve the different problems in closed form. The help of numerical methods could be necessary.
Appendix A. The generic local problem and its solving. Let us consider the following plane elastic problem which is set on the entire plane except a crack of length a centered at the origin in the x 1 direction: where T(x 1 ) represents the normal force distribution on the crack lips. The solution which is defined up to an arbitrary rigid displacement can be found by using complex potentials, cf [27]. Accordingly, the components of the displacement and of the stress are given in terms of the complex potential ϕ(z), z = x 1 + ix 2 . In particular, one has ϕ being holomorphic in the plane without the crack, the bar denoting the complex conjugate. By a standard procedure, we get from which one deduces the normal stress distribution along the axis x 2 = 0. Specifically, outside the crack lips, one gets Therefore the normal stress is in general singular at the crack tips with a singularity of the type σ 22 (a+r, 0) ∼ K I / √ 2πr for r close to 0. Accordingly, the relation between the mode I stress intensity factor K I at ±a (which is the same at the two tips by symmetry) and the normal stress distribution T is given by After an integration of (48), one obtains the normal displacement jump across the crack (the arbitrary rigid displacement does not play any role):