# An Accurate Thermoviscoelastic Rheological Model for Ethylene Vinyl Acetate Based on Fractional Calculus

Authors: Marco Paggi; Alberto Sapora
Journal: International Journal of Photoenergy
 (2015)

---

## Abstract

The thermoviscoelastic rheological properties of ethylene vinyl acetate (EVA) used to embed solar cells have to be accurately described to assess the deformation and the stress state of photovoltaic (PV) modules and their durability. In the present work, considering the stress as dependent on a noninteger derivative of the strain, a two-parameter model is proposed to approximate the power-law relation between the relaxation modulus and time for a given temperature level. Experimental validation with EVA uniaxial relaxation data at different constant temperatures proves the great advantage of the proposed approach over classical rheological models based on exponential solutions.


---

## Body

## 1. Introduction

Ethylene vinyl acetate (EVA) is the most used material for encapsulating solar cells in photovoltaic (PV) modules due to its low cost, high transmittance of light, good thermal conduction, and long operating range. This polymer is used for binding the glass cover, the solar cells, and the backsheet together in order to realize the layered stack of a typical PV module [1, 2]. Moreover, EVA provides protection of cells and interconnections from moisture, foreign impurities, and mechanical damage. The provision of such properties makes it a vital component on which the performance of a PV module depends [3].In this context, modeling accurately the thermoviscoelastic behavior of EVA over the whole thermal range used in applications, which varies from −40°C, the lowest temperature for standard qualification testing of PV modules [4], up to +150°C, corresponding to the lamination temperature [5], is a crucial and not trivial task. The elastic modulus may vary up to three orders of magnitude depending on the temperature and for temperatures less than zero, it may vary up to one order of magnitude depending on the relaxation time [2, 6]. Eva constitutive behavior significantly influences the stress and deformation states of solar cells [1, 2, 6–8]. It is responsible for the nonlinear variation of the gap between solar cells during temperature excursions experienced in the field, which may lead to the interconnection failure [9–13]. EVA is also subjected to degradation phenomena and discoloration due to aging [14–18] and, as shown very recently, its properties can have a role in the mechanism of backsheet delamination [19], influencing the toughness of the interface between backsheet and the polymer.EVA shows a relaxation behavior of power-law type, as many other natural and artificial materials [20–23]. Classical rheological models used so far generally provide exponential-type relations and, in order to approximate a power-law trend, a huge number of elements (and thus of model parameters) should be taken into account. Generalized Maxwell models or Prony series have been usually proposed to fit EVA experimental relaxation curves and have been implemented in commercial finite element software like Abaqus, Ansys, or Comsol multiphysics as user defined elements. A Prony series with up to 100 Maxwell arms in parallel was considered in [6], whereas models with up to 32 arms were used in [14].This gives rise to several drawbacks, including the fact that it is not always possible to provide a clear mechanical meaning to all the parameters and that the numerical procedure for their identification is not trivial, being these subjected to several physical constraints [24].On the other hand, power-law relations rise out naturally by assuming a material constitutive law of fractional type, that is, involving noninteger order derivatives of stress and/or strain [25, 26]. The ability of fractional calculus to model hereditary phenomena with long memory has now been widely accepted by the scientific community [27–29].In the present work, a two-parameter element consisting of a fractional dashpot (or springpot) is considered to model the viscoelastic behavior of EVA. The aim is to determine a simple constitutive model of easy numerical implementation able to capture the material mechanical behavior at all temperatures and with only two free parameters to be identified. The fractional approach will be successfully applied to relaxation tests carried out on EVA samples at different constant temperatures ranging from −35°C to 139°C [6] and used to assess the evolution of the gap between solar cells [1, 9]. The advantages with respect to Prony series in terms of both accuracy and model parameters identification are proven.The paper is organized as follows: classical rheological models are summarized in Section2, while the proposed fractional approach is described in Section 3. Section 4 includes details on the numerical procedure to fit the model parameters, the comparison between experimental data and theoretical results and comments on the relevance of the proposed constitutive model for computational simulations on PV modules. Concluding remarks are given in Section 5.


