Aurélien Doitrand

Dominique Leguillon

Eric Martin

Aurelien Doitrand
email: aurelien.doitrand@gmail.com

Computation of generalized stress intensity factors of 3D singularities

Keywords: 3D elasticity, Generalized stress intensity factors, Coupled fracture criterion

come

Introduction

Since the pioneering work of Irwin [START_REF] Irwin | Analysis of stresses and strains near the end of a crack traversing a plate[END_REF] in 1957, the role of the stress intensity factor (SIF) of the crack tip singularity has been enhanced. Griffith's energy release rate G (MPa m) [START_REF] Griffith | The phenomenon of rupture and flow in solids[END_REF] was directly expressed in terms of the SIF

G = 1 -ν 2 E K 2 I (1)
Where E (MPa) is the Young modulus and ν the Poisson ratio of the material, and K I (MPa m 1/2 ) the SIF of mode I crack tip singularity. The above relationship (1) holds true within the 2D plane strain framework for a pre-existing crack growing straight under an opening load. It is based on the Williams' expansion [START_REF] Williams | On the stress distribution at the base of a stationary crack[END_REF] of an elastic solution U in the vicinity of the crack tip

U (x 1 , x 2 ) = C + K I √ ru I (ϕ) + ... (2)
The Cartesian (x 1 , x 2 ) and polar (r, ϕ) coordinates emanating from the singular point are mixed without risk of confusion. The leading term C is an irrelevant constant related to the rigid body translation and √ ru I (ϕ) is the opening mode I.

Eqn. [START_REF] Apel | Edge singularities and structure of the 3-D Williams expansion[END_REF] leads to an analogous relationship between the critical values of G and K I , respectively the fracture energy G c (MPa m) and the fracture toughness K Ic (MPa m 1/2 ) considered as material parameters

G c = 1 -ν 2 E K 2 Ic ( 3
)
More recently, the interest was focused on the generalized stress intensity factor (GSIF), also called notch stress intensity factor, of more general singularities at re-entrant corners [START_REF] Carpinteri | Stress singularity and generalized fracture toughness at the vertex of re-entrant corners[END_REF] or interfaces breaking a free surface in composites for instance. GSIF computation is crucial for failure initiation prediction since it can, for instance, help determining the crack initiation location and loading level in the presence of several singular points. The Williams' expansion is now written in the vicinity of the singular point

U (x 1 , x 2 ) = C + k 1 r λ 1 u 1 (ϕ) + ... ( 4
)
The GSIF is k 1 (MPa m 1-λ 1 ), λ 1 is the singularity exponent (for a V-notch in a homogeneous material 0.5 ≤ λ 1 ≤ 1) and u 1 (ϕ) is the associated mode. The series in (4) is made of successive terms with increasing exponents.

In 1997, Dunn et al. [START_REF] Dunn | Fracture initiation at sharp notches: correlation using critical stress intensities[END_REF] carried out experiments on V-notched specimens of PMMA with various opening angles and depths. They showed that the relevant parameter to predict crack initiation at the root of the notch is the GSIF k 1 of the corner singularity, new crack onset occurs if the GSIF reaches a critical value k c (MPa m 1-λ 1 ).

Shortly after, Leguillon [START_REF] Leguillon | Strength or toughness ? A criterion for crack onset at a notch[END_REF], using matched asymptotic expansions, established a relationship, based on a twofold condition in energy and stress, between this critical value k c on the one hand and the fracture energy G c and the tensile strength σ c (MPa) on the other hand

k c = G c ξ 1-λ 1 σ 2λ 1 -1 c (5)
Here ξ (MPa -1 ) is a scaling coefficient that depends on the geometry, the material properties and the boundary conditions. In case of a re-entrant corner in an isotropic media, it depends only on the notch opening, it can also depend of the local elastic properties in case of a bi-material for instance. It can be computed by FE independently of σ c and G c , which thus allows determining the critical GSIF for any material parameter couples (σ c ,G c ) using Eqn. [START_REF] Doitrand | 3D application of the coupled criterion to crack initiation prediction in epoxy/aluminum bimaterial specimens under four point bending[END_REF]. This relationship is at the origin of the coupled criterion (CC) to predict crack initiation at any singular points which has proven its effectiveness in many situations [START_REF] Weissgraeber | A review of Finite Fracture Mechanics: crack initiation at singular and non-singular stress-raisers[END_REF].

Various methods have been proposed to compute the 2D GSIFs: a least square fitting [START_REF] Dunn | Fracture initiation at sharp notches: correlation using critical stress intensities[END_REF] or an extraction from the strain energy density [START_REF] Lazzarin | Rapid calculations of notch stress intensity factors based on averaged strain energy density from coarse meshes: Theoretical bases and applications[END_REF] for instance. But the most accurate method, which allows in addition to get rid of the mixed modes, is based on a path independent integral [START_REF] Labossiere | Calculation of stress intensities at sharp notches in anisotropic media[END_REF][START_REF] Labossiere | Stress intensities at interface corners in anisotropic bimaterials[END_REF][START_REF] Leguillon | Computation of singular solutions in elliptic problems and elasticity[END_REF]. To the authors' knowledge, 3D GSIF computation has not been implemented yet. It differs from the 2D case because of extra geometric complexity, and differences from a theoretical point of view for the eigenvalue problem that must be solved to determine the singularity exponents and the associated 3D primal and dual modes.

All this takes place within the framework of plane elasticity. Except for the pure experimental results from Labossiere and Dunn in 2001 [START_REF] Labossiere | Fracture initiation at three-dimensional bimaterial interface corners[END_REF], it was only in 2014 that a first attempt was made to extend these results and more especially the asymptotic approach of the CC to the 3D case [START_REF] Leguillon | An attempt to extend the 2d coupled criterion for crack nucleation in brittle materials to the 3d case[END_REF]. Further 3D applications of the criterion were performed recently [START_REF] Doitrand | 3D application of the coupled criterion to crack initiation prediction in epoxy/aluminum bimaterial specimens under four point bending[END_REF][START_REF] Doitrand | Comparison between 2D and 3D applications of the coupled criterion to crack initiation prediction in scarf adhesive joints[END_REF][START_REF] Doitrand | Experimental characterization and numerical modeling of crack initiation in rhombus hole PMMA specimens under compression[END_REF][START_REF] García | A numerical study of transverse cracking in cross-ply laminates by 3D finite fracture mechanics[END_REF]; they are based on full FE computations due in part to a lack of tools to calculate the 3D GSIFs. Indeed, 3D GSIF computation is essential for 3D CC applications in its asymptotic form to establish the link between the imposed loading to a structure and the critical initiation GSIF. Few 3D asymptotic approaches have also been developed in the case of V-notches [START_REF] Mittelman | Asymptotic analysis of the potential energy difference because of a crack at a V-notch edge in a 3D domain[END_REF][START_REF] Yosibash | A 3-D failure initiation criterion from a sharp V-notch edge in elastic brittle structures[END_REF] where the GSIFs have been calculated using a method very specifically dedicated to V-notches [START_REF] Yosibash | Edge stress intensity functions in polyhedral domains and their extraction by a quasidual function method[END_REF].

The aim of this work is to remedy to this deficiency and to develop a method to compute the 3D GSIFs of any singularity using a path independent integral. Some applications to the 3D CC in its asymptotic form will be briefly presented. However, in some cases the smallness assumptions of the asymptotic approach may not be fulfilled, which requires the use of full FE calculations.

It should be noted that all the results presented here have more or less straightforward generalizations to mixed and complex modes.

In the following, when possible the index 1 is omitted for simplicity and terms refer to the leading singularity.

3D GSIF computation

Obviously the main difference with the 2D case is based on the complexity, especially the geometric complexity, of the 3D case, surfaces become volumes and lines become surfaces. But there are also differences of a more theoretical nature as seen below. The eigenvalue problem to be solved to compute the 3D singular exponents and associated modes differs from the 2D case. Essential in the computation of the generalized stress intensity factors, the rule to associate 3D primal and dual singularities differs as well. The transition from 2D to 3D is not so straightforward and moreover the 3D case has never been addressed in this way.

Singularity exponent, primal and dual fields

The 3D singular exponent as well as its associated eigenmode can be computed in a similar way to 2D as detailed in [START_REF] Leguillon | Computation of singular solutions in elliptic problems and elasticity[END_REF][START_REF] Leguillon | Computation of 3D-singularities in elasticity[END_REF]. It is briefly recalled here. The local geometry and material properties around the singular point, chosen as the origin, have to depend only on the spherical angles θ and ϕ. The generic example is formed by Fichera's corner as illustrated in Figure 1a.

We are only interested in local solutions of the form r λ u(θ, ϕ) (if exist) disregarding remote boundary conditions. We assume that there is locally no body forces and that the neighboring faces undergo vanishing either Dirichlet (displacements) or Neumann (tractions) conditions. Non homogeneous conditions can be added afterward [START_REF] Leguillon | Computation of singular solutions in elliptic problems and elasticity[END_REF]. The calculation of λ and u are based on the variational formulation (VF)

of the elasticity problem. The above searched function is inserted in the bilinear form of the VF together with test functions of the form φ(r)v(θ, ϕ)

Ω C : ∇U : ∇V dV = 0 ∀V with U (r, θ, φ) = r λ u(θ, ϕ) and V (r, θ, φ) = φ(r)v(θ, ϕ) (6)
Where C is the stiffness tensor. All these functions are assumed to have the required smoothness in order to get a well-posed problem, in particular, φ(r) has a compact support in the vicinity of the singular point. Under this assumption, integrating by part, we obtain a non-linear 2D (in θ and ϕ)

eigenvalue problem of the form Find λ and u such that

-λ(λ + 1)a(u, v) -(λ + 1)b(u, v) + λc(u, v) + d(u, v) = 0 ∀v (7)
where a, b, c and d are bilinear forms, a(u, v) and d(u, v) are symmetric while b(u, v) = c(v, u). Their detailed form, which can be found in [START_REF] Leguillon | Computation of 3D-singularities in elasticity[END_REF], are recalled in the Appendix. The above VF (7) can be discretized using finite elements, leading to allows solving [START_REF] Doitrand | Experimental characterization and numerical modeling of crack initiation in rhombus hole PMMA specimens under compression[END_REF] for λ, taking precautions to avoid overflows. Then X is obtained by the inverse iterations method.

