{'fraction_non_alphanumeric': 0.03049207089552239, 'fraction_numerical': 0.012826492537313433, 'mean_word_length': 4.878341329677862, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Morphogenesis and propagation of complex cracks induced by thermal shocks

We study the genesis and the selective propagation of complex crack networks induced by thermal shock or drying of brittle materials. We use a quasi-static gradient damage model to perform large scale numerical simulations showing that the propagation of fully developed cracks follows Griffith criterion and depends only on the fracture toughness, while crack morphogenesis is driven by the material's internal length. Our numerical simulations feature networks of parallel cracks and selective arrest in two dimensions and hexagonal columnar joints in three dimensions, without any hypotheses on cracks geometry and are in good agreement with available experimental results.

We study the genesis and the selective propagation of complex crack networks induced by thermal shock or drying of brittle materials. We use a quasi-static gradient damage model to perform large scale numerical simulations showing that the propagation of fully developed cracks follows Griffith criterion and depends only on the fracture toughness, while crack morphogenesis is driven by the material's internal length. Our numerical simulations feature networks of parallel cracks and selective arrest in two dimensions and hexagonal columnar joints in three dimensions, without any hypotheses on cracks geometry and are in good agreement with available experimental results. Complex crack patterns are ubiquitous in nature and in technology applications. Yet the theoretical understanding and predictive numerical simulation of how and when complex crack patterns arise (nucleation) and how they evolve (crack propagation) is fraught with challenges. Although approaches based on phase fields [1] or variational regularizations [2] have led to significant advance in the numerical simulation of complex crack patterns, short of introducing initial flaws at the structural scale [3], prescribing ad-hoc stress criteria [4], or accepting global energy minimization arguments whose physical relevance is debated [5][6][7], the predictive understanding of crack nucleation is still an elusive goal.
It is well-accepted that while Griffith-like models are appropriate for crack propagation at the scale of a structure, they are inadequate for the modeling of crack nucleation in brittle materials. Arguably, finer models, where a microscopic (material) length scale plays a fundamental role, are necessary to determine the critical load and crack geometry at the onset, especially in situations where complex crack patterns arise straight from the nucleation. The consistent combined modeling and numerical simulation of crack nucleation and propagation from the material to the structural length-scale is a challenging and largely open issue.
In this Letter, we study the morphogenesis and the selective growth of complex crack patterns induced by material shrinking under thermal shock. We report unprecedented quantitative agreement between numerical simulations, a theoretical model, and experiments at scales spanning from the material internal length to the structural length-scale. Our numerical simulations predict key features of fracture patterns observed in experiments, such as the formation of periodic patterns and the scaling laws governing their selective propagation in two and three dimensions, and do not require any

Numerics
Semi-analytical scaling law from [9] (c) a priori hypotheses on cracks geometry. The method we use leverages recent progress in the understanding of the links between damage models [10,11] and the variational approach to fracture [2,12]. It is based on a rate-independent gradient damage model with stress-softening based on two material parameters: the fracture toughness which rules the evolution of fully developed cracks, and the material's internal length arXiv:1310.0501v2 [cond-mat.mtrl-sci] 28 Nov 2013 which controls the initial stages of crack nucleation.
We investigate the thermal shock of a brittle ceramics, a now classical experimental setup [3,13,14] where a sample initially at a uniform temperature T 0 is quenched in a cold bath at temperature T 0 − ∆T . We consider a rectangular slab Ω exposed to the thermal shock through its thin faces only. We focus first on very thin slabs, which we represent by a two-dimensional body in plane stress. We assume that within the range of temperatures involved, the material properties remain constant. Denoting by u the displacement field and ε = (∇u+∇ T u)/2 the linear strain tensor, we consider for the sound material a linear elastic behavior of energy density ψ t (ε) = A 0 (ε − ε th t ) · (ε − ε th t )/2 where A 0 is the isotropic elastic stiffness tensor. The inelastic deformation induced by the time-dependent temperature field T t is ε th t = β(T t − T 0 )I, where I is the identity matrix. The index t is meant to highlight the dependence on time. We neglect the cracks influence on heat transfer so that the temperature field T t solves the heat equation Phase changes, non-uniform convection and other nonlinear aspects of the heat exchange between the fluid and the sample are neglected by assuming that the temperature of the domain boundary exposed to the thermal shock is constant and equal to that of the water bath, i.e. T t = T 0 − ∆T . Inertial effects are not considered because the diffusion velocity of the temperature field is much slower than the wave speed in the material at the relevant scales in time and space. This hypothesis is universally accepted in the literature on thermal shock problems [3,4,9,14,15]. We model material failure using a gradient damage model characterized by the energy function where α is a scalar damage field varying between 0 (sound material) and 1 (fully damaged material), G c is the material's fracture toughness, an internal length, and c w = 1 0 w(s) ds a normalization constant. In a timediscrete setting, the quasi-static evolution is obtained by solving at time t i the following minimization problem min u,α≥αi−1 E ti (u, α), where the unilateral constraint on α enforces the irreversibility condition on the damage. The compliance function s and the energy dissipation function w should be chosen such that (1) converges as → 0 to a Griffith-like energy Ω\Γ ψ t (ε) dx + G c S(Γ), where S is the surface measure of the crack Γ [2,17,18]. In this model, material interpenetration in the fully damaged area is possible. In all the simulations presented here, it can be checked a posteriori that this issue does not present itself. Here, we use s(α) = 1/(1 − α) 2 and w(α) = α, a choice motivated by the convenience of its numerical implementation and specific analytical studies [11,19]. With this choice the damage model has a stress-softening behavior and remains purely elastic without damage until the stress reaches the critical value: The relation above may be used to determine the numerical value of the internal length for a specific material from the knowledge of its elastic limit σ c , Young modulus E, and fracture toughness G c [19]. The present model is in many aspects similar to the phase-field models of fracture developed independently [20]. Those with singlewell dissipation potentials [1,21] are in the form of (1) with One significant difference is that while phase-field models typically involve some form of viscous regularizations, our formulation is rate-independent. In addition, the current literature based on phase-field models is concerned only with the propagation of a pre-existing cracks and does not consider the initiation problem. The dimensional analysis of the energy (1) highlights three characteristic lengths: the geometric dimension of the domain L, the internal length and the Griffith length 0 = G c / Eβ 2 ∆T 2 . Using the material's internal length as the reference unit, the problem can be reformulated in terms of two dimensionless parameters, the dimension of the structure L/ (a geometric parameter) and the intensity of the thermal shock 0 / (a loading parameter). This is a significant departure from the classical Griffith setting where the only relevant parameter is L/ 0 [4,5,9].  (2) gives = 46 µm and 0 = 14 µm. As our model is rate independent its solution are independent of k c , up to a change of time scale. The numerical results are obtained through a finite element discretization and the approach of [2,19,22]. The main technical difficulties are the constrained minimization of a non-convex energy, and the need for a spatial discretization adapted to the material length-scale scale . Cracks correspond to the localized bands where α goes from 0 to 1 and back to 0. The qualitative agreement between experiments and simulation is very good. In particular, our simulations reproduce the key phenomenon: the emergence of a periodic array of parallel short cracks at the initiation and their selective propagation toward the interior of the slab. Figure 1(c) shows a quantitative comparison between the numerical simulation of FIG. 1(a) and experimental data from [8] by plotting the average crack spacing d as a function of the distance a to the edge exposed to the thermal shock for the final configuration, the agreement is striking. Note that in the experimental results shorter cracks are probably filtered out by the adopted experimental crack detection methods [19]. In a first regime, very short equi-distributed cracks nucleate (the plateaus of the crack spacing for short depth in the numerical experiments), followed by selective arrest and period doubling [23]. In the central region of the plot, we can compare experimental and simulation data with a scaling law obtained in [9] through linear fracture mechanics calculations by imposing a bifurcation condition between crack propagation modes with period doubling or not (solid line). For larger values of a, we observe the final crack arrest caused by the finite size of the sample, again in very good agreement with the experiments. Whereas classical theories can be applied in the second and third regimes consisting of fully developed cracks, they cannot properly account for the nucleation phenomenon observed here without preexisting flaws. Our simulations are initialized with a null damage field, an homogenous material, and an unflawed geometry. The crack nucleation is due to the softening character of the material behavior.
The second series of simulations focuses precisely on the crack nucleation process and hence on short times. In this setting, one can assume that the domain is semiinfinite so that the geometric parameter L/ is infinite and the only parameter is the intensity of the thermal shock 0 / . For an undamaged material, the stress is uniaxial and reaches its maximum value σ max = Eβ∆T at the surface of the thermal shock. Since (σ max /σ c ) 2 = 3 /8 0 , for mild-enough thermal shocks ( 0 / > 8/3), the critical stress is never reached and the solution remains elastic at all time. If 0 / < 8/3, damage takes place at t = 0, is homogeneous in the horizontal direction, and non-null in a band of finite thickness D, which penetrates progressively inside the body until a critical time t * . At t = t * the horizontally-homogeneous solution becomes unstable and the damage field develops oscillations of periodicity λ * (FIG. 2(a)).
An analytical solution for the damage field in the first stage of the evolution and its bifurcation and stability analysis is reported in [24], providing semi-analytical results for the periodicity λ * , the damage penetration D * , and the time t * at the bifurcation. Here we perform several simulations varying 0 / and detect the critical parameter at the bifurcation. In Figure 2(b) the numerical simulations (dots) are compared to [24] (solid lines). The good agreement provides an excellent verification of our numerical model. For severe shocks ( 0 ), the results disclose a welldefinite asymptotic behavior with λ * ∼ √ 0 , D * ∼ , and t * ∼ 0 /k c . In this regime we observe numerically that all oscillations at the bifurcation develop in fully formed cracks (max α = 1), which is not the case for milder shocks ( 0 ∼ ). However, the full post-bifurcation analysis remains an open problem at this time.
Experimental studies show that in three-dimensions cracks delineate cells with coarsening polygonal crosssections [25]. Because of the complexity of the problem, the few available theoretical and numerical studies are based either on simplified two-dimensional models [16,26] or on strong assumptions on the crack geometry [27]. The numerical simulation and analysis of the full three-dimensional problem is a major challenge for classical fracture mechanics tools and remains therefore largely unexplored.   the computational cost, the computation is performed on a fragment of width 5 mm and height 1 mm of the domain, and the temperature is assumed constant throughout the thickness of the sample. We observe a transition from transverse cracks to three-dimensional fracture patterns delimiting polygonal cells between 480 • C and 580 • C, which is consistent with [8, Fig 5]. Another series of simulations performed at constant temperature for increasing sample thickness (not shown here) highlight the same behavior: transverse cracks for thin domains, transitioning to three-dimensional cracks for thicknesses between 1 mm and 2 mm. This is also consistent with the observations in [8, Fig 5] and justifies the use of a two-dimensional model, a posteriori. Figure 4 shows a fully tri-dimensional crack pattern obtained for a domain of dimension 150 × 150 × 20 for 0 = 0.05 . During the simulation, a disordered pattern of small cells nucleates in the first time steps and propagates quasi-statically inside the domain. A selection mechanism leading to honeycomb patterns with increasingly large and regular cell arises from energy minimization. Tracking the propagation of three-dimensional crack front of Figure 4 using a classical Griffith-based model requiring an explicit description of the crack surface and its propagation criterion would be prohibitively complex. Instead, our three-dimensional computations are performed through a straightforward extension of the discretization and minimization algorithm for the energy (1). Obtaining an accurate scaling law for the cell diameter as a function of the depth as the one of the two-dimensional case of Figure 1 would require simulations on larger domains which would rapidly become computationally prohibitive. Hence, we perform a series of numerical experiments by fixing the loading 0 and by varying at the same time the internal length and the domain size L so as to keep their ratio equal to L/ = 40. The number of element is kept constant with a mesh size h = /5. For each computation, we compute the average cell diameter d as a function of distance from the bottom edge a using a post-processing software. This process does not involve any adjustable parameter, yet our results match the two-dimensional scaling law of [9] over several orders of magnitude, leading us to conclude that the scale selection mechanism in two and three dimensions are identical. In addition, while the initial phase of the evolution depends strongly on , later time evolution of fully developed cracks at the structural is unaffected by this parameter, matching the general scaling law for a Griffith-based model. This finding is consistent with the properties of the energy functional E t which is known to lead to a Griffith-type propagation criterion [18,21,28,29].
Our simulations show that a purely quasi-static model based on energy minimization can fully explain the formation of imperfect polygonal patterns and their selective coarsening as "maturation" mechanism during propagation, a phenomenon sometimes attributed to nonequilibrium processes [25]. We show that a carefully chosen gradient damage model can be used to account simultaneously for the nucleation of complex crack patterns and their propagation following Griffith criterion. Further works will be carried on to perform a careful statistical analysis of the geometry of the 3D crack patterns and further comparisons to experimental results. The present modeling framework has a general validity and can be applied to other domains including for example the formation of basalt columns with uniform cross-sectional diameters through the solidification of lava fronts [25] or shaping of biological systems as observed of the scales on the heads of crocodiles [30]. We are also considering stronger thermo-mechanical coupling including the effect of cracks on heat transfer as in [15,31].
The authors wish to thank Yingfeng Shao for providing the experimental data used in FIG. 1. B.B. work was supported in part by the National Science Foundation grant DMS-0909267. Some numerical experiments were performed using resources of the Extreme Science and Engineering Discovery Environment (XSEDE), which is supported by National Science Foundation grant number OCI-1053575 under the Resource Allocation TG-DMS060014. J.J.M. and C.M. gratefully acknowledge the funding of the ANR program T-Shock OTP J11R087.