Research Article  Open Access
An Accurate Thermoviscoelastic Rheological Model for Ethylene Vinyl Acetate Based on Fractional Calculus
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 twoparameter model is proposed to approximate the powerlaw 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.
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 powerlaw type, as many other natural and artificial materials [20–23]. Classical rheological models used so far generally provide exponentialtype relations and, in order to approximate a powerlaw 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, powerlaw 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 twoparameter 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 Section 2, 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 stateoftheart 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 stressstrain () constitutive relation writeswith 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: 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 relaxationbased formulation, it is now useful to introduce the relaxation modulus for the two elements considered so far. While for the spring it is simply , for the dashpot it takes the form , with being the Diracfunction.
2.2. Maxwell Model and Its Generalization
The Maxwell model consists of a spring and a dashpot in series (Figure 1(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): By introducing the relaxation time , the corresponding relaxation function isThe Maxwell element is therefore based just on two parameters, and (or ). Indeed, a more accurate but complicated model can be obtained by arranging Maxwell elements (also called arms) in parallel, leading to the generalized model depicted in Figure 1(b). Note that an isolated spring with stiffness 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 takes the following form: which is referred to as a Prony series [30]. The difficulties in fitting the related parameters will be discussed in Section 4.
(a)
(b)
(c)
3. A Rheological Model Based on Fractional Calculus
The fractional generalization of the models presented in Section 2 can be achieved by replacing the firstorder derivative in (2) with the derivative of order [27]. According to Caputo’s definition, the αderivative of a generic function is given by [31] where is Euler’s Gamma function. By reminding that and , according to (6) the function itself is recovered for , while the firstorder classical derivative is obtained for . For quiescent systems at , 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 firstorder derivative with the derivative of order is equivalent, from a mechanical point of view, to replace a dashpot with a fractional dashpot (or springpot) of order . The simplest fractional element (also known as ScottBlair element) consists of a single springpot (Figure 1(c)). The constitutive equation writes ()Equation (7) tends to (1), that is, the classic spring element, for , while it provides (2), that is, the dashpot, for . As varies, the mechanical meaning of the parameter (with SI units Pa ) changes accordingly, passing from a stiffness () to a viscosity (). For this element, the relaxation modulus assumes the following powerlaw form:
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 powerlaw 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 Section 3 is more efficient to handle, since a powerlaw relaxation modulus (8) comes naturally from the constitutive law (7). Furthermore, only two coefficients 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 (and vice 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 (Figure 2). Predictions according to two different generalized Maxwell models, one with 3 arms () and one with 5 arms (), are firstly considered. For the implementation, the values (s) are chosen equidistant to cover the whole time range; that is, for and for . Given the times , the remaining coefficients (Pa) () are obtained by fitting the Prony series to the experimental curve, resulting in E = (0.6557, 0.3541, 0.3548, 0.3216) for and E = (0.7774, 0.1253, 0.1744, 0.1888, 0.1746, 0.1052) for .
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 are a = 1.544 Pa , (Table 1).

A 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 Figure 2 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 up to 5, but with the necessity of identifying 6 parameters. The simpler 2parameter springpot model proposed in this work, leading to a straight line with negative slope −α and intercept a, 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 Figure 3. The identified values are collected in Table 1, together with the mean absolute percentage error of the correlations, which is always very low.
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 timetemperature 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 timescale. 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 WilliamsLandelFerry (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 timetemperature 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 arms that was deemed necessary in [14], it would be necessary to determine fitting parameters for each testing temperature and then construct a lookup 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.
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 ( = −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 , 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 , depending on temperature: A similar analysis can be performed for the coefficient , 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.
Therefore, the limited amount of identified data in Table 1 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 timetemperature 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 GrunwaldLetnikov 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 timetemperature 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.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgments
Professor Paggi would like to acknowledge funding from the European Research Council under the European Union’s Seventh Framework Programme (FP/2007–2013)/ERC Grant Agreement no. 306622 (ERC Starting Grant “MultiField and MultiScale Computational Approach to Design and Durability of PhotoVoltaic Modules” (CA2PVM)). The research activity by Dr. Sapora has been supported by the Italian Ministry of Education, University and Research in the project FIRB 2010 Future in Research “Structural Mechanics Models for Renewable Energy Applications” (RBFR107AKG), which is also gratefully acknowledged.
References
 M. Paggi, S. KajariSchröder, and U. Eitner, “Thermomechanical deformations in photovoltaic laminates,” Journal of Strain Analysis for Engineering Design, vol. 46, no. 8, pp. 772–782, 2011. View at: Publisher Site  Google Scholar
 M. Paggi, M. Corrado, and M. A. Rodriguez, “A multiphysics and multiscale numerical approach to microcracking and powerloss in photovoltaic modules,” Composite Structures, vol. 95, pp. 630–638, 2013. View at: Publisher Site  Google Scholar
 O. Hasan and A. F. M. Arif, “Performance and life prediction model for photovoltaic modules: effect of encapsulant constitutive behavior,” Solar Energy Materials and Solar Cells, vol. 122, pp. 75–87, 2014. View at: Publisher Site  Google Scholar
 K. Agroui, A. Maallemi, M. Boumaour, G. Collins, and M. Salama, “Thermal stability of slow and fast cure EVA encapsulant material for photovoltaic module manufacturing process,” Solar Energy Materials and Solar Cells, vol. 90, no. 15, pp. 2509–2514, 2006. View at: Publisher Site  Google Scholar
 International Electrotechnical Commission, “Crystalline silicon terrestrial photovoltaic (PV) modules—design qualification and type approval,” International Standard IEC 61215:2006, 2006. View at: Google Scholar
 U. Eitner, Thermomechanics of photovoltaic modules [Ph.D. thesis], Zentrum für Ingenieurwissenschaften der MartinLutherUniversität HalleWittenberg, 2011.
 S. Dietrich, M. Pander, M. Sander, S. H. Schulze, and M. Ebert, “Mechanical and thermomechanical assessment of encapsulated solar cells by finiteelementsimulation,” in Reliability of Photovoltaic Cells, Modules, Components, and Systems III, vol. 7773 of Proceedings of SPIE, August 2010. View at: Publisher Site  Google Scholar
 M. Koehl, M. Heck, and S. Wiesmeier, “Modelling of conditions for accelerated lifetime testing of Humidity impact on PVmodules based on monitoring of climatic data,” Solar Energy Materials and Solar Cells, vol. 99, pp. 282–291, 2012. View at: Publisher Site  Google Scholar
 U. Eitner, M. Köntges, and R. Brendel, “Use of digital image correlation technique to determine thermomechanical deformations in photovoltaic laminates: measurements and accuracy,” Solar Energy Materials and Solar Cells, vol. 94, no. 8, pp. 1346–1351, 2010. View at: Publisher Site  Google Scholar
 U. Eitner, S. KajariSchröder, M. Köntges, and R. Brendel, “Nonlinear mechanical properties of ethylenevinyl acetate (EVA) and its relevance to thermomechanics of photovoltaic modules,” in Proceedings of the 25th European Conference on Photovoltaic Solar Energy, pp. 4366–4368, 2010. View at: Google Scholar
 J. Schlothauer, S. Jungwirth, M. Köhl, and B. Röder, “Degradation of the encapsulant polymer in outdoor weathered photovoltaic modules: spatially resolved inspection of EVA ageing by fluorescence and correlation to electroluminescence,” Solar Energy Materials & Solar Cells, vol. 102, pp. 75–85, 2012. View at: Publisher Site  Google Scholar
 A. Sapora and M. Paggi, “A coupled cohesive zone model for transient analysis of thermoelastic interface debonding,” Computational Mechanics, vol. 53, no. 4, pp. 845–857, 2014. View at: Publisher Site  Google Scholar  MathSciNet
 F. J. Pern, “Factors that affect the EVA encapsulant discoloration rate upon accelerated exposure,” Solar Energy Materials and Solar Cells, vol. 4142, pp. 587–615, 1996. View at: Publisher Site  Google Scholar
 M. Pander, S. Dietrich, S.H. Schulze, U. Eitner, and M. Ebert, “Thermomechanical assessment of solar cell displacement with respect to the viscoelastic behaviour of the encapsulant,” in Proceedings of the 12th International Conference on Thermal, Mechanical and MultiPhysics Simulation and Experiments in Microelectronics and Microsystems (EuroSimE '11), pp. 1–6, Linz, Austria, April 2011. View at: Publisher Site  Google Scholar
 E. Kaplani, “Detection of degradation effects in fieldaged cSi solar cells through IR thermography and digital image processing,” International Journal of Photoenergy, vol. 2012, Article ID 396792, 11 pages, 2012. View at: Publisher Site  Google Scholar
 N. C. Park, W. W. Oh, and D. H. Kim, “Effect of temperature and humidity on the degradation rate of multicrystalline silicon photovoltaic module,” International Journal of Photoenergy, vol. 2013, Article ID 925280, 9 pages, 2013. View at: Publisher Site  Google Scholar
 T. Kojima and T. Yanagisawa, “The evaluation of accelerated test for degradation a stacked aSi solar cell and EVA films,” Solar Energy Materials and Solar Cells, vol. 81, no. 1, pp. 119–123, 2004. View at: Publisher Site  Google Scholar
 D. Berman and D. Faiman, “EVA browning and the timedependence of IV curve parameters on PV modules with and without mirrorenhancement in a desert environment,” Solar Energy Materials and Solar Cells, vol. 45, no. 4, pp. 401–412, 1997. View at: Publisher Site  Google Scholar
 F. D. Novoa, D. C. Miller, and R. H. Dauskardt, “Environmental mechanisms of debonding in photovoltaic backsheets,” Solar Energy Materials and Solar Cells, vol. 120, part A, pp. 87–93, 2014. View at: Publisher Site  Google Scholar
 P. G. Nutting, “A new general law of deformation,” Journal of the Franklin Institute, vol. 191, no. 5, pp. 679–685, 1921. View at: Publisher Site  Google Scholar
 R. Hilfer, Fractional Calculus in Physics, World Scientific, Singapore, 2000.
 R. L. Magin, Fractional Calculus in Bioengineering, Begel House, Redding, Calif, USA, 2006.
 J. T. Machado, V. Kiryakova, and F. Mainardi, “Recent history of fractional calculus,” Communications in Nonlinear Science and Numerical Simulation, vol. 16, no. 3, pp. 1140–1153, 2011. View at: Publisher Site  Google Scholar  MathSciNet
 J. Sorvari and M. Malinen, “On the direct estimation of creep and relaxation functions,” Mechanics of TimeDependent Materials, vol. 11, no. 2, pp. 143–157, 2007. View at: Publisher Site  Google Scholar
 G. Samko, A. A. Kilbas, and O. I. Marichev, Fractional Integrals and Derivatives, Gordon and Breach, Amsterdam, The Netherlands, 2000.
 H. Schiessel, R. Metzler, A. Blumen, and T. F. Nonnenmacher, “Generalized viscoelastic models: their fractional equations with solutions,” Journal of Physics A: General Physics, vol. 28, no. 23, pp. 6567–6584, 1995. View at: Publisher Site  Google Scholar
 F. Mainardi, Fractional Calculus and Waves in Linear Viscoelasticity. An Introduction to Mathematical Models, Imperial College Press, London, UK, 2010. View at: Publisher Site  MathSciNet
 M. di Paola, A. Pirrotta, and A. Valenza, “Viscoelastic behavior through fractional calculus: an easier method for best fitting experimental results,” Mechanics of Materials, vol. 43, no. 12, pp. 799–806, 2011. View at: Publisher Site  Google Scholar
 M. di Paola, F. P. Pinnola, and M. Zingales, “A discrete mechanical model of fractional hereditary materials,” Meccanica, vol. 48, no. 7, pp. 1573–1586, 2013. View at: Publisher Site  Google Scholar  MathSciNet
 N. W. Tschoegl, The Phenomenological Theory of Linear Viscoelastic Behavior: An Introduction, Springer, Berlin, Germany, 1989. View at: Publisher Site  MathSciNet
 I. Podlubny, Fractional Differential Equations, vol. 198 of Mathematics in Science and Engineering, Academic Press, San Diego, Calif, USA, 1999. View at: MathSciNet
 R. L. Bagley, “Power law and fractional calculus model of viscoelasticity,” AIAA Journal, vol. 27, no. 10, pp. 1412–1417, 1989. View at: Publisher Site  Google Scholar
 A. Carpinteri and F. Mainardi, Fractals and Fractional Calculus in Continuum Mechanics, Springer, Wien, Germany, 1997. View at: Publisher Site  MathSciNet
 R. D. Bradshaw and L. C. Brinson, “A sign control method for fitting and interconverting material functions for linearly viscoelastic solids,” Mechanics of TimeDependent Materials, vol. 1, no. 1, pp. 85–108, 1997. View at: Publisher Site  Google Scholar
 Y. Eom, L. Boogh, V. Michaud, P. Sunderland, and J.A. Månson, “Timecuretemperature superposition for the prediction of instantaneous viscoelastic properties during cure,” Polymer Engineering and Science, vol. 40, no. 6, pp. 1281–1292, 2000. View at: Publisher Site  Google Scholar
 W. Stark and M. Jaunich, “Investigation of Ethylene/Vinyl Acetate Copolymer (EVA) by thermal analysis DSC and DMA,” Polymer Testing, vol. 30, no. 2, pp. 236–242, 2011. View at: Publisher Site  Google Scholar
 S. Müller, M. Kästner, J. Brummund, and V. Ulbricht, “On the numerical handling of fractional viscoelastic material models in a FE analysis,” Computational Mechanics, vol. 51, no. 6, pp. 999–1012, 2013. View at: Publisher Site  Google Scholar  MathSciNet
 K. B. Oldham and J. Spanier, The Fractional Calculus, Academic Press, New York, NY, USA, 1974. View at: MathSciNet
 M. Simon and E. L. Meyer, “Detection and analysis of hotspot formation in solar cells,” Solar Energy Materials and Solar Cells, vol. 94, no. 2, pp. 106–113, 2010. View at: Publisher Site  Google Scholar
 M. Paggi and A. Sapora, “Numerical modelling of microcracking in PV modules induced by thermomechanical loads,” Energy Procedia, vol. 38, pp. 506–515, 2013. View at: Publisher Site  Google Scholar
Copyright
Copyright © 2015 Marco Paggi and Alberto Sapora. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.