Find (λ, X) ∈ R N +1 such that [-λ(λ + 1)A -(λ + 1)B + λC + D]X = 0 (8)
Remark 1: If λ is replaced by -λ -1 in (8) (or similarly in [START_REF] Doitrand | Numerical modeling of the nucleation of facets ahead of a primary crack under mode I+III loading[END_REF]) and taking into account the symmetry properties of the matrices A and D as well as the relationship B = C T , leads to conclude that if λ is a solution to [START_REF] Doitrand | Numerical modeling of the nucleation of facets ahead of a primary crack under mode I+III loading[END_REF][START_REF] Doitrand | Experimental characterization and numerical modeling of crack initiation in rhombus hole PMMA specimens under compression[END_REF] then -λ -1 is a solution too. The condition λ > -1/2 is necessary to have a solution with a finite energy, nevertheless, -λ -1 is a mathematical solution to [START_REF] Doitrand | Numerical modeling of the nucleation of facets ahead of a primary crack under mode I+III loading[END_REF][START_REF] Doitrand | Experimental characterization and numerical modeling of crack initiation in rhombus hole PMMA specimens under compression[END_REF] with its own eigenvector u -(θ, ϕ). The functions u(θ, ϕ) and u -(θ, ϕ) are called respectively primal and dual modes. This property will play an essential role in the calculation of the 3D GSIFs.

Remark 2: There is a major difference with the 2D case where the dual exponent to λ ≥ 0 is -λ [START_REF] Leguillon | Computation of singular solutions in elliptic problems and elasticity[END_REF].

2D and 3D dual modes are also baptized the complementary displacement fields by some authors [START_REF] Labossiere | Calculation of stress intensities at sharp notches in anisotropic media[END_REF][START_REF] Labossiere | Stress intensities at interface corners in anisotropic bimaterials[END_REF].

The exponents and the discretized values of the primal and dual modes (displacements and stresses) are tabulated to be used in the calculations described in the sequel.

The contour integral for GSIF calculation

In the vicinity of a singular point, the solution to an elastic problem U can be developed in terms of the primal modes already described in the above subsection. It is the so called 3D Williams' expansion

[27, 1] U (x 1 , x 2 , x 3 ) = C + kr λ u(θ, ϕ) + ... (10)
Where again Cartesian (x 1 , x 2 , x 3 ) and spherical (r, θ, ϕ) coordinates are mixed without risk of confusion.

Here, C is an irrelevant constant related to the rigid body translation and k (MPa m 1-λ 1 ) is the GSIF of the singular term. The goal of this section is to propose a method to extract this GSIF from U or more likely from a finite element (FE) approximation of it: U FE . To this aim, let us consider the following contour integral for two elastic solutions U and V fulfilling locally the homogeneous equilibrium equations and boundary conditions

Ψ(U , V ) = S [σ(U ).n.V -σ(V ).n.U ]dS (11)
Where S is any surface encompassing the singular point and finishing on the neighboring traction free (or vanishing displacements) faces, σ(U ) is the stress field associated to a displacement field U through Hooke's law, n is the normal to S pointing towards the origin (the singular point). This integral is similar to the integral employed for computing the 2D GSIFs [START_REF] Labossiere | Calculation of stress intensities at sharp notches in anisotropic media[END_REF][START_REF] Labossiere | Stress intensities at interface corners in anisotropic bimaterials[END_REF][START_REF] Leguillon | Computation of singular solutions in elliptic problems and elasticity[END_REF] except for the integration contour which is a line in 2D and a surface in 3D. The path independence of the integral ( 11) can be proved by considering two contours S 1 and S 2 encompassing the singular point and ending on the neighboring traction free faces. Defining D as the volume delimited by both contours S 1 and S 2 and the neighboring traction free faces Σ 1 and Σ 2 , the symmetry of the stiffness tensor C ensures:

D C : ∇U : ∇V dV = D C : ∇V : ∇U dV (12)
and therefore :

D σ(U ) : ∇V dV - D σ(V ) : ∇U dV = 0 (13)
Integrating by parts and using Green theorem it comes:

D ∇.σ(U ).V -∇.σ(V ).U dV + ∂D σ(U ).n.V -σ(V ).n.U dS = 0 (14)
The first integral in Eqn. ( 14) is null since U and V both satisfy equilibrium equations (∇. is the divergence operator). The contour ∂D is composed of the two contours S 1 and S 2 as well as the traction free faces Σ 1 and Σ 2 . Noticing that the second integral in Eqn. ( 14) is null over Σ 1 and Σ 2 (traction free faces) and that the S 1 and S 2 normal directions are opposite, it comes:

S 1 σ(U ).n.V -σ(V ).n.U dS - S 2 σ(U ).n.V -σ(V ).n.U dS = 0 ( 15
)
which proves the integral path independence. As a particular case, let us consider two singular exponents α and β solutions to [START_REF] Doitrand | Experimental characterization and numerical modeling of crack initiation in rhombus hole PMMA specimens under compression[END_REF], their eigenmodes u α (θ, ϕ) and u β (θ, ϕ) and the associated stress fields through Hooke's law r α-1 s α (θ, ϕ) and r β-1 s β (θ, ϕ). They locally fulfill the equilibrium equations, then considering the integral Ψ [START_REF] Griffith | The phenomenon of rupture and flow in solids[END_REF] where S is a sphere of radius R, it comes

Ψ(r α u α , r β u β ) = R α+β+1 θ φ [s α (θ, ϕ).n.u β (θ, ϕ) -s β (θ, ϕ).n.u α (θ, ϕ)] sin θdθdϕ (16)
As a consequence of the path independence, the above integral vanishes except for β = -α -1, i.e. if β is the dual exponent to α. Thus, noting

U = i k i r α i u α i (θ, ϕ), it comes: Ψ(U , r -α-1 u -α-1 ) = i k i Ψ(r α i u α i , r -α-1 u -α-1 ) = k α Ψ(r α u α , r -α-1 u -α-1 ) (17)
Therefore we obtain:

k = Ψ(U FE , r -λ-1 u -) Ψ(r λ u, r -λ-1 u -) (18)
Note that the denominator in ( 18) can be calculated once and for all independently of the contours as a consequence of ( 16). It can be emphasized that as long as the different modes (e.g. opening, symmetric and antisymmetric shear modes) correspond to different singularity exponents, the 3D GSIFs associated to each mode can be computed separately using Eqn. [START_REF] Lazzarin | Rapid calculations of notch stress intensity factors based on averaged strain energy density from coarse meshes: Theoretical bases and applications[END_REF].

In case of multiple modes (i.e. modes that share the same singularity exponent, for instance in the case of a crack), say m modes, ( 18) is replaced by a system m × m, since the dual modes have the same multiplicity m than the primal ones. Moreover, it must be emphasized that the approach described in this subsection is also appropriate for complex exponents and modes.

The coupled criterion

The coupled criterion basically states that crack initiation requires the simultaneous fulfillment of two conditions:

-The tensile stress must overcome the material strength all along the presupposed crack path prior to initiation.

-The potential energy (W p ) change due to crack initiation must be larger than the energy G c S required to create a crack of surface S.

Both conditions can be written in the following form:

   G inc = -δWp S G c
σ σ c throughout the presupposed crack surface [START_REF] Leguillon | Computation of singular solutions in elliptic problems and elasticity[END_REF] where G inc is called the incremental energy release rate. In practice, solving Eqn. ( 19) reverts to determine the minimum loading level for which both conditions are satisfied. This can be achieved by computing both the stress variation along the crack path and the potential energy variation due to crack initiation for several crack lengths.

The asymptotic form of the 3D coupled criterion

As already mentioned, the 3D CC is based as in 2D on the Williams' expansion [START_REF] García | A numerical study of transverse cracking in cross-ply laminates by 3D finite fracture mechanics[END_REF] of an elastic solution in the vicinity of the singular point. The challenge is now to estimate the shape of the crack at initiation. Three assumptions are made to overcome this difficulty:

1) In the first stage of initiation, the crack is plane, as observed experimentally for instance in [START_REF] Pham | On the growth of cracks under mixed-mode I + III loading[END_REF].

In any case, this can be considered as an acceptable simplification, the initial crack being assumed to be small.

2) The crack is embedded in the plane undergoing the maximum tensile stress (the crack angle can also be obtained as the one minimizing the initiation loading obtained using the CC).

3) The shape of the crack is defined by the isolines of the tensile stress acting on the plane (see Fig. 2), in the asymptotic case they are homothetic curves).

This allows defining the crack size through a single parameter s, which can be chosen for instance as the crack extension in a given direction [START_REF] Doitrand | Numerical modeling of the nucleation of facets ahead of a primary crack under mode I+III loading[END_REF][START_REF] Doitrand | Experimental characterization and numerical modeling of crack initiation in rhombus hole PMMA specimens under compression[END_REF]. The latter option introduces a coefficient b to derive a relation between s and S S = bs 2 (20)

Note that both definitions are strictly equivalent provided there exists a bijection between the crack surface and the crack extension along the given direction. Such a condition was verified in previous works [START_REF] Doitrand | 3D application of the coupled criterion to crack initiation prediction in epoxy/aluminum bimaterial specimens under four point bending[END_REF][START_REF] Doitrand | Comparison between 2D and 3D applications of the coupled criterion to crack initiation prediction in scarf adhesive joints[END_REF][START_REF] Doitrand | Numerical modeling of the nucleation of facets ahead of a primary crack under mode I+III loading[END_REF][START_REF] Doitrand | Experimental characterization and numerical modeling of crack initiation in rhombus hole PMMA specimens under compression[END_REF] where the possible initiation crack shapes were determined based on the stress isocontours.

In the vicinity of a singular point, the change in potential energy due to crack initiation over a surface S = bs 2 is:

-δW p = k 2 s 2λ ζ + ... (21)
and the stress along the presupposed crack path prior to crack initiation writes:

σ = ks λ-1 + ... (22)
The CC states that crack initiation may occur if the tensile stress is larger than the material strength and if the change in potential energy is larger the crack surface creation energy G c S:

   -δWp S = k 2 c s 2λ-1 ζ G c σ = k c s λ-1 σ c (23)
Then combining equations ( 21) the critical value k c of the GSIF k when crack onset occurs is

