Christian J Cyron

Jean-Jacques Marigo

Paul Sicsic
email: paul.sicsic@polytechnique.edu

Cohesive model approach to the nucleation and propagation of cracks due to a thermal shock

Keywords: Thermal Shock, Cohesive Surface Energy, Crack Initiation, Micro-cracks, Stability

This paper studies the initiation of cohesive cracks in the thermal shock problem through a variational analysis. A two-dimensional semi-infinite slab with an imposed temperature drop on its free surface is considered. Assuming that cracks are periodically distributed and orthogonal to the surface, at short times we show that the optimum is a distribution of infinitely close cohesive cracks. This leads us to introduce a homogenized effective behavior which reveals to be stable for small times, thanks to the irreversibility. At a given loading cracks with a non-cohesive part nucleate. We characterize the periodic array of these macro-cracks between which the micro-cracks remain. Finally, for longer times, the cohesive behavior converges towards that from Griffith's evolution law. Numerical investigations complete and quantify the analytical results.

Introduction

Specific topology of crack patterns arise in a various of situation and are often characteristic of the material, the loading and the structure. The nucleation of these specific topologies, especially away from singularities remains a difficult task for fracture mechan-ics. In particular, Griffith's theory of fracture [START_REF] Griffith | The phenomena of rupture and flow in solids[END_REF], until now the most used, is unable to account for crack initiation in a sound elastic body. The predictive understanding of crack nucleation is still an elusive goal. Short of introducing initial flaws at the structural scale [START_REF] Bahr | Thermal-shock crack patterns explained by single and multiple crack propagation[END_REF], prescribing ad-hoc stress criteria [START_REF] Jagla | Stable propagation of an ordered array of cracks during directional drying[END_REF], or accepting global energy minimization arguments whose physical relevance is debated [START_REF] Marigo | Initiation of cracks in griffith's theory: An argument of continuity in favor of global minimization[END_REF], initiation can not be addressed.

A typical situation of crack nucleation is that of an array of parallel cracks after shrinkage. In this paper, we focus on the thermal shock problem of a brittle slab, for which experimental results are reported in [START_REF] Geyer | Experimental investigation of thermally induced interacting cracks in brittle solids[END_REF]; [START_REF] Jiang | A study of the mechanism of formation and numerical simulations of crack patterns in ceramics subjected to thermal shock[END_REF]. The specimen is a thin slab, free at the boundary, composed of a homogeneous material without prestress in its initial configuration. It is uniformly heated and then quenched in a cold bath inducing a thermal shock on the exposed surfaces. A first mechanism of selection of materials towards thermal shock based on their surface energy was introduced by [START_REF] Hasselman | Elastic energy at fracture and surface energy as design criteria for thermal shock[END_REF]. The crack selective arrest has been explained by [START_REF] Bazant | Stability in postcritical growth of a cooling or shrinking cracks[END_REF] using a bifurcation analysis based on the change of sign of the second derivative of the strain energy with respect to the crack penetration. More recently, the spacing and initiation by global minimization of the Griffith energy were derived by [START_REF] Jenkins | Optimal spacing and penetration of cracks in a shrinking slab[END_REF]. [START_REF] Bahr | Scaling behavior of thermal shock crack patterns and tunneling cracks driven by cooling or drying[END_REF] derived a scale law relating the space between non cohesive cracks and their penetration.

The prediction of nucleation of cracks with only Griffith's evolution law leads necessarily to the initiation of cracks of finite length. In this sense it is not satisfactory. The introduction of a critical stress either with in the cohesive setting or with gradient damage models [START_REF] Sicsic | Initiation of a periodic array of cracks in the thermal shock problem: a gradient damage modeling[END_REF] is required. In a setting close to the thermal shock, i.e. the cooling of a concrete af-ter hydration, [START_REF] Bazant | Cohesive fracturing and stresses caused by hydration heat in massive concrete wall[END_REF] use cohesive surface energies to link the crack opening to the crack spacing.

The variational approach to fracture, in the spirit of [START_REF] Francfort | Revisiting brittle fracture as an energy minimization problem[END_REF] is based on the definition of an energy and a stability principle (see [START_REF] Bourdin | The variationnal approach to fracture[END_REF], for a broad description of the developments). Cohesive surface energies introduced by Dugdale (1960); [START_REF] Barenblatt | The mathematical theory of equilibrium cracks in brittle fracture[END_REF] have been reformulate into a full evolution law by [START_REF] Charlotte | Revisiting brittle fracture as an energy minimization problem: comparison of griffith and barenblatt surface energy models[END_REF]; [START_REF] Marigo | Initiation and propagation of fracture in the models of griffith and barenblatt[END_REF]; [START_REF] Charlotte | Initiation of cracks with cohesive force models: a variational approach[END_REF]; [START_REF] Ferdjani | Insensitivity to small defects of the rupture of materials governed by the dugdale model[END_REF].

The aim of this paper is to give further insight into the initiation phenomenon in thermal shock fracture. By doing this one provides a non-trivial example of the study of the evolution and bifurcation problem of cohesive models; which reveals a micro-layer behavior. The terminology of microcracking can cover many situations. Here, we do not refer to a process zone nor the weakening effect of micro-cracks ahead of a main crack [START_REF] Ortiz | Microcrack coalescence and macroscopic crack growth initiation in brittle solids[END_REF]. Rather we will account for micro-cracks as a constitutive behavior of the material. In this spirit, [START_REF] Pichler | Instability during cohesive zone growth[END_REF] study the stability of a set of cracks in a representative element volume. We focus on the thermal shock problem for a two-dimensional slab, in a quasi-static setting. By assuming a perfect conductivity at the surface of the thermal shock, we consider a Dirichlet boundary condition on the temperature and use the analytically calculated temperature field, function of space and time, to evaluate the mechanical loading in the form of thermally induced inelastic strains.

In view of the experimental results [START_REF] Jiang | A study of the mechanism of formation and numerical simulations of crack patterns in ceramics subjected to thermal shock[END_REF] only cracks which are orthogonal to the exposed surface are considered. A cohesive surface energy density κ of Barenblatt-type is introduced. The variational formulation, in the wake of [START_REF] Bourdin | The variationnal approach to fracture[END_REF], is based on this surface energy and a directional stability principle. If the thermal loading is not sufficiently severe the elastic response is stable at any time. Otherwise cracks nucleate as soon as the loading begins. In the sequel the focus is set on the latter case, assuming that all cracks have the same length and are periodically distributed. In order to study the singularity due to the thermal shock and to simplify the presentation a dimensionless setting is introduced. The problem is thus governed by a single parameter η, the ratio between the materials critical stress and the maximal stress induced by the thermal loading.

This total dimensionless energy is first minimized at short times, which gives rise to a "limit" which can be seen as a distribution of infinitely close cohesive cracks. This type of microstructuration is well known in phase transform problems where the energy is not convex [START_REF] Ball | Fine phase mixtures as minimizers of energy[END_REF][START_REF] Bhattacharya | Relaxation of some multi-well problems[END_REF]. Thus, a homogenized (effective) behavior ψ r associated with this optimal microcracking structuration is established. This solution is then confirmed by computer simulations. The irreversibility in the cohesive case, which can be defined in different ways, is a real issue for the definition of the effective behavior. The effective behavior resembles that of a perfectly plastic material (Fig. 7) but with no unloading part due to the irreversibility constraint. Then the stability of this micro-cracked layer is studied by testing the nucleation of macro-cracks with a non cohesive part. Due to the irreversivility the lost of the stability always arises at a finite time t c (Property 2), at the exception of Dugdale's cohesive model for which the micro-layer is always stable.

The lost of stability arises for a single crack in the center of the domain. Keeping the hypothesis of periodicity spacing we search for a periodic array of macrocracks whose optimal periodicity d opt is characterized numerically. This numerical investigation is carried on to reveal the crack selection mechanism where every other crack stops as in the Griffith case.

Specifically the paper is organized as follows. In the following section the problem statement is given, as well as a precise description of the cohesive crack model (surface energy and evolution principle). The elastic solution is given and its stability studied. Then the minimization of a parallel array of cracks is introduced (Section 3). In Section 4, the initiation of cracks for short times is studied. This allows us to construct the homogenized behavior of this microcrack state. In section 4.3 the stability of these homogenized micro-layer is established. Section 5 tackles the issue of the evolution of the non-cohesive cracks.

From the technical standpoint, we essentially use basic tools of the Calculus of Variations [START_REF] Dacorogna | Direct methods in the calculus of variations[END_REF]. Table 1 summarizes the main nomenclature used in this article. Vectors and second order tensors are denoted by boldface letters, e.g., u for the displacement and σ for the stress. Their components are denoted by the respective non-boldface letters by lower indices like u 1 or σ 11 . The magnitude of a vector is denoted by the respective non-boldface letter without index, e.g., u. The inner product between two vectors or two tensors of the same order is indicated by a dot. The identity tensor is denoted I. Infimum of the Rayleigh ratio minimization (37) dc = Eδc/(ωσc) Characteristic length of the material λ = t/L Slenderness of the micro-crack layer We consider the plane domain Ω = (-L, +L) × R + with homogeneous initial temperature θ(x) = θ 0 at time t < 0. At t = 0 the thermal shock θ 0 -ϑ with ϑ > 0 is imposed on its boundary x 2 = 0 (Fig. 1) whereas the lateral sides x 1 = ±L are thermally insulated, i.e. ∂θ/∂x 1 = 0 on x 1 = ±L. From the mechanical viewpoint, the upper side x 2 = 0 is force free whereas the lateral sides are shear free and their normal displacement is blocked. Thus, for all t > 0, the boundary conditions read as

θ = θ 0 -ϑ, σ 12 = σ 22 = 0 on x 2 = 0, ∂θ ∂x 1 = 0, u 1 = 0, σ 21 = 0 on x 1 = ±L.
The domain is composed of an isotropic brittle material with Young's modulus E, Poisson ratio ν = 0, thermal expansion coefficient a and thermal diffusivity k. The change of temperature is taken into account through the isotropic thermal strain field

