{'fraction_non_alphanumeric': 0.05396624138204295, 'fraction_numerical': 0.01022268008558523, 'mean_word_length': 4.130920918885953, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 17, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A quasi-static model for craquelure patterns

We consider the quasi-static evolution of a brittle layer on a stiff substrate; adhesion between layers is assumed to be elastic. Employing a phase-field approach we obtain the quasi-static evolution as the limit of time-discrete evolutions computed by an alternate minimization scheme. We study the limit evolution, providing a qualitative discussion of its behaviour and a rigorous characterization, in terms of parametrized balanced viscosity evolutions. Further, we study the transition layer of the phase-field, in a simplified setting, and show that it governs the spacing of cracks in the first stages of the evolution. Numerical results show a good consistency with the theoretical study and the local morphology of real life craquelure patterns.


Introduction
Craquelure and crazing usually denote network of cracks which appear on thin superficial layers of materials (see Figure 1). As a mathematical model for this kind of phenomena, we consider an elastic, brittle layer placed on a rigid adhesive substrate, which displaces the layer (and thus plays the role of a driving force). Having in mind the formation of patterns in a long time scale, we consider a quasi-static evolution driven by a time depending displacement of the substrate; in particular, we neglect diffusion of temperature and inertia. We employ a phase-field approach and we focus on a couple of complementary aspects: • characterization of quasi-static evolutions, obtained by time-discrete alternate minimization (staggered) schemes, • pattern formation and crack spacing as the result of different generations of cracks.
Other interesting aspects have been recently studied for similar mechanical systems, e.g. [5,1,8].  We assume that the reference configuration of the brittle film is an open bounded set Ω ⊂ R 2 . The physical parameters are the Lamè coefficients λ > 0 and µ > 0, the fracture toughness G c > 0 and the adhesion parameter β > 0 (which measures the elastic response of the adhesive between the brittle layer and the substrate). Finally, we consider an in-plane displacement g which gives, as a function of time, the displacement of the substrate (note that here g is a datum). For sake of simplicity we assume that g(t) is of the form tĝ for someĝ ∈ H 1 (Ω, R 2 ).
The sharp crack energy associated to the system is, in some sense, a combination of Griffith energy [10] and Mumford-Shah functional [15]; in the (weak) setting of SBD 2 spaces it reads where W (depending on λ and µ) denotes the linear elastic energy density. In this framework the set of discontinuity points J u represents the crack.
Actually, in view of the numerical simulations, it is far more convenient to work with the phase-field approximation of F given by [6] F ε (t, u, v) = 1 2 ∫ Ω (v 2 + η ε )W( (u)) dx + It well known [9] that the Γ-limit of the phase-field energy F ε (t, ·, ·) is indeed the sharp crack energy F (t, ·). Now, let us introduce the time-discrete evolution which is used in the numerical simulations. Let t k = kτ for τ > 0. Known the configuration (u k−1 , v k−1 ) at time t k−1 the updated configuration (u k , v k ) at time t k is given by v k = lim m→+∞ v k,m and u k = lim m→+∞ u k,m , where u k,m and v k,m are the steps of the staggered scheme: The theoretical question we are interested is the characterization of the time continuous limit. A precise statement is contained in §2.3, asuuming a stronger irreversibility condition, i.e. v ≤ v k,m−1 . Roughly speaking we will get an evolution t → (u(t), v(t)) which is possibly discontinuous and such that (u(t), v(t)) is a critical point of the energy, for every continuity point t. The technically difficult part, is instead the characterization of the behaviour in the discontinuity points, where the system makes an instantaneous (catastrophic) transition between (u(t − ), v(t − )) and (u(t + ), v(t + )) following some sort of staggered scheme.
We finally provide some numerical results, which shows the pattern formation for a prototype problem and for the real life specimen of Figure 1. It turns out that the pattern is formed by generations of cracks (which nucleate at different time) whose spacing follows a regular scheme: in the first generation, spacing of cracks is not predicted by minimizers (unless the domain is very short), it is instead dictated by the size of the boundary layer of the strain and of the phase-field function; further generations of cracks follow instead the dyadic behaviour of minimizers.
Now, let us turn to the Γ-limit of F ε (t, ·, ·) as ε → 0 + . To this end, it is convenient to define the extended functionalF ε : Then, by [9] we known that the Γ-limit of F (t, ·, ·), with respect to the topology of Let (u ε , v ε ) be a family of minimizers of F ε (t, ·, ·). Note that u ε is bounded in L 2 , then, by a fundamental result in the theory of Γ-convergence, see e.g. [7, Theorem 3.3], we know that u ε converge to a minimizer u of F (t, ·); more precisely, there exists a subsequence converging to u in the topology of L 2 (Ω; R 2 ). On the contrary, not much is known on the convergence of the critical points of F ε (t, ·, ·) which will appear in the evolution.

Time-discrete evolution
In this section we will describe the time-discrete evolution on which the numerical calculations are based. Here we can give only a brief description in the simplest possible setting, for complete proofs and generalizations the reader should make reference to [12,4]. For τ > 0 consider a time discretization t k = kτ. After setting the initial conditions u 0 and v 0 , at time t = 0, the configuration (u k , v k ), at time t k , is computed from (u k−1 , v k−1 ) by the following incremental scheme, known as alternate minimization algorithm. Let us introduce the further sequences u k,m and v k,m , with v k,0 = v k−1 and u k,0 = u k−1 . Then, for m ≥ 1 we define Note that this scheme takes full advantage of the fact that F (t k , ·, ·) is separately quadratic, and indeed its numerical implementation is very convenient, even if in practice the algorithm may converge quite slowly. Then, (up to subsequences) we let v k = lim m→+∞ v k,m and u k = lim m→+∞ u k,m . Note that the distance between v k−1 and v k could be either small or large, in other terms, there is no a priori control on the speed (v k − v k−1 )/τ since the system is rate-independent. Moreover, note that (u k , v k ) is an equilibrium point for F (t k , ·, ·), unilateral w.r.t. v. Finally, it is important to comment on the constraint v ≤ v k−1 which models irreversibility. The convergence analysis under the latter constraint in the "genuinely" rate-independent setting is still open. At the current stage a couple of alternative are feasible. The first consists in solving a discrete viscous parabolic evolution for the phase-field, with the contraint v ≤ v k−1 , followed by a vanishing viscosity procedure: by [3] and [2] this approach provides in the end a quasi-static evolution. The second option consists instead in replacing the constraint v ≤ v k−1 by the (stronger) constraint v ≤ v k,m−1 . In this way monotonicity is imposed at each alternate iteration. On the theoretical level, this assumption allows to give a full characterization of the evolution as the time step τ → 0, without passing through viscosity solutions. On the other hand, in some cases this assumption may be too strong in the numerical simulations; we will see an example in §5.3. Here we will follow the latter strategy.

Time-continuous evolution
The time discrete scheme of the previous subsection gives for every τ > 0 a finite sequence (v k , u k ) in the points t k . We will denote by (u τ , v τ ) : [0, T] → H 1 (Ω; R 2 ) × H 1 (Ω) an interpolation of (u k , v k ) in the points t k (here we will enter into the delicate technical issue about the choice of the interpolation). Our goal is to characterize the limit of (u τ , v τ ) as τ → 0.
First, we give a "qualitative" description of the limit evolution, as a function of time. We have already remarked that in this setting there is no a priori control of . In other terms the limit evolution t → v(t) could be discontinuous in time. Note that this is observed in the numerical simulations. For this reason, in the limit as τ → 0 we expect an evolution t → (u(t), v(t)) of class BV in time. Now, let us briefly describe the behaviour in continuity and discontinuity points (for more details we refer to [4,12]).
If t ∈ [0, T] is a continuity point for the evolution then (u(t), v(t)) is an equilibrium configuration for the system, i.e.
Equilibrium of v is unilateral because, by irreversibility, only negative variations are allowed. Note that the configuration (u(t), v(t)) is not necessarily a (joint) global or a local minimizer of F ε (t, ·, ·), it is actually a separate minimizer, by the separate convexity of F ε (t, ·, ·).
On the contrary, if t is a discontinuity point we expect an alternate "evolution", connecting (u(t − ), v(t − )) and (u(t + ), v(t + )). More precisely, in the limit the path between In particular the instantaneous transition between (u(t − ), v(t − )) and (u(t + ), v(t + )) is not simultaneous in u and v. Once again, this is confirmed by numerical results.
In order to give a rigorous description of the limit evolution, it is necessary to introduce the derivatives of the energy F ε , which will provide the driving forces for the evolution. We follow in particular [4]. Clearly, we can take the partial derivatives of F ε ; moreover, we can define the following slopes Note that the slope w.r.t. v is unilateral since, by irreversibility, we are interested only in negative variations. From the technical point of view it is fundamental that the slopes are weakly lower semicontinuous. Technically, we characterize the evolution in terms of a parametrized "balanced viscosity solution" [14]. First, in order to describe both the behaviour in continuity and discontinuity points we employ an arc-length parametrization s → (t(s), u(s), v(s)). In this setting, roughly speaking, discontinuity points (in time) correspond to intervals [s − , s + ] where t is constant, and thus t = 0; vice-versa, continuity points (in time) correspond to points s where t (s) > 0. Then, the evolution is characterized by the following set of conditions: (a) the map t : [0, S] → [0, T] is non-decreasing and surjective, (e) for every s ∈ [0, S] the following energy identity holds: For further details the reader may refer to [4].

A one-dimensional case study
As a preliminary study, it is interesting to focus on the behaviour and on the pattern generated by the minimizers of the sharp-crack energy, even if the evolution follows only partially this scheme. To this end it is useful to consider a simpler one-dimensional setting, where Ω = (−L, L) and g(t, x) = t x; the corresponding energy is then of the form We will assume that the initial configuration (at time t = 0) is g(0) = 0.
Then, a simple (but quite long) computation shows that In the sequel it will be important we understand the dependence of u 0 on L.
for different values of L (see Figure 4) shows that u L changes in the boundary layers while it is almost constant in the "interior". In order to provide some quantitative estimate of the size of the layer it is convenient to study how the derivative u L (x) scales with L. Thus, given a point x ∈ (0, L ) let us compute the (corresponding) point x ∈ (0, L) such that u L (x) = u L (x ). The latter identity reads A plot of the behaviour of x as a function of L (see Figure 5) shows that the distance between x and L is almost independent of L. Thus, we can expect that the length of the boundary layer remains almost constant with L. As we will see in the numerical simulations of § 5, the "scale invariance" of the boundary layer explains the periodic spacing of the first generations of cracks better than the dyadic structure (explained in the next subsections), which occurs only for small values of the size L (more precisely, when L is comparable with the size of layer). In order to better explain the fact that L − x is almost constant, we can consider the case when x and L are large enough to approximate cosh with exp/2, then the above identity can be approximated by which gives x ∼ L − b. Roughly speaking, the transition layer is of constant size b.

First generation of cracks by minimality
The content of this section follows the study performed in [5]. In the case of a single crack let us assume (for the moment) that J u = {0}, so that the interval (−L, L) splits into the subintervals (−L, 0) and (0, L); the minimizer, denoted by u 1 (t), can then be computed using the previous Euler-Lagrange equation in the subintervals (−L, 0) and (0, L) (see again Figure 3). Then, the energy reads F(t, u 1 (t)) = G c + t 2 2F(L/2). Now we will compare the energy F(t, u 0 (t)) of the continuous solution with the energy F(t, u 1 (t)) of the discontinuous solution. To this end, it is interesting to consider in the (t, L) plane the set of "critical transition times", satisfying It is easy to see that the function vanishes for L = 0 and is increasing for L > 0; hence it is positive and, given L, there exists a time for which the two energies coincide. Further, t 2 (F(L) − 2F(L/2)) < G c for t < t L while t 2 (F(L) − 2F(L/2)) > G c for t > t L ; in other terms, which justifies the name "critical transition time" at length L.
Before proceeding, let us check that for a single crack the least energy is always assumed when J u = {0}. For L 1 + L 2 = L, we have to check the energy inequality After some algebraic manipulations, the previous inequality boils down to 1 2 (tanh(λL 1 ) + tanh(λL 2 )) < tanh(λ(L 1 + L 2 )/2), which is true by the strict concavity of tanh.
Now, let us consider the general case of m cracks (for m ≥ 1). By the previous symmetry argument it is not restrictive to assume that J u = {−L + k L/(m + 1) : for 1 ≤ k ≤ m}. Arguing as above, the energy for m cracks reads We want to compare F(t, u 0 (t)) and F(t, u m (t)). This time the set of critical points is defined by t 2 F (L) − (m + 1)F(L/(m + 1)) = mG c , which gives the "critical transition time" .
Let us show that the sequence t m is monotone increasing w.r.t. m, independently of L. To this end it is convenient to study the function s → (s − 1) s tanh(λL/s) − tanh(λL) which coincide with t 2 m (up to a positive multiplicative constant) when s = m + 1. A simple calculation shows that this function is increasing and thus t m is monotone increasing as well. Note that t 1 = t L .
As a consequence, if the evolution is driven by energy minimization, in the first generation it appears a single crack in the center of the bar (−L, L). Indeed, for t ∈ [0, t 1 ) we have F(t, u 0 (t)) < F(t, u 1 (t)) and thus u(t) = u 0 (t). For t > t 1 we have F(t 1 , u 0 (t)) > F(t 1 , u 1 (t)), moreover, being t m > t 1 for every m > 1, we have F(t, u m (t)) > F(t, u 0 (t)), at least for t ∼ t 1 . Hence, we expect that u(t) = u 1 (t) for t ∈ (t 1 , t 2 ) where t 2 denotes the time when the second generation of crack will occur.

Second and further generations of cracks by minimality
When the first crack appears, at time t 1 , the interval (−L, L) is splitted into the subintervals (−L, 0) e (0, L). By periodicity, this is equivalent to consider the behaviour of solutions in the interval (−L/2, L/2). Arguing exactly as in the previous section, just replacing L with L/2, provides the existence of the critical transition time t L/2 when the bar will split again. Note that t L/2 > t L because the function t L in (8) is decreasing w.r.t. L (remember that ∆ 2 is increasing with respect to L).
Proceeding by induction, we find the times t L/2 k (see Figure 6) when the dyadic structure evolves; in the physical literature this is often called "halving" effect.
To conclude this section, it is fair to remark that in this one dimensional example the transition between u L and u L/2 cannot occur following a continuous energy decreasing path, because of the activation threshold G c which is payed as soon as the crack opens. In the numerical simulations, this topological problem is avoided by the phase-field regularization. Moreover, the uni-axial numerical results of §5 we show that the first generation of cracks does not follow this dyadic structure (at least for L large), rather the first crack pattern depends on the size of the boundary layer, described in the previous subsection.

Alternate minimization in the one dimensional setting
In this section we briefly study the evolution generated by the alternate minimization scheme in the one-dimensional setting of §3. For sake of simplicity we consider Dirichlet, instead of Neumann, boundary conditions. This example will be useful to understand the behaviour of the numerical results. We consider again the sharp crack energy F of (5) and the corresponding phase-field energy Here, we will further assume that u(t, ±L) = g(t, ±L) = ±t L. Fix t k = kτ assume that v k−1 = c k−1 is constant. We will show hereafter that the update v k is constant and that u k is constant, as well. Looking at this result in the time discrete scheme, it turns out that for Dirichlet boundary conditions there will be no nucleation of cracks if the initial phase-field v 0 is homogeneous. Other interesting results on homogeneous states are contained in [17]. For sake of simplicity, we will solve the following alternate scheme without the irreversibility constraint, which would actually not change the qualitative result, The Euler-Lagrange equation is of the form −a k−1 u +2βu = 2βg(t k ) in (−L, L) with Dirichlet boundary conditions. The explicit solution is simply u k,m (x) = g(t k , x) = −t k x independently of a k−1 . Since u k is constant the Euler-Lagrange equation for v is of the form −v + bv = c, whose solution v k,m is again a constant. In particular the staggered scheme finds a critical point after two iterations.

Numerical results for uni-axial problems
In the following subsections are reported the numerical results obtained for a bar Ω = (−L, L) × (−H, H) of variable length and fixed width 2H = 5. We have chosen fracture toughness G c = 1.0, Young modulus E = 1.0, Poisson ratio ν = 0.15, adhesive constant β = 0.15. Moreover, we assume that the displacement of the substrate is of the form g(t, x) = (t x 1 , 0). Since g is uni-axial we expect solutions in accordance with the theoretical arguments of § 3. Since the datum g is monotone we neglected the irreversibility constraint in the alternate minimization scheme (4). With this choice crack patterns look sharper. A comparison between the solutions obtained with and without irreversibility constraint (see § 5.3) show that the crack patterns behave in a similar way. Numerical results have been computed using FreeFem++ [11].

Dyadic structure for short bars
In this subsection we present and discuss the numerical results obtained with L = 6.5. The configuration of the phase field v shown in Figure 7 are computed with the unconstrained version of scheme (4), i.e. without any constraint on v.
First, it is interesting to comment on the very first image. In this case there are no pre-existing fractures, hence nucleation is a fundamental ingredient in the evolution. Figure 7 shows that nucleation requires first a diffuse damage, with the phase-field v decreasing in a wide region, and then concentration, in a second stage (see also [13,Fig 4.8]). Note that, in the first image, the profile of v is closely related to the profile of u L (the derivative of the one dimensional solution u L ) presented in §3.1. Indeed, since the datum is uniaxial and Poisson ratio ν is small, we can expect that u L provides a good approximation of the strain energy density; then, in phase-field models the larger is the elastic energy density the smallest is the phase-field. This is why the profile of v is consistent with that u L (compare the first plot of Figure 4 with the first image of Figure 7).
Then, the central crack nucleates and the domain splits. At that point, a similar evolution "restarts" in each subinterval, leading to the second generation of cracks.
Note that in this case the "transition layer" for u L is, roughly speaking, of the same size of the half interval (0, L). For this reason the evolution generates a dyadic fracture pattern, clearly visible in the last image of Figure 7. As we will see, the pattern of the first generation is not dyadic for larger L.

Periodic patterns for longer bars
In this subsection we present the numerical results obtained for L = 12.5. In the second picture, when the first generation of crack appears (at time t = 2.8) it is evident that the position of the first crack (from left to right) corresponds to the transition region of v, which, comparing with Figure 4, corresponds to the transition layer of u L . Moreover, all the cracks, apart from the innermost, are equally spaced and their distance is approximately 5.5, which is the value expected from the estimates of the size of the transition layer (see Figure 5). Note that, following the intermediate steps of the alternate minimization algorithm (at time t = 2.8), the outermost cracks nucleates firts, followed by the innermost. The reason behind the fact that cracks do not nucleate in the center of the bar is not fully clear, however, it should be found in the fact that the phase field v is almost constant in the inner part of the domain and thus its evolution does not promote concentration of strain and nucleation of cracks, as described in §4. The second generation of cracks follows instead the dyadic structure: this is due to the fact that the spacing of the cracks after the first generation is comparable with the size of the transition layer, as in the previous example. In conclusion, we could say that the first generation of cracks is determined by a characteristic length (the width of the transition layer) with cracks nucleating from the boundary to the interior, while the further generations evolve according to the dyadic scheme. Numerical results for longer bars confirm this behaviour. Finally we remark this type of evolution, triggered by the boundary layer, has been observed also in numerical simulations of brittle layered materials, see [1, Fig. 7], and soil drying (which share several features with our problem), see [18, Fig. 6].

Irreversibility
In this section we briefly report the numerical results obtained with the scheme (4) under the constraint v ≤ v k−1 in the case L = 6.5. Clearly, the main difference is the fact that here the phase-field v cannot increase after the nucleation of cracks, because of the irreversibility constraint. From the results, it is evident that this evolution is very similar to the evolution of Figure 7, at least as far as the structure of the crack pattern, the time and the nucleation sites.
In the case L = 12.5 the behaviour is similar to that of Figure 8.

Local craquelure patterns for a real life specimen
In order to produce some more realistic craquelure patterns and crack morphologies, we ran several numerical experiment changing both the value of the adhesive parameter β and the datum g, of the form g(x) = Ax for A ∈ R 2×2 . Even if the values are not provided by experimental measurements, it is clear from Figure 10 that this quasi-static phase-field model captures quite well the local morphological features of the patterns in Figure 1. Note that these plots do not show the entire evolution leading to the formation of the patterns but they show the final snapshot, at a certain time T. In the evolution the cracks actually appear at different times, following a scheme reminiscent of the uni-axial case.
In the bi-axial setting it is more difficult to explain the behaviour of the evolution. However, some observations are due as far as crack junctions, which of course did not occur in the uniaxial setting. First, remember that for local minimizers of the (scalar) Mumford-Shah functional only Y-junctions (i.e. triple junctions) at 2 3 π angles occur, see [15,Theorem 2.1]. We may expect a similar behaviour in our linear elasticity context. However, in the first image only T-junctions occur, this is simply due to the fact that horizontal cracks nucleate in the first generation and vertical cracks in the second generation, therefore, once horizontal straight cracks are formed it is no longer possible to have 2 3 π angles. In this example, it would be natural to see orthogonal crossing of cracks; on the contrary cracks seem to "shift" passing from one horizontal stripe to the other. This feaure could be due to remeshing and to the fact that horizontal cracks are not equally spaced, however, it is interesting to note that this feature occurs also in several points of the real life picture. In some cases (see the pictures on the right) the crack behaves instead as a local minimizer, splitting a T-or an X-junction into three (or more) Y-junctions, with angles close to the optimal value 2 3 π, compare with [15, Figure 14]. In this way, small triangular regions are formed. Once again, it is worth to remark that this the case both in the numerical experiment and in the real life specimen.