k c = G c ζ 1-λ σ 2λ-1 c (24)
In an analogous way to ξ in (5), ζ (MPa -1 ) is a scaling coefficient that depends on the local geometry and possibly on the local elastic properties in case of a bi-material for instance. The main difference between the 2D (Eqn. ( 5)) and 3D (Eqn. ( 24)) cases being the dependency on the coefficient b involved The isolines of the tensile stress acting on the presupposed crack surface (here the interface) and the parameter s in the case [START_REF] Labossiere | Fracture initiation at three-dimensional bimaterial interface corners[END_REF] (see Sec. 5). The lines were obtained numerically [START_REF] Doitrand | 3D application of the coupled criterion to crack initiation prediction in epoxy/aluminum bimaterial specimens under four point bending[END_REF], for small values of s compatible with the asymptotic hypothesis, they form homothetic lines.

in the relationship between the surface S of the new crack and the parameter s (see assumption 3 above and Eqn. ( 20)).

Implementation in a FE code

Once a FE calculation has been performed on a structure which geometry includes a 3D singularity, the computation of the GSIF only requires a post-processing of this calculation. The definition of a domain for computing the contour integral is first required. It can be done either analytically in simple cases where the contour is for instance, a quarter, half or a full sphere or cube. A simple discretization of such contours in (θ,ϕ) space can be easily obtained with a grid. It is thus straightforward to compute the weight of the integration points to estimate k based on [START_REF] Lazzarin | Rapid calculations of notch stress intensity factors based on averaged strain energy density from coarse meshes: Theoretical bases and applications[END_REF]. For more complex contour shapes, it is also possible to obtain the integration contour by performing Boolean operations between a sphere (or any other geometrical objects encompassing the singular point) and the geometry around the singular point. First, the complementary geometry of the singular point is generated. In the case of the Fichera's corner depicted in Fig. 1, the complementary geometry simply is a cube. Then, a difference between the sphere and this complementary geometry is applied so as to obtain a contour that strictly ends on the faces attached to the singular point. An example of a sphere-based contour obtained by this procedure is displayed in Fig. 1b-c. The contour can then be discretized by generating a surface mesh of the resulting contour geometry. Then, the weight of the integration points is obtained by computing the surface of the elements of this mesh. Once the contour geometry is defined, the second step consists in extracting the quantities of interest on the contour (displacements, stresses). FE displacements and stresses are obtained through post-processing the FE solution by interpolating the calculated quantities at the contour integration points using element shape functions and from the tables providing the primal and dual modes. The typical GSIF computation time is some seconds. In the following, all the FE calculations have been performed using Abaqus TM and the complete procedure is sketched in Figure 3.

Path independence verification on examples of 3D GSIF computation

Glued cube

This example is intended to be a preliminary study of the detachment of electronic chips under the effect of thermal stresses. However, we did not address the case of a thermal loading. It could be done and trigger a little more complexity but located upstream of the implementation at the asymptotic expansions level (see for instance [START_REF] Henninger | Adhesive fracture of an epoxy joint under thermal and mechanical loadings[END_REF] for thermo-elastic problems and [START_REF] Leguillon | Computation of singular solutions in elliptic problems and elasticity[END_REF] for a more general point of view on non-homogeneous problems). It plays no role on the implementation itself. Thus, to focus on our main concern we prefer considering a simplified mechanical loading acting on a 1 × 1 × 1 mm 3 cube.

A 0.1 mm displacement is prescribed on the upper face (z = 1 mm) whereas the bottom face (z = 0 mm) is clamped and the other faces are traction free. We consider a model material which parameters are E = 1 GPa, ν = 0.3. The 3D singular exponent is λ = 0.624 whereas the 2D singularity along the edges, which is also a 3D singularity [START_REF] Apel | Edge singularities and structure of the 3-D Williams expansion[END_REF], is weaker λ 2D = 0.711. Thus, failure will likely initiate from the corner. The contour encompassing the singular point for GSIF calculation can be chosen as one sphere eighth, corresponding to (θ,ϕ) ∈ [0, π 2 ]×[0, π 2 ] in spherical coordinates. In this case, the only point that can be checked is the independence of the integral Ψ (cf. Eqn. [START_REF] Griffith | The phenomenon of rupture and flow in solids[END_REF]) with respect to the integration path and as a consequence the independence of the calculated values of the GSIF k. Results are listed in Tab. 1 for ten contours with radius lying between 0.05 mm and 0.5 mm displayed in Fig. 4b using 30 × 30 integration points on the contour, resulting in differences on the computed GSIF lower than 0.1%. The same magnitude of differences is obtained if 10 × 10 integration points are used. Differences rises up to 0.5% and 1.3% respectively if only 5 × 5 or 3 × 3 integration points are used. The same order of magnitude is obtained for the differences between several contours using either linear or quadratic elements. It was also checked for this example and the following ones that differences lower than 1% were obtained with different mesh element densities.

Fichera's corner

The geometry of the Fichera's corner specimen studied herein is obtained by removing a cube of side 1 mm from a cube of side 5 mm (see Fig. 1a), it is made of the same model material than above (E = 1 GPa, ν = 0.3). Indeed, the mechanical interest is lower since the 3D singularity at the corner is λ = 0.776 whereas the 2D singularity along the edges is stronger λ 2D = 0.545, so that failure will likely initiate from the edges rather than at the corner. Nevertheless, the proposed approach allows computing the GSIF at the corner. The bottom face (z = 0 mm) is clamped whereas a 0.5 mm vertical displacement is prescribed on the upper face (z = 5 mm). The contour encompassing the singular point for GSIF calculation can be chosen as seven eighths of a sphere, corresponding to (θ,ϕ) 2: GSIF values computed for several contour radius between 0.05 mm and 0.5 mm in the case of the Fichera's corner.

∈ [0, π 2 ]×[ π 4 , 7π 4 ] ∪ [ π 2 ,π]×[0,

Application to V-notch specimen experiments

GSIF computation based on experimental data

This section addresses 3D GSIF computation in the case of PMMA V-notch specimens undergoing 3-point bending load. Dunn et al. [START_REF] Dunn | Fracture initiation at sharp notches: correlation using critical stress intensities[END_REF] performed such experiments on specimens with dimensions L=76.2 is the same whatever the V-notch depth). The failure load was recorded for all V-notch angle/depth configurations and the GSIFs were computed by least-squares fitting procedure of the asymptotic displacement fields along the V-notch flanks from 2D FE calculations of three-point bending of V-notch specimens using the loading level recorded experimentally. The material properties given in [START_REF] Dunn | Fracture initiation at sharp notches: correlation using critical stress intensities[END_REF] are E = 2.3 GPa, ν = 0.36. Three dimensional simulations of three-point bending of V-notch specimens with angles and depths described previously have been performed in order to compute the GSIFs at the middle of the V-notch edge. Meshes with 8-node linear elements are used with a refined zone along the V-notch singular edge. The contour encompassing the singular point (located at the center of the V-notch edge) for GSIF calculation is chosen as three quarters of a sphere (see Fig. 5), corresponding to (θ,ϕ)

∈ [0,π]×[- π 2 + β 2 , 3π 2 -β 2 ]
in spherical coordinates. First, the GSIF for each depth/angle V-notch configuration has been computed prescribing the failure load measured experimentally. Fig. 6a displays the computed GSIFs (plain and dashed line) and those obtained in [START_REF] Dunn | Fracture initiation at sharp notches: correlation using critical stress intensities[END_REF] by least square fitting (error bars) as a function of the V-notch depth to specimen height ratio a/h for several angles. It can be noted that for a given V-notch configuration (fixed value of depth a and angle β), three values are obtained for the numerical predictions corresponding to the minimum, maximum (dashed lines) and mean (plain line) failure forces measured experimentally. These values are displayed as error bars in Figs. 6 and7. It seems that an almost constant GSIF is obtained as a function of the V-notch depth for β =60 and 90 deg. angles, and a slight decrease with increasing V-notch depth is obtained for β =120 deg. angle. This could be explained in part by a larger scattering in the measures when the opening angle increases [START_REF] Leguillon | Crack initiation at a v-notch under complex loadings -Statistical scattering[END_REF]. The predicted GSIFs are in good agreement with data obtained from least-square fitting method based on experiments as shown in Fig. 6b which displays the mean GSIF obtained using the contour integral or least-square fitting method as a function of the V-notch angle. Fig. 7 shows the failure stress as a function of the Vnotch depth to specimen height ratio for several V-notch angles measured experimentally and obtained numerically for the GSIF level measured experimentally. Similarly to Fig. 6, three values are obtained for each numerical prediction corresponding to the minimum, maximum (dashed lines) and mean (plain line) GSIF levels measured experimentally. The failure stress decreases with increasing V-notch depth to specimen height ratio or increasing V-notch angle.

GSIF computation based on coupled criterion prediction

In this section, a full numerical procedure is employed so as to predict both the initiation failure stress and the corresponding GSIFs for all the studied V-notch configurations. First, the stress at failure is obtained numerically through a 3D application of the CC. The general steps that are followed so as to implement the CC solution are briefly described hereafter. For a more detailed description of the procedure, the reader can refer to [START_REF] Doitrand | Experimental characterization and numerical modeling of crack initiation in rhombus hole PMMA specimens under compression[END_REF] which includes a 3D application of the CC in a similar case as treated herein. For each configuration, a first FE calculation with no crack is performed so as to compute the stress field in the specimen middle plane at the V-notch. The stress isocontours in this plane are used so as to define the possible crack initiation shapes. This method is convenient because it allows to strictly fulfill the stress condition of the CC on the whole crack surface and also to define a three dimensional crack by a single parameter, for instance its surface or its extension along a given direction (see Subsection 2.4). Then, for each isoline (15 lines are enough in practice), a mesh of the V-notch specimen including a crack that follows this isoline is generated, which allows computing the potential energy released by such a crack and hence the energy condition of the CC. The CC is then solved by determining the minimum loading for which both the stress and the energy conditions are fulfilled, which is the initiation loading. Solving the CC requires the knowledge of the material fracture parameters. These parameters were taken within the bounds measured experimentally in [START_REF] Dunn | Fracture initiation at sharp notches: correlation using critical stress intensities[END_REF] so that a σ c =140 MPa strength and G c =0.53 MPa mm fracture energy were used in the calculations. Fig. 8 shows the predicted failure stress as a function of the V-notch depth to specimen height ratio for several V-notch angles. The same trends as observed experimentally are obtained using the CC as well as a satisfying agreement with the failure stress measured experimentally. Fig. 9 shows the GSIF computed numerically for the failure stress levels predicted by the CC. For these failure stress levels, a almost constant GSIF level is obtained whatever the V-notch depth to specimen height ratio. The condition that must be verified to ensure a constant GSIF whatever the notch depth is the crack smallness assumption of the asymptotic approach. It has been checked that the initiation crack extensions in the specimen middle plane predicted by the CC are lower than 0.02mm, therefore it is likely that this assumption is verified in this case. A good agreement with GSIFs obtained in [START_REF] Dunn | Fracture initiation at sharp notches: correlation using critical stress intensities[END_REF] based on experiments by least-square fitting is obtained for β =60 and 90 deg. angles whereas for β =120 deg. the predictions overestimate the GSIF level obtained using least-square fitting approach, which decreases with increasing V-notch depth. It is not straightforward to explain the decrease in the GSIF for a β = 120 deg. angle whereas it is not for the other angles. It has to be noted that in [START_REF] Dunn | Fracture initiation at sharp notches: correlation using critical stress intensities[END_REF], the GSIF is not obtained by a direct experimental measurement but requires a numerical least-squares procedure to be indirectly determined, which could explain the observed decrease.