ε th t (x) = a(θ t (x) -θ 0 )I.
To calculate ε th t (x) one first needs to know the solution for the temperature field in space and time. Assuming that the heat diffusion is not influenced by the possible presence of cracks, the solution of the heat equation gives the temperature field

θ t (x) = θ 0 -ϑf c x 2 2 √ kt
where f c denotes the complementary error function

f c (y 2 ) = 1 -erf(y 2 ) = 1 - 2 √ π y2 0 e -ξ 2 dξ.
Note that f c is monotonically decreasing with f c (0) = 1, f c (0) = -2/ √ π and f c (+∞) = 0. Recalling that the Poisson ratio is taken equal to 0, the thermoelastic energy density is given by

ψ(ε e ) := 1 2 E(ε e 11 ) 2 + E(ε e 12 ) 2 + 1 2 E(ε e 22 ) 2 (1)
where the elastic strain ε e is related to the total strain ε and the thermal strain ε th by

ε e = ε -ε th .
The total strain field ε is related to the displacement field by ε = ε(u) := 1 2 (∇u + ∇ T u), i.e. , ε is the symmetrized linearized gradient of u, while the thermal strain is given by

ε th t (x) = -aϑf c x 2 2 √ kt I. (2)
Accordingly, the stress-strain relation reads as σ = Eε e .

Fig. 1 Geometry and boundary conditions at t > 0.

The elastic response

In absence of any crack, the elastic response of the body is such that the displacement field at time t takes the form

u el t (x) = U x 2 2 √ kt e 2
and its unique non null component U is such that σ 22 = 0 everywhere and at every positive time. Therefore

U (y 2 ) = +∞ y2 f c (ξ)dξ
and the stress fied is given by

σ el t (x) = Eaϑf c x 2 2 √ kt e 1 ⊗ e 1 .
So, the material is submitted to a tensile stress σ 11 which is maximal at the upper side x 2 where it is equal to Eaϑ and is decreasing to 0 when x 2 goes to infinity (at given time). At a given point x 2 > 0, the tensile stress is monotonically increasing with time.

The cohesive crack surface energy

Because of the tensile stresses induced by the heat diffusion process and the boundary conditions, cracks can nucleate as soon as the thermal shock is applied and propagate throughout the process. We will only consider cracks which are parallel to the x 2 direction. This a reasonable assumption by virtue of the geometry and the loading (and that corresponds to what is observed in experimental tests [START_REF] Jiang | A study of the mechanism of formation and numerical simulations of crack patterns in ceramics subjected to thermal shock[END_REF]). From the variational viewpoint this condition will be considered as a constraint which is prescribed to the crack path. We assume that the lips of those cracks are submitted to normal cohesive forces which are given in terms of the crack opening u 1 through a surface energy function κ. (We will neglect tangential cohesive forces associated with possible tangential displacement jumps.) Specifically, the cohesive surface energy density κ is of Barenblatt-type Barenblatt (1962), u 1 → κ( u 1 ) is increasing, concave and tends to a finite limit G c when u 1 goes to infinity. Specifically, introducing dimensionless quantities, κ can read as

κ( u 1 ) = σ c δ c κ ( u 1 /δ c ) with κ(0) = 0, κ (0) = 1, κ (0) ≤ 0, κ ≥ 0, κ ≤ 0, κ(∞) < ∞.
Accordingly, the normal cohesive force σ 11 is given by the derivative of κ with respect to u 1 :

σ 11 = κ ( u 1 ) = σ c κ ( u 1 /δ c ) (3)
and hence σ c represents the cohesive force corresponding to an infinitesimal crack opening whereas δ c is a material characteristic length. As a particular case, we will sometimes consider the following family of cohesive models which depends on the parameter ω ∈ [0, 1]:

κ(ζ) = ζ -ω 2 ζ 2 if 0 ≤ ζ ≤ 1 1 -ω 2 if ζ ≥ 1 .
The associated relationship between the cohesive force and the crack opening is given by

σ nn = (1 -ωζ)σ c if 0 ≤ ζ < 1 0 if ζ > 1 .
The case ω = 0 corresponds to Dugdale's model while the case ω = 1 corresponds to a linear cohesive force model. Note that G c := κ(+∞) = (1 -ω 2 )σ c δ c . This family of models is illustrated in Fig. 2 where δ c is adjusted so that the values of σ c and G c be the same for all models. Remark 1 The cohesive law (3) is valid as long as the opening is an increasing function of time. When an unloading occurs, it is not reasonable from the physical viewpoint to assume that the behavior is reversible and that the relation between the cohesive force and the opening is still given by (3), one must consider an irreversibility condition. However, to simplify the presentation, that irreversibility condition will be only introduced in Section 4. Accordingly, we will work throughout this section without any irreversible condition and call the corresponding setting the reversible case.

Variational formulation of the problem in the reversible case

Let v be a displacement field which satisfies the boundary condition v 1 = 0 on x 1 = ±L and which is possibly discontinuous across a family of lines parallel to the x 2 direction. This set of points where v is discontinuous is called the jump set of v and denoted S v . In order that there is no interpenetration, the jump discontinu-ity must satisfy v 1 ≥ 0 on S v . Such a v is said kinematically admissible and one associates to it the total mechanical energy of the body at time t by

E t (v) = Ω\Sv ψ(ε(v) -ε th t ) dx + Sv κ( v 1 ) dx 2 , ( 4
)
where ε th t denotes the thermal strain field at time t given by (2). Following the variational approach proposed in [START_REF] Bourdin | The variationnal approach to fracture[END_REF], the real displacement field at time t, say u t , must be found among the kinematically displacement fields which are local minima of the energy at time t. Specifically, for any kinematically admissible displacement field v, it must exist h > 0 such that for all h ∈ (0, h) the following inequality holds

E t (u t + h(v -u t )) ≥ E t (u t ).
(5) (This definition of stability is in fact a condition for u t to be a local minimum in each given direction v, h depending in general on the direction v.) Expanding this inequality with respect to h up to the second order leads to

0 ≤ hE t (u t )(v -u t ) + 1 2 h 2 E t (u t )(v -u t ) + o(h 2 ), ( 6
)
where E t (u)(v) and E t (u)(v) denote respectively the first and second directional derivatives of E t at u in the direction v, i.e.

E t (u)(v) := Ω\(Su∪Sv) E(ε(u) -ε th t ) • ε(v) dx + Su∪Sv κ ( u 1 ) v 1 dx 2 E t (u)(v) := Ω\Sv Eε(v) • ε(v) dx + Su∪Sv κ ( u 1 ) v 1 2 dx 2 .
Dividing by h and passing to the limit when h goes to 0 in (6), one obtains the so-called first order stability condition,

E t (u t )(v -u t ) ≥ 0, ∀v kinematically admissible. ( 7
)
This condition is only a necessary stability condition which is not sufficient in general. Specifically, let u t be an admissible field which satisfies (7). In any direction v such that E t (u t )(v -u t ) > 0, (5) holds for h small enough. But, in the directions v such that E t (u t )(vu t ) = 0, one must study the sign of the second order term. That leads to the second order stability condition which reads as

E t (u t )(v -u t ) ≥ 0, ∀v kinematically admissible (8) such that E t (u t )(v -u t ) = 0.
The condition ( 8) is still a necessary condition of stability and, coupled with (7), it becomes a sufficient condition when the strict inequality holds true for all non trivial v. By standard arguments of Calculus of Variations [START_REF] Bourdin | The variationnal approach to fracture[END_REF], it turns out that the first order stability condition is satisfied at time t if and only if the displacement field u t and the associated stress field σ t are such that:

Stress-strain relation

σ t = E(ε(u t ) -ε th t ) in Ω \ S ut (9) Equilibrium div σ t = 0 in Ω \ S ut (10)
Limit stress condition

σ t11 ≤ σ c in Ω \ S ut (11)
Cohesive force relation

σ t e 1 = κ ( u t1 )e 1 on S ut (12)
Boundary conditions σ t e 2 = 0 on x 2 = 0,

σ t21 = 0, u t1 = 0 on x 1 = ±L (13)
Let us note in particular that the tensile normal stress σ 11 must be less than the critical value σ c everywhere and even in the uncracked part of the body. (In fact, this condition would be enlarged and would require that the maximal principal stress max i σ i must be less than σ c everywhere if we had considered arbitrary crack directions, see [START_REF] Charlotte | Initiation of cracks with cohesive force models: a variational approach[END_REF].) Note also that, since the cohesive surface energy density does not depend on the tangential displacement jump, the variational formulation automatically gives that there is no shear stress on the lips of the cracks.

2.5 Discussion on the existence and the stability of the elastic response By construction the jump set of the elastic displacement field u el t is empty and the associated stress field σ el t satisfies the equilibrium equation and the boundary conditions. Since σ el t11 is maximal at x 2 = 0 where it is equal to Eaϑ, the limit stress condition is satisfied if and only if the thermal shock is small enough, specifically iff ϑ ≤ σ c /Ea. Let us distinguish the two cases.

1. Case ϑ ≤ σ c /Ea. The elastic response satisfies the first order stability condition (at any time). It remains to check that it satisfies the complete stability condition (5). A direct calculation using the particular form of the elastic response gives, for any kinematically admissible field v and at any time:

E t (u el t )(v-u el t ) = Ω\Sv σ el t11 ∂v 1 ∂x 1 dx+ Sv σ c v 1 dx 2 .
Integrating by parts the first term in the right-hand side above and using the fact that the elastic stress field does not depend on x 1 lead to

E t (u el t )(v -u el t ) = Sv (σ c -σ el t11 ) v 1 dx 2 ≥ 0
where the inequality holds because σ el t11 < σ c and

