Fracture of solar-grade anisotropic polycrystalline Silicon: A combined phase field–cohesive zone model approach

Available online at www.sciencedirect.com ScienceDirect Comput. Methods Appl. Mech. Engrg. 330 (2018) 123–148 www.elsevier.com/locate/cma Fracture of solar-grade anisotropic polycrystalline Silicon: A combined phase field–cohesive zone model approach M. Paggia,∗, M. Corradob, J. Reinosoc aIMT School for Advanced Studies Lucca, Piazza San Francesco 19, 55100 Lucca, Italy bDepartment of Structural, Geotechnical and Building Engineering, Politecnico di Torino, Corso Duca degli Abruzzi 24, 10129, Torino, Italy cElasticity and Strength of Materials Group, School of Engineering, Universidad de Sevilla, Camino de los Descubrimientos s/n, 41092 Seville, Spain Received 26 May 2017; received in revised form 19 October 2017; accepted 20 October 2017 Available online 20 November 2017 Graphical Abstract Highlights •Combined phase-field and cohesive zone model for fracture. •Transgranular and intergranular fracture in solar-grade Silicon. •Anisotropic phase field model for fracture. •Identification of polycrystalline Silicon mechanical parameters. ∗Corresponding author. E-mail address: [email protected] (M. Paggi). https://doi.org/10.1016/j.cma.2017.10.021 0045-7825/ c 2017 The Author(s). Published by Elsevier B.V. This is an open access article under the CC BY-NC-ND license (http: //creativecommons.org/licenses/by-nc-nd/4.0/). 124 M. Paggi et al. / Comput. Methods Appl. Mech. Engrg. 330 (2018) 123–148 Abstract This work presents a novel computational framework to simulate fracture events in brittle anisotropic polycrystalline materials at the microscopical level, with application to solar-grade polycrystalline Silicon. Quasi-static failure is modeled by combining the phase field approach of brittle fracture (for transgranular fracture) with the cohesive zone model for the grain boundaries (for intergranular fracture) through the generalization of the recent FE-based technique published in [M. Paggi, J. Reinoso, Comput. Methods Appl. Mech. Engrg., 31 (2017) 145–172] to deal with anisotropic polycrystalline microstructures. The proposed model, which accounts for any anisotropic constitutive tensor for the grains depending on their preferential orientation, as well as an orientation-dependent fracture toughness, allows to simulate intergranular and transgranular crack growths in an efficient manner, with or without initial defects. One of the advantages of the current variational method is the fact that complex crack patterns in such materials are triggered without any user-intervention, being possible to account for the competition between both dissipative phenomena. In addition, further aspects with regard to the model parameters identification are discussed in reference to solar cells images obtained from transmitted light source. A series of representative numerical simulations is carried out to highlight the interplay between the different types of fracture occurring in solar-grade polycrystalline Silicon, and to assess the role of anisotropy on the crack path and on the apparent tensile strength of the material. c 2017 The Author(s). Published by Elsevier B.V. This is an open access article under the CC BY-NC-ND license (http: //creativecommons.org/licenses/by-nc-nd/4.0/). Keywords: Cohesive zone model; Phase field modeling of fracture; Anisotropic elasticity; Solar-grade polycrystalline silicon; Finite element method 1. Introduction Accurate prediction of failure in solids is a matter of a great relevance in different industrial sectors ranging from aerospace and aeronautics to renewable energy. At the macroscopic scale, fracture is generally originated by several microscopic defects, i.e. voids and cracks, which grow and coalesce leading to the deterioration of the overall material properties and the load carrying capacity of the structural component. In particular, at the microscopic scale, the interplay between intergranular (cracks propagating along grain boundaries) and transgranular (cracks propagating through the grains) fracture in polycrystalline materials is of major concern in material science and engineering applications. In the related literature, simplified models specific for intergranular fracture have been proposed in [1–5] based on the cohesive zone model (CZM) inserted along the internal grain boundaries. Cohesive crack formulations can be understood as phenomenological models in which fracture events are triggered by evaluating a particular traction– displacement law. This numerical technique inherently incorporates a characteristic length scale and has been widely used to simulate damage in polycrystalline materials [6], though crack paths are constrained along element edges, see a wide discussion in [7–9], among others. However, transgranular fracture is also a relevant failure mode in some materials and, in many cases, it is developed in competition with intergranular failure. Therefore, the development of numerical techniques that allow a reliable modeling of grains cracking with arbitrary failure paths are indeed necessary to provide a comprehensive simulation tool. In practical situations, for instance in micro-electro-mechanical systems (MEMS), both types of crack propagation have been reported in [10,11]. In those investigations, the authors comprehensively document the importance of considering the in plane grain anisotropy of Silicon which notably affects the fracture behavior. More specifically, regarding photovoltaics systems, solar-grade polycrystalline Silicon with grains randomly oriented in 3D are frequently used to manufacture thin solar cells. In such a case, both intergranular and transgranular fracture have been experimentally reported [12,13], see also some crack patterns in Fig. 1 experimentally observed by the present authors in solar cells embedded in photovoltaic modules, using the electroluminescence technique. A first attempt to model both types of fracture has been pursued in [14] by an intrinsic cohesive zone model approach. To avoid mesh dependency and un-physical material compliances, the extrinsic approach could be put forward, as in [15]. However, for very brittle materials like Silicon, cohesive zone model approaches suffer from the complexity in resolving the process zone which is very small as compared to the grain size or to the specimen size. In contrast to the previous numerical methods relying on cohesive zone models only, several studies have also been proposed to account for transgranular crack propagation in polycrystals by means of the extended finite element M. Paggi et al. / Comput. Methods Appl. Mech. Engrg. 330 (2018) 123–148 125 Fig. 1. Examples of crack patterns in polycrystalline Silicon solar cells embedded in photovoltaic modules. Dark areas in the images correspond to electrically inactive zones not contributing to energy production. method (X-FEM) [16], whereby the crack path does not explicitly depends upon the underlying spatial discretization of the physical domain into finite elements. In particular, Sukumar and co-authors [17] advocated a failure criterion based on the ratio between the fracture toughness of the grain boundaries Ggb cand the toughness of the grains Gi c, achieving a good qualitative agreement with respect to the experimental observations therein analyzed. Alternative computational techniques that feature explicit crack representations such as the enhanced finite element method (EFEM) [18,19], relying on the enhanced assumed variational methods [20], and the generalized finite element method (GFEM) [21] can also be used to model fracture events in such scenarios, though different operative difficulties such as the identification of crack initiation location and path might arise. Within the context of modeling sharp crack discontinuities, the recent phase field approaches for fracture [22,23] can be seen as a very promising numerical technique, especially suitable for quasi-brittle materials like Silicon (monocrystalline and polycrystalline solar-grade). Particularly, the explicit crack methods as those mentioned above suffer in situations with very complex fracture topologies including branching, kinking and coalescence. With the aim of overcoming such difficulties, grounded on the Griffith’s fracture theory [24], phase field models encompass a diffusive crack approach based on the definition of the crack phase field variable and its corresponding evolution equation. In the spirit of the Γ-convergence regularizations [25,26], discontinuities due to cracks are regularized using a characteristic length scale l. One of the most appealing aspects of phase field models of fracture is their inherent versatility, sharing several common aspects with gradient-enhanced damage formulations [27–30]. This feature has motivated the extension of the thermodynamic consistent formulation proposed by Miehe and co-authors [22,31] to dynamic [32,23] and ductile [33] fracture, anisotropic fracture in biomechanics [34], multi-physic systems [35,36], and shell structures [37–40], among many others. In this concern, recent investigations have applied the phase field approach to simulate failure in polycrystalline materials [41–43]. Moreover, its rational coupling with a novel cohesive formulation proposed in [44] allows fundamental Linear Elastic Fracture Mechanics (LEFM) predictions to be retrieved for the limit case of brittle cohesive interfaces in a consistent and robust manner. This formulation conceptually postulated the split of the dissipated energy into the bulk and interface contributions for the interplay between both failure mechanisms, i.e. crack growth in the bulk or debonding along the existing interfaces. Departing from the computational framework developed in [44], a novel formulation based on the phase field modeling of fracture for the grains combined with the cohesive zone model for the grain boundaries is proposed in this article. The principal motivation for the combination of both the cohesive zone and phase field models for triggering damage events in polycrystalline materials is to capture the fracture phenomena related to intergranular and transgranular crack growth and their potential interplay. Polycrystalline Silicon grains for solar-grade photovoltaics display a very brittle response, so that the use of the phase field approach to brittle fracture, herein extended to anisotropic materials with elastic and toughness properties dependent on the grain orientation, is a proper choice. Regarding the grain boundaries, on the other hand, a cohesive zone model approach is pursued to simulate their separation and allow for a richer description of their fracture behavior. In fact, although a tension cut-off cohesive zone model will be used in the present study, interface finite elements used to discretize the grain boundaries can accommodate more accurate traction–separation relations stemming from molecular dynamics or ab initio simulations, see e.g. [45,46], towards the development of a multi-scale computational model of technological impact for photovoltaic applications. Moreover, the proposed interface finite element formulation can be extended to take into account the additional thermal resistance of the grain boundaries for multi-field problems, which in turn does depend on the crack opening at the sub-micro scale, following the approach put forward in [47]. Within such a framework, 126 M. Paggi et al. / Comput. Methods Appl. Mech. Engrg. 330 (2018) 123–148 further developments might also regard modeling of the optoelectronic properties of grain boundaries [48], which is an important concern for the solar energy conversion efficiency. Based on the previous considerations, the main goal of the current approach concerns the simulation of fracture events in anisotropic polycrystalline materials with high technological impact, as is the case of solar-grade polycrystalline Silicon. First, the proposed formulation extends the concepts introduced in [44] for anisotropic materials via the definition of the fracture toughness of the grains as a function of the grain orientation. This modification will be used in the unified variational principle for the competition between intergranular and transgranular crack propagation. Second, due to the practical importance of polycrystalline Silicon in solar systems, special attention is devoted for identification of the model parameters for such applications through the use of solar cells images from transmitted light source. The manuscript is organized as follows. In Section 2, the principal aspects of the proposed modeling framework to simulate intergranular and transgranular failure in polycrystalline microstructures are presented, where special attention is paid to the modeling assumptions of anisotropic behavior of solar-grade Silicon. Section 3outlines the determination of the actual material properties of solar cells obtained from Silicon wafers in solar energy systems. Algorithmic details regarding the insertion of interface models within the grain boundaries are addressed in Section 4. The performance of the proposed computational method is assessed in Section 5by means of several representative examples. Finally, the main conclusions of the present study are drawn in Section 6. Additional electronic material collects MATLAB routines for the generation of finite element meshes for polycrystalline materials and for the computation of the anisotropic constitutive tensor of the solar-grade Silicon grains depending on their orientation, adopting the conventional Miller indices representation used in materials science. 2. Phase field approach for transgranular crack growth combined with the cohesive zone model for intergranular fracture This section revisits the main aspects of the consistent modeling framework for the combination of the phase field approach of fracture for transgranular crack growth and the cohesive zone model used to simulate intergranular fracture. Section 2.1 outlines the variational framework and the interface model, whereas the corresponding weak form and finite element formulation is presented in Section 2.2. Subsequently, Section 2.3 addresses the main aspects concerning the material modeling for anisotropic Silicon grains. As stated above, from the computational standpoint, the main novelty of the current study concerns the incorporation of the anisotropic character of polycrystalline materials –particularized to Silicon photovoltaics– in terms of elastic and fracture properties within the original modeling framework proposed in [44]. Therefore, this technique allows for the simulation of the competition between intergranular and transgranular crack propagation in an efficient and robust manner. 2.1. Variational framework and interface formulation Let Ω∈Rndim be an arbitrary cracked solid in the Euclidean space of dimension ndim , whose boundary is denoted by ∂Ω. Let Γ∈Rndim −1the crack surface within Ω, see Fig. 2.a. This system also includes the presence of an interface Γi. To clarify the notation, in the sequel, the position vectors of the grains are denoted by x, whereas xcstands for the position vectors referred to the interface Γi. The variational approach to brittle fracture governing crack nucleation, propagation and branching is set up through the definition of the following functional [22,23]: Π(u,Γ)=ΠΩ(u,Γ)+ΠΓ(Γ)=∫Ω\Γ ψe(ε) dΩ+∫Γ GcdΓ,(1) where ψe(ε) is the elastic energy density that depends upon the strain field ε, and Gcis the fracture energy. In Eq. (1), the term ΠΩ(u,Γ) 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 ΠΓ(Γ). Following [44,49], the fracture energy function is split into the corresponding counterparts associated with the dissipated energy in the anisotropic polycrystalline grains Ωdue to transgranular fracture (governed by the phase field approach of brittle fracture for the prospective discontinuities Γt) and along the grain boundaries (Γi) as a manifestation of intergranular fracture: ΠΓ=ΠΓt+ΠΓi=∫Γt G⟨abc⟩(u,d) dΓ+∫Γi Gi(u,d) dΓ.(2) M. Paggi et al. / Comput. Methods Appl. Mech. Engrg. 330 (2018) 123–148 127 Fig. 2. (a) The multi-field problem for cracked bodies based on the phase field approach of fracture with the presence of interfaces Γi. Left: Mechanical problem. Right: Crack phase field problem. (b) Traction–separation law of the interface model compatible with the phase field approach of fracture: representation for a generic fracture Mode that releases the energy Gi C:τ,τc,0and τcare the traction component at the interface, its critical value for undamaged and partially damaged surrounding bulk, whilst g,gc,0and gcare their respective relative displacements at the interface. In the grains, the fracture energy is dissipated according to the Griffith hypothesis [50]. In the present anisotropic model, the fracture toughness of each grain, G⟨abc⟩, depends on the grain orientation ⟨abc⟩, where a,b, and cdenote the standard Miller indices. The phase field variable in Eq. (2) has the physical meaning of an internal state damage variable (d∈[0,1], where d=0 represents an intact material, while d=1 identifies the fully damaged state), which is an implicit function of the displacement field and of a regularization parameter related to the smeared crack width, l. When the characteristic regularization parameter tends to zero (l→0), then the formulation outlined in Eq. (2) tends to Eq. (1) in the sense of the so-called Γ-convergence (see [50] for more details). For finite non-vanishing values of l, this parameter controls the value of the apparent material strength in uniaxial tensile tests, see e.g. [44]. At this point, it is worth mentioning that recent studies have investigated the relation between the parameter land the failure stress σc=√Gb cE/L, where Lis the characteristic size of the specimen and Gb cand Estand for the toughness and the Young’s modulus of the bulk, respectively. In particular, Tann´ e and co-authors [51] proposed that, once these material properties are known, i.e. σc, Gb cand E, then the parameter lcan be set as l=27 256 (Gb cE σ2 c ). Therefore, based on this consideration, the regularization parameter can be considered as a material internal length scale, in line with the discussion in [52,53]. For the grain boundaries, the corresponding fracture energy is released according to a cohesive zone model formulation. In particular, in the following we assume that the interface behavior is ruled by a linear traction– separation relation with tension cut-off, to simulate the brittle response of solar-grade polycrystalline Silicon. According to [44], the energy dissipation for intergranular fracture is characterized by the fracture energy function Gi, which is related to the displacement discontinuities at the interface, g, and also to the average phase field degradation variable of the adjacent continuum, d: Gi=G(g,d).(3) In the present study, a linear tension cut-off cohesive zone model (CZM) is adopted for Mode I and Mode II interface fracture, see [44] for more details omitted here for the sake of conciseness. The traction–separation relations are defined by the peak tractions σmax and τmax, and by the areas below them, which respectively correspond to the 128 M. Paggi et al. / Comput. Methods Appl. Mech. Engrg. 330 (2018) 123–148 grain boundary Mode I and Mode II fracture toughnesses, see Fig. 2.b for a schematic representation of the interface response. The particular forms of the traction–separation relation of the current cohesive model, which depends on the damage state of the surrounding bulk relying on the phase field variable, are given by: σ=⎧ ⎪ ⎨ ⎪ ⎩ kn gn gnc ,if 0 <gn gnc <1; 0,if gn gnc ≥1, τ=⎧ ⎪ ⎨ ⎪ ⎩ kt gt gtc ,if 0 <gt gtc <1; 0,if gt gtc ≥1. (4) It is worth mentioning that, differing from [43], the shape of the proposed cohesive model is consistently coupled with the phase field model in the bulk, but preserving the energy dissipated along the interface. The corresponding energy release rates for Modes I and II read: Gi I(d)=1 2kn,0g2 n g2 nc,0 [(1 −d)gnc,0+dgnc,1]2;Gi I I (d)=1 2kt,0g2 t g2 tc,0 [(1 −d)gtc,0+dgtc,1]2,(5) where kn,0and kt,0stand for the interface stiffness values for d=0, with their corresponding critical relative displacements gnc,0and gtc,0. The respective values for d=1 are denoted by kn,1,kt,1,gnc,1and gtc,1. Finally, for a generic value of the crack phase field variable in the bulk, d, the corresponding interface parameters read: kn,kt, gnc and gtc, whose definitions can be found in [44]. Interface failure is attained through the adoption of the standard fracture criterion based on the quadratic composition of the Mode I and Mode II energy release rates: (Gi I Gi IC )2 +(Gi I I Gi I IC )2 =1.(6) Based on the previous modeling assumptions, the functional in Eq. (1) can be recast as: Π(u,Γt,Γi)=ΠΩ+ΠΓt+ΠΓi=∫Ω\Γ ψe(ε) dΩ+∫Γt G⟨abc⟩(u,d) dΓ+∫Γi Gi(g,h,d)dΓ,(7) where the functional corresponding to the grains is: Πt(u,Γt)=ΠΩ(u,Γt)+ΠΓt(Γt)=∫Ω\Γ ψe(ε) dΩ+∫Γt G⟨abc⟩(u,d) dΓ.(8) Within the regularized framework of the phase field approach in continuous bodies [54], given a crack surface topology Γc⊂R2, the regularized crack functional adopts the form Γc(d):= ∫Ω γ(d,∇xd) dΩ,with γ(d,∇xd)=1 2ld2+l 2|∇xd|2.(9) In line with [22,31], the minimization problem associated with the diffusive crack topology reads d(x)=Arg{inf d∈WΓ(t) Γc(d)}(10) subject to the Dirichlet-type constraint: WΓ(t)= {d|d(x)=1 at x∈Γc}. Correspondingly, the potential energy of the system is decomposed into two terms: Πt(u,d)=∫Ω ψ(ε,d) dΩ+∫Ω G⟨abc⟩γ(d,∇xd) dΩ.(11) As a result, the total free energy density of the bulk, ˆ ψ, reads: ˆ ψ(ε,d)=ψ(ε,d)+G⟨abc⟩γ(d,∇xd).(12) The corresponding Euler equations associated with the phase field problem within the grains are given by d−l2∇2 xd=0 in Ωand ∇xd·n=0 in ∂Ω,(13) where ∇2 xdstands for the Laplacian of the phase field variable. M. Paggi et al. / Comput. Methods Appl. Mech. Engrg. 330 (2018) 123–148 129 Regarding the strain energy density of the bulk, ψ(ε,d), the following expression is herein considered, accounting for the anisotropic constitutive response of the grains: ψ(ε,d)=g(d)ψe(ε)=[(1 −d)2+K]ψe(ε).(14) In the previous equation, ψe(ε) is the elastic strain energy density of each undamaged anisotropic grain under plane stress assumptions, whereas g(d) stands for the degradation function, which satisfies g(0)ψe(ε)=ψe(ε) for undamaged states, and g(1)ψe(ε)=Kψe(ε) for fully deteriorated states, being Ka residual positive parameter that prevents numerical instabilities. Hence, the resulting form of the functional corresponding to the current modeling framework reads Π(u,d)=∫Ω g(d)ψe(ε) dΩ+∫Ω G⟨abc⟩γ(d,∇xd) dΩ+∫Γi Gi(g,d)dΓ.(15) Finally, note also that a suitable decomposition of the free-energy function into positive and negative counterparts as proposed in [22] is recalled to preclude crack propagation under compression. 2.2. Weak form of the variational problem and finite element formulation The construction of the variational formulation is carried out through the exploitation of the variation of the functional given in Eq. (15) with respect to the displacements and the crack phase field variable. For transgranular failure, we postulate kinematic and traction boundary conditions, which are defined along the disjoint parts ∂Ωu⊂∂Ωand ∂Ωt⊂∂Ω, respectively, with ∂Ωt∪∂Ωu=∂Ωand ∂Ωt∩∂Ωu= ∅, yielding: u=uon ∂Ωuand t=σ·non ∂Ωt,(16) where nis the outward normal unit vector to the body, and σis the Cauchy stress tensor. Recalling standard arguments, the following weak form for the variational problem for the grains is obtained: δΠt(u, δu,d, δd)=∫Ω σ:δεdΩ−∫Ω 2(1 −d)δdψe(ε) dΩ +∫Ω G⟨abc⟩l[1 l2dδd+∇xd·∇x(δd)]dΩ+δΠt,ext(u, δu),(17) where δucollects the test functions referred to the displacement field, which satisfy (Vu={δu|u=uon ∂Ωu,u∈ H1}); δdidentifies the phase field test functions (Vd={δd|δd=0 on Γb,d∈H0}), and finally, the external contribution to the variation of the bulk functional reads δΠt,ext(u, δu)=∫∂Ω t·δud∂Ω+∫Ω fv·δudΩ.(18) The irreversible character of the fracture process within the grains is guaranteed by means of the incorporation of a penalty term accounting for the history of the local damage variable [22,55]. The thermodynamic consistency according to the Clausius–Planck inequality of the present formulation has been comprehensively addressed in [22], and consequently specific details are omitted here for the sake of brevity. Similarly, with regard to the interface contribution to the complete functional of the system, its virtual variation with respect to the two primary fields takes the form δΠΓi(u, δu,d, δd)=∫Γi(∂Gi(u,d) ∂uδu+∂Gi(u,d) ∂dδd)dΓ,(19) where the displacement and crack phase field test functions are defined in accordance to the transgranular failure problem. The finite element discretization of the weak forms given in Eqs. (17)–(19) is performed complying with the isoparametric concept. Through the definition of the Lagrangian shape functions NI(ξ), which are defined in the parametric space ξ= {ξ1, ξ2}, the following discretization scheme for the displacements and the crack phase variable at the grains are introduced: x∼ =N˜ x;u∼ =Nd;δu∼ =Nδd;∆u∼ =N∆d,(20) 130 M. Paggi et al. / Comput. Methods Appl. Mech. Engrg. 330 (2018) 123–148 d∼ =Nd;δd∼ =Nδd;∆d∼ =N∆d,(21) where the interpolation functions are arranged in the operator N; the nodal reference positions are collected in the vector ˜ x, whereas the nodal displacements and crack phase field variable are denoted by the vectors dand d, respectively. Moreover, the nodal trial and incremental displacements are denoted by δdand ∆d, respectively, whereas the corresponding vectors for the crack phase field variable are identified by δdand ∆d. Note that the current modeling framework adopts the same interpolation scheme for the displacement and the phase field variables, as given in Eqs. (20)–(21). In particular, with respect to the phase field interpolation, the average value of the phase field at the two flanks of the grain boundary is computed, to be associated to the midline which is taken as the usual reference in CZM formulations (see [44]). This leads to a straightforward implementation within a monolithic fully implicit solution algorithm, as herein pursued. Nevertheless, it is worth mentioning that different interpolation schemes could be exploited and easily implemented within a staggered (sequential) approach, where the phase field Poisson-type partial differential equation is solved separately from the nonlinear mechanical equilibrium problem and an overall convergence loop is introduced. Based on the authors’ experience, the monolithic approach was found to be faster to converge than the staggered one, see the specific convergence study for fracture in solid shells reported in [40]. Within a monolithic strategy, coupling of different discretization schemes for the phase field and for the displacement field could be achieved by a Lagrange multiplier approach, such as the mortar method, which is a solution worth exploring in future developments. The interpolation of the strain field (ε) and the gradient of the phase field (∇xd), their variations (δε) and (∇xδd), and their increments (∆ε) and (∇x∆d) read ε∼ =Bdd;δε∼ =Bdδd;∆ε∼ =Bd∆d(22) ∇xd∼ =Bdd; ∇x(δd)∼ =Bdδd; ∇x(∆d)∼ =Bd∆d,(23) where Bddenotes the displacement–strain compatibility operator, and Bdidentifies the suitable matrix for the interpolation of the gradient of the phase field variable. After some straightforward operations, the consistent linearization of the discrete residual equations leads to the following coupled system: [Kt dd Kt dd Kt ddKt dd][∆d ∆d]=[ft d,ext 0]−[ft d,int ft d],(24) where Kt ab, with a,b=d,dare the tangent matrices to trigger transgranular failure within the grains, ft d,ext is the external force vector, whereas the internal residuals corresponding to the displacements and the crack phase field variable are denoted by ft d,int and ft d, respectively. In an analogous manner, for the intergranular failure along the interface, the local gap vector gloc at any point inside Γel ican be interpolated as follows gloc ∼ =ˆ Rˆ Bbd(25) where ˆ Rdenotes the standard rotation matrix of the interface traction–separation relation, whilst ˆ Bd=NL is the interface compatibility operator. For the crack phase field variables at the interface, the interpolation of the average value defined on the middle interface surface renders d∼ =NdMd¯ d=ˆ Bd¯ d,(26) where Mdis a suitable average operator and ˆ Bd=NdMdstands for the compatibility operator corresponding to the phase field. The particular expressions of the previous interpolation operators are detailed in [56,57] and are omitted here for brevity reasons. Finally, the discrete system of equations involving the displacement and the phase fields for intergranular failure renders [Ki dd Ki dd Ki ddKi dd][∆d ∆d]=[fi d fi d].(27) M. Paggi et al. / Comput. Methods Appl. Mech. Engrg. 330 (2018) 123–148 131 The particular expressions for tangent operators Ki ab, with a,b=d,d, and for the internal force vectors fi dand fi dare comprehensively derived in [44]. The numerical solution procedure for the systems Eqs. (24) and (27) is performed in a monolithic manner. In particular, two user-defined elements for the intergranular and transgranular failure are implemented into the FE code FEAP [58]. 2.3. Anisotropic polycrystalline Silicon: grains material model With regard to the material anisotropy, the stress–strain relation for a 3D Silicon polycrystal is usually specified in a frame with coordinate axes aligned along with its symmetry planes, i.e., using Miller’s indices, ⟨100⟩,⟨010⟩, and ⟨001⟩. In this material frame, the stress–strain relation can be symbolically written as σlocal i j =Clocal i jkl εlocal kl , where Clocal i jkl is the fourth order anisotropic constitutive tensor. For the particular case of Silicon, due to its cubic symmetry, we have C1111 =C2222 =C3333 =165.7 GPa, C1122 =C1133 =C2211 =C2233 =C3311 =C3322 =63.9 GPa, C1313 =C3113 =C1331 =C3131 =79.6 GPa, C2323 =C3223 =C2332 =C3232 =79.6 GPa, and C1212 =C2112 = C1221 =C2121 =79.6 GPa. To effectively account for the material properties in the general global setting, a generic orientation can be characterized by the three Euler angles defining the intersection between the surface plane [abc] of the grain defined by its normal vector ⟨abc⟩and the symmetry planes [100], [010], and [001]. We shall denote these angles as θ100,θ010, and θ001, for simplicity of the notation. They can be easily computed from the Miller’s indices a,b,c, as follows: θ100 =acos (a √a2+b2+c2), θ010 =acos (b √a2+b2+c2), θ001 =acos (c √a2+b2+c2). Thus, the stress–strain constitutive relation for each arbitrarily oriented grain in the global reference system reads σglobal i j =Cglobal i jkl εglobal kl , where Cglobal i jkl =Rpi Rq j Rrk Rsl Clocal pqrs (repeated indices denote summation from 1 to 3). The rotation matrix Rappearing in this transformation formula can be computed as the product R001R010R100, whose terms are defined as follows: R100 =⎡ ⎣ 1 0 0 0 cos(θ100)−sin(θ100) 0 sin(θ100) cos(θ100) ⎤ ⎦,(28a) R010 =⎡ ⎣ cos(θ010) 0 sin(θ010) 0 1 0 −sin(θ010) 0 cos(θ010) ⎤ ⎦,(28b) R001 =⎡ ⎣ cos(θ001)−sin(θ001) 0 sin(θ001) cos(θ001) 0 0 0 1 ⎤ ⎦.(28c) The derivation of a closed-form expression for Cglobal is very lengthy and impractical to display, so that it is more convenient to proceed numerically based on the above transformation. Moreover, since the computational model of the solar cell addressed in this work complies with plane stress modeling assumptions, the determination of the 2D constitutive stress–strain tensor C2D of the anisotropic grain is therefore required. This can be derived by imposing that the stress components σzz,σxz, and σyz are equal to zero. From the computational point of view, it is convenient to compute the tensor D=C−1 global relating the strain components to the stresses as a first operation, and then to take the sub-matrix D3×3whose components D3×3(i,j) are those contained in D(i,j) with i,j=1,2,6. The desired expression for C2D is finally obtained by inverting D3×3, i.e., C2D =D−1 3×3. The present method for the computation of the constitutive tensors has been implemented in the MATLAB routine Cmat.m, which is provided open source as a complementary electronic material of this article. 3. Anisotropic elastic parameters and fracture toughness depending on grain orientation The elastic properties of Silicon crystals are strongly dependent on their orientation, as carefully pointed out in [59], and their variations can be quite significant. Hence, stochastic models accounting for the random orientation of the 138 M. Paggi et al. / Comput. Methods Appl. Mech. Engrg. 330 (2018) 123–148 (a) v=16.5µm. (b) v=18.0µm. (c) v=18.5µm. (d) v=20.0µm. Fig. 7. Anisotropic material model. Contour lines of the inelastic vertical displacement field identifying the crack pattern, for Gi/G⟨100⟩=0.75. is composed of an initial transgranular path followed by a sequence of intergranular paths which are localizing the deformation along the grain boundaries and are inducing a stress relief in the continuum. A further increase of the toughness of the grain boundaries, Gi/G⟨100⟩=1.00, leads to a prevailing transgranular crack path over the whole specimen, completed by a short crack path due to intergranular crack growth for v≥20.0µm, see Fig. 8. Finally, considering grain boundaries tougher than the grains, Gi/G⟨100⟩=2.00, only transgranular crack growth is predicted to occur, see Fig. 9. Near the right hand side of the specimen, the crack path approaches the boundary between two adjacent grains, while remaining in the finite elements of the grains where the phase field variable attains unity. The average vertical stress, σ, computed as the sum of the nodal reaction forces in the vertical direction acting on the upper side, divided by the length of the same side and the unitary out-of-plane thickness, is plotted vs. the imposed vertical displacement vand the result is shown in Fig. 10 for the different values of the toughness ratio Gi/G⟨100⟩ previously considered. The case with Gi/G⟨100⟩=0.25, leading to a crack pattern at failure of pure intergranular type, is characterized by the minimum value of the peak stress. In all the other cases, the value of the peak stress is dictated by the transgranular crack growth into the first grain where the defect is located, which is a common feature for all the cases with Gi/G⟨100⟩≥0.75. The post-peak branch of the stress–displacement curve, on the other hand, is influenced by the last part of the crack pattern, which can be intergranular or transgranular depending on the grain boundary toughness. In the case of intergranular crack growth, a more brittle post-peak branch is noticed, see the curve in Fig. 10 corresponding to Gi/G⟨100⟩=0.75. M. Paggi et al. / Comput. Methods Appl. Mech. Engrg. 330 (2018) 123–148 139 (a) v=17.5µm. (b) v=19.0µm. (c) v=20.0µm. (d) v=21.0µm. Fig. 8. Anisotropic material model. Contour lines of the inelastic vertical displacement field identifying the crack pattern, for Gi/G⟨100⟩=1.00. 5.2. Isotropic model with an initial defect In this test, the same geometry and boundary conditions as for the anisotropic case are considered, while all the grains are assumed to be ⟨100⟩equi-oriented, as a possible model for isotropic monocrystalline Silicon. This example encompasses a simple isotropic model, so that anisotropic effects in terms of stiffness and toughness properties are herewith precluded. In line with the previous application, considering a grain boundary toughness Gi/G⟨100⟩=2.00, transgranular cracking is prevailing over intergranular one, as for the anisotropic case, see the crack pattern shown in Fig. 11. However, as compared to the anisotropic case, the crack path at failure is simply straight through the specimen, rather than curved as in Fig. 9. Therefore, these results pinpoint that elastic anisotropy and the dependence of the fracture toughness on the grain orientation play a decisive role on the shape of the crack path at failure. Although the crack pattern of the isotropic material model is different from the anisotropic one, the curves relating the average vertical stress, σ, to the imposed vertical displacement, v, are similar, see Fig. 12. 5.3. Anisotropic model without defects Finally, in this subsection, the same anisotropic model as in the test problem whose mechanical response was illustrated in Fig. 9 is examined, removing the initial defect inside the grain on the left hand side of the microstructure. Considering the high toughness of the interface as compared to that of the grains, Gi/G⟨100⟩=2.00, transgranular 140 M. Paggi et al. / Comput. Methods Appl. Mech. Engrg. 330 (2018) 123–148 (a) v=16.8µm. (b) v=18.0µm. (c) v=19.2µm. (d) v=20.0µm. Fig. 9. Anisotropic material model. Contour lines of the inelastic vertical displacement field identifying the crack pattern, for Gi/G⟨100⟩=2.00. Fig. 10. Average vertical stress σvs. imposed vertical displacement vfor the anisotropic model and for different values of the toughness ratio, Gi/G⟨100⟩. M. Paggi et al. / Comput. Methods Appl. Mech. Engrg. 330 (2018) 123–148 141 (a) v=16.8µm. (b) v=18.0µm. (c) v=20.0µm. (d) v=22.0µm. Fig. 11. Isotropic material model. Contour lines of the inelastic vertical displacement field identifying the crack pattern, for Gi/G⟨100⟩=2.00. Fig. 12. Average vertical stress σvs. imposed vertical displacement vfor the isotropic and the anisotropic material models, for Gi/G⟨100⟩=2.00. 142 M. Paggi et al. / Comput. Methods Appl. Mech. Engrg. 330 (2018) 123–148 (a) v=25.2µm. (b) v=26.0µm. (c) v=26.8µm. (d) v=27.6µm. Fig. 13. Anisotropic material model without defects. Contour lines of the inelastic vertical displacement field identifying the crack pattern, for Gi/G⟨100⟩=2.00. fracture is expected to prevail over the intergranular one. The evolution of the crack pattern shown in Fig. 13 confirms the trend observed in the presence of a defect, i.e., transgranular fracture only. In the absence of defects, however, the grain boundaries near the upper side of the microstructure contribute to the localization of deformation, with a phase field crack propagating very close to the material interfaces, thus leading to the occurrence of a sub-interfacial crack growth. In terms of the overall mechanical response, see Fig. 14, the presence of an initial defect significantly reduces the value of the peak stress as compared to the un-notched case, as expected also from classical theoretical considerations. In the latter situation, the value of the peak stress significantly increases towards the nominal strength of monocrystalline Silicon, i.e., 190 MPa, but it is still lower due to strain localization effects induced by the presence of grain boundaries and material anisotropy. 6. Conclusions In this work, a new computational framework to simulate transgranular and intergranular crack growth in anisotropic polycrystalline microstructures has been proposed. Inspired by [44], the proposed methodology enables the incorporation of the anisotropic character of solar-grade Silicon in a robust and straightforward manner. The current approach overcomes the difficulties associated with the competition between both dissipative phenomena. To do so, the interplay between the phase field approach for transgranular cracking and the cohesive zone model M. Paggi et al. / Comput. Methods Appl. Mech. Engrg. 330 (2018) 123–148 143 Fig. 14. Average vertical stress σvs. imposed vertical displacement vfor the anisotropic model with or without an initial defect, for Gi/G⟨100⟩=2.00. for intergranular failure is carried out by stating an interface response whose stiffness depends upon the damage accumulation in the surrounding bulk. This approximation leads to a fully coupled system of equations that is solved by means of a monolithic Newton–Raphson solution scheme for the two primary fields, i.e. the displacements and the crack phase field variable. In addition, due to its practical importance, special attention has been devoted to the identification of the anisotropic elastic parameters of the grains and the fracture toughness in solar-grade Silicon cells, which depend on grain orientation. The proposed computational method and its predictive capabilities have been assessed by means of a series of representative tests concerning anisotropic polycrystalline Silicon microstructures. These applications show that the proposed approach is able to handle complex crack patterns and the combination of both fracture modes mentioned above. In order to draw more general conclusions, the proposed model can be applied in the future within a statistical framework, for instance Monte Carlo-based simulations, to analyze the effects of the variability of the material properties on the resulting crack pattern (as it was done, for instance, in [76] for polysilicon in MEMS). On the other hand, it is relevant to recall that recent experimental studies on full size solar-grade polycrystalline Silicon [77] have highlighted a very good reproducibility of the crack pattern in the case of twin specimens produced with exactly the same grain geometry. This evidence suggests that the scatter in the elastic and fracture properties from a specimen to another is almost negligible and that only the grain-size distribution plays a key role. Besides a statistical study, also the size of the microstructural model should be object of further studies, since the apparent strength and the crack pattern are indeed affected by the ratio between the process zone size of the intergranular cohesive cracks and the specimen size, as shown in previous studies on intergranular fracture of polycrystals [1,4]. Acknowledgments The research presented in this article has received funding from the European Research Council through 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 no. 306622. JR is also grateful to the support to the projects funded by the Spanish Ministry of Economy and Competitiveness (Projects MAT2015–71036-P and MAT2015–71309-P) and the Andalusian Government (Projects of Excellence No. TEP-7093 and P12-TEP-1050). 144 M. Paggi et al. / Comput. Methods Appl. Mech. Engrg. 330 (2018) 123–148 Appendix A. Anisotropic elastic coefficients of the grains This Appendix collects the 3D and 2D plane stress anisotropic elastic coefficients of the grains with the most common orientations observed in solar cells. These coefficients have been computed following the procedure presented in Section 3, by preliminary evaluating the three Euler angles for each orientation using the MATLAB script provided as additional electronic material. The unit of measure of all the coefficients is GPa. For Silicon grains whose surface is aligned along the [100] plane (⟨100⟩normal vector), the three Euler angles formed with the [100], [010], and [001] planes are θ100 =0◦,θ010 =90◦, and θ001 =90◦, respectively. As a result, the 3D anisotropic constitutive tensor reads: C100 = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 165.7 63.9 63.9000 63.9 165.7 63.9000 63.9 63.9 165.7000 0 0 0 79.600 0 0 0 0 79.6 0 0 0 0 0 0 79.6 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ (A.1) and its plane stress counterpart is: C100 2D =⎡ ⎣ 141.1 39.3 0.0 39.3 141.1 0.0 0.0 0.0 79.6 ⎤ ⎦.(A.2) For Silicon grains whose surface is aligned along the [110] plane (⟨110⟩normal vector), the three Euler angles are θ100 =45◦,θ010 =45◦, and θ001 =90◦. The 3D anisotropic constitutive tensor reads: C110 = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 194.4 49.5 49.5−14.300 49.5 201.5 42.4 7.200 49.5 42.4 201.6 7.200 −14.3 7.2 7.2 58.100 0 0 0 0 65.2−14.3 0 0 0 0 −14.3 65.2 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ (A.3) and the plane stress tensor is: C110 2D =⎡ ⎣ 177.7 40.7 0 40.7 192.1 0 0 0 62.1 ⎤ ⎦.(A.4) For Silicon grains whose surface is aligned along the [111] plane (⟨111⟩normal vector), the three Euler angles are θ100 =θ010 =θ001 =54.7◦. The 3D anisotropic constitutive tensor reads: C111 = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 194.1 45.5 53.9−4.5 8.9 9.6 45.5 186.2 61.8 1.4 5.1−13.4 53.9 61.8 177.8 3.1−14.1 3.7 −4.5 1.4 3.1 77.5 3.7 5.1 8.9 5.1−14.1 3.7 69.6−4.5 9.6−13.4 3.7 5.1−4.5 61.2 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ (A.5) and the plane stress tensor is: C111 2D =⎡ ⎣ 174.6 24.8 9.7 24.8 163.3−14.0 9.7−14.0 60.5 ⎤ ⎦.(A.6) M. Paggi et al. / Comput. Methods Appl. Mech. Engrg. 330 (2018) 123–148 145 Table A.1 Eigenvalues of the constitutive tensor and Young’s modulus, E, depending on the grain orientation. Grain orientation λ1(GPa) λ2(GPa) λ3(GPa) E(GPa) ⟨100⟩180.3 101.8 79.6 130.1 ⟨110⟩226.3 143.6 62.1 169.1 ⟨111⟩194.4 146.8 57.1 187.8 ⟨211⟩214.2 130.4 66.8 169.1 ⟨344⟩198.9 150.2 56.3 185.4 For Silicon grains whose surface is aligned along the [211] plane (⟨211⟩normal vector), the three Euler angles are θ100 =35.26◦,θ010 =65.90◦, and θ001 =65.90◦. The 3D anisotropic constitutive tensor reads: C211 = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 182.1 55.5 55.9−5.9 13.0 8.0 55.5 193.1 44.9−6.9−7.6−9.4 55.9 44.8 192.7 12.9−5.4 1.4 −5.9−6.9 12.9 60.6 1.4−7.6 13.0−7.6−5.4 1.4 71.6−5.9 8.0−9.4 1.4−7.6−5.9 71.2 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ (A.7) and the plane stress tensor is: C211 2D =⎡ ⎣ 161.2 42.1 7.5 42.1 180.5−11.5 7.5−11.5 69.8 ⎤ ⎦.(A.8) For Silicon grains whose surface is aligned along the [344] plane (⟨344⟩normal vector), the three Euler angles are θ100 =62.06◦,θ010 =51.34◦, and θ001 =51.34◦. The 3D anisotropic constitutive tensor reads: C344 = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 197.2 43.8 52.6−3.1 6.6 9.3 43.8 190.3 59.4 5.9 9.1−8.6 52.6 59.4 181.5−2.8−15.7−0.8 −3.1 5.9−2.8 75.1−0.8 9.1 6.6 9.1−15.7−0.8 68.3−3.1 9.3−8.6−0.8 9.1−3.1 59.5 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ (A.9) and the plane stress tensor is: C344 2D =⎡ ⎣ 180.0 24.4 10.3 24.4 167.2−8.5 10.3−8.5 58.2 ⎤ ⎦.(A.10) To appreciate the mismatch in the elastic constitutive properties, the eigenvalues of the above constitutive tensors are collected in Table A.1, together with the corresponding Young’s modulus. In particular, the Young’s moduli for the different orientations of the surface of the polycrystalline grains, defined by the normal vector ⟨abc⟩, have been computed from the coefficients of the compliance tensor D100, by applying the following equation [78]: 1 Eabc =D11 −2(D11 −D12 −1 2D66)(cos(θ100)2cos(θ010)2+cos(θ010)2cos(θ001)2 +cos(θ100)2cos(θ001)2)(A.11) where cos(θ100), cos(θ010) and cos(θ001) are the direction cosines as defined in Section 2.3. Appendix B. Supplementary data Supplementary electronic material related to this article can be found online at http://dx.doi.org/10.1016/j.cma. 2017.10.021. 146 M. Paggi et al. / Comput. Methods Appl. Mech. Engrg. 330 (2018) 123–148 References [1] M. Paggi, P. Wriggers, A nonlocal cohesive zone model for finite thickness interfaces–Part I: mathematical formulation and validation with molecular dynamics, Comput. Mater. Sci. 50 (5) (2011) 1625–1633. [2] M. Paggi, P. Wriggers, A nonlocal cohesive zone model for finite thickness interfaces–Part II: FE implementation and application to polycrystalline materials, Comput. Mater. Sci. 50 (5) (2011) 1634–1643. [3] M. Paggi, E. Lehmann, C. Weber, A. Carpinteri, P. Wriggers, M. Schaper, A numerical investigation of the interplay between cohesive cracking and plasticity in polycrystalline materials, Comput. Mater. Sci. 77 (2013) 81–92. [4] M. Paggi, P. Wriggers, Stiffness and strength of hierarchical polycrystalline materials with imperfect interfaces, J. Mech. Phys. Solids 60 (4) (2012) 557–572. [5] H.D. Espinosa, P.D. Zavattieri, A grain level model for the study of failure initiation and evolution in polycrystalline brittle materials. Part II: Numerical examples, Mech. Mater. 35 (3–6) (2003) 365–394. [6] M. Ortiz, A. Pandolfi, Finite deformation irreversible cohesive elements for three-dimensional crack-propagation analysis, Internat. J. Numer. Methods Engrg. 44 (1999) 1267–1282. [7] T. Luther, C. Könke, Polycrystal models for the analysis of intergranular crack growth in metallic materials, Eng. Fract. Mech. 76 (15) (2009) 2332–2343. [8] C.V. Verhoosel, M.A. Gutiérrez, Modelling inter- and transgranular fracture in piezoelectric polycrystals, Eng. Fract. Mech. 76 (6) (2009) 742–760. [9] M. Paggi, M. Corrado, M.A. Rodriguez, A multi-physics and multi-scale numerical approach to microcracking and power-loss in photovoltaic modules, Compos. Struct. 95 (2013) 630–638. [10] A. Corigliano, A. Ghisi, G. Langfelder, A. Longoni, F. Zaraga, A. Merassi, A microsystem for the fracture characterization of polysilicon at the micro-scale, Eur. J. Mech. A Solids 30 (2) (2011) 127–136. [11] R. Vayrette, M. Galceran, M. Coulombier, S. Godet, J.-P. Raskin, T. Pardoen, Size dependent fracture strength and cracking mechanisms in freestanding polycrystalline silicon films with nanoscale thickness, Eng. Fract. Mech. 168 (Part A) (2016) 190–203. [12] M. Sander, S. Dietrich, M. Pander, M. Ebert, J. Bagdahn, Systematic investigation of cracks in encapsulated solar cells after mechanical loading, Sol. Energy Mater. Sol. Cells 111 (2013) 82–89. [13] L. Zhao, A. Maynadier, D. Nelias, Stiffness and fracture analysis of photovoltaic grade silicon plates, Int. J. Solids Struct. 97–98 (2016) 355–369. [14] A. Infuso, M. Corrado, M. Paggi, Image analysis of polycrystalline solar cells and modelling of intergranular and transgranular cracking, J. Eur. Ceram. Soc. 34 (11) (2014) 2713–2722. [15] M. Corrado, J.-F. Molinari, Effects of residual stresses on the tensile fatigue behavior of concrete, Cem. Concr. Res. 89 (2016) 206–219. [16] N. Moës, J. Dolbow, T. Belytschko, A finite element method for crack growth without remeshing, Internat. J. Numer. Methods Engrg. 46 (1) (1999) 131–150. [17] N. Sukumar, D.J. Srolovitz, T.J. Baker, J.-H. Prévost, Brittle fracture in polycrystalline microstructures with the extended finite element method, Internat. J. Numer. Methods Engrg. 56 (14) (2003) 2015–2037. [18] J. Oliver, A.E. Huespe, S. Blanco, D.L. Linero, Stability and robustness issues in numerical modeling of material failure with the strong discontinuity approach, Comput. Methods Appl. Mech. Engrg. 195 (52) (2006) 7093–7114. [19] C. Linder, F. Armero, Finite elements with embedded strong discontinuities for the modeling of failure in solids, Internat. J. Numer. Methods Engrg. 72 (12) (2007) 1391–1433. [20] J.C. Simo, J. Oliver, F. Armero, An analysis of strong discontinuities induced by strain-softening in rate-independent inelastic solids, Comput. Mech. 12 (5) (1993) 277–296. [21] A. Simone, C.A. Duarte, E. Van der Giessen, A generalized finite element method for polycrystals with discontinuous grain boundaries, Internat. J. Numer. Methods Engrg. 67 (8) (2006) 1122–1145. [22] C. Miehe, M. Hofacker, F. Welschinger, A phase field model for rate-independent crack propagation: Robust algorithmic implementation based on operator splits, Comput. Methods Appl. Mech. Engrg. 199 (45) (2010) 2765–2778. [23] M.J. Borden, C.V. Verhoosel, M.A. Scott, T.J.R. Hughes, C.M. Landis, A phase-field description of dynamic brittle fracture, Comput. Methods Appl. Mech. Engrg. 217 (2012) 77–95. [24] A.A. Griffith, The phenomena of rupture and flow in solids, Phil. Trans. R. Soc. Lond. Ser. A, Contain. Pap. Math. Phys. Character 221 (1921) 163–198. [25] L. Ambrosio, V.M. Tortorelli, Approximation of functional depending on jumps by elliptic functional via t-convergence, Comm. Pure Appl. Math. 43 (8) (1990) 999–1036. [26] G.A. Francfort, J.-J. Marigo, Revisiting brittle fracture as an energy minimization problem, J. Mech. Phys. Solids 46 (8) (1998) 1319–1342. [27] M. Frémond, B. Nedjar, Damage, gradient of damage and principle of virtual power, Int. J. Solids Struct. 33 (8) (1996) 1083–1103. [28] R.H.J. Peerlings, M.G.D. Geers, R. De Borst, W.A.M. Brekelmans, A critical comparison of nonlocal and gradient-enhanced softening continua, Internat J. Solids Struct. 38 (44) (2001) 7723–7746. [29] C. Comi, U. Perego, Fracture energy based bi-dissipative damage model for concrete, Int. J. Solids Struct. 38 (36) (2001) 6427–6454. [30] B.J. Dimitrijevic, K. Hackl, A regularization framework for damage–plasticity models via gradient enhancement of the free energy, Int. J. Numer. Methods Biomed. Eng. 27 (8) (2011) 1199–1210. [31] C. Miehe, F. Welschinger, M. Hofacker, Thermodynamically consistent phase-field models of fracture: Variational principles and multi-field FE implementations, Internat. J. Numer. Methods Engrg. 83 (10) (2010) 1273–1311. [32] M. Hofacker, C. Miehe, A phase field model of dynamic fracture: Robust field updates for the analysis of complex crack patterns, Internat. J. Numer. Methods Engrg. 93 (3) (2013) 276–301. M. Paggi et al. / Comput. Methods Appl. Mech. Engrg. 330 (2018) 123–148 147 [33] M. Ambati, T. Gerasimov, L. De Lorenzis, Phase-field modeling of ductile fracture, Comput. Mech. 55 (5) (2015) 1017–1040. [34] O. Gültekin, Hgf. Dal, G.A. Holzapfel, A phase-field approach to model fracture of arterial walls: Theory and finite element analysis, Comput. Methods Appl. Mech. Engrg. 312 (2016) 542–566 Phase Field Approaches to Fracture. [35] C. Miehe, L.M. Schänzel, H. Ulmer, Phase field modeling of fracture in multi-physics problems. Part I. Balance of crack surface and failure criteria for brittle crack propagation in thermo-elastic solids, Comput. Methods Appl. Mech. Engrg. 294 (2015) 449–485. [36] X. Zhang, A. Krischok, C. Linder, A variational framework to model diffusion induced large plastic deformation and phase field fracture during initial two-phase lithiation of silicon electrodes, Comput. Methods Appl. Mech. Engrg. 312 (2016) 51–77. [37] F. Amiri, D. Millán, Y. Shen, T. Rabczuk, M. Arroyo, Phase-field modeling of fracture in linear thin shells, Theor. Appl. Fract. Mech. 69 (2014) 102–109. [38] P. Areias, T. Rabczuk, M.A. Msekh, Phase-field analysis of finite-strain plates and shells including element subdivision, Comput. Methods Appl. Mech. Engrg. 312 (2016) 322–350. [39] J. Kiendl, M. Ambati, L. De Lorenzis, H. Gomez, A. Reali, Phase-field description of brittle fracture in plates and shells, Comput. Methods Appl. Mech. Engrg. 312 (2016) 374–394. [40] J. Reinoso, M. Paggi, C. Linder, Phase field modeling of brittle fracture for enhanced assumed strain shells at large deformations: formulation and finite element implementation, Comput. Mech. 59 (2017) 981–1001. [41] J.D. Clayton, J. Knap, Phase field modeling of directional fracture in anisotropic polycrystals, Comput. Mater. Sci. 98 (2015) 158–169. [42] J.D. Clayton, J. Knap, Phase field modeling and simulation of coupled fracture and twinning in single crystals and polycrystals, Comput. Methods Appl. Mech. Engrg. 312 (2016) 447–467 Phase Field Approaches to Fracture. [43] T.-T. Nguyen, J. Réthoré, J. Yvonnet, M.-C. Baietto, Multi-phase-field modeling of anisotropic crack propagation for polycrystalline materials, Comput. Mech. (2017) 1–26. [44] M. Paggi, J. Reinoso, 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, Comput. Methods Appl. Mech. Engrg. 321 (2017) 145–172. [45] P.R. Budarapu, J. Reinoso, M. Paggi, Concurrently coupled solid shell-based adaptive multiscale method for fracture, Comput. Methods Appl. Mech. Engrg. 319 (2017) 338–365. [46] Y. Umeno, A. Kushima, T. Kitamura, P. Gumbsch, Ab initio study of the surface properties and ideal strength of (100) silicon thin films, Phys. Rev. B 72 (2005) 165431. [47] A. Sapora, M. Paggi, A coupled cohesive zone model for transient analysis of thermoelastic interface debonding, Comput. Mech. 53 (2014) 845–857. [48] J.-W. Lee, S.-H. Bae, N. De Marco, Y.-T. Hsieh, Z. Dai, Y. Yang, The role of grain boundaries in perovskite solar cells, Mater. Today Energy, (in press). [49] C.V. Verhoosel, R. de Borst, A phase-field model for cohesive fracture, Internat. J. Numer. Methods Engrg. 96 (1) (2013) 43–62. [50] H. Amor, J.-J. Marigo, C. Maurini, Regularized formulation of the variational brittle fracture with unilateral contact: numerical experiments, J. Mech. Phys. Solids 57 (8) (2009) 1209–1229. [51] E. Tanné, T. Li, B. Bourdin, J.-J. Marigo, C. Maurini, Crack nucleation in variational phase-field models of brittle fracture, J. Mech. Phys. Solids 110 (2017) 80–99. [52] C. Comi, U. Perego, Criteria for mesh refinement in nonlocal damage finite element analyses, Eur. J. Mech. A Solids 23 (4) (2004) 615–632. [53] C. Comi, S. Mariani, U. Perego, An extended FE strategy for transition from continuum damage to mode I cohesive crack propagation, Int. J. Numer. Anal. Methods Geomech. 31 (2007) 213–238. [54] B. Bourdin, G.A. Francfort, J.-J. Marigo, The variational approach to fracture, J. Elasticity 91 (1) (2008) 5–148. [55] M.-A. Msekh, J.-M. Sargado, M. Jamshidian, P.M. Areias, T. Rabczuk, Abaqus implementation of phase-field model for brittle fracture, Comput. Mater. Sci. 96 (2015) 472–484. [56] J. Reinoso, M. Paggi, A consistent interface element formulation for geometrical and material nonlinearities, Comput. Mech. 54 (6) (2014) 1569–1581. [57] M. Paggi, J. Reinoso, An anisotropic large displacement cohesive zone model for fibrillar and crazing interfaces, Int. J. Solids Struct. 69 (2015) 106–120. [58] O.C. Zienkiewicz, R.L. Taylor, The Finite Element Method: Solid Mechanics, Vol. 2, Butterworth-Heinemann, 2000. [59] M.A. Hopcroft, W.D. Nix, T.W. Kenny, What is the Young’s modulus of Silicon? IEEE J. Micromechanical Syst. 19 (2010) 229–238. [60] A. Masolin, P.O. Bouchard, R. Martini, M. Bernacki, Thermo-mechanical and fracture properties in single-crystal silicon, J. Mater. Sci. 48 (3) (2013) 979–988. [61] A. Ghisi, F. Fachin, S. Mariani, S. Zerbini, Multi-scale analysis of polysilicon MEMS sensors subject to accidental drops: Effect of packaging, Microelectron. Reliab. 49 (3) (2009) 340–349. [62] A. Ghisi, S. Kalicinski, S. Mariani, I. De Wolf, A. Corigliano, Polysilicon MEMS accelerometers exposed to shocks: numerical-experimental investigation, J. Micromech. Microeng. 19 (3) (2009) 035023. [63] K. Fujiwara, Crystal growth behaviors of Silicon during melt growth processes, Int. J. Photoenergy (2012) 169829. [64] T. Lehmann, M. Trempa, E. Meissner, M. Zschorsch, M. Reimann, J. Friedrich, Laue scanner: A new method for determination of grain orientations and grain boundary types of multicrystalline silicon on a full wafer scale, Acta Mater. 69 (2014) 1–8. [65] D. Lausch, M. Glaeser, C. Hagendorf, Determination of crystal grain orientations by optical microscopy at textured surfaces, J. Phys. D: Appl. Phys. 114 (2013) 194509. [66] B. Sopori, D. Guhabiswas, P. Rupnowski, S. Shet, S. Devayajanam, H. Moutinho, A new method for rapid measurement of orientations and sizes of grains in multicrystalline silicon wafers, in: 37th IEEE Photovoltaic Specialists Conference, June 2011, pp. 001680–001685.