Crack channelling mechanisms in brittle coating systems under moisture or temperature gradients

Crack channelling is predicted in a brittle coating-substrate system that is subjected to a moisture or temperature gradient in the thickness direction. Competing failure scenarios are identified, and are distinguished by the degree to which the coating-substrate interface delaminates, and whether this delamination is finite or unlimited in nature. Failure mechanism maps are constructed, and illustrate the sensitivity of the active crack channelling mechanism and associated channelling stress to the ratio of coating toughness to interfacial toughness, to the mismatch in elastic modulus and to the mismatch in coefficient of hygral or thermal expansion. The effect of the ratio of coating to substrate thickness upon the failure mechanism and channelling stress is also explored. Closed-form expressions for the steady-state delamination stress are derived, and are used to determine the transition value of moisture state that leads to unlimited delamination. Although the results are applicable to coating-substrate systems in a wide range of applications, the study focusses on the prediction of cracking in historical paintings due to indoor climate fluctuations, with the objective of helping museums developing strategies for the preservation of art objects. For this specific application, crack channelling with delamination needs to be avoided under all circumstances, as it may induce flaking of paint material. In historical paintings, the substrate thickness is typically more than ten times larger than the thickness of the paint layer; for such a system, the failure maps constructed from the numerical simulations indicate that paint delamination is absent if the delamination toughness is larger than approximately half of the mode I toughness of the paint layer. Further, the transition between crack channelling with and without delamination appears to be relatively insensitive to the mismatch in the elastic modulus of the substrate and paint layer. The failure maps developed in this work may provide a useful tool for museum conservators to identify the allowable indoor humidity and temperature fluctuations for which crack channelling with delamination is prevented in historical paintings.


Introduction
The prediction of the fracture response of a bilayer composed of a brittle coating adhering to a substrate is highly relevant to the structural integrity of microelectronics components (van Gils et al. 2004), thermal barrier coatings (Choi et al. 1999;Hille et al. 2009), ceramic multi-layers (Sørensen et al. 2012), road pavements (Timm et al. 2003), natural phenomena (surface layers in dry soils (Velde 1999)), and also to the failure and preservation of cultural heritage objects including historical oak wood cabinets (Luimes et al. 2018a) and historical paintings (Giorgiutti-Dauphiné and Pauchard 2016; Léang et al. 2017). This study is mainly motivated from observations of surface crack development in historical paintings, as driven by indoor climate variations.
Historical paintings typically consist of one or more layers of paint adhering to a canvas or wooden substrate, and their degradation is sensitive to environmental changes associated with moisture fluctuations, and to a lesser extent to temperature changes (Mecklenburg and Tumosa 1991). Degradation originates from differential swelling and shrinkage across the layers: indoor climate variations, along with differential mechanical and hygro-thermo-expansive properties of the paint and substrate materials, result in tensile self-stresses in the paint layer. If the tensile stress at some location of the paint layer attains the fracture strength, then paint cracking will occur. The most common mechanism of paint cracking is crack channelling, as observed at the surface of panel and canvas paintings, and is usually denoted as craquelure. The network of these cracks is complex, and generally depends on the geometry and materials used and on the environmental conditions to which the painting is exposed (Bucklow 1997). Moisture-driven channelling cracks in wooden panel paintings have been analysed (Bratasz 2013;Krzemień et al. 2016;Bratasz and Vaziri Sereshk 2018), and the sensitivity of crack stability to micro-climate variations, and to the thickness of layers have been reported. A channelling crack may kink at the interface between the paint layer and the substrate, generating interfacial delamination. This mechanism has also been identified in paintings (Crook and Learner 2000;Young 2007), and its occurrence under the influence of relative humidity cycles has been investigated in Tantideeravit et al. (2013), Wood et al. (2018). If the delamination is unstable such that it grows without limit, it may ultimately lead to spallation or flaking of the paint material, as reported in Brewer and Forno (1997), O'Donoghue et al. (2006). Flaking degrades the visual appearance of the painting, and needs to be understood in detail for museums to adequately regulate their indoor climate control and safely preserve their objects.
An illustrative overview of the three crack channelling mechanisms as described above is given in Fig. 1. In the lower left corner of the painting displayed in Fig. 1a, a fine network of channelling cracks can be observed at the paint surface, see Fig. 1b. This may lead to crack deflection at the paint-substrate interface, characterized by two opposite delaminations, see Fig. 1c. A delamination crack may advance without limit, thereby triggering paint flaking, see Fig. 1d.
From a more general perspective, the study of cracking in coating-substrate systems has been the topic of various investigations, as follows. The channelling behaviour of a single mode I crack in a brittle, isotropic elastic layer subjected to residual tension and bonded either to a semi-infinite or finite-thickness substrate has been analysed in Beuth (1992) and Vlassak (2003), respectively. These analyses show that the reduction in film stress associated with elastic deformation of the substrate may have a significant effect on the energy release rate for crack propagation. The channelling of a series of parallel cracks in a surface layer has been addressed in Thouless (1990), Thouless et al. (1992), and the crack spacing has been predicted as a function of the stress in the surface layer, its thickness, and the toughness of the layer. A similar problem has been studied for the case of thermal loading in Timm et al. (2003). Several studies have focussed on crack path selection (He and Hutchinson 1989a, b), whereby conditions were identified such that a brittle mode I crack in the surface layer penetrates the underlying substrate, or kinks along the layer interface and generates delamination. Various aspects of crack deflection at the interface between the coating and the substrate have been investigated, including the formation of multi-fissures and debonding (León Baldelli et al. 2013), the origin of crack spacing (Lazarus and Pauchard 2011;Lazarus 2017;Ma et al. 2019), the contribution of residual stresses in relation to bending (Charalambides et al. 1990;Delette et al. 2012;Forschelen et al. 2016), and the influence of a thermal misfit between the coating and substrate (Lee et al. 2006). When delamination is absent, brittle mode I cracks may channel through both the coating and substrate; depending upon the stiffness mismatch, the depth of the channelling crack may be an order of magnitude larger than the coating thickness (Thouless et al. 2011). An extensive overview of studies on crack channelling, interfacial delamination and substrate damage is provided in the reference works Hutchinson and Suo (1992), Freund and Suresh (2004). of cracks channelling through the paint layer; c channelling crack that has kinked at the interface between the substrate and the paint layer, thereby introducing two opposite delaminations; d unlimited delamination that has led to flaking of the paint material. Pictures courtesy of Matteo Rossi Doria Inspired by the failure mechanisms observed in historical paintings, see Fig. 1, the present study addresses the behaviour of a brittle crack channelling in a bilayer under a moisture (or temperature) gradient in the thickness direction. The moisture gradient originates from a jump in humidity across the thickness of the system, with diffusion in steady-state, such that time transients do not play a role. The mode I channelling crack may deflect at the coating-substrate interface into two opposite delaminations of equal length, leading to the doubly deflected crack as sketched in Fig. 2.
Since the delamination length can vary, in principle, over the full range from zero to infinity, this configu-ration is representative of all three crack channelling scenarios as illustrated in Fig. 1b-d for a historical painting. The coating and substrate are characterized by two isotropic, linear elastic solids of dissimilar hygroscopic and mechanical properties. The assumption of elastic isotropy is acceptable for most paints; the substrate, however, is anisotropic in the case of wood or canvas, and this will lead to some inaccuracy in the predicted fracture response as computed herein. However, as demonstrated in a recent study on the fracture of both historic and new oak wood samples (Luimes et al. 2018b), the error is minor if the crack face normal aligns with a specific material direction, and the assumed value of Young's modulus for the isotropic elastic solid equals the anisotropic stiffness in the direction of the crack face normal.
Crack channelling is investigated by considering the three distinct fracture scenarios as summarized in Fig.  3a-c. These scenarios are: (i) channelling of a mode I crack in the coating with delamination absent (mechanism 1), (ii) channelling of a doubly deflected crack with finite delamination length (mechanism 2), and (iii) channelling of a mode I crack with unlimited delamination in all directions (mechanism 3). It is emphasised that these three mechanisms are representative of the cracking scenarios that are regularly observed in historical paintings, recall Fig. 1b The 3D failure mechanisms sketched in Fig. 3 are analysed by making use of the results from finite element calculations of steady-state crack channelling and plane-strain delamination. The critical remote stress for steady-state crack channelling may be determined from a plane-strain elasticity solution of a doubly deflected crack, whereby the difference in strain energy upstream and downstream of the channelling crack front is equated to the work of fracture of the doubly deflected crack (Hutchinson and Suo 1992;Beuth 1992;Suiker and Fleck 2004). The delamination toughness is taken to be constant and independent of the mode-mix; this assumption is acceptable, as the modemix attains a steady-state value at a relatively short delamination length for all configurations investigated. This modelling strategy for the analysis of 3D crack channelling mechanisms is similar to that used in Fleck and Zhao (2000) for microbuckle tunnelling in fibre composites and in Suiker and Fleck (2004), Cox and Marshall (1996), Suiker and Fleck (2006) for crack tunnelling in layered solids. The numerical results are used in the construction of failure mechanism maps, which illustrate the sensitivity of the active fracture mechanism and corresponding critical channelling stress to the ratio of layer interface to coating toughness, and to the mismatches in stiffness and in coefficient of hygral expansion. In the analysis of the results, the moisture content profile driving the fracture process is idealised by the superposition of a uniform contribution and a contribution that varies linearly across the thickness of the bilayer. In accordance with this decomposition, the remote stresses generated in the coating and substrate may be decomposed into constant and linear parts. These two stress contributions can be conveniently combined to construct the overall critical remote stress for crack channelling. Results are presented in terms of the hygral boundary conditions, but, from the analogy between hygral diffusion and thermal conduction, they can be immediately applied to channelling cracks under thermal boundary conditions.
The study considers a coating-substrate system containing a single crack. In historical paintings, however, multiple cracks can develop, with the typical crack spacing ranging between 10 and 50 times the coating thickness (Giorgiutti-Dauphiné and Pauchard 2016). It has been verified from the results of the numerical simulations that, for almost all the configurations analysed (with the exception of systems with a thick but compliant substrate), the stress field in the coating approaches the value of the remote stress at a distance from the delamination tip of less than six times the coating thickness. Since the crack spacing in historical painting typically exceeds the coating thickness by an order of magnitude, the cracks may be treated as isolated defects. Hence, the results of this study can also be applied to systems containing multiple parallel cracks. This conclusion is consistent with that reported in Thouless (1990). Although the present study has been motivated from observations of crack patterns in historical paintings, the modelling approach applies to the general case in which a surface layer on a substrate experiences channelling cracks under moisture and/or temperature gradients.
The paper is organized as follows. The problem is defined in Sect. 2. In Sect. 3 the governing equations for steady-state crack channelling and plane-strain delamination are presented. Section 4 gives the numerical results for a bilayer system composed of layers of equal thickness. The influence of the relative substrate thickness upon the failure response is then studied in Sect. 5. The main conclusions of the study are summarised in Sect. 6.