v 1 ≥ 0. Moreover E t (u el t )(v -u el t ) = 0 iff the di- rection v is such that v 1 = 0 on S v .
For such a direction, the cohesive surface energy vanishes and one gets

E t (u el t +hv)-E t (u el t ) = h 2 2 Ω\Sv Eε(v)•ε(v) dx ≥ 0.
Therefore the elastic response is really stable and no crack should occur. 2. Case ϑ > σ c /Ea. The elastic response is never stable and can never be observed. In such a case, cohesive cracks necessarily exist at any t > 0.

3 Energy minimization in the case of periodic crack distribution

Main assumptions

Knowing that cracks nucleate at t = 0 when ϑ > σ c /Ea, the scenario of crack nucleation for t close to 0 is investigated in this section. To this purpose, we consider the Hypothesis 1 All cracks are parallel to the x 2 -direction, have the same length and are periodically spaced.

We are seeking out their optimal distribution at a time t close to 0, i.e. the configuration which leads to the least energy.

In this scenario the displacement field is periodic and its jump across the crack is purely normal. Accordingly we can only consider kinematically admissible displacements of the same form and, still by symmetry, one has to study only one half of a periodic cell. Without loss of generality, we can suppose that the unit (half-)cell is the semi-infinite strip Y d = (0, d) × R + , d being the half-period, see Fig. 3. The side where the crack is located is Γ 0 = {0} × (0, ∞). The opposite side of the cell is Γ d = {d} × (0, ∞). The total mechanical energy (4) becomes a functional of d and of the (virtual) displacement field v at a given time t. Accordingly, taking into account the periodicity, one sets

E t (d, v) := 2L d Y d ψ(ε(v)-ε th t ) dx+ L d Γ0 κ(2v 1 ) dx 2 . ( 14
)
Owing to the symmetry, the kinematical conditions now read

v 1 ≥ 0 on Γ 0 , v 1 = 0 on Γ d ,
and the jump set of v is

S v = {x = (0, x 2 ) : v 1 (x 2 ) > 0} ⊂ Γ 0 .
Fig. 3 The unit cell and its boundary conditions .

Dimensionless problem statement

In order to simplify theoretical considerations, we first formulate the above described problem in a dimensionless setting. To this end we introduce the following transformations

Y d → Y d, x → y = x 2 √ kt , d → d = d 2 √ kt (15)
for coordinates and distances in space, at a given time whereas the displacements are mapped by

v(x) → v(y) = v(x) 2 √ kt aϑ . ( 16
)
Accordingly, this leads to define the normalized strain field by

ε(v) = 1 2 ∇ y v + ∇ T y v = 1 aϑ ε(v),
the normalized elastic strain energy function by

ψ(ε e ) = ψ(ε e ) a 2 ϑ 2 E = 1 2 εe 11 (v) 2 + εe 12 (v) 2 + 1 2 εe 22 (v) 2 (17)
and the normalized stresses by σij = σ ij Eaϑ .

The intensity of the thermal shock is characterized by the dimensionless parameter:

η = σ c Eaϑ ,
whereas the time t is replaced by the following dimensionless parameter:

t = 2 √ kt aϑ δ c .
Using the above definitions, the elastic response can be written in a dimensionless form as

σel 11 = f c (y 2 ), σel 12 = 0, σel 22 = 0, ( 18
)
εel 11 = 0, εel 12 = 0, εel 22 = -f c (y 2 ), ( 19
)
which corresponds to the following dimensionless displacement field

ūel (y) = +∞ y2 f c (ξ)dξ e 2 . ( 20
)
Finally the total mechanical energy functional (14) becomes

E t (d, v) = 4EL √ kt a 2 ϑ 2 Ēt( d, v) with Ēt( d, v) = 1 2 d Y d (ε(v) + f c (y 2 )I) • (ε(v) + f c (y 2 )I) dy + η d Γ0 κ(2 tv 1 ) 2 t dy 2 . ( 21
)
Therefore the problem consists now in minimizing Ēt with respect to d and v at given t. This dimensionless formulation shows that the problem essentially depends on two parameters: the dimensionless intensity of the thermal shock η and the dimensionless time t. It is important to see that the time has disappeared from the bulk energy to reappear in the surface energy. This will simplify the analysis.

Minimization of the energy at short times

To study the behavior for t close to 0, one formally passes to the limit in (21) when t goes to 0, at given ( d, v). Since limt →0 κ(2 tv)/2 t = κ (0)v 1 = v1 , that leads to consider the following energy functional

Ē0 ( d, v) = 1 2 d Y d (ε(v) + f c (y 2 )I) • (ε(v) + f c (y 2 )I) dy + η d Γ0 v1 dy 2 . ( 22
)
This energy is a quadratic strictly convex functional of v at given d. Therefore, there exists a unique displacement field which minimizes the energy over the convex set C( d) of kinematically displacement fields, with

C( d) = {v : v1 ≥ 0 on Γ 0 = {0} × (0, ∞), v1 = 0 on Γ d = { d} × (0, ∞)}.
As far as the minimization of Ē0 ( d, v) with respect to ( d, v) is concerned, we will use the following property:

Property 1 Let d, η and f be positive constants. Then

min ϕ∈L 2 (0, d) d 0 ϕ(y1) dy1≤0 1 2 d d 0 (ϕ(y 1 ) + f ) 2 -2ηϕ(y 1 ) dy 1 = f 2 /2 -(f -η) 2 /2 if f ≥ η f 2 /2 if f < η .
Proof Decomposing ϕ into its mean value ϕ and a function with zero mean value, i.e.

ϕ(y

1 ) = ϕ + φ(y 1 ), ϕ ≤ 0, d 0 φ(y 1 ) dy 1 = 0,
the functional to minimize reads as

1 2 d d 0 (ϕ(y 1 ) + f ) 2 -2ηϕ(y 1 ) dy 1 = 1 2 ( ϕ + f ) 2 -η ϕ + 1 2 d d 0 φ(y 1 ) 2 dy 1 .
Minimizing with respect to φ gives φ = 0 as the minimizer. Therefore, the problem consists in minimizing

1 2 ( ϕ + f ) 2 -η ϕ with respect to ϕ ≤ 0. One eas- ily obtains that the minimizer is ϕ = 0 if η > f and ϕ = η -f if η ≤ f . The property follows.
Let us return to the minimization problem of Ē0 ( d, v) by remarking first that the energy can read as

Ē0 ( d, v) = 1 2 d Y d v1,1 + f c (y 2 ) 2 -2ηv 1,1 + v2,2 + f c (y 2 ) 2 + 1 2 (v 1,2 + v2,1 ) 2 dy, ( 23
)
where vi,j stands for ∂v i /∂y j . Moreover, since v1 = 0 on Γ d, the condition v1 ≥ 0 on Γ 0 is equivalent to

d 0 v1,1 (y 1 , y 2 ) dy 1 ≤ 0, ∀y 2 > 0.
In the case where v = ūel an easy calculation based on (20) gives

Ēel := Ē0 ( d, ūel ) = ∞ 0 1 2 f c (y 2 ) 2 dy 2 , ∀ d > 0.
Let us now distinguish the two cases η ≥ 1 and η < 1.

1. Case η ≥ 1, i.e. ϑ ≤ σ c /Ea. Since f c (y 2 ) ≤ η for all y 2 > 0, Property 1 gives 1 2 d Y d v1,1 + f c (y 2 ) 2 -2ηv 1,
. Case η < 1, i.e. ϑ > σ c /Ea. Defining ¯ c = f c -1 (η) one has f c (y 2 ) > η if y 2 < ¯ c , f c (y 2 ) < η if y 2 > ¯ c .
Then, for any d > 0 and any v ∈ C d, Property 1 gives

1 2 d Y d v1,1 + f c (y 2 ) 2 -2ηv 1,1 dy ≥ Ēel -1 2 ¯ c 0 (f c (y 2 ) -η) 2 dy 2 ,
and the equality holds if and only if

v1 (y) = (f c (y 2 ) -η) + ( d -y 1 ),
(where x + denotes the positive part of x). Inserting into (23), one deduces that

∀ d > 0, ∀v ∈ C d: Ē0 ( d, v) ≥ Ēel -1 2 ¯ c 0 (f c (y 2 ) -η) 2 dy 2 .
Therefore we have obtained that Ēel -

1 2 ¯ c 0 (f c - η) 2 dy 2 is a lower bound for Ē0 . Let us now prove that it is the infimum, i.e. inf d>0 min v∈C d Ē0 ( d, v) = Ēel -1 2 ¯ c 0 (f c (y 2 ) -η) 2 dy 2 .
Let us first remark that this lower bound cannot be a minimum, i.e. it cannot be reached by any pair ( d, v) with d > 0. Indeed, in order that the lower bound be reached, one should have both v1 (y) = (f c (y 2 )-η) + ( d-y 1 ), v2,2 +f c (y 2 ) = 0 and v1,2 +v 2,1 = 0, which is impossible. In order to prove that the lower bound is the infimum, one must construct a sequence whose energy converges to the lower bound. So, let us consider d = 1/n with n ∈ N * and vn

(y) = (f c (y 2 ) -η) + 1 n -y 1 e 1 + ∞ y2 f c (ξ) dξ e 2 , for y ∈ (0, 1/n) × (0, ∞). Inserting into (23) gives Ē0 (1/n, vn ) = Ēel -1 2 ¯ c 0 (f c (y 2 ) -η) 2 dy 2 + 1 12n 2 ¯ c 0 f c (y 2 ) 2 dy 2 .
Passing to the limit when n goes to infinity gives the desired result.

Conclusion

The analysis of the previous subsection leads to the following conclusion concerning the minimization problem of the limit energy at short times:

1. If the thermal shock is small enough, i.e. if η := σ c /(Eaϑ) ≥ 1, then the elastic response minimizes the (limit) energy (at least in the class of displacement fields considered in Hypothesis 1), result which confirms the analysis made in Subsection 2.2. 2. On the other hand, if the thermal shock is large enough, i.e. if η < 1, then the infimum of the (limit) energy cannot be reached by any finite spacing of the cracks but is approached by the following minimizing sequence ( d, ū d) when the crack spacing d goes to 0:

ū d(y) = (f c (y 2 ) -η) + ( d -y 1 ) e 1 + ∞ y2 f c (ξ) dξ e 2 ,
where x + denotes the positive part of x. Accordingly, the "limit" when d tends to 0 can be seen as a distribution of infinitely close cohesive cracks of length ¯ c = f c -1 (η) and across which the displacement jump is infinitely small. Specifically, one gets

lim d→0 ū d (y 2 ) d = 2(f c (y 2 ) -η)e 1 if y 2 ≤ ¯ c 0 otherwise .
Let us recall that these results have been obtained by minimizing not the true dimensional energy functional Ēt but only its limit Ē0 when t goes to 0.

Numerical confirmation of the minimization of the total mechanical energy

The minimum of the total mechanical energy predicted in Section 3.3 is confirmed by finite element simulations. The dependence of the total mechanical energy of the system on the periodicity d (Fig. 3) is studied with the commercial code COMSOL for different values of η. In the finite element model the height H of the domain in y 2 -direction is chosen large enough to make sure that at y 2 > H the changes in the stress field due to crack formation are negligible. A structured mesh of quadrilateral linear finite elements with edge length min( d, H)/30 is used. The cohesive length ¯ c for a certain d is such that

σ11 (0, y 2 ) = η if y 2 ≤ ¯ c , σ11 (0, y 2 ) < η if y 2 > ¯ c with ū1 (0, y 2 ) = 0 for y 2 > ¯ c .
Owing to the non-linearity, ¯ c ( d) cannot be computed a priori and the boundary condition on y 1 = 0 is not straightforward. To deal with this difficulty we prescribe a symmetry boundary condition on {y 1 = 0, y 2 > ¯ max c } where ¯ max c is chosen large enough to ensure that it is always larger than ¯ c . On {y 1 = 0, y 2 ≤ ¯ max c } we impose as a nonlinear Neumann boundary condition in y 1 -direction the stress

η reg = η if ū1 ≥ ūtol 1 ηū 1 /ū tol 1 otherwise (24)
where ūtol 1 d. Choosing ūtol 1 sufficiently small ensures that the boundary is either subject to the cohesive stress η or satisfies ū1 < ūtol 1 ≈ 0. It is made sure that u tol 1 is chosen small enough to introduce a numerical error not greater than the one resulting from other sources such as the finite discretization length of the finite element mesh. For the simulations we chose in (24) u tol 1 = 0.25 d × 10 -4 for η = 0.9 and u tol 1 = 1.25 d × 10 -4 for η = 0.5. In Fig. 4 the convergence of min v∈C d Ē0 ( d, v) in finite element simulations towards the infinum predicted in Section 3.3 for d → 0 is depicted for η = 0.9 and η = 0.5. Excellent agreement between the numerical results and the analytical predictions is observed underlining especially that for d → 0 indeed a unique minimum of total mechanical energy is reached.

Not only do these computer simulations confirm the theoretical analysis in Section 3.3, but they can also be exploited to illustrate the mechanical meaning of the results. Thus in Fig. 5 the stress field σ11 observed in the finite element simulations with η = 0.9 is plotted for different values of d. -For d = 0.01 obviously σ11 ≈ η in the whole part above the crack tip, i.e. for y 2 < ¯ c = f c -1 (.9) = 0.089 as expected.

-For increasing values of d, however, we observe in the upper part of the domain that σ11 ≈ η only for small y 1 , whereas for increasing y 1 the effect of the cohesive crack diminishes smoothly. -The stress increases and becomes more and more similar to the thermal stress field from ( 18). -For d = 1 the cohesive stress is observed only in a comparatively small part of the domain whereas for y 1 → d the stress field becomes approximately the one of an uncracked body.

Precisely, for any finite d we expect from Saint Venant's principle that in the upper part of the domain there are sections parallel to the y 2 -axis along which σ11 > η. Obviously, the creation of a new cohesive cracks along such sections is energetically favorable which means that for any finite d we expect that a subdivision of the considered unit cell by an additional crack should happen. This immediately leads to the conclusion that d → 0 is the energetically optimal periodicity between adjacent cracks.

4 Crack initiation by formation of a micro-cracked layer

Effective behavior of micro-cracks

