Phase field modeling of anisotropic silicon crystalline cracking in 3D thin-walled photovoltaic laminates

Int J Fract https://doi.org/10.1007/s10704-024-00821-3 RESEARCH Phase field modeling of anisotropic silicon crystalline cracking in 3D thin-walled photovoltaic laminates Z. Liu ·P. Lenarda ·J. Reinoso ·M. Paggi Received: 26 May 2024 / Accepted: 15 November 2024 © The Author(s) 2025 Abstract Anovelcomputationalframeworkintegrating the phase field approach with the solid shell formulation at finite deformation is proposed to model the anisotropic fracture of silicon solar cells in the thin-walled photovoltaic laminates. To alleviate the locking effects, both the enhanced assumed strain and assumed natural strain methods are incorporated in the solid shell element formulation. Aiming at tackling the poor convergence performance of standard Newton schemes, the efficient and robust quasi-Newton scheme is adopted for the solution of phase field modeling with enhanced shell formulation in a monolithic manner. Due to fracture anisotropy of the brittle silicon solar cells, the second-order structural tensor that is defined by the normal of preferential crack plane is introduced into the crack energy density function in Z. Liu ·P. Lenarda ·M. Paggi IMT School for Advanced Studies, Piazza San Francesco 19, 55100 Lucca, Italy e-mail: [email protected] Z. Liu·J. Reinoso (B) Grupo de Ingeniería de Estructuras y Materiales. Departamento de Mecánica de Medios Continuos y Teoría de Estructuras, ETS Ingeniería, Universidad de Sevilla, Camino de los Descubrimientos S/N, 41092 Seville, Spain e-mail: [email protected] J. Reinoso Escuela Politécnica Superior, Universidad de Sevilla, C/ Virgen de África, 7, Sevilla 41011, Spain J. Reinoso ENGREEN - Laboratory of Engineering for Energy and Environmental Sustainability, Universidad de Sevilla, Seville, Spain the phase field modeling. On the other hand, to efficiently predict the crack growth of silicon solar cells, a global–local approach in the 3D setting proposed in the previous work is adopted here for the fracture modeling. In this approach, both mechanical deformation and phasefieldfractureareaccountedfor atthelocalmodel, while only mechanical deformation is addressed at the globallevel.Ateachtimestep,thesolutionoftheglobal model is used to drive the local model, which correspondsto the one-way coupling inline with experimental evidence that the silicon cell cracking has negligible influence on the stiffness of photovoltaic modules. The capability of the modeling framework is demonstrated through numerical simulation of silicon solar cell cracking in the photovoltaic modules when subjected to different loading cases. Keywords Global–local approach ·Phase field fracture ·Solid shell element ·Quasi–Newton scheme 1 Introduction Nowadays, crystalline silicon photovoltaic (PV) modules have been massively deployed all over the world since the last century, and according to Xu et al. (2021), the installed capacity of PV approximately increases from 40 GW to 715 GW in the recent 10 years. Among the major components of PV products, including tempered glass, ethylene-co-vinyl acetate (EVA), silicon cells, and backsheet, the silicon wafer accounts for 123 Z. Liu et al. Fig. 1 Different crack patterns in silicon solar cells (Papargyri et al. 2020) more than 40% of the manufacturing cost of crystalline modules (Papargyri et al. 2020). It was proposed in the 2010 International Technology Roadmap for Photovoltaicsthatthethicknessofsiliconcellshouldbesignificantly decreased so as to save cost of cystalline PV modules. However, the decrease of silicon cell thicknesscouldreduceitsrobustnessundermechanicalloading and potentially leading to the formation of microcracking events within the conductor layer. Hence, the identification of crack initiation and propagation that may trigger performance degradation of PV modules has received a lot of attention in the recent decades. During the fabrication process of silicon solar cells, permanent deformations are induced from the thermomechanicalloadings, which correspondsto the residual stress that leads to cracks. It was pointed out in Rupnowski and Sopori (2009) that 2% of the silicon wafers present defects, causing the increase of production cost and material losses. Even though the fabrication process can be optimized, the imperfections inside the silicon cells during production are unavoidable, especially when the wafer thickness is reduced (Pingel et al. 2009; Liu et al. 2023). Besides, crack formation of silicon cells can also take place during the transportation and installationof PV modules in the field, as wellas during the operation when subjected to the harsh environmental conditions such as wind, snow loading, hail impacts and so on Köntges et al. (2016), Assmus et al. (2011), Liu et al. (2023), Liu et al. (2022). The cracks of silicon solar cells are usually invisible, but could produce over time due to ageing electrically disconnected regions that significantly increase the electrical resistance, and hence reduce the power output of PV modules. Overall, the length, width and orientation of cracks in the silicon wafers directly influence the output of solar panels as pointed out in Munoz et al. (2011), Javvaji et al. (2018), Buerhop et al. (2018). According to the survey of over 200 PV modules, more than 20% of power loss has beendetecteddue to the cellcracksincombinationwith EVA delamination and degradation (Käsewieter et al. 2014). Experimental studies in Chaturvedi et al. (2013) reported that approximately 4% power degradation are caused by cracks of solar cells in the mechanical load tests. When the disconnected cell area caused by crack formation and propagation is greater than 8% of the total cell layer, the associated power loss was found to be roughly proportional to the electrically inactive area (Köntges et al. 2017). Regarding the crack direction in solar cells, several types of cracks can be observed from previous experimental studies, including parallel and perpendicular cracks, +45 and −45 cracks, and multiple direction cracks (Kajari-Schröder et al. 2011), as shown in Fig. 1. Different cracks lead to different power degradation of PV modules. Therefore, it is of utmost importance to develop reliable simulation tools for the fracture modeling of brittle solar cells in order to understand the impact of cracks on the performance of PV modules. In order to accurately predict the failure and crack propagation, several formulations on the basis of different numerical techniques have been developed in the past decades. Simplified phenomenological cohesive zone models specific for the crack modeling along internal boundaries have been proposed in Paggi and Wriggers (2011a), Paggi and Wriggers (2011b), Paggi et al. (2013), Paggi and Wriggers (2012), Infuso et al. (2014), in which cracking events are triggered by the evaluationoftraction-separationconstitutivelaws.This modeling technique has been widely employed to identify crack paths along the element edges by incorporating a characteristic length scale, see the discussions in Ortiz and Pandolfi (1999), Paggi et al. (2013), Mergheim et al. (2005), Mergheim and Steinmann (2006), Spannraft et al. (2023), among many others. In contrast to the cohesive zone models, the extended finiteelement methodhasalsobeendevelopedtomodel crack propagation, and it does not explicitly depend on the finite element discretization corresponding to the physical domain (Moës et al. 1999). Other alternative computational techniques such as enhanced finite ele123 Phase field modeling of anisotropic silicon crystalline cracking ment method (Simo et al. 1993;Oliver et al. 2006) and generalizedfiniteelementmethod (Simone et al. 2006), can also be used to model crack events although possible operative difficulties to identify crack initiation and paths might arise. In order to overcome the disadvantages of the above-mentioned explicit modeling method for complex crack topologies, the phase field approach (Miehe et al. 2010a;Borden et al. 2012) that is based on the Griffith’s theory has been proposed to address the crack modeling of quasi-brittle materials such as silicon solar cells (Paggi et al. 2018;Carollo et al. 2017). In this approach, the sharp crack is generally diffused through the definition of so-called phase field variable and the crack propagation is characterized by the evolution of the corresponding governing equations. According to the -convergence theory in Francfort and Marigo (1998), the crack discontinuities are regularized by a characteristic phase field length scale. Note that the phase field formulation shares several common aspects with the gradient enhanced damage formulations, which is one of the most appealing advantages of this crack modeling method (Frémond and Nedjar 1996;Peerlings et al. 2001;Comi and Perego 2001). In view of this feature, the phase field approach accounting for thermodynamic consistency (Miehe et al. 2010a;Hofacker and Miehe 2013) has been extended to the modeling of ductile fracture (Ambati et al. 2015), anisotropic fracture (Gültekin et al. 2016), multiphysics fracture (Miehe et al. 2015),andpolycrystallinematerials(ClaytonandKnap 2015,2016), among many others. Given the promising aspects of phase field approach, this technique has been employed to model the crack propagation of brittle silicon solar cells in this work. As mentioned above, to reduce the production cost of silicon solar panels, the average thickness of solar cells has decreased from 300 μm to 150 μm during the recent decades (Wohlgemuth et al. 2008;Terheiden et al. 2015;Wang 2006), which increases the breakage rates and crack events. Based on the computational framework developed in Liu et al. (2022), the solid shell formulation combined with anisotropic phase field model at finite deformation is proposed to model the preferential crack propagation in thinfilm solar cells. In the literature, phase field models have been coupled with shell kinematics for the fracture modeling of thin-walled structures through many different ways, see Liu et al. (2022), Reinoso et al. (2017) and references therein. Previous attempts in Ulmer et al. (2012) have been made to model geometrically linear problems according to the Reissner-Midlin theory, but the developments are limited to standard finite elements. Alternative investigation adopting the Local Maximum-Entropy approximation is proposed in Amiri et al. (2014), making it difficult to be implemented into the standard finite element method. The tension-compression split method by integrating the energycontributionthroughtheshellelementthickness has been proposed in Paul et al. (2020), Kiendl et al. (2016), and later extended to the isogeometric modeling of multipatch structures (Proserpio et al. 2020;Paul et al. 2020). Recently, this approach has been adopted in Kikis et al. (2021), Pillai et al. (2020) for phase field fracture modeling of the Reissner-Mindlin shells and plates. Besides, other models have been developed in combination of phase field theory with solid shell and Kirchhoff-Love formulations to model brittle and ductilefracture(Reinosoetal.2017;AmbatiandDeLorenzis 2016;Areias et al. 2016;Proserpio et al. 2021). Hence, the attempt to employ the phase field solid shell formulation, extended to account for the anisotropic fracture orientation, is a proper choice for the crack modeling of solar cells in thin-walled PV laminates. The structure of this work is organized as follows. In Sect.2, the primary aspects of phase field approach for fracture modeling and solid shell kinematics at finite deformation are described in detail. The numericalimplementation and quasi-Newtonmonolithic solution scheme are outlined in Sect.3. Then the numerical examplesare presented in Sect.4,andfinallysomeconcluding remarks are summarized in Sect.5. 2 Phase field solid shell formulation In this section, the phase field solid shell formulation is presented for the modeling of anisotropic crack propagataion in the thin-film brittle solar cells. Firstly, the basics of phase field approach are revisited in Sect.2.1. For a more comprehensive description, the detailed derivation of this theory can be found in the significant work (Miehe et al. 2010a;Miehe et al. 2010b;Bourdin et al. 2000). In Sect.2.2, the solid shell kinetics at finite deformation are described in the convective curvilinear coordinate system. 123 Z. Liu et al. 2.1 Basics of the phase field approach to fracture In the 3D setting, let B0⊂R3and Bt⊂R3denote the reference configuration and current configuration of an arbitrary solid body with existing crack represented by , and Xand xstand for the position vectors in the corresponding configurations, respectively. Furthermore, let ∂B0and ∂Btbe the exterior boundary of the solid bodyin the reference and current configuration, respectively. The motion of the material point Xinside the body is denoted by ϕ(X,t):B0×[0,t]→R3that maps into the corresponding position xin the current configuration during the time interval [0,t]. The deformation gradient Fuis defined as Fu:= ∂Xϕ(X,t), where ∂Xrepresents the partial derivative with respect to the position Xin the reference configuration, and its determinant Ju=det[Fu]denotes the Jacobian. The phase field approach is conceived as the regularization of sharp crack topology within a diffusive crack zone characterized by the scalar-valued function d:B0×[0,t]→[0,1]to model brittle fracture. The phase field variable dsmears out the sharp crack by a diffusive crack area of width l,asshownin Fig. 2. It is noted that the width of regularization area depends on the parameter l, which is called the phase field length scale that controls the transition between intact and damaged parts. In the modeling framework, dis regarded as a smooth function of (X,t), and d=0 and d=1 denote intact and cracked states, respectively. According to the variational theory of fracture (Bourdin et al. 2000), the total energy functional that governs the cracked body under external loadings takes the form of (u,)=int (u,d,)+ext(u)(1) where udenotes the displacement field, and int and ext represent the internal and external energy functionals, respectively. Note that the internal energy functional term int is defined as the sum of elastic energy stored in the solid body Band the energy dissipated through crack propagation . Based on the Griffith’s theory, the internal energy functional int is defined as int (u,d,) =B(u,d)+() =B\ (u,d)d+ Gcd(2) where denotes the specific elastic energy, and Gc is the critical energy release rate. The dissipated fracture energy during the crack propagation is evaluated through the Griffith theory. It should be pointed out that the competition between the elastic energy in the solid body and fracture energy is directly defined in the context of minimization problem. The evaluation of this competition is computationally challenging when the space discretization methods are used to track the crack propagation. To circumvent the use of tracking algorithms, the phase field approach that smears the crack over the whole domain of the body is employed in the same line with gradient damage formulations (Francfort and Marigo 1998;Bourdin et al. 2000;Forest 2009). Inthecontextofthisapproach,thedissipated surface energy can be approximated by a volume integral  Gcd≈B0 Gcγ(d,∇Xd)d(3) where γ(d,∇Xd)is the crack surface density per unit volume, and ∇Xdis the gradient of phase field variable. In this way, the sharp crack is regularized over the body inducing the diffusive representation, as shown in Fig. 2. The crack surface density γ(d,∇Xd)is defined according to the modified Ambrosio-Tortorelli functional (Gültekin et al. 2018) γ(d,∇Xd,ω)=1 2ld2+l 2ω:(∇Xd⊗∇ Xd)(4) where ωis the second order structural tensor characterising the material anisotropy. In order to make the energy release rate orientation-dependent (Clayton and Knap 2015;Nguyen et al. 2017), the structural tensor can be defined as ω=I2+αp(I2−Np⊗Np)(5) where I2is the second-order identity tensor, Nprepresents the unit vector normal to the preferential cleavage plane, and αpis the penalty parameter that is used to prevent damage to develop on planes not normal to Np. According to the previous study (Clayton and Knap 2014), this parameter should be greater than 1.0, and in case of isotropic material, the value is equal to 0.0. The specific elastic energy is influenced by the phase field variable d, which is motivated by the fact that this energy in the damage transition zone has to decrease so as to ensure the thermodynamic equilibrium. Furthermore, to avoid the asymmetric damage behaviour, the elastic energy density can be split into the active contribution +and the passive contribution −(Miehe et al. 2010b). During the reversal loading process, the crack closing precludes the damage evolution,andalsoprovokesthestiffnessdiscovery(Cavuoto 123 Phase field modeling of anisotropic silicon crystalline cracking Fig. 2 Schematic diagram of solid body with (a) sharp crack topology, and (b) diffusive crack topology et al. 2022). Based on the aforementioned viewpoints, the elastic energy can be defined as =g(d)++−(6) where g(d)=(1−d)2+Kdenotes the monotonic degradation function, K≈0 is a positive parameter to ensure numerical stability in case of fully material degradation, and is set to 1.0e−7 in this work. 2.2 Solid shell kinematics In the concept of solid shell theory, the position vectors Xand xin the reference and current configurations, and the phase field variable dat any material point can be approximated by the corresponding vectors on the bottom and top surfaces of solid-like shell element, and can be defined as Xξ1,ξ2,ξ3=1 21+ξ3Xtξ1,ξ2 +1 21−ξ3Xbξ1,ξ2(7a) xξ1,ξ2,ξ3=1 21+ξ3xtξ1,ξ2 +1 21−ξ3xbξ1,ξ2(7b) dξ1,ξ2,ξ3=1 21+ξ3dtξ1,ξ2 +1 21−ξ3dbξ1,ξ2(7c) where the parametric space is identified as: A:= ξ=ξ1,ξ2,ξ3∈R3|−1≤ξi≤+1;i=1,2,3, the subscript b and t represent bottom and top surfaces, respectively, and ξ1,ξ2,ξ3represent the coordinates in the parametric space. The kinematics of solid shell element can be described by the use of the convective curvilinear systemsothattheANSinterpolationforthetransversenormal and shear strain components can be implemented. The covariant tangent vectors Gi(ξ)and gi(ξ)in the reference and current configuration are defined as the partial derivatives of corresponding position vectors X and xwith respect to the convective coordinates ξi Gi(ξ):= ∂X(ξ) ∂ξi,gi(ξ):= ∂x(ξ) ∂ξi,i=1,2,3(8) The contravariant basis vectors can be determined in a standard manner by Gi·Gj=δj iand gi·gj=δj i, and metric tensors are defined as G=GijGi⊗Gj= GijGi⊗Gj,g=gijgi⊗gj=gijgi⊗gj. In curvilinear setting, the deformation gradient Fu is given by Fu=∂x ∂X=gi⊗Gi(9) where the Einstein summation convention on repeated indices is adopted here. Through the definition of metric tensor components Gij =Gi·Gjand gij =gi·gj in the reference and current configuration, the displacement derived Green-Lagrange strain tensor is defined as Eu:= 1 2FuTFu−I2=1 2gij −GijGi⊗Gj (10) The energetically conjugated second Piola-Kirchhoff stress tensor is defined as S=SijGi⊗Gj(11) where Sij represents the contravariant component. 3 Numerical implementation In this section, the numerical strategy rooted in the use of finite element method for the spatial approximation is described briefly. Since this work is restricted to quasi-static analysis, no temporal integration scheme is required, which leads to an equilibrium problem at 123 Z. Liu et al. each pseudo-time step. Firstly, the variational formulation and finite element interpolations are derived in Sect.3.1.Morespecific details regarding the discretizations are ignored here for brevity, but can be found in the previous work (Liu et al. 2022). Secondly, the solution schemes proposed for the phase field solid shell formulation are depicted in Sect.3.2. 3.1 Finite element interpolation ThemixedHu-Washizu variationalprinciple is adopted for the derivation of phase field solid shell formulation incorporatingtheEnhancedAssumedStrain(EAS)and Assumed Natural Strain (ANS) methods to alleviate the locking pathologies. It should be pointed out that the EAS method is employed to remedy volumetric and Poisson thickness locking, while the membrane and inplane locking effects are tackled by the ANS method (Simo and Rifai 1990;Dvorkin and Bathe 1984). Given the enhancement based on the EAS method at finite deformation, the Green-Lagrange strain consists of two parts following the approach proposed in Bischoff and Ramm (1997), including the displacement derived compatible strain Euand incompatible strain ˜ E, and its complete form reads: E=Eu+˜ E. The multi-field variational functional of the solid body takes the form of (S,˜ E,u,d)=B0 Eu,˜ E,dd +B0 Gcl 2d2 l2+ω:(∇Xd⊗∇ Xδd)d−ext ss (12) where ext identifies the external energy functional. Note that the displacement u, the phase field variable d, and the imcompatible strain ˜ Eare the independent variables. Given the orthogonality condition between the stress and strain fields, the stress field is ignored in line with the previous work (Simo and Armero 1992). The first variation of Eq. (12) with respect to the independent fields is given by Ru(u,δu,˜ E,d)=B0 ∂ ∂E:∂Eu ∂uδud −δext(u)=0,∀δu∈Vu(13a) R˜ E(u,˜ E,δ˜ E,d)=B0 ∂ ∂E:δ˜ Ed=0,∀δ˜ E∈V˜ E(13b) Rd(u,˜ E,o,δd)=B0 −2(1−d)δd+d +B0 Gcl1 l2dδd+ω:(∇Xd⊗∇ Xδd)d=0, ∀δd∈Vd(13c) where Vu,V˜ Eand Vdare the admissible spaces of independent variables. According to the concept of isoparametric interpolation, the position vectors Xand xin the reference and current configurations, the displacement vector uand its variation δu, and the phase field variable dand its variation δdcan be approximated by X=N˜ X,x=N˜ x(14a) u=Nd,δu=Nδd(14b) d=Nd¯ d,δd=Ndδ¯ d(14c) where ˜ Xand ˜ xare the corresponding nodal position vectors, ddenotes the nodal displacement vector, and ¯ drepresents the nodal phase field vector. The shape function matrix Ndis defined as Nd=[N1,N2,N3,N4,N5,N6,N7,N8](15) where its component NIis given by NI=1 81+ξ1 Iξ11+ξ2 Iξ21+ξ3 Iξ3(16) with ξ1 I,ξ2 I,ξ3 I=±1. The material gradient of phase field ∇Xdand its variation ∇Xδdcan be defined as ∇Xd=G−T∇ξ¯ d=Bd(ξ)¯ d, ∇Xδd=G−T∇ξδ¯ d=Bd(ξ)δ¯ d(17) where ∇ξdenotes the gradient with respect to the natural coordinates. The vector form of the Green-Lagrange strain tensorisgivenbyE=[E11,2E12,2E13,E22,2E23,E33]T. To alleviate the curvature thickness locking effects, the ANS method proposed in Betsch and Stein (1995)is adopted to modify the transverse normal strain component E33 by four collocation points defined in convective coordinates ξCias ξC1=(−1,−1,0),ξC2= (1,−1,0),ξC3=(1,1,0), and ξC4=(−1,1,0),see Fig. 3. Besides, to prevent transverse shear locking, the ANS method proposed in Dvorkin and Bathe (1984) is also employed in this work. The four collocation points for the enhancement of transverse shear strain components are ξA1=(0,−1,0),ξA2=(0,1,0), ξB1=(−1,0,0), and ξB2=(1,0,0), see Fig. 3.Given 123 Phase field modeling of anisotropic silicon crystalline cracking Fig. 3 Position of collocation points in the element parametric space for the ANS method Fig. 4 Schematic of the benchmark problem for anisotropic fracture modeling the aforementioned ANS interpolations, the strain vector is given by Eu= ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 1 2(g11 −G11) (g12 −G12) 1−ξ2gA1 13 −GA1 13 +1+ξ2gA2 13 −GA2 13  1 2(g22 −G22) 1−ξ1gB1 23 −GB1 23 +1+ξ1gB2 23 −GB2 23  4 i=11 41+ξ1 iξ11+ξ2 iξ21 2gCi 33 −GCi 33 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ (18) where the superscripts A1,A 2,B 1,B 2, and Ciwith i=1,2,3,4 stand for the corresponding collocation points. The approximation of virtual strain is given by δEu=Bδd with B=[B1,B2,B3,B4,B5,B6,B7,B8](19) where BIis the interpolation matrix for each node I. According to the previous work (Klinkel and Wagner 1997;Vu-Quoc and Tan 2003), the interpolations of incompatible strain vector ˜ Eand its variation δ˜ Eare given by ˜ E≈M(ξ)ζ,δ ˜ E≈M(ξ)δζ(20) whereM(ξ) istheinterpolationmatrixoftheenhancing modes ζ, and is defined as M(ξ) =det J0 det JT−T 0˜ M(ξ) (21) whereJ=[G1,G2,G3]T,J0isitsevaluationattheelement center, T0is the tranformation matrix, and ˜ M(ξ) is the interpolation matrix in the parametric space (Liu et al. 2022). 123 Z. Liu et al. Fig. 5 The phase field contour plots of fully cracked specimen with the preferential crack orientation of Grain 1 equal to 0 degree, 22.5 degree, 45 degree, and 67.5 degree Fig. 6 The obtained force vs. displacement curves with the preferential crack orientation of Grain 1 equal to 0 degree, 22.5 degree, 45 degree, and 67.5 degree Inserting the aforementioned interpolations into Eqs. (13), the discretes residual equations are given by Rd(d,δd,¯ d,ζ)=B0 g(d)B(d)TSd−Rd ext (22a) Rζ(d,¯ d,ζ,δζ)=B0 g(d)M(ξ)TSd(22b) Rd(d,d,δd,ζ)=B0 −2(1−d)N(ξ)T+d +B0 GclBdTW∇Xd+1 l2N(ξ)Tdd (22c) To avoid the irreversible growth of the fracture process, a history variable His introduced to modify the residual vector Rd(Msekh et al. 2015), which is transformed into 123 Phase field modeling of anisotropic silicon crystalline cracking Fig. 7 Schematic diagram of the local fracture modeling of one single solar cell when the photovoltaic module is subjected to tensile loading Rd=B0 −2(1−d)N(ξ)THd +B0 GclBdTW∇Xd+1 l2N(ξ)Tdd (23) Inviewoftheirrevisibility(Mieheetal.2010b),thehistory variable Hmust obey the following Kuhn-Tucker conditions +−H⩽0,˙ H⩾0,˙ H(+−H)=0(24) At the certain time point t, the history variable Hcan be expressed as H=max τ∈[0,t]+(τ) (25) Tosolvethe multi-field problem, an iterativescheme is adopted, and the consistent linearization of Eqs. (22) obtained from the concept of Gateaux directional derivative can be derived as ⎡ ⎣ kdd kdζ0 kζdkζζ 0 00k dd ⎤ ⎦⎡ ⎣ d ζ ¯ d ⎤ ⎦=⎡ ⎣ Rd ext 0 0 ⎤ ⎦−⎡ ⎣ Rd Rζ Rd ⎤ ⎦(26) It is worth mentioning that the derivation of stiffness terms are omitted here for brevity, see (Liu et al. 2022). Since inter-element continuity of the enhanced strain field is not required, the second equation of Eqs. (26) can be condensed out in the element level in line with (Bischoff and Ramm 1997), and the condensed equations can be expressed as k∗ dd 0 0k dd d d=Rd ext 0−Rd∗ Rd(27) Table 1 Mechanical properties for photovoltaic modules (Paggi et al. 2011;Corrado et al. 2017) E(GPa)Density (kg/m3) Backsheet 2.8 1200 EVA 0.01 1180 Glass 73 2500 where the modified stiffness k∗ dd and residual vector Rd∗are defined as k∗ dd =kdd −kdζk−1 ζζ kζd(28a) Rd∗=Rd−kdζk−1 ζζ Rζ(28b) 3.2 Solution schemes In order to solve the nonlinear iterative equations, two different schemes are widely employed for the coupled phase field-displacement problem, including the monolithic and staggered schemes. Monolithic scheme retains unconditional stability, but suffers from poor convergence performance that hinders its wide application. On the other hand, the staggered scheme is very robust and can overcome the convergence issue. However, very small time increment must be adopted to prevent the solution deviating from the equilibrium, so this solution scheme is very time-consuming. In the following, the quasi-Newton monolithic scheme with improved performance in terms of both convergence and computational efficiency is introduced here for the 123 Z. Liu et al. field contour plots of the local model with the preferential crack orientation angle equal to 45 degree, which presents completely different crack patterns compared to the former. Note that when the preferential crack orientation is set to 0 degree, the silicon solar cell is fully cracked when the loading displacement imposed on the global photovoltaic module reaches 76.5 mm, while in the case with the preferential crack orientation equal to 45 degree, the silicon solar cell is not fully cracked until when the loading displacement reaches 109.5 mm. It can be concluded here that the orientation of silicon cell will significantly influence the crack growth when the photovoltaic module is subjected to bending loading, which could provide guidance to the photovoltaic industry for the better design of photovoltaic products with excellent crack-resistant performance. 5 Concluding remarks The present modeling strategy constitutes a major progress with respect to the state of the art, since it provides the first proof of concept of a computational methodology integrating structural mechanics considerations for real photovoltaic installations and advanced fracture mechanics models for the assessment of damage distributions in solar cells. To accurately capture the structural deformation of very thin silicon solar cells, the solid shell element is formulated at large deformation, which is incorporated into the phase field approach for the cracking modeling using the efficient and robust quasi-Newton solution. Note that the fracture anisotropy is also taken into account in the phase field solid shell formulation so that experimental crack patterns with varying orientation can be reproduced. For the sake of reducing computational cost of phase field fracture modeling, a global local approach suitable for the local fracture modeling of solar cells in the global level of photovoltaic module is explored, which is demonstrated by the simulation of severaldifferentloadingcases.Giventhe complexityof loading scenarios of photovoltaic modules in the outdoor environment, this global local fracture modeling method can be very promising for the possible realistic prediction of crack growth of silicon solar cells. Acknowledgements This paper is dedicated to the memory of Prof. Dominique Leguillon, who was a close friend and an extraordinary researcher Funding Funding for open access publishing: Universidad de Sevilla/CBUA The authors acknowledge funding received from the European Union’s H2020-MSCA-ITN-2019 research and innovation program under the Marie Skłodowska-Curie grant agreement no. 861061 - Project NEWFRAC “New strategies for multifield fracture problems across scales in heterogeneous systems for Energy, Health and Transport”. J. Reinoso acknowledges the support of Ministerio de Ciencia e Innovación de España through the Project TED2021-131649B-I00 Data Availability No datasets were generated or analysed during the current study. Declarations Conflict of interest The authors declare no conflict of interest. Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or otherthirdpartymaterialinthisarticleareincludedinthearticle’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/ by/4.0/. References Ambati M, De Lorenzis L (2016) Phase-field modeling of brittle and ductile fracture in shells with isogeometric NURBSbased solid-shell elements. Comput Methods Appl Mech Eng 312:351–373 Ambati M, Gerasimov T, De Lorenzis L (2015) Phase-field modeling of ductile fracture. Comput Mech 55:1017–1040 Amiri F, Millán D, Shen Y, Rabczuk T, Arroyo M (2014) Phasefield modeling of fracture in linear thin shells. Theoret Appl Fract Mech 69:102–109 Areias P, Rabczuk T, Msekh M (2016) Phase-field analysis of finite-strain plates and shells including element subdivision. Comput Methods Appl Mech Eng 312:322–350 Assmus M, Jack S, Weiss K-A, Koehl M (2011) Measurement and simulation of vibrations of PV-modules induced by dynamic mechanical loads. Prog Photovolt Res Appl 19(6):688–694 Betsch P, Stein E (1995) An assumed strain approach avoiding artificial thickness straining for a non-linear 4-node shell element. Commun Numer Methods Eng 11(11):899–909 Bischoff M, Ramm E (1997) Shear deformable shell elements for large strains and rotations. Int J Numer Meth Eng 40(23):4427–4449 Borden MJ, Verhoosel CV, Scott MA, Hughes TJ, Landis CM (2012) A phase-field description of dynamic brittle fracture. Comput Methods Appl Mech Eng 217:77–95 123 Phase field modeling of anisotropic silicon crystalline cracking Bourdin B, Francfort GA, Marigo J-J (2000) Numerical experiments in revisited brittle fracture. J Mech Phys Solids 48(4):797–826 Buerhop C, Wirsching S, Bemm A, Pickel T, Hohmann P, Nieß M, Vodermayer C, Huber A, Glück B, Mergheim J et al (2018) Evolution of cell cracks in PV-modules under field and laboratory conditions. Prog Photovolt Res Appl 26(4):261–272 Carollo V, Reinoso J, Paggi M (2017) A 3d finite strain model for intralayer and interlayer crack simulation coupling the phase field approach and cohesive zone model. Compos Struct 182:636–651 Cavuoto R, Lenarda P, Misseroni D, Paggi M, Bigoni D (2022) Failurethroughcrack propagationin componentswithholes and notches: an experimental assessment of the phase field model. Int J Solids Struct 257:111798 Chaturvedi P, Hoex B, Walsh TM (2013) Broken metal fingers in silicon wafer solar cells and PV modules. Sol Energy Mater Sol Cells 108:78–81 Clayton JD, Knap J (2014) A geometrically nonlinear phase field theory of brittle fracture. Int J Fract 189(2):139–148 ClaytonJ, KnapJ(2015) Phasefieldmodeling of directionalfracture in anisotropic polycrystals. Comput Mater Sci 98:158– 169 Clayton J, Knap J (2016) Phase field modeling and simulation of coupled fracture and twinning in single crystals and polycrystals. Comput Methods Appl Mech Eng 312:447–467 Comi C, Perego U (2001) Fracture energy based Bi-dissipative damage model for concrete. Int J Solids Struct 38(36– 37):6427–6454 Corrado M, Infuso A, Paggi M (2017) Simulated hail impacts on flexible photovoltaic laminates: testing and modelling. Meccanica 52:1425–1439 Dvorkin EN, Bathe K-J (1984) A continuum mechanics based four-node shell element for general non-linear analysis, Eng Comput Forest S (2009) Micromorphic approach for gradient elasticity, viscoplasticity, and damage. J Eng Mech 135(3):117–131 Francfort GA, Marigo J-J (1998) Revisiting brittle fracture as an energy minimization problem. J Mech Phys Solids 46(8):1319–1342 Frémond M, Nedjar B (1996) Damage, gradient of damage and principle of virtual power. Int J Solids Struct 33(8):1083– 1103 Gültekin O, Dal H, Holzapfel GA (2016) A phase-field approach to model fracture of arterial walls: theory and finite element analysis. Comput Methods Appl Mech Eng 312:542–566 Gültekin O, Dal H, Holzapfel GA (2018) Numerical aspects of anisotropic failure in soft biological tissues favor energybased criteria: a rate-dependent anisotropic crack phasefield model. Comput Methods Appl Mech Eng 331:23–52 Hofacker M, Miehe C (2013) A phase field model of dynamic fracture: Robust field updates for the analysis of complex crack patterns. Int J Numer Meth Eng 93(3):276–301 Infuso A, Corrado M, Paggi M (2014) Image analysis of polycrystalline solar cells and modelling of intergranular and transgranular cracking. J Eur Ceram Soc 34(11):2713–2722 Javvaji B, Budarapu PR, Paggi M, Zhuang X, Rabczuk T (2018) Fracture properties of graphene-coated silicon for photovoltaics. Adv Theor Simul 1(12):1800097 Kajari-Schröder S, Kunze I, Eitner U, Köntges M (2011) Spatial and orientational distribution of cracks in crystalline photovoltaic modules generated by mechanical load tests. Sol Energy Mater Sol Cells 95(11):3054–3059 Käsewieter J, Haase F, Larrodé MH, Köntges M (2014) Cracks in solar cell metallization leading to module power loss under mechanical loads. Energy Procedia 55:469–477 Kiendl J, Ambati M, De Lorenzis L, Gomez H, Reali A (2016) Phase-fielddescriptionofbrittlefractureinplatesandshells. Comput Methods Appl Mech Eng 312:374–394 Kikis G, Ambati M, De Lorenzis L, Klinkel S (2021) Phasefield model of brittle fracture in reissner-mindlin plates and shells. Comput Methods Appl Mech Eng 373:113490 Klinkel S, Wagner W (1997) A geometrical non-linear brick element based on the EAS-method. Int J Numer Meth Eng 40(24):4529–4545 Köntges M, Siebert M, Morlier A, Illing R, Bessing N, Wegert F (2016) Impact of transportation on silicon wafer-based photovoltaic modules. Prog Photovolt Res Appl 24(8):1085– 1095 Köntges M, Oreski G, Jahn U, Herz M, Hacke P, Weiß K.- A. (2017) Assessment of Photovoltaic Module Failures in the Field: International Energy Agency Photovoltaic Power Systems Programme: IEA PVPS Task 13, Subtask 3: Report IEA-PVPS T13-09: 2017, International Energy Agency Liu Z, Reinoso J, Paggi M (2022) Phase field modeling of brittle fracture in large-deformation solid shells with the efficient quasi-newton solution and global-local approach. Comput Methods Appl Mech Eng 399:115410 Liu Z, Reinoso J, Paggi M (2022) A humidity dose-CZM formulation to simulate new end-of-life recycling methods for photovoltaic laminates. Eng Fract Mech 259:108125 Liu Z, Marino M, Reinoso J, Paggi M (2023) A continuum largedeformation theory for the coupled modeling of polymersolvent system with application to PV recycling. Int J Eng Sci 187:103842 Liu Z, Lenarda P, Reinoso J, Paggi M (2023) A multifield coupled thermo-chemo-mechanical theory for the reactiondiffusion modeling in photovoltaics. Int J Numer Meth Eng 124(12):2876–2901 Liu Z, Reinoso J, Paggi M (2022) Hygro-thermo-mechanical modeling of thin-walled photovoltaic laminates with polymeric interfaces. J Mech Phys Solids 169:105056 Matthies H, Strang G (1979) The solution of nonlinear finite element equations. Int J Numer Meth Eng 14(11):1613– 1626 Mergheim J, Steinmann P (2006) A geometrically nonlinear Fe approach for the simulation of strong and weak discontinuities. Comput Methods Appl Mech Eng 195(37–40):5037– 5052 Mergheim J, Kuhl E, Steinmann P (2005) A finite element method for the computational modelling of cohesive cracks. Int J Numer Meth Eng 63(2):276–289 Miehe C, Hofacker M, Welschinger F (2010) A phase field model for rate-independent crack propagation: Robust algorithmic implementation based on operator splits. Comput Methods Appl Mech Eng 199(45–48):2765–2778 Miehe C, Welschinger F, Hofacker M (2010) Thermodynamically consistent phase-field models of fracture: variational principles and multi-field Fe implementations. Int J Numer Meth Eng 83(10):1273–1311 123 Z. Liu et al. Miehe C, Schaenzel L.-M, Ulmer H (2015) 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, Computer Methods in Applied Mechanics and Engineering 294 449–485 Moës N, Dolbow J, Belytschko T (1999) A finite element method for crack growth without remeshing. Int J Numer Meth Eng 46(1):131–150 Msekh MA, Sargado JM, Jamshidian M, Areias PM, Rabczuk T (2015) Abaqus implementation of phase-field model for brittle fracture. Comput Mater Sci 96:472–484 Munoz M, Alonso-García MC, Vela N, Chenlo F (2011) Early degradation of silicon PV modules and guaranty conditions. Sol Energy 85(9):2264–2274 Nguyen T-T, Réthoré J, Yvonnet J, Baietto M-C (2017) Multiphase-field modeling of anisotropic crack propagation for polycrystalline materials. Comput Mech 60:289–314 Oliver J, Huespe AE, Blanco S, Linero D (2006) Stability and robustness issues in numerical modeling of material failure with the strong discontinuity approach. Comput Methods Appl Mech Eng 195(52):7093–7114 Ortiz M, Pandolfi A (1999) Finite-deformation irreversible cohesive elements for three-dimensional crack-propagation analysis. Int J Numer Meth Eng 44(9):1267–1282 Paggi M, Wriggers P (2011) A nonlocal cohesive zone model for finite thickness interfaces-Part I: mathematical formulation and validation with molecular dynamics. Comput Mater Sci 50(5):1625–1633 Paggi M, Wriggers P (2011) A nonlocal cohesive zone model for finite thickness interfaces-Part II: Fe implementation and application to polycrystalline materials. Comput Mater Sci 50(5):1634–1643 PaggiM,WriggersP (2012) Stiffnessandstrengthof hierarchical polycrystalline materials with imperfect interfaces. J Mech Phys Solids 60(4):557–572 Paggi M, Kajari-Schröder S, Eitner U (2011) Thermomechanical deformations in photovoltaic laminates. J Strain Anal Eng Design 46(8):772–782 Paggi M, Lehmann E, Weber C, Carpinteri A, Wriggers P, SchaperM (2013)Anumerical investigationof theinterplay between cohesive cracking and plasticity in polycrystalline materials. Comput Mater Sci 77:81–92 Paggi M, Corrado M, Rodriguez MA (2013) A multi-physics and multi-scale numerical approach to microcracking and power-loss in photovoltaic modules. Compos Struct 95:630–638 Paggi M, Corrado M, Reinoso J (2018) Fracture of solar-grade anisotropic polycrystalline silicon: a combined phase fieldcohesive zone model approach. Comput Methods Appl Mech Eng 330:123–148 PaggiM,CorradoM, BerardoneI(2016)Aglobal/localapproach for the prediction of the electric response of cracked solar cellsinphotovoltaicmodulesundertheactionofmechanical loads. Eng Fract Mech 168:40–57 Papargyri L, Theristis M, Kubicek B, Krametz T, Mayr C, Papanastasiou P, Georghiou GE (2020) Modelling and experimental investigations of microcracks in crystalline silicon photovoltaics: a review. Renew Energy 145:2387– 2408 Paul K, Zimmermann C, Mandadapu KK, Hughes TJ, Landis CM, Sauer RA (2020) An adaptive space-time phase field formulation for dynamic fracture of brittle shells based on LR NURBS. Comput Mech 65(4):1039–1062 Paul K, Zimmermann C, Duong TX, Sauer RA (2020) Isogeometric continuity constraints for multi-patch shells governed by fourth-order deformation and phase field models. Comput Methods Appl Mech Eng 370:113219 Peerlings R, Geers M, De Borst R, Brekelmans W (2001) A critical comparison of nonlocal and gradient-enhanced softening continua. Int J Solids Struct 38(44–45):7723–7746 Pillai U, Triantafyllou SP, Ashcroft I, Essa Y, de la Escalera FM (2020) Phase-field modelling of brittle fracture in thin shell elements based on the MITC4+ approach. Comput Mech 65(6):1413–1432 Pingel S, Zemen Y, Frank O, Geipel T, Berghold J (2009) Mechanical stability of solar cells within solar panels, Proc. of 24th EUPVSEC 3459–3464 Proserpio D, Ambati M, De Lorenzis L, Kiendl J (2020) A framework for efficient is geometric computations of phase-field brittle fracture in multipatch shell structures. Comput Methods Appl Mech Eng 372:113363 Proserpio D, Ambati M, De Lorenzis L, Kiendl J (2021) Phasefield simulation of ductile fracture in shell structures. Comput Methods Appl Mech Eng 385:114019 Reinoso J, Blázquez A (2016) Application and finite element implementation of 7-parameter shell element for geometrically nonlinear analysis of layered CFRP composites. Compos Struct 139:263–276 Reinoso J, Paggi M, Linder C (2017) Phase field modeling of brittle fracture for enhanced assumed strain shells at large deformations: formulation and finite element implementation. Comput Mech 59(6):981–1001 RupnowskiP,SoporiB(2009)Strengthofsiliconwafers:fracture mechanics approach. Int J Fract 155:67–74 Sander M, Dietrich S, Pander M, Ebert M, Bagdahn J (2013) Systematic investigation of cracks in encapsulated solar cells after mechanical loading. Sol Energy Mater Sol Cells 111:82–89 Simo J-C, Armero F (1992) Geometrically non-linear enhanced strain mixed methods and the method of incompatible modes. Int J Numer Meth Eng 33(7):1413–1449 Simo JC, Rifai M (1990) A class of mixed assumed strain methods and the method of incompatible modes. Int J Numer Meth Eng 29(8):1595–1638 Simo JC, Oliver J, Armero F (1993) An analysis of strong discontinuities induced by strain-softening in rate-independent inelastic solids. Comput Mech 12(5):277–296 Simone A, Duarte CA, Van der Giessen E (2006) A generalized finite element method for polycrystals with discontinuous grain boundaries. Int J Numer Meth Eng 67(8):1122–1145 Spannraft L, Steinmann P, Mergheim J (2023) A generalized anisotropicdamageinterfacemodelforfinitestrains.JMech Phys Solids 174:105255 Terheiden B, Ballmann T, Horbelt R, Schiele Y, Seren S, Ebser J, Hahn G, Mertens V, Koentopp MB, Scherff M et al (2015) Manufacturing 100-μm-thick silicon solar cells with efficiencies greater than 20% in a pilot production line. Phys Status Solidi (a) 212(1):13–24 Ulmer H, Hofacker M, Miehe C (2012) Phase field modeling of fracture in plates and shells. PAMM 12(1):171–172 123 Phase field modeling of anisotropic silicon crystalline cracking Vu-Quoc L, Tan X (2003) Optimal solid shells for non-linear analyses of multilayer composites. I. statics. Comput Methods Appl Mech Eng 192(910):975–1016 Wang PA (2006) Industrial challenges for thin wafer manufacturing, in: 2006 IEEE 4th World Conference on Photovoltaic Energy Conference, Vol. 1, IEEE, pp. 1179–1182 WohlgemuthJH,CunninghamDW,PlacerNV, KellyGJ,Nguyen AM, The effect of cell thickness on module reliability, in, (2008) 33rd IEEE Photovoltaic Specialists Conference. IEEE 2008:1–4 Wu J-Y, Huang Y, Nguyen VP (2020) On the BFGS monolithic algorithmfor the unified phase field damage theory. Comput Methods Appl Mech Eng 360:112704 Xu X, Lai D, Wang G, Wang Y (2021) Nondestructive silicon wafer recovery by a novel method of solvothermal swelling coupled with thermal decomposition. Chem Eng J 418:129457 Publisher’s Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations. 123