Definition of the problem
The three channelling mechanisms of Fig. 3 shall be predicted by considering the plane-strain crack as detailed in Fig. 4. Consider a bilayer system composed of a coating of thickness h 1 bonded to a substrate of thickness h 2 , with h = h 1 + h 2 the total thickness of the bilayer. Assume that a crack-like flaw develops across the coating thickness and subsequently deflects at the interface with the substrate into two opposite delaminations, each of length , see Fig. 4. The location of the crack is described by means of a Cartesian coordinate system (x, y), with the x-axis located along the interface between the coating and the substrate, and the y-axis coinciding with the vertical symmetry axis at the centre of the system. The coating and substrate are isotropic elastic solids, of Young's modulus E i , Poisson's ratio ν i , hygral expansion coefficient β i , and moisture diffusion coefficient D i , with i = 1 designating the coating and i = 2 denoting the substrate. There is choice in the assumed boundary condition on the bottom of the substrate. Here, roller supports are imposed to allow for horizontal displacements and to prevent vertical displacements (and consequently impose zero curvature on the system). The bilayer is subjected to plane-strain conditions, and experiences a steady-state hygroscopic loading as defined by a moisture content m that varies as a function of vertical coordinate y in the thickness direction.
The moisture content variation depends on the moisture diffusion coefficients of the coating, D 1 , and substrate, D 2 , and on the changes in moisture content at the top, m 1 = m 1 − m 0 , and at the bottom bottom, m 2 = m 2 − m 0 , where m 0 is a reference value of Fig. 4 Plane-strain cracking in a bilayer system. Doubly deflected crack at the interface between the coating and the substrate as a result of the applied moisture content profile m(y).
Bending and vertical deformations of the system are prevented through the fixed and roller supports at the bottom of the substrate. The remote stress profile depends on the system geometry, the moisture content profile, and the mismatch in hygromechanical properties of the individual layers moisture content. For a steady-state diffusion process, the moisture content m(y) across the thickness is bilinear, 1 and each linear branch can be split into a uniform part (with magnitude at y = 0) and a part that scales linearly with y as given by The in-plane, axial stresses σ 1 and σ 2 that are generated in the coating and substrate by the moisture content profile follow from the constitutive equations as are the plane-strain elastic modulus and coefficient of hygral expansion of layer i ∈ {1, 2}, and ε i (x, y) and β i m(y) are the total strain and hygroscopic strain, respectively.
Closed-form expressions for the remote, moistureinduced layer stresses that drive fracture can be derived 1 The bi-linear moisture profile given by expression (1) is obtained by solving the steady-state moisture diffusion equations for the coating and the substrate. The boundary conditions are imposed such that at the top and bottom surfaces of the system the moisture content variation is equal to m 1 and m 2 , respectively. Further, at the interface between the two layers, continuity conditions are imposed for the moisture content variation and moisture flux. from force equilibrium on an upstream cross-section remote from the crack. The total force on the crosssection vanishes. Upon expressing the upstream remote stresses in the coating and substrate as σ 1 (x → ∞, y) = σ ∞ 1 (y) and σ 2 (x → ∞, y) = σ ∞ 2 (y), respectively, this leads to Now substitute the moisture content profile (1) into the constitutive relations (2) 1 and (2) 2 , and insert the result into (3). Also, enforce the compatibility statement that the total axial strain is uniform across the upstream cross-section: Consequently, the remote stress distributions σ ∞ 1 (y) and σ ∞ 2 (y) become wherē Bilayer system subjected to moisture content variation. Remote stress componentsσ 1 and 1 presented by Eqs. (6) 1 and (6) 2 , as a function of the relative thickness h 2 / h 1 , for a selection of stiffness mismatchesĒ 2 /Ē 1 . The stress componentσ 1 is associated with a uniform moisture content profile equal to the moisture content at the coating-substrate interface, and the stress component 1 is associated with a linear moisture content profile that vanishes at the coating-substrate interface. The layers have the same diffusion coefficient, D 2 /D 1 = 1. The mismatch in the coefficient of hygroexpansion is aβ 2 /β 1 = 0.1 and b β 2 /β 1 = 10 Note from Eq. (5) that, in each layer i, the remote stress σ ∞ i is characterized by the sum of constant and linear stress contributions across the thickness. The constant stress contributionσ i is defined by Eqs. (6) 1 and (6) 4 , and is due to the first term on the right end side of Eq. (1): it is the contribution to σ ∞ i from a uniform moisture content m = m(0) = ( m 1 D 1 h 2 + m 2 D 2 h 1 )/(D 1 h 2 + D 2 h 1 ) that equals the moisture content value at the coatingsubstrate interface. The linear stress is quantified by the stress measure i as defined by Eqs. (6) 2 and (6) 5 , and is due to the linear moisture profile as stated by the second term on the right end side of Eq.(1). It is emphasised that this linear moisture content variation vanishes at the coating-substrate interface.
The coating stress contributionsσ 1 and 1 as given by Eqs. (6) 1 and (6) 2 are depicted in Fig. 5 as a function of the relative thickness h 2 / h 1 , for the case of vanishing moisture content at the bottom face of the bilayer, m 2 = 0. It is further assumed that the layers have equal diffusion coefficients, D 2 /D 1 = 1. The stresses in the coating are given in dimensionless form for selected stiffness mismatchesĒ 2 /Ē 1 = [0.1, 1, 10], and for mismatches in hygroexpansion coefficient of β 2 /β 1 = 0.1, see Fig. 5a, and ofβ 2 /β 1 = 10, see Fig. 5b. It is clear from Fig. 5a that, for a low coefficient of hygroexpansion of the substrateβ 2 /β 1 = 0.1, the stressesσ 1 and 1 asymptote under increasing substrate thickness to a compressive stress that is independent of stiffness mismatchĒ 2 /Ē 1 . For a coating and substrate of equal thickness, h 2 / h 1 = 1, a coating of moderate to high stiffness, 0.1 ≤Ē 2 /Ē 1 ≤ 1, experiences tensile stresses 1 that may induce cracking. Conversely, the constant stress contributionσ 1 is compressive, irrespective of the value of the stiffness mis-matchĒ 2 /Ē 1 and the relative coating thickness h 2 / h 1 . Figure 5b illustrates that for a relatively high coefficient of hygroexpansion of the substrateβ 2 /β 1 = 10, the stressσ 1 unconditionally lies in the tensile regime, while 1 lies in the compressive regime. Both stress parameters increase in magnitude with increasing substrate thickness, eventually approaching a limit value. It is emphasised that the above discussion relates to the case of hygroscopic swelling, with the imposed change in moisture content being positive, m 1 > 0. In the case of hygroscopic shrinkage such that the moisture content variation is negative, m 1 < 0, the stresses along the vertical axes of Fig. 5a, b change sign, and consequently "tension" and "compression" in the above discussion are interchanged. It is concluded that the constant and linear stresses in the coating may either lay in the compressive or tensile regime, depending on the moisture profile, stiffness mismatch, hygroexpansion coefficient mismatch and substrate-to-coating thickness ratio.
Consider first the case of a uniform moisture content profile that generates a uniform stress distribution within the coating, such that 1 = 0. Mode I fracture in the coating can only occur if the constant remote stressσ 1 in the coating is tensile, i.e., Referring to the remote stress definition (6) 1 , under hygroscopic swelling, m = ( m 1 D 1 h 2 + m 2 D 2 h 1 ) /(D 1 h 2 + D 2 h 1 ) > 0, this condition is met for a hygroexpansion coefficient mismatchβ 2 /β 1 > 1, independent of the value of the stiffness mismatch and the relative substrate thickness. Conversely, for the case of hygroscopic shrinkage, m < 0, the coating is in tension for a mismatch in hygroexpansion coefficient β 2 /β 1 < 1. Next consider the case of a linear moisture content variation that vanishes at the layer interface, along with σ 1 = 0. For simplicity, it is assumed that the linear stress distribution in the coating that follows from Eq. (5) 1 as 1 (1 + y D 2 /(η 1 (D 1 h 2 + D 2 h 1 ))), is tensile across the entire coating thickness. This condition is satisfied for the choice: Note from the remote stress definition (6) 2 that requirement (8) 1 is met if the function η 1 and the moisture content difference m 2 − m 1 have the same sign. For the choice m 1 − m 2 > 0, this implies that η 1 < 0, and, upon referring to the definition (6) 3 , this is equivalent toĒ 2 D 1β2 h 2 2 /(Ē 1 D 2β1 h 2 1 ) < 1. Conversely, for m 1 − m 2 < 0, the remote stress 1 is tensile if η 1 > 0, leading toĒ 2 D 1β2 h 2 2 /(Ē 1 D 2β1 h 2 1 ) > 1. Additionally, condition (8) 2 holds when it is satisfied at its most critical location, which is at the top surface of the coating, y = h 1 . Note that (8) 2 depends solely on the hygro-mechanical and geometrical properties of the bilayer system, and is independent of the moisture content profile imposed. The material parameters and geometrical properties of the bilayer systems analysed in Sects. 4 and 5 are selected such that condition (8) 2 is satisfied at y = h 1 , when condition (8) 1 holds. Specifically, the selected elastic mismatches between the coating and substrate areĒ 2 /Ē 1 = [0.1, 0.3, 1, 3, 10], with their Poisson's ratios being taken equal, ν 1 = ν 2 = 0.3. Further, the diffusion coefficient is chosen to be equal for the two layers, D 2 /D 1 = 1. If the thicknesses of the coating and substrate are in the same range, i.e., h 2 / h 1 < 10, for the range of stiffness ratios mentioned above, condition (8) 2 is unconditionally satisfied at y = h 1 for hygroexpansion coefficient mismatchesβ 2 /β 1 ≥ 10. Hence, as an illustrative case, in the numerical analyses the mismatch in hygroexpansion coefficients will be taken asβ 2 /β 1 = 10.
Note finally that the above model may be also used for thermo-mechanical loading, whereby Eq. (1) refers to a temperature variation and the coefficients D 1 and D 2 are substituted by the thermal conductivities of the coating and substrate, respectively. Further, the coefficientsβ 1 andβ 2 in the constitutive equations (2) 1 and (2) 2 are replaced by the thermal expansion coefficients of the coating and substrate, respectively.