Application to mixed mode loading

The example chosen to illustrate the computation of GSIFs in the case of mixed mode is 3-point bending of PMMA specimens containing a slanted V-notch, the V-notch orientation angle being denoted γ. It can be noted that for γ = 0 deg., the configuration studied in Section 4 is retrieved. The specimen length, width and height are the same as those described in Section 4: L=76.2 mm × l=17.8 mm × the specimen center, a mixed mode loading is generated along the slanted V-notch. This leads to crack nucleation of a network of facets that grow and coalesce, then resulting in a progressive rotation of the crack so as to reach pure mode I loading from a macroscopic point of view, as observed experimentally [START_REF] Buchholz | Fracture analyses and experimental results of crack growth under general mixed mode loading conditions[END_REF][START_REF] Lazarus | Comparison of predictions by mode II or mode III criteria on crack front twisting in three or four point bending experiments[END_REF][START_REF] Citarella | Comparison of crack growth simulation by DBEM and FEM for SEN specimens undergoing torsion or bending loading[END_REF], see also [START_REF] Doitrand | Numerical modeling of the nucleation of facets ahead of a primary crack under mode I+III loading[END_REF] for a theoretical analysis based on the CC of this mechanism in case of a blunted notch.

The singularity exponents corresponding to mode I, mode II and mode III are computed following the method presented in Section 2:

-mode I : λ I =0.545 -mode III : λ III =0.668 -mode II : λ II =0.908 Remark 3: Note that, it is through a slight misuse of writing that the modes are still referred to as I, II and III, a notation a priori reserved for a crack. Mode I means symmetric opening mode, mode II antisymmetric mode in the y direction and mode III antisymmetric mode in the z direction (see Fig. 10). To avoid any ambiguity, GSIFs are noted as lower-case k with arabic indices. It is emphasized that the singularity exponents are the same all along the V-notch except at the two end points where the V-notch edge reaches the specimen free surfaces.

An example of contour used to compute the GSIF is depicted in Fig. orientation γ. For γ =0 deg., mode II (k 2 ) and III (k 3 ) GSIFs are zero and the pure mode I loading configuration presented in Section 4 is retrieved. With increasing V-notch orientation, k 2 remains zero and k 1 decreases while k 3 increases. The ability to compute the mode mixity for a given configuration and loading allows making the link with the initiation GSIF computed using the matched asymptotic approach of the CC [START_REF] Doitrand | Numerical modeling of the nucleation of facets ahead of a primary crack under mode I+III loading[END_REF][START_REF] Leguillon | An attempt to extend the 2d coupled criterion for crack nucleation in brittle materials to the 3d case[END_REF].

The proposed method allows computing the GSIFs all along the V-notch edge. Fig. 11b shows the mode I, II and III GSIFs variations along the V-notch edge as a function of the relative position z/h (the V-notch edge center is located at z/h=0.5 whereas the edge end is located at z/h=1). It can be seen that k 1 and k 3 are almost constant near the center of the V-notch edge, and slightly decrease or increase when approaching the specimen traction free faces, whereas k 2 increases from 0 in the center to its maximum value at the specimen edge. It means that crack initiation is more likely to begin from the edge center. However, experimental evidences of crack initiation starting from the edge center would require the use of a rapid camera and the ability to capture the very moment of fracture initiation since at initiation the crack nucleates quasi-instantaneously over a finite surface [START_REF] Doitrand | Experimental characterization and numerical modeling of crack initiation in rhombus hole PMMA specimens under compression[END_REF] .

Application to a bimaterial interface corner

Examples of GSIF calculation in previous sections involve a 3D singular point in homogeneous materials. This section addresses GSIF computation for a singular point located at a bimaterial interface.

The configuration under investigation is aluminum/epoxy bimaterial specimens with a square cross section (width w, cf. Fig. 12). Labossiere and Dunn [START_REF] Labossiere | Fracture initiation at three-dimensional bimaterial interface corners[END_REF] performed four point bending tests on several specimens with different widths to study crack initiation, a 45 deg. rotation of the specimens being adopted so that crack initiation occurs at the bimaterial interface corner under an opening mode. The singular exponent corresponding to this mode, computed using the method described in Section 2, is λ =0.640. The material properties taken from [START_REF] Labossiere | Fracture initiation at three-dimensional bimaterial interface corners[END_REF] are E a =70 GPa, ν a =0.33 for aluminum and E e =2.98

GPa and ν e =0.38 for epoxy. The meshes consisting of linear hexahedra are refined in the vicinity of the bimaterial interface corner. The contour encompassing the singular point for GSIF calculation can be chosen as one quarter of a sphere (Fig. 12), corresponding to (θ,ϕ)

∈ [0,π]×[ π 4 , 3π 4
] in spherical coordinates. For each specimen width, ten different contours have been used to compute the GSIF, an example of values computed on the different contours in the case of w=12.5 mm specimen width is shown in Tab. 3. It can be noted that in the case of a singular point located at a bimaterial interface, the GSIF variations obtained using several contours are larger than in the case of homogeneous materials.

The GSIFs computed as a function of the specimen width for the load at failure measured by Labossiere and Dunn [START_REF] Labossiere | Fracture initiation at three-dimensional bimaterial interface corners[END_REF] are presented in Tab. 4. The GSIF obtained in [START_REF] Labossiere | Fracture initiation at three-dimensional bimaterial interface corners[END_REF] by a least-squares fitting procedure are also recalled. For the largest specimen width, similar GSIF as computed in [START_REF] Labossiere | Fracture initiation at three-dimensional bimaterial interface corners[END_REF] is obtained with the present approach. However, it is found that the GSIF decreases with decreasing specimen width whereas an almost constant value was obtained in [START_REF] Labossiere | Fracture initiation at three-dimensional bimaterial interface corners[END_REF]. Since whatever the specimen width the 3D singularity is the same, according to [START_REF] Labossiere | Fracture initiation at three-dimensional bimaterial interface corners[END_REF] and to the 3D CC, one could expect that the same GSIF at crack initiation is obtained.

However, this is only the case if the assumptions of the asymptotic approach are verified as explained in Section 4. The decrease in GSIF with decreasing specimen size may thus be explained by the fact that as w decreases, the asymptotic assumption of crack initiation smallness with respect to other dimensions does not hold true. The initiation crack extension as a function of the specimen width was already computed in a previous work using the CC [START_REF] Doitrand | 3D application of the coupled criterion to crack initiation prediction in epoxy/aluminum bimaterial specimens under four point bending[END_REF]. The crack extension at initiation to specimen size along the corner bisector ratio s c /( √ 2w) is provided for the different specimen widths in Tab. 4.

For small specimen widths, it is clear that the assumption of a small initiation crack with respect to other dimensions is not respected. Obviously, the larger the specimen the truer the assumption and the better the agreement with experimental results. Contour radius (mm) 0.3 0.37 0.43 0.5 0.57 0.63 0.70 0.77 0.83 0.9 GSIF (MPa.mm 0.376 ) 7.43 7.30 7.49 7.47 7.60 7.70 7.85 7.70 7.82 7.94

Conclusion

Similarly to the 2D case, the GSIF of a 3D singularity can be calculated using a path independent integral, using the so-called dual mode to the primal singular mode. The transition from 2D to 3D presents some changes from a theoretical point of view related to the singularity exponent and primal/dual mode calculation, but, above all, it is the implementation that constitutes the main obstacle.

Calculations must be carried out on various contours and the independence of the results to the contours is an element that makes it possible to rule out implementation errors, without obviously also guaranteeing the correctness.

The knowledge of an accurate value of a GSIF allows using the CC in its asymptotic form [START_REF] Mittelman | Asymptotic analysis of the potential energy difference because of a crack at a V-notch edge in a 3D domain[END_REF] to predict crack nucleation at singular points, provided the asymptotic assumption of smallness (i.e. s (see Fig. 2) small compared to any dimensions of the structure) holds true.

A point that has not been considered but that can be addressed without any additional difficulty is the study of the singularities at the ends of a V-notch at a traction free surface (see Sect. 5). This is a crucial information for studying the initiation of a crack either from an inner point of the notch root or from the ends as in Section 5.

Among the immediate perspectives, we plan to study the singularities around a Vickers indenter, which presents both singular vertices and edges. This can be achieved quite quickly by neglecting plastic deformations or compaction and then apply the CC to predict the formation of cracks and compare these predictions to to the rich amount of data, in particular experimental micrographs for crack shape analysis. Future work will also cover the experimental GSIF determination using Digital Image Correlation..

Appendix

We recall here the expression of the bilinear forms a, b, c and d defined in Eqn. [START_REF] Doitrand | Numerical modeling of the nucleation of facets ahead of a primary crack under mode I+III loading[END_REF]

Figure 1 :

 1 Figure 1: (a-b) Geometry of Fichera's corner and (b-c) integration contour obtained using Boolean operation between a sphere and the Fichera's corner.

  where A, B, C and D are N × N matrices, A and D are symmetric and B = C T (transposed). The equation Det[-λ(λ + 1)A -(λ + 1)B + λC + D] = 0 (9)

