Fracture Mechanics Solutions for Interfacial Cracks between Compressible Thin Layers and Substrates

The decohesion of coatings, thin films, or layers used to protect or strengthen technological and structural components causes the loss of their functions. In this paper, analytical, computational, and semi-analytical 2D solutions are derived for the energy release rate and mode-mixity phase angle of an edge-delamination crack between a thin layer and an infinitely deep substrate. The thin layer is subjected to general edge loading: axial and shear forces and bending moment. The solutions are presented in terms of elementary crack tip loads and apply to a wide range of material combinations, with a large mismatch of the elastic constants (isotropic materials with Dundurs’ parameters − 1 ≤ α ≤ 1 and − 0.4 ≤ β ≤ 0.4 ). Results show that for stiff layers over soft substrates ( α → 1 ), the effects of material compressibility are weak, and the assumption of substrate incompressibility is accurate; for other combinations, including soft layers over stiff substrates ( α → − 1 ), the effects may be relevant and problem specific. The solutions are applicable to edge- and buckling-delamination of thin layers bonded to thick substrates, to mixed-mode fracture characterization test methods, and as benchmark cases.


Introduction
The decohesion of thin layers from thicker substrates is a well-known failure mechanism which may occur at very different scales in structural and technological components due to various loading conditions, and is controlled by the interfacial fracture toughness, a material property which depends on the relative magnitude of shear and normal tractions on the interface [1][2][3][4]. Delamination (or debonding) of the face sheets from the core is one of the dominant failure mechanisms of thick sandwich structures used in naval, aeronautical, automotive, and wind energy applications, and can be due to direct loadings, stress waves, or internal pressure cycles [4][5][6][7][8]. In these systems, the face sheets are typically laminate composite materials, and the core is made up of cellular foams, honeycomb or balsa; exemplary applications are naval bulkheads and deck panels, primary aerostructures, e.g., wing leading edges and rudders, and wind turbine blades. The debonding of Fiber Reinforced Polymer (FRP) sheets used to repair or strengthen reinforced concrete members in civil engineering applications prevents load transfer and nullifies the reinforcing action [9,10]. In thermal barrier coatings designed to reduce the temperature of the underlying metallic substrates for aircraft and power-generating turbines, delamination and spalling are common failure modes [11][12][13][14]. The cracking of thin films used in electronic packages or of metal alloy coatings on elastomers, used in soft electronics, is often driven in soft electronics, is often driven by residual stresses [2,3,[15][16][17][18]. The cracking and spalling of brittle layers may be a consequence of edge loadings [19]. The decohesion of polymer coatings used to protect metals for different applications, e.g., in beverage or food cans, biomedical and mechanical devices, or with anti-fouling functionality, leads to the loss of the protective and aesthetic properties of the coatings [20,21]. Many other thin film applications depend on the adhesion of the film to the substrate [22][23][24].
This paper presents fracture mechanics solution methods and results for the fundamental 2D problem in Figure 1, which describes an edge delamination crack between a thin layer and a semiinfinite substrate. The thin layer is subjected to general edge loadings and mode mixity phase angle and energy release rate associated with the collinear propagation of the crackare obtained in terms of crack tip force and moment resultants, Figure 1b [25][26][27][28] to provide analytical solutions for axial, bending, and shear loading with 0 β = , an efficient and accurate numerical method for the derivation of the fracture mechanics parameters in the numerically challenging film/substrate problem, and numerical coefficients for the application of the semi-analytic solutions in [25,26] to bi-materials with a general mismatch of the elastic constants and For its technological importance, the problem in Figure 1 has been studied and used extensively in the literature since the early work by Obreimoff on the splitting strength of mica [1][2][3][29][30][31][32]. Most studies and applications are limited to conventional bi-material systems with 0.8 0.8 α − ≤ ≤ and neglect the effects of the mismatch of the volumetric stiffness of the layers; this is done by assuming the second Dundurs' parameter 0 β = "because of its secondary role". In addition, most studies neglect the effects of crack tip shear, V in Figure 1b, because of its limited importance when dealing with long cracks. However, new material combinations with a larger mismatch of the elastic constants are increasingly being used [2,3,8,11,[33][34][35], where the effects of 0 β ≠ on the fracture parameters might be important. Table 1 lists examples of thin-layer/substrate material combinations used in electronics, structural composite sandwiches, and thermal barrier coatings. In addition, shear does For its technological importance, the problem in Figure 1 has been studied and used extensively in the literature since the early work by Obreimoff on the splitting strength of mica [1][2][3][29][30][31][32]. Most studies and applications are limited to conventional bi-material systems with −0.8 ≤ α ≤ 0.8 and neglect the effects of the mismatch of the volumetric stiffness of the layers; this is done by assuming the second Dundurs' parameter β = 0 "because of its secondary role". In addition, most studies neglect the effects of crack tip shear, V in Figure 1b, because of its limited importance when dealing with long cracks. However, new material combinations with a larger mismatch of the elastic constants are increasingly being used [2,3,8,11,[33][34][35], where the effects of β = 0 on the fracture parameters might be important. Table 1 lists examples of thin-layer/substrate material combinations used in electronics, structural composite sandwiches, and thermal barrier coatings. In addition, shear does contribute to the fracture parameters in many applications and thin-layer tests, such as the blister/bulge test and the peel test for thin films [26,31], or the Single Cantilever Beam (SCB) test for sandwich structures [36]. The analytical and semi-analytical expressions presented in this paper can be used for the solution of various film/substrate problems, for instance for interface toughness measurements or to study buckle-driven delamination. They can be used to validate and support the simplifying assumptions. They also provide a useful reference limiting solutions for problems where the substrate is not semi-infinite, e.g., in SCB sandwich tests. Table 1. Dundurs' parameters for bi-material interfaces (Mat.#1/Mat.#2) with a large mismatch of the elastic constants (plane strain conditions) [2,3,8,11,[33][34][35]; β is expressed here as a function of α and the Poisson ratios. The paper is organized as follows. Analytical solutions are presented in Section 2 and Appendix A for the interface crack in Figure 1 with a → ∞ , β = 0 and a general mismatch of the effective Young's moduli, −1 ≤ α ≤ 1. A numerical model which implements a crack surface displacement extrapolation method into the finite element code ANSYS for the accurate definition of the fracture parameters is presented in Section 2 and Appendix B; the model is verified in Section 3 using the analytical results. In Section 2, the semi-analytic solutions in References [25,26] are particularized to the problem in Figure 1, and the numerical coefficients are derived and presented in tabular form in Section 3 for bi-materials with a general mismatch of the elastic constants and Dundurs' parameters −0.99 ≤ α ≤ 0.998 and −0.4 ≤ β ≤ 0.4. Extreme asymptotics and intermediate expansions of the fracture parameters are derived in Section 3. A discussion on the results is presented in Section 4.

Materials and Methods
The problem in Figure 1a, representing an edge-delamination crack of length a between a layer of thickness h 1 and a semi-infinite plane h 1 /h 2 → 0 , subjected to general edge loadings per unit width, P, M , V with M = M − Va, is examined. The crack tip loads are defined by the elementary, self-equilibrated systems in Figure 1b: pure bending moments, M, axial forces, P (plus the compensating moment, M * = M + Ph 2 /2, since h 1 /h 2 → 0 ) and shear forces, V [25,26].
The layer, material #1, and the half plane, material #2, are linearly elastic and isotropic, with Young's moduli and Poisson ratios E 1 , ν 1 and E 2 , ν 2 ; the interface between the two materials is abrupt, smooth, and adherent for x ≥ 0. The mismatch of the elastic constants is described by the Dundurs' parameters [37]: where E i = E i and κ i = (3 − ν i )/(1 + ν i ) for plane-stress, and E i = E i /(1 − ν i 2 ) and κ i = 3 − 4ν i for plane-strain; G i = E i /(2 + 2ν i ) is the shear modulus. The ratio of the effective Young's moduli is Σ = E 1 /E 2 and is related to α through α = (Σ − 1)/(Σ + 1). The parameters α, β define measures of the mismatch in uniaxial and bulk compliances of the two materials.
The crack length a is sufficiently long to ensure that the crack tip fields are unaffected by the actual distribution of the end tractions; since the layer is homogeneous and isotropic, this condition is satisfied if a is a few times the layer thickness, h 1 (Saint Venant's principle). Details on the applicability of these solutions to systems where h 1 /h 2 = 0, namely with finite substrates or finite thickness layers, are given in Section 4.
The stress and displacement fields in the singularity dominated zone at the crack tip are described by the stress and displacement components, σ yy , σ xy and u y , u x . For general interfaces, they are oscillating [38,39], and the oscillations are described by the oscillatory index ε given in Equation (1): with r is the distance from the crack tip (ahead for θ = 0 or behind for θ = ±π), Figure 1b, E * = (1 − α)E 1 , and K = K 1 + iK 2 is the complex stress intensity factor [1,40]. The oscillations imply crack surface interpenetration at the crack tip, which can be discounted if contact occurs in a small region near the crack tip (small-scale contact) [40]. For β = 0, the index ε = 0, the fields are similar to those in a homogeneous material and K I = K 1 and K I I = K 2 define the mode I and mode II components of the stress intensity factor. The energy release rate for the collinear propagation of the crack is related to the modulus of the complex stress intensity factor, |K| = K 1 2 + K 2 2 1/2 , by: [41]. The crack tip conditions are generally mixed mode and defined by the mode-mixity phase angle.
For bi-material interfaces with β = 0, K = K I + iK I I = Ke iψ , and the mode mixity phase angle is uniquely defined by: For bi-material interfaces with β = 0, the mode-mixity phase angle is not uniquely defined due to the oscillating fields. Here, following [40], the mode-mixity angle is related to the ratio between the shear and normal stresses at a reference distance,r, ahead of the crack tip: Results will be presented in this paper assumingr = h 1 . The choice ofr is arbitrary, since the phase angle at a different reference length is given by ψr = ψ h1 + ε ln(r/h 1 ) [40]; as a consequence, the interfacial fracture toughness defined through fracture experiments will depend on the mode mixity and on the value ofr used to define ψr [35,42].
For most applications of interface mechanics in the literature, the effect of the second Dundurs' parameter on the solution is neglected by assuming β = 0. This is done to simplify the treatment of the problem and avoid dealing with the oscillatory fields. This assumption is valid for many material combinations, especially when α is small, e.g., −0.4 ≤ α ≤ 0.4, since β would then vary in the approximate range −0.1 < β < 0.1, which is exact for ν 1 = ν 2 = 1/3 under plane-strain conditions [34]. However, the assumption is questionable for other cases and larger values of α (see Table 1).

Analytical Model: Materials with β = 0
An analytical solution of the problem in Figure 1 for β = 0 has been derived in [28] using the Laplace transform and reducing the problem to a homogeneous matrix Riemann-Hilbert problem following [46,47]; see Appendix A for more details. The solution has been obtained under the assumption that the loads are applied at an infinite distance from the crack tip. The crack is in mixed-mode conditions with the Mode I and mode II stress intensity factors given by: where and ω and δ are given in Equations (A1) and (A2); they depend on Σ = E 1 /E 2 (or α) and have been tabulated for ease of application in Table 2. Table 2. Selected values of ω (in radians) and δ on varying α for materials with β = 0. The mode-mixity phase angles associated with the elementary loads in Figure 1b are derived from Equations (4), (7) and (8): (9) and the dimensionless functions which define the elementary energy release rates in Equation (6) are obtained as: where f M and f P coincide with those derived in [25]. The analytical results for ψ P = ω, ψ V , and f V will be presented in Section 3.3 for selected values of α; extreme asymptotics and intermediate expansions, which describe very soft and stiff substrates and near homogeneous problems, will be derived in Section 3.3. The global mode mixity angle when P, M, V = 0 in Figure 1, is obtained through Equations (4) and (7), and the global energy release rate in Equation (3) becomes: which coincides with results in [25] when V = 0.

Computational Model
The finite element model used to analyze the problem in Figure 1 is presented in Appendix B. The model implements the Crack Surface Displacement Extrapolation (CSDE) technique presented in [48] to define the energy release rate and mode-mixity phase angle. The CSDE is based on a linear extrapolation technique, where the energy release rate and mode-mixity phase angle calculated from relative nodal pair displacements at locations behind the crack tip are extrapolated linearly to the crack tip location. The energy release rate and mode-mixity phase angle are defined in terms of relative crack opening, ∆u y , and shearing, ∆u x , displacements at a distance r behind the crack tip [40]: The displacements are calculated using the finite element model presented in Appendix B.
The CSDE technique was originally developed to solve convergence problems for high stiffness mismatch bi-material interfaces, e.g., sandwich face/core interfaces, encountered with other mode separation methods, such as the Virtual Crack Closure Technique (VCCT), where nodal forces at or close to the crack tip are utilized to calculate the mode-mixity. Beuth [49] showed that the traditional VCCT by Rybicki and Kanninen [50] is not able to generate converged results for the mode-mixity due to the near-tip stress and displacement oscillations. Therefore, a modified version of the VCCT was presented in [49] where the near-tip oscillations are neutralized analytically during the post-processing of the near tip nodal force and displacement results to calculate the mode-mixity, and in this way achieve converged mode-mixity results. However, in most numerical finite element representations for applied bi-material crack problems, nodal results at or close to the crack tip are often associated with high numerical errors. This is both due to high stiffness mismatches over the interface, and thus the inability of the elements to accurately capture the complex stress field near the crack tip, but also to the inability of most models to capture the oscillations close to the crack tip. As a result, even the modified VCCT presented in [49] is unable to accurately predict the mode-mixity for high stiffness mismatch bi-material problems, like sandwich face/core interfaces [51]. The CSDE method, due to its linear extrapolation technique, bypasses the numerical error-zone close to the crack tip, and therefore constitutes a very stable and robust mode separation technique which is able to calculate the mode-mixity phase angle for even very challenging extreme stiffness mismatch bi-material crack problems [51].

Semi-Analytical Model-General Material Mismatch β = 0
Following the formulations in [25,26], the energy release rate and mode mixity phase angle for the bi-material system in Figure 1 are defined in terms of crack tip force and moment resultants, P, M, V, by operating in the complex plane of the stress intensity factors and noting that for each elementary load, e.g., V, the following relations hold: This and Equation (6) where ψ M = ω − π/2. For β = 0, f V , ω and ψ V , are given in Equations (9) and (10), and G coincides with Equation (11); the expressions in Equation (13) are then fully analytical. For β = 0, ω(α, β), f V (α, β) and ψ V (α, β) have to be defined numerically; this will be done in Section 3.3. A discussion on the applicability of the analytical and semi-analytical solutions to problems where substrate and crack have finite sizes is presented in Section 4.

Results
In this section, the accuracy of the finite element model and the CSDE technique for the challenging problem in Figure 1 is verified using the analytical results for ω, ψ V and f V in materials with β = 0. The finite element model is then applied to derive the dimensionless coefficients of the semi-analytic solutions in Equation (13) for β = 0. Extreme asymptotics and intermediate expansions for the coefficients of the analytical solutions are derived to facilitate the application of the analytical results to systems with no mismatch in the volumetric stiffnesses.

Finite Element Results for Materials with
The finite element model, Appendix B, and the CSDE technique implemented in ANSYS are first applied to define the energy release rate, stress intensity factors, and mode mixity angles associated with the elementary loads in Figure 1b for a bi-material interface with Dundurs' parameter β = 0. The phase angle ψ P = ω associated to the elementary axial load P (plus the compensating bending moment) are presented in Table 3 and Figure 2 on varying α and compared with the analytical solutions in Equations (7)- (9). Phase angle, energy release rate, and stress intensity factors for the elementary load of pure crack tip shear, V, are presented in Tables 4 and 5 and Figure 3.
The absolute error between the finite element based CSDE results and the analytical results is always below 0.04 • on the phase angle, and below 0.001 on the dimensionless stress intensity factors for −0.8 ≤ α ≤ 0.8; errors are slightly larger for −0.99 ≤ α ≤ −0.9 and 0.9 ≤ α ≤ 0.99, but always below 0.2 • and 0.01. The results have also been successfully compared with solutions in the literature [25,26] (a discrepancy is observed in the solutions for pure shear loading and α ≤ 0 in [26]; given the agreement with the analytical solution, the present results appear to be more accurate). The results in Figures 2 and 3 show large variations for α > 0.8, which corresponds to stiff layers on soft material substrates. These cases will be examined later in Section 3.3, where extreme asymptotics of the analytical solutions will be derived. (13)-Materials with β = 0 The finite element model and the CSDE technique have been applied to the geometry in Figure 1 to define energy release rate and mode mixity phase angle for the elementary loads in bi-material interfaces with general elastic mismatch and β = 0.

Tabulated Coefficients for the Semi-Analytic Solutions in Equation
The dimensionless function f V = E 1 hG V /V and the mode-mixity phase angle ψ V are presented in Tables 4 and 5 for −0.99 ≤ α ≤ 0.998 and −0.4 ≤ β ≤ 0.4. The phase angle ψ P = ω for −0.99 ≤ α ≤ −0.9 and 0.9 ≤ α ≤ 0.998 and −0.4 ≤ β ≤ 0.4, which has not been presented before in the literature, has also been derived, and results are listed in Table 3; these results can be combined with those tabulated in [25] for −0.8 ≤ α ≤ 0.8 to cover all material combinations (a verification performed on the results in [25] using the analytical model, Equation (7), and the finite element model formulated here indicates an uncertainty of ±0.5 • ). The uncertainties are given in the tables and have been defined through comparison with the analytical results and results in [25,26]. They are mainly determined by cases with α → −1 , for which the numerical results are less accurate. However, as it evident also from Figures 3b, 4a, 5 and 6, the coefficients approach a plateau for large negative values of α and the numerical derivation for α → −1 is not necessary.

Explicit Solutions for Soft/Stiff Material Substrates and Nearly Homogeneous Systems-β = 0
To facilitate the application of the analytical results presented in Section 2, extreme asymptotics and intermediate expansions of ω, ψ V and f V are derived for stiff layers over soft substrates, where Σ = E 1 /E 2 >> 1 is very large and α → 1 , for soft layers over stiff substrates, where Σ = E 1 /E 2 << 1 is very small and α → −1 , and near homogeneous cases, where Σ = E 1 /E 2 ≈ 1 and α ≈ 0. They are listed below.

•
Phase angle associated with the axial force, ψ P = ω (in radians) [28]: The uncertainties are below 0.2 • for E 1 /E 2 ≥ 330 (α ≥ 0.993) in Equation (14), for E 1 /E 2 ≤ 0.4 (α ≤ −0.43) in Equation (15), and for 0. (16). These uncertainties imply a relative percent error always lower than 0.4%. The equations can be used with confidence for all values of α outside the interval 0.51 < α < 0.993 where the exact analytical solutions in Equation (A1) or interpolation of the results in Table 3 should be applied. Table 3. Mode-mixity phase angle associated with the axial force, Figure 1a, with P = 0, M = V = 0.  • Phase angle associated with shear, ψ V (in radians): The uncertainties are below 0.2 • for E 1 /E 2 ≥ 330 (α ≥ 0.994) in Equation (17), for E 1 /E 2 ≤ 0.37 (α ≤ −0.46) in Equation (18), and for 0 ≤ E 1 /E 2 ≤ 6.2 (−1 ≤ α ≤ 0.72) in Equation (19). These uncertainties imply a relative percent error higher than that on ω, since ψ V is much smaller, which can be as high as 9% (Figures 2 and 3). However, this occurs when ψ V is closed to zero, and the problem is essentially mode I, and the error is therefore irrelevant for engineering applications. The equations can be used with confidence for all values of α outside the interval 0.72 < α < 0.994, where the exact analytical solutions in Equations (9) and (A2), or interpolation of the results in Table 4, should be applied.
The uncertainties are below 0.2° for 1 Equation 19. These uncertainties imply a relative percent error higher than that on ω , since V ψ is much smaller, which can be as high as 9% (Figures 2 and 3). However, this occurs when V ψ is closed to zero, and the problem is essentially mode I, and the error is therefore irrelevant for engineering applications. The equations can be used with confidence for all values of α outside the interval 0.72 0.994 α < < , where the exact analytical solutions in Equations (9) and (A2), or interpolation of the results in Table  4, should be applied.
The uncertainties are below ±0.5 • . Analytical results are shown in brackets.

•
The dimensionless function which defines the energy release rate associated with shear, f V : The relative error is below 1% for E 1 /E 2 ≥ 200 (α ≥ 0.99) in Equation (20), for E 1 /E 2 ≤ 0.6 (α ≤ −0.25) in Equation (21), and for 0.4 ≤ E 1 /E 2 ≤ 2.1 (−0.44 ≤ α ≤ 0.36) in Equation (22). The equations can be used with confidence for all values of α outside the interval 0.36 < α < 0.99 where the exact analytical solutions in Equations (10) and (A2), or the results in Table 5, should be applied. If a maximum relative error of 5% is acceptable, the validity of Equation (20) can be extended to all E 1 /E 2 ≥ 0.36 (α ≥ −0.47). Table 5. The dimensionless function defining the energy release rate associated with the elementary shear load, Figure 1, with V = 0, P = M = 0. The asymptotics for small Σ = E 1 /E 2 and for Σ → 1 have been obtained by developing the integrands of the corresponding values into series followed by integration. The asymptotics for large Σ = E 1 /E 2 have been obtained by first substituting s = s Σ 1/3 into Equations (A1) and (A2) and then developing them into series. The asymptotics (20) which define the energy release rate associated with shear in the limit of a soft substrate coincides with that obtained by approximating the thin layer as a beam (plate) [52]. An approximate expression for ω, valid in the range −0.8 ≤ α ≤ 0.8, and obtained through a numerical fit of the results in [25] is presented in [3]. Figure 4 shows results for ψ P = ω as a function of E 1 /E 2 in a logarithmic scale to highlight the behaviors. The analytical results are compared with the numerical results obtained here and presented in Table 3 for −0.99 ≤ α ≤ −0.8 and 0.8 ≤ α ≤ 0.999, and with those [25] for −0.8 ≤ α ≤ 0.8. The figure shows the two asymptotics in Equations (14) and (15) for small and large values of E 1 /E 2 and the asymptotics presented in [2] for large E 1 /E 2 . The latter has been obtained using the scaling, K I /K I I ≈ (E 1 /E 2 ) −1/3 , derived from the model of a beam on a semi-infinite substrate, and a numerical fit for α = 0.95; comparison with the analytical results and the analytical asymptotic in Equation (14) shows that the approximation is quite accurate for large values of Σ = E 1 /E 2 .
Coatings 2019, 9, x FOR PEER REVIEW 12 of 20 , derived from the model of a beam on a semi-infinite substrate, and a numerical fit for 0.95 α = ; comparison with the analytical results and the analytical asymptotic in Equation (14) shows that the approximation is quite accurate for large values of (a) (b)    Figure 5 shows results for ψ V as function of E 1 /E 2 on a logarithmic scale to highlight behavior in the extreme cases. The analytical results are compared with the numerical results in Table 4, for −0.99 ≤ α ≤ 0.998. The figure also shows the two asymptotics in Equations (17) and (18), for small and large values of E 1 /E 2 , and the intermediate expansion in Equation (19).

Discussion
Analytical expressions for the energy release rate and mode-mixity phase angle associated with the elementary loads in Figure 1b, Equations (9) and (10), have been derived from the results in Reference [28] for materials with 0 β = , and depend on

Discussion
Analytical expressions for the energy release rate and mode-mixity phase angle associated with the elementary loads in Figure 1b, Equations (9) and (10), have been derived from the results in Reference [28] for materials with β = 0, and depend on Σ = E 1 /E 2 or α = (E 1 − E 2 )/(E 1 + E 2 ). Useful asymptotics solutions have been derived in Equations (14)- (22) for stiff layers over soft substrates, where α → 1 or E 1 /E 2 → ∞ , e.g., metal coatings on elastomers, and for soft layers over stiff substrates, where α → −1 or E 1 /E 2 → 0 , e.g., polymer coatings over metals; furthermore, intermediate expansions have been derived for near-homogeneous problems.
The effects of β = 0 on the fracture parameters are presented in Tables 3-5. In the absence of shear forces, for V = 0, M, P = 0, the energy release rate, Equation (13), does not depend on β and the effects on the mode mixity phase angle ψ P = ω are synthesized in Figure 6. For stiff layers over soft substrates, where α → 1 or E 1 /E 2 → ∞ , the second Dundurs parameter β = 1/2(1 − 2ν 2 )/(1 − ν 2 ) does not depend on the Poisson ratio of the layer, which can be approximated as a beam [2,52,53]. The effects of β = 0 in this limit are not very pronounced, and this implies that (i) the assumption of incompressible substrate β = 0 is accurate and (ii) the fracture mechanics parameters could be obtained using the approximation of a beam on a Winkler foundation, for M, or the membrane approximation, for P. For soft layers over stiff substrates, where α → −1 or E 1 /E 2 → 0 , the second Dundurs parameter β = 1/2(1 − 2ν 1 )/(1 − ν 1 ) depends on the Poisson ratio of the layer only, and the beam/membrane approximations are not acceptable. The effects of β = 0 are important and controlled by the near-tip deformations in the layer ahead of the crack tip. These effects should not be discounted a priori.  Table 3 and [25]. The angle has been calculated at a reference distance Similar conclusions apply to the mode mixity phase angle associated with pure shear loading 0, , The energy release rate associated with pure shear, Table 4, is also affected by  [2,35,42]. Material systems with 0 β ≠ can be studied using the semi-analytical solutions in Equation (13) and the numerical coefficients in Tables 3-5. They can be applied also to problems with large deformations, e.g., buckle-driven delamination or peeling problems if the crack tip forces are previously calculated using a model which accounts for the geometrical nonlinearities [16,26]. Both analytical and semi-analytical solutions presented in this paper account for the presence of crack tip shear. This contribution is often neglected under the assumption that the effect of the shear becomes negligible for long cracks. This is indeed correct. However, a simple example may clarify the importance of shear in problems where the layer is subjected to a transverse load, e.g., in a blister test, peel test, or SCB face/core interface characterization sandwich testing. Consider the problem in Figure 1 and assume that only a transverse load 0 V ≠ acts at a distance a from the crack tip (similar to a shaft-loaded blister test but with free rotations at the layer end). The transverse load will induce a bending moment, M Va = , and shear, V , at the crack tip, Figure 1b. The energy release rate of Equation (13) then simplifies: The dominant term in G is related to the crack tip bending moment and depends quadratically on the crack length. The term which is directly related to the crack tip shear (the second on the righthand side) is independent of the crack length, and its effect is negligible apart from the case of  Table 3 and [25]. The angle has been calculated at a reference distancer = h 1 .
Similar conclusions apply to the mode mixity phase angle associated with pure shear loading V = 0, M, P = 0. The energy release rate associated with pure shear, Table 5, is also affected by β = 0, even if the effect is weaker. The effects are limited for intermediate values of α and become negligible in the limit α → 1 where the solution does not depend on the Poisson ratios of the substrate and layer and the approximation of a beam over a Winkler foundation is accurate. In the limit α → −1 , neglecting the compressibility of the layer could lead to errors as high as 10%.
The assumption β = 0 should therefore be used with care. This conclusion is supported by experimental results on bi-material systems which show an important variation of the interfacial fracture toughness on varying the mode mixity, especially when crack tip conditions are mode II dominated [2,35,42].
Material systems with β = 0 can be studied using the semi-analytical solutions in Equation (13) and the numerical coefficients in Tables 3-5. They can be applied also to problems with large deformations, e.g., buckle-driven delamination or peeling problems if the crack tip forces are previously calculated using a model which accounts for the geometrical nonlinearities [16,26]. Both analytical and semi-analytical solutions presented in this paper account for the presence of crack tip shear. This contribution is often neglected under the assumption that the effect of the shear becomes negligible for long cracks. This is indeed correct. However, a simple example may clarify the importance of shear in problems where the layer is subjected to a transverse load, e.g., in a blister test, peel test, or SCB face/core interface characterization sandwich testing. Consider the problem in Figure 1 and assume that only a transverse load V = 0 acts at a distance a from the crack tip (similar to a shaft-loaded blister test but with free rotations at the layer end). The transverse load will induce a bending moment, M = Va, and shear, V, at the crack tip, Figure 1b. The energy release rate of Equation (13) then simplifies: The dominant term in G is related to the crack tip bending moment and depends quadratically on the crack length. The term which is directly related to the crack tip shear (the second on the righthand side) is independent of the crack length, and its effect is negligible apart from the case of very short cracks. The effects of shear are important through the third term, which depends linearly on the crack length and has a non-negligible effect also for intermediate crack lengths. The mechanical significance of this term has been first explained in [27]. It describes the effect of crack tip shear on the near-tip deformations generated by the bending moment. These deformations are typically modeled as a crack tip root-rotation, which is the rotation the cross-section experiences at the crack tip (here assumed as positive if clockwise) due to the elastic deformability. If the layer is sufficiently long, the crack tip root rotations are related to the dimensionless crack tip bending moment, (Va)/(E 1 h 1 2 ), through a compliance coefficient, a M 1 . For the problem in Figure 1, the compliance coefficient is ) and the third term on the righthand side of Equation (23) defines the work done by the crack tip shear force on the crack tip root rotation. Equation (23) then shows that, in the presence of shear, the crack tip root rotation due to the crack tip bending moment does contribute to the energy release rate. The effect is important for stiff layers over soft substrates, where the function f V becomes very large (Figure 3b and Table 5). In the absence of shear, the root rotations due to pure bending would not have any effects on the fracture parameters [25], but only on the deformability of the system, e.g., in buckle-driven delamination [15,16,54]. These concepts have recently been discussed in [55] for sandwich beams and in [56] for bi-material DCB specimens.
The solutions derived here have been obtained assuming the crack faces to be traction free and the loads to act at a sufficient distance from the crack tip. However, they can be used to obtain approximate solutions for problems where the loads act as distributed tractions applied along the crack faces, e.g., the blister test or pressure-loaded debonds, using the crack tip force and moment resultants produced by the loads at the crack tip [44,45,57].
The results obtained in this paper can be applied to problems where the geometry has finite dimensions h 2 × 2c and h 1 /h 2 = 0 (with h 2 and 2c depth and width, see Figure A1). This is possible if a, h 2 , and 2c are above minimum values, which must ensure that (i) the crack tip fields are unaffected by the actual distribution of the end loads applied to the thin layer and by the presence of boundaries (free or clamped) and (ii) the stresses generated at the crack tip do not interact with the boundaries. For Saint Venant's principle, the crack length should be just a few times the thickness of the layer h 1 . The analyses in [2] using the model of a beam on a semi-infinite half-plane indicate that the width of the domain should be c/h 1 > 3E 1 /E 2 to ensure that the interfacial tractions ahead of the crack tip drop to almost zero. This would imply c/h > 3000,30000 for E 1 /E 2 = 1000,10000 or α = 0.998,0.999. The calculations performed here and the results in Tables 3-5 show that good engineering accuracy is achieved on the fracture parameters, for −0.99 ≤ α ≤ 0.999, with much smaller widths (c/h = 1000). This conclusion is supported by analyses performed in [8,55,58] on three-layers structures. According to [2], the depth of the substrate should be h 2 /h 1 > 10E 1 /E 2 ; again, the convergence analyses performed here indicate that a smaller depth is necessary to obtain engineering accuracy for large E 1 /E 2 (h 2 /h 1 = 1000). The problem of a corner interface crack, where the crack emerges at the left substrate edge, imposes different restrictions on the crack length a, which are similar to those for c [12,59].
The results and methods in this paper can be extended to orthotropic materials using an orthotropic rescaling technique [43][44][45]60]. For the problem in Figure 1, this extension might be important when dealing with debonds in sandwich structures where the core is often orthotropic. The models applied here to the thin-layer/substrate system can be extended to general bi-material interfaces in beam-like structures. The extension allows for accurate and efficient calculation of the fracture parameters, avoiding case specific numerical calculations, in a large number of fracture tests used for material characterization, e.g., Double-Cantilever Beam tests, mixed-mode bending tests, end-notched flexural tests, four-point bending, and inverted four-point bending tests. Analytical solutions for bi-material beams with h 1 /h 2 = 1 and β = 0 have been recently derived in [61], and semi-analytical solutions for beams with h 1 /h 2 = 1 and β = 0 are presented in [62].

Conclusions
Analytical and semi-analytical solutions have been derived for the energy release rate and mode mixity phase angle of an edge-delamination between a thin layer and an infinitely deep substrate. The 2D solutions are in terms of three elementary crack tip load systems, axial forces, shear forces, and bending moment, and can be used to define the fracture parameters for general loading conditions. Isotropic materials with an exceptionally large mismatch of the elastic constants have been analyzed, and solutions derived for Dundurs' parameters −1 ≤ α ≤ 1 and −0.4 ≤ β ≤ 0.4. The solutions account for the effects of shear, which can be important also for intermediate crack lengths, especially when the near-tip deformations are large as in the presence of soft substrates.
Extreme asymptotics and intermediate expansions have been presented for materials with β = 0, which highlight the scaling of the fracture parameters on varying E 1 /E 2 and provide useful closed-form solutions applicable to stiff layers over soft substrates, to soft layers over stiff substrates and to near homogeneous problems. Results show that for stiff layers over soft substrates ( α → 1 ), the effects of material compressibility, β = 0, are weak; the assumption of substrate incompressibility is accurate as well as the approximation of the layer as a beam on an elastic foundation and also as a beam on a Winkler foundation. For other combinations, including soft layers over stiff substrates ( α → −1 ), the effects may be relevant and problem specific.
A computational model, previously formulated and verified for the derivation of the fracture parameters in bi-material beam-like structures, which uses a crack surface displacement extrapolation technique for mode partitioning, has been verified, and its accuracy proven also for this numerically challenging problem (long slender layer + thick substrate + bi-materials with large mismatch of the elastic constants).