Steady-state crack channelling and plane-strain delamination
The governing equations for steady-state crack channelling and plane-strain delamination in a coatingsubstrate system under a moisture gradient can be derived from a framework similar to that employed for microbuckle tunnelling in fibre composites in Fleck and Zhao (2000) and that developed for crack tunnelling in layered solids (Cox and Marshall 1996;Fleck 2004, 2006). Accordingly, it is assumed that the channelling crack has nucleated from an initial flaw in the coating and develops as a result of a moistureinduced tensile remote stress σ ∞ 1 (y) given by Eq. (5) 1 . During steady-state channelling, the channelling front has a constant shape, whereby the energy release rate does not depend on the channelling length in the outof-plane direction of the crack. Hence, the energy released per unit advance of crack channelling can be calculated as the difference in the elastic strain energy W upstream and downstream of the channelling front (Hutchinson and Suo 1992;Beuth 1992;Fleck and Zhao 2000;Suiker and Fleck 2004), which equals the difference between the strain energy stored in the uncracked plane-strain solid and in the cracked plane-strain solid. For the channelling crack shown in Fig. 4, this energy difference equals where δ(y) is the opening displacement across the crack faces within the coating. For convenience, this crack opening displacement δ(y) is decomposed as where the contributionδσ (y) represents the crack opening displacement when the remote stress in the coating is constant and equal toσ 1 , whileδ (y) is the crack opening displacement associated with a linear remote stress, characterized by the stress measure 1 , see Eq. (5) 1 . Insert Eqs.
(2) 1 and (10) into Eq. (9), and apply Betti's reciprocity principle; then, the elastic strain energy can be decomposed into three contributions, in a similar way to that done in Vlassak (2003). The first contribution Wσ relates to the stressσ 1 that is constant across the coating thickness, the second contribution W is associated with a linear stress profile described by the stress measure 1 , and the third contribution Wσ results from the coupling between the constant stress and linear stress distributions, 2 i.e., where with the corresponding average displacement parameters as 2 Note that the elastic strain energy Wσ associated with the coupling between the constant and linear stress distributions has been defined here by relations (12) 3 and (13) 3 . By referring to Betti's reciprocity principle, it can be equivalently written as 3.1 Steady-state crack channelling with a finite delamination length In accordance with the decomposition of the strain energy as stated by Eq. (11), steady-state channelling of the doubly deflected crack can be analysed by separately considering the constant and linear stress distributions in the coating. Hence, the contribution Wσ associated with the constant stressσ 1 alone, from a uniform moisture content variation, is considered first. Based on dimensional considerations, the corresponding average displacementδσ is expressed in the form where the dimensionless function f depends on the aspect ratio / h 1 of the doubly deflected crack, the stiffness mismatchĒ 2 /Ē 1 , the mismatch in hygroexpansion coefficientsβ 2 /β 1 and the ratio between the moisture diffusion coefficients D 2 /D 1 . Upon substituting Eq. (14) into Eq. (12) 1 , the energy drop per unit crack length becomes Under a uniform moisture content, the moisture content profile (1) 1 becomes m(y) = m 1 = m 2 , and the remote stress σ ∞ 1 coincides with the constant stress contributionσ 1 , see Eq. (5) 1 . The corresponding average energy release rate per unit advance of a channelling crack G c,σ is directly related to the elastic energy drop via the expression (Hutchinson and Suo 1992;Beuth 1992;Fleck and Zhao 2000;Suiker and Fleck 2004) The energy drop Wσ reflects the energy necessary to form a mode I crack of length h 1 in the coating and two delamination cracks of length . Write the mode I toughness of the coating as I and the delamination toughness as d ; then, the energy drop can be written as In general, the delamination toughness d is a function of the mode-mix of the interfacial crack, as will be discussed in more detail in Sect. 3.2. The critical remote stress for steady state channellingσ 1 =σ c can be computed by combining Eqs. (17) and (15), leading tō Now consider the second term W of the energy drop given by Eq. (11), which is related to the linear stress measure 1 . Analogous to Eq. (14), the average dis-placementD associated with 1 can be written as Here, the non-dimensional function g characterizes the doubly deflected crack, and is a dimensionless function of the aspect ratio / h 1 , the stiffness mismatchĒ 2 /Ē 1 , the mismatch in hygroexpansion coefficientβ 2 /β 1 and the ratio between the moisture diffusion coefficients D 2 /D 1 . Substituting Eq. (19) into Eq. (12) 2 , the energy drop W can be expressed as For a linear profile in moisture content, with a vanishing value at the coating-substrate interface, Eq.
(1) 1 reduces to m(y) = y D 2 ( m 1 − m 2 )/(D 1 h 2 +D 2 h 1 ), and the remote stress in Eq. (5) 1 reduces to the linear distribution, σ ∞ 1 (y) = 1 (1+ y D 2 /(η 1 (D 1 h 2 + D 2 h 1 ))). The corresponding average energy release rate for unit advance of a channelling crack G c, is expressed as The critical remote stress for steady-state channelling, 1 = c , is obtained by combining Eqs. (21) and (20), and by assuming that the energy drop W is related to the mode I toughness I and the delamination toughness d in a similar manner to that given in (17), such that Consider finally the energy drop contribution Wσ in Eq. (11) that is associated with the coupling between the constant and linear parts of the stress field. From dimensional considerations, the corresponding displacement Dσ defined by relation (13) 3 can be written in the form Here, k is a dimensionless function that depends upon the aspect ratio of the crack / h 1 , the stiffness mis-matchĒ 2 /Ē 1 , the mismatch in hygroexpansion coefficientβ 2 /β 1 and the ratio between the moisture diffusion coefficients D 2 /D 1 . Substituting Eq. (23) into Eq. (12) 3 , the energy drop Wσ follows as The remote channelling stress σ ∞ c (y) associated with the general moisture content profile as given by (1) can be computed from the total energy drop W in Eq. (11) for the doubly deflected crack, such that where G c is the average energy release rate for crack channelling. Inserting the energy components Wσ , W and Wσ as given by Eqs. (15), (20) and (24), respectively, into Eq. (11), and combining the result with Eq. (25), leads to the general result The above relation gives a failure locus inσ 1 versus 1 space for any assumed delamination length . In order to determine the operative mechanism 1, 2 or 3 of Fig. 3, assume thatσ 1 is specified, and solve (26) for 1 as a function of . If 1 is smallest for = 0 then mechanism 1 operates, if 1 is smallest for a finite but nonvanishing value of then mechanism 2 operates and if In particular, to calculate 1 from relation (26), first writeσ c as the uniform remote stressσ 1 for steady-state channelling such that 1 = 0, recall (18). Next, assume thatσ 1 is prescribed asσ 1 = ασ c , where the weighting factor α lies within the range 0 ≤ α ≤ 1. Since the sign ofσ c is positive, the range chosen for α ensures that the uniform remote stressσ 1 corresponds to a tensile stress, consistent with condition (7) that enables mode I crack development in the coating. The quadratic equation (26) is then solved for 1 , giving a positive and a negative solution for the corresponding remote stress measure for steady-state crack channelling. In accordance with condition (8) 1 , the positive solution is assumed for the critical channelling stress, 1 = α c . Then, upon substitutingσ 1 = ασ c and 1 = α c into Eqs. (6) 1 and (6) 2 , the critical moisture content values m 1 and m 2 at the top and bottom surfaces of the system are determined, and the corresponding moisture profile for steady-state crack channelling via Eq. (1), for the assumed value of . Furthermore, the selected value for ασ c and the associated value of α c can be substituted into Eq. (5) 1 in order to obtain the remote stress profile for crack channelling, σ ∞,α c (y). Note that the limit α = 1 results in the uniform remote channelling stress, σ ∞,α c (y) =σ c , withσ c given by Eq. (18), and the limit α = 0 leads to the linear remote channelling stress, σ ∞,α c (y) = c (1 + y D 2 /(η 1 (D 1 h 2 + D 2 h 1 ))), with c following from (22).

