Revisiting the problem of a crack impinging on an interface: a modeling framework for the interaction between the phase field approach for brittle fracture and the interface cohesive zone model

TL;DR

This paper introduces a combined phase field and cohesive zone model to study crack interaction with interfaces, revealing complex phenomena like simultaneous penetration and deflection, and transitions influenced by model parameters.

Contribution

A novel integrated modeling framework that couples phase field brittle fracture with cohesive zone interface modeling, providing new insights into crack-interface interactions.

Findings

The model captures classical fracture mechanics criteria for crack behavior.

Increasing process zone size leads to complex crack phenomena.

Parameter variations induce transitions between crack penetration and deflection.

Abstract

The problem of a crack impinging on an interface has been thoroughly investigated in the last three decades due to its important role in the mechanics and physics of solids. In this investigation, this problem is revisited in view of the recent progresses on the phase field approach to brittle fracture. In this concern, a novel formulation combining the phase field approach for modeling brittle fracture in the bulk and a cohesive zone model for pre-existing adhesive interfaces is herein proposed to investigate the competition between crack penetration and deflection at an interface. The model, implemented within the finite element method framework using a monolithic fully implicit solution strategy, is applied to provide a further insight into the understanding of the role of model parameters on the above competition. In particular, in this study, the role of the fracture toughness ratio between the interface and the adjoining bulks and the characteristic fracture-length scales of the dissipative models are analyzed. In the case of a brittle interface, the asymptotic predictions based on linear elastic fracture mechanics criteria for crack penetration, single deflection or double deflection are fully captured by the present method. Moreover, by increasing the size of the process zone along the interface, or by varying the internal length scale of the phase field model, new complex phenomena are emerging, such as simultaneous crack penetration and deflection and the transition from single crack penetration to deflection and penetration with subsequent branching into the bulk.

Equations

$$\mathbf{u}=\overline{\mathbf{u}}\text{ on }\partial {\Omega }_u\text{ and }\overline{\mathbf{t}}=\mathbf{\sigma }\cdot \mathbf{n}\text{ on }\partial {\Omega }_t,$$

$$\Pi (\mathbf{u},\Gamma )={\Pi }_{\Omega }(\mathbf{u},\Gamma )+{\Pi }_{\Gamma }(\Gamma )={\int }_{\Omega \backslash \Gamma }{\psi }^e(\mathbf{\epsilon })\mathrm{d}\Omega +{\int }_{\Gamma }{\mathcal{G}}_c\mathrm{d}\Gamma ,$$

$${\Pi }_{\Gamma }={\Pi }_{{\Gamma }_b}+{\Pi }_{{\Gamma }_i}={\int }_{{\Gamma }_b}{\mathcal{G}}_c^b(\mathbf{u},\mathfrak{d})\mathrm{d}\Gamma +{\int }_{{\Gamma }_i}{\mathcal{G}}^i(\mathbf{u},\mathfrak{d})\mathrm{d}\Gamma .$$

$${\mathcal{G}}^i=\mathcal{G}(\mathbf{g},\mathfrak{h},\mathfrak{d}).$$

$$\Pi (\mathbf{u},{\Gamma }_b,{\Gamma }_i)={\Pi }_{\Omega }+{\Pi }_{{\Gamma }_b}+{\Pi }_{{\Gamma }_i}={\int }_{\Omega \backslash \Gamma }{\psi }^e(\mathbf{\epsilon })\mathrm{d}\Omega +{\int }_{{\Gamma }_b}{\mathcal{G}}_c^b(\mathbf{u},\mathfrak{d})\mathrm{d}\Gamma +{\int }_{{\Gamma }_i}{\mathcal{G}}^i\left(\mathbf{g},\mathfrak{h},\mathfrak{d}\right)\mathrm{d}\Gamma ,$$

$${\Pi }_b(\mathbf{u},{\Gamma }_b)={\Pi }_{\Omega }(\mathbf{u},{\Gamma }_b)+{\Pi }_{{\Gamma }_b}({\Gamma }_b)={\int }_{\Omega \backslash \Gamma }{\psi }^e(\mathbf{\epsilon })\mathrm{d}\Omega +{\int }_{{\Gamma }_b}{\mathcal{G}}_c^b(\mathbf{u},\mathfrak{d})\mathrm{d}\Gamma .$$

$${\Pi }_b(\mathbf{u},\mathfrak{d})={\int }_{\Omega }\psi (\mathbf{\epsilon },\mathfrak{d})\mathrm{d}\Omega +{\int }_{\Omega }{\mathcal{G}}_c^b\gamma (\mathfrak{d},{\nabla }_{\mathbf{x}}\mathfrak{d})\mathrm{d}\Omega \text{,}$$

$$\hat{\psi }(\mathbf{\epsilon },\mathfrak{d})=\psi (\mathbf{\epsilon },\mathfrak{d})+{\mathcal{G}}_c^b\gamma (\mathfrak{d},{\nabla }_{\mathbf{x}}\mathfrak{d}).$$

$$\gamma (\mathfrak{d},{\nabla }_{\mathbf{x}}\mathfrak{d})=\frac{1}{2l}{\mathfrak{d}}^2+\frac{l}{2}|{\nabla }_{\mathbf{x}}\mathfrak{d}{|}^2$$

$$\mathfrak{d}-l^2{\nabla }_{\mathbf{x}}^2\mathfrak{d}=0\text{ in Ω}\text{ and }{\nabla }_{\mathbf{x}}\mathfrak{d}\cdot \mathbf{n}=0\text{ in ∂Ω},$$

$$\psi (\mathbf{\epsilon },\mathfrak{d})$$

$${\psi }_{+}^e(\mathbf{\epsilon })$$

$${\psi }_{-}^e(\mathbf{\epsilon })$$

$$\mathfrak{g}(\mathfrak{d})={\left(1-\mathfrak{d}\right)}^2+\mathcal{K},$$

$$\mathbf{\sigma }:=\frac{\partial \hat{\psi }}{\partial \mathbf{\epsilon }}=\mathfrak{g}(\mathfrak{d}){\mathbf{\sigma }}_{+}+{\mathbf{\sigma }}_{-};\text{ with }{\mathbf{\sigma }}_{\pm }=\lambda \left(⟨\text{tr}[\mathbf{\epsilon }]{⟩}_{\pm }\right)1+2\mu {\mathbf{\epsilon }}_{\pm },$$