## 2. Classical Rheological Models

Basic models in the classical rheology are herein briefly summarized, since they constitute the actual state-of-the-art modelling of the viscoelastic behavior of EVA.


### 2.1. Spring and Dashpot Elements

The material elastic behavior is classically represented in mechanics by a spring element, whose stress-strain (σ-ε) constitutive relation writes(1)σt=Eεt,with E being the spring stiffness.On the other hand, the limit case of a fluid viscous response is described by a dashpot element, whose constitutive law can be expressed:(2)σt=ηdεtdt,where η is the dynamic viscosity.Viscoelastic materials exhibit both elastic and viscous behaviors, which are then modeled by properly combining springs and dashpots. Since the viscoelastic constitutive equation will be expressed by using a relaxation-based formulation, it is now useful to introduce the relaxation modulus for the two elements considered so far. While for the spring it is simplyE(t)=E, for the dashpot it takes the form E(t)=ηδ(t), with δ(t) being the Dirac-function.


### 2.2. Maxwell Model and Its Generalization

The Maxwell model consists of a spring and a dashpot in series (Figure1(a)). Since the two elements are subjected to the same stress, while the total strain is the sum of two contributions—one of the spring and one of the dashpot, respectively—the constitutive relation is obtained by differentiating (1) and adding the result to (2): (3)σt+ηEdσtdt=ηdεtdt.By introducing the relaxation time τ=η/E, the corresponding relaxation function is(4)Et=Eexp-tτ.The Maxwell element is therefore based just on two parameters, E and η (or τ). Indeed, a more accurate but complicated model can be obtained by arranging n Maxwell elements (also calledarms) in parallel, leading to the generalized model depicted in Figure 1(b). Note that an isolated spring with stiffness E0 is also introduced, in order for the relaxation function to achieve an asymptotic value different from zero. This choice is supported by several experimental data available in the literature [10]. The function E(t) takes the following form: (5)Et=E0+∑i=1nEiexp-tτi,which is referred to as a Prony series [30]. The difficulties in fitting the related 2n+1 parameters will be discussed in Section 4.Figure 1
Maxwell model (a); generalized Maxwell model (b); Scott-Blair element or springpot (c).
(a)
(b)
(c)


## 2.1. Spring and Dashpot Elements

The material elastic behavior is classically represented in mechanics by a spring element, whose stress-strain (σ-ε) constitutive relation writes(1)σt=Eεt,with E being the spring stiffness.On the other hand, the limit case of a fluid viscous response is described by a dashpot element, whose constitutive law can be expressed:(2)σt=ηdεtdt,where η is the dynamic viscosity.Viscoelastic materials exhibit both elastic and viscous behaviors, which are then modeled by properly combining springs and dashpots. Since the viscoelastic constitutive equation will be expressed by using a relaxation-based formulation, it is now useful to introduce the relaxation modulus for the two elements considered so far. While for the spring it is simplyE(t)=E, for the dashpot it takes the form E(t)=ηδ(t), with δ(t) being the Dirac-function.


## 2.2. Maxwell Model and Its Generalization

The Maxwell model consists of a spring and a dashpot in series (Figure1(a)). Since the two elements are subjected to the same stress, while the total strain is the sum of two contributions—one of the spring and one of the dashpot, respectively—the constitutive relation is obtained by differentiating (1) and adding the result to (2): (3)σt+ηEdσtdt=ηdεtdt.By introducing the relaxation time τ=η/E, the corresponding relaxation function is(4)Et=Eexp-tτ.The Maxwell element is therefore based just on two parameters, E and η (or τ). Indeed, a more accurate but complicated model can be obtained by arranging n Maxwell elements (also calledarms) in parallel, leading to the generalized model depicted in Figure 1(b). Note that an isolated spring with stiffness E0 is also introduced, in order for the relaxation function to achieve an asymptotic value different from zero. This choice is supported by several experimental data available in the literature [10]. The function E(t) takes the following form: (5)Et=E0+∑i=1nEiexp-tτi,which is referred to as a Prony series [30]. The difficulties in fitting the related 2n+1 parameters will be discussed in Section 4.Figure 1
Maxwell model (a); generalized Maxwell model (b); Scott-Blair element or springpot (c).
(a)
(b)
(c)