Plane-strain crack with a finite delamination length
In addition to analysing the problem of crack channelling with a finite delamination length, it is of interest to consider delamination growth from the tip of a plane-strain crack that exists across the thickness of the coating, recall Fig. 4. The plane-strain crack problem is closely related to the channelling problem, and helps in the identification of mechanisms 1, 2 or 3 for crack channelling. The details are as follows.
First, characterize the plane-strain delamination crack between two elastic isotropic but dissimilar materials, by introducing the Dundurs' parameters A and B that quantify the elastic mismatch between the two solids (Dundurs 1969). They are defined by where G j = E j /(2(1 + ν j )) is the shear modulus of layer j. In general, the singular stress field at the tip of an interfacial crack is characterized by a complex stress intensity factor K = K 1 + i K 2 , with i = √ −1 and K 1 , K 2 the real and the imaginary parts of the stress intensity factor, respectively. The singular normal σ yy and shear σ xy stress components at a distance r from the tip in the asymptotic limit read (Hutchinson et al. 1987;Rice 1988) where r = e i lnr = cos( lnr ) + i sin( lnr ), and the oscillatory index is defined as in terms of the second Dundurs' parameter B as introduced in Eq. (27) 2 . The phase angle is directly related to the ratio of shear stress to normal stress on the crack plane immediately ahead of the crack tip, as follows. This ratio is dependent on the distance r from the tip, and consequently it is necessary to evaluate the modemix at an arbitrary, but specified, reference lengthl, ahead of the crack tip. As a consequence of the oscillatory stress behaviour at the tip of an interfacial crack between two dissimilar solids, the classical definition of the mode-mix is extended to the form (Rice 1988) The choice of the reference lengthl has only a minor effect on the value of the angle for commonly encountered values of the oscillatory index, | | 1 (Rice 1988).
The plane-strain energy release rate per unit advance of delamination along a bi-material interface is given by (Hutchinson and Suo 1992) where E * = 2(Ē −1 1 +Ē −1 2 ) −1 . Delamination will occur when the energy release rate as defined in (31) equals the delamination toughness at the appropriate modemix d ( (l)), i.e., The plane-strain energy release rate for delamination can be further related to the energy drop W , given by relation (9). For the doubly deflected crack sketched in Fig. 4, the energy release rate for each delamination can be calculated as where the factor of 2 is due to the fact that there are two delamination tips. Similar to the decomposition of the energy drop W as given by Eq. (11), the energy release rate may be decomposed into Upon inserting expression (15) into (33), the first term in the right-hand side of Eq. (34) becomes with f representing the partial derivative f = ∂ f /∂ . Hence, if the coating-substrate system is subjected to a uniform moisture content with m(y) = m 1 = m 2 , then it follows from Eq. (34) that G d = G d,σ , so that the critical remote stress for plane-strain delamination,σ 1 =σ d , is obtained by combining Eq. (32) with Eq. (35), to givē In an analogous fashion, the energy release rate per unit advance of delamination G d, , which is related to the stress measure 1 , is derived from expression (20), such that where g represents the partial derivative g = ∂g/∂ . Accordingly, in the case of a linear moisture content profile of vanishing value at the coating-substrate interface, m(y) = y D 2 ( m 1 − m 2 )/(D 1 h 2 + D 2 h 1 ), the energy release rate for plane-strain delamination follows from relation (34) as G d = G d, . Now combine Eq. (32) with Eq. (37); then, the critical remote stress for delamination, 1 = d , follows as Consider finally the energy drop Wσ associated with the coupling between the constant and the linear parts of the stress field, as given by Eq. (24). The corresponding energy release rate G d,σ for plane-strain delamination is computed as in terms of the partial derivative k = ∂k/∂ . Substituting expressions (35), (37) and (39) into relation (34) then leads to the energy release rate for plane-strain delamination under the arbitrarily linear moisture content profile given by relation (1). Now equate this value of energy release rate to the delamination toughness, upon recalling Eq. (32), to obtain The remote delamination stress under the moisture content profile (1) 1 is finally found in a similar way to that demonstrated at the end of Sect. 3.1 for the channelling stress. Relation (40) represents a failure locus inσ 1 versus 1 space. It is assumed thatσ 1 is specified, and Eq. (40) is solved with respect to 1 . In particular, first takeσ d as the uniform remote stressσ 1 for plane-strain delamination, with 1 = 0, in accordance with (36). Next, the value ofσ 1 is prescribed asσ 1 = ασ d , with α a weighting factor in the range 0 ≤ α ≤ 1. The positive sign ofσ d , along with the selected range for α, guarantees that the uniform remote stressσ 1 is tensile, in accordance with condition (7). Further, Eq. (40) is solved for 1 , providing a positive and a negative solution for the corresponding remote stress for planestrain delamination. The positive solution is taken as the critical delamination stress, 1 = α d , consistent with requirement (8) 1 . Next, the critical moisture content values m 1 and m 2 at the top and bottom surfaces of the system can be determined, by insertinḡ σ 1 = ασ d and 1 = α d into Eqs. (6) 1 and (6) 2 . Using the values of m 1 and m 2 , the corresponding moisture profile for plane-strain delamination can be calculated through Eq. (1). Finally, by substituting into Eq. (5) 1 the assumed value for ασ d and the corresponding value of α d , the remote stress distribution for plane-strain delamination σ ∞,α d (y) can be determined. Analogous to the mechanism of crack channelling, for α = 1 the remote delamination stress becomes uniform, σ ∞,α d (y) =σ d , withσ d in accordance with (36), whereas for α = 0 it reduces to the linear stress profile σ ∞,α d (y) = d (1 + y D 2 /(η 1 (D 1 h 2 + D 2 h 1 ))), with d given by Eq. (38).

Unlimited delamination
When the delamination length exceeds the thickness of the coating, / h 1 > 1, the energy release rate for delamination typically approaches an asymptotic steady-state value, G d,ss = G d ( / h 1 → ∞), see also Hutchinson and Suo (1992), Freund and Suresh (2004), Suiker and Fleck (2004), Forschelen et al. (2016). This value can be derived analytically by computing the difference in the elastic strain energy downstream and upstream of the delamination tip. For this purpose, the geometry depicted in Fig. 4 is considered, whereby the system is subjected to the arbitrarily linear moisture content profile (1). Recall that the remote stresses σ ∞ 1 (y) and σ ∞ 2 (y), as specified by (5) 1 and (5) 2 , respectively, have been computed via relations (2) 1 and (2) 2 in combination with force equilibrium of the intact, upstream cross-section via Eq. (3). In a similar fashion, force equilibrium at the fractured, downstream cross-section x = 0 may be expressed as where, in correspondence with Eq.
(2) 2 , the stress σ 0 2 (y) = σ 2 (x = 0, y) is evaluated at x = 0. Inserting the moisture content profile (1) 2 into the constitutive relation (2) 2 , and substituting the result into relation (41) allows one to solve for the strain ε 2 (x = 0, y). Substituting this result back into Eq. (2) 2 yields for the stress σ 0 2 : Note that σ 0 2 is described solely by the stress contribution associated with the linear moisture content variation in the substrate, and the contribution from the uniform moisture content variation vanishes, i.e.,σ 0 2 = 0. The steady-state energy release rate at an individual delamination tip can be obtained from the difference between the energies upstream and downstream of the delamination tip, i.e., with ε e,∞ i (y) the remote, upstream elastic strain of layer i and ε e,0 2 (y) the downstream elastic strain in the substrate. Now insert expressions (5) 1 , (5) 2 and (42) for the upstream and downstream stresses into (43), and invoke the constitutive relation σ i =Ē i ε e i with i ∈ {1, 2}, to obtain G d,ss = G ss,σ + G ss, + G ss,σ , where G ss,σ = 1 2σ Here, the sub-index "d" has been dropped in the above energy release rate, as the steady-state value is representative of both plane-strain delamination and crack channelling, i.e., G ss = G d ( / h 1 → ∞) = G c ( / h 1 → ∞). Specifically, the equality between the crack channelling and delamination energy release rates occurs at any point of the G c versus / h 1 curve for which ∂G c /∂ = 0 (Beuth 1992;Fleck and Zhao 2000;Suiker and Fleck 2004). It will be demonstrated in Sect. 4 that the finite element results confirm this condition when / h 1 → ∞.
During unlimited delamination the energy release rate equals the delamination toughness d at the appropriate mode-mix, see Eq. (32). Upon substituting expression (45) 1 into relation (32) allows to solve for the steady-state stress valueσ ss =σ d ( / h 1 → ∞) =σ c ( / h 1 → ∞) corresponding to a uniform moisture content profile m(y) (= m 1 = m 2 ), with 1 = 0. Similarly, substituting expression (45) 2 into (32) leads to the steady-state stress Hence, the expressions for the stressesσ ss and ss in the coating arē with the parameters η 1 and η 2 in Eq. (45) 2 provided by Eqs. (6) 3 and (6) 6 , respectively. Note that expression (46) 1 for the critical steady-state stress under a uniform moisture content distribution is independent of both the mismatch in coefficient of hygral expansion,β 2 /β 1 and the mismatch in diffusion coefficient, D 2 /D 1 . The remote critical steady-state stress value corresponding to an arbitrary moisture profile can be finally obtained, as follows. For a given value of α, the contributionσ 1 isσ 1 = ασ ss , and upon equating the energy Fig. 6 Steady-state critical stress σ ∞,α ss (0) in the coating versus stiffness mismatchĒ 2 /Ē 1 for different weighting factors α, see Eq. (47). The stress is computed for a bilayer system with layers of equal thickness, h 2 / h 1 = 1, equal diffusion coefficient, D 2 /D 1 = 1, and for a mismatch in coefficient of hygral expansionβ 2 /β 1 = 10 release rate to the delamination toughness d in (44), the steady-state stress 1 = α ss is calculated. The steady-state stress profile σ α ss (y) follows immediately as σ ∞,α ss (y) = ασ ss + α ss 1 + y D 2 η 1 (D 1 h 2 + D 2 h 1 ) .
(47) Figure 6 depicts the steady-state stress, Eq. (47), for arbitrary moisture profiles, evaluated at the interface between the coating and substrate with y = 0, i.e. σ ∞,α ss (0), as a function of the stiffness mismatch, E 2 /Ē 1 , for a selection of weighting factors, α = [0, 0.5, 0.7, 0.9, 1]. The stress is computed for a system characterized by equal layer thicknesses, h 1 = h 2 , by equal diffusion coefficients, D 1 = D 2 , and by a mismatch in the coefficients of hygral expansion β 2 /β 1 = 10, in order to satisfy condition (8) 2 . The stress is obtained from expression (47), for which the limit values α = 0 and α = 1 respectively reduce to the linear σ ∞,α ss (0) = ss and uniform σ ∞,α ss (0) =σ ss steady-state stress measures given in Eq. (46). It can be seen that for all values of α the steady-state stress grows monotonically with an increasing stiffness ratio. The response envelope is indeed provided by the cases α = 0 and α = 1, from which it is concluded that the corresponding stress measuresσ ss and ss may respec-tively serve as lower and upper bounds in the practical design against plane-strain delamination.