The theoretical analysis above suggests that, at a macroscopic point and for a given local macroscopic strain, a configuration consisting of a layer of cohesive microcracks could be energetically favored over the elastic solution. That depends on the intensity of the macroscopic strain. Specifically, the best deformation in order to minimize the (total) energy consists in: (i) either a pure elastic response without microcrack, if the strain component ε e 11 is small enough; (ii) or a piecewise elastic response intersected by a fine arrangement of vertical cohesive micro-cracks across which the displacement field suffers an infinitesimal jump, otherwise. This type of microstructuration is well known in phase transform problems where the the energy is not convex [START_REF] Ball | Fine phase mixtures as minimizers of energy[END_REF]; [START_REF] Bhattacharya | Relaxation of some multi-well problems[END_REF]. The associated procedure of energy relaxation can be transposed to the present context [START_REF] Bourdin | The variationnal approach to fracture[END_REF]; [START_REF] Bouchitté | A global method for relaxation[END_REF]. From the physical viewpoint, this consists in replacing the original energy functional by that corresponding to the homogenized (effective) behavior associated with this optimal microcracking structuration. Of course, the construction of this effective behavior must take into account the history of the macroscopic strain because the current response depends on the presence of the previously created micro-cracks. However, since the irreversibility condition into a cohesive force model can be introduced in many different ways (whereas there exists a quite natural way in Griffith's theory setting where the cohesive forces are neglected), the definition of the effective behavior will be made in two steps: first, by considering the so-called reversible case where there is no condition of irreversibility; then, by introducing a specific irreversibility condition. In both cases, since we are merely interested here by the thermal shock problem, we have only considered the possibility of micro-cracks parallel to the x 2 direction. Accordingly, the behavior in the x 1 direction only is affected by the micro-cracks and we can consider that the representative volume element is in fact a one-dimensional interval which after a suitable rescaling is assumed to be Y = (-1/2, 1/2).

-The reversible case. Let ε e be the given macroscopic elastic strain tensor. Then, by definition, the effective strain energy ψ r is defined by minimizing, at the microscale, the total mechanical energy with respect to all possible micro-crack distributions at that given macroscopic strain. Specifically, that leads to the following definition:

ψ r (ε e ) = min v1∈C Y \Sv 1 ψ(ε e + v 1 (y 1 )e 1 ⊗ e 1 ) dy 1 + y1∈Sv 1 σ c v 1 (y 1 ) , (25)
where

C = {v 1 : v 1 (±1/2) = 0, v 1 ≥ 0 on S v1 }.
In (25), ψ is the original elastic potential given by (1) and S v1 denotes the set of points of the interval (-1/2, 1/2) where the (dimensionless) scalar field v 1 is discontinuous. This set can be empty and one has to find the (possibly discontinuous) optimal field v 1 such that the total mechanical energy be minimal, at given macroscopic strain. Note that the effective energy is the sum of the bulk elastic energy and the surface energy due to the micro-cracks and this latter energy only involves σ c = κ (0) because the jump of the displacement is necessarily infinitesimal. Indeed, finite jumps would give rise to an infinite energy when summed over all the micro-cracks. This minimization problem is essentially the same as that considered in Subsection 3.3 and hence we simply give the final result. The effective potential ψ r is given by:

ψ r (ε e ) =                    E 2 ε e 11 2 + ε e 22 2 + 2ε e 12 2 if ε e 11 ≤ σ c /E σ c ε e 11 - σ 2 c 2E + E 2 ε e 22 2 + 2ε e 12 2 if ε e 11 > σ c /E .
Thus, the dependence on ε e 11 only is affected. The quadratic dependence is replaced by an affine dependence when ε e 11 becomes larger than σ c /E so that the effective stress σ 11 remains constant and equal to σ c . Note that the minimizer is v 1 = 0 when ε e 11 ≤ σ c /E and hence is not discontinuous. On the other hand, when ε e 11 > σ c /E, there exist an infinite number of minimizers which differ only by their jump set. Specifically, in this case, any minimizer has to satisfy

v 1 (y 1 ) = σ c E -ε e 11 in Y \ S v1 , v 1 > 0 on S v1 , y1∈Sv 1 v 1 (y 1 ) = ε e 11 - σ c E .
Thus, the number and the position of the discontinuity points are arbitrary (at the micro-scale). For instance, in the case where the minimizer contains a single discontinuity at y 1 = 0 then v 1 reads as

v 1 (y 1 ) = σc E -ε e 11 (y 1 + 1 2 ) if -1 2 < y 1 < 0 σc E -ε e 11 (y 1 -1 2 ) if 0 < y 1 < 1 2 .
-The irreversible case. In the construction above, the effective behavior only depends on the current strain tensor ε e and not on its history. In particular, the opening and the closure of the micro-cracks are considered as a reversible process. It is not realistic and one has to introduce some irreversibility conditions. We adopt here the conditions proposed in Jaubert and Marigo (2006); Abdelmoula et al. ( 2010) and [START_REF] Talon | A model of adhesion coupled to contact and friction[END_REF]. Specifically, the rate of the displacement jump across the micro-cracks is governed by the following flow rule:

v1 ≥ 0 if σ 11 = σ c , v1 = 0 if 11 ∈ (0, σ c ) v1 ≤ if σ 11 = 0. ( 26
)
Thus, the displacement jump increases only when the stress is at the critical value σ c and decreases only when the stress vanishes. Consequently, the total (dimensionless) opening δ, (sum of the) displacement jump over the (set of) microcrack(s), can be seen as an internal variable δ which enters in the stress-strain relation

σ 11 = E(ε e 11 -δ), δ := y1∈Sv 1 v 1 (y 1 ) ≥ 0.
Using (26), the evolution law of δ can be formulated in terms of the following Kuhn-Tucker relations:

when δ = 0,      δ ≥ 0, σ 11 -σ c ≤ 0 (σ 11 -σ c ) δ = 0 ; (27) when δ > 0, 0 ≤ σ 11 ≤ σ c σ 11 δ = σ c δ+ . (28)
In ( 27), δ+ denotes the positive part of δ, i.e. δ+ = max{0, δ}. Accordingly, the effective behavior of the material can be seen as that of an elastic-plastic material where δ represents the (unique non vanishing) plastic strain component. This plastic strain exists only when the material is not in compression and it evolves only when the tensile stress is at the ends of the interval [0, σ c ]. That leads to a behavior which is schematized on Figure 7 where the arrows indicate the direction in which the paths must be followed (a double arrow means that the path is reversible). The effective strain energy ψ r is now function not only on the current strain tensor ε e , but also on the current total opening δ and even on the cumulated opening ∆ defined at time t by

∆(t) = t 0 δ(s) + ds.
Specifically, the effective strain energy ψ r reads as

ψ r (ε e , δ, ∆) = E 2 (ε e 11 -δ) 2 + ε e 22 2 + 2ε e 12 2 + σ c ∆. (29)
For any accessible state (ε e , δ, ∆), ψ r enjoys the following inequality which plays a fundamental role to obtain stability properties:

ψ r (ε e , δ + hδ * , ∆ + hδ + * ) ≥ ψ r (ε e , δ, ∆), (30)
∀h > 0 and

∀δ * ≥ 0 if δ = 0 ∀δ * if δ > 0 .
Indeed, an easy calculation gives

ψ r (ε e , δ + hδ * , ∆ + hδ + * ) -ψ r (ε e , δ, ∆) = h(σ c δ + * -σ 11 δ * ) + h 2 2 Eδ 2 * .
If δ = 0, then σ 11 ≤ σ c and the inequality (30) follows. If δ > 0, then 0 ≤ σ 11 ≤ σ c and the inequality (30) still holds.

Thus in (30) a second order term is introduced which will allow the lost of stability.

The response corresponding to the propagation of a micro-cracked layer

Throughout the present subsection we only consider the case where Eaϑ > σ c . The previous analysis shows that a possible solution of the evolution problem corresponds to the growth, from the upper boundary, of a micro-cracked layer the depth of which continuously increases with time. Let us precisely define such an evolution by discriminating between the reversible and the irreversible cases.

1. The reversible case. The total mechanical energy of the body at time t associated with the kinematically admissible displacement field v is defined by replacing the elastic potential ψ by the effective strain energy ψ r . Therefore it now reads

E t (v) = Ω\Sv ψ r (ε(v) -ε th t ) dx + Sv κ( v 1 ) dx 2 .
The real displacement and stress fields u t and σ t at time t must satisfy the first order stability condition which differ only from ( 9)-( 13) in the stress-strain relation which becomes

σ t = ∂ψ r ∂ε (ε(u t ) -ε th t ).
It is easy to check that the elastic displacement field u el t given in Subsection 2.2 satisfies the new first order stability conditions. The associated stress field now reads

σ t (x) = σ c if x ∈ Ω c t Eaϑf c x2 √ 2kt if x ∈ Ω e t ,
where Ω c t and Ω e t denote respectively the microcracked layer and the uncracked part of the body at time t, i.e.

Ω c t = (0, L) × (0, t ), Ω e t = (0, L) × [ t , ∞),
t being the depth of the micro-cracked layer

t = ¯ c √ 2kt, ¯ c = f c -1 (η).
2. The irreversible case. The displacement and stress fields are the same as in the reversible case. The unique difference is the existence of two internal variables: the opening field δ (which plays the role of a residual strain) and the cumulated opening field ∆. At each time t these fields are equal because there is no unloading, and vanish outside the microcracked layer. Inside the layer, they are given by

δ t (x) = ∆ t (x) = aϑf c x 2 √ 2kt - σ c E if x ∈ Ω c t .

Stability of the micro-cracked layer

We have shown in the previous subsection that growth of a micro-cracked layer without any macro-crack is a possible solution to the evolution problem. It remains to check the stability of this micro-cracked layer, which we will do in the following subsection, where for the sake of simplicity only the irreversible case is considered.

Definition of the stability in the irreversible case

Since the stability condition has only been defined in the reversible case, one must redefine this condition in the irreversible case, where the principle remains the same: comparing the energy of the body in its tested state with the energy that the body would have if one perturbs this state in any possible direction. Specifically, let t > 0 be a given time and (u el t , δ t , ∆ t ) be the state of the body at this time corresponding to the micro-cracked layer. Recalling that ∆ t = δ t , the energy of the body in this state reads as

E t (u el t , δ t , ∆ t ) = Ω ψ r ε(u el t ) -ε th t , δ t , δ t dx,
where ψ r is the effective energy density given by ( 29).

Let v be a kinematically admissible displacement field (which can be discontinuous across its jump set S v ) and let δ * be a direction of perturbation of the opening field, δ * can be any square integrable function over Ω with the constraint that δ * must be non-negative outside the micro-cracked layer. If one perturbs the body at time t by changing u el t to u el t + hv, δ t to δ t + hδ * and hence ∆ t = δ t to δ t + hδ + * where h is a positive number, then the energy of the body becomes

E t (u el t + hv, δ t + hδ * , ∆ t + hδ + * ) = = Ω\Sv ψ r ε(u el t + hv) -ε th t , δ t + hδ * , δ t + hδ + * ) dx + Sv κ(h v 1 ) dx 2 .
The condition for stability is that for any admissible direction of perturbation (v, δ * ), the following inequality

E t (u el t + hv, δ t + hδ * , ∆ t + hδ + * ) ≥ E t (u el t , δ t , ∆ t )
holds if h is sufficiently small. Expanding with respect to h up to the second order gives

0 ≤ h Ω\Sv σ t • (ε(v) -δ * e 1 ⊗e 1 ) + σ c δ + * dx + h Sv σ c v 1 dx 2 + h 2 2 Ω\Sv E(ε(v) -δ * e 1 ⊗e 1 )•(ε(v) -δ * e 1 ⊗e 1 ) dx + h 2 2 Sv κ (0) v 1 2 dx 2 + o(h 2 ). ( 31
)
Using the fact that σ t = σ t11 e 1 ⊗ e 1 , the first order stability condition reads as

0 ≤ Sv σ c -σ t11 ) v 1 dx 2 + Ω\Sv (σ c δ + * -σ t11 δ * ) dx.
The inequality above is automatically satisfied, because 0 ≤ σ t11 ≤ σ c and v 1 ≥ 0 everywhere. Moreover, the inequality becomes an equality if and only if the direction of perturbation is such that v 1 = δ * = 0 outside the micro-cracked layer and δ * ≥ 0 inside the microcracked layer. Considering such directions of perturbation, the second order stability condition reads as

0 ≤ Ω\Sv E(ε(v) -δ * e 1 ⊗e 1 )•(ε(v) -δ * e 1 ⊗e 1 ) dx + Sv κ (0) v 1 2 dx 2 . ( 32
)
Since the above inequality must hold for any δ * such that δ * = 0 in Ω e t and δ * ≥ 0 in Ω c t , it must hold for the field δ * which minimizes the right-hand side of (32) at given v. An elementary calculation gives that this minimizer is

δ * (x) = v + 1,1 (x) if x ∈ Ω c t \ S v 0 if x ∈ Ω e t ,
where the plus still denotes the positive part. Inserting into (32) the inequality becomes

Ω c t \Sv E (v - 1,1 ) 2 + (v 2,2 ) 2 + 1 2 (v 1,2 + v 2,1 ) 2 dx + Ω e t Eε(v) • ε(v) dx ≥ -κ (0) Sv v 1 2 dx 2 , ( 33
)
where a -denotes the negative part of a, i.e. a -= max{0, -a}. The inequality (33) must hold for any v such that S v ⊂ Ω c t and v 1 ≥ 0 on S v . It is automatically satisfied in the case of Dugdale's model because κ (0) = 0. Therefore, the response associated with the growth of a micro-cracked layer is always stable for Dugdale's model.

Let us examine now the general case where κ (0) < 0. Introducing the Rayleigh ratio

R t (v) := E |κ (0)| N t (v) Sv v 1 2 dx 2 , ( 34
)
where the quadratic form N t reads

N t (v) := Ω c t \Sv ε(v) • ε(v) -(v + 1,1 ) 2 dx + Ω e t ε(v) • ε(v) dx.
The stability condition consists in comparing the infimum of the Rayleigh ratio (over the set of admissible v) with 1. Specifically, setting

R t := inf v : v1 ≥0 on Sv⊂Ω c t , v1=0 on x1=±L R t (v), ( 35
)
one gets if R t > 1 the micro-cracked state is stable at t, if R t < 1 the micro-cracked state is unstable at t.

Minimization of the Rayleigh ratio

The Rayleigh ratio depends on time t only through the depth t of the micro-cracked layer, and so does also its infimum. Moreover, the infimum depends also on the width L of the body and on material parameters. By simple arguments of physical dimension, it is easy to check that R t can read as

R t = E |κ (0)|L R t L . (36)
In ( 35), R is a dimensionless positive function associated with the following rescaled Rayleigh ratio minimization problem:

R(λ) = inf v∈C(λ) R(λ, v), (37)
which is posed on a domain of width 2 containing a layer of depth λ = t /L. Specifically, one sets

C(λ) := {v : v 1 ≥ 0 on S v ⊂ Ω c (λ), v 1 = 0 on x 1 = ±1}, R(λ, v) := N (λ, v) Sv v 1 2 dx 2 , ( 38
) N (λ, v) := Ω c (λ)\Sv ε(v) • ε(v) -(v + 1,1 ) 2 dx + Ω e (λ) ε(v) • ε(v) dx
and