Figure 2 :

 2 Figure2: The isolines of the tensile stress acting on the presupposed crack surface (here the interface) and the parameter s in the case[START_REF] Labossiere | Fracture initiation at three-dimensional bimaterial interface corners[END_REF] (see Sec. 5). The lines were obtained numerically[START_REF] Doitrand | 3D application of the coupled criterion to crack initiation prediction in epoxy/aluminum bimaterial specimens under four point bending[END_REF], for small values of s compatible with the asymptotic hypothesis, they form homothetic lines.

Figure 3 :

 3 Figure 3: Steps for the computation of integral (18).

Figure 4 :

 4 Figure 4: (a) FE mesh of the glued cube under investigation refined in the vicinity of the singular point. (b) Contours used for the calculation of the GSIF.

Figure 5 :

 5 Figure 5: V-notch specimen and example of contour used for GSIF computation.

Figure 6 :

 6 Figure 6: (a) Critical GSIF k c as a function of normalized V-notch depth a/h obtained in [9] by least-square fitting (error bars) and using the contour integral for the failure stress levels σ f measured experimentally (plain and dashed lines respectively corresponding to the predictions obtained using the mean or the minimum/maximum failure stresses measured experimentally) for several V-notch opening angles β. (b) Critical GSIF k c as a function of V-notch angle β obtained by least square fitting and using the contour integral.

Figure 7 :

 7 Figure7: Failure stresses σ f as a function of normalized V-notch depth a/h obtained experimentally (error bars) and numerically for the GSIF levels obtained in[START_REF] Dunn | Fracture initiation at sharp notches: correlation using critical stress intensities[END_REF] by least-square fitting (plain and dashed lines respectively corresponding to the predictions obtained using the mean or the minimum/maximum GSIF measured experimentally) for several V-notch opening angles β.

Figure 8 :

 8 Figure 8: Failure stresses σ f obtained experimentally (error bars) and predicted numerically by the coupled criterion (plain line) for several V-notch opening angles β using G c =0.53 MPa.mm and σ c =140 MPa.

Figure 9 :

 9 Figure 9: Critical GSIF k 1c obtained from experiments by least-square fitting (error bars) and numerically (plain line) for the failure stress level σ f predicted by the coupled criterion for several V-notch opening angles β using G c =0.53MPa.mm and σ c =140MPa.

Figure 10 :

 10 Figure 10: Specimen with a slanted V-notch and example of contour used for mode I and mode III GSIF computation.

  h=12. mm. A a =3.56 mm deep V-notch with β =90 deg. angle is studied. The material properties are E=2.3 GPa and ν=0.36. Linear hexahedra are used, the meshes being refined in the V-notch singular edge center. Fig. 10 depicts the specimen geometry and an example of contour used for GSIF calculations. Contrary to the example presented in Section 4 corresponding to a pure mode I loading at

10

 10

Figure 11 :

 11 Figure 11: Mode I, II and III GSIFs (a) computed at the V-notch in the middle of the specimen as a function of the V-notch orientation γ and (b) computed along the V-notch edge (z/h=0.5 corresponds to the specimen middle and z/h=1 corresponds to the specimen edge) for a γ = 30 deg. V-notch orientation.

Figure 12 :

 12 Figure 12: Epoxy-Aluminum bimaterial specimen under investigation and example of contour used for the GSIF computation.

Table 1 :

 1 GSIF values computed for several contour radius between 0.05 mm and 0.5 mm in the case of the glued cube with 10 × 10 integration points.

	Contour radius (mm) 0.05	0.1	0.15	0.2	0.25	0.3	0.35	0.4	0.45	0.5
	GSIF (MPa.mm 0.376 ) 4.769 4.771 4.771 4.772 4.772 4.772 4.772 4.773 4.773 4.773

Table

  2π] in spherical coordinates. Differences on GSIF values are lower than 1.% for 10 contours with radius lying between 0.05 mm and 0.5 mm as can be observed in Tab. 2.

	Contour radius (mm) 0.05	0.1	0.15	0.2	0.25	0.3	0.35	0.4	0.45	0.5
	GSIF ( MPa.mm 0.224 ) 2.374 2.388 2.392 2.395 2.395 2.395 2.395 2.395 2.395 2.395

Table 3 :

 3 GSIF values computed for several contour radius between 0.3 mm and 0.9 mm in the case of the bimaterial specimen with 4.0 mm width. MPa.mm 0.360 ) 4.98 + 0.38 6.41 + 0.60 6.78 + 0.25 7.62 + 0.32 8.32 + 0.21 GSIF from[START_REF] Labossiere | Fracture initiation at three-dimensional bimaterial interface corners[END_REF] (MPa.mm 0.360 ) 8.83 + 0.66 8.57 + 0.66 8.39 + 0.72 8.45 + 0.64 8.42 + 0.59 s c /(

	Specimen width w (mm)	4.0	6.2	8.9	12.5	17.8
	Force at failure (N)	554.8	212.1	84.4	32.3	10.2
	GSIF (√	2w)	0.24	0.15	0.08	0.04	0.02

Table 4 :

 4 GSIF as a function of the specimen width obtained numerically and experimentally from[START_REF] Labossiere | Fracture initiation at three-dimensional bimaterial interface corners[END_REF]

  . The C ijkh are the components of the stiffness tensor C. ijkh (A j B h u i ∂v k ∂θ + A j C h u i

					a(u, v) =	C ijkh A j A h u i v k sin(θ)dθdϕ		(25)
							S
			b(u, v) =	S	C ijkh (B j A h	∂u i ∂θ	v k + C j A h	∂u i ∂ϕ	v k ) sin(θ)dθdϕ	(26)
			c(u, v) =								∂v k ∂ϕ	) sin(θ)dθdϕ	(27)
	d(u, v) =	S	C ijkh (B j B h	∂u i ∂θ	∂v k ∂θ	+ B j C h	∂u i ∂θ	∂v k ∂ϕ	+ C j B h	∂u i ∂ϕ	∂v k ∂θ	+ C j C h	∂u i ∂ϕ	∂v k ∂ϕ	) sin(θ)dθdϕ	(28)

S

C


Luiz Felipe

Faria Ricardo

Dominique Leguillon

Guillaume Parry

Aurelien Doitrand
email: aurelien.doitrand@insa-lyon.fr

Modelling the thermal shock induced cracking in ceramics

Keywords: Quenching, Thermal shock, Coupled criterion, Cohesive zone

Prediction of surface cracking in ceramics due to quenching is performed numerically using either the coupled criterion or a cohesive zone model. Under such a thermal shock, a network of short cracks with minimal spacing between them initiate and propagate until some of them stop while the others continue propagating. The numerical implementation consists of a periodic array of cracks modeled by a representative volume element. It allows crack initiation, simultaneous propagation and period doubling to be predicted. The investigation of the crack period doubling allows a precise determination of the optimal crack spacing, which decreases with an increasing thermal shock amplitude. The predicted crack spacing results are in agreement with experimental measurements.

Introduction

Industrial ceramics are known for their low coefficient of thermal expansion (CTE) which gives them high thermal stability, hence their use in thermal engines and in the aeronautics and aerospace industry. However, it is observed that surface cracking, an undesired phenomenon, occurs when they are sharply cooled down [START_REF] Bahr | Scaling behavior of thermal shock crack patterns and tunneling cracks driven by cooling or drying[END_REF]. Usually, in a metal, quenching triggers a brutal shrinkage of the surface which deforms plastically. When cooling down goes on, the core of the material in turn shrinks, putting the surface under compression and leading to a strengthening effect. Unfortunately, being brittle, ceramics cannot deform plastically, so if the surface shrinkage is too large, cracks will appear in order to accommodate the deformation. To have experimental evidence of this phenomenon, quenching experiments have been conducted by various authors on ceramic bars (Al 2 O 3 ) using a specific setup [START_REF] Jiang | A study of the mechanism of formation and numerical simulations of crack patterns in ceramics subjected to thermal shock[END_REF][START_REF] Shao | Crack Patterns in Ceramic Plates after Quenching[END_REF][START_REF] Shao | Effect of Crack Pattern on the Residual Strength of Ceramics After Quenching[END_REF][START_REF] Wu | Size effect of thermal shock crack patterns in ceramics andnumerical predictions[END_REF] (see Fig. 1).

Figure 1: One side of a specimen after a thermal shock about 400 K (reprinted from [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 thermal shock caused by a sudden cooling plays an important role in the fracture pattern. Indeed, in case of a slow cooling or drying process of a thin film on a rigid substrate which prevents effective shrinkage, a first group of cracks is formed with a given (rather large) spacing and new cracks are continually added as the cooling or drying process goes on, ultimately leading to saturation [START_REF] Jenkins | Optimal spacing and penetration of cracks in a shrinking slab[END_REF][START_REF] Leguillon | A simple model of thermal crack pattern formation using the coupled criterion[END_REF][START_REF] Leguillon | Multi-cracking in brittle thin layers and coatings using a FFM model[END_REF][START_REF] Leguillon | Prediction of multi-cracking in sub-micron films using the coupled criterion[END_REF]. On the contrary, when cooling down is sudden, it seems that a set of cracks with minimal spacing appears. These cracks are constrained to grow with the minimal spacing until an instability takes place leading to several successive period doubling events [START_REF] Bahr | Scaling behavior of thermal shock crack patterns and tunneling cracks driven by cooling or drying[END_REF][START_REF] Jenkins | Optimal spacing and penetration of cracks in a shrinking slab[END_REF][START_REF] Jenkins | Determination of crack spacing and penetration due to shrinkage of a solidifying layer[END_REF].

The numerical simulation of such a phenomenon has proved difficult to carry out. In parallel with the simple model proposed by Bahr [START_REF] Bahr | Scaling behavior of thermal shock crack patterns and tunneling cracks driven by cooling or drying[END_REF] other approaches have emerged. Jenkins [START_REF] Jenkins | Optimal spacing and penetration of cracks in a shrinking slab[END_REF][START_REF] Jenkins | Determination of crack spacing and penetration due to shrinkage of a solidifying layer[END_REF] and Jiang et al. [START_REF] Jiang | A study of the mechanism of formation and numerical simulations of crack patterns in ceramics subjected to thermal shock[END_REF] used a method based on finite elements (FE) computations and an energy minimization principle to predict crack spacing. To this aim Li et al. [START_REF] Li | Direct numerical simulations on crack formation in ceramic materials under thermal shock by using a non-local fracture model[END_REF] proposed a convincing, but not so easy to implement, method relying on a non-local damage model. More recently a phase field approach was developed [START_REF] Bourdin | Morphogenesis and Propagation of Complex Cracks Induced by Thermal Shocks[END_REF][START_REF] Sicsic | Initiation of a periodic array of cracks in the thermal shock problem: A gradient damage modeling[END_REF]. In this latter approach, a thin layer of homogeneously damaged material first appears on the surface before an instability takes place leading to the formation of a first periodic network of close cracks. The subsequent study of the growth of these cracks, which is allowed by the phase field method, describes the instabilities leading to period doubling.

The phenomenon of multi-cracking resulting of surface shrinking is in its first stage more akin to a crack initiation problem than to a crack growth. The traditional tools of brittle fracture mechanics are therefore ineffective. Specific models must be used, two of them are proposed hereafter. The coupled criterion (CC) combines an energy and a stress condition [START_REF] Leguillon | Strength or toughness? A criterion for crack onset at a notch[END_REF][START_REF] Weissgraeber | A review of Finite Fracture Mechanics: Crack initiation at singular and non-singular stress-raisers[END_REF] and the cohesive zone model (CZM) relies on a traction-separation relationship across fracture surfaces [START_REF] Park | Cohesive Zone Models: A critical review of traction-separation relationships across fracture surfaces[END_REF]. They allow determining, for a given temperature drop ∆T and a given crack spacing d:

1) crack initiation time and length, 2) crack period doubling point time and length, 3) crack length variation as a function of time.