Geometry and modelling features
Crack channelling and plane-strain delamination are studied first for a system in which the coating and substrate have the same thickness, i.e., h 1 / h 2 = 1. The stress state is analysed using the commercial finite element program ABAQUS Standard. 3 Due to symmetry, only half of the domain is considered. Roller and fixed supports impose the required symmetry and prevent rigid body motions, and vertical supports at the bottom edge prevent bending of the bilayer. The geometry is taken to be sufficiently long in order for boundary effects to be negligible, in accordance with a horizontal system length equal to 200h 1 . The fracture response has been analysed for a crack with a delamination length within the range / h 1 ∈ [0.015, 20]. For each delamination length a different finite element mesh was used. The finite element configurations contain 2500 to 3600 plane-strain 8-nodes iso-parametric elements, equipped with 3 × 3 Gauss quadrature. At the delamination tip, the mid-side nodes on the crack faces are moved to the 1/4 point nearest to the crack tip, in order to represent the square root singularity of the stress field. Further, for each tip element, three neighbouring nodes are collapsed to the crack tip. As already discussed in Sect. 2, the elastic mismatch between the coating and substrate is taken to bē E 2 /Ē 1 = [0.1, 0.3, 1, 3, 10], and the Poisson's ratio of the coating and substrate is ν 1 = ν 2 = 0.3. Consequently, the Dundurs' parameters as given by Eq. (27) take the values A = [0.81, 0.53, 0, −0.5, −0.82] and B = [0.23, 0.15, 0, −0.14, −0.23]. The hygroexpansion coefficient mismatch is chosen asβ 2 /β 1 = 10, and for all analyses performed the coating is fully loaded in tension. The diffusion coefficients of the layers are considered to be equal, i.e., D 2 /D 1 = 1. The system is subjected to the moisture profile as given by Eq. (1), whereby the moisture content at the top surface of the system equals m 1 and at the bottom surface is m 2 = 0. Although the coating-substrate sys-tem is subjected to a specific moisture content profile, the FEM results for steady-state crack channelling and plane-strain delamination can be generalized to a wide range of critical moisture content profiles, which are obtained by varying the value of the weighting factor α introduced in relation to the solution of Eqs. (26) and (40), see also the discussion below.
The FEM data required for the calculation of the critical remote stresses for steady-state crack channelling and plane-strain delamination are (i) the displacement profile δ(y) of the nodes at the mode I crack faces, (ii) the energy release rate per unit advance of delamination G d at the delamination tip, which for a crack in an elastic solid equals the path-independent J -integral, G d = J (Rice 1968), and (iii) the complex stress intensity factor at the delamination tip, K = K 1 + i K 2 . The real and imaginary parts of the stress intensity factor, K 1 and K 2 , are calculated via an embedded routine within ABAQUS using an interaction integral method: the mode I and mode II stress intensity factors K 1 and K 2 are extracted from the energy release rate G d for plane-strain delamination by combining the solution of the actual crack tip field with that of an auxiliary field (Matos et al. 1989). In accordance with the modelling framework presented in Sect. 3, the critical remote stresses for steady-state crack channelling and plane-strain delamination are characterized by constant (σ c andσ d ) and linear ( c and d ) stress contributions, which can be derived from two separate FEM simulations that respectively consider (i) only the constant part of a simulated moisture content profile, and (ii) the total moisture content profile. Considering a system with specific stiffness and toughness properties, first a FEM simulation is performed in which only a constant moisture content part is applied. For the moisture content profile and geometry considered, with D 2 /D 1 = 1, the value of the constant moisture content follows from Eq. (1) as m = m(0) = m 1 /2. The constant moisture content induces a constant remote stress in the coating, σ ∞ 1 (y) =σ 1 , and in the substrate, σ ∞ 2 (y) =σ 2 , see Eq. (5). The generated displacements δσ (y) of the mode I crack faces provide the average crack opening displacementδσ via relation (13) 1 , and its non-dimensional value f by means of Eq. (14). Subsequently, the drop in strain energy Wσ follows from Eq. (12) 1 , and this is substituted into (16) to give the energy release rate G c,σ for steady-state crack channelling. The function f further results in the constant remote stressσ c for crack channelling via Eq. (18). In addition, the J -integral leads to the energy release rate for plane-strain delamination, G d,σ = J , and thereby to the derivative f through relation (35). The constant remote stressσ d for plane-strain delamination follows from Eq. (36).
As a next step, the same coating-substrate system is subjected to the total moisture content profile, which, for the moisture content profile and geometry considered, with D 2 /D 1 = 1, follows from Eq. (1) as m(y) = m 1 /2+ m 1 y/ h. The mode I crack opening displacement δ(y) calculated for the total moisture content profile is combined in Eq. (10) with the displacement function δσ (y) for a constant moisture content to obtain the mode I crack opening displacement under the linear moisture content profile, δ (y) = δ(y) − δσ (y). Subsequently, the average crack opening displacement parametersD andDσ related to the linear and constant moisture contributions are calculated from Eqs. (13) 2 and (13) 3 , and the non-dimensional values g and k are thereby obtained via relations (19) and (23), respectively. The corresponding drops in elastic energy, W and Wσ , result from Eqs. (12) 2 and (12) 3 , respectively, and consequently the total drop in elastic energy W is determined via (11), and the energy release rate G c, per unit advance of steadystate crack channelling through (21). Subsequently, by choosing a value for α, the fraction of the constant channelling stress experienced by the coating becomes σ 1 = ασ c . This stress value is inserted into (26) to give the corresponding linear stress measure for crack channelling, 1 = α c , which via (5) 1 leads to the total stress for steady-state crack channelling, σ ∞,α c (y). Furthermore, the energy release rate for plane-strain delamination generated under the total moisture content profile is obtained from the J -integral as G d = J . From the corresponding values of the mode I and mode II stress intensity factors at the delamination tip, K 1 and K 2 , the stress intensity factors K 1, and K 2, due to a linear moisture content profile are calculated as where K 1,σ and K 2,σ are the delamination stress intensity factors for a constant moisture content profile. The values of K 1, and K 2, are inserted into Eq. (31) to give the delamination energy release rate G d, under a linear moisture content profile. In accordance with (34), the coupling term for the delamination energy release rate can now be calculated as G d,σ = G d − G d,σ − G d, , and the corresponding derivative k follows from Eq. (39). The value chosen for α specifies the fraction of the constant delamination stress present in the coating, σ 1 = ασ d , which, after substitution into (40), leads to the critical linear stress measure for plane-strain delamination, 1 = α d . The total stress for plane-strain delamination, σ ∞,α d (y), is obtained from Eq. (5) 1 . As pointed out above, the value of α specifies the relative contributions of the constant and linear remote stress measuresσ 1 and 1 to the activation of steadystate crack channelling and plane-strain delamination. For any given value of α, the critical moisture contents m 1 and m 2 for steady-state crack channelling and plane-strain delamination at the top and bottom surfaces of the coating-substrate system can be obtained from the corresponding constant and linear stress measures via relations (6) 1 and (6) 2 . By varying the factor α in the range 0 ≤ α ≤ 1, the two FEM simulation results above thus provide a family of critical linear moisture content profiles for both steady-state crack channelling and plane-strain delamination, with the limit value α = 1 reflecting the uniform moisture content distribution and α = 0 representing the linear moisture content distribution with a vanishing value at the coating-substrate interface.
4.2 Cracking due to a uniform moisture content distribution (α = 1) For the coating-substrate system described in Sect. 4.1, the failure mechanisms of steady-state crack channelling and plane-strain delamination are first examined for a uniform moisture content distribution, m(y) = m 1 = m 2 . Accordingly, the weighting factor α is set to unity, α = 1: the remote stresses for crack channelling and plane-strain delamination are uniform, such that σ ∞,α c (y) =σ c and σ ∞,α d (y) =σ d . All numerical results will be presented using appropriate dimensionless parameters. For the case of a uniform moisture content contribution, the obtained modemix and critical stresses turn out to be independent of both the mismatch in the coefficient of hygral expansion,β 2 /β 1 , and the mismatch in diffusion coefficient, D 2 /D 1 . The normalization selected for the energy release rates makes these results also independent of the hygroscopic coefficient mismatch,β 2 /β 1 , and diffusion coefficient mismatch, D 2 /D 1 .

Mode-mix and energy release rate
For the considered configuration, the mode-mix σ associated with a uniform moisture content is illustrated in Fig. 7a as a function of the relative delamination length / h 1 , for selected values of the stiffness mismatchĒ 2 /Ē 1 . In correspondence with Eq. (30), the phase angle of mode-mix is calculated as σ = atan(Im(Kσl i )/Re(Kσl i )), where the complex stress intensity factor Kσ is defined by Kσ = K 1,σ +i K 2,σ , in terms of the mode I and mode II stress intensity factors for interfacial delamination, K 1,σ and K 2,σ , respectively. The reference length at which the mode-mix is evaluated is set equal to the thickness of the coating,l = h 1 . For all stiffness mismatches considered, the mode-mix σ varies over a relatively small range, and a steady-state value is attained when the delamination length approaches the coating thickness, i.e., when / h 1 ≥ 0.8. The values of the phase angle at steadystate lies within the range 52 • to 57 • . The evolution of the mode-mix with delamination length is comparable to that computed in Suiker and Fleck (2004) for a doubly deflected crack generated under a remote tensile stress in the two outer layers of a layered solid, see Figure 16a in Suiker and Fleck (2004). Figure 7b, c depict the energy release rates for planestrain delamination G d,σ and steady-state crack channelling G c,σ , respectively, as a function of the delamination length / h 1 , generated under a constant moisture content. The energy release rate for plane-strain delamination is strongly sensitive to the delamination length for / h 1 < 1, whereby in the asymptotic limit / h 1 → 0 it becomes equal to zero for stiffness mis-matchesĒ 2 /Ē 1 < 1, and goes to infinity forĒ 2 /Ē 1 > 1. This behaviour has also been reported in other studies on interfacial delamination (He and Hutchinson 1989b;Ye et al. 1992;Lu 1996;Suiker and Fleck 2004). For E 2 /Ē 1 = 1 the energy release rate at zero delamination is finite, and is slightly less than the corresponding value for a doubly-deflected crack in the two outer layers of a multi-layered solid (Suiker and Fleck 2004). For l/ h 1 > 1 the energy release rate G d,σ for planestrain delamination has almost attained a steady state, and the values for the range of stiffness mismatches are in close agreement with the analytical steady-state energy release rate, G ss,σ = G d,σ (l/ h 1 → ∞), as given by Eq. (45) 1 .
The energy release rate G c,σ for crack channelling, shown in Fig. 7c, decreases monotonically with increasing delamination length / h 1 for all values of stiffness mismatchĒ 2 /Ē 1 . Compared to the energy release rate for plane-strain delamination, the steadystate energy release rate G ss,σ is reached at a substantially larger delamination length. Moreover, the steadystate value is reached at a smaller delamination length with increasing stiffness mismatchĒ 2 /Ē 1 . These two features have also been observed for the crack tunnelling mechanisms studied in Suiker and Fleck (2004).