Ω c (λ) := (-1, 1)×(0, λ), Ω e (λ) := (-1, 1)×(λ, +∞).
Let us prove that λ → R(λ) is not increasing, i.e.

R(λ

1 ) ≥ R(λ 2 ) if λ 1 < λ 2 . (39) First, since Ω c (λ 1 ) ⊂ Ω c (λ 2 ), one gets C(λ 1 ) ⊂ C(λ 2 ) and R(λ 1 , v) ≥ R(λ 2 , v) for all v ∈ C(λ 2 ). Therefore, inf v∈C(λ2) R(λ 2 , v) ≤ inf v∈C(λ2) R(λ 1 , v) ≤ inf v∈C(λ1) R(λ 1 , v)
and ( 39) follows.

Let us construct an upper bound for R(λ). Considering the following displacement field

v(x) = sin π 2 1 - x 2 λ + sign(x 1 ) sinh π(1 -|x 1 |) 2 √ 2 λ e 1
which belongs to C(λ) and whose jump set is S v = {0} × (0, λ), a tedious calculation eventually leads to the following estimate

R(λ) ≤ R(λ, v) = π 4 √ 2 λ tanh π 2 √ 2 λ . ( 40
)
Fig. 8 Normalized displacement test field for the evaluation of the upper bound.

To obtain a lower bound, one uses the positivity of each term of the bulk energy to get

R(λ) ≥ inf v∈C(λ) Ω c (λ)\Sv (v - 1,1 ) 2 dx Sv v 1 2 dx 2 . ( 41
)
The infimum above is equal to 1/2 (the proof is left to the reader) and is reached by the field v which vanishes in Ω e (λ) and is given by

v(x) = sign(x 1 )(1 -|x 1 |)e 1 in Ω c (λ).
(This displacement field is not in C(λ) because it is discontinuous at x 2 = λ, but it can be approached by elements of C(λ).) Finally, we have obtained the following bounds for the infimum of the Rayleigh ratio:

Eδ c 2ωσ c L ≤ R t ≤ π 4 √ 2 tanh πL 2 √ 2 t Eδ c ωσ c t , ( 42
)
where ω = -κ (0) > 0. Let us remark that the lower and upper bounds are equal when t (and hence t ) go(es) to infinity. This mathematical result can be interpreted mechanically as pointed out in remark 2. These estimates are sufficient to discuss the stability of the micro-cracked layer at large times, but must be improved at short times. Indeed, let us define the characteristic length of the material

d c = Eδ c ωσ c . ( 43
)
If the body width L is small enough, i.e. 2L ≤ d c , one deduces from ( 42) that R t > 1 for all t and hence the micro-cracked layer is always stable. On the other hand, if 2L > d c , then R t < 1 for sufficiently large times t so that the micro-cracked layer becomes unstable. Still in case that 2L > d c , to conclude on the stability of the micro-cracked layer at small times requires to find lim λ→0 R(λ). The determination of this limit is a real issue and we make the following conjecture which is based on some arguments and the computations presented on the next subsection:

Hypothesis 2 The function λ → R(λ) goes to infinity like 1/λ when λ goes to 0. Specifically, there exists R 0 > 0 such that lim

λ→0 λ R(λ) = R 0 .
Note that the constant R 0 is necessarily independent of the material parameters and of the thermal shock.

It can only depend on the shape of the body and on the type of boundary conditions which are prescribed.

Next, we assume that R(λ) is strictly monotonically decreasing in order to ensure the existence of the inverse function R -1 for all λ (we have already proved that R(λ) is not increasing):

Hypothesis 3 The function λ → R(λ) is strictly monotonically decreasing.

Finally, we assume similarly as already for the microcracks Hypothesis 4 Macro-cracks are parallel to the x 2 -direction and form a periodic pattern in x 1 -direction.

From Hypothesis 2, 3 and 4 as well as (36) we conclude that Property 2 As far as the stability of the micro-cracked layer is concerned, one has 1. If 2L ≤ d c , then the micro-cracked layer is always stable; 2. If 2L > d c , then there exists a critical time t c before which the micro-cracked layer is stable and after which the micro-cracked layer is unstable. This critical time is given by

t c = R -1 (L/d c ) erfc -1 (η) 2 L 2 4k . (44)
Therefore, the stronger the thermal shock (i.e. the smaller η) or the smaller the material characteristic length d c , the shorter the critical times.

3. If 2L d c , then the critical time is approximately given by

t c ≈ t ∞ c = R 0 erfc -1 (η) 2 d 2 c 4k
and corresponds to the time when the depth of the micro-cracked layer is tc ≈ R 0 d c .

Note that for large L, t c is a finite, positive value independent on L, which means that the critical time for destabilization of the micro-crack layer never drops under a minimal value. Therefore, regardless of the body width the micro-crack layer will always be stable. Figure 9 plots the critical time versus the load parameter η, and thus t c is increasing with the thermal load parameter η. The results from Section 4.3.2 are checked numerically by finite element computations of the rescaled Rayleigh coefficient R(λ) for different values of the ratio λ between thermal penetration depth and cell width.

This Rayleigh coefficient in the irreversible case can be determined only numerically. From ( 37) and ( 38) we conclude that for a given ratio λ between thermal penetration depth and domain width, the rescaled Rayleigh ratio R(λ) is the smallest real value satisfying for any

variation v the condition 0 = Ω c (λ)\Sv (v - 1,1 ) 2 + (v 2,2 ) 2 + 1 2 (v 1,2 + v 2,1 ) 2 dx + Ω e (λ)\Sv ε(v) • ε(v) dx -R(λ) Sv v 1 2 dx 2 . ( 45
)
To compute this value, we perform a finite element computation for reasons of symmetry only on the half-cell with 0 ≤ x 1 ≤ 1. This half-cell is discretized by a rectangular grid and we use bilinear basis functions and a constant discretization length h 2 in y 2 -direction on the boundary Γ . The column vector V I concatenates all x 1 -displacements of nodes situated on the crack boundary S v , which forms a part of the boundary of the half-cell considered. Note that the nodal displacements V I on the crack boundary in our finite element computation distinguish from the jump v 1 in (45) by a factor of two. Therefore (45) can be written in a discrete manner as

0 = V I V II T      K el I,I K el I,II K el II,I K el II,II K el -2R(λ)h 2 1 0 0 0      K el,crack V I V II (46)
where 1 denotes the identity matrix, h 2 the discretization length in x 2 -direction. K el is the elastic stiffness matrix of the considered half-cell with zero Dirichlet boundary conditions everywhere except for S v and the body surface at x 2 = 0, which are free boundaries. K el,crack can be seen as the total stiffness of the complete system consisting of an elastic continuum and a cohesive crack. As ( 46) can be satisfied if and only if V I and V II form an eigenvector of K el,crack with eigenvalue zero,

V II = -(K el II,II ) -1 K el II,I V I . (47)
Let us denote µ i the eigenvalues of

K el,ef f = K el I,I -K el I,II (K el II,II ) -1 K el II,I
Recalling that R(λ) is the smallest real value satisfying (46), we arrive at

R(λ) = 1 2h 2 µ 1 (48)
Mechanically, in (48) K el,ef f describes the effective stiffness of the elastic material in the periodic cell towards a crack opening in direction V I and its smallest eigenvalue normalized by the discretization length h 2 describes the minimal stiffness of the continuum towards a crack opening in any direction, which is exactly the mechanical meaning of the Rayleigh-like coefficient R(λ). Using (48), R(λ) can be computed conveniently in M AT LAB for any given λ. The results of these computations are depicted in Fig. 10. By parameter studies with both different domain depths in x 2 -direction and different discretization lengths it is verified that the numerical error by discretization and finiteness of the simulated domain in x 2 -direction is negligible in all simulations.

Obviously, in the limit 1/λ → 0, (i.e. t → ∞) we arrive at λR(λ) ∼ λ/2, i.e. at lim λ→∞ R(λ) = 1/2, which is exactly what we expect from ( 41) and ( 42) for large λ where upper and lower bound converge to the same value. On the other hand, for λ → 0 we obtain the constant lim λ→0 λR(λ) = R 0 ≈ 0.232, which backs up numerically Hypothesis 2. Note that this numerical result is also in accordance with the upper bound from (40) for small λ, which requires lim λ→0 λR(λ) ≤ π/(4 √ 2) ≈ 0.555. As can be seen in Fig. 10, R(λ) is indeed not only monotonically decreasing as pointed out already in (39), but even strictly monotonically decreasing as assumed in Hypothesis 3. Thus both the Hypothesis 2 and 3 are confirmed by the numerical computation of R(λ).

Remark 2 (Interpretation of R(λ) for λ → ∞) For very large values of λ the domain becomes a quasi onedimensional strip stretched out in x 2 -direction. In such a continuum shear deformation and strain in x 2 -direction are expected to be negligible so that v 2,1 = v 2,2 = 0 in (38) and the strain in x 1 -direction is expected to be constant. Geometric compatibility requires that the strain v - 1,1 sums up over the cell of width 2 to the jump v 1 :

Ω c (λ)\Sv v - 1,1 dx = 2 Sv v - 1,1 dx 2 = Sv v 1 dx 2 , ( 49
)
and thus in Ω c (λ), v - 1,1 = v 1 /2 so that (38) gives

lim λ→∞ R(λ) = 1/2. ( 50
)
Thus the fact that in (42) the upper and lower bounds converge for t → ∞, i.e. , λ → ∞ to the same value which is equivalent to R(λ) = 1/2, can be understood mechanically simply from the fact that in this limit the domain becomes quasi one-dimensional so that all strain components become negligible except for the one in direction of the width of the body, which becomes constant and thus directly related to the jump at the crack.

Comparison with Griffith's surface energy and gradient based damage models

