A geometrical multi-scale numerical method for coupled hygro-thermo-mechanical problems in photovoltaic laminates

A comprehensive computational framework based on the finite element method for the simulation of coupled hygro-thermo-mechanical problems in photovoltaic laminates is herein proposed. While the thermo-mechanical problem takes place in the three-dimensional space of the laminate, moisture diffusion occurs in a two-dimensional domain represented by the polymeric layers and by the vertical channel cracks in the solar cells. Therefore, a geometrical multi-scale solution strategy is pursued by solving the partial differential equations governing heat transfer and thermo-elasticity in the three-dimensional space, and the partial differential equation for moisture diffusion in the two dimensional domains. By exploiting a staggered scheme, the thermo-mechanical problem is solved first via a fully implicit solution scheme in space and time, with a specific treatment of the polymeric layers as zero-thickness interfaces whose constitutive response is governed by a novel thermo-visco-elastic cohesive zone model based on fractional calculus. Temperature and relative displacements along the domains where moisture diffusion takes place are then projected to the finite element model of diffusion, coupled with the thermo-mechanical problem by the temperature and crack opening dependent diffusion coefficient. The application of the proposed method to photovoltaic modules pinpoints two important physical aspects: (i) moisture diffusion in humidity freeze tests with a temperature dependent diffusivity is a much slower process than in the case of a constant diffusion coefficient; (ii) channel cracks through Silicon solar cells significantly enhance moisture diffusion and electric degradation, as confirmed by experimental tests.