Failure mechanisms
The crack channelling stressσ c and plane-strain delamination stressσ d associated with a uniform moisture content profile in the coating-substrate system are shown in dimensionless form in Fig. 8 as a function of delamination length / h 1 . Figures 8a-c contain the numerical results for stiffness mismatchesĒ 2 /Ē 1 = 0.1,Ē 2 /Ē 1 = 1 andĒ 2 /Ē 1 = 10, respectively. The channelling stressσ c has been computed using Eq. (18) and is indicated by the dark blue lines for a selection of toughness ratios d / I = [0.05, 0.1, 0.2, 0.5, 1, 10].
Although the delamination toughness d is, in principle, a function of the mode-mix σ , for the calculation of the channelling stressσ c a constant delamination toughness is assumed. This assumption is reasonable since the mode-mix typically has only a small variation and reaches a steady-state value at relatively short delamination length, / h 1 ≈ 0.8, see Fig. 7a. Since the stresses are normalized with respect to the delamination toughness d , a channelling stressσ c associated with a larger toughness ratio d / I corresponds to a smaller mode I toughness of the coating material. The curve for the delamination stressσ d , plotted by the light blue line, is obtained from Eq. (36), and may be interpreted as a crack-growth resistance curve (R-curve). The depicted analytical steady-state stress, σ d ( / h 1 → ∞) =σ c ( / h 1 → ∞) =σ ss , follows from Eq. (46) 1 , and adequately reflects the asymptotic limit to which theσ c -andσ d -curves converge under increasing delamination.
In accordance with the procedure presented in Cox and Marshall (1996), Fleck and Zhao (2000), Suiker and Fleck (2004), from the results in Fig. 8 the three main failure scenarios sketched in Fig. 3 can be distinguished, which depend on the toughness ratio d / I and the elastic mismatchĒ 2 /Ē 1 . Consider first the crack channelling and delamination stresses for the case of a stiffness mismatchĒ 2 /Ē 1 = 0.1 as depicted in Fig. 8a. The delamination stress has an initial, rising branch (R-curve behaviour) for relatively short delamination lengths, / h 1 ≤ 0.45; subsequently, after having reached its maximum value, it drops slightly under growing delamination in order to reach the steadystate stressσ ss . Since stable plane-strain delamination is characterized by an increasing stressσ d versus delamination length / h 1 , channelling of a doubly deflected crack with stable delamination corresponds to a crack channelling curve intersecting with the rising branch of the delamination curve, which occurs for d / I ≥ 0.58, up to some maximum value that is discussed below. At the intersection between the actualσ c ( / h 1 ) curve and theσ d ( / h 1 ) curve, the corresponding delamination length / h 1 can be read off along the horizontal axis of Fig. 8a. Further, the corresponding value of the channelling stress is represented by a minimum (extremum), which becomes clear when applying the condition ∂σ c /∂ = 0 to expression (18) for the channelling stress, leading toσ c,min =σ d with σ d given by Eq. (36). The above fracture scenario is denoted in Fig. 3 as mechanism 2. (a) (b) (c) Fig. 8 Crack channelling and plane-strain delamination in a bilayer system with layers of equal thickness (h 2 / h 1 = 1). Remote stressσ 1 versus delamination length / h 1 for a constant moisture content profile, considering a stiffness mismatch of aĒ 2 /Ē 1 = 0.1, bĒ 2 /Ē 1 = 1 and cĒ 2 /Ē 1 = 10. The dark blue lines represent the crack channelling stressσ 1 =σ c for selected toughness ratios d / I . The light blue line reflects the plane-strain delamination stressσ 1 =σ d . The results are independent of the mismatch in coefficient of hygral expansion, β 2 /β 1 , and the mismatch in diffusion coefficient, D 2 /D 1 . The depicted steady-state valueσ ss is computed with Eq. (46) 1 In contrast, for toughness ratios d / I < 0.58, the channelling stressσ c decreases monotonically with increasing delamination / h 1 , as a result of which the critical, minimum channelling stress coincides with the steady-state value at infinite delamination,σ c,min = σ c ( / h → ∞) =σ ss . Hence, crack channelling is characterized by unlimited delamination in all directions, as visualized in Fig. 3 by mechanism 3. Finally, it can be observed from Fig. 8a that for values of d / I towards 10 the crack channelling curveσ c ( / h 1 ) drops below the delamination curveσ d ( / h 1 ), as a result of which there are no intersection points between the two curves. The precise value of the toughness ratio at which this transition happens appears to be d / I = 2.16. Accordingly, for d / I ≥ 2.16 delamination remains absent during crack channelling, which is in agreement with the observation that the minimum value of the crack channelling curve in this range occurs at zero delamination, / h 1 = 0, see Fig. 8a. In Fig. 3 this fracture scenario is designated as mechanism 1. In summary, from the numerical results presented in Fig. 8a, it is concluded that crack channelling with delamina-tion absent (mechanism 1) is operational for toughness ratios in the range d / I ≥ 2.16, crack channelling with constant delamination (mechanism 2) occurs for 0.58 ≤ d / I < 2.16, and crack channelling with unlimited delamination in all directions (mechanism 3) is operational for d / I < 0.58. Note that the corresponding value of the minimum crack channelling stress in the coating,σ c,min , can be translated to a critical value for the uniform moisture content via Eq. (6) 1 , upon substituting for the condition m 1 = m 2 .
Consider now a stiffness mismatchĒ 2 /Ē 1 = 1 for which the numerical results are illustrated in Fig. 8b. The peak value for the delamination stressσ d is found at relatively small delamination length / h 1 ≈ 0.15, after which the delamination curve rapidly decreases towards its steady-state valueσ ss . Due to the small rising part of the delamination curveσ d ( / h 1 ), mechanism 2 only appears for a small range of toughness ratios, 0.44 ≤ d / I < 0.49, and is characterized by small delamination lengths, / h 1 ≤ 0.15. For d / I ≥ 0.49 and d / I < 0.44, mechanism 1 and mechanism 3 are operational, respectively. Figure 8c finally refers to a stiffness mismatch E 2 /Ē 1 = 10. Observe that, after a minor oscillation at / h 1 = 0.03, the delamination stressσ d drops monotonically towards its steady-state valueσ ss upon increasing delamination. Neglecting the initial, local oscillation in the delamination stress, from the decreasing trend of the delamination curveσ d ( / h 1 ) it follows that mechanism 2 does not occur. Accordingly, for low toughness ratios d / I ≤ 0.43 mechanism 3 is operational, while the bilayer system may fail in accordance with mechanism 1 when d / I > 0.43.
From the above identification procedure of the three crack channelling mechanisms, a failure mechanism map can be constructed in which the minimum crack channelling stressσ c,min in the coating is plotted as a function of the toughness ratio d / I for a selection of stiffness mismatchesĒ 2 /Ē 1 , see Fig. 9. This failure map thus provides an estimate of the critical crack channelling stress, given the values of the toughness ratio d / I and the elastic stiffness mis-matchĒ 2 /Ē 1 . The three failure mechanisms are represented by the regions defined by the dotted orange lines. For coating-substrate systems characterized by a relatively low toughness ratio d / I ≤ 0.43, mechanism 3 prevails for any of the stiffness mismatches considered. For higher toughness ratios d / I > 0.43, mechanism 1 operates at higher stiffness mismatches, Fig. 9 Failure mechanism map for a bilayer system with layers of equal thickness (h 2 / h 1 = 1). Minimum crack channelling stressσ c,min in the coating versus toughness ratio d / I for a constant moisture content profile, considering a broad selection of stiffness mismatchesĒ 2 /Ē 1 . The dotted orange lines define the regions corresponding to the three failure mechanisms presented in Fig. 3. The results are independent of the mismatch in coefficient of hygral expansion,β 2 /β 1 , and the mismatch in diffusion coefficient, D 2 /D 1 E 2 /Ē 1 ≥ 3, while for moderate to low stiffness mis-matchesĒ 2 /Ē 1 < 3 mechanism 3 gradually turns into mechanism 2, and finally into mechanism 1 under increasing toughness ratio d / I . Note that the transitions in failure mechanism occur at a larger toughness ratio d / I if the stiffness ratioĒ 2 /Ē 1 is lower.

Influence of the moisture profile on the cracking response via the parameter α
This section discusses the influence of the considered moisture profile on the predicted fracture response, by exploring the effect of the weighting factor α. For the bilayer system described in Sect. 4.1, the fracture scenarios are first analysed corresponding to a linear moisture content distribution with a vanishing moisture content value at the layer interface, which, for the case D 2 /D 1 = 1, is defined via Eq. (1) 1 as m(y) = ( m 1 − m 2 )y/ h. Correspondingly, the weighting factor α is set to zero, α = 0, yielding the remote stresses for crack channelling and planestrain delamination as σ ∞,α c (y) = c (1 + y/(hη 1 )) and σ ∞,α d (y) = d (1 + y/(hη 1 )), respectively. Next, an arbitrary linear moisture content profile, as represented by relation (1), is considered. The correspond- Fig. 10 Failure mechanism map for a bilayer system with layers of equal thickness (h 2 / h 1 = 1). Minimum crack channelling stress c,min in the coating versus toughness ratio d / I for a linear moisture content profile related to a weighting factor α = 0, considering a broad selection of stiffness mismatches E 2 /Ē 1 . The dotted orange lines define the regions corresponding to the three failure mechanisms presented in Fig. 3. The mismatch in coefficient of hygral expansion equalsβ 2 /β 1 = 10, and the ratio of diffusion coefficients is D 2 /D 1 = 1 ing remote channelling and plane-strain delamination stresses σ ∞,α c (y) = ασ c + c (1 + y/(hη 1 )) and σ ∞,α d (y) = ασ d + d (1 + y/(hη 1 )), respectively, are computed by assuming a weighting factor α = 0.5, and are evaluated at the interface between the coating and substrate, i.e. at y = 0. Note that the numerical results are now dependent on the mismatch in the coefficient of hygral expansion,β 2 /β 1 , and the mismatch in the diffusion coefficient, D 2 /D 1 , which quantify the remote stress measure 1 , and thus σ ∞ 1 (y), via the parameter η 1 , see Eqs. (6) 2,3 . The mismatch in coefficient of hygral expansion has been chosen asβ 2 /β 1 = 10, so that condition (8) 2 is satisfied for the full range of selected stiffness mismatches,Ē 2 /Ē 1 ∈ [0.1, 10].

Failure mechanisms for a linear moisture content distribution with a vanishing value at the layer interface (α = 0)
Similarly to that done in Sect. 4.2, the mode-mix , the energy release rates G c, , G d, and the channelling c and plane-strain delamination d stresses are predicted in the case of a linear moisture content distribution with a vanishing value at the layer interface (α = 0), as a function of the delamination length / h 1 . The obtained results present trends that are qual- Fig. 11 Cracking in a bilayer system with layers of equal thickness (h 2 / h 1 = 1). Possible failure mechanisms as a function of the elastic mismatchĒ 2 /Ē 1 and toughness ratio d / I . The light blue dotted lines and the dark blue dashed lines define the regions corresponding to the three failure mechanisms depicted in Fig. 3 under linear (α = 0) and uniform (α = 1) moisture content profiles, respectively. The results obtained under a linear moisture content profile relate to a mismatch in coefficient of hygral expansion ofβ 2 /β 1 = 10, and a ratio of diffusion coefficients D 2 /D 1 = 1 itatively comparable to those discussed for a uniform moisture profile (α = 1), and are presented in detail in Appendix 1. These results can be summarized in the failure mechanism map shown in Fig. 10, illustrating the regions in which the three crack channelling mechanisms depicted in Fig. 3 are active. In Fig. 10, the minimum channelling stress c,min in the coating is plotted as a function of the toughness ratio, d / I , for selected values of the stiffness mismatch,Ē 2 /Ē 1 . In the range of toughness ratios 0.4 ≤ d / I < 0.52, mechanism 1 directly converts into mechanism 3 for moderate to large values of the stiffness mismatchĒ 2 /Ē 1 ≥ 1, while for low stiffness mismatch valuesĒ 2 /Ē 1 < 1 mechanism 2 intervenes. However, since the corresponding values of constant delamination length are relatively small, with a maximum value of / h 1 = 0.3 forĒ 2 /Ē 1 = 0.1, in practice it may be difficult to distinguish mechanism 2 from mechanism 1.