Following Property 2, when 2L

d c t c ≈ R 0 erfc -1 (η) 2 d 2 c 4k , tc ≈ R 0 d c . (51)
It is instructive to compare the current results to those from a gradient-damage model for the same problem [START_REF] Sicsic | Initiation of a periodic array of cracks in the thermal shock problem: a gradient damage modeling[END_REF]. Qualitatively, in both cases no macroscopic cracks are initiated at t = 0. Rather for the damage model one observes a homogeneous (with respect to y 2 -direction) damaged layer growing in time which can be considered the counterpart to the micro-crack layer in the cohesive model. At a critical time t c for both models a loss of stability of this damaged layer is observed. In the damage model the critical time (which corresponds in fact to a bifurcation point) is according to (Sicsic et al., 2013, Proposition 10) given by

t b = τ 2 c (η) η 2 ¯ 2 dam 4k ( 52
)
where τ c (η) is a function of the thermal load parameter η determined numerically in [START_REF] Sicsic | Initiation of a periodic array of cracks in the thermal shock problem: a gradient damage modeling[END_REF], ¯ dam is a characteristic length of the material in the damage model. Obviously, this length affects t c the same way as the characteristic material length d c in the cohesive force model. Furthermore, the thermal diffusivity k plays the same role in ( 51) and ( 52). As to the influence of the thermal load there, interestingly, these both expressions are numerically very similar (Fig. 9).

5 From the loss of stability of the micro-crack layer to a parallel array of macro-cracks

For 2L > d c the micro-crack layer loses its stability at some time t c > 0. In the sequel we will restrict to case where the slab is such that 2L d c . In the following we will discuss the post-critical development following the principle of local minimum of total mechanical energy under the irreversibility constraint. Studying the full complexity of the mixture of micro-cracks and macrocracks goes beyond the scope of this paper. Therefore we limit the discussion of the post-critical developments to some simple qualitative arguments revealing the general lines along which the system will evolve after the loss of stability of the micro-cracked layer in the irreversible case.

Initiation of macro-crack pattern

In the previous section, the micro-layer is unstable at t = t c where a single macrocrack appears in the center of the domain. To understand how the pattern of macro-cracks will evolve in time, we assume for the moment that the macro-cracks are not formed directly at t = t c , but at some later time. We seek those of a periodicity d opt that minimize the energy per unit length.

From Hypothesis 1 and 4, the energy saved in a periodic cell of width d opt by formation of macro-cracks can be expressed by the second order term in h in (31) interpreting Ω and S v as domain and crack boundary in this periodic cell. As the term h 2 2 Sv κ (0) v 1 2 dx 2 does not depend on d opt , the minimization of this energy is with (34) and ( 35) equivalent to the minimization of R d t + 1 where R d t is referring to a periodic cell of width d instead of the whole domain of width L as R t . Then the energy per unit length and thus the total mechanical energy of the whole system is minimized by the crack spacing

d opt = arg min d R d t + 1 d ( 53
)
where we neglect the requirement that the body width L has to be a multiple of the optimal crack spacing d opt .

From the strict monotonicity of R d t in d we conclude that this minimum is obtained under the condition