Introduction and motivations
Photovoltaic modules are composite structures obtained by laminating layers of various materials. Some have the role to guarantee protection from the environment by a suitable sealing of the device, others have specific electric features to produce energy. The durability of these devices is a serious concern due to fracture events promoted by the mismatch between the thermo-mechanical properties of the material constituents, often amplified by the severe working conditions they are exposed to. Moreover, their modelling is changelling and it requires a multi-physics framework to gain an accurate picture of their overall working conditions, performance and degradation [1].
Typical photovoltaic (PV) modules are laminates made of a thick glass superstrate, an encapsulating polymer (usually ethylene vinyl acetate, EVA), a layer of interconnected Silicon solar cells separated by few centimeters of EVA in their plane, another layer of EVA, and finally a polymeric protective backsheet. EVA provides protection of cells and interconnections but it is permeable to moisture, which diffuses from the backsheet and percolates along the surface of the solar cells. In turn, moisture induces chemical oxidation of the grid line deposited on the surface of the solar cells, giving rise to electrically inactive areas and power-loss. This phenomenon has been reported in damp heat tests in [2] prescribed by the international qualification standards [4], where PV modules were exposed to a very aggressive environment at constant 85 • C temperature and 85% of air humidity. In particular, it has been shown a progressive increase of dimmer areas in time in the electroluminescence (EL) images starting from the edges of the solar cells towards their center (see Fig. 1a). Correspondingly, the current-voltage of the PV module degrades, with a significant power-loss (see Fig. 1b).
For these reasons, an accurate modelling of the EVA is crucial to determine the lifetime of a PV module, moisture diffusion, chemical reactions induced by moisture, as well as its reduction of cohesive energy promoting delamination of the layers (see Fig. 2).
The EVA polymer displays a strong thermo-visco-elastic constitutive response, as experimentally reported in [5,6], with a variation of the elastic modulus of up to three orders of magnitude depending on temperature. Generalized Maxwell rheological models used so far generally provide exponential type relations for the relaxation modulus and, in order to approximate the experimentally observed power-law trend, a huge number of elements (and thus of model parameters) has to be taken into account. To significantly simplify the Fig. 2 Example of delamination of the glass from the solar cells, adapted from [3] Fig. 3 Multi-physics modelling showing the proposed solution schemes and the interactions between the various fields task of parameters identification, modelling the visco-elastic behaviour via fractional derivatives has been proved to be very effective [8][9][10]. For rheologically complex polymers as EVA, whose microstructure changes with temperature, the fractional calculus formulation in [11] allows the use of only two temperature dependent parameters for its complete description.
Another complexity regards the moisture diffusion properties of EVA, strongly temperature dependent as experimentally reported in [12], which implies coupling between moisture and temperature fields. In return, moisture is degrading the cohesive energy of the EVA encapsulant, promoting decohesion of the backsheet or separation of the glass cover from the solar cells [13,14]. This mathematically corresponds to coupling between the mechanical and the diffusive fields, with a further acceleration of moisture diffusion and degradation as reported in [15], which is a feedback coupling from diffusion to mechanics.
A proper modelling of these coupled nonlinear phenomena, whose conceptual interactions are sketched in Fig. 3, requires a comprehensive computational framework where coupled thermo-elastic and heat conduction problems are accurately solved at the module level in the three-dimensional space. Then, moisture diffusion inside the EVA layers has to be simulated by considering the dependency on the temperature and the thermo-elastic fields via the diffusive constitutive equations. So far, the state-of-the-art simulations on moisture diffusion in [12] consider the EVA layer only and treat diffusion as a one dimensional problem without updating the diffusivity of the material based on the actual temperature of the system. The former approximation of considering moisture diffusion as a pure 1D problem fails when channel cracks in Silicon solar cells are present, since they can also be a source of moisture percolation from the backsheet to the front side of the solar cells. The latter assumption of using a constant diffusivity (uncoupling with the thermal field) is an approximation valid only in the steady state temperature regime. While this is the case of the damp heat test, its validity in the case of a cyclic variation of temperature from −40 up to 85 • C as in the humidity freeze test is highly questionable.
To shed light into the above issues, and provide a comprehensive physical modelling and computational framework for the study of these phenomena, a geometrical multiscale approach is herein proposed by following the seminal work in [16] for biophysical systems. Starting from the evidence that moisture diffusion takes place in a physical domain with a lower dimension with respect to that of the thermo-mechanical and heat conduction problems, two different finite element models are used in parallel. The coupled thermo-mechanical and heat conduction problems are solved in the three-dimensional setting (or in the two-dimensional one in the case of a cross-section of the PV module). As a further simplification, the EVA encapsulant layers are modelled as zero-thickness interfaces, whose thermo-viscoelastic constitutive response is taken into account by a novel thermo-viscoelastic cohesive zone model. As compared to other cohesive zone model formulations available in the literature [17,18], the present formulation is based on fractional calculus and it is able to simulate rheologically complex materials.
The thermo-mechanical problem, which is much faster than moisture diffusion, is solved first via a fully implicit solution scheme in space and time, see Fig. 2. Temperature and relative displacements computed in the Gauss points along the encapsulant interfaces are then projected to the nodes of another finite element model specific for the solution of moisture diffusion. This second model is used to discretize the domain where moisture diffusion takes place. In particular, it is represented by the mid-surfaces of the encapsulant layers and of the channel cracks through Silicon, see Fig. 4.
This article is structured as follows. In Sect. 2, the variational framework for thermo-mechanics and heat conduction for the layers is presented, along with the interface model for the thermo-visco-elastic encapsulant, as well as for moisture diffusion. The weak forms of the partial differential equations are established in Sect. 3 and the finite element