Crack channelling mechanisms under linear
(α = 1) versus uniform (α = 0) moisture profiles The regimes of dominance of the three crack channelling scenarios depicted in Figs. 9 and 10 differ somewhat under uniform and linear moisture content profiles. In order to clearly illustrate these differences, Fig. 11 summarizes the failure scenarios activated under the two types of moisture variations in a single failure map in material property space, by plotting on a log-linear scale the elastic stiffness mismatch E 2 /Ē 1 as a function of the toughness ratio d / I . Here, the light blue dotted lines delineate the failure regions under the presence of a linear moisture content profile (whereby the weighting factor α = 0), and the dark blue dashed lines delineate failure regions under a uniform moisture content profile (whereby the weighting factor α = 1). It can be clearly seen that for stiffness mismatchesĒ 2 /Ē 1 ≥ 1 mechanism 2 is (virtually) absent, with the transition from mechanism 1 to mechanism 3 under uniform and linear moisture content profiles occurring at similar values of the toughness ratio d / I . Conversely, forĒ 2 /Ē 1 < 1 mechanism 2 intervenes, although, as discussed previously, the corresponding constant delamination length remains limited to a fraction of the coating thickness, / h 1 < 1. At the lowest value of the stiffness mismatch,Ē 2 /Ē 1 = 0.1, the range of toughness ratios across which mechanism 2 operates under a linear moisture content profile is substantially smaller than under a uniform moisture content profile. Hence, for the experimental identification of mechanism 2 it is recommended to use a bilayer specimen with a relatively compliant substrate that is exposed to a uniform moisture content profile.

Failure mechanisms for an arbitrarily linear moisture content profile (0 < α < 1)
The failure scenario corresponding to an arbitrarily linear moisture content profile is next analysed. Observe first from Fig. 6, illustrating the dependence of the steady-state stress σ ∞,α ss (0) on the parameter α, that the failure response for arbitrarily linear moisture profiles (0 < α < 1) is bounded by the cases α = 0 and α = 1. In the following, the crack channelling stress σ ∞,α c (0) and delamination stress σ ∞,α d (0) are thus calculated based on an intermediate value of the weighting factor α, i.e., α = 0.5. This is done by repectively substituting the results from Eqs. (26) and (40) in Eq. (5) 1 . The obtained plots for the remote crack channelling stress and delamination stress as a function of the delamination length result to be qualitatively comparable to those shown in Fig. 8 for a constant moisture profile Fig. 12 Failure mechanism map for a bilayer system with layers of equal thickness (h 2 / h 1 = 1). Minimum crack channelling stressσ ∞,α c,min (0) in the coating versus toughness ratio d / I for a linear moisture content profile related to a weighting factor α = 0.5, considering a broad selection of stiffness mismatches E 2 /Ē 1 . The dotted orange lines define the regions corresponding to the three failure mechanisms presented in Fig. 3. The mismatch in coefficient of hygral expansion equalsβ 2 /β 1 = 10, and the ratio of the diffusion coefficients is D 2 /D 1 = 1 related to α = 1, and therefore are omitted here. The failure mechanism map depicted in Fig. 12 for α = 0.5 also is comparable to those in Fig. 9 for α = 1 and in Fig. 10 for α = 0. However, for relatively low stiffness ratiosĒ 2 /Ē 1 < 0.3, the transition from mechanism 1 to mechanism 2 in Fig. 12 with α = 0.5 and in Fig. 9 for α = 1 clearly occurs at a larger toughness ratio d / I than in Fig. 10 with α = 0.

Steady-state delamination under varying relative substrate thickness
The influence of the relative layer thickness h 2 / h 1 on the steady-state delamination stress is studied by varying the substrate thickness h 2 in the range h 2 = [h 1 , 1000h 1 ]. The steady-state stressesσ ss , ss and σ ∞,α ss (0) generated in the coating under uniform and linear moisture profiles are respectively shown in Fig. 13a-c as a function of the relative layer thickness h 2 / h 1 . The constant and linear stress measures, σ ss (α = 1) and ss (α = 0), are calculated from Eq. (46), and the total stress σ ∞,α ss (0) is computed (a) (b) (c) Fig. 13 Influence of the relative substrate thickness h 2 / h 1 on the steady-state coating stress for unlimited plane-strain delamination, considering a broad selection of stiffness mismatches E 2 /Ē 1 . a Constant delamination stressσ ss (α = 1), Eq. (46) 1 , whereby the results are independent of the ratiosβ 2 /β 1 and D 2 /D 1 ; b Linear delamination stress ss (α = 0), Eq. (46) 2 , forβ 2 /β 1 = 10 and D 2 /D 1 = 1, and c Total delamination stress σ ∞,α ss (0) (with α = 0.5), Eq. (47), forβ 2 /β 1 = 10 and D 2 /D 1 = 1 from expression (47), with the weighting factor set equal to α = 0.5. The constant steady-state delamination stressσ ss in the coating is associated with a uniform moisture content distribution, and is independent of the mismatch in the coefficient of hygral expansion,β 2 /β 1 , see also expression (46) 1 . Figure 13a illustrates that its magnitude increases with increasing stiffness mismatchĒ 2 /Ē 1 . In addition, the delamination stress grows monotonically with increasing relative layer thickness h 2 / h 1 , and asymptotes to a limit valuē σ ss (h 2 / h 1 → ∞) = 2Ē 1 d / h 1 that is independent of the stiffness mismatch. Note that under increasing relative layer thickness h 2 / h 1 the limit value is reached earlier for a larger stiffness ratioĒ 2 /Ē 1 . Figure 13b, c shows that the steady-state delamination stress ss (α = 0) and the total stress σ ∞,α ss (0) (α = 0.5) in the coating follow a similar trend as the constant stress depicted in Fig. 13a. Although ss and σ ∞,α ss (0) are computed for a mismatch in the coefficient of hygral expansion ofβ 2 /β 1 = 10 and a ratio of the diffusion coefficients of D 2 /D 1 = 1, the effect of these specific values on the fracture response vanishes for a large substrate thickness. Specifically, at large values of the relative thickness, h 2 / h 1 ≥ 5, it can be observed that the failure stressesσ ss , ss and σ ∞,α ss (0) are almost identical, irrespective of the value of the stiffness mismatchĒ 2 /Ē 1 . As can be confirmed from the second term in Eq. (5) 1 , this implies that for a large substrate the linear variation in stress across the coating becomes relatively small, such that the stress within the coating is approximately uniform and thus independent of the values ofβ 2 /β 1 and D 2 /D 1 . Conversely, when the coating and substrate have a comparable thickness, h 2 / h 1 < 5, the linear variation in stress has a substan-tial effect on the total stress σ ∞,α ss (0), and needs to be explicitly accounted for in the failure analysis. Note again that under these circumstances the constant and linear stress measuresσ ss and ss respectively serve as upper and lower bounds for the total steady-state stress σ ∞,α ss (0) for any value of the stiffness mismatch E 2 /Ē 1 . Observe further that the linear failure stress ss approaches zero for a bilayer system with layers of equal thickness, h 2 / h 1 = 1, and a stiffness mismatch equal toĒ 2 /Ē 1 = 0.1, see Fig. 13b. This means that under a linear moisture content profile with a vanishing value at the layer interface this system has no resistance against delamination, and will fail spontaneously. However, the practical relevance of this particular case is small, since a linear moisture content profile in general will also induce a constant stress contribution (i.e., in practice α = 0), for which the resistance against delamination is finite, see Fig. 13a. 5.2 Failure mechanisms for a bilayer system with a thick substrate (h 2 / h 1 = 10) Consider now a bilayer system with a relatively thick substrate characterized by h 2 / h 1 = 10. Under the application of a linear moisture profile associated with α = 0.5, the generated remote crack channelling and delamination stresses σ ∞,α , as obtained by respectively inserting the results from Eqs. (26) and (40) into Eq. (5) 1 , are shown in Fig. 14a-c as a function of the delamination length / h 1 for stiffness mismatchesĒ 2 /Ē 1 = 0.1, 1 and 10, respectively.
As argued in the previous section, for such a thick substrate the normalized failure stresses are mainly governed by the uniform part of the moisture content distribution. It is interesting to observe that for the stiffness mismatchesĒ 2 /Ē 1 = 0.1 andĒ 2 /Ē 1 = 1 the rising part of the delamination curve σ ∞,α d (l/ h 1 ) is more prominent than for the case of a bilayer system with layers of equal thickness discussed in Sect. 4; for both stiffness mismatches the increase of the delamination stress continues up to a stable delamination length of / h 1 ≈ 4, which is sufficiently large to adequately identify mechanism 2 in practical applications. The range of toughness mismatches for which mechanism 2 is operational is 0.10 ≤ d / I < 0.73 for E 2 /Ē 1 = 0.1 and 0.22 ≤ d / I < 0.41 forĒ 2 /Ē 1 = 1. Above and below these ranges crack channelling occurs by means of mechanism 1 and mechanism 3, respectively. ForĒ 2 /Ē 1 = 10 the toughness range related to mechanism 2 is minor and corresponds to relatively small delamination lengths, / h 1 < 0.1; hence, under an increasing toughness ratio mechanism 3 more or less directly converts into mechanism 1 at a value d / I ≈ 0.4. The above-mentioned regions of the three failure mechanisms are designated in more detail in the failure mechanism map presented in Fig. 15; the figure clearly shows that mechanism 2 indeed becomes Fig. 15 Failure mechanism map for a bilayer system with a relatively thick substrate (h 2 / h 1 = 10). The mismatch in coefficient of hygral expansion equalsβ 2 /β 1 = 10 and the ratio of the diffusion coefficients is D 2 /D 1 = 1. Minimum crack channelling stress σ ∞,α c,min (0) in the coating versus toughness ratio d / I for a linear moisture content profile related to a weighting factor α = 0.5, considering a broad selection of stiffness mismatches E 2 /Ē 1 . The dotted orange lines define the regions corresponding to the three failure mechanisms presented in Fig. 3 a distinguished failure mechanism at lower values of the stiffness ratioĒ 2 /Ē 1 , i.e., for a bilayer system with a relatively compliant substrate.

