Modeling complex crack paths in ceramic laminates: A novel variational framework combining the phase field method of fracture and the cohesive zone model

Contents lists available at ScienceDirect Journal of the European Ceramic Society journal homepage: www.elsevier.com/locate/jeurceramsoc Original Article Modeling complex crack paths in ceramic laminates: A novel variational framework combining the phase field method of fracture and the cohesive zone model V. Carollo a,⁎ , J. Reinoso b,a , M. Paggi a a IMT School for Advanced Studies Lucca, Piazza San Francesco 19, 55100 Lucca, Italy b Elasticity and Strength of Materials Group, School of Engineering, Universidad de Sevilla, Camino de los Descubrimientos s/n, 41092 Seville, Spain ARTICLE INFO Keywords: Phase field model of fracture Cohesive zone model Crack deflection Crack branching Laminates ABSTRACT The competition between crack penetration in the layers and cohesive delamination along interfaces is herein investigated in reference to laminate ceramics, with special attention to the occurrence of crack deflection and crack branching. These phenomena are simulated according to a recent variational approach coupling the phase field model for brittle fracture in the laminae and the cohesive zone model for quasi-brittle interfaces. It is shown that the proposed variational approach is particularly suitable for the prediction of complex crack paths involving crack branching, crack deflection and cohesive delamination. The effect of different interface properties on the predicted crack path tortuosity is investigated and the ability of the method to simulate fracture in layered ceramics is proven in relation to experimental data taken from the literature. 1. Introduction Ceramic materials are largely used in technological applications, especially with the aim of achieving a desired resistance to severe wear or corrosive phenomena at high temperatures. However, the main drawback of ceramics regards their brittleness and, to increase their toughness, laminates are often used alternating ceramic and metallic layers. For instance, in [1,2], Al/SiC and Al/TiN laminates have been explored and tested. The metallic layers make the composite able to withstand higher deformations by means of the development of plasticity at several locations within the specimen, and therefore increasing the overall toughness of the laminate. The same toughening process has been achieved by alternating ceramic layers with polymeric layers in [3]. The main drawback of these solutions is that metals and polymers loose their mechanical properties at high temperatures and have a low wear resistance. A possible way to enhance the toughness of ceramics is to introduce quasi-brittle interfaces [4,5]. Then, a stack of ceramic layers alternated by thin layers of a very brittle ceramic is a possible effective technological solution. Such brittle layers act as a quasi-brittle interfaces which make the crack path very complex, thus increasing the overall material toughness by acting on the crack tortuosity. This mechanism has been firstly theorized in the framework of linear elastic fracture mechanics (LEFM) by the so called Cook–Gordon mechanism [6]. It is the result of crack branching and crack deflection typical of the competition between crack penetration in the layers and delamination along the existing interfaces. Another approach to foster complex crack paths is to introduce porous layers between the ceramic ones [7]. However, the drawback of porous materials is their low wear resistance. Alternatively, tough interfaces with preexisting defects could be inserted among the layers [8]. The tough interface can be made of a ceramic material which guarantees the resistance to wear and corrosion, also at high temperatures. The defects in the tough interfaces guarantee the development of crack deflection with a consequent toughening of the material. Considering the technological strategies herein described, the introduction of an interface with tailored properties clearly emerges as a strategy to create a complex crack pattern and consequently enhance the apparent material toughness. In this article, we investigate these possibilities by examining and modelling the interaction between crack penetration in the layers and cohesive delamination at the interfaces between materials with different elastic and fracture properties. This problem has been mathematically addressed in many publications using different methodologies, cf. [9–11]. Differing from previous investigations, this problem is herein analysed by means of the numerical method published in [12], which has been recently extended to a 3D finite deformation framework in [13] and applied to fracture of anisotropic polycrystalline Silicon in [14]. The novelty of the current approach relies upon the innovative variational framework combining the https://doi.org/10.1016/j.jeurceramsoc.2018.01.035 Received 30 November 2017; Received in revised form 21 January 2018; Accepted 22 January 2018 ⁎ Corresponding author. E-mail address: [email protected] (V. Carollo). Journal of the European Ceramic Society 38 (2018) 2994–3003 Available online 02 February 2018 0955-2219/ © 2018 The Author(s). Published by Elsevier Ltd. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/BY-NC-ND/4.0/). T phase field method of fracture and the cohesive zone model. In particular, the role played by the internal characteristic length scales of the two approaches is rigorously analysed in order to understand their effect on the resulting crack path and its tortuosity as a way to enhance the overall composite toughness. The manuscript is organized as follows. Section 2summarizes the principal features of the variational framework herein employed. In Section 3, the numerical method is applied to predict the crack path in laminates. Special attention is devoted to examining the effect of tough and quasi-brittle interfaces and reproducing also experimental results related to a ceramic laminate taken from literature. Finally, the main conclusions of the current investigation are drawn in Section 4. 2. Variational model In this section, the coupled phase field and cohesive zone model formulation developed in [12] is resumed. We first present the fundamental hypothesis of the current coupling approach in Section 2.1. Later, the phase field method for brittle fracture formulation and the cohesive zone model compatible with the phase field model are outlined in Sections 2.2 and 2.3, respectively. Finally, the corresponding finite element formulation within the infinitesimal deformation setting is derived in Section 2.4. 2.1. Fundamental hypothesis The formulation herein presented is developed in the general Euclidean space of dimension n dim under infinitesimal deformation setting. Let us consider a body ∈Ωℝ ndi m with a generic shape, where the boundaries of the body are denoted by ∂ ∈ − Ωℝ n 1 dim (Fig. 1). Kinematic and traction boundary conditions can be respectively prescribed on the disjointed parts of the boundaries ∂Ω u and ∂Ω t (with ∂Ω t ∪∂Ω u =∂Ωand ∂Ω t ∩∂Ω u =∅). Then, the prescribed displacements and tractions are denoted by: =∂ = ∂ u utσnon Ω and · on Ω , ut (1) where nis the outward normal unit vector to the body, and σis the Cauchy stress tensor. The body forces are represented by the function → f :Ω ℝ vndi m . The composite is characterized by quasi-brittle interfaces, Γ i , and a crack in the layer represented as an internal discontinuity, Γ b (see Fig. 1(a)). A generic point in the bulk of the body is denoted by the vector of its Cartesian coordinates x, while a generic point on the interface Γ i is denoted by the vector x c . The free energy functional which governs the mechanics of the body Ωis defined as [15,16]: ∫∫ =+= + ∖ψuu ε Π (,Γ) Π(,Γ) Π(Γ) ()dΩ dΓ, ecΩΓ ΩΓ Γ G (2) where ψ e (ε) is the elastic energy density, εis the strain field, and c G is the fracture energy. The main idea to couple the phase field approach for brittle fracture and the cohesive zone model is to split the fracture energy function c G in two parts. One part ( c b G ) describes fracture in the layers and it is modelled by the phase field approach. The second part ( i G ) describes the cohesive fracture of the interfaces and it is modelled by the cohesive zone approach. Then, the free energy functional in Eq. (2) can be rewritten as: ∫∫ ∫ =++= + +∖ψuεu g Π( , Γ , Γ) Π Π Π ( )dΩ ( , )d Γ (,, )dΓ, bi ec b i ΩΓ ΓΩΓ Γ Γ bib i G G d hd (3) where gdenotes the vector of displacement discontinuities at the interface, h is an history parameter as in [17] to avoid re-healing of the material and the non uniqueness of the solution, and d is the phase field degradation variable which will be further detailed in the next section. 2.2. Phase field approach for brittle fracture The phase field approach for brittle fracture [15,18] is a variational approach which considers a crack as a diffuse damage instead of a sharp discontinuity (Fig. 1(a)). Within this framework, the potential energy of the bulk is formulated as follows: ∫∫ =+∇ψγuε Π (, ) (, )dΩ (, )dΩ, bc bx ΩΩ G dd dd (4) Fig. 1. Schematic representation of an arbitrary body with a discontinuity in the domain and an interface: (a) Left: discrete discontinuity in the domain. Right: smeared discontinuity in the domain based on the phase field concept. (b) Diffusive crack modeling solution for the one-dimensional crack problem. V. Carollo et al. Journal of the European Ceramic Society 38 (2018) 2994–3003 2995 where ψ ε(, ) dis the elastic energy stored in the bulk, the symbol ∇ x • denotes the spatial gradient operator, and ∇γ(, ) x dd is the so called crack density functional which reads [15]: ∇= +∇γl l (, ) 1 22 || xx 22 dd d d (5) where lis the phase field internal length parameter governing the sharpness of the crack according to the equation presented in Fig. 1(b) for the mono-dimensional case. The elastic energy in Eq. (4) takes the following form: =+ +− ψ ψψεεε(, ) () () () , ee dgd (6a) =+ +++ ψ λμεε ε() 2(tr[] ) tr[ ] , e22 (6b) =+ −−− ψ λμεε ε() 2(tr[] ) tr[ ] , e22 (6c) where λand μare the Lamé constants, tr[•] denotes the trace operator, ϵ + and ϵ − denote the positive and negative counterpart of the strain tensor, respectively, and g()dis a degradation function: =− +() (1 ) . 2 K gd d (7) In Eq. (6a), the elastic energy has been split in its positive and negative counterparts according to the formulation in [19,20]. The positive counterpart of the elastic energy is produced by the tensile stresses while the negative counterpart is produced by compression. Positive and negative stresses can be computed by the derivative of the elastic energy with respect to the strain tensor. Then, the introduction of the positive and negative split of the strain tensor leads to: =∂ ∂=+ = + +− ± ±± ψλμσεσσ σ ε1 ε:ˆ() ; with (tr[] ) 2 , gd (8) Then, since the degradation function in Eq. (7) affects only the positive counterpart, damage can only develop when the material is under tension, avoiding damage growth in compression. 2.3. Cohesive zone model compatible with the phase field In this section, the classical linear cohesive zone model with tension cut-off[21] is particularized in order to take into account the effect of the bulk damage d . First of all, the cohesive counterpart of the fracture energy in Eq. (3) is decomposed in the sum of the Mode I and Mode II fracture energies, I G and II G , respectively. Based on the formulation outlined in [12], the critical crack opening displacement (g c ) depends on the bulk damage d according to the linear relation =− +ggg() (1 ) ccc,0 , 1 ddd , where ==gg(0 ) cc,0 d and ==gg(1 ) cc,1 d. Then, the cohesive traction vs. relative displacement laws for Mode I and Mode II take the form shown in Fig. 2, and are described by the following equations: =⎧ ⎨ ⎩ << ≥=⎧ ⎨ ⎩ << ≥ σ k τ k,if0 1; 0, if 1, ,if0 1; 0, if 1. n g g g g g g t g g g g g g nn n tt t nc nc nc tc tc tc (9) where σand τare the Mode I and Mode II tractions, respectively, gis the relative displacement, and the subscript nand trefers to opening and sliding, respectively. The stiffness of the cohesive relation, k, depends on damage d according to the formulae: == () () kk kk,. nn g gtt g g ,0 2 ,0 2 nc,0 nc tc,0 tc (10) where k 0 and g 0 are, respectively, the stiffness and critical relative displacements at =0 d . In this formulation, among the different potential alternatives, we have adopted the hypothesis of keeping constant the critical energy release rate of the interface, c i G , for many reasons. Firstly, the model is formulated starting from the Griffith energy balance criterion, then in Eq.(3) we set up a clear and explicit split between the dissipated energy due to the bulk fracture and to the interface delamination. Correspondingly, c i G becomes the material parameter that governs the interface failure. This choice endows a clear characterization of the interface fracture energy. In fact, among the four parameters entering the cohesive law, c i G can be obtained in a straightforward manner from experimental tests [22]. Finally, the mixed mode failure criterion proposed in [23] is herewith considered to trigger the interface failure: ⎜ ⎟⎜⎟ ⎛ ⎝ ⎞ ⎠+⎛ ⎝⎞ ⎠=1 , I i i i i IC 2 II IIC 2 G G G G (11) where I i G and i II G are the dissipated fracture energies which take the form: == −+ −+ ng kg() , () . I itn g gg itt g gg 1 2,0 2 [(1 ) ] II 1 2,0 2 [(1 ) ] nc,0 2 nc,0 nc,1 2 tc,0 2 tc,0 tc,1 2 GG dd dd dd (12) The critical fracture energies i IC G and i IIC G are: ==gk gk,. init IC 1 2nc,0 2,0 IIC 1 2tc,0 2,0 GG (13) 2.4. Finite element formulation The finite element formulation of the previous fracture mechanics models is herein derived. First of all, the weak form of the free energy functional in Eq. (3) is deduced using the standard Galerkin procedure. Then, the variation of the bulk energy functional (Eq. (4)) with respect to the displacements uand the phase field variable d takes the form: Fig. 2. Schematic representation of the cohesive zone model coupled with the phase field variable for brittle fracture in the bulk. (a) Mode I CZM traction σvs. g n . (b) Mode II CZM traction τvs. g t . V. Carollo et al. Journal of the European Ceramic Society 38 (2018) 2994–3003 2996 ∫∫ ∫ =−− + ⎡ ⎣+∇ ∇ ⎤ ⎦ +∀ + ψ δδ uu σε ε uu u δΠ ( , δ , d, δd) :δ dΩ 2(1 d)δd ()dΩ l·()dΩ δΠ ( , δ ) δ , δd, be c b l b xx ΩΩ Ω 1 ,ext 2 G dd d d (14) where ∈= = ∂ ∈δδuuuuu{| onΩ, } uu1 H V is the vector of the displacement test function, and ∈= = ∈δδδ{| 0onΓ, } b0 H dV dd d dis the damage test function. The contribution of the external forces in the variation of the bulk energy functional is defined as follow: ∫∫ =∂+ ∀ ∂ δδ δ δδδuu tu fu uΠ(,) ·dΩ ·dΩ , . bv,ext ΩΩ d (15) Finally, the variation of the interface energy functional Π Γ iin Eq. (3) is defined as: ∫ ⎜⎟ =⎛ ⎝∂∂+∂∂⎞ ⎠∀δδδ δ δ δδuu u uuuuΠ(, , , ) (, ) (, ) dΓ , . ii ΓΓ ii GG dd dd ddd (16) The phase field model has been implemented within a 4-node isoparametric finite element. The cohesive zone model has been implemented using a 4-node interface finite element. The detailed finite element implementation in the software FEAP [24] and all the related operators can be found in [12] and are omitted here for the sake of brevity. 3. Simulation of complex crack paths in layered ceramics The prediction of the crack path in ceramic laminates is a very challenging issue. Various factors such as the elastic mismatch of the constituent materials and the properties of the interfaces [23,25] make the crack propagating along tortuous paths which are very difficult to be simulated using previous numerical methods. The present approach resolves most of their drawbacks primarily associated to remeshing or corner-case problems. 3.1. Complex crack paths for laminates under tensile loading: analysis of failure patterns with quasi-brittle or tough interfaces In this section we study the effect of the interface toughness on the resulting crack path. The case of a single-edge notched bi-material laminate under tension (Fig. 3) is herein considered. The laminate has been modelled using the phase field finite elements for the bulk and the interface finite elements compatible with the phase field between each layer. The materials which compose the laminate are: a soft material 1 with high fracture toughness; a stiffmaterial 2 with low fracture toughness. The material parameters are collected in Table 1. Three cases have been examined: (1) laminate with perfectly bonded layers; (2) laminate with tough interfaces; (3) laminate with brittle interfaces. Note that the toughness of the interfaces is always larger than the toughness of the materials composing the laminae, to simulate configurations consistent with technological solutions for ceramics laminates. In the first simulation, no interface elements are introduced in order to simulate fully bonded layers. The phase field internal length lis set very small for both layers (see Table 1), to reproduce LEFM predictions as discussed in [26]. The results of the simulations are shown in Fig. 4. The first layer that shows crack nucleation is the second one (material 2). In such a layer, two parallel cracks are predicted to propagate simultaneously (Fig. 4(a)) and can be considered as two branches of the initial notch in the layer 1. After increasing the applied load, each crack in the second layer further branches in the next brittle layer (Fig. 4(b)). The same process continues for the next brittle layer but only two branches are now developed (Fig. 4(c)). At this stage, the cracks in the brittle layers start connecting through the material 1. Finally, failure of the specimen is achieved (Fig. 4(d)) and it is the result of a complex crack path mostly localized in the mid-span cross-section of the specimen. In the second simulation, we introduce interface elements between each layer. The parameters used for the interface are σ c,0 =τ c,0 = 100 MPa and k 0 = 2000 MPa/mm to model a stiffquasibrittle interface. As expected, the evolution of the predicted crack path is very different from that of the previous simulation. The first layer where cracks nucleate is the second one (Fig. 5(a)) with the same pattern as in the first simulation. Immediately after that, delamination along the interface between the first and the second layer takes place (Fig. 5(b)). Then, the first layer is cracked by the propagation from the Fig. 3. Specimen geometry. Table 1 Geometry and material/interface parameters. Geometry parameters L6 mm Specimen length l y 0.25 mm Layer thickness h0.005 mm Interface thickness Mechanical parameters material 1 E 1 70,000 MPa Material 1 Young modulus v 1 0.34 Material 1 Poisson ratio 1 G 0.025 N/mm Material 1 fracture energy l 1 0.0075 mm Material 1 phase field length scale parameter Mechanical parameters material 2 E 2 300,000 MPa Material 2 Young modulus v 2 0.14 Material 2 Poisson ratio 2 G 0.005 N/mm Material 2 fracture energy l 2 0.0075 mm Material 2 phase field length scale parameter Mechanical parameters interface k 0 2000 MPa/mm Initial stiffness of interface σ c,0 ,τ c,0 100 MPa Initial peak stress of the tough interface c i G 2.5 N/mm Critical energy release rate for tough interface σ c,0 ,τ c,0 1 MPa Initial peak stress of the brittle interface c i G 0.025 N/mm Critical energy release rate for brittle interface σ c,0 /σ c ,τ c,0 / τ c 1 Ratio between the initial and final value of the peak stress of the interface V. Carollo et al. Journal of the European Ceramic Society 38 (2018) 2994–3003 2997 notch (Fig. 5(c)). Continuing with the simulation, cracking proceeds in the material 2 layers together with the development of delamination at interfaces (Fig. 5(d)). At failure, delamination makes the crack pattern distributed along the whole specimen, as a primary difference from the results of the first simulation. Another important aspect is that the majority of the material 1 layers are not cracked, apart from the first layer containing the notch. In the third simulation we introduce a more brittle interface, setting σ c,0 =τ c,0 = 1 MPa and k 0 = 2000 MPa/mm. The evolution of the crack path is again quite different from the previous cases (Fig. 6). Thus, in the current case, first, delamination is predicted to occur between the first and the second layer (Fig. 6(a)). Then, the crack starts propagating from the notch until it impinges onto the delaminated interface (Fig. 6(b)). Subsequently, branching is predicted to take place in the second layer (Fig. 6(c)), and each branched crack is developed starting from the points where delamination was arrested. Furthermore, crack Fig. 4. Crack evolution in the simulation with fully bonded layers. Fig. 5. Crack evolution in the simulation with tough interface (σ c,0 = 100 MPa). V. Carollo et al. Journal of the European Ceramic Society 38 (2018) 2994–3003 2998 penetration triggers delamination along the second interface. Finally, the sample failure is the result of a sudden delamination along all the interfaces and cracks in the brittle layers (Fig. 6(d)). Again, the material 1 layers have been preserved from cracking. The force–displacement curves of the above three simulations are shown in Fig. 7. All of them present an initial nonlinearity due to the formation of damage in the bulk. After that, the curve of the first simulation starts loosing the load-bearing capacity slowly with a smooth softening. This is due to the gradual crack propagation within the specimen. The curves corresponding to the second and third simulations show a completely different crack pattern and the nonlinear effects are much more pronounced in the load–displacement curve due to the occurrence of delamination. The development of cracking and delamination events cause multiple drops in the resistant force. The correspondence between crack/delamination events and the drops in the force can be examined closely in Fig. 8. In these graphs, the force–displacement curves are plotted together with other two quantities, the total crack propagation length and the portion of crack propagation length in just the bulk. Both quantities are normalized with respect to the layer thickness, l y . These plots show that, in general, delamination and cracking occur simultaneously. This is in agreement with the patterns in Figs. 5 and 6, where in some cases delamination triggers crack formation in the layers or vice-versa. Another important consideration is that the delamination length in the simulation with quasi-brittle interfaces is bigger than the in simulation with tougher interfaces. Moreover, we compare the current predictions with the theoretical results presented in [6], which complied with the so-called Cook–- Gordon mechanism. According to that theory, the introduction of a brittle interface should increase the apparent material strength. The results reported in Fig. 7 seem to challenge this theory, since a brittle interface reduces the apparent material strength as compared a tougher one. This contradiction can be explained by noting that the assumptions of the Cook–Gordon model are not fully satisfied in the present computational setting. One of the most relevant differences regards the fact that in the Cook–Gordon model the material is considered homogeneous and linear elastic until the condition of fracture of the interface is achieved. After delamination, discontinuities are considered along the interface and, consequently, the interface starts interfering with the linear elastic fracture mechanics crack tip stress distribution. In the present model, on the other hand, interfaces are introduced as an additional compliant material from the very beginning of the simulation. In [12] it has been shown that the behaviour of such a system depends on the ratio between the interface process zone size and the bulk energy dissipation zone size (l CZM /l). Depending on this ratio, the apparent strength of the overall material ranges between the strengths ruled by the following limit models: a model with prefect bonded interfaces and finite l; a model with elastic bulk material and finite l CZM . In our simulations, the ratio l CZM /lis kept constant. In fact, the bulk properties do not change (then lis constant), and l CZM is also constant since kis fixed. This latter aspect is a consequence of dimensional analysis considerations in [12] from which we deduced that ∝=lEσEk//(2) c i CZM max 2 G . The consequence of having a constant ratio l CZM /lis that delamination is triggered earlier in the presence of a brittle interface. Due to this, we have also longer delamination paths for a brittle interface. These conclusions lead to the second important difference with respect to the Cook–Gordon model, where the load is supposed to be applied as a remote tensile stress at infinity. The assumption of infinite plane plays an important role, since there is no constrain on the size of the delamination path. As a result, a brittle interface could lead to a longer delamination with a consequent delayed crack propagation, increasing the elastic energy release rate. The specimens that we simulated are not long enough to be considered as infinite. Due to that, when the interface is brittle, delamination reaches the boundaries of the specimen causing the failure of the whole interface. This is also what accelerates the final failure of the specimen. As a result, for a brittle interface, the apparent strength of the specimen results lower than for a more ductile one. Finally, there are other assumptions made in the Cook–Gordon model that are not in line with the hypotheses of the present approach. For instance, Cook and Gordon neglected the effect of the vertical stress σ y in their interface fracture criterion. They made this simplification Fig. 6. Crack evolution in the simulation with brittle interface (σ c,0 = 1 MPa). V. Carollo et al. Journal of the European Ceramic Society 38 (2018) 2994–3003 2999 because their system is an homogeneous material with an interface, while here we have a laminate with different elastic properties for the laminae. This means that σ y is not constant across the laminae and there are jumps in correspondence of the interfaces. We cannot say quantitatively how much is the effect of the variation of σ y on the interface delamination and on the overall strength. Nevertheless, we have good motivations to believe that σ y cannot be neglected and could be another important source for possible discrepancies with respect to the Cook–Gordon model. 3.2. Crack propagation in Si 3 N 4 /BN micro-laminate In this section, we reproduce the experimental results in [4] concerning the 4-point bending test of a silicon nitride/boron nitride (Si 3 N 4 /BN) micro-laminate. Both constituent materials are brittle ceramics. The laminate is structured with layers of Si 3 N 4 of thickness between 40 μm and 60 μm, alternated by layers of BN of variable thickness between 2 μm and 10 μm(Fig. 9(a)). The BN layers act as quasi-brittle interface between the Si 3 N 4 layers. The 4-point bending test geometry is shown in Fig. 9(b). The dimensions of the specimen are: total span L= 5 mm, thickness T= 3 mm, width W= 4 mm. The outer and inner span of the 4-point bending test are, respectively, S 1 = 4 mm and S 2 = 2 mm. This geometry is discretized with phase field finite elements for the Si 3 N 4 layers, representing the bulk material, and cohesive interface finite elements compatible with phase field for the BN layers. Due to the variable thickness of the Si 3 N 4 layers, in our simulation the thickness associated to the layers of this material is set equal to 40, 50 or 60 μm. The assignment of the Si 3 N 4 layer thickness is randomly chosen according to a uniform distribution. The interface thickness, on the other hand, is set Fig. 7. Force–displacement curve of the three simulations. Fig. 8. Force–displacement curve compared with the crack and delamination length for the case of σ max = 100 MPa (left) and σ max = 1 MPa (right) normalized by the layer thickness l y . The letter refer to the corresponding image in Fig. 5 for the case of σ max = 100 MPa and to Fig. 6 for the case of σ max = 1 MPa. Fig. 9. (a) Si 3 N 4 /BN structure of the microlaminate; (b) 4-point bending experimental test geometry. Image (a) reprints from [4]. V. Carollo et al. Journal of the European Ceramic Society 38 (2018) 2994–3003 3000 Fig. 10. Numerically predicted vs. experimental force–displacement curves. Fig. 11. Fracture propagation and delamination resulting from the numerical simulation. For each subfigure there is the damage contour plot on top and the x-displacement contour plot on bottom. V. Carollo et al. Journal of the European Ceramic Society 38 (2018) 2994–3003 3001 constant and equal to 5 μm. The material and fracture parameters of the bulk are: Young modulus E = 310 GPa, Poisson ratio ν= 0.27, fracture energy s G =9μN/μm and phase field internal length scale l=2μm. The interface parameters are: σ c,0 =τ c,0 = 32 MPa and k 0 = 70 MPa/ μm. The experimental force–displacement curve (Fig. 10) shows an initial linear behaviour until the peak load of 475 N is reached. After this point, the load-carrying capacity of the specimen drops down to around 40% of the peak load. Then, the load continue increasing until a second drop is observed when the value of 240 N is reached. After this second drop, the load-carrying capacity is reduced to 10% of that at peak load. The specimen maintains this level until final failure. Fig. 11 shows the evolution of the predicted crack pattern. Due to symmetry of the geometry and of the boundary conditions, only half of the domain has been simulated, using around 380.000 finite element nodes and a finer mesh in the region where the crack is expected to propagate. The simulation shows an initial linear behaviour until the peak load of 475N is reached (Fig. 11(a)). Then, an interface situated in the middle of the specimen thickness fails due to delamination (Fig. 11(b)). After this first delamination event, the specimen continues gaining load-carrying capacity until the Si 3 N 4 layers start failing. According to the numerical predictions, the layers that first fail are that at the intrados and the one immediately next to the delaminated interface (Fig. 11(c)). The crack is predicted to continue its propagation towards the extrados of the specimen until final failure (Fig. 11(d)). The result of the simulation in terms of force–displacement curve shows a good agreement with the experimental one (Fig. 10). The initial part of the curve, the peak load and the drop in load-carrying capacity are very well predicted, capturing the critical load at which damage events occur. The final part of the simulation shows a more progressive failure evolution until collapse. This is not in perfect agreement with experiments, since we suppose that the second drop in the force noticed in experiments is caused by a second severe delamination event which we were not able to reproduce numerically. One possible source for this mismatch can be attributed to the interface parameters, since no experimental characterization was available. Another important result of the simulation is the very satisfactory prediction of the crack pattern features conforming to the experimental evidences. The experimental image in Fig. 12(a) shows the strong crack deflection and branching developed due to delamination. The same behaviour has been reproduced by the numerical simulation (Fig. 12(b)) where both phenomena are present in the simulated crack pattern. 4. Conclusions In this work, a novel variational framework combining the phase field method of fracture and the cohesive zone model approach has been applied for modelling complex crack paths in ceramic laminates within the context of the finite element method. In particular, the predictive capabilities of the proposed modelling tools have been exploited in relation to ceramic laminates. The effect of tough or quasi-brittle interfaces has been quantified and compared to the case of fully bonded layers. Based on the numerical predictions, we have numerically quantified that introducing interfaces is a way to increase the tortuosity of the crack path increasing the energy dissipated before failure. It has been also identified that the peak cohesive tractions govern the development of crack penetration, crack branching and crack deflection phenomena. Moreover, the applicability of the current computational method has been examined through its application to real laminates structures composed by Si 3 N 4 /BN layers. The numerical predictions showed that the proposed framework enabled reproducing the key features of the crack pattern and the force–displacement curve observed in the experiments, proving the capabilities of the proposed approach to be used as a tool for the design and characterization of ceramic laminates. In light of the previous arguments, this work is expected to provide a suitable modelling framework for more application-based problems concerning crack pattern predictions in ceramic composites, among other alternative multilayer or heterogeneous materials. Moreover, as a research perspective, a comparison of the proposed approach with the method provided in the theory of structured deformations [27–29] could be of interest for further investigations, to provide a re-interpretation of the phase field approach to fracture in terms of structured deformations. Acknowledgments The authors acknowledge funding received from the European Research Council under the European Union's H2020 Programme ERC Grant Agreement No. 737447 (ERC Proof of Concept 2016 “Photovoltaic with superior crack resistance”–PHYSIC). JR acknowledges the support of the project funded by the Spanish Ministry of Economy and Competitiveness/FEDERMAT2015-71036-P and the Andalusian Government Project of Excellence No. P12-TEP-1050. The authors would like to thank the anonymous reviewer for his useful suggestions. References [1] M. Sebastiani, K.-E. Johanns, E.-G. Herbert, F. Carassiti, G.-M. Pharr, A novel pillar indentation splitting test for measuring fracture toughness of thin ceramic coatings, Philos. Mag. 95 (16–18) (2015) 1928–1944. [2] N. Li, H. Wang, A. Misra, J. Wang, In situ nanoindentation study of plastic codeformation in Al-TiN nanocomposites, Sci. Rep. 4 (2014). Fig. 12. (a) Experimental crack pattern; (b) numerical predictions of the crack pattern (magnified). Image (a) adapted from [4]. V. Carollo et al. Journal of the European Ceramic Society 38 (2018) 2994–3003 3002