Variational framework
In this section, the variational framework describing moisture diffusion and the thermo-visco-elastic response of a laminate made of linear elastic homogeneous and isotropic layers separated by polymeric thermo-visco-elastic laminae is presented. In these laminates, moisture diffusion takes place in the polymeric layers, which progressively percolates from the free edges and from the permeable backsheet towards the centre of the solar cells. EVA layers used to protect Silicon solar cells are permeable to moisture, which is one of the major concerns for the degradation of the electrical output of the PV module over time. Due to moisture diffusion, the adhesive properties of EVA progressively degrade and the corresponding layers may experience a lack of cohesion leading to separation of the backsheet from the Silicon cells, or between the Silicon cells from the glass superstrate. In order to efficiently simulate cohesive degradation of EVA and delamination, we propose to treat the polymeric layers as zero-thickness internal interfaces with suitable tractionseparation relations accounting for their thermo-visco-elastic properties.
Let the laminate occupy a volume Ω = ∪ n m=1 R(m) ⊂ R 3 in the reference undeformed configuration, where each layer is geometrically identified by h m for m = 1, . . . , n. Moreover, let model a generic polymeric layer of thickness h between laminae R( p) = R − and R( p The position of each material point inside Ω is identified by the coordinate vector x = (x 1 , x 2 , x 3 ) T in a three-dimensional cartesian orthonormal frame {e 1 , e 2 , e 3 }. Let u I (x 1 , x 2 , x 3 , t) (I = 1, 2, 3) and θ(x 1 , x 2 , x 3 , t) = T (x 1 , x 2 , x 3 , t) − T 0 be the displacement field and temperature variation from a reference one, T 0 , inside the material during the time interval 0 ≤ t ≤ t f . The index I is used to denote the component of the displacement field along the corresponding coordinate.

Thermo-mechanical formulation of the layers
The global dynamics of each material layer R(m) obeys the equations of coupled linear isotropic thermo-elasticity (see e.g. [19,20]). The Cauchy thermal stress tensor is defined for each layer m as: where the Einstein summation notation has been adopted.
Here, C m I J K L is the fourth order elastic constitutive tensor, and β m is the coupling thermal factor related to the thermal expansion coefficient α m . The infinitesimal strain tensor is: The balance of linear momentum is given by: where ρ m is the density of the m-th material. Let q I be the heat flux inside each layer R(m), and assume that it is related to the temperature variation θ by the Fourier law: where k m > 0 is the thermal conductivity of the m-th material. Hence, the heat transfer partial differential equation is given by: where c m is the heat capacity of the m-th material.

Thermo-visco-elastic polymeric interfaces
Under the assumption of replacing each polymeric layer by a zero-thickness imperfect interface, we allow u I and the temperature θ to be discontinuous across the interface separating R − from R + . We therefore define the jumps on S( p) as: The average temperature across the interface is given by: Hence, the coupled thermo-elastic model given in the previous section for the continuum layers is enriched by adding the presence of a cohesive traction field and a heat flux normal to S( p). The relations between those fields and [[u I ]], [[θ ]] are provided by the constitutive relations for the interface. By assuming the continuity of the traction vector components t I , the in-plane sliding and out-of-plane tearing components of the traction vector read: where I = 1, 2 and J I = (−δ c I , +δ c I ), while the opening traction component is: where J 3 = (0, δ c 3 ). Equation (9) corresponds to a tension cut-off tractionseparation curve, suggesting that the interface is not able to transfer tractions when the critical opening displacement δ c 3 is overcome. A similar brittle behavior is assumed for sliding and tearing fracture modes (see Fig. 5). In compression, a penalty formulation with penalty parameter is adopted as in [7].
In order to obtain a structural response of the interface equivalent to that of the real EVA layers, the modulus K 3 of the zero-thickness interface is related to the actual stiffness of the layer in the direction n, i.e., it can be evaluated as the ratio between the EVA Young's modulus, E EVA , and its thickness, h EVA , i.e., Since polymers have a thermo-visco-elastic constitutive behavior, the stiffnesses K I have to depend both on the average temperature θ and time t. Instead of using a Prony series representation, a fractional calculus approach [8,9] is herein adopted to synthetically characterize those dependencies. This approach has been proved in [11] to be very effective for parameters identification. Accordingly, E EVA is defined as follows: for functions 0 < a, α < 1 and Γ (t) is the Euler gamma function Function h(t, T ) is a history dependent function of time and temperature used to model thermo-rheologically complex materials where the principle of time-temperature superposition does not apply. This is due to a modification of the internal microstructure of the polymer driven by a temperature change above a threshold. Hence, h(t, T ) is equal to the current time t minus the time t 0 corresponding to such a microstructure modification. At the interface, we remark that cohesive tractions are continuous and opposing to each other, viz.
As far as heat conduction is concerned, we assume that the heat flux across the interfaces is oriented along the direction orthogonal to the surface S( p), which is a reasonable approximation for thin polymeric layers. Hence, q 1 = q 2 = 0 and q = q 3 is given by where κ 0 is the thermal conductivity of the interface without decohesion, i.e., for [[u 3 ]] = 0. Note that the heat flux is assumed to be a decreasing function of the normal gap, in order to account for partial heat transfer in the case of a damaged interface, see also [21]. Continuity of heat flux at the interface is preserved, viz.

Moisture diffusion along polymeric interfaces
Durability tests of PV modules inside a climate chamber are characterized by temperature and moisture dependent on time according to specified ramps. In these composites, moisture diffusion takes place along the layers of the polymeric encapsulant, or along channel cracks in Silicon. Since their thickness is very small, it is possible to neglect moisture flux in the direction orthogonal to the EVA layers. Under these assumptions, moisture diffusion can be modelled as a diffusion process taking place over the mid-surface of a generic encapsulant layer, which corresponds to S( p).
Hence, the aim of the numerical method reduces to simulate and predict diffusion of water content c(x 1 , x 2 , t) along the encapsulant mid-surface S( p) for each point and time.
The following initial and boundary value problem can be considered, where an imposed water content c * is imposed on the boundary: where D is the encapsulant diffusivity. In addition to a reduced dimensionality of the domain where moisture diffusion takes place with respect to thermoelasticity, it is also remarkable to notice that these physical problems are characterized by very different time scales. The characteristic velocity of thermal diffusion is in fact ruled by the ratio k m /(ρ m c m ), while that of moisture diffusion is related to the diffusivity D. Considering characteristic values for EVA, the ratio between the velocities of the two phenomena is: so that heat transfer is about six order of magnitude faster than moisture diffusion. From this observation, we can state that moisture diffusion is dependent on heat transfer and not viceversa. Hence, the diffusivity of the encapsulant has to be considered as temperature dependent and, based on the experimental evidence [12], of Arrhenyus type: In order to take into account the effect of a possible debonding of the encapsulant, which would enhance moisture diffusion, D is assumed to be a linear increasing function of the normal gap

Weak forms
The partial differential equations governing the dynamic equilibrium of the body, Eq. (3), and heat conduction, Eq. (5), for each layer R(m) (m = 1, . . . , n) and the constitutive relations for the interfaces, Eqs. (8), (9) and (13) for each S( p), define an initial boundary value problem describing the debonding of a thermo-mechanical layered PV panel with thermo-visco-elastic polymeric interfaces.
Let t * I and u * I be the surface traction and the prescribed boundary displacement such that: and q * = q * I n I , and θ * be the imposed normal heat flux and the imposed temperature at the boundaries such that: where the indices ∂ D and ∂ N denote the Dirichlet and the Neumann portions of the boundary ∂ R.
The weak form corresponding to Eq. (3) is obtained by multiplying it by a virtual displacement δv I having a virtual gap [[δv I ]] along S( p) and by integrating the result on each domain R(m). After applying the divergence theorem as customary and by dropping the index m to simplify notation, we obtain: where δ I J is the Kronecker operator, and δ denotes a virtual variation of the variable to which is applied. Analogously, the weak form corresponding to Eq. (5) is obtained by multiplying it by a test function δθ having a gap [[δθ]] on S( p) and by integrating the result on each domain R(m). After some calculation and dropping the index m to simplify notation, we get: As far as moisture diffusion is concerned, the corresponding weak form is constructed by multiplying Eq. (15) by a test function δc. After integration by parts we have that the concentration c(x 1 , x 2 , t) solves the following equation for all the admissible δc and ∀t ∈ [0, t f ]: 4 Finite element discretization

Discretized weak forms for the thermo-elastic and heat conduction problems
According to the finite element method, the domain R is discretized into a finite number of bulk R e and interfaceS e elements so that: We also introduce for the purpose of numerical integration the mid-plane surface S( p) ≈ ∪ e S e , where S e is the middle surface of each interface elementS e .
The class of interface elements considered here consists of two surface elements coincident with the facets of the bulk elements used to discretize the continuum that are bricks or tetrahedra. For consistency between interfaces and bulk, the same order of interpolation is used. In the case of 2D simulations on cross-sections of a laminate, the present formulation still holds, provided that bulk elements are represented by quadrilateral or triangular plane strain finite elements and interface elements are given by two opposing lines. Again, the same interpolation order has to be used.
By introducing the shape functions, the finite element approximation for the bulk reads: a=1 defined in the natural reference system −1 ≤ ξ 1 , ξ 2 , ξ 3 ≤ +1, where N (e) is the number of element nodes, which is equal to 8 for a 3D linear brick element, or 4 for a 2D linear 4-node plane strain element.
Similarly, for the interface elements, the gaps are approximated as: where the shape functions {Ψ 1 (ξ 1 , ξ 2 )} S(e) a=1 are defined along the mid-surface plane in the natural reference system −1 ≤ ξ 1 , ξ 2 ≤ +1 and 2S(e) is the number of nodes of the interface element which is equal to 8 for a 3D interface element compatible with bricks, or 4 for a 2D interface element compatible with plane strain elements. The nodal displacement vector is: Similarly, the operator [ e ] ab applied to the nodal temperatures of the interface element leads to the temperature jumps between the (+) and the (−) interface flanks, i.e., After introducing these expressions in the weak form (17), its discretized version reads: Similarly, the discretized weak form (18) is: hence, Eq. (25) can be recast as: where the global mass, damping and stiffness matrices and the load vector are assembled as follows:

Discretized weak form of moisture diffusion
The discretization of the weak form for moisture diffusion is derived by introducing a finite element mesh of the midsurface S( p) of the encapsulant layer. In principle, since a staggered geometrical multiscale solution scheme is adopted, the spacing of the mesh used to solve moisture diffusion can be different from that used for the discretization of the thermo-mechanical problem. In that case, a projection of the nodal temperatures from the discretized thermo-mechanical problem to the nodes of the mesh used to solve moisture diffusion has to be performed via a suitable interpolation scheme.
In the sequel, without any loss of generality, we consider a finite element discretization for moisture diffusion coincident with the middle surface discretization of each interface element, i.e., S( p) ≈ e S e . By introducing the shape functions Ψ a , the water concentration in a generic point of coordinates (x 1 , x 2 ) and at a time t is Introducing Eq. (29) into (19), we obtain: the assembling of the global matrices leads to: providing a linear system of ordinary differential equations: Note that the matrix [B] is time dependent because it contains the diffusivity D which changes with time according to (16).

Proposed numerical solution scheme
To solve numerically the system (28) resulting from the thermo-mechanical problem, we adopt an Euler backward implicit scheme so that: where n = 1, 2, . . . , N and t n = t n−1 + t. Hence, Eq. (28) becomes: which is a nonlinear system of equations in the unknown { } n+1 , where the nonlinearity relies in the load vector due to the nonlinear relations between the cohesive tractions and the relative displacements, and between the heat flux and the temperature jump at the polymeric interfaces. This problem is solved iteratively using a Newton-Raphson scheme. At the iteration k + 1 we have an approximation { } n+1 (k) for { } n+1 and we seek for a better approximation { } n+1 (k+1) such that: By substituting this expression back to the system (35) and rearranging the various terms, leads: where the right-hand side is the so-called residual vector, {R} n+1 (k) , so that we can write: where: and [T e ] n+1 (k) is given by [P] The various terms entering Eq. (42) are reported in the Appendix.
Once temperature and displacements are computed at a given time step n + 1, these nodal results are transferred to the discretized moisture diffusion problem. Its solution is then performed by using an Euler backward time integration scheme with the same partition of the temporal interval [0, t f ] as for the thermo-mechanical problem. For n = 1, . . . , N time steps t, we solve the linear system of equations: are computed using the values Θ n+1 , [[U 3 ]] n+1 obtained from the thermo-mechanical problem. The algorithm for the proposed time stepping method with a staggered scheme is detailed in Algorithm 1. The Newton-Raphson iteration is performed until machine precision is achieved, i.e., up to a tolerance in the norm of the residual vector tol = 1 × 10 −15 .
It has to be remarked that the time-dependency of the visco-elastic constitutive equation (10) requires the use of a history variable {h v } for all the nodes of the finite element mesh for the thermo-mechanical problem. To model relaxation, this variable is set to zero at any change of temperature (state variable), while it is updated by the current time increment if the temperature remains constant with respect to the previous time step, within a given tolerance tol 2 . This method allows the simulation of the thermo-viscoelastic behavior of polymeric materials when the temperature-time superposition principle does not apply, for instance due to a change of microstructure by varying temperature as it happens in the case of epoxy-vinil-acetate used in photovoltaics.
Algorithm 1: Numerical scheme for the solution of the hygro-thermo-mechanical problem with a staggered approach.

Application to photovoltaics
In this section we propose the simulation of the two tests prescribed by international standars [4], namely the damp heat test and the humidity freeze test. While the former allows the complete uncoupling between the solution of the thermomechanical problem from moisture diffusion and allows the derivation of a closed form solution useful for benchmarking, the latter requires the present fully-coupled solution scheme. Moreover, the role of a temperature-dependent diffusion coefficient and the role of cohesive cracks in Silicon are investigated, in comparison to experimental results.

Damp heat test
Let us consider a laminate of span L = 125 cm and made of a Glass-Glass structure separated by EVA as in Fig. 6. The thickness of each glass is 3 mm, while the thickness of the EVA is 0.5 mm. In this laminate, moisture is diffusing from the free edges towards the centre, since glass is not permeable to moisture. As far as the initial and boundary conditions are concerned, let us consider the prescriptions by international standards [4] for the damp heat test, that is, a constant temperature of 85 • C and an air relative humidity of 85%. This relative humidity corresponds to a moisture content c * = 0.55 g/cm 3 imposed at the free edges of the laminate, i.e., for x 1 = 0 and x 1 = L, where x 1 is the distance from the left edge of the PV module.
Since temperature is held constant, the problem of diffusion can be solved independently from the thermomechanical problem, considering a constant diffusivity D = 5×10 −5 cm 2 /s corresponding to 85 • C. For this special case, the analytical solution was obtained in [12] and it is used as a benchmark for our computational scheme: The relative humidity in the air is kept constant at 85% for the range of temperatures where its control is thermodynamically feasible.
Due to a non constant temperature boundary condition, this problem is much more difficult to be simulated as compared to the damp heat test. In particular, the cohesive properties of EVA have to be updated during the simulation, as well as its diffusivity. More specifically, as far as the Young's modulus of EVA is concerned, the parameters α(T ) and a(T ) are herein considered to be temperature-dependent as experimentally evaluated in [6] and interpreted via a fractional calculus model in [21], see the plot for α(T ) and E(T ) in Fig. 9.
Regarding the diffusive properties, we consider the expression of D(T ) for EVA as reported in [12] and shown in Fig. 9 a EVA relaxation modulus versus time at various temperatures in a double logarithmic scale; b temperature dependent fractional exponent α(T ), adapted from [21] Fig. 10. Such trends can be fitted according to the Arrhenius type equation (16).
The critical crack opening, δ 3 , to be assigned to the cohesive zone model can be estimated from published experimental data in [13] reporting the variation of the Mode I fracture energy with temperature. Since the fracture energy G is the area below the traction-separation curve, the following relation holds: Hence, G(T ) experimental data can be converted in δ c 3 (T ) data based on the known temperature dependency of the Young's modulus of EVA as in Fig. 8a, evaluated for the asymptotic condition of an infinite time. Based on these data,  we obtain the correlations shown in Fig. 11 and used as input for the numerical simulations.
A sketch of the cross-section of a PV mini-module simulated in the present study and containing 3 solar cells is shown in Fig. 10. The lateral size of each Silicon cell is 125 mm and the the interspace between two cells is 2 mm. The module is made of a glass superstrate with thickness h G = 4 mm, an encapsulating polymer layer (EVA) with thickness h EVA = 0.5 mm, the Silicon solar cell with thickness h SI = 0.166 mm, another layer of EVA with the same thickness as the previous one, and finally a thin backsheet made of an ethylene tetrafluoroethylene core with silicon nitride coating (isovolta Icosolar T 2754), with thickness h BS = 0.1 mm. Thermal and mechanical parameters of each layer are collected in Table 1 and are taken from [5,6]. Since moisture penetrates from the backsheet and percolates along the interspace between solar cells, in the numerical simulations it is possible to impose a constant value of moisture concentration, c * , directly at the boundary of each solar cell embedded in the PV module (Fig. 12).
The temperature distribution inside a portion of the PV module cross-section near one of the free boundaries is shown in the contour plots in Fig. 13 for selected time steps. After the first ramp from 0 to 85 • C completed after 0.5 hours, heat has diffused inside the panel and temperature is almost uniform everywhere and equal to 85 • C. During the subsequent decreasing ramp from 85 to −40 • C, the Silicon cells and the EVA around them remain warmer than the other component. This temperature mismatch progressively shrinks during the further stage of the simulation at constant temperature θ * . This trend is quantified in Fig. 14 by plotting the temperature along a vertical line at x 1 = 2 mm from the free edge of the laminate, through the panel thickness.
The evolution of moisture concentration in the encapsulant vs. time by using a time-dependent diffusivity is shown in Fig. 15a. The same simulation with a constant diffusivity D = 5 × 10 −5 cm 2 /s corresponding to 85 • C is shown for comparison in Fig. 15b. As it can be noticed, the proper update of the diffusivity based on the actual temperature of EVA during the simulation provides very different results from those based on a constant diffusivity. In particular, moisture diffusion is a much slower phenomenon than what expected by the approach presented in [12], which is based on the use of a constant diffusivity equal to that at the maximum temperature experienced during the test.
The time step used for the simulation of the humidity freeze test is t = 18 s. At each time step, the Newton-Raphson iterative scheme used to solve the nonlinear thermo-mechanical block has a quadratic convergence and it requires a maximum of four iterations to satisfy the condition of a relative residual error norm less than the machine precision. For given variables computed from the solution of the thermo-mechanical block, the moisture diffusion block is linear and therefore it converges in one single iteration. In terms of CPU time, results are shown in Fig. 16, where the plot of the logarithm of CPU time for the solution of the nonlinear thermo-mechanical block, the moisture diffusion block, and the total time of the staggered scheme is shown for the first 90 time steps of the simulation of the humidity freeze test, corresponding to the first temperature ramp from 0 up to 85 • C. Simulations have been run on the server Proliant DL585R07 (128 GB Ram, 4 processors AMD Opteron 6282 SE 2.60 GHz, 16 cores). In the case of delamination,  The proposed numerical method can be also employed to simulate the effect of cracks in Silicon on the evolution of moisture concentration. To have a benchmark experimental result to compare with, minimodules composed of 3 × 3 multicrystalline solar cells with glass cover were realized in the Institute for Solar Energy Research of Hamelin, Germany (courtesy of Dr. M. Köentges) and were subjected to a threepoint bending test to induce cracks in the central solar cell, see Fig. 17.
This minimodule was then subjected to the humidity freeze test inside a climate chamber at Politecnico di Torino, Italy. Electroluminescence images taken regularly during the test (the reader is referred to [22] for more details about this nondestructive technique) show electrically damaged areas (black region) near the crack, see the EL image after 2400 h of testing in Fig. 18a. This electric degradation is presumably induced by chemical oxidation of the grid line deposited on the surface of the solar cell, enhanced by moisture. The presence of a crack appears to be relevant, since all the other solar cells show much lower degradation, mostly in form of dimmer areas at the edges of the solar cells due to moisture diffusion from the backsheet through the EVA interspace.
To prove that cracks can enhance electric degradation due to moisture diffusion through them, a numerical simulation is performed according to the present computational framework by imposing the value of c * = 0.055 g/cm 3 along the back side of the PV module. Diffusion can take place now not only along the EVA interspaces, but also along the two channel cracks, where 3D diffusive finite elements are inserted. Moisture concentration inside the EVA layer above the solar cells is shown in Fig. 18b after 2400 h of humidity freeze test. The presence of cracks enhances moisture diffusion in the central cell, with a contour plot of moisture concentration that correlates very well with the EL image of the electrically damaged areas (compare the red areas in Fig. 18b having c = c * with the black areas in Fig. 18a, confirming that electric degradation can be significantly enhanced by the presence of cracks.

Conclusions and outlook
A comprehensive finite element computational framework for the simulation of coupled hygro-thermo-mechanical problems in photovoltaic laminates has been proposed. To achieve the computational efficiency required to simulate large scale commercial PV modules consisting of up to 60 solar cells, the EVA encapsulant layers have been modelled as zero-thickness interface elements whose traction-separation relations take into account the complex thermo-visco-elastic rheological response of the polymer as per experiments. Moreover, a staggered solution scheme has been proposed and implemented in FEAP [23] to solve the partial differential equations governing heat transfer and thermo-elasticity in the three-dimensional domain of the laminate, and then predict moisture diffusion in the two dimensional domains represented by the EVA layers and channel cracks in Silicon by considering a material diffusivity dependent on temperature and crack opening. Fractional calculus has been used to model the visco-elastic behaviour of the EVA material layer, while moisture diffusion has been assumed to be of classical type. The use of fractional derivatives in partial differential equations (PDEs) related to diffusion is also a viable modelling approach and it is usually related to a diffusion process taking place across a non-Euclidean domain, see e.g. [24]. However, due to a lack of specific experimental evidence of fractal diffusive domains, the classic PDE for diffusion has been considered in the present study.
The proposed methodology has been successfully applied to the simulation of the qualification tests required by the International Electrotechnical Commission [4], namely the damp heat test and the humidity freeze test. In the latter, due to a continuous variation of temperature during the test, we have shown that coupling between the thermo-mechanical field and moisture diffusion has to be taken into account to correctly predict the spatio-temporal evolution of moisture in the PV module. Moreover, percolation through channel cracks in Silicon solar cells has been found to significantly enhance moisture diffusion. This trend has been confirmed by experimental tests performed by the present authors and showing an increased electric degradation in cracked solar cells with respect to the intact ones.
With the present computational tool available, further developments may regard the simulation of other moisture diffusion phenomena observed in experiments, in addition to percolation through Silicon channel cracks and to diffusion along the EVA layers. As shown in [2], in fact, imperfect sealing of the region near the main electric conductors (busbars) soldered on the surface of Silicon can further enhance percolation along them, creating preferred streams for moisture diffusion and chemical oxidation, see Fig. 19. To account for this type of degradation, a possible method is to model Fig. 19 Electroluminescence image of a PV module showing moisture degradation along the busbars, due to imperfect sealing (adapted from [2]) inhomogeneous diffusivity properties inside the EVA layer above the busbars by specifying an initial separation of the EVA interface. where: being the characteristic function such that [M] ab for (1 ≤ a ≤ S(e), 1 ≤ b ≤ 2S(e)) is the mean operator computing the average temperature across the EVA layer, i.e.