Hygro-thermo-mechanical modeling of thin-walled photovoltaic laminates with polymeric interfaces

J. Mech. Phys. Solids 169 (2022) 105056 Available online 5 September 2022 0022-5096/© 2022 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/). Contents lists available at ScienceDirect Journal of the Mechanics and Physics of Solids journal homepage: www.elsevier.com/locate/jmps Hygro-thermo-mechanical modeling of thin-walled photovoltaic laminates with polymeric interfaces Z. Liu a,b,∗, J. Reinosob, M. Paggi a aIMT School for Advanced Studies Lucca, Piazza San Francesco 19, 55100, Lucca, Italy bDepartamento de Mecánica de Medios Continuos y Teoría de Estructuras, School of Engineering, Universidad de Sevilla, Camino de los Descubrimientos s/n, 41092, Seville, Spain ARTICLE INFO Keywords: Hygro-thermo-mechanical coupling Thermo-visco-elastic behavior Moisture diffusion Finite element method ABSTRACT A three-dimensional hygro-thermo-mechanical computational framework for photovoltaic (PV) laminates as well as its numerical implementation are established in this work. Aiming at an efficient thermo-mechanical modeling of thin-walled structures with polymeric interfaces, the solid shell element, which incorporates the enhanced assumed strain (EAS) method and the assumed natural strain (ANS) method for the alleviation of locking pathologies, and the interface element with thermo-visco-elastic cohesive zone model using fractional calculus are formulated. Besides, the finite element (FE) implementation of moisture diffusion in the 3D setting along the polymeric interfaces is also derived with the consideration of spatial and temporal variation of diffusivity due to its temperature and material decohesion dependencies. Given the difference between the time scales of moisture diffusion and thermo-mechanical problems, a staggered scheme is proposed for the solution of the coupled hygro-thermo- mechanical governing equations. Specifically, the relative displacement and temperature fields are firstly solved from the thermo-mechanical analysis, and then projected to the FE model of moisture diffusion to determine the diffusion coefficient for its subsequent analysis. The proposed method is applied to the simulation of three international standard tests of PV modules, namely the damp heat test, the humidity freeze test, and the thermal cycling test, and numerical predictions are compared with analytical solution for the damp heat case with a constant temperature boundary condition, as well as experimental electroluminescence (EL) images obtained from the thermal cycling test with the cyclic temperature boundary condition. A very satisfactory consistency demonstrates the effectiveness and reliability of this modeling framework. 1. Introduction Typical photovoltaic modules are laminates consisting of different layers, including tempered glass, encapsulant layers (usually ethylene vinyl acetate, EVA), a layer of interconnected silicon solar cells, and backsheet, see Fig. 1. The layer of silicon solar cells is used to produce electric energy, while the other layers like glass and backsheet play the role of protecting the modules from the outdoor harsh environments (Deng et al.,2021). The durability analysis of PV modules has received a great deal of attention in the past decades (Liu et al.,2022a;Gagliardi and Paggi,2019;Antartis and Chasiotis,2014;Lenarda and Paggi,2016), and the modeling of polymeric materials is quite challenging as it requires a multiphysics framework to accurately predict their overall performance and degradation (Paggi et al.,2013;Bahrololoumi et al.,2020,2021;Silberstein and Boyce,2010,2011). ∗Corresponding author at: IMT School for Advanced Studies Lucca, Piazza San Francesco 19, 55100, Lucca, Italy. E-mail address: [email protected] (Z. Liu). https://doi.org/10.1016/j.jmps.2022.105056 Received 25 April 2022; Received in revised form 25 July 2022; Accepted 30 August 2022 Journal of the Mechanics and Physics of Solids 169 (2022) 105056 2 Z. Liu et al. Fig. 1. The general structure of a typical PV module (Liu et al.,2022a). Fig. 2. EL images of PV modules with moisture induced dimmer areas during the damp heat test (Herrmann and Bogdanski,2011). The EVA layer of PV modules is usually made of polymeric materials, and permeable to moisture, which diffuses from the edges, backsheet and interspaces between different solar cells. This moisture diffusion can induce chemical oxidation of grid lines deposited on the surface of solar cells and subsequently lead to degradation and power loss, which has been reported in the damp heat test of PV modules exposed to the very aggressive environment with 85 ◦C temperature and 85% relative humidity (Iseghem et al.,2012; Köntges et al.,2017,2014). As shown in Fig. 2, the dark bands can be observed at the edges of solar cells in the captured EL image after different elapsed time, which were not detected in the beginning but became evident during ageing. It is important to remark that in photovoltaics, the larger the dark – electrically inactive – area, the higher the electrical power loss of the module, see Paggi et al. (2016). Particularly, in the damp heat accelerated ageing test, the dark bands relate to the moisture ingress into the module from the edges of the solar cells over time. Since moisture can induce chemical reactions inside the modules and delamination failure between different laminae, it is crucial to establish a reliable and accurate modeling strategy of moisture diffusion along the EVA layer. The complexity regards the strong dependency of moisture diffusion coefficient of EVA on the temperature, which indicates the coupling between the moisture and thermal fields, as shown in Kempe (2006). On the other hand, moisture diffusion also degrades the cohesive energy of EVA layer, giving rise to the delamination failure between the different laminae (Liu et al., 2019), such as the decohesion between solar cell and backsheet or tempered glass (Novoa et al.,2014,2016), which corresponds to the coupling between the mechanical and moisture fields (Dunlop,2003). Regarding the thermo-mechanical behavior, the EVA polymer shows a strong thermo-visco-elastic constitutive response, and its Young’s modulus varies within a range by three orders of magnitude depending on temperature, as experimentally reported in Eitner et al. (2011), Paggi et al. (2011). In order to approximate the power-law trend from experimental observation, the generalized Maxwell rheological models can be used to determine the relaxation modulus with exponential type equations. However, a great number of elements as well as model parameters have to be taken into account, which requires laborious calibration work. To simplify the procedure of parameter identification, the fraction calculus method has been proved to be very effective for the modeling of visco-elastic constitutive behavior (Mainardi,2010;Di Paola et al.,2013,2011). Thus this formulation has been adopted to describe the rheologically complex behavior of polymeric EVA, whose microstructure changes with temperature, and, according to Paggi and Sapora (2015), only two temperature dependent parameters are required for its complete description. To model these coupled nonlinear hygro-thermo-mechanical problems in PV modules, a comprehensive computational framework in the three-dimensional space, where the coupled thermo-mechanical problem and the moisture diffusion are solved in a staggered manner, is required so that the dependency of diffusion properties of EVA on the thermo-mechanical fields can be accounted for Journal of the Mechanics and Physics of Solids 169 (2022) 105056 3 Z. Liu et al. Fig. 3. The proposed solution scheme for the 3D hygro-thermo-mechanical modeling framework. properly. Analytical solution for the moisture diffusion along EVA layer proposed in Kempe (2006) assumes constant diffusivity without any update based on the actual temperature, and the spatial variation due to temperature dependency is ignored. Although it might be feasible in the steady state temperature case, such as in the damp heat test, its validity for cyclic temperature boundary conditions in the humidity freeze and in the thermal cycling tests is hard to guarantee. Besides, since moisture can also diffuse from the channel cracks and interspaces between the solar cells, it is almost impossible to obtain analytical solution in these complex cases. To overcome the limitations of analytical solution methods, the finite element modeling methodology in the three-dimensional setting is proposed to simulate the coupled problems in this work. The thermo-mechanical phenomena in the PV laminate, which is much faster than moisture diffusion, is firstly solved through a monolithic fully implicit solution scheme, see Fig. 3. The calculated temperature and displacement fields inside the EVA layer are then projected to the nodes of another finite element model specific for moisture diffusion, so that the diffusion properties at the corresponding time and space can be determined accordingly. With updated diffusion coefficient from thermo-mechanical solution at each time increment, the moisture diffusion can be solved by the Newton–Raphson scheme in a standard way. For the sake of computational efficiency, the EVA layers are modeled using zero-thickness interface finite elements with a thermovisco-elastic cohesive zone model based on fractional calculus, and a thermo-mechanical solid shell formulation incorporating the EAS and ANS methods to remedy potential locking pathologies, which is kinematically compatible with the 3D interface element, is established here for the modeling of the different thin-walled laminae in PV modules such as glass layer, backsheet and silicon solar cell layer. Modeling thermo-mechanical responses in thin-walled structures necessitates the consistent derivation of formulations accounting for the stress and heat transfer across the internal interfaces. To use the 3D constitutive laws without further modifications, recent advances in shell element formulation aim at incorporating the three-dimensional effects into the corresponding implementation. To achieve this, it is pointed out in Rah et al. (2013), Vu-Quoc and Tan (2003) that a linear normal strain distribution through the thickness direction has to be embodied in the shell element formulation. In this regard, two paths have been followed: (i) shell formulations that include the linear distribution using either quadratic displacement distribution or enhanced strain methods through the reference surface of body in the thickness direction (Başar and Ding,1997;Betsch et al.,1996; Brank,2005;Büchter et al.,1994;Parisch,1995;Simo et al.,1990;Liu and Xia,2021), and (ii) shell formulations relying on the solid shell concept through the parametrization of the top and bottom surfaces of the body (Miehe,1998;Reinoso and Blázquez, 2016;Harnau and Schweizerhof,2002;Klinkel and Wagner,1997;Schwarze and Reese,2009;Liu and Xia,2019;Liu et al.,2022b). The latter has been extensively developed in the past decades since the complex update procedure regarding the rotation tensor can be completely avoided. However, the extension of solid shell formulation to thermo-mechanical applications has received very limited attention, see Braun et al. (1994), Bathe and Dvorkin (1985), Bischoff and Ramm (1997), Hauptmann et al. (2000). Besides, to alleviate the locking pathologies in shell formulations complying with the low-order kinematic interpolation, different numerical strategies have been proposed such as the EAS method (Korelc and Wriggers,1996;Simo et al.,1993;Simo and Rifai,1990;Valente et al.,2003), the ANS method (Betsch and Stein,1995;Dvorkin and Bathe,1984), and combination of them (Cardoso et al.,2008; Gruttmann and Wagner,2006;Kasper and Taylor,2000). This work is structured as follows. In Section 2, the primary aspects of coupled thermo-mechanical analysis along with the definition of kinematics and constitutive formulation, the thermo-visco-elastic cohesive zone model through the fractional calculus method, as well as temperature and gap dependent moisture diffusion are presented in detail. The weak forms of the governing equations for the hygro-thermo-mechanical problems and the corresponding finite element discretization are given in Section 3. In particular, the Hu–Washizu variational principle with thermo-mechanical solid shell formulation accounting for the EAS and ANS methods to alleviate locking effects is outlined in Section 3.1, while Section 3.2 and Section 3.3 cover the 3D thermomechanical interface and moisture diffusion finite element implementation, respectively. Details of the numerical algorithm for the implementation of the hygro-thermo-mechanical framework with staggered solution scheme are provided in Section 4. Then Journal of the Mechanics and Physics of Solids 169 (2022) 105056 4 Z. Liu et al. Fig. 4. The finite thermo-mechanical deformation of a three-dimensional body from the reference configuration to the current configuration. this modeling methodology is applied to simulate the three international standard tests of PV laminates, and validated with the analytical solution as well as experimental data, which is presented in Section 5. Finally, some concluding remarks of the proposed modeling framework are drawn in Section 6. 2. Hygro-thermo-visco-elastic modeling framework for the PV laminate In this section, the modeling framework for the thermo-mechanical coupling in thin-walled PV laminae separated by thermovisco-elastic interfaces and moisture diffusion along the polymeric layers is presented. Moisture diffusion takes place in the polymeric encapsulant layers of laminate, which usually percolates from the free edges and towards the center of solar cell panels. The encapsulant layers are made of visco-elastic polymer materials and permeable to moisture, which is one of the dominant factors leading to the electrical degradation of PV modules during the service life. Besides, moisture diffusion can also degrade the adhesion strength of encapsulant layers, causing delamination failure between silicon solar cells and backsheet or the glass superstrate. In order to effectively model the cohesive mechanical behavior of encapsulant layers, a 3D cut-off traction–separation law is formulated. 2.1. Coupled thermo-mechanical kinematics and formulation For thermo-mechanical boundary value problems in solids, displacement and absolute temperature fields are usually chosen as the independent fields in the coupled governing equations. In the 3D setting, let 0⊂R3denotes the reference configuration, while 𝑡⊂R3denotes the current configuration. The position vectors of material points in the reference and current configurations are represented by 𝐗and 𝐱, respectively. The displacement field is denoted as 𝐮(𝐗, 𝑡) ∶ 0×[0, 𝑡]→R3, which is a vector valued function during the time interval [0, 𝑡]. The absolute temperature field 𝑇(𝐗, 𝑡) ∶ 0× [0, 𝑡]→R+is a smooth scalar-value function. As shown in Fig. 4, the body motion denoted by 𝝋(𝐗, 𝑡) ∶ 0× [0, 𝑡]→R3maps the material point 𝐗in the reference configuration onto its corresponding point 𝐱in the current configuration. The mechanical boundary conditions are imposed on the boundary of deformable body in the reference configuration 𝜕0, which is divided by disjoint sets 𝜕0,𝑢 ⊂ 𝜕0and 𝜕0, 𝑡⊂ 𝜕0with 𝜕0,𝑢 ∩𝜕0, 𝑡= ∅. The Dirichlet boundary condition 𝐮= 𝐮is applied on 𝜕0,𝑢, while Neumann boundary condition characterized by prescribed tractions 𝝈⋅𝐧= 𝐭(𝐗, 𝑡)is applied on 𝜕0, 𝑡. Similarly, the temperature boundary condition 𝑇=𝑇0is prescribed over 𝜕0,𝑇 ⊂ 𝜕0, while the heat flux boundary condition 𝑞=𝑞 is prescribed over 𝜕0, 𝑞 ⊂ 𝜕0, such that 𝜕0,𝑇 ∩𝜕0, 𝑞 = ∅. Besides, the deformable body can also be divided into parts 𝑖 0(𝑖 > 0) such that 0= ∪𝑖𝑖 0, which could be assigned with different mechanical and thermal constitutive behaviors. The deformation gradient, which is defined as the gradient of deformation map with regard to the reference configuration, can be expressed as 𝐅𝑢∶= 𝜕𝐗𝝋(𝐗, 𝑡)=∇𝐗𝝋(𝐗, 𝑡),(1) where ∇𝐗[∙] denotes the Lagrangian gradient operator. This operator represents the map of line element d𝐗in the reference configuration onto the corresponding line element d𝐱=𝐅𝑢d𝐗in the current configuration. The Jacobian of the deformation is defined as the determinant of deformation gradient 𝐽𝑢∶= det[𝐅𝑢]>0, where det[∙] represents the determinant operator. The covariant tangent vectors 𝐆𝑖(𝝃)and 𝐠𝑖(𝝃)are defined as the partial derivatives of position vectors with respect to the convective coordinates 𝜉𝑖in the reference and current configurations, respectively 𝐆𝑖(𝝃) ∶= 𝜕𝐗(𝝃) 𝜕𝜉𝑖,𝐠𝑖(𝝃) ∶= 𝜕𝐱(𝝃) 𝜕𝜉𝑖, 𝑖 = 1,2,3.(2) Journal of the Mechanics and Physics of Solids 169 (2022) 105056 5 Z. Liu et al. The contravariant vectors are defined by 𝐆𝑖⋅𝐆𝑗=𝛿𝑗 𝑖and 𝐠𝑖⋅𝐠𝑗=𝛿𝑗 𝑖, and the metric tensors are defined as 𝐆=𝐺𝑖𝑗 𝐆𝑖⊗𝐆𝑗=𝐺𝑖𝑗 𝐆𝑖⊗𝐆𝑗, 𝐠=𝑔𝑖𝑗 𝐠𝑖⊗𝐠𝑗=𝑔𝑖𝑗 𝐠𝑖⊗𝐠𝑗. The left and right Cauchy–Green deformation tensors, which are derived from displacement field vectors, are given by 𝐂𝑢∶= [𝐅𝑢]T𝐠[𝐅𝑢],𝐛𝑢∶= [𝐅𝑢]𝐆−1[𝐅𝑢]T.(3) The Green–Lagrange strain tensor 𝐄𝑢is defined as 𝐄𝑢∶= 1 2[𝐂𝑢−𝐆].(4) In order to alleviate the locking effects according to Bischoff and Ramm (1997), an incompatible strain tensor  𝐄is incorporated into the displacement derived quantities, and the total Green–Lagrange strain by additive decomposition takes the form of 𝐄∶= 𝐄𝑢+ 𝐄.(5) This decomposition constitutes the fundamentals of assumed strain method, and thus the enhanced right Cauchy–Green tensor 𝐂is given by 𝐂∶= 𝐂𝑢+ 𝐂= 2(𝐄𝑢+ 𝐄) + 𝐆.(6) Then the enhanced deformation gradient 𝐅can be derived by the polar decomposition theorem. The displacement derived deformation gradient can be decomposed as 𝐅𝑢∶= 𝐑𝐔𝑢, where 𝐑and 𝐔𝑢represents the rotation and right-stretch tensors, respectively. Based on Eq. (6), the modified right-stretch tensor 𝐔that accounts for the incompatible strain is calculated as 𝐔∶= √𝐂.(7) Since the calculation of rotation tensor is straightforward, the modified deformation gradient 𝐅can be computed as 𝐅∶= 𝐑𝐔.(8) The modified corresponding Jacobian is given by 𝐽=det[𝐅]. Consider 0⊂0from the continuum body 0in the reference configuration with the delimiting boundary 𝜕0, and its spatial counterpart 𝑡⊂𝑡, with the boundary 𝜕𝑡in the current configuration. It is postulated in Cauchy stress theorem that there is a linear dependency between the normal vector 𝐧of 𝜕𝑡and traction 𝐭through Cauchy stress tensor 𝝈, which is given by 𝐭=𝝈⋅𝐧.(9) In line with this theorem, the first Piola–Kirchhoff stress tensor 𝐏, which is also known as nominal stress tensor, is defined as 𝐏=𝐽𝝈⋅𝐅−1.(10) The first Piola–Kirchhoff traction vector  𝐓is defined as  𝐓=𝐏⋅𝐍, which satisfies the force equality  𝐓d𝑆=𝐭d𝑠, and 𝐍denotes the normal direction in the reference configuration. Accordingly, the second Piola–Kirchhoff stress tensor 𝐒, which is a symmetrical stress tensor, is given by 𝐒=𝐅−1 ⋅𝐏=𝐽𝐅−1 ⋅𝝈⋅𝐅−1.(11) Analogously, Stokes heat flux theorem in the current configuration reads 𝑞𝑛=𝐪⋅𝐧,(12) where 𝑞𝑛and 𝐪represent the scalar and vector heat flux, respectively, and 𝐧denotes the normal direction in the current configuration. The material heat flux 𝐐can be determined from the equality condition 𝐪⋅𝐧d𝑠=𝐐⋅𝐍d𝑆, and thus its definition can be expressed as 𝐐=𝐽𝐅−1 ⋅𝐪.(13) Recalling the theoretical developments outlined in Weber and Anand (1990), Bargmann and Steinmann (2006), the Helmholtz free energy function for the thermo-elastic Kirchhoff–Saint-Venant constitutive material law is given by 𝛹(𝐄, 𝑇 ) = 1 2𝜆(tr[𝐄])2+𝜇tr[𝐄2]−3𝜅𝛼tr[𝐄](𝑇−𝑇0)+𝜌0𝑐𝑝[(𝑇−𝑇0)−𝑇log 𝑇 𝑇0],(14) where 𝜅is the bulk modulus, 𝜆and 𝜇are the Lame constants, 𝜌0is the mass density, 𝛼represents the thermal expansion coefficient and 𝑇0stands for the reference temperature. The second Piola–Kirchhoff stress tensor can be obtained from the constitutive relationship, which reads 𝐒∶= 𝜕𝐄𝛹=𝜆(tr[𝐄])𝟏+ 2𝜇𝐄− 3𝜅𝛼 (𝑇−𝑇0)𝟏,(15) where 𝟏is the second-order identity tensor. The constitutive operators in the curvelinear setting take the form C=𝜕2 𝐄𝐄𝛹=[𝜆𝐺𝑖𝑗 𝐺𝑘𝑙 +𝜇(𝐺𝑖𝑘𝐺𝑗𝑙 +𝐺𝑖𝑙𝐺𝑗𝑘)]𝐆𝑖⊗𝐆𝑗⊗𝐆𝑘⊗𝐆𝑙,(16a) Journal of the Mechanics and Physics of Solids 169 (2022) 105056 6 Z. Liu et al. Fig. 5. The traction versus separation curves of the polymeric interface: (a) sliding and shearing modes, and (b) the opening mode. 𝐙=𝜕2 𝑇𝐄𝛹= −3𝜅𝛼𝐺𝑖𝑗 𝐆𝑖⊗𝐆𝑗.(16b) where Cis the fourth-order material stiffness tensor, and 𝐙is the second-order tensor associated with the coupling term of the Helmholtz function. Assuming isotropic conductivity 𝐤=𝑘𝑔𝑖𝑗 𝐠𝑖⊗𝐠𝑖(Kuhl et al.,2004), the material heat flux vector in the curvilinear setting can be expressed as 𝐐= −𝐽𝐅−1𝐤𝐅−T∇𝐗𝑇= −𝐽𝑘 (𝐆𝑖⊗𝐠𝑖)(𝑔𝑘𝑙𝐠𝑘⊗𝐠𝑙)(𝐠𝑗⊗𝐆𝑗)∇𝐗𝑇= −𝐽𝑘𝐂−1∇𝐗𝑇 . (17) where 𝐂−1 stands for the inverse of right Cauchy–Green tensor. The formulation can be further simplified as 𝐐= −𝑘0𝐆∇𝐗𝑇under isotropic assumption, where 𝑘0represents the thermal conductivity in the reference configuration. 2.2. Thermo-visco-elastic cohesive interface model of polymeric layers Assuming that the polymeric layers can be treated as zero-thickness imperfect interfaces, the displacement and temperature fields inside the encapsulant layers are allowed to be discontinuous through the thickness direction. The cohesive zone model assumes the existence of free energy density per unit undeformed area, and thermo-mechanical coupling can be achieved by making conductivity properties dependent on the cohesive damage and interface tractions dependent on the temperature field. Let define the displacement gaps along the interfaces in three-dimensional setting as: ∆= (𝛥𝑛, 𝛥𝑡, 𝛥𝑠), where 𝛥𝑛, 𝛥𝑠, 𝛥𝑡stand for normal, tangential and shear displacement gaps, respectively, and temperature gap 𝛥𝑇along the interface, which play the role of internal variables for the description of the debonding process along polymeric interfaces. Hence, the coupled thermo-mechanical model for deformable laminae body is enriched due to the presence of cohesive traction field and heat flux through the thickness direction of polymeric layer. Under the assumption of continuity along the interface, the out-of-plane tearing and in-plane sliding traction vectors are defined as 𝜏𝐼={𝐾𝐼(𝑡, ⟨𝑇⟩)𝛥𝐼,if 𝛥𝐼∈𝐽𝐼 0,if 𝛥𝐼∉𝐽𝐼 (18) where 𝐽𝐼=(−𝛿𝑐 𝐼,+𝛿𝑐 𝐼),⟨𝑇⟩is the average temperature along the interface and 𝐼=𝑡, 𝑠, while the opening traction component 𝜎is defined as 𝜎=⎧ ⎪ ⎨ ⎪ ⎩ 𝜖𝛥𝑛,if 𝛥𝑛<0 𝐾𝑛(𝑡, ⟨𝑇⟩)𝛥𝑛,if 𝛥𝑛∈𝐽𝑛 0,if 𝛥𝑛∉𝐽𝑛and 𝛥𝑛≥0 (19) where 𝐽𝑛=(0, 𝛿𝑐 𝑛), and 𝜖is the penalty parameter in compression (Corrado and Paggi,2015). This corresponds to a tension cut-off traction separation cohesive law, in which the interface cannot transfer tractions when the critical opening gap 𝛿𝑐 𝑛is reached. The similar brittle behavior is assumed for the tearing and sliding modes, see Fig. 5. Journal of the Mechanics and Physics of Solids 169 (2022) 105056 7 Z. Liu et al. To obtain the structural response of encapsulant layers, the stiffness 𝐾𝑛is related to the actual stiffness of polymeric materials in the normal direction, which is calculated as the ration between the Young’s modulus 𝐸𝑝𝑜𝑙𝑦 and its thickness ℎ𝑝𝑜𝑙𝑦, i.e., 𝐾𝑛= 𝐸𝑝𝑜𝑙𝑦∕ℎ𝑝𝑜𝑙𝑦. Analogously, the tearing and sliding stiffness can be expressed as: 𝐾𝑠=𝐾𝑡=𝐸𝑝𝑜𝑙𝑦∕[2ℎ𝑝𝑜𝑙𝑦(1 + 𝜈𝑝𝑜𝑙𝑦)]. Since polymeric materials have thermo-visco-elastic constitutive behaviors, the Young’s modulus 𝐸𝑝𝑜𝑙𝑦 depends both on the average temperature ⟨𝑇⟩and time history 𝑡. To synthetically characterize these dependencies, a fractional calculus approach proposed in Di Paola et al. (2013) is adopted here instead of the use of Prony series representation, which has been proved to be very effective for parameters identification (Paggi and Sapora,2015;Lenarda and Paggi,2022). Accordingly, the modulus 𝐸𝑝𝑜𝑙𝑦 is given by 𝐸𝑝𝑜𝑙𝑦(𝑡, 𝑇 ) = 𝑎(𝑇)ℎ(𝑡, 𝑇 )−𝛼(𝑇) 𝛤(1 − 𝛼(𝑇)) (20) where 𝑎and 𝛼are two temperature dependent functions, such that 0< 𝑎, 𝛼 < 1, and 𝛤(𝑡)is the Euler gamma function 𝛤(𝑡) = ∫∞ 0 𝑒−𝑥𝑥𝑡−1 d𝑥(21) The function ℎ(𝑡, 𝑇 )is a time history and temperature dependent function, which is used to model the rheologically complex polymeric materials when the time–temperature superposition principle is not applicable. This can be ascribed to the change of polymer microstructure driven by temperature above a threshold. Hence, this history function ℎ(𝑡, 𝑇 )is equal to the difference between the current time 𝑡and 𝑡0corresponding to the microstructure modification. With regard to the heat conduction, it is assumed that the heat flux across the interface is oriented in the direction orthogonal to the thin polymeric surface. Thus, 𝑞1=𝑞2= 0 and 𝑞=𝑞3is expressed as 𝑞=⎧ ⎪ ⎨ ⎪ ⎩ −ℎ0(1 − 𝛥𝑛 𝛿𝑐 𝑛)⟨𝑇⟩if 𝛥𝑛∈𝐽𝑛 0if 𝛥𝑛∉𝐽𝑛 (22) where ℎ0is the thermal conductivity of interface without crack opening, i.e., 𝛥𝑛= 0. Note that the heat conductivity is assumed to be a decreasing function of gap opening so that the partial heat transfer in case of damaged interface can be taken into account properly (Sapora and Paggi,2014). 2.3. Moisture diffusion along the polymeric interfaces Durability tests of photovoltaic laminates, including damp heat test, humidity freeze test, and thermal cycling test, are characterized by time-dependent temperature and moisture conditions in accordance with prescribed ramps inside a climate chamber. Moisture diffusion mainly takes place along the encapsulant polymeric layers between the solar cell and tempered glass or backsheet. Generally, the aim of the numerical method is to predict the moisture content 𝑐(𝑥1, 𝑥2, 𝑥3, 𝑡)inside the polymeric layer for each material point and time. The boundary value problem for moisture diffusion, in which an imposed moisture content 𝑐∗is applied to the boundary, can be described as follows ⎧ ⎪ ⎨ ⎪ ⎩ 𝜕𝑐 𝜕𝑡 (𝑥1, 𝑥2, 𝑥3, 𝑡)−𝐷∇2𝑐(𝑥1, 𝑥2, 𝑥3, 𝑡)= 0 in 𝑝𝑜𝑙𝑦 ×[0, 𝑡𝑓] 𝑐(𝑥1, 𝑥2, 𝑥3,0)= 0 in 𝑝𝑜𝑙𝑦 𝑐(𝑥1, 𝑥2, 𝑥3, 𝑡)=𝑐∗in 𝜕𝑝𝑜𝑙𝑦 ×(0, 𝑡𝑓](23) where 𝐷is the moisture diffusion coefficient. It should be pointed out that the moisture diffusion is characterized by a different time scale from that of the thermo-mechanical problem. The characteristic velocity of moisture diffusion is related to the diffusion coefficient 𝐷, while that of temperature diffusion is governed by the ratio 𝑘0∕𝜌𝑐𝑝. Given the characteristic values for polymeric layers, the ratio between these two physical phenomena is [𝑘0∕(𝜌𝑐𝑝)]∕𝐷≈ 106. Hence, the moisture diffusion is dependent on the coupled thermo-mechanical problem and not viceversa since heat transfer is about six order faster than moisture diffusion. Based on the experimental evidence reported in Kempe (2006), the moisture diffusion coefficient should be considered as temperature and gap dependent based on an Arrhenius type function 𝐷=⎧ ⎪ ⎨ ⎪ ⎩ 𝐴exp (−𝐸𝑎 𝑇 𝑅 ),if 𝛥𝑛≤𝛿𝑐 𝑛 𝐴exp (−𝐸𝑎 𝑇 𝑅 )𝛥𝑛 𝛿𝑐 𝑛 ,if 𝛥𝑛> 𝛿𝑐 𝑛 (24) For the purpose of accounting for the debonding effect of the encapsulant layer on the moisture diffusion, 𝐷is assumed to be a linear increasing function of the interfacial gap 𝛥𝑛when it overcomes the critical value 𝛿𝑐 𝑛. Given the different time scales of moisture diffusion and heat transfer, a staggered scheme is adopted for the solution of this computational framework, where the diffusion coefficient 𝐷is determined by the average temperature and gap displacement at the previous time increment computed from the coupled thermo-mechanical problem. Journal of the Mechanics and Physics of Solids 169 (2022) 105056 8 Z. Liu et al. 3. Variational form and finite element interpolation In this section, the derivations for the variational form and finite element formulation of the coupled thermo-mechanical problem with cohesive interface and 3D moisture diffusion are outlined. Section 3.1 introduces the multi-field Hu–Washizu variational principle and discretization for the coupled thermo-mechanical problem, which is particularized for solid shell formulation incorporating the EAS and ANS method to alleviate the locking effects. Subsequently, Section 3.2 presents the consistent linearization of the governing equations for the thermo-mechanical cohesive interface. Finally, the finite element implementation of the 3D moisture diffusion along the encapsulant polymeric layers is detailed in Section 3.3. 3.1. Implementation of the coupled thermo-mechanical problem and solid shell formulation In the following, the variational basis and finite element interpolation of the initial boundary value problem with the coupled thermo-mechanical solid shell formulation are presented. The laminae in photovoltaic modules are usually made of thin-walled structures, including tempered glass, silicon solar cell, and backsheet, and to accurately simulate their mechanical behaviors, a solid shell formulation incorporating the EAS and ANS methods for the alleviation of locking effects is derived from the mixed Hu–Washizu variational principle. The weak form of energy balance equation in absence of heat sources and dissipative mechanism in the reference configuration is given by 𝑇(𝐮, 𝐄, 𝑇 , 𝛿𝑇 ) = ∫0 𝜌0𝑐𝑝 𝑇 𝛿𝑇 d𝛺−∫0 𝑇𝐙∶ 𝐄𝛿𝑇 d𝛺+∫0 ∇⋅𝐐𝛿𝑇 d𝛺= 0, (25) where 𝛿𝑇 is the virtual temperature field. Invoking the Gauss’s theorem, the third term relevant to the divergence of heat flux can be reformulated as ∫0 ∇⋅𝐐𝛿𝑇 d𝛺=∫𝜕0, 𝑞 𝑄𝑁𝛿𝑇 d𝜕𝛺 −∫0 𝐐∇𝐗𝛿𝑇 d𝛺, (26) where 𝑄𝑁=𝐐⋅𝐍denotes the Neumann boundary condition on 𝜕0, 𝑞. Subsequently, Eq. (25) can be rewritten as 𝑇(𝐮, 𝐄, 𝑇 , 𝛿𝑇 ) = ∫0 𝜌0𝑐𝑝 𝑇 𝛿𝑇 d𝛺−∫0 𝑇𝐙∶ 𝐄𝛿𝑇 d𝛺+∫𝜕0, 𝑞 𝑄𝑁𝛿𝑇 d𝜕𝛺 −∫0 𝐐∇𝐗𝛿𝑇 d𝛺= 0. (27) By inserting the Duhamel’s law, i.e., simplified Eq. (17) under isotropic assumption, the previous expression can be formulated as 𝑇=∫0 𝜌0𝑐𝑝 𝑇 𝛿𝑇 d𝛺−∫0 𝑇𝐙∶ 𝐄𝛿𝑇 d𝛺+∫𝜕0, 𝑞 𝑄𝑁𝛿𝑇 d𝜕𝛺 +∫0[∇𝐗𝛿𝑇 ]T𝑘0𝐆∇𝐗𝑇d𝛺= 0,(28) with the internal residual 𝑇 int and external counterpart 𝑇 ext being identified as 𝑇 int =∫0 𝜌0𝑐𝑝 𝑇 𝛿𝑇 d𝛺−∫0 𝑇𝐙∶ 𝐄𝛿𝑇 d𝛺+∫0[∇𝐗𝛿𝑇 ]T𝑘0𝐆∇𝐗𝑇d𝛺= 0,(29) 𝑇 ext =∫𝜕0, 𝑞 𝑄𝑁𝛿𝑇 d𝜕𝛺. (30) Analogously, the weak form of the linear momentum balance governing can be derived in this procedure. However, low-order shell elements usually suffer from various locking pathologies, such as volumetric locking, Poisson thickness locking, transverse shear locking and so on. To alleviate such deficiencies, the mixed formulation proposed in Gruttmann and Wagner (2006), Klinkel et al. (2006), Simo and Armero (1992) is adopted here, which incorporates the ANS and EAS methods through the enhancement of the displacement derived strain field by several collocation points and a set of incompatible strains. In this three-field Hu–Washizu variational principle, the displacement 𝐮, the incompatible strain  𝐄, and stress 𝐒are the independent unknown variables for this numerical strategy. Given the use of EAS method, the strain field can be decomposed into a displacement derived compatible part 𝐄𝑢and an enhanced incompatible part  𝐄(Bischoff and Ramm,1997;Hauptmann et al.,2000), i.e., 𝐄=𝐄𝑢+ 𝐄. Accounting for the orthogonality condition between the stress spaces and enhanced strain fields, the stress field can be removed from the subsequent derivations. Hence, the weak form of the linear momentum balance equation can be expressed as 𝑢(𝐮, 𝐄, 𝑇 , 𝛿𝐮) = ∫0 𝐒∶𝛿𝐄𝑢d𝛺−∫0 𝜌0 𝜸𝛿𝐮d𝛺−∫𝜕0, 𝑡  𝐭𝛿𝐮d𝜕𝛺 = 0,(31)  𝐸(𝐮, 𝐄, 𝑇 , 𝛿  𝐄) = ∫0 𝐒∶𝛿 𝐄d𝛺= 𝐸 int = 0,(32) where  𝐸 int represents the internal contribution from the enhanced strain field to the functional. Here the internal residual 𝑢 int and external residual 𝑢 ext corresponding to the kinematic compatible part can be defined as 𝑢 int =∫0 𝐒∶𝛿𝐄𝑢d𝛺, (33) Journal of the Mechanics and Physics of Solids 169 (2022) 105056 9 Z. Liu et al. Fig. 6. The position of collocation points in the parametric space for the ANS method. 𝑢 ext =∫0 𝜌0 𝜸⋅𝛿𝐮d𝛺−∫𝜕0, 𝑡  𝐭⋅𝛿𝐮d𝜕𝛺. (34) The initial configuration 0is discretized into 𝑛𝑒non-overlapping finite elements, i.e., 0≈⋃𝑛𝑒 𝑒=1 (𝑒) 0. For the thermo-mechanical modeling framework stated above, it is particularized for the so-called solid shell parametrization, which assumes the approximation of vector at any material point in terms of the counterparts on the corresponding top and bottom surfaces of the shell element. In line with this parametrization, the position vectors and temperature of any material point in the reference and current configurations can be computed as 𝐗(𝝃) = 1 2(1 + 𝜉3)𝐗𝑡(𝜉1, 𝜉2) + 1 2(1 − 𝜉3)𝐗𝑏(𝜉1, 𝜉2),(35a) 𝐱(𝝃) = 1 2(1 + 𝜉3)𝐱𝑡(𝜉1, 𝜉2) + 1 2(1 − 𝜉3)𝐱𝑏(𝜉1, 𝜉2),(35b) where the subscripts 𝑡and 𝑏stand for the top and bottom surfaces, respectively. The parametric space is given by: ∶= {𝝃= (𝜉1, 𝜉2, 𝜉3) ∈ R3|− 1 ≤𝜉𝑖≤+1; 𝑖= 1,2,3}, where (𝜉1, 𝜉2) represent the in-plane directions and 𝜉3identifies the thickness direction. Based on the isoparametric concept, the approximation of position vectors 𝐗and 𝐱are interpolated through standard trilinear shape functions 𝐍𝐼as follows 𝐗≈ 𝑛𝑛 ∑ 𝐼=1 𝐍𝐼(𝝃)𝐗𝐼=𝐍 𝐗,𝐱≈ 𝑛𝑛 ∑ 𝐼=1 𝐍𝐼(𝝃)𝐱𝐼=𝐍 𝐱,(36) where 𝐗𝐼and 𝐱𝐼are the nodal position vectors in the reference and current configurations, respectively, with number of nodes 𝑛𝑛= 8, and  𝐗and  𝐱are the respective nodal position vectors in the element level. The shape functions 𝐍𝐼is defined as 𝐍𝐼=𝐝𝐢𝐚𝐠 [𝑁𝐼, 𝑁𝐼, 𝑁𝐼](37) where 𝑁𝐼=1 8(1 + 𝜉1 𝐼𝜉1)(1 + 𝜉2 𝐼𝜉2)(1 + 𝜉3 𝐼𝜉3),𝐼= 1,2,…,8. The interpolation of displacement and temperature fields (𝐮, 𝑇 ), the variations (𝛿𝐮, 𝛿𝑇 ), and the increments (𝛥𝐮, 𝛥𝑇 ) are given by 𝐮≈ 𝑛𝑛 ∑ 𝐼=1 𝐍𝐼(𝝃)𝐝𝐼=𝐍𝐝, 𝛿𝐮≈ 𝑛𝑛 ∑ 𝐼=1 𝐍𝐼(𝝃)𝛿𝐝𝐼=𝐍𝛿𝐝, 𝛥𝐮≈ 𝑛𝑛 ∑ 𝐼=1 𝐍𝐼(𝝃)𝛥𝐝𝐼=𝐍𝛥𝐝,(38) 𝑇≈ 𝑛𝑛 ∑ 𝐼=1 𝑁𝐼(𝝃)𝑇𝐼= 𝐍 𝑻, 𝛿𝑇 ≈ 𝑛𝑛 ∑ 𝐼=1 𝑁𝐼(𝝃)𝛿𝑇𝐼= 𝐍𝛿 𝑻, 𝛥𝑇 ≈ 𝑛𝑛 ∑ 𝐼=1 𝑁𝐼(𝝃)𝛥𝑇𝐼= 𝐍𝛥 𝑻, (39) where 𝐝𝐼and 𝑇𝐼represent the nodal displacement and temperature, respectively, and 𝐝and  𝑻are the corresponding vectors at the element level. The interpolation of the spatial temperature gradient in the reference configuration ∇𝐗𝑇, its variation ∇𝐗𝛿𝑇 and its increment ∇𝐗𝛥𝑇 reads ∇𝐗𝑇=𝐆−T∇𝝃𝑇≈𝐆−T∇𝝃 𝐍 𝑻,∇𝐗𝛿𝑇 ≈𝐆−T∇𝝃 𝐍𝛿 𝑻,∇𝐗𝛥𝑇 ≈𝐆−T∇𝝃 𝐍𝛥 𝑻,(40) where ∇𝝃represents the gradient with respect to the natural coordinates in the curvilinear setting. The Green–Lagrange strain vector is expressed as 𝐄𝑢=[E11,2E12,2E13,E22,2E23,E33]T. To overcome the curvature thickness locking, the ANS method proposed in Betsch and Stein (1995) is adopted to modify the strain component E33, which requires four collocation points defined in convective coordinates 𝝃C𝑖as 𝝃C1= (−1,−1,0),𝝃C2= (1,−1,0),𝝃C3= (1,1,0), and 𝝃C4= (−1,1,0), see Fig. 6. Besides, to prevent transverse shear locking, the ANS method proposed in Dvorkin and Bathe (1984) is also adopted in this work. The corresponding collocation points are 𝝃A1= (0,−1,0),𝝃A2= (0,1,0),𝝃B1= (−1,0,0), and 𝝃B2= (1,0,0), see Fig. 6. Journal of the Mechanics and Physics of Solids 169 (2022) 105056 16 Z. Liu et al. Fig. 13. (a) Relaxation modulus of EVA vs. time curves at different temperatures, (b) temperature dependent fractional calculus coefficient 𝑎(𝑇), (c) temperature dependent exponent 𝛼(𝑇). Source: Adapted from Paggi and Sapora (2015). Journal of the Mechanics and Physics of Solids 169 (2022) 105056 17 Z. Liu et al. Fig. 14. The temperature contour plot of the cross section of the double glass laminate at time point A in the first cycle during humidity freeze simulation. Fig. 15. The temperature contour plot inside EVA at six different time points in the first cycle during humidity freeze simulation. Journal of the Mechanics and Physics of Solids 169 (2022) 105056 18 Z. Liu et al. Fig. 16. The temperature versus distance from edge curves at six different time points in the first cycle during humidity freeze simulation. where 𝐍𝑐= [𝑁1, 𝑁2, 𝑁3, 𝑁4, 𝑁5, 𝑁6, 𝑁7, 𝑁8]is the shape function matrix, and  𝐜is the moisture concentration nodal vector. Besides, the interpolation of the spatial concentration gradient ∇𝐗𝑐and its variation ∇𝐗𝛿𝑐 can be expressed as ∇𝐗𝑐≈𝐆−T∇𝝃𝐍𝑐 𝐜=𝐁𝑐 𝐜,∇𝐗𝛿𝑐 ≈𝐆−T∇𝝃𝐍𝑐𝛿 𝐜=𝐁𝑐𝛿 𝐜,(77) Introducing Eqs. (76) and (77) into Eq. (75), the variational form of moisture diffusion is given by 𝛿 𝐜T∫0(𝐷𝐁T 𝑐∇𝐗𝑐+𝐍T 𝑐𝑐)d𝛺=𝛿 𝐜T𝐟𝑐 int (78) where 𝐟𝑐 int is the internal residual vector, which is defined as 𝐟𝑐 int =∫0(𝐷𝐁T 𝑐∇𝐗𝑐+𝐍T 𝑐𝑐)d𝛺. (79) The linearization for moisture diffusion at the (k+ 1)𝑡ℎ iteration in the Newton–Raphson iterative solution scheme reads 𝐊𝑐,𝑘𝛥 𝐜𝑘+1 = −𝐟𝑐,𝑘 int (80a)  𝐜𝑘+1 = 𝐜𝑘+𝛥 𝐜𝑘+1 (80b) where 𝐊𝑐is the element stiffness for moisture diffusion, which is given by 𝐊𝑐=∫0(𝐷𝐁T 𝑐𝐁𝑐+1 𝛥𝑡 𝐍T 𝑐𝐍𝑐)d𝛺(81) 4. Staggered solution scheme for the multi-field framework In this section, the staggered computational procedure for the solution of the thermo-mechanical problem with solid shell elements and interface finite elements, along with temperature-dependent 3D moisture diffusion, is outlined. With regard to the thermo-mechanical governing equations incorporating the solid shell formulation, let consider the time increment [𝑡𝑛, 𝑡(𝑘) 𝑛+1], where 𝑡𝑛and 𝑡(𝑘) 𝑛+1 stand for the previous converged increment and prospective current increment at iteration k, respectively. Given the data {𝐝𝑛,𝝇𝑛, 𝐓𝑛}at the previous converged increment, the nonlinear incremental solution requires the Newton iterations at the intermediate state {𝐝(𝑘) 𝑛+1,𝝇(𝑘) 𝑛+1, 𝐓(𝑘) 𝑛+1}. Note that {𝐝𝑛,𝝇𝑛, 𝐓𝑛}and {𝐝(𝑘) 𝑛+1,𝝇(𝑘) 𝑛+1, 𝐓(𝑘) 𝑛+1}denote the nodal displacement vector, the enhancing vector, and the nodal temperature vector at the previous converged increment and prospective current increment at (k)𝑡ℎ iteration, respectively. Based on the static condensation described in Section 3.1, the nodal displacement and temperature vectors are defined as unknowns in the element level, and the increment of enhancing vector in the (k)𝑡ℎ iteration of next time increment, 𝛥𝝇(𝑘) 𝑛+1, should be determined. According to the procedure proposed in Reinoso and Blázquez (2016), 𝛥𝜻(𝑘) 𝑛+1 is given by 𝛥𝝇(𝑘) 𝑛+1 = − [𝐤𝜍𝜍,𝑛]−1 [𝐑𝜍 int,𝑛 +𝐤𝜍𝑑,𝑛𝛥𝐝(𝑘) 𝑛+1 +𝐤𝜍𝑇 ,𝑛𝛥 𝐓(𝑘) 𝑛+1](82) It should be pointed out that the increments 𝛥𝐝(𝑘) 𝑛+1 and 𝛥 𝐓(𝑘) 𝑛+1 are provided by the solver, while the element matrices at the previous increment [𝐤𝜍𝜍,𝑛]−1,𝐑𝜍 int,𝑛,𝐤𝜍𝑑,𝑛,𝐤𝜍𝑇 ,𝑛, and 𝝇𝑛are all stored as internal variables. The numerical algorithm for the finite element Journal of the Mechanics and Physics of Solids 169 (2022) 105056 19 Z. Liu et al. Fig. 17. The contour plot of moisture diffusion along the EVA layer during the humidity freeze simulation. implementation of the thermo-mechanical problem with solid shell and interface element formulation is detailed in Algorithm 1. The Newton–Raphson iteration is performed until the machine precision is achieved, i.e., the tolerance of residual vector is up to 10−15. It is worth noting that the history variables ℎat the integration points need to be determined due to the time dependency of visco-elastic constitutive equation, Eq. (20). To model relaxation, this history variable is set to zero in case of any temperature change, which is stored as a state variable in the interface element implementation, while updated by the current time increment if the temperature value remains constant compared with that of previous increment, see Algorithm 1. Hence, the drawback of temperature–time superposition principle for thermo-visco-elastic modeling of polymeric materials is overcome. Once the temperature and displacement field in the coupled thermo-mechanical problem are solved in a given time increment, these nodal values are transferred to the moisture diffusion problem, which is solved by the Euler backward time integration scheme using the same temporal interval with that of the thermo-mechanical analysis. As stated in Section 2.3, given the time scale difference between the moisture diffusion and thermo-mechanical phenomena, the coefficient in the current increment of moisture diffusion Journal of the Mechanics and Physics of Solids 169 (2022) 105056 20 Z. Liu et al. analysis is determined from the nodal temperature and gap values at the previous time increment. The detailed computational procedure of moisture diffusion is outlined in Algorithm 2. Algorithm 1: Numerical implementation procedure of the coupled thermo-mechanical formulation if Solid shell element then Data: d𝑛, 𝐓𝑛,𝛥d(𝑘) 𝑛+1,𝛥 𝐓(𝑘) 𝑛+1 Result: d𝑛+1, 𝐓𝑛+1 Initialization of 𝝇𝑛,𝐑𝜍 int,𝑛,𝐤𝜍𝑑,𝑛,[𝐤𝜍𝜍,𝑛]−1; while ||R∗ 𝑑||>tolerance do Compute 𝛥𝝇(𝑘) 𝑛+1 = − [𝐤𝜍𝜍,𝑛]−1 [𝐑𝜍 int,𝑛 +𝐤𝜍𝑑,𝑛𝛥𝐝(𝑘) 𝑛+1 +𝐤𝜍𝑇 ,𝑛𝛥 𝐓(𝑘) 𝑛+1]; Update the enhancing vector 𝝇(𝑘) 𝑛+1 =𝝇𝑛+𝛥𝝇(𝑘) 𝑛+1; for n←1to 8 integration points do Compute the curvilinear basis 𝐆(𝑘) 𝑛+1 and 𝐠(𝑘) 𝑛+1; Compute the B matrices 𝐁(𝑘) 𝑛+1 and 𝐁(𝑘) 𝑇 ,𝑛+1; Modify the 𝐁(𝑘) 𝑛+1 matrix according to the ANS method; Compute C(𝑘) 𝑛+1,𝐒(𝑘) 𝑛+1, and 𝐙(𝑘) 𝑛+1; Compute the EAS operator 𝐌(𝑘) 𝑛+1; end Compute the element stiffness matrices 𝐤(𝑘) 𝑑𝑑,𝑛+1,𝐤(𝑘) 𝑑𝜍,𝑛+1,𝐤(𝑘) 𝜍𝑑,𝑛+1,𝐤(𝑘) 𝜍𝜍,𝑛+1, and 𝐤(𝑘) 𝑇 𝑇 ,𝑛+1; Compute the internal force vectors 𝐑𝑑(𝑘) int,𝑛+1,𝐑𝜍(𝑘) int,𝑛+1, and 𝐑𝑇(𝑘) int,𝑛+1; Perform the static condensation and final assembly; end end else if Cohesive interface element then Data: d𝑛, 𝐓𝑛,𝛥d(𝑘) 𝑛+1,𝛥 𝐓(𝑘) 𝑛+1 Result: d𝑛+1, 𝐓𝑛+1 while ||R∗ 𝑑||>tolerance do Compute the operators 𝐋𝑐𝑧 and 𝐌𝑐𝑧; for n←1to 4 integration points do Compute the shape function 𝐍(𝑘) 𝑐𝑧,𝑛+1 ; Compute the matrix 𝐁(𝑘) 𝑐𝑧,𝑛+1 =𝐍(𝑘) 𝑐𝑧,𝑛+1𝐋𝑐𝑧; Compute the rotation matrix in the curvilinear system 𝐑(𝑘) 𝑐𝑧,𝑛+1; Compute the local gap 𝐠(𝑘) loc,𝑛+1 =𝐑(𝑘) 𝑐𝑧,𝑛+1𝐁(𝑘) 𝑐𝑧,𝑛+1  𝐝(𝑘) 𝑛+1; Compute the temperature value at the integration point 𝑇(𝑘) 𝑖𝑛𝑡,𝑛+1; if |𝑇(𝑘) 𝑖𝑛𝑡,𝑛+1 −𝑇𝑖𝑛𝑡,𝑛|>tolerance then Compute the history variable ℎ(𝑘) 𝑖𝑛𝑡,𝑛+1 =ℎ𝑖𝑛𝑡,𝑛; else Compute the history variable ℎ(𝑘) 𝑖𝑛𝑡,𝑛+1 =ℎ𝑖𝑛𝑡,𝑛 +𝛥𝑡; end Compute the stiffness matrix 𝐂(𝑘) 𝑐𝑧,𝑛+1 and the traction vector 𝐓(𝑘) 𝑛+1; end Compute the element stiffness matrix 𝐊e(𝑘) int,𝑛+1; Compute the internal residual vector 𝐟(𝑘) int,𝑛+1; end end Journal of the Mechanics and Physics of Solids 169 (2022) 105056 21 Z. Liu et al. Algorithm 2: Numerical implementation procedure of the moisture diffusion along the interface Data:  𝐜𝑛,𝛥 𝐜(𝑘) 𝑛+1 Result:  𝐜𝑛+1 Initialization of 𝐝𝑛and  𝐓𝑛; while ||𝐟𝑐 𝑖𝑛𝑡||>tolerance do for n←1to 8 integration points do Compute the shape function 𝐍(𝑘) 𝑐,𝑛+1; Compute the B matrix 𝐁(𝑘) 𝑐,𝑛+1; Compute the temperature value at the integration point 𝑇𝑖𝑛𝑡,𝑛 with  𝐓𝑛; Compute the normal gap at the integration point 𝛥𝑖𝑛𝑡,𝑛 with 𝐝𝑛; if 𝛥𝑖𝑛𝑡,𝑛 ≤𝛿𝑐 𝑛then Compute the coefficient 𝐷(𝑘) 𝑖𝑛𝑡,𝑛+1 =𝐴exp (−𝐸𝑎 𝑇𝑖𝑛𝑡,𝑛𝑅); else Compute the coefficient 𝐷(𝑘) 𝑖𝑛𝑡,𝑛+1 =𝐴exp (−𝐸𝑎 𝑇𝑖𝑛𝑡,𝑛𝑅)𝛥𝑖𝑛𝑡,𝑛 𝛿𝑐 𝑛; end end Compute the element stiffness matrices 𝐊(𝑘) 𝑐,𝑛+1; Compute the internal force vectors 𝐟𝑐(𝑘) 𝑖𝑛𝑡,𝑛+1; Perform the final assembly; end 5. Numerical applications to photovoltaics In this section, three international standard tests for photovoltaics, including damp heat test, humidity freeze test, and thermal cycling test, are simulated by the proposed computational framework. Note that the damp heat test requires constant temperature and humidity testing environments, which uncouples the moisture diffusion from the thermo-mechanical problem and thus allows the derivation of analytical solution useful for benchmark targets, while the latter require the fully coupled solution scheme due to the spatial variation of temperature. Besides, the role of cracks in silicon on moisture diffusion pattern is also investigated and compared with experimental electroluminescence images. 5.1. Damp heat test As prescribed by the international standard, the damp heat test is conducted in the environmental chamber with the constant temperature of 85 ◦C and relative humidity of 85% RH. This humidity environmental condition corresponds to the initial boundary condition with a moisture concentration 𝑐∗= 0.0056 g/cm3imposed on the free edges of the photovoltaic laminate. Since temperature is kept constant in this test, the moisture diffusion problem can be solved independently from the coupled thermo-mechanical problem with a constant diffusivity. In this work, the specimen is a double-glass PV laminate of span 𝑙= 40 mm separated by EVA. The thickness of glass laminae is 3 mm, and the thickness of EVA is 0.5 mm. As glass is not permeable, the moisture will only diffuse from the free edges of EVA towards the central part in this test, see Fig. 9. The analytical solution was obtained in Crank (1979), which will be treated as a benchmark of the computational scheme proposed above. Assuming the initial laminate is dry, the spatial and temporal solution 𝑐 can be obtained as 𝑐=𝑐∗+4𝑐∗ 𝜋 ∞ ∑ 𝑚=0 1 2𝑚+ 1 sin [(2𝑚+ 1)𝜋𝑥 𝑙]e(−𝐷(2𝑚+1)2𝜋2𝑡∕𝑙2)(83) The mesh size of finite element model is 1 mm, and there are 1600 solid diffusion elements in total. The contour plot of moisture concentration inside the EVA is shown in Fig. 10 and, as clearly seen, water concentration gradually diffuses from the edges into the middle area of EVA during simulation. At 100 h, the most area of EVA is still dry with almost zero concentration, see Fig. 10(a), while after 1000 h, the dry area becomes obviously smaller. To quantitatively validate the computational diffusion model, the predicted moisture concentration with respect to the distance from the edge after 1000 h (the black solid line in Fig. 11) is compared with that obtained from the analytical solution (the black dotted line in Fig. 11), and they agree with each other very well. Besides, this model is also used to simulate the experimental results obtained from Miami as reported in Kempe (2006), and numerical predictions also agree very well with the data from the impermeable double-glass laminate exposed to Miami environment for 1 year, 2 years, and 3 years, which further proves the reliability of the diffusion modeling method proposed in this work. Journal of the Mechanics and Physics of Solids 169 (2022) 105056 22 Z. Liu et al. Fig. 18. The comparison of moisture concentration versus distance from edge curves after 700 h and 1000 h between the damp heat and humidity freeze simulations. Fig. 19. The comparison of time history of moisture concentration at the positions 3 mm, 5 mm, and 8 mm away from the edge between the damp heat and humidity freeze simulations. 5.2. Humidity freeze test In the humidity freeze test of PV laminates, as requested by the international standard, the modules are subjected to the cycling temperature condition from −40 ◦C to 85 ◦C with the constant humidity condition of 85% RH. The temperature ramp can be expressed as 𝑇∗(𝑡) = ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ 𝑡 𝑡∗ 1 𝑇∗ 10≤𝑡<𝑡∗ 1 𝑇∗ 1𝑡∗ 1≤𝑡<𝑡∗ 2 𝑡∗ 3−𝑡 𝑡∗ 3−𝑡∗ 2 𝑇∗ 1𝑡∗ 2≤𝑡<𝑡∗ 3 𝑡−𝑡∗ 3 𝑡∗ 4−𝑡∗ 3 𝑇∗ 2𝑡∗ 3≤𝑡<𝑡∗ 4 𝑇∗ 2𝑡∗ 4≤𝑡<𝑡∗ 5 𝑡∗ 6−𝑡 𝑡∗ 6−𝑡∗ 5 𝑇∗ 2𝑡∗ 5≤𝑡<𝑡∗ 6 (84) where 𝑇∗ 1= 85 ◦C,𝑇∗ 2= −40 ◦C, and 𝑡∗ 1= 1.0h, 𝑡∗ 2= 21.0h, 𝑡∗ 3= 22.0h, 𝑡∗ 4= 22.5h, 𝑡∗ 5= 23.5h, 𝑡∗ 6= 24.0h, as shown in Fig. 12. Compared with the damp heat test, this case is much more complex and almost impossible to obtain the analytical solution due to the spatial and temporal temperature variation inside the EVA, to which a non-constant temperature boundary condition is applied. Particularly, the cohesive properties and diffusivity of EVA need to be updated during simulation, and to be specific, the viscoelastic parameters 𝛼(𝑇)and 𝑎(𝑇)for the calibration of Young’s modulus of EVA are temperature dependent as experimentally Journal of the Mechanics and Physics of Solids 169 (2022) 105056 23 Z. Liu et al. Fig. 20. EL images of the minimodule taken at different cycles during the thermal cycling test. Source: Adapted from Berardone and Paggi (2019). Fig. 21. Temperature profile imposed inside the chamber during the thermal cycling test. Journal of the Mechanics and Physics of Solids 169 (2022) 105056 24 Z. Liu et al. Fig. 22. The sketch of the PV minimodule in the thermal cycling test. Fig. 23. The temperature contour plot of the cross section of the minimodule at time points B and E in the first cycle during the thermal cycling test. Table 1 Mechanical and thermal properties of the PV materials. 𝐸(GPa)𝜌(kg∕m3)𝛼 𝑐𝑝(J/kg K)𝑘0(W/m K) Backsheet 2.8 1000 5.04e-5 300 0.36 Glass 73 2300 8e-6 500 0.8 Silicon 130 2500 2.49e-6 715 148 evaluated in Paggi et al. (2011) and interpreted through the fractional calculus method (Sapora and Paggi,2014), see the plot of modulus 𝐸,𝑎(𝑇), and 𝛼(𝑇)in Fig. 13. Regarding the properties of moisture diffusion, the coefficient 𝐷is determined according to the Arrhenius type equation Eq. (24) reported in Kempe (2006). The thermal and mechanical properties of different PV materials are taken from Eitner et al. (2011), Paggi et al. (2011), which are listed in Table 1. Besides, the critical opening of polymeric interface 𝛿𝑐 𝑛can be estimated from the experimental data of variation of Mode 1 fracture energy with respect to temperature in Novoa et al. (2014), since the area below the traction versus separation curve is the fracture energy. The temperature variation inside the double glass laminate is strongly dependent on the heat conduction properties of different layers, as shown by the temperature contour of the cross-section at time point A in Fig. 14. The contour plot of temperature distribution inside the EVA layer in the first cycle during simulation is shown in Fig. 15. At time point A, the boundary temperature reaches 85 ◦C, and is kept constant for 20 h. As shown in Fig. 15(a)–(c), heat gradually diffuses from the edges into the middle area of EVA during this period. At time point D, the boundary temperature drops to −40 ◦C, but the temperature in the middle area of EVA layer is still higher than 65 ◦C due to the inertial effect of heat conduction, see Fig. 15(d)–(f). Remarkably, this spatial and temporal variation of temperature inside EVA leads to the difference of moisture diffusion properties from that of EVA subjected to constant condition in the damp heat test. To quantitatively describe the spatial variation of temperature inside EVA, the temperature versus distance from edge curves at the six different time points in the first cycle is plotted in Fig. 16. As can be seen, the plateau value, which corresponds to the central area temperature of EVA, gradually goes up from time point A to C with the evolution of thermal conduction. At time C, the temperature at almost all the area of EVA from the edge to the center reaches around 85 ◦C, which means the EVA layer is fully and evenly heated after 20 h. Besides, the plateau at the time period when boundary temperature changes to −40 ◦C, which can be ascribed to the hysteresis of heat conduction, gradually goes down from time point D to F. It is worth noting that the temperature at the central area is still higher than 30 ◦C at the end of one complete cycle due to the relatively short duration of constant low-temperature boundary loading. The contour plot history of moisture diffusion during the humidity freeze simulation is shown in Fig. 17. The moisture gradually diffuses into the central area, which is very similar to the phenomena obtained from previous damp heat simulation. However, the Journal of the Mechanics and Physics of Solids 169 (2022) 105056 25 Z. Liu et al. Fig. 24. The temperature contour plot inside EVA at six different time points in the first cycle during thermal cycling simulation. amount of moisture concentration is totally different in two cases due to the spatial variation of temperature inside EVA in the humidity freeze simulation. As shown in Fig. 18, the difference of moisture concentration spatial variation after 700 h and 1000 h between the damp heat and humidity freeze simulation can be observed. Note that the amount of moisture concentration in all the area of EVA after 700 h and 1000 h in the damp heat case is higher than that in the humidity freeze case, since the moisture diffusion coefficient in the former case is much higher at the constant high-temperature condition (85 ◦C). Furthermore, It can be seen from Fig. 18 that the amount of moisture concentration at all the areas of EVA after 700 h in the damp heat case is even higher than that in the humidity freeze case after 1000 h, which further demonstrates the effect of temperature on moisture diffusion. Besides the difference of spatial variation, the temporal variation between the damp heat and humidity freeze simulation is also shown and compared in this work, see Fig. 19. The time history of moisture concentration at three positions 3 mm, 5 mm, and 8 mm away from the edge is also plotted, and the curves obtained from the damp heat simulations (solid lines) are higher than all the counterparts obtained from the humidity freeze simulations (dotted lines). Thus, conclusions can be drawn that temperature plays an import role on the moisture diffusion problem as the coefficient is determined from the Arrhenius equation, and this coupling can be described with the proposed modeling method, which can be demonstrated by the example of humidity freeze simulation.