Comparison of failure response for bilayer
systems with h 2 / h 1 = 10 and h 2 / h 1 = 1 The results obtained in Figs. 12 and 15 are summarized in Fig. 16, illustrating the operative failure mechanisms as a function of the elastic stiffness mismatch E 2 /Ē 1 and toughness ratio d / I . Dark blue dashed lines refer to a system with layers of equal thickness (h 2 / h 1 = 1) and light blue dotted lines refer to a system characterized by a relatively thick substrate (h 2 / h 1 = 10). It can be seen that the toughness range d / I within which mechanism 2 is active generally corresponds to lower values when the substrate is thicker. Also, at larger stiffness ratiosĒ 2 /Ē 1 > 1 (i.e., for systems with a relatively stiff substrate) mechanism 2 vanishes for a bilayer with h 1 = h 2 , while for the case with h 2 = 10h 1 it remains present. It is further instructive to display the critical total stress for crack channelling σ ∞,α c,min (0) (using a weighting factor α = 0.5) for mechanisms 1 and 3 as a function of the stiffness mismatchĒ 2 /Ē 1 , considering the Fig. 16 Cracking in a bilayer system for different relative layer thicknesses and a linear moisture content profile related to a weighting factor α = 0.5. The mismatch in coefficient of hygral expansion equalsβ 2 /β 1 = 10 and the ratio of the diffusion coefficients is D 2 /D 1 = 1. Possible failure mechanisms as a function of the elastic mismatchĒ 2 /Ē 1 and toughness ratio d / I . The dark blue dashed lines and the light blue dotted lines define the regions corresponding to the three failure mechanisms depicted in Fig. 3 for h 2 / h 1 = 1 and h 2 / h 1 = 10, respectively two relative layer heights, h 2 / h 1 = 1 and h 2 / h 1 = 10, see Fig. 17.
The critical channelling stress in the coating is normalized by the term (Ē 1 c )/ h 1 , where c = I when referring to mechanism 1 (indicated by the solid lines) and c = d when denoting mechanism 3 (indicated by the dashed lines). The critical channelling stress for mechanism 1 is calculated from the FEM results presented in Sects. 4.3.3 and 5.2, by taking the limit σ ∞,α c,min = σ ∞,α c ( / h 1 → 0), while for mechanism 3 it follows from insertingσ 1 = ασ ss , together with the fracture criterion (32), into Eq. (44), and solving for 1 , after which the solution is substituted into expression (47). In all cases the critical channelling stress monotonically increases with increasing value of the stiffness mismatch. It can be further observed that for both mechanisms the corresponding critical channelling stress is lower for the system with the smaller relative substrate thickness, h 2 / h 1 = 1. The figure may serve as a design graph provided the failure mechanism is known. Note that mechanism 2 can intervene, depending on the value of the toughness mismatch d / I , as demonstrated by the failure mechanism maps depicted in Figs. 12 and 15. Fig. 17 Cracking in a bilayer system for different relative layer thicknesses and a linear moisture content profile related to a weighting factor α = 0.5. The mismatch in coefficient of hygral expansion equalsβ 2 /β 1 = 10 and the ratio of the diffusion coefficients is D 2 /D 1 = 1. Critical channelling stress σ ∞,α c,min (0) in the coating for mechanism 1 and mechanism 3 as a function of the elastic mismatchĒ 2 /Ē 1 . The dark blue thin lines and the light blue thick lines refer to the relative layer heights h 2 / h 1 = 1 and h 2 / h 1 = 10, respectively. In the normalization of the minimum channelling stress along the vertical axis, for mechanisms 1 and 3 the toughness parameter is taken as c = I and c = d , respectively

Conclusions
Crack channelling is addressed for a brittle coatingsubstrate system subjected to a moisture (or temperature) gradient in the thickness direction. Three distinct failure scenarios have been identified: (i) channelling of a mode I crack in the coating with delamination absent, (ii) channelling of a doubly deflected crack with constant delamination length, and (iii) channelling of a mode I crack with unstable, unlimited delamination in all directions. Failure mechanism maps have been constructed, illustrating the dependence of the active crack channelling mechanism and the corresponding critical channelling stress to the ratio of layer interface to coating toughness, and to the mismatches in stiffness and in coefficient of hygral expansion. The influence of the thickness ratio of the coating and substrate on the critical channelling stress and failure mechanism has also been explored. Closed-form expressions are derived for the steady-state delamination stresses, which allow one to determine the critical moisture conditions that lead to unlimited delamination. Due to the general character of the study, the results can be applied to coating-substrate systems in various applications. Specifically, the study serves to predict cracking phenomena generated in historical paintings under indoor climate fluctuations, thereby helping museums to preserve their art objects. In particular, failure mechanisms 2 and 3 (crack channelling with finite and unlimited delamination, respectively) are the most critical for a painting and must be avoided, while mechanism 1 (crack channelling with no delamination) is the preferential failure regime. From the proposed failure maps, it is possible to derive (i) the material property range in which mechanism 1 is active; (ii) the stress level at which it occurs and (iii) the corresponding critical moisture content values. This information may thus provide museums guidelines on appropriate conservation interventions.
Acknowledgements E.B. is grateful for the support by the Netherlands Organization for Scientific Research (NWO), Project 15873, "Engineering goes Beauty-A computational multiphysics modelling approach towards the preservation of historical oil paintings", within the funding scheme "NWO Veni Award". A.S.J.S. has received funding from the European Union's Horizon 2020 research and innovation programme under grant agreement N. 814624. The authors would like to thank Matteo Rossi Doria of CBC Conservazione Beni Culturali Società Cooperativa in Rome for providing images of cracking mechanisms in historical paintings.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/ by/4.0/. Appendix 1: Cracking due to a linear moisture content distribution with a vanishing value at the layer interface (α = 0) The results for the bilayer system presented in Sect. 4.1, characterized by equal diffusion coefficients D 1 = D 2 , and subjected to a linear moisture content distribution with a vanishing moisture content value at the layer interface, i.e., m(y) = ( m 1 − m 2 )y/(h 1 + h 2 ), are presented in detail in the following. A1.1 Mode-mix and energy release rate The phase angle generated by the linear moisture content variation is depicted in Fig. 18a as a function of the delamination length / h 1 . In alignment with Eq. (30), the mode-mixity is calculated as = atan(Im(K l i )/Re(K l i )), whereby the complex stress intensity factor K is defined by K = K 1, + i K 2, , with K 1, and K 2, as given by relation (48). The reference length at which the modemix is evaluated is taken equal to the layer height l = h 1 . The phase angle attains a steady state at a delamination length / h 1 ≈ 0.8 for all stiffness mismatches selected, with the asymptotic value lying within the range [46 • , 54 • ]. The phase angle at steadystate typically decreases with increasing stiffness mis-matchĒ 2 /Ē 1 , with the exception being the slightly lower steady-state value forĒ 2 /Ē 1 = 0.3 compared tō E 2 /Ē 1 = 0.1.
The normalized energy release rates for plane-strain delamination G d, and crack channelling G c, are shown respectively in Fig. 18b, c, as a function of the delamination length / h 1 . The results forĒ 2 /Ē 1 = 0.1 have been omitted in this figure, due to the relatively large values found for this case. The trends observed for the energy release rates are comparable to those calculated for the constant moisture profile, see Fig. 7b, c, and asymptote to the analytical steady-state value G ss, as given by Eq. (45) 2 A1.2 Failure mechanisms Figure 19 contains plots of the stresses d for planestrain delamination and c for steady-state crack channelling as a function of the delamination length / h 1 , for a linear moisture distribution. Figure 19a, b respectively refer to stiffness mismatchesĒ 2 /Ē 1 = 1 and 10. The results for a stiffness mismatch ofĒ 2 /Ē 1 = 0.1 are not depicted here, since the remote stress values calculated for this case almost vanishes. The channelling stress c , indicated in Fig. 19 by the dark blue lines, has been computed using Eq. (22) for the following selected values of toughness ratio, d / I = [0.05, 0.1, 0.2, 0.5, 1, 10]. The delamination stress d , designated by the light blue line, has been calculated by means of Eq. (38). For a system characterized by a mismatch in elastic stiffness ofĒ 2 /Ē 1 = 1, mechanism 2 operates along the small rising part of the (a) (b) (c) Fig. 18 Cracking in a bilayer system with layers of equal thickness (h 2 / h 1 = 1). The mismatch in coefficient of hygral expansion equalsβ 2 /β 1 = 10 and the ratio of diffusion coefficients is D 2 /D 1 = 1. a Mode-mix = atan(Im(K h i 1 )/Re(K h i 1 )) as a function of the delamination length / h 1 , for a linear moisture content profile with a vanishing value at the layer interface.
Energy release rate for b plane-strain delamination G d, , Eq. (37), and c steady-state crack channelling G c, , Eq. (21), as a function of the delamination length / h 1 , for a linear moisture content profile that vanishes at the layer interface. The depicted steady-state values G ss, are computed with Eq. (45) 2 c ( / h 1 ) curve, see Fig. 19a. The toughness range over which this occurs is limited, 0.40 ≤ d / I < 0.44, whereby the constant delamination length remains relatively small, / h 1 ≤ 0.1. For toughness values above and below this range mechanism 1 and mechanism 3 become active, respectively. At a large value of the stiffness mismatch,Ē 2 /Ē 1 = 10, apart from a minor, local oscillation the delamination curve is monotonically decreasing, see Fig. 19b, so that a transition from mechanism 1 to mechanism 3 takes place under an increasing toughness ratio at a value d / I = 0.45. Crack channelling and plane-strain delamination in a bilayer system with layers of equal thickness (h 2 / h 1 = 1). The mismatch in coefficient of hygral expansion equalsβ 2 /β 1 = 10, and the ratio of diffusion coefficients is D 2 /D 1 = 1. Remote stress 1 versus delamination length / h 1 for a linear moisture content profile, considering a stiffness mismatch of aĒ 2 /Ē 1 = 1 and bĒ 2 /Ē 1 = 10. The dark blue lines represent the crack channelling stress 1 = c , equation (22), for selected toughness ratios d / I . The light blue line reflects the plane-strain delamination stress 1 = d , equation (38). The depicted steady-state value ss is computed with Eq. (46) 2