## 3. A Rheological Model Based on Fractional Calculus

The fractional generalization of the models presented in Section2 can be achieved by replacing the first-order derivative in (2) with the derivative of order α∈(0,1) [27]. According to Caputo’s definition, the α-derivative of a generic function f is given by [31] (6)dαfxdxα=1Γ1-α∫0xf′yx-yαdy,where Γ is Euler’s Gamma function. By reminding that Γ(0+)=+∞ and Γ(1)=1, according to (6) the function f itself is recovered for α=0, while the first-order classical derivative is obtained for α=1. For quiescent systems at x=0, the definition (6) coincides with that provided by Riemann and Liouville [31]. We will thus refer with the same symbol to both definitions, without losing of generality.To substitute the first-order derivative with the derivative of orderα∈(0,1) is equivalent, from a mechanical point of view, to replace a dashpot with a fractional dashpot (orspringpot) of order α. The simplest fractional element (also known as Scott-Blair element) consists of a single springpot (Figure 1(c)). The constitutive equation writes (α∈(0,1))(7)σt=adαεtdtα.Equation (7) tends to (1), that is, the classic spring element, for α→0+, while it provides (2), that is, the dashpot, for α→1. As α varies, the mechanical meaning of the parameter a (with SI units Pa sα) changes accordingly, passing from a stiffness (α=0) to a viscosity (α=1). For this element, the relaxation modulus E(t) assumes the following power-law form:(8)Et=at-αΓ1-α.Equation (8) will be sufficient for the purposes of the present study, although the singular nature of the power law’s relationship, that is, infinite modulus at time zero, is inappropriate for the goal of portraying the phenomenon at times close to the onset of the strain. This point was deeply discussed in [32].The analysis of fractional viscoelastic models is limited to what is presented above, for the sake of simplicity. Interested readers can refer to [27, 33] for a deeper insight on this topic. It is however important to realize that the substitution of a dashpot by a springpot involves an additional parameter in the modeling analysis, for a total of two parameters to be identified.


## 4. Experimental Validation of the Proposed Approach for EVA and Critical Comparison with Classical Methods