$$\sigma=\left\{\begin{array}[]{ll}k_{n}\dfrac{g_{n}}{g_{nc}},&\hbox{if $0<\dfrac{g_{n}}{g_{nc}}<1$;}\\ 0,&\hbox{if $\dfrac{g_{n}}{g_{nc}}\geq 1$,}\end{array}\right.$$

$${\mathcal{G}}_{IC}=\frac{1}{2}{k}_n{g}_{nc}^2.$$

$$k_n=k_{n,0}{\left(\frac{g_{nc,0}}{g_{nc}}\right)}^2,$$

$${\mathcal{G}}_I^i(\mathfrak{d})=\frac{1}{2}{k}_{n,0}{g}_n^2\frac{g_{nc,0}^2}{{\left[(1-\mathfrak{d})g_{nc,0}+\mathfrak{d}{g}_{nc,1}\right]}^2}.$$

$$\tau=\left\{\begin{array}[]{ll}k_{t}\dfrac{g_{t}}{g_{tc}},&\hbox{if $0<\dfrac{g_{t}}{g_{tc}}<1$;}\\ 0,&\hbox{if $\dfrac{g_{t}}{g_{tc}}\geq 1$.}\end{array}\right.$$

$$k_t=k_{t,0}{\left(\frac{g_{tc,0}}{g_{tc}}\right)}^2.$$

$${\mathcal{G}}_{II}^i(\mathfrak{d})=\frac{1}{2}{k}_{t,0}{g}_t^2\frac{g_{tc,0}^2}{{\left[(1-\mathfrak{d})g_{tc,0}+\mathfrak{d}{g}_{tc,1}\right]}^2}.$$

$${\left(\frac{{\mathcal{G}}_I^i}{{\mathcal{G}}_{IC}^i}\right)}^2+{\left(\frac{{\mathcal{G}}_{II}^i}{{\mathcal{G}}_{IIC}^i}\right)}^2=1,$$

$${\mathcal{G}}_{IC}^i=\frac{1}{2}{g}_{nc,0}^2{k}_{n,0};{\mathcal{G}}_{IIC}^i=\frac{1}{2}{g}_{tc,0}^2{k}_{t,0}.$$

$$\delta {\Pi }_b(\mathbf{u},\delta \mathbf{u},\mathfrak{d},\delta \mathfrak{d})={\int }_{\Omega }\mathbf{\sigma }:\delta \mathbf{\epsilon }\mathrm{d}\Omega -{\int }_{\Omega }2(1-\mathfrak{d})\delta \mathfrak{d}{\psi }_{+}^e(\mathbf{\epsilon })\mathrm{d}\Omega +{\int }_{\Omega }{\mathcal{G}}_c^{b}l\left[\frac{1}{l^2}\mathfrak{d}\delta \mathfrak{d}+{\nabla }_{\mathbf{x}}\mathfrak{d}\cdot {\nabla }_{\mathbf{x}}(\delta \mathfrak{d})\right]\mathrm{d}\Omega +\delta {\Pi }_{b,\text{ext}}(\mathbf{u},\delta \mathbf{u}),$$

$$\delta {\Pi }_{b,\text{ext}}(\mathbf{u},\delta \mathbf{u})={\int }_{\partial \Omega }\overline{\mathbf{t}}\cdot \delta \mathbf{u}\mathrm{d}\partial \Omega +{\int }_{\Omega }{\mathbf{f}}_v\cdot \delta \mathbf{u}\mathrm{d}\Omega .$$

$$\delta {\Pi }_{{\Gamma }_i}(\mathbf{u},\delta \mathbf{u},\mathfrak{d},\delta \mathfrak{d})={\int }_{{\Gamma }_i}\left(\frac{\partial {\mathcal{G}}^i(\mathbf{u},\mathfrak{d})}{\partial \mathbf{u}}\delta \mathbf{u}+\frac{\partial {\mathcal{G}}^i(\mathbf{u},\mathfrak{d})}{\partial \mathfrak{d}}\delta \mathfrak{d}\right)\mathrm{d}\Gamma ,$$

$$\delta {\tilde{\Pi }}_{{\Gamma }_i}^{el}(\mathbf{d},\delta \mathbf{d},\bar{\mathfrak{d}},\delta \bar{\mathfrak{d}})={\int }_{{\Gamma }_i^{el}}\left(\frac{\partial {\mathcal{G}}^i(\mathbf{d},\bar{\mathfrak{d}})}{\partial \mathbf{d}}\delta \mathbf{d}+\frac{\partial {\mathcal{G}}^i(\mathbf{d},\bar{\mathfrak{d}})}{\partial \bar{\mathfrak{d}}}\delta \bar{\mathfrak{d}}\right)\mathrm{d}\Gamma ,$$

$$\mathbf{g}=\mathbf{NLd}={\hat{\mathbf{B}}}_{\mathbf{d}}\mathbf{d},$$

$${\mathbf{g}}_{\text{loc}}\cong \mathbf{Rg}=\mathbf{R}{\hat{\mathbf{B}}}_{\mathbf{b}}\mathbf{d}.$$

$$\mathfrak{d}\cong {\mathbf{N}}_{\mathfrak{d}}{\mathbf{M}}_{\mathfrak{d}}\bar{\mathfrak{d}}={\hat{\mathbf{B}}}_{\mathfrak{d}}\bar{\mathfrak{d}},$$

Full text

††journal: ———-

1 Introduction

The problem of a crack impinging on an interface plays a major role in the performance of many modern multi-component structures, especially involving composite systems. The mechanical responses of such systems strongly depend on the capabilities of an interface to deflect a crack, relying on the competition between deflection and penetration events.

From the theoretical standpoint, it is well established that fracture of bi-material systems is strongly governed by the energy dissipation of adhesives and interfaces. In the last decades, a significant effort has been devoted to the investigation of the fundamental competition between crack penetration into the bulk and deflection along the interface (see Fig.1 for a range of possible problems). From the pioneering work by Zak and Williams [1], it is known that the power of the stress-singularity for a crack meeting an interface between two bonded linear elastic materials is influenced by the elastic mismatch. In this context, He and Hutchinson [2] have demonstrated that the competition between crack penetration and deflection depends on the ratio between the toughness of the interface and that of the bulk, with also the possibility for an asymmetric single-sided deflection. Moreover, again based on linear elastic fracture mechanics arguments, this competition is found to depend also on the angle of the crack which is impinging on the interface. This problem has been further re-examined by various authors again for linear elastic problems [3, 4, 5, 6], and then within the framework of nonlinear fracture mechanics using cohesive zone models [7]. Specifically, focusing on layered composite materials, several studies have deeply investigated the role of the elastic mismatch of the components, see [8, 6, 4] and the references therein given. Mentionable contributions analyzing the singularity of the stress field have been carried out by using asymptotic methodologies for perfectly bonded [9, 10, 11, 12, 13] and cohesive interfaces [15]. The effect of the mismatch in the plastic behaviour of two bonded similar elastic materials on crack-tip shielding and amplification for fracture perpendicular to a bi-material interface has also been experimentally analyzed in [14]. In this setting, current experimental studies [16] have confirmed the linear elastic criterion for deflection and propagation by He and Hutchinson [2], opening also new perspectives of research for wavy interfaces.

Experimental investigations provide more complex physical phenomena with regard to crack propagation in engineering systems, which have not yet been tackled today. One of the most challenging scenarios, reported in [17], regards the simultaneous occurrence of crack penetration and deflection in a bi-material system. In this particular application, the observed crack pattern is not expected according to the existing linear elastic fracture mechanics criteria and, to the best of the authors’ knowledge, no computational model has been able to simulate its occurrence so far. Moreover, Parab and Chen [18] reported a different and very complex crack pattern in a borosilicate glass/borosilicate glass bi-layered system showing an initial crack arrest at the adhesive interface, followed by penetration into the second layer after some delay. Such a penetration can be straight with a single crack, or with multiple branching, depending on the thickness of the epoxy adhesive. At present, this problem is still unsolved due to the fact that most of the current numerical methodologies suffer from serious operative drawbacks for modeling such complex scenarios.

In addition to the previous considerations, the development of numerical methods (especially finite element (FE)-based formulations) to predict fracture onset, propagation and branching in engineering components has been a matter of intensive research during the last decades, to tackle problems that cannot be solved by analytical methods. Most of the extensively used techniques to trigger quasi-brittle and ductile fracture events fall into the following general categories: $(i)$ Continuum Damage Mechanics (CDM) models accounting for a smeared crack representation [19], which in their local version suffer from mesh dependency that has been partially alleviated by using integral-based non-local and gradient enhanced procedures [20, 21, 22, 23, 24]; $(ii)$ extended FE strategies with nodal kinematic enrichment (extended-FEM, X-FEM) that rely on Partition of Unity Methods (PUM) [25, 26, 27] and element enrichment formulations (enhanced-FEM, E-FEM) [28, 29, 30, 31]; $(iii)$ adaptive insertion of cohesive interface elements during the computation or their prior embedding along all the finite element edges [32, 33, 34, 35, 36, 37, 38]; $(iv)$ thick-level set approaches [39, 40]. Although these strategies have been successfully applied to many different fracture mechanics problems, they all present limitations with regard to predicting crack initiation, crack branching, and crack coalescence for multiple fronts.

To overcome these shortcomings, multi-field variational formulations (usually denominated phase field methods), which account for a nonlocal phase variable governed by a Poisson-type partial differential equation to model fracture events, have recently been proposed in the related literature, see the pioneering studies by Francfort and Marigo [41] and by Amor et al. [42]. These approaches share some mathematical and modeling aspects with CDM models but incorporate a non-local formulation. The foundations of phase field approaches for brittle fracture can be traced back to the classical energy-based Griffith criterion [43] through the introduction of a total energy functional that is the sum of the fracture and elastic energy contributions. The minimization of this functional allows triggering crack nucleation, propagation and coalescence in the continuum. In this regard, remarkable contributions are the seminal formulations in [41, 44, 45], whereas the comprehensive treatment of the so-called $\Gamma$-convergence concept has been addressed in [46, 47, 48] and in the references therein given. Quasi-static phase field formulations for brittle fracture have been proposed by Bourdin et al. [44, 49] and the thermodynamically consistent framework has been extensively developed by Miehe and coworkers [50, 51], Kuhn and Müller [52], and Borden et al. [53]. Recent studies have further extended this modeling strategy to shell structures [54], ductile fracture [55, 56], cohesive-based failure [57], dynamic fracture [58, 59], and multi-physics applications [60, 61], to quote some of the most notable contributions.

Alternative numerical procedures to the previous FE-based approaches based on meshfree techniques have been extensively developed in the last decades [62, 63]. As a consequence of their features, these computational strategies offer several appealing aspects, which are especially suitable for modeling initiation and propagation of crack events in solids or any other source of displacement discontinuity [64, 65, 66]. With the aim of exploiting such capabilities, meshfree techniques have been recently combined with variational approaches for modeling fracture in solids based on local maximum entropy approximants [67].

However, at present, the majority of the investigations within the context of the phase field approach to fracture have been devoted to the analysis of continuous and homogeneous bodies [50], and reinforced-composites [60]. A first attempt to model cohesive fracture in the bulk using the phase field approach was proposed in [57] via a suitable modification of the variational formulation to account for the displacement jumps. In that framework, modeling of the displacement discontinuities is found to be a significant complication requiring an additional constraint to be imposed on the auxiliary field that must be constant in the direction orthogonal to the crack. A recent modeling scheme within the context of the phase field formulation, which accounts for both bulk brittle fracture and interfacial damage has been recently proposed in [68, 69]. This alternative approach relies on the definition of a new energy formulation mixing bulk damageable energy and cohesive surface energy, which is activated based on the level set method. This methodology [68, 69], although very promising for stiff interfaces, does not allow for the consideration of pre-existent discontinuities, which is on the other hand the case of adhesive layers.

In the present study, a modeling framework which combines the phase field model for brittle fracture in the bulk and the cohesive zone model for a pre-existing interface is developed. Similarly to [68, 69], the current modeling strategy is based on the definition of a single functional, which accounts for the energy dissipation of the two aforementioned fracture models. However, instead of using the level set method to trace the jump discontinuity at crack faces, a new interface finite element fully compatible with the phase field approach is proposed. Within this formulation, cohesive tractions are computed based on the relative displacements at the interface as in classical interface elements [72, 73, 35]. In this concern, a possible coupling between the interface fracture energy $\mathcal{G}^{i}$ and the non-local damage in the surrounding bulk is postulated via a dependency of the interface stiffness on the average phase field variable evaluated at the interface flanks. This formulation allows to univocally distinguish between the forms of dissipation at the interface and in the bulk, and can also treat complex situations where the amount of damage in the bulk affects the interface response as well, for instance by a degradation of the interface strength.

The proposed formulation, which is implemented in the finite element analysis program FEAP [74] using a monolithic full implicit solution scheme, is then applied to the study of crack propagation at an interface. First, the classical problem of competition between crack penetration and deflection in homogeneous systems at a brittle interface is re-examined, obtaining results in close agreement with linear elastic fracture mechanics predictions. Moreover, these results are extended to the complex case of a cohesive interface where analytical solutions are not available. Finally, the problem of a crack meeting perpendicularly a bi-material (heterogeneous systems) interface is re-examined, providing for the very first time a plausible explanation with regard to the complex branching phenomena observed in previous experimental investigations [17, 18].

2 Fundamental aspects of the proposed interface model compatible with the phase field approach to brittle fracture in the bulk

This section outlines the fundamentals of the proposed consistent interface formulation to be used in combination with the phase field model for brittle fracture in the bulk. Sects. 2.1-2.3 address the theoretical formulation of the current modeling framework, whereas Sects. 2.4 and 2.5 are devoted to the variational formulation and finite element dicretization of the proposed interface model compatible with the phase field approach to brittle fracture in the bulk, respectively.

2.1 Fundamental hypothesis

Let consider an arbitrary body $\Omega\in\mathbb{R}^{n_{dim}}$ in the Euclidean space of dimension $n_{dim}$, in which the existence of an interface111Note that, without loss of generality, the existing interface herewith assumed can be also placed between two adjoining bodies as is usually idealized in cohesive interface formulations. $\Gamma_{i}$ and an evolving internal discontinuity $\Gamma_{b}$ is postulated, see Fig.2(a). The position of a material point is denoted by the vector $\mathbf{x}$ in the global Cartesian frame within the bulk, whereas $\mathbf{x}_{c}$ identifies an arbitrary point of $\Gamma_{i}$. The body forces are denoted by $\mathbf{f}_{v}:\Omega\rightarrow\mathbb{R}^{n_{dim}}$. The boundary of the body is denoted by $\partial\Omega\in\mathbb{R}^{n_{dim}-1}$. Kinematic and traction boundary conditions are prescribed along the disjoining parts $\partial\Omega_{u}\subset\partial\Omega$ and $\partial\Omega_{t}\subset\partial\Omega$, respectively, with $\partial\Omega_{t}\cup\partial\Omega_{u}=\partial\Omega$ and $\partial\Omega_{t}\cap\partial\Omega_{u}=\emptyset$, yielding:

[TABLE]

where $\mathbf{n}$ denotes the outward normal unit vector to the body, and $\boldsymbol{\sigma}$ is the Cauchy stress tensor.

The variational approach to brittle fracture governing the crack nucleation, propagation and branching is set up through the definition of the following free energy functional [50, 58]:

[TABLE]

where $\psi^{e}(\boldsymbol{\varepsilon})$ is the elastic energy density that depends upon the strain field $\boldsymbol{\varepsilon}$, and $\mathcal{G}_{c}$ is the fracture energy. In Eq.(2), the term $\Pi_{\Omega}(\mathbf{u},\Gamma)$ identifies the elastic energy stored in the damaged body, while the energy required to create the crack complying with the Griffith criterion is denoted by $\Pi_{\Gamma}(\Gamma)$.

The central idea of the present formulation regards the split of the fracture energy function into the corresponding counterparts associated with the dissipated energy in the bulk $\Omega$ (governed by the phase field approach of brittle fracture for the prospective discontinuities $\Gamma_{b}$) and along the existing interface ($\Gamma_{i}$) as follows:

[TABLE]

Therefore, while in the bulk the fracture energy $\mathcal{G}_{c}^{b}$ is dissipated according to the Griffith hypothesis [42], at the interface the corresponding fracture energy is released according to a cohesive zone formulation. In particular, in the following we assume that the interface behavior is ruled by a linear cohesive zone model with tension cut-off, though any other cohesive zone model can be easily incorporated into the present framework.

The energy dissipation at the interface is characterized by the fracture energy function $\mathcal{G}^{i}$, which can be related to the displacement discontinuities at the interface, $\mathbf{g}$, a history parameter, $\mathfrak{h}$, as in [57], but also on the phase field degradation variable of the bulk, $\mathfrak{d}$:

[TABLE]

The phase field variable in Eq.(3) has the physical meaning of an internal state damage variable ($\mathfrak{d}\in[0,1]$, where $\mathfrak{d}=0$ represents an intact material, while $\mathfrak{d}=1$ identifies the fully damaged state), and $l$ stands for a regularization parameter related to the smeared crack width (see Fig.2(b) for an illustration of the effect of this regularization length). Thus, when the characteristic regularization parameter (used for the description of the actual width of the smeared crack) tends to zero $(l\to 0)$, then the formulation outlined in Eq.(3) tends to Eq.(2) in the sense of the so-called $\Gamma$-convergence.

Based on the previous modeling assumptions, the functional in Eq.(2) can be recast as:

[TABLE]

where the functional corresponding to the bulk is:

[TABLE]

2.2 Phase field approach for brittle fracture in the bulk

Within the regularized framework of the phase field approach [49, 50], the potential energy of the system is decomposed into two terms:

[TABLE]

where $\psi(\boldsymbol{\varepsilon},\mathfrak{d})$ is the energy density of the bulk for the damaged state, and $\gamma(\mathfrak{d},\nabla_{\mathbf{x}}\mathfrak{d})$ stands for the so-called crack density functional, with $\nabla_{\mathbf{x}}\bullet$ denoting the spatial gradient operator. As a result, the total free energy density of the bulk $\hat{\psi}$ reads:

[TABLE]

According to [50], the functional $\gamma(\mathfrak{d},\nabla_{\mathbf{x}}\mathfrak{d})$, which is a convex function composed by a quadratic term of $\mathfrak{d}$ and another quadratic term involving its gradient, is given by

[TABLE]

and the corresponding Euler equations associated with the phase-field problem take the form

[TABLE]

where $\nabla_{\mathbf{x}}^{2}\mathfrak{d}$ stands for the Laplacian of the phase field variable.

Regarding the energy density in the bulk $\psi(\boldsymbol{\varepsilon},\mathfrak{d})$, the following positive-negative decomposition is assumed [75, 76]:

[TABLE]

where $\lambda$ and $\mu$ are the Lamé constants, $\text{tr}[\bullet]$ denotes the trace operator, and $g(\mathfrak{d})$ is a degradation function that takes the form

[TABLE]

being $\mathcal{K}$ a parameter that defines a residual stiffness to prevent numerical instabilities in the computational implementation, and simultaneously preventing that the resulting system of equations becomes ill-conditioned. In the above equations, the decomposition of the strain tensor into its positive and negative counterparts, $\boldsymbol{\varepsilon}=\boldsymbol{\varepsilon}_{+}+\boldsymbol{\varepsilon}_{-}$, is exploited in order to account for damage under tensile loading only. The spectral decomposition of the positive part of the strain tensor reads $\boldsymbol{\varepsilon}_{+}=\sum\limits_{i=1}^{n_{dim}}\langle\varepsilon^{i}\rangle_{+}\mathbf{n}_{\varepsilon}^{i}\otimes\mathbf{n}_{\varepsilon}^{i}$, where $\varepsilon^{i}$ and $\mathbf{n}_{\varepsilon}^{i}$ identify the eigenvalues and the eigenvectors of the strain tensor and $\langle\bullet\rangle_{+}=(\bullet+|\bullet|)/2$.

Relying on standard arguments [77], the Cauchy stress tensor is defined as:

[TABLE]

where $\mathbf{1}$ denotes the second-order identity tensor.

The irreversibility of the fracture process is guaranteed by means of the incorporation of a penalty term accounting for the history of the local damage variable [50, 87]. The thermodynamic consistency according to the Clausius-Plank inequality of the present formulation has been been comprehensively addressed in [50], and consequently specific details are omitted here for the sake of brevity.

2.3 Cohesive zone model for interface delamination coupled with the phase field

Particularizing the formulation for two-dimensional applications, the interface fracture energy function introduced in Eq.(4) is assumed to be decomposed in the sum of the Mode I and Mode II energy release rates, $\mathcal{G}_{I}$ and $\mathcal{G}_{II}$, based on the considered cohesive zone model. In the present study, without loss of generality, we adopt a linear Mode I cohesive zone model with tension cut-off upon failure, see previous applications in [78, 79, 80]. Moreover, the same traction-separation profile is used for the cohesive zone relation corresponding to Mode II fracture, see Fig.3. The stiffness $k_{n}$ is usually proportional to the ratio between the Young module of the adhesive and its thickness, see [78].

To propose a formulation as general as possible, a dependency of the CZM description on the phase field variable in the surrounding bulk is herein postulated. In this concern, the critical opening displacement can be considered as a function of the phase field variable $\mathfrak{d}$ that triggers fracture events in the adjoining continuum body. In particular, a linear dependency is herein adopted. The critical opening displacement can be reduced or increased depending on the value of $\mathfrak{d}$ which ranges from zero to unity. The former situation (the reduction of the interface stiffness) can be representative of a damaged interface induced by the growth of damage in the adjoining material. The latter scenario, which considers the reduction of $g_{nc}$ by increasing $\mathfrak{d}$, can be related to structured biological interfaces where fibrils are progressively activated by the increase of deformation in the surrounding material (fiber recruitment) [81, 82, 38].

It should be kept in mind that the previous scenarios are defined by considering a constant interface fracture energy with respect to $\mathfrak{d}$. Therefore, there is no modification with respect to the energy dissipation according to the characteristic properties of the interface. Based on this consideration, an increase of $g_{nc}$ by increasing $\mathfrak{d}$ due to the effect of damage in the surrounding bulk implies a simultaneous reduction of the stiffness $k_{n}$ and of the peak traction $\sigma_{c}$. The same behavior is assumed for the Mode II cohesive tractions, see Fig.3.

By stating a linear relation between the Mode I critical opening displacement and the phase field variable $\mathfrak{d}$, the following governing equation for $g_{nc}$ can be defined: $g_{nc}(\mathfrak{d})=(1-\mathfrak{d})g_{nc,0}+\mathfrak{d}g_{nc,1}$, where $g_{nc,0}=g_{nc}(\mathfrak{d}=0)$ and $g_{nc,1}=g_{nc}(\mathfrak{d}=1)$. Hence, for the Mode I cohesive traction we deduce:

[TABLE]

where $\sigma$ denotes the normal traction component of the interface, being $\sigma_{c}$ its corresponding critical value.

The corresponding Mode I interface fracture energy reads

[TABLE]

Note that, by imposing the condition that $\mathcal{G}_{IC}^{i}$ is constant with respect to the phase-field variable $\mathfrak{d}$, the expression for $k_{n}$ is derived by equating the generic value of the interface fracture energy $\mathcal{G}_{IC}^{i}$ to the value corresponding to the absence of damage in the bulk $(\mathfrak{d}=0)$:

[TABLE]

where $k_{n,0}$ is the interface stiffness for $\mathfrak{d}=0$.

Moreover, due to the above constraint, the following closed-form expression for the Mode I energy release rate is deduced:

[TABLE]

The same functional dependencies are herein proposed for the fracture Mode II:

[TABLE]

where $\tau$ identifies the tangential traction component along the interface, whose critical value is $\tau_{c}$, and $g_{t}$ denotes the relative sliding displacement. Its critical value, $g_{tc}$, also obeys $g_{tc}(\mathfrak{d})=(1-\mathfrak{d})g_{tc,0}+\mathfrak{d}g_{tc,1}$ as for Mode I. In order to provide a Mode II interface fracture energy independent of $\mathfrak{d}$, the stiffness $k_{t}$ of the traction-sliding relation has to satisfy the following condition:

[TABLE]

Hence, the Mode II energy release rate reads:

[TABLE]

Finally, to treat Mixed Mode fracture conditions, the use of a standard quadratic criterion is adopted:

[TABLE]

where:

[TABLE]

Nevertheless, it should be pointed out that the present formulation allows also the use of any Mixed Mode fracture criteria available in the literature that are usually specifically tailored based on the technological application.

2.4 Weak form of the variational problem

In this section, the weak forms corresponding to the phase field model for brittle fracture in the bulk and to the cohesive zone model for the interface are derived.

Following a standard Galerkin procedure, the weak form of the coupled displacement and fracture problem in the bulk according to Eq.(7) reads:

[TABLE]

where $\delta\mathbf{u}$ is the vector of the displacement test functions ($\mathfrak{V}^{u}=\left\{\delta\mathbf{u}\,|\,\mathbf{u}=\overline{\mathbf{u}}\text{ on }\partial\Omega_{u},\mathbf{u}\in\mathcal{H}^{1}\right\}$), and $\delta\mathfrak{d}$ stands for the phase field test function ($\mathfrak{V}^{\mathfrak{d}}=\left\{\delta\mathfrak{d}\,|\,\delta\mathfrak{d}=0\text{ on }\Gamma_{b},\mathfrak{d}\in\mathcal{H}^{0}\right\}$). Eq.(23) holds for any trial functions $\delta\mathbf{u}$ and $\delta\mathfrak{d}$. The external contribution to the variation of the bulk functional in Eq.(23) is defined as follows:

[TABLE]

Regarding the interface contribution to the functional of the system corresponding to the term $\Pi_{\Gamma_{i}}$ in Eq.(5), its virtual variation reads:

[TABLE]

where the displacement test functions corresponding to the displacement field and to the phase field variable are defined in close analogy with the formulation for the bulk.

2.5 Finite element formulation

This section details the numerical strategy pursued to solve the simultaneous quasi-static evolution problems for brittle fracture in the bulk and cohesive fracture along the pre-existing interfaces according to the formulation outlined in Section 2. Standard low-order finite elements are used for the spatial discretization, where a monolithic fully coupled solution scheme for the displacement and the phase field nodal variables is considered.

The principal aspects of the finite element discretization for the phase field approach of brittle fracture for the bulk is addressed in A.

With regard to the developed interface, in line with the discretization of the bulk (see A), $\mathbf{d}$ denotes the vector of nodal unknown displacements, and $\bar{\mathfrak{d}}$ stands for the vector of nodal unknown phase field values of the interface element. Accordingly, Eq.(25) can be recast in a discretized form for each interface finite element $\Gamma_{i}^{el}$ ($\Gamma_{i}\sim\bigcup\Gamma_{i}^{el}$):

[TABLE]

where $\mathcal{G}^{i}=\mathcal{G}_{I}^{i}+\mathcal{G}_{II}^{i}$, whose expressions are reported in the previous section.

The gap vector $\mathbf{g}$ at any point inside $\Gamma_{i}^{el}$ is the result of the difference between the displacements of the opposing points at the interface flanks, which is obtained via the interpolation of the nodal displacements $\mathbf{d}$ multiplied by the matrix operator $\mathbf{L}$:

[TABLE]

where $\mathbf{N}$ denotes a matrix collecting the standard Lagrangian shape functions of the element and $\hat{\mathbf{B}}_{\mathbf{d}}=\mathbf{N}\mathbf{L}$ identifies the interface compatibility operator.

To apply the CZM relation, which is expressed in a local reference setting defined by the normal and tangential unit vectors at the interface [35, 38], the global gap vector in Eq.(27) is multiplied by the standard rotation matrix $\mathbf{R}$ for the computation of the gap $\mathbf{g}_{\text{loc}}$ in the local reference system:

[TABLE]

Similarly, the following expressions and operators are introduced to compute the average phase field variable $\mathfrak{d}$ across the interface $\Gamma^{el}_{i}$ at the element level:

[TABLE]

where $\mathbf{M}_{\mathfrak{d}}$ is an average operator and $\hat{\mathbf{B}}_{\mathfrak{d}}=\mathbf{N}_{\mathfrak{d}}\mathbf{M}_{\mathfrak{d}}$ is the compatibility operator corresponding to the phase field.

Accordingly, the discretized weak form reads:

[TABLE]

which leads to the residual vector components:

[TABLE]

Through the consistent linearization of the residual vectors, the tangent operators of the proposed interface finite element for the fully-coupled implicit solution scheme are derived:

[TABLE]

where the tangent constitutive operators of the interface assume the following form for the present CZM traction-separation relation:

[TABLE]

In the expressions above, the terms $\hat{\alpha}$ and $\hat{\beta}$ read:

[TABLE]

Analogously to Eq.(A.9), the coupled system of equations involving the displacement and the phase fields for the interface element takes the form

[TABLE]

3 Competition between penetration and deflection of a crack impinging on an interface with an inclined angle

The first application herein investigated is concerned with the competition between penetration and deflection of a crack impinging on an inclined interface inside a homogeneous system. In particular, let us consider the square domain under plain strain conditions sketched in Fig.4, containing an initial horizontal notch meeting a cohesive interface at the angle $\vartheta$ from the horizontal axis. The dimensions of the system are $L=B=1$ mm and the initial notch length is set equal to $B/2$. The system is subjected to a uniform displacement $\Delta$ applied along its lower and upper boundaries. The Lamé coefficients of the bulk are $\lambda=121.15$ GPa and $\mu=80.77$ GPa, and $l=0.015$ mm as in a similar case study without the interface discussed in [50].

According to linear elastic fracture mechanics, the ratio between the energy release rate for crack deflection along the interface, $\mathcal{G}^{i}$, and the energy release rate for crack penetration into the bulk, $\mathcal{G}^{b}$, depends on the inclination angle $\vartheta$ as addressed in [2]:

[TABLE]

To assess whether the crack either deflects along the interface or propagates into the bulk, the ratio $\dfrac{\mathcal{G}^{i}}{\mathcal{G}^{b}}$ has to be compared with the ratio between the corresponding critical values (fracture toughnesses) of the interface, $\mathcal{G}^{i}_{c}$, and of the bulk, $\mathcal{G}^{b}_{c}$. The condition for crack deflection reads [2, 3]:

[TABLE]

otherwise the crack penetrates into the bulk.

The critical curve separating these two possible scenarios is shown in Fig.5. In the case of $\vartheta=30^{\circ}$, the threshold value of $\mathcal{G}^{i}_{c}/\mathcal{G}^{b}_{c}$ distinguishing between penetration and deflection is approximately equal to $0.87$.

Considering the total force acting on the system, corresponding to the sum of the tractions on the lower or on the upper boundaries, as the representative mechanical response of the specimen subjected to imposed displacements, the following most general functional dependency on the material and geometrical properties can be stated:

[TABLE]

Assuming in the present problem the same properties for Mode I and Mode II CZM relations, i.e., $\sigma_{c}=\tau_{c}$ and $\mathcal{G}^{i}_{c}=\mathcal{G}^{i}_{Ic}=\mathcal{G}^{i}_{IIc}$, the previous functional dependency can be reduced to:

[TABLE]

According to the $\Pi$-theorem of dimensional analysis [83], the following dimensionless representation is derived by selecting $\sigma_{c}$ and $L$ as the physical independent quantities:

[TABLE]

where $\Phi_{0}$ is a dimensionless function.

The first dimensionless number into the previous parentheses can be replaced without any loss of generality by a linear combination of the first two dimensionless numbers, obtaining $\Pi_{1}=\mathcal{G}^{b}_{c}/\mathcal{G}^{i}_{c}$. This yields:

[TABLE]

Moreover, the second and the third dimensionless numbers in Eq.(41) can be replaced by their combination as done in [35]:

[TABLE]

where we recognize that the second dimensionless number $\Pi_{2}\sim l_{\mathrm{CZM}}/L$ is proportional to the ratio between the process zone size along the interface, $l_{\mathrm{CZM}}\sim\left(\mathcal{G}^{i}_{c}E\right)/\sigma_{c}^{2}$, and the sample size, $L$. This number rules the size-scale effects which are typical of nonlinear fracture mechanics in the presence of a cohesive interface [84, 35].

In the case of a very small value of $\Pi_{2}$ $(\Pi_{2}\to 0)$, the interface is very brittle and linear elastic fracture mechanics is expected to be retrieved as a limit scenario. In this situation, the competition between crack deflection and propagation is solely ruled by $\Pi_{1}$ according to the well-known criterion previously recalled in Eq.(37) [2, 3].

To assess this argument, let us consider a set of material parameters leading to $\Pi_{2}=1.25\times 10^{-7}$. In such a case, the present model provides numerical results which are in very good agreement with the analytic linear elastic fracture mechanics predictions reported in [2], see Fig.6 for two values of the dimensionless number $\Pi_{1}$ which is associated to the ratio between the fracture energies of the bulk and the interface. Specifically, for $\Pi_{1}=0.70<0.87$, crack deflection is predicted to occur, while a prevailing crack penetration is estimated for $\Pi_{1}=1.00>0.87$, see the contour plots of the phase field variable in Fig.6 illustrating the numerically predicted crack paths.

For a cohesive interface with a finite process zone size, the competition between crack deflection and penetration is much more complex and cannot be predicted analytically according to linear elastic fracture mechanics. In general, $\Pi_{2}$ is expected to come into play in addition to $\Pi_{1}$. By selecting material parameters yielding to $\Pi_{1}=1.00$, which would correspond to a value leading to crack penetration according to linear elastic fracture mechanics, parametric simulations are performed by varying $\Pi_{2}$ over two orders of magnitude (from $1.25\times 10^{-7}$ to $6.23\times 10^{-5}$). The corresponding numerical predictions are depicted in Fig.7. In these contour plots of the phase field variable it can be seen that the crack penetrates for the lowest value of $\Pi_{2}$. Conversely, a longer deflection path along the interface is predicted by increasing $\Pi_{2}$. For each of these cases, a subsequent branching into the bulk is also observed. This trend is motivated by the increase in the size of the process zone along the cohesive interface. Based on these results it can be observed that the predicted position of the branching point $\eta_{p}=\eta/L_{\rm{int}}$ is an increasing function of $\Pi_{2}$, where $L_{\rm{int}}$ denotes the length of the interface, see Fig.8.

Examining the effect of the interface angle $\vartheta$, in the limit case of a brittle interface $(\Pi_{2}=1.25\times 10^{-7})$, and setting $\Pi_{1}=0.50$, an increase of the interface inclination is expected to promote the transition from deflection to penetration according to linear elastic fracture mechanics reasonings, see Fig.9(a). This theoretical trend is also captured by the present model, see the contour plots of the phase field variable showing the crack path for the cases labeled $A$ $(\vartheta=30^{\circ})$, $B$ $(\vartheta=45^{\circ})$ and $C$ $(\vartheta=60^{\circ})$ in Fig.9, with a progressive reduction of the length of the delamination path before penetration into the bulk by increasing $\vartheta$.

4 The role of the internal fracture-length scales: competition and interplay between the phase field approach of brittle fracture and the cohesive zone model

The main objective of this section concerns providing a thorough explanation with regard to the competition and interplay between the two fracture mechanics models herein considered, namely the phase field approach for brittle fracture to trigger damage in the bulk and the cohesive zone approach to model crack propagation along a pre-existing interface. In particular, the present discussion analyzes the role of the intrinsic fracture-length scales of both methodologies.

With reference to the phase field approach, following [71, 70], the regularizing parameter $l$ can be interpreted as characteristic fracture-length scale of the bulk, which influences the apparent failure stress $\sigma_{\mathrm{PF}}$ of the system. In particular, these authors established the following relationship between such a length scale, $l$, and the material properties of the bulk and $\sigma_{\mathrm{PF}}$ from numerical tests:

[TABLE]

On the other hand, the fracture process zone size of the cohesive zone approach is also affecting the apparent tensile strength of the mechanical system, say $\sigma_{\mathrm{CZM}}$, and it is related to the CZM parameters as pinpointed in previous studies [35, 7]:

[TABLE]

According to dimensional analysis considerations (see the results derived in Sect. 3), the transition between crack deflection and crack penetration, and the competition between both dissipative phenomena, is expected to depend on the value of these internal fracture-length scales.

To quantitatively investigate this issue, a set of parametric simulations has been carried out by considering one of the problems discussed in Sect. 3, namely the competition between penetration and deflection for a crack impinging on an inclined interface at 30*∘* with respect to the horizontal axis in a square specimen with lateral size $L$. In particular, we consider the scenario where the ratio between the fracture toughness corresponding to the bulk and that of the interface is equal to 500, i.e., the interface is much tougher than the bulk (however, note that the conclusions stemming from the current analysis are of general validity). Three different configurations are examined. ($i$) An asymptotic model with a perfectly bonded cohesive interface, where straight crack propagation in the bulk is the only possible failure mode ($l_{\mathrm{CZM}}$ tends to zero and $l$ has a finite value). ($ii$) An asymptotic model with the cohesive interface embedded into an elastic continuum that is characterized by different process zone sizes. In such a case, interface decohesion is the only potential failure mode since the parameter $l$ tends to zero and $l_{\mathrm{CZM}}$ has a finite value. ($iii$) A coupled problem where both phase field and interface cohesive fracture might take place, with the corresponding characteristics fracture-length scales both finite valued.

As far as the first asymptotic model is concerned, the bulk fracture energy, $\mathcal{G}_{c}^{b}$, is set constant and equal to 0.0054 N/mm, while different values of the characteristic length scale $l$ are examined in order to assess the effect of this phase field parameter on the apparent strength of the system. The results of the numerical simulations are shown in Fig.10 in terms of average stress vs. average strain. The average stress $\bar{\sigma}$ has been obtained by computing the sum of the reaction forces acting on the upper boundary, and dividing it by the lateral size $L$ of the specimen and its unit out-of-plane thickness. The average strain $\bar{\epsilon}$ is given by the imposed vertical displacement $\Delta$ divided by the specimen lateral size. The system response is almost linear till brittle crack growth takes place, inducing a post-peak softening branch. Consistently with previous results reported in the literature [71, 70], the apparent strength $\sigma_{\mathrm{PF}}$ evaluated as the maximum of the average stress-strain curves is increasing by reducing the length scale $l$.

Regarding the second asymptotic model, a set of simulations with different maximum cohesive tractions are carried out, setting $\sigma_{c}=\tau_{c}$ for simplicity, whilst the interface fracture toughness $\mathcal{G}_{c}^{i}$ is kept constant and equal to 2.7 N/mm. Under these conditions, the normal and tangential stiffness of the interface are increasing functions of the maximum cohesive tractions and, correspondingly, the process zone size is diminishing. The mechanical response of the system is almost linear until interface crack growth takes place, leading to a very brittle post-peak softening branch, see Fig.11. The apparent strength of the system is affected by the change in $\sigma_{c}$, (as mentioned before, whose increase is reducing the process zone size $l_{\mathrm{CZM}}$), see [85, 84, 35, 7], among others, for a series of fracture mechanics problems involving single material or bi-material systems. Moreover, the apparent stiffness strongly depends on the value of $\sigma_{c}$. This behavior is due to the presence of its cohesive zone with a finite process zone size, which is contributing with a compliance to the system in addition to the compliance of the linear elastic bulk.

Finally, in the coupled problem, model labelled as ($iii$), the fracture-length scale $l$ of the phase field model of fracture is set equal to 0.015 mm and different values of the cohesive maximum tractions, $\sigma_{c}=\tau_{c}$ are explored in order to vary the process zone size $l_{\mathrm{CZM}}$ and investigate the interplay between the two failure modes. For all the cases herein examined, the obtained FE results can be expressed by a relationship of the type $\bar{\sigma}=k\bar{\epsilon}^{\bar{\alpha}}$, with an initial stage almost linear ($\bar{\alpha}=1$) till the onset of softening. To compare the predictions of the asymptotic models and those of the coupled simulations within a single chart, a bi-logarithmic diagram is preferred over a bi-linear one due to the very different average strains experienced by the system in the simulations. In the bi-logarithmic diagram, the stress-strain relation assumes the form $\log\bar{\sigma}=\log k+\bar{\alpha}\bar{\epsilon}$. Therefore, the apparent stiffness of the system can be quantitatively assessed by the value of the intercept of the curves, $\log k$. For visual comparison, the current FE predictions corresponding to the asymptotic model ($i$) are shown in Fig.12 with blue dots (refer to the online version of the article for colors), while the predictions corresponding to the asymptotic model ($ii$) are shown with red line in the same diagram. The predictions corresponding to the coupled problem, model ($iii$), are also superimposed to the same chart with black dots.

The emergent mechanical response of the system with coupling between the phase field approach of brittle fracture and the CZM interface delamination shows a very complex trend which lies in between the two asymptotic models ($i$) and ($ii$), depending on the ratio $l_{\mathrm{CZM}}/l$. The apparent strength of the system, $\bar{\sigma}_{f}$, results to be the minimum between the apparent strength of the model ($i$), $\sigma_{\mathrm{PF}}$, and of the model ($ii$), $\sigma_{\mathrm{CZM}}$, corresponding to the results of the asymptotic models for the same values of the variables $l$ and $l_{\mathrm{CZM}}$, i.e., $\bar{\sigma}_{f}\sim\text{min}\{\sigma_{\mathrm{PF}},\sigma_{\mathrm{CZM}}\}$. The value of $\sigma_{\mathrm{PF}}$ corresponding to $l$=0.015 mm is marked in Fig. 12 on the curve corresponding to the predictions of the asymptotic model ($i$). The curves A, B, C, D and E correspond to $\sigma_{c}$ equal to 10000 MPa, 1000 MPa, 500 MPa, 300 MPa, and 20 MPa, respectively. For these cases A, B, C and D, the corresponding $\sigma_{\mathrm{CZM}}$ is higher than $\sigma_{\mathrm{PF}}$ and therefore the apparent strength is limited by the phase field model. On the other hand, for the case E, the situation is opposite and the coupled model predicts an apparent strength of the system closer to that of the asymptotic model ($ii$) for the same value of $\sigma_{c}$.

Examining in the detail the contour plots of the phase field variable (Figs.13(a)-(c)) and of the corresponding vertical displacements at failure (Figs. 13(d)-(f)) depending on the ratio $l_{\mathrm{CZM}}/l$, the final crack pattern can be visualized. Crack propagation into the bulk (due to the phase field) is eventually prevailing over an initial interface decohesion for $l_{\mathrm{CZM}}/l<1$ (Fig.13(a) corresponding to the case labeled A in Fig.12), with the crack pattern clearly defined by the level set of the phase field variable equal to unity, see the corresponding vertical displacements in Fig.13(d). On the other hand, interface decohesion is the predominant dissipative mechanism for $l_{\mathrm{CZM}}/l>1$ (Fig.13(c) corresponding to the case labeled E in Fig.12), with a phase field variable generally less than unity over the whole domain, and triggering significant interface relative opening displacements (Fig.13(f)). Interestingly, for $l_{\mathrm{CZM}}/l\sim 1$, interface decohesion is predicted to take place and the phase field variable is also reaching unity (Fig.13(b) corresponding to the case labeled C in Fig.12), with a level set $\mathfrak{d}=1$ coincident with the interface trajectory, whose vertical displacements are depicted in Fig.13(e).

Theoretical results of the asymptotic models in Eqs.(43) and (44) suggest that the ratio $l_{\mathrm{CZM}}/l$ is proportional to $\mathcal{G}_{c}^{i}/\mathcal{G}_{c}^{b}(\sigma_{\mathrm{PF}}/\sigma_{c})^{2}$, i.e., to $\Pi_{1}(\sigma_{\mathrm{PF}}/\sigma_{c})^{2}$. Considering $\Pi_{1}=500$ and $\sigma_{\mathrm{PF}}$ given by the asymptotic model ($i$) for $l$=0.015 mm, the value of $\sigma_{max}$ to reach $l_{\mathrm{CZM}}/l\sim 1$ is predicted to be about 300 MPa, which is exactly the value used to obtain the above transitional configuration C. Finally, we notice that for $l_{\mathrm{CZM}}/l\rightarrow 0$, the apparent stiffness of the system asymptotically approaches the value corresponding to that provided by the phase field model of fracture.

5 Competition between penetration and deflection for a crack perpendicular to a bi-material interface

In this section, the well-known problem of a crack perpendicular to a bi-material interface is re-examined according to the current framework. The bi-material specimen is a square domain with lateral side $L=1$ mm and with an edge crack, see Fig.14 for the geometry and the boundary conditions. The system is subjected to uniform tensile loading by applying imposed displacements at the lower and upper sides. The initial edge crack triggers crack propagation in Mode I into material 1, until the crack meets the bimaterial interface. While this stage is mostly governed by the phase field approach to fracture in the bulk (through material 1), the subsequent crack pattern may involve deflection or penetration into material 2 since it is affected by the interface and by the elastic mismatch between materials 1 and 2. Hence, the Dundurs’ parameters $\alpha$ and $\beta$ are introduced to characterize the elastic mismatch of the bi-material system as in [2]:

[TABLE]

where $\mu_{i}$, $\nu_{i}$ $(i=1,2)$ denote the Lamé constant and the Poisson ratio of the two materials under consideration.

He and Hutchinson [2] found that three different mechanisms can take place in linear elasticity when the crack meets the bimaterial interface depending on the number $1/\Pi_{1}$, which is equal to ratio between the fracture toughness of the interface over the corresponding value for the bulk: $(i)$ double deflection for small values of $1/\Pi_{1}$, which corresponds to a situation where material 2 is much tougher than the interface; $(ii)$ single deflection for larger values of $1/\Pi_{1}$; $(iii)$ penetration in the bulk for very large values of $1/\Pi_{1}$. The curves separating these three scenarios are shown in Fig.15(a) for a bi-material system with $\beta=0$ and different $\alpha$. Three material configurations labeled as $A$, $B$ and $C$ have been selected from this diagram as representative for these three situations in order to assess the capability of the proposed numerical method to capture these theoretical trends for a brittle interface $(\Pi_{2}\to 0)$. The corresponding contour plots of the predicted vertical displacement field at failure for such three cases, to highlight the crack pattern and the displacement discontinuities at the interface by the abrupt change of color from blue to red, are shown in Figs.15(b)-15(d). Analyzing the results shown in these graphs it can be seen that the obtained numerical predictions are in very good agreement with theoretical results based on linear elastic fracture mechanics. This again pinpoints the predictive capability of the proposed methodology without any kind of numerical perturbation to capture the single sided deflection.

Moreover, to assess the role played by a not negligible size of the process zone along the cohesive interface, let us focus on a crack perpendicular to the interface with vanishing elastic mismatch $(\alpha=0)$ and setting $\lambda_{1}=\lambda_{2}=121.15$ GPa, $\mu_{1}=\mu_{2}=80.77$ GPa and $\Pi_{1}=1$. According to the diagram shown in Fig.15(a) based on linear elastic fracture mechanics arguments, crack penetration is expected to occur. However, this is valid for a brittle interface having a vanishing process zone size, i.e., for $\Pi_{2}\to 0$. To highlight the role of $\Pi_{2}$ on the competition between delamination and penetration, a constant interface fracture energy $\mathcal{G}_{c}^{i}=0.0054$ N/mm is considered and the maximum peak stress $\sigma_{c}$ of the cohesive zone model is progressively reduced. Based on the definition of $\Pi_{2}$, this aspect implies an increase of this dimensionless number and, consequently, of the size of the process zone along the interface. In turn, this corresponds to a reduction of the stiffness of the interface. The computed dimensionless force-displacement curves for three different values of $\Pi_{2}$ obtained by the application of the current numerical method are shown in Fig.16(a), while the corresponding contour plots of the vertical displacements at failure are displayed in Figs.16(b)-16(d). Analyzing these results, it can be readily observed that, by increasing $\Pi_{2}$, a transition from penetration to delamination is predicted to take place. For $\Pi_{2}=0.125$, the interface has a very high maximum traction and a tiny process zone size. This configuration leads to crack penetration into the bulk and the dimensionless force-displacement response is almost unaffected by the presence of the interface, see the red curve in Fig.16(a). This prediction is in excellent qualitative agreement with previous linear elastic fracture mechanics considerations. Indeed, the contour plots of the phase field variable for different levels of the imposed displacement show that there is no build up of sliding along the interface and failure is due to straight crack growth (see Fig.16(b)). Therefore, crack penetration into the second layer is estimated to be propagated without any delay with respect to the pseudo-time of the quasi-static computation (see the labels in Fig.17 reporting the value of the dimensionless imposed displacement $\Delta/L$ in correspondence of different propagation steps).

For a mid value of $\Pi_{2}=0.498$ (Fig.16(c)), partial interface delamination is predicted to occur. This is observable from the discontinuity in the contour plot representing the dimensionless vertical displacement field in the domain. This discontinuity can be interpreted as a measure of the sliding taking place along the interface. In any case, the final crack pattern is still represented by a subsequent penetration of the crack into the adjacent bulk. This phenomenon of simultaneous delamination and penetration has been experimentally observed in [17], see Fig.18.

For a large value of $\Pi_{2}=12.5$ (Fig.16(d)), the significant drop in the load carrying capacity of the system observed in the dimensionless force-displacement diagram in Fig.16(a) is due to the development of crack deflection along the interface, which is again clearly visible by the discontinuity in the depicted contour plot. By further increasing the imposed displacement $\Delta$, crack penetration is finally predicted to take place again in the bulk after some pseudo-time delay. In such a situation, this pattern is followed by crack branching in the bulk. The corresponding evolution of the phase field variable for all these stages is shown in Fig.19 for different values of $\Delta/L$ to appreciate the delay taking place before crack penetration caused by a consistent build up of sliding at the interface (see the difference between the dimensionless imposed displacement in the labels of Figs.19(a) and 19(b)).

It is remarkable to note that the above numerical predictions (Figs. 13 and 15) provide a mechanical interpretation to the complex crack pattern observed in the experimental results reported in [18] in which the morphology of the crack pattern varies from straight crack penetration in the case of a thin interface to simultaneous delamination followed by penetration and branching into the second material layer, see Fig.20. Moreover, Parab and Chen [18] also reported a very different time delay before the occurrence of crack penetration, depending on the adhesive thickness, see the values given in the caption of Fig.20.

To provide a mechanical interpretation to such complex crack patterns, it is remarkable to note that the previous numerical predictions leading to the crack patterns depicted in Figs.17 and 19 have been obtained through the variation of the peak traction value $\sigma_{c}$ for a constant interface fracture energy. Considering that $k_{n}=\sigma_{c}/g_{nc}=E/t$ for an adhesive [78, 86, 35], where $E$ and $t$ are the Young modulus and thickness of the interface, respectively, it is possible to state that $\Pi_{2}=\mathcal{G}^{i}_{c}E/(\sigma_{c}^{2}L)=0.5t/L$ for the present cohesive zone model. Hence, the larger the value of $\Pi_{2}$, the thicker the interface for this case problem.

In this concern, consistently with the experimental trends [18], in the present simulations crack penetration is observed for small values of $\Pi_{2}$ (thin brittle adhesive), without appreciable accumulation of sliding along the interface (Fig. 13). On the other hand, for large values of $\Pi_{2}$ (thick adhesive with a larger process zone size), significant sliding takes place along the interface before crack penetration and the subsequent branching (Fig. 15). Furthermore, the present model is also able to explain the reason for the increasing time delay before penetration by increasing the interface thickness observed in the experiments ($\Delta\tau=8.2$ $\mu$s for $t=0.2$ mm; 60.4 $\mu$s for $t=1.0$ mm; 83.9 $\mu$s for $t=2.7$ mm). According to the present predictions, it is possible to state that such a time delay is due to the accumulation of sliding along the adhesive interface, which is a phenomenon that is promoted by configurations with large values of $\Pi_{2}$. Although crack opening data vs. time are not provided in [18] and therefore it is not possible to translate the time delay in the displacement delay, a qualitative satisfactory agreement is achieved by the present model. Thus, it is worth noting that, for a thin interface, crack penetration is predicted to occur as soon as the crack meets the interface (see the values of the dimensionless imposed displacements in the captions of Figs. 17 before and after penetration). Conversely, a significant increase in the applied remote displacements is required to trigger crack penetration for a thicker interface. This latter event can stem from the result of the energy dissipation provoked by interface sliding, see the analogous captions in Figs. 19.

6 Conclusion

In this study, a novel modeling framework which combines the phase field approach for brittle fracture and a cohesive formulation for a pre-existing interface has been proposed.

Several benchmark test problems inspired by well-known linear elastic fracture mechanics methodologies regarding a crack impinging on an interface with or without elastic mismatch between the two jointed layers have been fully retrieved by the present approach in the case of a quasi-brittle interface with a small process zone size. Further numerical predictions for scenarios in which the interface has a finite process zone size have been carried out to characterize configurations that cannot be analyzed using analytical methods based on linear elastic fracture mechanics.

The proposed methodology has been also validated by the comparison with relevant experimental results showing very complex crack patterns that were unexplained so far in the related literature. In particular, a thorough analysis of the effect of the adhesive interface thickness on the crack pattern has been conducted. Results provide a clear mechanical explanation to the time delay experimentally observed for thick adhesives caused by the pile up of sliding along the interface and the transition from straight penetration to simultaneous delamination and penetration followed by branching by increasing the adhesive thickness.

Additionally, a comprehensive analysis with regard to the competition and interplay between both fracture mechanics models relying on their internal length scales has been conducted. In particular, such a competition and interplay is dependent on the ratio $l_{\mathrm{CMZ}}/l$ between the characteristic length scales of both methodologies. This ratio determined the apparent strength of the specimen under analysis, $\bar{\sigma}_{f}$, being possible to state that $\bar{\sigma}_{f}\sim\min\{\sigma_{\mathrm{PF}},\sigma_{\mathrm{CZM}}\}$ when both fracture-length scales are finite valued.

The developed strategy constitutes a very promising simulation tool that allows modeling complex crack patterns in a wide range of engineering systems involving interfaces.

Future developments might regard further analysis on the interaction between bulk damage and interface degradation or strengthening. The latter possibility can be very important for biological applications in which fibre recruitment phenomena take place by increasing the deformation level. Hence, this consideration is motivating further research within the framework of finite elasticity and dynamics.

Acknowledgments

This article is dedicated to the memory of Professor Christian Miehe (University of Stuttgart) for his eminent contributions to the phase field approach to fracture that have inspired the present authors in their research. MP and JR would like to thank the European Research Council for supporting the ERC Starting Grant “Multi-field and multi-scale Computational Approach to Design and Durability of PhotoVoltaic Modules” - CA2PVM, under the European Union’s Seventh Framework Programme (FP/2007-2013) / ERC Grant Agreement n. 306622. JR acknowledges the support of the Spanish Ministry of Economy and Competitiveness (DPI2012-37187) and the Andalusian Government (Project of Excellence No. TEP-7093).

Appendix A Finite element formulation and implementation aspects of the phase field model of brittle fracture for the bulk

This appendix outlines the finite element formulation for the phase field model of brittle fracture for the bulk. Complying with the isoparametric concept, standard Lagrangian shape functions $N^{I}(\boldsymbol{\xi})$, which are defined in the parametric space $\boldsymbol{\xi}=\{\xi^{1},\xi^{2}\}$, are used for the interpolation of the geometry ($\mathbf{x}$), the displacement field ($\mathbf{u}$), its variation ($\delta\mathbf{u}$) and its linearization ($\Delta\mathbf{u}$):

[TABLE]

where $n$ identifies the number of nodes at the element level, and $\mathbf{x}_{I}$ and $\mathbf{d}_{I}$ denote the discrete nodal coordinates and displacements values, respectively, which are collected in the corresponding vectors $\widetilde{\mathbf{x}}$ and $\mathbf{d}$. The interpolation functions are arranged in the operator $\mathbf{N}$ as usual. The strain field ($\boldsymbol{\varepsilon}$), its variation ($\delta\boldsymbol{\varepsilon}$) and its linearization ($\Delta\boldsymbol{\varepsilon}$) are interpolated through the displacement-strain $\mathbf{B}_{\mathbf{d}}$ operator as follows:

[TABLE]

The phase field variable interpolation, its variation and linearization read:

[TABLE]

where $\mathfrak{\overline{d}}_{I}$ stands for the nodal phase field values, which are collected in the vector $\mathfrak{\overline{d}}$. Note that the same shape functions $\mathbf{N}$ are considered for the interpolation of the kinematics and of the phase field variable.

The gradient of the phase field ($\nabla_{\mathbf{x}}\mathfrak{d}$), its variation ($\nabla_{\mathbf{x}}\delta\mathfrak{d}$) and linearization ($\nabla_{\mathbf{x}}\Delta\mathfrak{d}$) are interpolated via the $\mathbf{B}_{\mathfrak{d}}$ operator:

[TABLE]

With the previous interpolation scheme at hand, the discretized version of Eq.(23) at the element level (denoted by the superscript $el$) reads:

[TABLE]

where

[TABLE]

[TABLE]

[TABLE]

where $\mathbf{f}_{\mathbf{d},\text{int}}^{b}$ and $\mathbf{f}_{\mathbf{d},\text{ext}}^{b}$ denote the internal and external residual vectors of the displacement field with $\mathbf{f}_{\mathbf{d}}^{b}=\mathbf{f}_{\mathbf{d},\text{ext}}^{b}-\mathbf{f}_{\mathbf{d},\text{int}}^{b}$, and $\mathbf{f}_{\mathfrak{d}}^{b}$ is the residual vector associated with the phase field variable.

Due to the strong nonlinearities involved into the proposed modeling framework which combines fracture in the bulk and along the interfaces, a fully coupled monolithic solution scheme is herein preferred over staggered methods used in [68, 69, 56] with the aim of preventing numerical instabilities.

The consistent linearization of the resulting nonlinear system of equations yields to the definition of the following coupled system:

[TABLE]

The specific form of the element stiffness matrices $\mathbf{K}_{\mathbf{d}\mathbf{d}}^{b}$, $\mathbf{K}_{\mathbf{d}\mathfrak{d}}^{b}$, $\mathbf{K}_{\mathfrak{d}d}^{b}$ and $\mathbf{K}_{\mathfrak{d}\mathfrak{d}}^{b}$ are omitted here for the sake of brevity. The reader is referred to [87] for further details on the derivation.