0 = ∂ ∂d R d t + 1 d d=dopt . ( 54
)
Using the numerical data from Fig. 10, we can plot d opt on the basis of (54) as depicted in Fig. 11. The limit for t → ∞ observed in this figure can be derived as follows: for t → ∞ we conclude from ( 42) and ( 43 In Fig. 11 we can clearly see that the optimal crack spacing decreases. The requirement that the body width has to be a multiple of the crack spacing is neglected. If this constraint is additionally taken into account, the graph for d opt becomes a staircase-shaped approximation of the graph in Fig. 11.

Remark 3 This results is different from the one given by global minimization with Griffith surface energy [START_REF] Jenkins | Optimal spacing and penetration of cracks in a shrinking slab[END_REF] where the optimal spacing at the surface depends on the loading 0 . In the case of gradient based damage models [START_REF] Sicsic | Initiation of a periodic array of cracks in the thermal shock problem: a gradient damage modeling[END_REF] the optimal spacing is proportional to a material length ¯ dam but also depends on the loading. Here it only depends on the material parameter d c and not on the loading.

Let us now consider a body with a finite, but very large width L. The changes of the stress field by the creation of a macro-crack at some time t c (L) (where the width is made explicit in (44)) can be expected to mainly affect the neighborhood of the crack, i.e. the region around x 1 = x1 (the center of the periodic cell considered, i.e. where the macro-crack forms). Far away from this first macro-crack the micro-crack layer will therefore persist and evolve nearly independently on macro-crack. Therefore, as time progresses and the optimal spacing for macro-cracks has decreased sufficiently in the way illustrated in Fig. 11, the two halfdomains in positive and negative direction of the first, central macro-crack are expected to be subdivided by another macro-crack. This process of continuous subdivision is actually expected to continue as time progresses although under no circumstances the macrocrack spacing is expected to drop below the minimal value d c from (56). To back up these assumptions about the post-critical behavior we performed finite element simulations (described in Appendix A) whose results are illustrated in Fig. 12.

The stress is observed in these simulations to drop with the formation of macro-cracks under the critical stress level η in the whole periodic cell except on the tips of the macro-cracks. However, if the crack spacing is large compared to the thickness of the microcrack layer, and thus to the length of the newly formed macro-crack, the stress field close to the boundary of the periodic cell remains almost unaffected by the formation of the macro-crack. Since σ 11 /σ c ≤ 1 before formation of the macro-crack we thus expect σ 11 /σ c < 1 afterwards everywhere in the domain. Recalling (31), the periodic cell is thus directly after the formation of the macrocrack in a stable state. If crack spacing d is small, the stress level does not recover back to σ c , but even slightly further decreases in the region close to x 2 = 0 where σ 11 /σ c = 1 previously (cf. Fig. 12 and13). This prevents from any further destabilization once the crack spacing has dropped under a certain level, which we expect from our simulations to be even significantly larger than the minimal periodicity d c discussed above. The situation is different, however, for a large crack spacing as illustrated in Fig. 12 and 13. In this case, at least far from the macro-crack, the stress σ 11 has never dropped much, but only slightly below σ c , and as time progresses, stress rises due to temperature diffusion and the critical stress level is quickly recovered in a large region as illustrated in the left plot of Fig. 13. This allows for further destabilization, i.e. the formation of new macro-cracks with a smaller crack spacing. The monotonicity of the Rayleigh ratio with respect to the width of the periodic cells considered with the fact that the critical stress level is expected to be recovered first at the maximal distance from the already existing macro-cracks allows to assume that the reduction of crack-spacing happens by successive bisections of the periodic cells. This way crack spacing is expected to reduce until it has become too small for further destabilization of the existing periodic cells.

Remark 4 Although Fig. 13 provides a general understanding of how crack spacing reduces over time as long as the distance between the macro-cracks is still large enough, it is not appropriate for quantitatively predicting when the existing crack spacing is replaced by a smaller one. The reason is that it is based on the assumption of a micro-crack layer without preexisting macro-cracks whereas the subsequent stepwise reduction of crack spacing discussed above happens in an environment consisting of a mixture of micro-cracked regions and already existing macro-cracks. Our simulations assume that both cases differ even far away from the macro-cracks by certain subtle features such as the fact that the critical stress level is seemingly never fully recovered at any point with x 2 = 0 after the formation of the first macro-cracks. In all simulations it is observed that the macrocracks feature, directly after their formation, a large non-cohesive zone and only a small cohesive zone. This can be understood immediately from the variational approach used for the computation of the critical time for destabilization. Macro-cracks can form if there exists at least one deformation mode v for opening the macrocrack where the reduction of cohesive surface energy thanks to κ by an increased crack opening surpasses the increase in elastic energy by this deformation mode. If such a deformation mode v exists, the more energy can be released the further the crack opens as long as neither the stiffness of the continuum nor κ change. The former remains constant in linear continuum mechanics so that a maximal release of energy necessarily goes along with an opening of the macro-crack until κ changes at least on a part of the crack lip, i.e. until a non-cohesive zone exists (where κ drops to zero). Therefore, once the micro-crack layer loses stability and macro-cracks appear, these always appear by a "brutal" crack opening creating immediately a non-cohesive zone.

Propagation of macro-cracks

Once a minimal periodicity has been reached, which is no longer reduced by further bisections of the periodic cells, the critical stress will appear only at the tips of the macro-cracks. Therefore at this time, microcracks cannot form or grow any longer whereas the existing macro-cracks will keep growing as thermal diffusion progresses. The size of the micro-cracked regions will therefore become smaller and smaller and finally negligible relative to the length of the macro-cracks. Together with the negligibly small cohesive zone at the tip of the macro-cracks the situation becomes almost identical to a pattern of Griffith cracks and thus for large t, one may expect Griffith-like crack propagation, in which according to [START_REF] Bahr | Scaling behavior of thermal shock crack patterns and tunneling cracks driven by cooling or drying[END_REF] periodicity of crack propagation gradually doubles.

This way, we obtain for the irreversible case on the basis of the principle of local minimum of total mechanical energy a crack initiation and propagation: (1) initiation of a micro-cracked layer at t = 0; (2) propagation of this layer; (3) loss of stability at some critical time t c by a single macro-crack; (4) formation of a macro-crack pattern with very large periodicity; (5) this periodicity reduces gradually by bisection of the periodic cells until a minimal periodicity of the macrocrack pattern is reached; (6) propagation of the array; (7) selective arrest, the periodicity of the macro-crack pattern gradually doubles.

For the long term behavior of the array of cracks, the irreversibility conditions is neglected for simplification of the computations. In Fig. 14 the optimal periodicity for the reversible case is plotted for different values of the load parameter η versus time t, these results are compared with those of (Bahr et al., 2010, Fig. 4). For large t, the simulations performed in order to plot Fig. 14 revealed cracks where the length of the cohesive zone is negligible compared to the length of the non-cohesive zone (e.g., for η = 0.1 and 4kt/(R 0 d c ) 2 = 3.1 × 10 6 , the ratio of both length is 1.1 × 10 -3 ). Thus, in this case the cracks formed according to the cohesive model are largely similar to Griffith cracks and one would thus expect a similar influence of equivalent model parameters in both cases. This is indeed the case. The logarithmic curves depicted in (Bahr et al., 2010, Fig. 4) can be expressed as log(d/ 0 ) = A η + k log(a/ 0 ) (57)

with crack periodicity d, crack length a, a real parameter k. The normalization constant reads 0 = G c /(Eaϑ) 2 , which is related to d c by 0 = (2 -ω)ωd c η 2 /2. A η denotes a real constant for a given loading η. According to (Bahr et al., 2010, Eq. ( 18)) we furthermore have

a = 4kt 1.735 d , (58)
which together with (57) gives

d = A η × 1-k 1+k 0 = A η × η 2 1-k 1+k .
In Fig. 14 we observe for the largest value of t the relation d = A η × η 0.38 , which is equivalent to k = 0.68, which is exactly the slope in Fig. 4b of [START_REF] Bahr | Scaling behavior of thermal shock crack patterns and tunneling cracks driven by cooling or drying[END_REF]. If the exponent in the power law relating optimal periodicity and crack length in our simulations is almost identical to the one proposed in [START_REF] Bahr | Scaling behavior of thermal shock crack patterns and tunneling cracks driven by cooling or drying[END_REF], the prefactor in the power law differs by roughly a factor of two. The exact reason for this deviation has not yet been analyzed. Possible explanations are deviations between the cohesive and Griffith model in the regime close to the critical time t c where the cohesive cracks are initiated. It should also not be forgotten that this regime is close to the point a/d ∼ 1 where the scaling assumptions on which the analysis of [START_REF] Bahr | Scaling behavior of thermal shock crack patterns and tunneling cracks driven by cooling or drying[END_REF] is based are no longer satisfied.

Conclusion and Perspectives

In this paper we studied crack initiation in a half-plane cooling by means of a cohesive force model. A temperature diffusion independent of the crack set and a perfect conductivity at the surface are assumed. If the loading is sufficiently severe (η ≤ 1), cracks nucleate at t = 0 whereas if η > 1 the elastic solution remains stable for any time. For short times a layer of micro-cracks with infinitesimal periodicity and opening is formed. This allows to define an effective behavior which reveals to be very similar to the one of an ideally plastic constitutive law.

Assuming crack propagation on the basis of the principle of local minimum of the total mechanical energy subject to a irreversibility condition this layer is shown to remain stable and grow until the critical point of time t c . The critical stress σ c and the toughness G c are not enough to model crack initiation and propagation as our paper clearly demonstrates that the shape of the cohesive model, through κ (0) as a third free parameter, has an influence on the stability. Indeed, for Dugdale (ω = -κ (0) = 0) surface energies, this homogenized behavior remains stable for any time t > 0. Furthermore the irreversibility assumption is fundamental, as with no irreversibility Dal [START_REF] Maso | Gradient bounds for minimizers of free discontinuity problems related to cohesive zone models in fracture mechanics[END_REF] the micro-crack layer is found to always be unstable. Irreversibility was already found to be a crucial modeling assumption in order to reproduce fatigue effects [START_REF] Abdelmoula | Construction and justification of paris-like fatigue laws from dugdale-type cohesive models[END_REF]. The micro-cracked layer we observe initially may be viewed as the counterpart of a homogeneous damaged layer found for the same problem at small times theoretically and numerically on the basis of a gradient damage model in [START_REF] Sicsic | Initiation of a periodic array of cracks in the thermal shock problem: a gradient damage modeling[END_REF].

At t c a single macro-crack appears at the center of the domain. The post-critical development after t c of the system subject to an irreversibility constraint depends significantly on the domain width L and is hard to describe quantitatively in a rigorous and general manner. Assuming that the macro-cracks appearing at and after t c form a periodic pattern and are of equal length, their optimal spacing d opt (t) quickly decreases, which gives rise to rapid bisections of the periodic cells by the formation of new macro-cracks.

The latter always nucleate brutally and exhibit a very small cohesive zone. These bisections occur due to a loss of stability of the periodic cells similar to the one observed at t c and stop once a stable cell periodicity is reached. As soon as this is the case, a minimal periodicity is reached. At this point the growth of the micro-cracked zones is stopped. Since, however, the macro-cracks keep propagating, the importance of the micro-cracks will finally become negligible.

For large times t our results are in accordance with the scale law found in [START_REF] Bahr | Scaling behavior of thermal shock crack patterns and tunneling cracks driven by cooling or drying[END_REF] assuming a Griffith model.

Fig. 2

 2 Fig. 2 Cohesive surface energy density κ (top) and cohesive force κ (bottom) depending on the crack opening u for the family of models parameterized by ω. Dugdale's model corresponds to ω = 0.

Fig. 5 Fig. 4

 54 Fig.5Stress field σ11 in finite element simulations with η = 0.9 for d = {0.01, 0.1, 1}, all stresses below 0.7 plotted in the same color in order to improve resolution for higher stresses, y 2 = ¯ c marked as horizontal continued line.

Fig. 6

 6 Fig.6The body with a microcrack-layer near the upper side (left); the local structuration of the cracks at the micro-scale (right).

Fig. 7

 7 Fig.7Effective behavior of the microcracked material when the irreversibility condition is introduced.

Fig. 9

 9 Fig. 9 Black curve: critical time corresponding to the loss of stability of the microcrack layer versus the load parameter η.Red curve: the time t b at which appears another branch of solution than the homogeneous damage layer versus the load parameter η in the case of a gradient damage model[START_REF] Sicsic | Initiation of a periodic array of cracks in the thermal shock problem: a gradient damage modeling[END_REF]

Fig. 10

 10 Fig. 10 Rayleigh-like coefficient versus normalized domain width

  Fig. 11 Optimal macro-crack spacing as a function of time.

  Fig. 13 Normalized stress field σ 11 /σc for η = 0.4 on boundary x 1 = d + x1 of periodic half-cell at the time tc(d) of macro-crack initiation at x 1 = x1 (black dashed line) and short time later at t = 1.02 × tc(d) (gray line).

Fig. 12

 12 Fig. 12 Normalized stress field σ 11 /σc for η = 0.4 in periodic half-cell with d/(2 ktc(d)) = 40 (top) and d/(2 ktc(d)) = 2 (bottom) at the moment tc(d) directly after the formation of a macro-crack pattern with periodicity d in a only micro-cracked body (left) and short time later at t = 1.02 × tc(d) (right). At tc(d) the critical stress is reached exactly only at the tip of the macro-crack on the left boundary of the domain although for d/(2 ktc(d)) = 40 there is a broad layer in the upper right part, where it is almost reached. The thickness of micro-cracked layer before formation of macro-cracks at tc(d) is depicted by horizontal black line and the center of full periodic cell is denoted by x1 .

Fig. 14

 14 Fig.14Optimal periodicity in the reversible case in the physical space.

Table 1

 1 Main nomenclature

	2 Problem statement
	2.1 Geometry, thermoelastic behavior and thermal
	shock

Acknowledgements C. Cyron and J.-J. Marigo were partially supported by the French Agence Nationale de la Recherche (ANR) under Grant T-Shock (OTP J11R087).

A Appendix: Numerical implementation for macro-crack

In this appendix we briefly describe how the results depicted in Fig. 13 are obtained. The general setting is similar to the one discussed in Section 3.5 with the following two important differences:

1. a non-zero second derivative κ (0) < 0 of the cohesive surface energy is taken into account; 2. the constitutive behavior is described by ψ r instead of ψ.

The first difference is a natural consequence of carrying out simulations for a finite time in the physical rather than the normalized space. The second one allows us to account for a preexisting layer of micro-cracks in numerical computations addressing the possibility of a formation of macro-cracks and in fact even for the formation of new micro-cracks during the formation of the macro-cracks. Thus we admit under the constraint of irreversibility an arbitrary mixture of micro-and macro-cracks when looking for the optimal spacing of macrocracks. In COM SOL the second feature can be implemented by means of an elasto-plastic material behavior for the stress component σ 11 with yield stress σc. To speed up and stabilize the computations we made use of such a constitutive law only in parts of the domain such as the micro-cracked zone 0 ≤ x 2 ≤ t, because far from the micro-cracked layer the elasto-plastic behavior does not appear.

The difficulty in simulating the formation of a macrocrack in the finite element simulations is that the preexisting micro-crack layer is actually in equilibrium, albeit not a stable one, and the seeked configuration with macro-crack is only a second equilibrium configuration which has to be found. To this end, we define in the computation a fixed length Sv over which we allow the macro-crack to form and apply on this crack-boundary a cohesive load σ num c = τ num σc with a continuation parameter τnum ∈ {0; 1} for path continuation, which is continuously increased from zero to one. This way COM SOL will provide among the two equilibrium configurations -the one with and without a macro-crack -the one with larger opening of the macro-crack, as the path continuation starts from a configuration where the potential crack boundary of length Sv is not supported by any boundary stress.

A disadvantage of this method compared to the regularized boundary condition ( 24), which cannot directly be applied here, is that the crack length Sv is no longer calculated automatically, but that it is an initially unknown input parameter for the simulation. The proper Sv can then be found either by comparing the total mechanical energy obtained for different Sv or by searching the case where σ 11 = σc at the crack tip, because both criteria are equivalent according to [START_REF] Ferdjani | Insensitivity to small defects of the rupture of materials governed by the dugdale model[END_REF]. Practically, the second criterion turned out to be easier to apply. The temporal development between subsequent points in time which is depicted in Fig. 13 is simply simulated by imposing an increasing penetration depth of the thermal shock on the domain.