Note that, the further crack growth is either simulated using the classical Griffith criterion or by using the cohesive zone model beyond the crack initiation stage.

It can be pointed out that the proposed model allows varying the crack spacing d which is not known a priori. A method aiming at determining the optimal crack spacing is presented in the following. In both approaches, the repeated crack pattern is described using a periodicity assumption. The width of the spacial period corresponds to the crack spacing.

The thermal conduction problem

To model the surface cracking phenomenon occurring due to surface shrinking in the very early times following the thermal shock and solve it numerically, it is important to have a sufficiently refined thermal conduction model, especially with regard to boundary conditions. It is briefly recalled below.

We consider the 2D problem of an infinitely long along Ox specimen with width 2w so that the heat equation is a 1D equation in Oy direction. The specimen is symmetric with respect to the Ox axis, T , t and y respectively denotes the temperature, time and space variable.

The heat equation writes:

ρ ∂T ∂t = a ∂ 2 T ∂y 2 (1)
where

a = k ρcp (m 2 s -1 ) is the diffusivity coefficient, k (W m -1 kg -1
) the thermal conductivity, ρ (kgm -3 ) the density and c p (J kg -1 K -1 ) the specificic heat capacity. We define the dimensionless reduced time t, space variable y and temperature T as:

t = at w 2 ; y = y w ; T = T -T 0 T ext -T 0 (2)
where T 0 is the initial temperature in the specimen and T ext is the outer temperature. In the present case of a thermal shock leading to shrinking and cracking T 0 >> T ext . As a consequence, Eqn. (1) rewrites:

∂T ∂t = ∂ 2 T ∂y 2 (3)
The Fourier condition on the bottom face of the specimen (y = -w) is :

k ∂T (-w, t) ∂y = -h(T (-w, t) -T ext ) (4)
where h (W m 2 K -1 ) is the convective heat transfer coefficient. Introducing the dimensionless Biot number B = hw k , Eqn. (4) becomes:

∂T (-1, t) ∂y = -BT (-1, t) (5)
In addition the symmetry condition along the (Ox) axis is

∂T ∂y = 0 i.e. ∂T ∂y = 0 (6)
The solution to this system is [START_REF] Carslaw | Conduction of heat in solids[END_REF]:

T (y, t) = i>0 A i exp(-c 2 i t)cos(c i y) -1 (7)
with

A i = 2sin(c i ) c i + sin(c i )cos(c i ) ( 8
)
where the c i 's are solutions to the transcendent equation:

c i sin(c i ) -Bcos(c i ) = 0 (9)
In the sequel, for simplicity, the material coefficients are assumed to be independent of the temperature although this is only an approximation [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 temperature variation in time and space is modeled by implementation of the solution [START_REF] Carslaw | Conduction of heat in solids[END_REF] to the heat equation in a user subroutine with 30 terms. Results are illustrated in Figure 2. Small oscillations can be seen in Figure 2b for t = 0.0005 due to the limited number of terms used in [START_REF] Carslaw | Conduction of heat in solids[END_REF]. It must be pointed out that the above equations and solution remain unchanged if cracks develop parallel to the 0y axis, since there is no thermal exchange in the Ox direction.

The coupled criterion to predict crack nucleation

The CC combines an energy and a stress condition. This approach is particularly well suited to simulate the initiation of new cracks at stress concentration points in brittle materials, such as ceramics [START_REF] Leguillon | Application of the coupled stress-energy criterion to predict the fracture behaviour of layered ceramics designed with internal compressive stresses[END_REF][START_REF] Leguillon | What is the tensile strength of a ceramic to be used in numerical models for predicting crack initiation?[END_REF][START_REF] Martin | Understanding the tensile strength of ceramics in the presence of small critical aws[END_REF].

The criterion

The CC extends in a way the Griffith criterion and allows to predict crack onset at stress concentration points [START_REF] Weissgraeber | A review of Finite Fracture Mechanics: Crack initiation at singular and non-singular stress-raisers[END_REF]. It meets two conditions: an energy one, close to the Griffith condition although written in an incremental form instead of a differential one, and a stress condition. It does not question Griffith's criterion, it coincides with it at the tip of a crack and with the maximum tension that a material can undergo (strength) in the absence of stress concentration. The major property that emerges from this criterion is that the initiation of a new crack is sudden with a jump over a given length (often small) before continuing to grow or stop.

Its implementation can be done in two ways: (i) the asymptotic approach, not used here for reasons that will be explained below, assumes that the initial crack jump is small compared to all the dimensions of the studied structure; (ii) the full FE approach which has fewer restrictive conditions but requires a higher numerical effort.

Let us consider an uncracked structure and solve the elastic or thermo-elastic problem using FE. First, the potential energy W P (0) can be computed as well as the tension σ(s) acting on the presupposed crack path (thus prior to crack initiation) emanating from the stress concentration point under consideration (s is the abscissa along the presupposed crack path). Then, if double nodes have been inserted along this presupposed path, each pair can be released step by step to simulate a crack of increasing length l. At each step the corresponding potential energy W P (l) is computed and the incremental energy condition can be written:

-∆W P = W P (0) -W P (l) ≥ G c l ⇒ G inc ≥ G c with G inc = W P (0) -W P (l) l ( 10
)
Where G c is the fracture energy of the material identified with the material toughness through Irwin's relationship. Here G inc is the so-called incremental energy release rate. It differs from the usual definition of the energy release rate G which is of a differential form

G = - ∂W P ∂l (11)
On the other hand, the stress condition states that the tensile stress prior to crack formation must exceed the tensile strength of the material σ c all along the initial crack jump:

σ(s) ≥ σ c for 0 ≤ s ≤ l (12)
These two inequalities allow to compute the crack jump l c at initiation and the critical prescribed load triggering the new crack formation. This full FE approach is adopted for two reasons. First, as observed in the experiments, cracks that appear are close to each other and the spacing is of the same order of magnitude than cracks length, violating the asymptotic assumption that the crack length is far smaller than any other dimension of the structure. Second, in the early times of the quenching the temperature gradient is very large which can lead to the definition of a second small parameter that may interact with the crack length.

The CC representative volume element

As already mentioned, a large number of cracks initially appear close to one another and it would be extremely cumbersome to model them all. To avoid this problem, calculations are carried out on an elementary volume representative of the entire structure by periodically repeating the pattern. This means that we study the central part of the specimen, far from the ends. The width d of this representative volume element (RVE) is a parameter of the model that will be used to represent the crack spacing (Figure 3, not at the right scale). This approach has already been used successfully in [START_REF] Leguillon | Multi-cracking in brittle thin layers and coatings using a FFM model[END_REF][START_REF] Leguillon | Prediction of multi-cracking in sub-micron films using the coupled criterion[END_REF].

Because of symmetries, only 1/4 of the RVE (Figure 3) is meshed, its height is 1/2 of the specimen width w. Along the presupposed crack path the mesh size increases upward from 1.7 µm to 8 µm at 0.2 mm from the bottom (it is checked afterward that this is big enough to cover the actual crack length). It then continues to increase to reach 0.2 mm on the upper part.

Due to the symmetry of the specimen with respect to Ox there is no bending, the actual displacement field in the RVE U is taken as the superposition of an elementary thermally activated field U th (baptized shrinking in Figure 4) and an elementary mechanical traction U el acting on the left and right faces (stretching in Figure 4): Where ∆T 0 = T ext -T 0 is the prescribed temperature drop. The unknown weight z, function of t, is linked to the average horizontal strain εxx

U (x, y, t) = ∆T 0 U th (x, y, t) + z U el (x, y) (13)
z = εxx d 2 ( 14
)
It is obtained writing that the resultant force R(U ) acting on the lateral faces of the specimen vanishes

∆T 0 R(U th ) + z R(U el ) = 0 ( 15
)
The two elementary solutions fulfil the balance equation

-∇.σ th = 0; -∇.σ el = 0 ( 16
)
Where the stress tensors σ th and σ el are defined by the constitutive equations

σ th = C : ∇U th -αI (17)
and

σ el = C : ∇U el (18)
Here C is the Hooke tensor depending on Young's modulus E and Poisson's ratio ν of the material (isotropy is assumed). In [START_REF] Leguillon | Multi-cracking in brittle thin layers and coatings using a FFM model[END_REF], α is the CTE and I the second order identity tensor (isotropy is assumed as well). Eqn. ( 17) is the thermoelastic constitutive law for a temperature change equal to unity. The bottom face Γ b of the RVE is traction free

σ th xy = σ th yy = 0 σ el xy = σ el yy = 0 (19)
Horizontal displacements are prescribed on the left face Γ l of the RVE

U th x = 0; U el x = -1 (20)
For symmetry reasons calculations are carried out only on 1/4 of the RVE (see Figure 3) and the two elementary solutions fulfil the symmetry conditions along Ox U th y = 0; σ th xy = 0; U el y = 0; σ el xy = 0 [START_REF] Li | Direct numerical simulations on crack formation in ceramic materials under thermal shock by using a non-local fracture model[END_REF] and along Oy U th x = 0;

σ th xy = 0; U el x = 0; σ el xy = 0 ( 22
)
Once U th and U el have been numerically computed using FE for various crack length l, the change in potential energy between the uncracked and cracked states can be calculated [START_REF] Leguillon | Prediction of multi-cracking in sub-micron films using the coupled criterion[END_REF] -

∆W P = 1 2 ∆T 2 0 ∆(KX th X th ) + ∆T 0 z∆(BX el ) - 1 2 z 2 ∆(KX el X el ) ( 23
)
Where K is the stiffness matrix that results from the discretization of the two problems (the same in both cases), X th and X el are the nodal vectors solution to the discretized systems and B the vector of nodal forces involved in the thermoelastic problem. Note that in (23) X th and z both depend on t and that this calculation must be carried out at each time step.

Once the cracks have initiated, the further growth can be studied using Griffith's criterion, the derivative being approximated by a difference in potential energy at each step of unbuttoning:

- ∂W P ∂l ≥ G c (24)
It is thus possible to compare from an energetic point of view, on the one hand the simultaneous crack growth using the RVE of Figure 3 and on the other hand the growth of one crack out of two by growing only the left crack in the RVE of Figure 5.

CC Results

For convenience, although it is a cooling process, the values of ∆T 0 = T ext -T 0 and ∆T = T -T 0 are noted as positive in the following.

The simulations have been carried out on the Alumina already described in [START_REF] Li | Direct numerical simulations on crack formation in ceramic materials under thermal shock by using a non-local fracture model[END_REF]. The material parameters are: Young's modulus E = 370 GPa, Poisson's ratio ν = 0.3, tensile strength σ c = 180 MPa, fracture energy G c = 24.3 J m -2 (toughness

K Ic = 3.14 MPa m 1/2 ), density ρ = 3980 kg m -3 , thermal conductivity k = 31 W m -1 K -1 , specific heat capacity c p = 880 J kg -1 K -1 , convective heat transfer coefficient h = 5 10 4 W m -2 K -1 (Biot's number B ≈ 8).
The way the CC is applied is illustrated in Figures 6, 7 and 8 for t = 1.6 ms and a given crack spacing 2 mm. It is written in a dimensionless way so that the two conditions 10 and 12 can be plotted on the same graph. Obviously, in Figure 6 the stress condition holds true for l ≤ 0.06 mm but the energy condition nowhere. No crack initiation is possible after t = 1.6 ms for a temperature drop ∆T 0 = 500 K. Under the same conditions, but with ∆T 0 = 550 K (Figure 7), again no possible crack initiation, the two conditions holds true but on disjointed segments. Initiation can occur for ∆T 0 = 600 K, with l c = 0.075 mm. This length is the initiation length but in general not an arrest length as seen further. Finally, from Figure 9 for ∆T 0 = 650 K, it is clear that it still exists a solution at an earlier time than t = 1.6 ms.

The elapsed time before crack nucleation depends on the temperature drop as shown in Figure 10. Moreover, for temperature drops smaller than ∆T 0 = ∆T min = 120 K no crack nucleation is predicted. The actual crack spacing is determined using the method developed in [START_REF] Doitrand | Numerical modeling of the nucleation of facets ahead of a primary crack under mode I+III loading[END_REF][START_REF] Leguillon | Multi-cracking in brittle thin layers and coatings using a FFM model[END_REF][START_REF] Leguillon | Prediction of multi-cracking in sub-micron films using the coupled criterion[END_REF]. At a given initiation time t, various decreasing spacings are tested. In each case, the temperature drop likely to cause crack onset is determined based on the CC. It is observed that this temperature drop keeps constant for large spacings and then starts to increase as the crack spacing becomes smaller than a threshold (Figure 11). The actual crack spacing is taken as the smallest spacing prior to the increase in the temperature drop, it corresponds to the maximum dissipation that can be reached at a given temperature drop. According to Figure 11 it is located between 0.8 and 1.4 mm, seemingly independent of the intensity of the thermal shock. It is the order of magnitude found in [START_REF] Jiang | A study of the mechanism of formation and numerical simulations of crack patterns in ceramics subjected to thermal shock[END_REF], however, it is obviously difficult to capture accurately the end of the plateau. This may not be very surprising taking into account the large scattering seen in the experiments, spacings of 0.4 mm to 1.6 mm are observed for instance for ∆T 0 = 400 K (see Figure 21 further on). It can be noted that the difficulty to accurately compute the optimal crack spacing is not due to a CC limitation but rather to the plateau beginning determination method itself since the same difficulty was met using the above method with CZM, therefore a different strategy will be proposed in Section 4.3.

To study the crack growth after initiation, the energy release rate G (Eqn. ( 11)) is superimposed in Figure 12. Then, according to [START_REF] Doitrand | Numerical modeling of the nucleation of facets ahead of a primary crack under mode I+III loading[END_REF][START_REF] Leguillon | Multi-cracking in brittle thin layers and coatings using a FFM model[END_REF][START_REF] Leguillon | Prediction of multi-cracking in sub-micron films using the coupled criterion[END_REF], it is clear in figure 12 that following the initiation, the crack grows almost instantaneously in an unstable manner (G > G c ). The "arrest length"has two bounds (quotation marks are used because it would be a real being entirely consumed in rupture (refer to [START_REF] Leguillon | Prediction of multi-cracking in sub-micron films using the coupled criterion[END_REF]). These bounds, together with the initiation length, are plotted in Figure 13 as functions of the temperature drop ∆T . It should be noted that the experimental measurement of the initiation and arrest lengths is not straightforward since it requires the ability to capture the very moment of crack initiation, for instance using a high-speed camera [START_REF] Shao | Observation of ceramic cracking during quenching[END_REF].

In the following, the lower bound will be used, this feature has been observed in many other cases except when the excess energy is extremely large [START_REF] Leguillon | Application of the coupled stress-energy criterion to predict the fracture behaviour of layered ceramics designed with internal compressive stresses[END_REF]. ∆T Once initiated and "stopped", the subsequent growth of the crack can be studied using Griffith's criterion. We propose to examine what happens immediately after this first phase, when time continues to elapse, by comparing the growth of the crack in the RVE depicted in Figure 3 (corresponding to the simultaneous growth of all cracks), to the growth of a single crack in the RVE depicted in Figure 5 (corresponding to the growth of one crack out of two, i.e. the doubling of the period). More precisely, twice the energy released by the growth of the crack in case 1 is compared to the energy released by the single crack in case 2 and the most favourable case from the point of view of the energy (i.e. the case leading to the maximum dissipation) is selected [START_REF] Nguyen | Energy methods in fracture mechanics: stability, bifurcation and second variations[END_REF]. Results are summarized in Table 1. Results of Table 1 are qualitatively in agreement with a comment in [START_REF] Jiang | A study of the mechanism of formation and numerical simulations of crack patterns in ceramics subjected to thermal shock[END_REF] stating that the short cracks become shorter as ∆T 0 increased.

The cohesive zone model

The model

A second approach allowing quenching-induced crack initiation and growth to be modeled is cohesive zone modeling (CZM) [START_REF] Barenblatt | The formation of equilibrium cracks during brittle fracture: general ideas and hypotheses[END_REF][START_REF] Dugdale | Yielding of steel sheets containing slits[END_REF]. Such an approach has already been employed to model cracking in ceramics materials for instance under mechanical [START_REF] Vargas | Analysis of a castable refractory using the wedge splitting test and cohesive zone model[END_REF] or thermal [START_REF] Buljak | Numerical modeling of thermally induced microcracking in porous ceramics: An approach using cohesive elements[END_REF] loading. It defines a relation between the cohesive traction T n and the opening displacement δ n of a process zone. As soon as the traction reaches a critical value (taken herein as the material strength), debonding occurs which triggers a decrease in the cohesive traction (damage) until the opening displacement reaches a critical value δ c for which the traction becomes null, which corresponds to the nucleation of a crack. The area of the traction-separation profile is the work for separation, i.e. G c for an elastic material with no dissipation effects. For a given traction-separation profile, a relation between G c , σ c and δ c can be straightforwardly derived so that only two of these parameters are required as input. The CZM approach thus requires the same ingredients as for the CC, i.e. the material strength σ c and fracture energy G c . In addition, the traction separation profile can also be varied. However, all the cohesive zone models are equivalent for studying stable rectilinear crack propagation [START_REF] Acary | Nonsmooth fracture dynamics using a cohesive zone approach[END_REF] and the traction-separation profile only has an influence on crack initiation [START_REF] Alfano | On the influence of the shape of the interface law on the application of cohesive-zone models[END_REF][START_REF] Doitrand | Comparison between cohesive zone and coupled criterion modeling of crack initiation in rhombus hole specimens under quasi-static compression[END_REF] or branching. In the present work, the influence of the traction-separation profile is not studied and a bilinear profile is adopted throughout the analysis (Figure 14).

The CZM representative volume element

As above and for analogous reasons, the FE model consists of an infinite specimen modeled by repeating a representative volume element (RVE) with periodic boundary conditions, which is thus representative of the specimen center far for the lateral edges. For practical reasons, the present RVE contains two parallel crack paths on which cohesive zone elements are inserted, which allows modeling initiation of a periodic array of cracks, the simultaneous growth of these cracks as well as the possible subsequent crack period doubling if one out of two cracks stops propagating (Fig. 15).

A material characteristic length defined by

L mat = EGc σ 2 c
allows estimating the process zone extent, which is actually a fraction of this length [START_REF] Turon | An engineering solution for mesh size effects in the simulation of delamination using cohesive zone models[END_REF], its value also depending on the prescribed boundary conditions. It was checked that a typical process zone size is around 0.12 mm, which corresponds to 0.44 L mat . The mesh is refined near the cracks so as to ensure that the process zone contains enough elements to well reproduce crack initiation and propagation. The mesh size along the crack path is set to 5 × 10 -3 mm.

Symmetry boundary condition U y (x, y 0 ) = 0 is applied on the top face of the RVE, moreover a displacement U x (x 0 , y 0 ) = 0 at the top left node (see Fig. 15)) is also prescribed to prevent rigid body motions. Periodic boundary conditions are applied to the RVE so as to model an infinite array of cracks. This means that, as in periodic homogenization, the strain field is periodic whereas the traction exerted on the lateral faces are antiperiodic, as a consequence the displacement field is the sum of a periodic function and a linear one

U x (x 0 + 2d, y) -U x (x 0 , y) = 2d εxx U y (x 0 + 2d, y) -U y (x 0 , y) = 2d εxy ( 25
)
for anypoint (x, y) belonging to the left boundary. Here εxx and εxy are the mean strain tensor elongation and shear components (see ( 14)), which are unknown a priori since a thermal loading is applied. However, according to Hooke's law σ xy = 2µε xy , then it comes εxy = 0 since σxy = 0 because (i) no efforts are applied on the lateral faces of the specimen, (ii) the thermal expansion coefficients are isotropic. In addition, according to the condition that prevents rigid body motion, U x (x 0 + 2d, y 0 ) = 2dε xx .

To summarize, the boundary conditions can be written (x varies from 0 to 2d and y from 0 to w/2)

U y (x, y 0 ) = 0 U y (x 0 + 2d, y) = U y (x 0 , y) U x (x 0 , y 0 ) = 0 U x (x 0 + 2d, y) = U x (x 0 , y) + 2d εxx with εxx = Ux(x0+2d,y0) 2d (26)
Due to the choice of a non-symmetrical RVE, writing boundary conditions is slightly more complicated than in Section 3.2 but the above conditions [START_REF] Shao | Effect of Crack Pattern on the Residual Strength of Ceramics After Quenching[END_REF] can be taken into account by most FE codes. In this case, it is essentially the number of cohesive zone elements and the total number of degrees of freedom that determine the cost of the computation.

Results with CZM

The model presented in Section 4.1 is now applied so as to determine the optimal crack spacing as a function of the temperature drop. Fig. 16 shows the crack length variation as a function of time for the two cracks in the RVE for a temperature drop ∆T 0 = 480 K and a crack spacing d = 1 mm. It can be seen that, first, both cracks initiate and propagate simultaneously until reaching a certain length which corresponds to the crack period doubling point, beyond which only one crack out of two propagates while the other one stops. This phenomenon can be explained by the fact that propagating only one crack from this point is a more energetically favourable configuration. This is not contradictory with results in Table 1 where for such a temperature drop, it is concluded that after initiation (∆T 0 = 500 K and l arrest 0.185 mm) there is simultaneous growth of the whole set of cracks, without auguring of further evolution.

We performed simulations with RVE's (Figure 15) containing either one crack (RVE width d) or two cracks (RVE width 2d) so as to compare models in which the crack period doubling is likely to occur (RVE with two cracks) with a case in which both cracks would continue propagating simultaneously (RVE with one crack). For crack lengths smaller than the one corresponding to the crack period doubling point, the energy of both configurations (i.e. twice the energy in the single crack case) are exactly the same, which is consistent with the simultaneous crack propagation observed prior to the doubling point in Fig. 16. For larger crack lengths, the energy of the configuration with only one out of two cracks propagating is smaller than if both cracks continue propagating (i.e. there is a larger dissipation). Therefore, the former is more energetically favorable than the latter and will be observed in practice. The crack period doubling point can be determined as a function of the crack spacing for a given temperature drop. Fig. 17 shows the crack length variation as a function of time for several crack spacings for a ∆T 0 = 380 K temperature drop. It can be observed that the smaller the crack spacing, the smaller the time and the length of the crack at the period doubling point.

The crack length at the doubling point as a function of the crack spacing, depicted in Fig. 18, exhibits a linear variation. Therefore, it is possible to determine the minimum spacing d opt for which crack period doubling occurs, which corresponds to the intersection between the line and the abscissa axis. It means that below this value, only one out of two cracks will initiate in the RVE, which finally corresponds to a configuration of a RVE with twice the crack spacing. Therefore, the minimum spacing for which crack period doubling occurs corresponds to the optimal spacing between cracks.

The same procedure has been repeated for several temperature drops. Fig. 19 shows the crack length at the crack period doubling point variation as a function of the crack spacing obtained for temperature drops between 280 and 580 K.

In all cases, a linear variation is obtained which allows determining the optimal crack spacing. The higher the temperature drop, the lower the optimal crack spacing. 20 shows the optimal crack spacing as a function of the temperature drop. A regression of the computed optimal spacings d opt as a function of the temperature drop ∆T 0 has been performed The chosen regression form has the advantage of involving a single parameter ∆T min which has a physical meaning since it corresponds to the minimum temperature drop below which no cracking occurs during the ceramic quenching. A least-squares fitting of this parameter leads to ∆T min = 146 K. It can be seen that the chosen function allows a satisfying representation of the optimal spacing variation as a function of the temperature drop and also give access to ∆T min , which is an essential data for manufacturing ceramic parts that may undergo temperature drops. Experimental results are taken from [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 shaded area corresponds to the bounds determined in Figure 11.

A comparison between the experimental results from [START_REF] Jiang | A study of the mechanism of formation and numerical simulations of crack patterns in ceramics subjected to thermal shock[END_REF] and the predictions of the FEM calculations using CZM is displayed in Figure 21. A very large scattering in the crack spacing measured experimentally can be observed. The error bars in the numerical approach correspond to d opt determined by [START_REF] Shao | Observation of ceramic cracking during quenching[END_REF] and twice this value. Indeed, we assume that the crack network will subdivide until it falls within this interval for which no further subdivision can occur since d opt is the smallest possible spacing. The shaded area corresponds to the two bounds determined using the CC (see Section 3.3 and Figure 11).

Conclusion

The thermal shock induced by a sudden cooling of ceramics results in the nucleation of a set of cracks more or less regularly distributed along the specimen edge. Both CC and CZM allow the determination of the initiation crack length and temperature drop corresponding to a given crack spacing by combining a stress and an energy conditions. Both approaches are complementary, they basically yield similar crack initiation loading and spacing, and their numerical efficiency depends on the searched purpose. A rough estimate of the initial crack spacing can also be obtained by detecting the beginning of a plateau in the temperature drop variation as a function of the crack spacing. However, although it is illusory trying to give a precise spacing so large is the experimental scattering, this method lacks accuracy.

At this stage, the biggest difference lies in the computation time for obtaining the results. The CZM requires solving nonlinear problems which calculation cost depends essentially on the number of nodes located in the cohesive zones. This cost is necessarily high to get a satisfying accuracy. On the contrary, the CC only requests linear resolutions, as many as there are nodes on the presupposed crack path, with a substantial saving in the calculation cost. However, the advantage of the CC stops here, it is only dedicated to the initiation phase, studying the subsequent growth of cracks and the period doubling are beyond its capabilities. This can be done using the Griffith criterion but the prodedure is difficult because it is necessary to consider numerous crack configurations. On the contrary, the same CZM approach can be continued for analyzing the further crack growth without introducing new complications. This is particularly interesting for investigating in depth the doubling of the period. This cohesive zone approach together with the period doubling made it possible to develop a strategy to determine the initial crack spacing with a higher precision, highlighting a trend in relation to the amplitude of the thermal shock, results are in agreement with the experimental measurements.

Figure 2 :

 2 Figure 2: Reduced temperature (a) drop T and (b) gradient as a function of the reduced distance to the edge at different reduced time t.

Figure 3 :

 3 Figure 3: Left: the RVE used to study the simultaneous onset and growth of a periodic set of cracks. Right: 1/4 of the RVE, the dashed lines are symmetry axes, the red solid line at the bottom left is the presupposed crack path.

Figure 4 :

 4 Figure 4: 1/4 RVE stretching due to a unit traction and shrinking due to a unit cooling.

Figure 5 :

 5 Figure 5: Left: the RVE used to study a possible doubling of the period. Right: 1/4 of the RVE, the dashed lines are symmetry axes, the red solid line at the bottom left is the growing crack and the one at the bottom right is an arrested one.

Figure 6 :

 6 Figure 6: The stress and energy conditions at t = 1.6 ms and ∆T 0 = 500 K for d = 2 mm. The stress condition is fulfilled for l ≤ 0.06 mm but the energy condition nowhere.

Figure 7 :

 7 Figure 7: The stress and energy conditions at t = 1.6 ms and ∆T 0 = 550 K for d = 2 mm. Ranges on which the stress criterion and the energy criterion are disjointed.

Figure 8 :

 8 Figure 8: The stress and energy conditions at t = 1.6 ms and ∆T 0 = 600 K for d = 2 mm. There is crack initiation with lc = 0.075 mm.

Figure 9 :

 9 Figure 9: The stress and energy conditions at t = 1.6 ms and ∆T 0 = 650 K for d = 2 mm. There exists a solution at an earlier time.

Figure 10 :Figure 11 :

 1011 Figure10: Actual temperature drop ∆T = T -T 0 (as already mentioned, although this is a cooling process, the ∆T values are noted positive) at different times t corresponding to crack nucleation ∆T 0 = 600 K, t = 1.6 ms, ∆T 0 = 500 K, t = 2.2 ms, ∆T 0 = 400 K, t = 3.3 ms, ∆T 0 = 300 K, t = 5.8 ms, ∆T 0 = 200 K, t = 14.8 ms (following the arrow).

Figure 12 :Figure 13 :

 1213 Figure 12: The crack initiation length and the two bounds for the arrest length derived from the CC.

Figure 14 :

 14 Figure 14: The bilinear cohesive zone model.

Figure 15 :

 15 Figure 15: Left: the RVE used with cohesive zones. Right: 1/2 of the RVE, the dashed line is a symmetry axis. Presupposed crack paths where cohesive zone elements are inserted are vertical lines embedding the red solid lines. The red solid lines are an example of cracks after doubling the period.

Figure 16 :

 16 Figure 16: Cracks extensions as a function of time for a RVE containing two cracks.

Figure 17 :Figure 18 : 1 ∆T0Figure 19 :Figure 20 :

 171811920 Figure 17: Crack extensions as a function of time for several crack spacings d.

Figure 21 :

 21 Figure 21: Crack spacing determined experimentally (red dashed line) and numerically (blue solid line) as a function of the temperature drop.Experimental results are taken from[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 shaded area corresponds to the bounds determined in Figure11.

Table 1 :

 1 Temperature drop, elapsed time at initiation, initiation length, "arrest length"after initiation, further growth. Values differ slightly from those used previously to take into account results presented in Figure11.

	growth

0 (K) t (ms) d (mm) l init (mm) l arrest (mm) Further