Several experimental tests on viscoelastic materials show relaxation functions with a power-law behavior [20, 28]. This is also the case of EVA; see the recent experimental campaign in [6]. If a classical rheological model is used, since an exponential function (4) is obtained through the Maxwell element, this implies the need of implementing generalized Prony series (5). Indeed, several drawbacks emerge from this approach. First of all, it is not always clear what the physical meaning of all the arms is. Secondly, several restrictions must be imposed for physically realistic materials; for example, the coefficients have to be strictly positive [34]. The fitting algorithm results therefore into a constrained least squares problem, and sophisticated numerical methods must be implemented; see, for example, [24, 30].Eventually, notice that a creep experiment is usually easier to perform than a relaxation one, so most of the experimental data are available in the form of creep compliance versus time. Since there is no known closed form for the (creep) compliance in terms of the coefficients of the Prony series, once the creep data are known, it is not straightforward at all to get the coefficients of the (relaxation) Prony series.On the other hand, the fractional model discussed in Section3 is more efficient to handle, since a power-law relaxation modulus (8) comes naturally from the constitutive law (7). Furthermore, only two coefficients (a,α) need to be fitted. Observe also that fractional operators have a simple definition in terms of the Laplace transforms, so it is not difficult to obtain the creep compliance starting from the relaxation modulus (andvice versa), with these functions being strictly related to each other in the Laplace domain [20, 25].Generalized Maxwell models and the proposed fractional element are now applied to the results of the uniaxial relaxation tests on EVA specimens carried out in [6] at different temperatures. The reader is referred to [6, 10] for a description on the testing procedure, which is standard, and it is not reported here to avoid duplication of information.Let us start by considering the relaxation data of EVA tested at 60°C (Figure2). Predictions according to two different generalized Maxwell models, one with 3 arms (n=3) and one with 5 arms (n=5), are firstly considered. For the implementation, the values τi (s) are chosen equidistant to cover the whole time range; that is, τ=(100,102,104) for n=3 and τ=(100,101,102,103,104) for n=5. Given the times τi, the remaining n+1 coefficients Ej (Pa) (j=0,…,n) are obtained by fitting the Prony series to the experimental curve, resulting inE = (0.6557, 0.3541, 0.3548, 0.3216) for n=3 andE = (0.7774, 0.1253, 0.1744, 0.1888, 0.1746, 0.1052) for n=5.Figure 2
Relaxation modulus versus time for EVA at 60°C: experimental data (circles), predictions by the generalized Maxwell model with 3 arms (dashed line), by the generalized Maxwell model with 5 arms (dotted line), and by the springpot model based on fractional calculus (continuous line).A similar fitting is carried out for the fractional model (8), where a procedure to evaluate the coefficients of a nonlinear regression function, using a least squares estimate, is adopted. The identified parameters area = 1.544 Pa sα, α=0.06542 (Table 1).Table 1
Identified parameters of the proposed rheological model based on fractional calculus used to fit experimental uniaxial relaxation tests at different temperatures.
Temperature(°C)
α
a (Pa sα)
Mean absolutepercentage error (%)
−35
0.22600
814.7
3.728
−28
0.16810
182.7
2.702
−18
0.10150
52.63
1.823
0
0.05566
23.55
1.851
20
0.04227
11.04
0.3044
40
0.07417
4.668
2.977
49
0.08634
4.116
5.467
60
0.06542
1.544
0.9898
80
0.05117
1.049
1.110
100
0.04179
0.9276
0.8064
119
0.03610
0.7965
0.9627
139
0.03311
0.8228
0.3811A comparison between the generalized Maxwell models with 3 or 5 arms, the proposed rheological model based on fractional calculus, and the experimental results is shown in Figure2 on a bilogarithmic plane. The generalized Maxwell model with 3 arms, often available in the majority of finite element programs, like the Finite Element Analysis Program (FEAP) by Zienkiewicz and Taylor, without the need of implementing new user element subroutines, generates a wavy curve with very low accuracy. More accurate predictions are achieved by increasing n up to 5, but with the necessity of identifying 6 parameters. The simpler 2-parameter springpot model proposed in this work, leading to a straight line with negative slope −α and intercepta, is indeed reasonably accurate and simple.With the efficiency of the proposed approach over the classic one being proven, the identification procedure for the fractional model is repeated here for relaxation data related to different testing temperatures; see Figure3. The identified values (a,α) are collected in Table 1, together with the mean absolute percentage error of the correlations, which is always very low.Figure 3
Relaxation modulus versus time for EVA specimens at different temperatures: experimental data (circles) and predictions by the springpot model based on fractional calculus (continuous lines).It is remarkable to note here that this further identification procedure is not usually done in the framework of the classic approaches which rely on the time-temperature superposition principle [35]. It states that the material response at a given temperature can be related to that at another temperature by a suitable change in time-scale. Actually, when this principle is checked against experimental results [14], its validity for EVA was found to be questionable and leading to poor predictions if tacitly assumed to be valid (see the comparison in Figure 4 between the actual shift factor deduced by experiments and the Williams-Landel-Ferry (WLF) theoretically expected shift function). This behavior is imputable to a modification of the EVA microstructure due to melting of the crystalline part, as reported in [14]. Therefore, the use of the classical principle of time-temperature superposition should not be applied to model the thermoviscoelastic behavior of EVA over the whole temperature range. In spite of this evidence, the problem was left unresolved up to now in the framework of a Prony series representation. In fact, for a model with n=25 arms that was deemed necessary in [14], it would be necessary to determine n+1=26 fitting parameters for each testing temperature and then construct a look-up table from where different model parameters can be picked up during numerical simulations. Clearly, this leads to a complex numerical implementation that cannot be handled by commercial finite element software. Moreover, interpolation between Prony series approximations corresponding to different temperatures is not an obvious procedure as well.Figure 4
Comparison between the actual shift factor determined from EVA experiments (dots) and the expected Williams-Landel-Ferry WLF shift function (solid curve) for the application of the time-temperature superposition principle (adapted from [14]). The reference temperature is Tref = −20°C.By using the proposed approach, only 2 model parameters have to be determined for each testing temperature, which makes it possible to determine useful correlations of much easier implementation in finite element software. In particular, by plotting the exponentα versus ΔT=T-Tref (Tref = −35°C) in Figure 5 (using the data collected in Table 1), we note that the transition in the EVA microstructure due to melting of the crystalline phase is well evidenced by this method due to the deviation from a single trend and it occurs at ΔT~80°C, consistently with that reported in [14] (see also [36]). Therefore, within the present approach, an accurate modelling of the rheological behavior of EVA suggests the use of the following two distinct correlations for the parameter α=α(T), depending on temperature: (9)αT=-6×10-7ΔT3+1×10-4ΔT2-0.0093ΔT+0.225,ΔT≤80°C,-1×10-7ΔT3+5×10-5ΔT2-0.0083ΔT+0.474,ΔT>80°C.A similar analysis can be performed for the coefficient a=a(T), which generally decreases by increasing temperature; see Table 1. This physically implies that the stiffness of the material increases as temperature decreases, consistently with the experimental evidence.Figure 5
Variation of the identified model parameterα versus temperature. Two different correlations are suggested due to a modification in the material microstructure of EVA for ΔT~80°C. The reference temperature is Tref = −35°C.R2 represents the coefficient of determination of the best-fit.Therefore, the limited amount of identified data in Table1 or, alternatively, correlations of the type (9) can be efficiently passed as input to a new user material element based on the proposed approach for accurate finite element simulations, without the need of relying on the time-temperature superposition principle. The implementation of the new constitutive equation (7) in the finite element method [37] via a user defined material will require the numerical approximation of the fractional derivative (6). Starting from the simple Grunwald-Letnikov definition, different algorithms based on fractional finite differences can be considered [38] and their comparison is left for further research.The rheological model proposed in the present work can be very useful for improving the accuracy of numerical simulations of coupled problems in multiphysics whenever time and temperature are varying at the same time [1, 2]. Examples of this type where the proposed methodology is expected to be beneficial are the simulation of the process of cooling down of PV modules to the ambient temperature from the lamination temperature in order to assess residual stresses; the simulation of aging of PV modules and crack propagation in Silicon during standard tests inside a climate chamber; temperature degradation of EVA in case of hot spots or defects [39, 40]. In all of these cases, in fact, the temperature is varying with position and time inside the PV module and an accurate instantaneous value of EVA Young’s modulus has to be used in each Gauss point and at each integration time step.


## 5. Conclusions

In this study, a rheological constitutive relationship based on fractional calculus has been proposed to model the viscoelastic behavior of EVA. The formulation involves just two parameters, the identification of which is much easier to be made than classical Prony series related to standard approaches. The relaxation modulus, as function of time, has been evaluated over the whole temperature range important in PV applications: theoretical results have been compared with experimental data, showing a very satisfactory agreement.The present modeling allows bypassing the problem of inapplicability of the time-temperature superposition principle for EVA noticed in the literature, since a very limited number of identified model parameters (or just correlations, see (9)) can be passed as input to a user defined material subroutine for finite element software.



---

*Source: 252740-2015-06-15.xml*