Multi-field and multi-scale Modeling of Fracture for renewable energy Applications
Abstract
This thesis is mainly focused on the computational modeling of solar cell cracking, multiphyics phenomena, and recycling of photovoltaic (PV) modules through the finite element method. Specifically, it consists of three parts. In the first part, a comprehensive hygro-thermo-mechanical computational framework in the 3D setting is proposed to model the coupled degradation phenomena in the PV modules for the durability analysis, and it 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. The second part is focused on the crack modeling of very thin and brittle silicon solar cells in the PV modules, and a reliable computational framework integrating solid shell element formulation with phase field fracture modeling is developed using the efficient quasi-Newton solution scheme and global local approach. The excellent performance is showcased through the simulation of different boundary value problems, and then applied to predict the crack growth of silicon solar cells when the PV modules are subjected to different external loadings. The third part addresses the efficient recycling of PV modules through the numerical modeling method by the development of 3D interface finite element with humidity-dose enhanced cohesive zone model for the peeling simulation to separate different layers, and diffusion-swelling large deformation continuum theory for the nondestructive recovery of silicon cells in the PV recycling using the solvent method. With these tools at hand, it is possible to design suitable virtual testing procedures for PV durability and recyclability analysis.
Full text
IMT School for Advanced Studies, Lucca Lucca, Italy Multi-field and multi-scale modeling of fracture for renewable energy applications PhD Program in Systems Science Track in Computer Science and Systems Engineering XXXVI Cycle By Zeng Liu 2023
The dissertation of Zeng Liu is approved. PhD Program Coordinator: Alberto Bemporad, IMT School for Advanced Studies Lucca Advisor: Prof. Marco Paggi, IMT School for Advanced Studies Lucca Co-Advisor: Prof. Jose Antonio Reinoso Cuevas, Universidad de Sevilla The dissertation of Zeng Liu has been reviewed by: Professor 1, University 1 Professor 2, University 2 IMT School for Advanced Studies Lucca 2023
TBD
Contents List of Figures x List of Tables xx Acknowledgements xxi Abstract xxiii 1 Introduction 1 1.1 Introduction to photovoltaic modules . . . . . . . . . . . . 1 1.2 Durability of photovoltaic modules: Accelerated hygrothermo-mechanical ageing . . . . . . . . . . . . . . . . . . . 4 1.3 Durability of photovoltaic modules: Silicon cell cracking . 8 1.4 Sustainability of photovoltaic modules: Recycling of the end-of-life PV modules . . . . . . . . . . . . . . . . . . . . . 11 1.5 Outline of the thesis . . . . . . . . . . . . . . . . . . . . . . . 13 2 Multiphysics modeling for the durability analysis of photovoltaic laminates 15 2.1 Hygro-thermo-mechanical modeling framework for the PV laminate ............................. 19 2.1.1 Coupled thermo-mechanical kinematics formulation 19 2.1.2 Thermo-visco-elastic cohesive interface model for polymericlayers .................... 24 2.1.3 Moisture diffusion along polymeric interfaces . . . 26 2.1.4 Variational form and finite element approximation 27 vii
2.1.5 Staggered solution scheme for the multi-field framework ........................... 41 2.2 Numerical applications to photovoltaics . . . . . . . . . . . 44 2.2.1 Damp heat test . . . . . . . . . . . . . . . . . . . . . 45 2.2.2 Humidity freeze test . . . . . . . . . . . . . . . . . . 54 2.2.3 Thermal cycling test . . . . . . . . . . . . . . . . . . 59 3 Global-local phase field modeling of silicon cell cracking in photovoltaic laminates 73 3.1 Phase field solid shell formulation . . . . . . . . . . . . . . 76 3.1.1 Basics of phase field approach . . . . . . . . . . . . 76 3.1.2 Solid shell kinematics . . . . . . . . . . . . . . . . . 79 3.2 Numerical implementation . . . . . . . . . . . . . . . . . . 81 3.2.1 Finite element interpolation . . . . . . . . . . . . . . 81 3.2.2 Solution schemes . . . . . . . . . . . . . . . . . . . . 87 3.3 The global-local fracture submodeling approach . . . . . . 89 3.4 Numerical examples . . . . . . . . . . . . . . . . . . . . . . 92 3.4.1 Single edge notched tension . . . . . . . . . . . . . . 93 3.4.2 Single edge notched shear . . . . . . . . . . . . . . . 95 3.4.3 Cylindrical structure under mixed tension and bending ............................ 98 3.4.4 Phase field fatigue modeling using solid shell formulation.........................102 3.4.5 Global-local fracture modeling of thin-walled structure with curved geometry . . . . . . . . . . . . . . 107 3.4.6 Benchmark example of anisotropic fracture . . . . . 111 3.4.7 Silicon solar cell cracking in the photovoltaic module114 4 Modeling and simulation of the recycling process for end-of- life photovoltaic laminates 124 4.1 Coupled modeling of polymer-solvent system in PV recycling................................128 4.1.1 Continuum kinematics . . . . . . . . . . . . . . . . . 128 4.1.2 Balancelaws.......................130 4.1.3 Thermodynamics and constitutive laws . . . . . . . 131 viii
4.1.4 Free energy definition . . . . . . . . . . . . . . . . . 134 4.1.5 Numerical implementation . . . . . . . . . . . . . . 136 4.1.6 Numerical examples . . . . . . . . . . . . . . . . . . 143 4.2 Peeling simulation of PV recycling . . . . . . . . . . . . . . 163 4.2.1 Humidity dose-based cohesive zone model . . . . . 163 4.2.2 3D large deformation interface element and FE implementation ......................167 4.2.3 Numerical examples . . . . . . . . . . . . . . . . . . 174 5 Conclusion and future developments 187 5.1 Summary and further developments for the multiphysics modeling of photovoltaic modules . . . . . . . . . . . . . . 187 5.2 Concluding remarks and further developments for the globallocal crack modeling of solar cells in the photovoltaic modules ................................189 5.3 Final conclusions and futher developments for the recycling of photovoltaic modules . . . . . . . . . . . . . . . . . 190 ix
69 The predicted phase field contour plots of the local model with no fracture anisotropy at the loading displacement 1.24 mm, 1.34 mm, 1.42 mm, and 1.48 mm. . . . . . . . . . 118 70 The predicted phase field contour plots of the local model with the preferential crack orientation angle equal to 45 degree at the loading displacement of 1.24 mm, 1.34 mm, 1.42 mm, and 1.48 mm. . . . . . . . . . . . . . . . . . . . . . 119 71 Schematic diagram of the local fracture modeling of one single solar cell when the photovoltaic module is subjected to three point bending. . . . . . . . . . . . . . . . . . . . . . 120 72 The predicted displacement contour plots of the whole photovoltaic module when subjected to three point bending at the loading displacement 67.5 mm, 70.5 mm, 73.5 mm, and 76.5mm. .............................121 73 The predicted phase field contour plots of the local model with no fracture anisotropy at the loading displacement 67.5 mm, 70.5 mm, 73.5 mm, and 76.5 mm. . . . . . . . . . 122 74 The predicted phase field contour plots of the local model with the preferential crack orientation angle equal to 45 degree at the loading displacement 100.5 mm, 103.5 mm, 106.5 mm, and 109.5 mm. . . . . . . . . . . . . . . . . . . . 123 75 Cylindrical EVA samples with three different aspect ratios (0.5, 1.0, and 2.0) between radius and height. . . . . . . . . 144 76 Comparison of the contour plots of the polymer volume fraction in Tetrahydrofuran after (a) 15 min, (b) 30 min, and (c) 45 min for the three different models with the ratio between radius and height equal to 0.5, 1.0, and 2.0. . . . . 146 77 The time history plots of swelling degree in Tetrahydrofuran for the three different models with the ratio between radius and height equal to 1.0, 0.5, and 2.0. . . . . . . . . . 147 78 Comparison of the time history curves of swelling degree in the three different solvents between experiments and numerical simulation. . . . . . . . . . . . . . . . . . . . . . 148 xvi
79 Schematic for the transient swelling of (a) pure EVA sample and (b) EVA/aluminum laminate. . . . . . . . . . . . . 149 80 Comparison of the contour plots of the polymer volume fraction at (a) 15 min, (b) 30 min, and (c) 45 min between the pure EVA samples with crosslinking degree α= 0.0 and α= 1.0. ...........................151 81 Comparison of the contour plots of the polymer volume fraction at (a) 15 min, (b) 30 min, and (c) 45 min between the EVA/aluminum laminate samples with crosslinking degree α= 0.0and α= 1.0(aluminum in grey color). . . . 152 82 Comparison of the time history curves of swelling degree between the pure EVA samples and EVA/aluminum laminates with crosslinking degree α= 0.0and α= 1.0. . . . . 154 83 The sectional sketch of the PV module in the recycling process.................................154 84 Finite element model of a quarter of the complete PV laminate. ...............................154 85 The contour plots of the polymer volume fraction ϕof the first EVA layer (EVA1) in the three different solvents Toluene, Tetrahydrofuran, and Octane after (a) 5 h, (b) 10 h,and(c)15h...........................156 86 The spatial variation of the polymer volume fraction ϕof the first EVA layer (EVA1) in the three different solvents Toluene, Tetrahydrofuran, and Octane after (a) 1 h, (b) 5 h, (c) 10 h, and (d) 15 h. . . . . . . . . . . . . . . . . . . . . . . 157 87 The contour plots of the polymer volume fraction ϕof the second EVA layer (EVA2) in the three different solvents Toluene, Tetrahydrofuran, and Octane after (a) 5 h, (b) 10 h,and(c)15h...........................158 88 The spatial variation of the polymer volume fraction ϕof the second EVA layer (EVA2) in the three different solvents Toluene, Tetrahydrofuran, and Octane after (a) 1 h, (b) 5 h, (c) 10 h, and (d) 15 h. . . . . . . . . . . . . . . . . . . . . . . 159 xvii
89 The time history evolution of the EVA swelling degree in the three different solvents Toluene, Tetrahydrofuran, and Octane...............................161 90 The maximum principal stress plots of the silicon cell layer in the three different solvents Toluene, Tetrahydrofuran, and Octane after (a) 5 h, (b) 10 h, and (c) 15 h. . . . . . . . . 162 91 The damage plots of the silicon cell layer in the three different solvents Toluene, Tetrahydrofuran, and Octane after (a) 5 h, (b) 10 h, and (c) 15 h. . . . . . . . . . . . . . . . . 164 92 A sketch of schematic and kinematic definitions of the interface between two bodies along the deformation process 168 93 A sketch of three-dimensional interface element with element nodes and integration points . . . . . . . . . . . . . . 171 94 A sketch of one single element simulation . . . . . . . . . . 175 95 Traction separation curves of one single element simulation with different ultimate displacement values . . . . . . 176 96 (a) 90-degree peeling test, and (b) discretization . . . . . . 177 97 Displacement contour in the peeling direction along the simulationprocess .......................178 98 Peeling force vs. displacement curves obtained from simulation at the condition 85 ◦C/85%RH. ...........178 99 Comparison of adhesion strength between experiment and simulation with different exposure duration at the same condition 85 ◦C/85%RH. ...................179 100 Peeling force vs. displacement curves obtained from simulation at the condition 85 ◦C..................180 101 Comparison of adhesion strength between experiment and simulation for different relative humidity values at the same condition 85 ◦C and exposure duration of 72 h. . . . . . . . 180 102 Peeling force vs. displacement curves obtained from simulation at the condition 85%RH. ...............181 103 Comparison of adhesion strength between experiment and simulation for different temperature values at the same condition 85%RH and exposure duration of 72 h. . . . . . 182 xviii
104 Comparison of adhesion strength versus humidity dose for different conditions between experimental fitting and numericalresults ........................183 105 Comparison of the required specific energy for different relative humidity and temperature conditions . . . . . . . 185 106 Relationship of average particle size with the required specific energy for crushing in recycling of PV [262] . . . . . . 186 xix
List of Tables 1 Mechanical and thermal properties of the PV materials. . . 56 2 Mechanical properties for Photovoltaic modules [87, 53]. . 115 3 Model parameters of EVA copolymer for the coupled modeling. ...............................144 4 Solvent properties for the swelling simulation [245, 247]. . 145 5 Mechanical properties for PV modules [87, 53]. . . . . . . . 155 6 Testing conditions for Photovoltaic modules. . . . . . . . . 175 xx
Acknowledgements Firstly, I am very grateful to my supervisors Prof. Marco Paggi and Prof. Jose Reinoso for their tremendous guidance and support in terms of both research and life during my PhD study. I enjoy a lot during this journey. This thesis is based on the articles listed in the section of Publication. The work in Chapter 2 is taken from Publication 1 and 2, the work in Chapter 3 is taken from Publication 3 and 6 (In Preparation), and the work in Chpater 4 is taken from Publication 4 and 5. I also would like to acknowledge the funding received from the European Union’s H2020-MSCA-ITN-2019 research and innovation program under the Marie Skłodowska-Curie grant agreement no. 861061 – Project NEWFRAC ”New strategies for multifield fracture problems across scales in heterogeneous systems for Energy, Health and Transport”. Finally, I would like to extend my heartfelt appreciation to all my friends and family members. xxi
Publications 1. Liu, Z., Reinoso, J. and Paggi, M., 2022. Hygro-thermo-mechanical modeling of thin-walled photovoltaic laminates with polymeric interfaces. Journal of the Mechanics and Physics of Solids, 169, p.105056. 2. Liu, Z., Lenarda, P., Reinoso, J. and Paggi, M., 2023. A multifield coupled thermo-chemo-mechanical theory for the reaction-diffusion modeling in photovoltaics. International Journal for Numerical Methods in Engineering, 124(12), pp.2876-2901. 3. Liu, Z., Reinoso, J. and Paggi, M., 2022. Phase field modeling of brittle fracture in large-deformation solid shells with the efficient quasi-Newton solution and global–local approach. Computer Methods in Applied Mechanics and Engineering, 399, p.115410. 4. Liu, Z., Marino, M., Reinoso, J. and Paggi, M., 2023. A continuum largedeformation theory for the coupled modeling of polymer–solvent system with application to PV recycling. International Journal of Engineering Science, 187, p.103842. 5. Liu, Z., Reinoso, J. and Paggi, M., 2022. A humidity dose-CZM formulation to simulate new end-of-life recycling methods for photovoltaic laminates. Engineering Fracture Mechanics, 259, p.108125. 6. Liu, Z., Lenarda, P., Reinoso, J. and Paggi, M, 2023. Phase field modeling of silicon crystalline cracking in 3D thin-walled photovoltaic laminates. (In Preparation) xxii
Abstract This thesis is mainly focused on the computational modeling of solar cell cracking, multiphyics phenomena, and recycling of photovoltaic (PV) modules through the finite element method. Specifically, it consists of three parts. In the first part, a comprehensive hygro-thermo-mechanical computational framework in the 3D setting is proposed to model the coupled degradation phenomena in the PV modules for the durability analysis, and it 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. The second part is focused on the crack modeling of very thin and brittle silicon solar cells in the PV modules, and a reliable computational framework integrating solid shell element formulation with phase field fracture modeling is developed using the efficient quasi-Newton solution scheme and global local approach. The excellent performance is showcased through the simulation of different boundary value problems, and then applied to predict the crack growth of silicon solar cells when the PV modules are subjected to different external loadings. The third part addresses the efficient recycling of PV modules through the numerical modeling method by the development of 3D interface finite element with humidity-dose enhanced cohesive zone model for the peeling simulation to separate different layers, and diffusion-swelling large deformation continuum theory for the nondestructive recovery of silicon cells in the PV recycling using the solvent method. With these tools at hand, it is possible to design suitable virtual testing procedures for PV durability and recyclability analysis. xxiii
Chapter 1 Introduction 1.1 Introduction to photovoltaic modules Nowadays renewable energy plays a key role in the socio-economic development, and different types of solutions through various technologies rapidly emerge due to the increasing demand of reducing the energy deficit. Photovoltaic (PV) technology is the dominant choice in the renewable energy market, and it is not only widely used for the industry, but also for individual needs. A PV module is usually a laminate composed of tempered glass, solar cells, backsheet and encapsulant layers, see Fig. 1. Different components are adhesively bonded to each other. This multilayer system can ensure the safety and performance of PV modules but forms a number of interfaces that are potential paths for contaminant ingress, causing delamination [2], which eventually promotes many other further ageing mechanisms such as polymer degradation, corrosion of metalic materials inside the system, loss of light transmission, and so on [3, 4]. The reliability of encapsulation is therefore crucial for the durability and the service life of PV modules [5, 6]. The most commonly used encapsulant material in the existing PV installations is Epoxy Vynil Acetate (EVA) that, interposed in the form of layers between the glass and the silicon solar cells, and between the solar cells and the backsheet, holds the components together 1
Frame Glass EVA Solar Cell EVA Backsheet Figure 1: The general structure of a typical PV module [1]. and protects the module from outside harsh environments [7]. The conversion of solar energy into thermal energy for the various industrial processes not only reduces the dependency on fossil resources but only minimizes the harmful emissions [9]. Photovoltaic technologies are widely used in the industrial processes, such as steam demand processes, industrial space heating, buildings and so on [10]. The majority of solar energy applications lies in the solar water heating industry, which is regarded as the most effective alternative among all the available solar thermal technologies. Solar water heating system usually consists of solar collectors and storage space, which works on the basis of the density difference between hot and cold water. Another application of solar energy is the refrigeration and air-conditioning system, and generally there are two different solar air-conditioning systems including the closed and open systems. Closed air-conditioning solar systems are capable of being installed with high temperature collectors, while open air-conditioning solar systems can heat the ambient air to certain temperature and regenerate the sorbent. Fig. 2 shows the rooftop airconditioning solar system [8]. Due to the increasing energy demand in agriculture, photovoltaic technologies have also been widely applied to 2
transportation and installation of PV modules in the field, as well as during the operation when subjected to the harsh environmental conditions such as wind, snow loading, hail impacts and so on [40, 41]. The cracks of silicon solar cells are usually invisible, but could produce electrically disconnected regions that significantly increase the electrical resistance, and hence reduce the power output of PV modules. Overall, the length, width and orientation of cracks in the silicon wafers directly influence the output of solar panels as pointed out in [42, 43, 44]. According to the survey of over 200 PV modules, more than 20 % of power loss has been detected due to the cell cracks in combination with EVA delamination and degradation [45]. Experimental studies in [46] reported that approximately 4 % power degradation are caused by cracks of solar cells in the mechanical load tests. When the disconnected cell area caused by crack formation and propagation is greater than 8 % of the total cell layer, the associated power loss was found to be roughly proportional to the electrically inactive area [28]. Regarding the crack direction in solar cells, several types of cracks can be observed from previous experimental studies, including parallel and perpendicular cracks, +45 and -45 cracks, and multiple direction cracks [47], as shown in Fig. 7. Different cracks lead to different power degradation of PV modules. Therefore, it is of utmost importance to develop reliable simulation tools for the fracture modeling of brittle thin-film solar cells in order to understand the impact of cracks on the performance of PV modules. Microcracks of silicon solar cells can be formed at different stages, including the cutting of an ingot, the lamination of PV modules, the transportation and installation, and the environmental loadings during the service lifetime. The silicon substrates for the manufacturing of PV modules are typically produced by casting the silicon ingots and subsequent wafer sawing, and the wire sawing of silicon ingots can result into the cracking of wafers [49]. According to the previous study [48], the risks of crack formation increase as the thickness of the silicon solar cells has been drastically decreased. Besides, due to the difference of thermal expansion properties, thermo-mechanical residual stress can be generated during the soldering process, leading to the silicon solar cell crackings. 9
Figure 7: Different crack patterns in silicon solar cells [48]. On the other hand, the residual stress can also be generated during the lamination process of PV modules, which increases the propensity of crackings. Besides the factors during production, cracks of silicon solar cells can be produced due to improper transportation as reported in [50]. Currently there is no standardard on the way of PV transportation, and it is pointed out in [40] that transporting the PV products vertically can reduce the risks of breaking silicon solar cells compared to the way of transporting the modules horizontally. Finally, cracks of silicon solar cells and power degradation may be caused by the environmental loadings during the service lifetime, such as wind [51], snow loading [52], and hail impact [53]. Such environmental loadings increase the stress beyond the residual stress from the production processes, and lead to higher crack probability of silicon solar cells [37]. Although cracks of silicon solar cells are unavoidable as mentioned above, it is of significant importance to identify the crack characteristics so as to improve the durability of the photovoltaic modules. 10
1.4 Sustainability of photovoltaic modules: Recycling of the end-of-life PV modules Massive PV modules have been deployed in the recent decades, and given the service lifetime of approximately 30 years, the pressure to handle the end-of-life modules starts to emerge. The PV modules that approach the end of the service lifetime can be reused or repaired so as to lower the environmental impacts without additional energy and material investment. Although numerous economic and regulatory challenges must be addressed, the reusing strategies of PV modules are worth investigating from the circular-economy point of view. The key is that the burden from the repair or reuse of end-of-life PV modules should be lower than the environmental benefit from the subsequent reuse. Note that the reuse of PV modules requires the necessary repair for safety and performance. Due to the shorter service lifetime and lower power efficiency of reused PV modules, the system owners may not gain much economic values from the second-hand products [54]. However, to avoid the material loss of end-of-life modules, the landfills have been discouraged by the European policy since the beginning of this century. The values of end-of-life PV materials are estimated to be approximately 15 billion US dollars and 2 billion new PV products can be manufactured by the reuse of the reclaimed materials, which accounts for around 630 GW of new capacity [55]. Furthermore, there is a strong incentive to avoid environmental pollution for the recycling of end-of-life PV modules [56, 57]. In general, the environmental requirement and material recovery necessitate the development of reliable PV recycling techniques for the circular source of PV materials. The improper disposal of PV panels by landfilling is intolerable nowadays as the hazardous metal components such as lead and cadmium, among others, would cause pollution to the environment [58]. Besides, it is of great importance to recover the precious silicon wafers and other metals for the manufacturing of new products to reduce the cost, supporting the sustainable development of energy production [59, 60, 61, 62]. Normally, the PV module is a laminate structure composing of dif- 11
ferent plies including tempered glass, silicon solar cells, and backsheet, which are bonded to each other by ethylene-co-vinyl acetate (EVA) layers [20]. Therefore, many approaches have been developed in the relatd literature in order to separate the different layers by the removal of EVA and then reclaim valuable materials in the end-of-life PV modules [1, 63, 64, 65]. The use of nitric acid to remove EVA layers was proposed in [66], but the resulting NOxemissions are harmful to the environments. Thermal decomposition is another alternative method to remove EVA layers but leads to the release of harmful gas because of the existence of fluorinated compounds [67, 68, 69]. Also, some thermal or mechanical pretreatments are designed to dissemble the backsheet before the further recycling process [70, 71, 72]. In addition, recycling of PV modules using the organic solvent method has received a great deal of attention in the recent years [73, 74, 75, 76, 77]. In spite of these efforts, it is very challenging to reclaim unbroken silicon wafers, which is a key problem in the PV recyling as the thickness of silicon cell layers was reduced to approximately 0.1 mm after 2006 [78]. The major factor responsible for the degradation of PV modules is the aging of EVA layers, while the service lifetime of silicon wafers is much longer than that of PV modules [73]. In this regard, much cost can be saved with the reclaim of structurally intact silicon wafers for the manufacturing of new PV products, since the silicon ingot manufacturing process can be omitted. High environment pollution and energy consumption can be avoided with such an ideal reusing strategy that also reduces around 40%cost of PV production [79]. Unfortunately, very scarce research activities have been devoted to the nondestructive reclaim of intact silicon wafers in PV recycling, such as thermal decomposition at specific heating rates [80] and breakage of glass layer before thermal treatment [81]. Although it is possible to recover the intact silicon wafers using the these methods, the thermal decomposition is actually difficult to control and the cracking of glass during the treatment tends to break the the silicon wafers. In addition to this, results in [73] have indicated that the swelling of EVA during the recycling process using the solvent method can lead to the cracking of silicon cell layer. Therefore, there is 12
an urgent need to develop a suitable computational framework to help control the recycling process for the recovery of structurally intact silicon wafers. In the recent decade, recycling of photovoltaic modules have received a lot attention, and different organizations are established to address this topic. PHOTORAMA is an EU funded innovation action striving to improve recycling of solar panels, and it is a consortium consisting of 13 organizations from 2021 to 2024. Another interesting project is the CABRISS, which is also an EU funded consortium with 11 companies and 5 research institutes for pioneering a circular economy on photovoltaics. A France-based start-up named ROSI has claimed to develop an innovative and viable process to recover silicon from the end-of-life photovoltaic modules. Besides, there is another European project ReProSolar led by Veolia Germany that designs a new delamination technology which is able to separate solar cells efficiently from the glass plate for the recovery of all materials without shredding the photovoltaic modules. 1.5 Outline of the thesis This work is mainly focused on the computational modeling of thinwalled photovoltaic laminates with application to the durability and recycling analysis. In this section, each chapter in the following is outlined. Chapter 2 presents the comprehensive computational framework at finite deformation for the multiphysics modeling of the degrdation phenomena in the photovoltaic modules. The developed framework is then applied to simulate the three standard international tests of PV modules, including the damp heat test, the humidity freeze test, and the thermal cycling test. Chapter 3 addresses the crack modeling of brittle silicon solar cells in the photovoltaic modules using the global-local approach with solid shell element formulation. For the efficient modeling, the quasi-Newton monolithic solution scheme is proposed, and its excellent performance is demonstrated by the comparison with the standard staggered Newton solution scheme through the different boundary value problems. This 13
computational framework is then applied to the modeling of silicon solar cell crackings in the PV modules under different boundary conditions, which showcases its capacity to predict the crack growth in the thin-film silicon solar cells. Chapter 4 outlines the numerical simulation of PV recycling problems. It consists of two parts. In the first part, a large deformation theory is proposed to model the swelling-diffusion problem in the PV recycling using the solvent method, which aims at predicting the EVA swelling induced cracking of silicon solar cells. In the second part, given the degradation of EVA in the end-of-life PV modules, a computational framework integrating large deformation interface element with the polynomial cohesive zone model incorporating the moisture and temperature effects is developed to simulate the peeling tests of PV panels. Chapter 5 addresses the main conclusions and future developments of this thesis. 14
Chapter 2 Multiphysics modeling for the durability analysis of photovoltaic laminates The durability analysis of PV modules has received much attention in the past decades [82, 83, 20, 84], and modeling polymeric materials EVA is quite challenging as it requires a multiphysics framework to predict their overall performance and degradation accurately [85]. This chapter establishes a comprehensive computational framework for the coupled modeling of hygro-thermo-chemo-mechanical degradation phenomena of PV laminates for durability analysis. The main content of this chapter is taken from [20, 84]. 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 [86, 87]. 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 proce- 15
dure of parameter identification, the fraction calculus method has been proved to be very effective for the modeling of visco-elastic constitutive behaviour [88, 89, 90]. Thus this formulation has been adopted to describe the rheologically complex behaviour of polymeric EVA, whose microstructure changes with temperature, and according to [91], 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 properly. Analytical solution for the moisture diffusion along EVA layer proposed in [30] 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 damp heat test, its validity in the cases of cyclic temperature boundary loading in the humidity freeze and 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. 8. 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-Rapson scheme in a standard way. 16
Nonlinear heat conduction Nonlinear thermoelasticity Thermo-mechanical coupling: •Thermo-mechanical shell •Thermo-visco-elastic CZM Nonlinear moisture diffusion Monolithic fully implicit solution scheme for thermomechanical problem Implicit Newton-Rapson solution scheme for moisture diffusion problem Figure 8: The proposed solution scheme for the 3D hygro-thermo- mechanical modeling framework. To save the computational efficiency, the EVA layers are modelled using zero-thickness interfaces elements with thermo-visco-elastic cohesive zone model based on the fractional calculus method, and a thermomechanical 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 [92, 93] 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 fol- 17
lowed: (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 [94, 95, 96, 97, 98, 99], and (ii) shell formulations relying on the solid shell concept through the parametrization of the top and bottom surfaces of the body [100, 101, 102, 103, 104]. 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 [105, 106, 107, 108]. 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 [109, 110, 111, 112], the ANS method [113, 114], and combination of them [115, 116, 117]. This section is structured as follows. In Section 2.1, 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 is 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 2.1.4. 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, and this section also covers the 3D thermo-mechanical 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 2.1.5. Then this modeling methodology is applied to simulate the three international standard tests of PV laminates, 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 the constant temperature boundary condition, as well as experimental electroluminescence images obtained from the thermal cycling test with 18
𝜹𝒔(𝒕) 𝒄 −𝜹𝒔(𝒕) 𝒄 ∆𝒔(𝒕) 𝝉𝒔(𝒕) 𝒄 −𝝉𝒔(𝒕) 𝒄 𝟏 𝑲𝒔(𝒕) 𝝉𝒔(𝒕) 𝜹𝒏 𝒄∆𝒏 𝝈𝒏 𝒄 𝟏 𝑲𝒏 𝝈𝒏 𝝐 𝟏 (a) (b) 𝜹𝒏 −𝝈𝒏 𝑲𝒔(𝒕) 𝟏 Figure 10: The traction versus separation curves of the polymeric interface: (a) sliding and shearing modes, and (b) the opening mode. To obtain the structural response of encapsulant layers, the stiffness Knis related to the actual stiffness of polymeric materials in the normal direction, which is calculated as the ration between the Young’s modulus Epoly and its thickness hpoly, i.e., Kn=Epoly/hpoly. Analogously, the tearing and sliding stiffness can be expressed as: Ks=Kt= Epoly/[2hpoly(1 + νpoly)]. Since polymeric materials have thermo-visco- elastic constitutive behaviors, the Young’s modulus Epoly depends both on the average temperature ⟨T⟩and time history t. To synthetically characterize these dependencies, a fractional calculus approach proposed in [89] is adopted here instead of the use of Prony series representation, which has been proved to be very effective for parameters identification [91, 122]. Accordingly, the modulus Epoly is given by Epoly(t, T) = a(T)h(t, T )−α(T) Γ(1 −α(T)) (2.20) where aand αare two temperature dependent functions, such that 0< a, α < 1, and Γ(t)is the Euler gamma function Γ(t) = Z∞ 0 e−xxt−1dx(2.21) 25
The function h(t, T)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 h(t, T)is equal to the difference between the current time tand t0corresponding 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, q1=q2= 0 and q=q3is expressed as q=(−h01−∆n δc n⟨T⟩if ∆n∈Jn 0if ∆n/∈Jn (2.22) where h0is the thermal conductivity of interface without crack opening, i.e., ∆n= 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 [123]. 2.1.3 Moisture diffusion along 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 c(x1, x2, x3, t)inside the polymeric layer for each material point and time. The boundary value problem for moisture diffusion, in which an imposed moisture content c∗is applied to the boundary, can be described as follows ∂c ∂t (x1, x2, x3, t)−D∇2c(x1, x2, x3, t)=0 in Bpoly ×[0, tf] c(x1, x2, x3,0) = 0 in Bpoly c(x1, x2, x3, t) = c∗in ∂Bpoly ×(0, tf] (2.23) 26
where Dis the moisture diffusion coefficient. It should be pointed out that the moisture diffusion is characterized with a different time scale from that of thermo-mechanical problem. The characteristic velocity of moisture diffusion is related to the diffusion coefficient D, while that of temperature diffusion is governed by the ratio k0/ρcp. Given the characteristic values for polymeric layers, the ratio between these two physical phenomena is [k0/(ρcp)] /D ≈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 [30], the moisture diffusion coefficient should be considered as temperature and gap dependent Arrhenyus type function D=(Aexp −Ea T R ,if ∆n≤δc n Aexp −Ea T R ∆n δc n,if ∆n> δc n (2.24) For the purpose of accounting for the debonding effect of the encapsulant layer on the moisture diffusion, Dis assumed to be a linear increasing function of the interfacial gap ∆nwhen it overcomes the critical value δc n. 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 Dis determined by the average temperature and gap displacement at the previous time increment computed from the coupled thermo-mechanical problem. 2.1.4 Variational form and finite element approximation 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. Firstly, this part 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, the consistent linearization of 27
the governing equations for the thermo-mechanical cohesive interface is presented. Finally, the finite element implementation of the 3D moisture diffusion along the encapsulant polymeric layers is detailed. In the following, the variational basis and finite element interpolation of the initial boundary value problem with the coupled thermomechanical 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 RT(u,˜ E, T, δT) = ZB0 ρ0cp˙ TδT dΩ−ZB0 TZ:˙ EδT dΩ +ZB0∇·QδT dΩ=0, (2.25) where δT is the virtual temperature field. Invoking the Gauss’s theorem, the third term relevant to the divergence of heat flux can be reformulated as ZB0∇·QδT dΩ = Z∂B0,¯q QNδT d∂Ω−ZB0 Q∇XδT dΩ, (2.26) where QN=Q·Ndenotes the Neumann boundary condition on ∂B0,¯q. Subsequently, Eq. (2.25) can be rewritten as RT(u,˜ E, T, δT) = ZB0 ρ0cp˙ TδT dΩ−ZB0 TZ:˙ EδT dΩ +Z∂B0,¯q QNδT d∂Ω−ZB0 Q∇XδT dΩ=0. (2.27) By inserting the Duhamel’s law, i.e., simplified Eq. (2.17) under isotropic 28
assumption, the previous expression can be formulated as RT=ZB0 ρ0cp˙ TδT dΩ−ZB0 TZ:˙ EδT dΩ +Z∂B0,¯q QNδT d∂Ω + ZB0 [∇XδT ]Tk0G∇XTdΩ=0, (2.28) with the internal residual RT int and external counterpart RT ext being identified as RT int =ZB0 ρ0cp˙ TδT dΩ−ZB0 TZ:˙ EδT dΩ +ZB0 [∇XδT ]Tk0G∇XTdΩ=0, (2.29) RT ext =Z∂B0,¯q QNδT d∂Ω.(2.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 [116, 124, 125] 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 u, the incompatible strain ˜ E, and stress Sare 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 Euand an enhanced incompatible part ˜ E[107, 108], i.e., E=Eu+˜ 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 Ru(u,˜ E, T, δu) = ZB0 S:δEudΩ−ZB0 ρ0¯ γδudΩ−Z∂B0,¯ t ˆ tδud∂Ω=0, (2.31) 29
R˜ E(u,˜ E, T, δ ˜ E) = ZB0 S:δ˜ EdΩ = R˜ E int = 0,(2.32) where R˜ E int represents the internal contribution from the enhanced strain field to the functional. Here the internal residual Ru int and external residual Ru ext corresponding to the kinematic compatible part can be defined as Ru int =ZB0 S:δEudΩ,(2.33) Ru ext =ZB0 ρ0¯ γ·δudΩ−Z∂B0,¯ t ˆ t·δud∂Ω.(2.34) The initial configuration B0is discretized into nenon-overlapping finite elements, i.e., B0≈Sne e=1 B(e) 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 X(ξ) = 1 21 + ξ3Xt(ξ1, ξ2) + 1 21−ξ3Xb(ξ1, ξ2),(2.35a) x(ξ) = 1 21 + ξ3xt(ξ1, ξ2) + 1 21−ξ3xb(ξ1, ξ2),(2.35b) where the subscripts tand bstand for the top and bottom surfaces, respectively. The parametric space is given by: A:= {ξ= (ξ1, ξ2, ξ3)∈ R3|−1≤ξi≤+1; i= 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 Xand xare interpolated through standard trilinear shape functions NIas follows X≈ nn X I=1 NI(ξ)XI=Ne X,x≈ nn X I=1 NI(ξ)xI=Nex,(2.36) 30
where XIand xIare the nodal position vectors in the reference and current configurations, respectively, with number of nodes nn= 8, and e X and exare the respective nodal position vectors in the element level. The shape functions NIis defined as NI=diag [NI, NI, NI](2.37) where NI=1 81 + ξ1 Iξ11 + ξ2 Iξ21 + ξ3 Iξ3,I= 1,2, ..., 8. The interpolation of displacement and temperature fields (u, T), the variations (δu, δT ), and the increments (∆u,∆T) are given by u≈ nn X I=1 NI(ξ)dI=Nd, δu≈ nn X I=1 NI(ξ)δdI=Nδd, ∆u≈ nn X I=1 NI(ξ)∆dI=N∆d, (2.38) T≈ nn X I=1 NI(ξ)TI=ˆ Nˆ T, δT ≈ nn X I=1 NI(ξ)δTI=ˆ Nδˆ T, ∆T≈ nn X I=1 NI(ξ)∆TI=ˆ N∆ˆ T, (2.39) where dIand TIrepresent the nodal displacement and temperature, respectively, and dand ˆ Tare the corresponding vectors at the element level. The interpolation of the spatial temperature gradient in the reference configuration ∇XT, its variation ∇XδT and its increment ∇X∆T reads ∇XT≈G−T∇ξˆ Nˆ T, ∇XδT ≈G−T∇ξˆ Nδˆ T, ∇X∆T≈G−T∇ξˆ N∆ˆ T, (2.40) where ∇ξrepresents the gradient with respect to the natural coordinates in the curvilinear setting. 31
Midsurface Element node Collocation point Midsurface A1 A2 B2 C1 B1 C2 C3 C4 ANS for transverse shear locking ANS for trapezoidal locking 𝒆𝟑𝒆𝟐 𝒆𝟏 Figure 11: The position of collocation points in the parametric space for the ANS method. The Green-Lagrange strain vector is expressed as Eu=[E11, 2E12, 2E13, E22, 2E23, E33]T. To overcome the curvature thickness locking, the ANS method proposed in [113] is adopted to modify the strain component E33, which requires four collocation points defined in convective coordinates ξCias ξC1= (−1,−1,0),ξC2= (1,−1,0),ξC3= (1,1,0), and ξC4= (−1,1,0), see Fig. 11. Besides, to prevent transverse shear locking, the ANS method proposed in [114] 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. 11. Accounting for the ANS interpolations, the strain vector is given by Eu= 1 2(g11 −G11) (g12 −G12) 1−ξ2gA1 13 −GA1 13 +1 + ξ2gA2 13 −GA2 13 1 2(g22 −G22) 1−ξ1gB1 23 −GB1 23+1 + ξ1gB2 23 −GB2 23 P4 i=1 1 41 + ξ1 iξ11 + ξ2 iξ21 2gCi 33 −GCi 33 (2.41) The approximation of strain variation and increment is interpolated by δEu≈ nn X I=1 BI(ξ)δdI=Bδd,∆Eu≈ nn X I=1 BI(ξ)∆dI=B∆d(2.42) 32
where BIis defined as BI= NI,1gT 1 NI,1gT 2+NI,2gT 1 1−ξ2 NA1 I,1gA1 3T+NA1 I,3gA1 1T!+1 + ξ2 NA2 I,1gA2 3T+NA2 I,3gA2 1T! NI,2gT 2 1−ξ1 NB1 I,2gB1 3T+NB1 I,3gB1 2T!+1 + ξ1 NB2 I,2gB2 3T+NB2 I,3gB2 2T! P4 i=1 1 41 + ξ1 iξ11 + ξ2 iξ2NI,3gi 3T (2.43) Similarly, the interpolation of enhanced incompatible strain field ˜ E, its variation δ˜ Eand its increment ∆˜ Etake the form of ˜ E≈M(ξ)ς, δ ˜ E≈M(ξ)δς,∆˜ E≈M(ξ)∆ς, (2.44) where ςis the enhancing modes vector to prevent locking pathologies as pointed out in [124, 93], and Mdenotes the interpolation matrix of incompatible strain [99], which is given by M(ξ) = det J0 det JT−T 0˜ M(ξ)(2.45) where J= [G1,G2,G3]T,J0is its evaluation at the element center (ξ1= 0,ξ2= 0,ξ3= 0), and the transformation matrix T0takes the form T0= J2 110J2 210J2 3102J110J2102J110J3102J210J310 J2 120J2 220J2 3202J120J2202J120J3202J220J320 J2 130J2 230J2 3302J130J2302J130J3302J230J330 J110J120J210J220J310J320J110J220+J210J120J110J320+J310J120J210J320+J310J220 J110J130J210J230J310J330J110J230+J210J130J110J330+J310J130J210J330+J310J230 J120J130J220J230J320J330J120J230+J22J130J120J330+J320J130J220J330+J320J230 (2.46) where JIJ0in the transformation matrix T0are the components of J0. To alleviate the membrane, volumetric and Poisson thickness locking effects, the interpolation matrix ˜ M(ξ)of the enhancing modes defined in the parametric space ξ={ξ1, ξ2, ξ3}reads ˜ M= ξ1000ξ1ξ20 0 0 0 0 0 0 0 ξ1ξ20 0 ξ1ξ20 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0ξ20 0 0 ξ1ξ20 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 ξ1ξ2ξ1ξ2 .(2.47) 33
By inserting the interpolation formulae, the discrete forms of internal energy balance residual Eq. (2.29) is given by ˆ RT int(d,ς,ˆ T, δ ˆ T) =δˆ TT"ZB0 ˆ NTρ0cp˙ TdΩ−ZB0 ˆ NTZT˙ ETdΩ +ZB0 k0BT TG∇XTdΩ#, (2.48) where BTis the interpolation operator for the temperature gradient. The discrete internal residual terms with regard to the kinematic compatible part, Eq. (2.33), and enhanced strain field, Eq. (2.32), can be expressed as ˆ Ru int(d,ς,ˆ T, δd) = δdTZB0 BTSdΩ,(2.49) ˆ R˜ E int(d,ς,ˆ T, δς) = δςTZB0 MTSdΩ.(2.50) For subsequent developments, the residual vectors associated with the displacement field, the incompatible strain field, and the temperature field, respectively, are defined as RT int =ZB0 ˆ NTρ0cp˙ TdΩ−ZB0 ˆ NTZT˙ ETdΩ + ZB0 k0BT TG∇XTdΩ, (2.51a) Rd int =ZB0 BTSdΩ,(2.51b) Rς int =ZB0 MTSdΩ.(2.51c) To solve the set of nonlinear residual equations, the iterative scheme is adopted for the multi-field coupled thermo-mechanical problem, and the consistent linearization of the system derived from the concept of Gateaux directional derivative in matrix form can be expressed as kdd kdς kdT kςd kςς kςT kT d kT ς kT T ∆d ∆ς ∆ˆ T = Rd ext 0 RT ext − Rd int Rς int RT int (2.52) 34
where Nc= [N1, N2, N3, N4, N5, N6, N7, N8]is the shape function matrix, and ˜ cis the moisture concentration nodal vector. Besides, the interpolation of the spatial concentration gradient ∇Xcand its variation ∇Xδc can be expressed as ∇Xc≈G−T∇ξNc˜ c=Bc˜ c,∇Xδc ≈G−T∇ξNcδ˜ c=Bcδ˜ c,(2.77) Introducing Eq. (2.76) and Eq. (2.77) into Eq. (2.75), the variational form of moisture diffusion is given by δ˜ cTZB0DBT c∇Xc+NT c˙cdΩ = δ˜ cTfc int (2.78) where fc int is the internal residual vector, which is defined as fc int =ZB0DBT c∇Xc+NT c˙cdΩ.(2.79) The linearization for moisture diffusion at the (k+ 1)th iteration in the Newton-Rapson iterative solution scheme reads Kc,k∆˜ ck+1 =−fc,k int (2.80a) ˜ ck+1 =˜ ck+ ∆˜ ck+1 (2.80b) where Kcis the element stiffness for moisture diffusion, which is given by Kc=ZB0DBT cBc+1 ∆tNT cNcdΩ (2.81) 2.1.5 Staggered solution scheme for the multi-field framework In this section, the staggered computational procedure for the solution of thermo-mechanical problem with solid shell element and interface element formulation and 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 htn, t(k) n+1i, 41
where tnand t(k) n+1 stand for the previous converged increment and prospective current increment at iteration k, respectively. Given the data {dn,ςn,ˆ Tn} at the previous converged increment, the nonlinear incremental solution requires the Newton iterations at the intermediate state {d(k) n+1,ς(k) n+1,ˆ T(k) n+1}. Note that {dn,ςn,ˆ Tn}and {d(k) n+1,ς(k) n+1,ˆ T(k) n+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)th iteration, respectively. Based on the static condensation described in Section 2.1.4, the nodal displacement and temperature vectors are defined as unknowns in the element level, and the increment of enhancing vector in the (k)th iteration of next time increment, ∆ς(k) n+1, should be determined. According to the procedure proposed in [101], it is calculated as ∆ζ(k) n+1 is given by ∆ς(k) n+1 =−[kςς,n]−1hRς int,n +kςd,n∆d(k) n+1 +kςT,n∆ˆ T(k) n+1i(2.82) It should be pointed out that the increments ∆d(k) n+1 and ∆ˆ T(k) n+1 are provided by the solver, while the element matrices at the previous increment [kςς,n]−1,Rς int,n,kςd,n,kςT,n, and ςnare all stored as internal variables. The numerical algorithm for the finite element implementation of 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 hat the integration points need to be determined due to the time dependency of visco-elastic constitutive equation, Eq. (2.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 thermomechanical problem are solved in a given time increment, these nodal 42
Algorithm 1: Numerical implementation procedure of the coupled thermo-mechanical formulation if Solid shell element then Data: dn,ˆ Tn,∆d(k) n+1,∆ˆ T(k) n+1 Result: dn+1,ˆ Tn+1 Initialization of ςn,Rς int,n,kςd,n,[kςς,n]−1; while ||R∗ d|| >tolerance do Compute ∆ς(k) n+1 =−[kςς,n]−1hRς int,n +kςd,n∆d(k) n+1 +kςT,n∆ˆ T(k) n+1i; Update the enhancing vector ς(k) n+1 =ςn+ ∆ς(k) n+1; for n←1to 8 integration points do Compute the curvilinear basis G(k) n+1 and g(k) n+1; Compute the B matrices B(k) n+1 and B(k) T,n+1; Modify the B(k) n+1 matrix according to the ANS method; Compute C(k) n+1,S(k) n+1, and Z(k) n+1; Compute the EAS operator M(k) n+1; end Compute the element stiffness matrices k(k) dd,n+1,k(k) dς,n+1,k(k) ςd,n+1,k(k) ςς,n+1, and k(k) T T,n+1; Compute the internal force vectors Rd(k) int,n+1,Rς(k) int,n+1, and RT(k) int,n+1; Perform the static condensation and final assembly; end end else if Cohesive interface element then Data: dn,ˆ Tn,∆d(k) n+1,∆ˆ T(k) n+1 Result: dn+1,ˆ Tn+1 while ||R∗ d|| >tolerance do Compute the operators Lcz and Mcz; for n←1to 4 integration points do Compute the shape function N(k) cz,n+1 ; Compute the matrix B(k) cz,n+1 =N(k) cz,n+1Lcz; Compute the rotation matrix in the curvilinear system R(k) cz,n+1; Compute the local gap g(k) loc,n+1 =R(k) cz,n+1B(k) cz,n+1 ˆ d(k) n+1; Compute the temperature value at the integration point T(k) int,n+1; if |T(k) int,n+1 −Tint,n|>tolerance then Compute the history variable h(k) int,n+1 =hint,n; else Compute the history variable h(k) int,n+1 =hint,n + ∆t; end Compute the stiffness matrix C(k) cz,n+1 and the traction vector T(k) n+1; end Compute the element stiffness matrix Ke(k) int,n+1; Compute the internal residual vector f(k) int,n+1; end end 43
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.1.3, given the time scale difference between the moisture diffusion and thermo-mechanical phenomena, the coefficient in the current increment of moisture diffusion 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 2: Numerical implementation procedure of the moisture diffusion along the interface Data: ˜ cn,∆˜ c(k) n+1 Result: ˜ cn+1 Initialization of dnand ˆ Tn; while ||fc int|| >tolerance do for n←1to 8 integration points do Compute the shape function N(k) c,n+1; Compute the B matrix B(k) c,n+1; Compute the temperature value at the integration point Tint,n with ˆ Tn; Compute the normal gap at the integration point ∆int,n with dn; if ∆int,n ≤δc nthen Compute the coefficient D(k) int,n+1 =Aexp −Ea Tint,nR; else Compute the coefficient D(k) int,n+1 =Aexp −Ea Tint,nR∆int,n δc n; end end Compute the element stiffness matrices K(k) c,n+1; Compute the internal force vectors fc(k) int,n+1; Perform the final assembly; end 2.2 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 en- 44
Glass EVA Glass H2O (O2) H2O (O2) Figure 14: The sketch of double-glass structure in the damp heat test. vironments, which uncouples the moisture diffusion from the thermomechanical 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 to the experimental electroluminescence images. 2.2.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. Since temperature is kept constant in this test, the moisture diffusion problem can be solved independently from the coupled thermo-mechanical problem with the constant diffusivity. The specimen is a double-glass PV laminate of span l= 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. 14. The analytical solution was obtained in [128], which will be treated as a benchmark of the computational scheme proposed above. Assuming the initial laminate is dry, the spatial and temporal solution c can be obtained as c=c∗+4c∗ π ∞ X m=0 1 2m+ 1 sin (2m+ 1)πx le(−D(2m+1)2π2t/l2)(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. 15, and as clearly seen, water 45
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. 15(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. 16) is compared with that obtained from the analytical solution (the black dotted line Fig. 16), 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 [30], 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. A series of damp heat tests were also examined with the double glass modules exposed to four different hot and humid conditions, namely 85 ◦C/85%RH, 65 ◦C/85%RH, 45 ◦C/85%RH, and 85 ◦C/40%RH, which correspond to the absolute humidity of 2.74 kg/cm3, 1.25 kg/cm3, 0.508 kg/cm3, and 1.29 kg/cm3, respectively [129]. The contour plots of moisture diffusion in the 85 ◦C/85%RH damp heat test obtained from simulation are shown in Fig. 17 and compared with the corresponding experimental EL images at different time points [129]. As can be seen from both predictions and EL images, moisture diffuses from the free edges of the EVA layer towards the central area during aging, and they agree with each other very well. The captured EL images show the dark bands at the side area of EVA, which are not detected at the beginning (Fig. 17(a)), and get thicker after 43 h (Fig. 17(b)) and 81 h (Fig. 17(c)). This trend that indicates the moisture ingress into the module over time [129] is well reproduced by the modeling method proposed in this work. To see the temperature effect on the diffusion phenomena, the contour plots of moisture diffusion and experimental EL images in the 65 ◦C/85%RH and 45 ◦C/85%RH damp heat tests are also compared and shown in Fig. 18 and Fig. 19, respectively. Due to the temperature-dependent diffusion properties, the absolute humidity in these two cases differs from that in the 85 ◦C/85%RH case. The dif- 46
(a) 100 h (g) 800 h (d) 500 h (b) 300 h (c) 400 h (e) 600 h (f) 700 h (h) 900 h (i) 1000 h 0.0000 g/cm30.0056 g/cm3 Figure 15: The contour plot of moisture diffusion along the EVA layer inside the double-glass laminate during the damp heat simulation. 47
0 4 8 1 2 1 6 2 0 0.000 0.001 0.002 0.003 0.004 0.005 0.006 M o i s t u r e c o n c e n t r a t i o n ( g / c m 3) D is t a n c e f r o m e d g e ( m m ) D a m p H e a t [ Kem pe, 2006] S im u l a t i o n _ D a m p H e a t M i a m i _ 1 y e a r [ Kem pe, 2006] S im u l a t i o n _ M i a m i _ 1 y e a r M i a m i _ 2 y e a r s [ Kem pe, 2006] S im u l a t i o n _ M i a m i _ 2 y e a r s M i a m i _ 4 y e a r s [ Kem pe, 2006] S im u l a t i o n _ M i a m i _ 4 y e a r s Figure 16: The comparison of moisture concentration vs. distance from the edge of laminate curves between diffusion simulation and damp heat test after 1000 h as well as experimental data in Miami. fusion phenomena are much slower in low-temperature conditions since the diffusion coefficients are smaller. It is clearly observed that the process for dark bands in the 65 ◦C/85%RH case (see Fig. 18) or 45 ◦C/85% RH case (see Fig. 19) to become thick is slower than that in the 85 ◦C/85% RH damp heat case. Besides, it is worth noting that the absolute humidity in the 65 ◦C/85%RH (see Fig. 18) and 85 ◦C/40%RH (see Fig. 20) are almost equal to each other (1.25 kg/cm3vs. 1.29 kg/cm3). However, due to the temperature difference, the diffusion phenomena are still different in these two cases. The predicted time history evolution of moisture concentration at the positions with different distances away from the edge in the four different cases is shown in Fig. 21. The moisture concentration at different positions increases over time as the diffusion phenomena proceed during the damp heat tests. Although the absolute humidity at the 65 ◦C/85%RH and 85 ◦C/40%RH cases are almost equal to each other, the diffusion time history demonstrates the significant difference, especially 48
(a) 0 h (b) 43 h (c) 81 h 2.74 kg/cm3 0.00 kg/cm3 Figure 17: Comparison between the moisture diffusion contour plots obtained from simulation and the EL images taken from [129] in the 85 ◦C/85%RH damp heat test after 0 h, 43 h, and 81 h. 49
(a) 0 h (b) 130 h (c) 279 h 1.25 kg/cm3 0.00 kg/cm3 Figure 18: Comparison between the moisture diffusion contour plots obtained from simulation and the EL images taken from [129] in the 65 ◦C/85%RH damp heat test after 0 h, 130 h, and 279 h. 50
100101102103104 10- 1 100 101 102 103 R e l a x a t i o n m o d u l u s ( P a ) T i m e ( s ) T = -3 5 ° C T = -2 8 ° C T = -1 8 ° C T = 0 ° C T = 20 ° C T = 40 ° C T = 49 ° C T = 60 ° C T = 80 ° C T = 119 ° C T = 139 ° C Figure 24: Relaxation modulus of EVA vs. time curves at different temperatures (adapted from [91]). hesive properties and diffusivity of EVA need to be updated during simulation, and to be specific, the viscoelastic parameters α(T)and a(T)for the calibration of Young’s modulus of EVA are temperature dependent as experimentally evaluated in [87] and interpreted through the fractional calculus method [123], see the plot of modulus Ein Fig. 24. Regarding the properties of moisture diffusion, the coefficient Dis determined according to the Arrhenius type equation Eq. (2.24) reported in [30]. The thermal and mechanical properties of different PV materials are taken from [86, 87], which are listed in Table 1. Besides, the critical opening of polymeric interface δc ncan be estimated from the experimental data of variation of Mode 1 fracture energy with respect to temperature in [32], 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. 25. 57
85.0 C 8.93 C Figure 25: 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. The contour plot of temperature distribution inside the EVA layer in the first cycle during simulation is shown in Fig. 26. At time point A, the boundary temperature reaches 85 ◦C, and is kept constant for 20 hours. As shown in Fig. 26(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. 26(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. 27. As can be seen, the platform, which denotes 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 hours. Besides, the platform 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 58
loading. The contour plot history of moisture diffusion during the humidity freeze simulation is shown in Fig. 28. The moisture gradually diffuses into the central area, which is very similar to the phenomena obtained from previous damp heat simulation. However, the 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. 29, the difference of moisture concentration spatial variation after 700 hours and 1000 hours 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 hours and 1000 hours 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. 29 that the amount of moisture concentration at all the areas of EVA after 700 hours in the damp heat case is even higher than that in the humidity freeze case after 1000 hours, which further demonstrates the effect of temperature on the 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. 30. 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. 2.2.3 Thermal cycling test In this section, as shown in Fig. 31, a minimodule composed of 3×3 multicrystalline solar cells with a crack in the middle one, which was 59
85.0 C 8.93 C 85.0 C 67.8 C 85.0 C 81.1 C (a) Time point A 65.9 C -40 C -40 C 58.2 C -40 C 50.9 C (b) Time point B (c) Time point C (d) Time point D (e) Time point E (f) Time point F Figure 26: Temperature contour plot inside EVA at six different time points in the first cycle during humidity freeze simulation. 60
0 4 8 1 2 1 6 2 0 - 6 0 - 3 0 0 3 0 6 0 9 0 T e m p e r a t u r e ( °C ) D is t a n c e f r o m e d g e ( m m ) T i m e P o i n t A T i m e P o i n t B T i m e P o i n t C T i m e P o i n t D T i m e P o i n t E T i m e P o i n t F Figure 27: Temperature versus distance from edge curves at six different time points in the first cycle during humidity freeze simulation. tested under the thermal cycling environmental condition with EL images taken regularly as reported in [130], was simulated using the proposed modeling framework. In this case, the temperature ramp is given by T∗(t) = t∗ 1−t t∗ 1(T∗ 1−T∗ 2) + T∗ 20≤t<t∗ 1 T∗ 2t∗ 1≤t<t∗ 2 t−t∗ 2 t∗ 3−t∗ 2(T∗ 3−T∗ 2) + T∗ 2t∗ 2≤t<t∗ 3 T∗ 3t∗ 3≤t<t∗ 4 t∗ 5−t t∗ 5−t∗ 4(T∗ 3−T∗ 1) + T∗ 1t∗ 4≤t<t∗ 5 (2.85) where T∗ 1= 25 ◦C, T∗ 2=−40 ◦C, T∗ 3= 85 ◦C, and t∗ 1= 0.5h, t∗ 2= 1.5h, t∗ 3= 2.5h, t∗ 4= 3.5h, t∗ 5= 4.0h, as shown in Fig. 32. The sketch of the cross-section of the PV minimodule is shown in Fig. 33. The module consists of a tempered glass with the thickness of 4 mm, an encapsulanting EVA with the thickness of 0.5 mm, the silicon solar cell with the thickness of 0.166 mm, another layer of EVA, and finally a thin backsheet made of an ethylene tetrafluoroethylene core and silicon nitride coating with the 61
(a) 100 h (g) 800 h (d) 500 h (b) 300 h (c) 400 h (e) 600 h (f) 700 h (h) 900 h (i) 1000 h 0.0000 g/cm30.0056 g/cm3 Figure 28: Contour plot of moisture diffusion along the EVA layer during the humidity freeze simulation 62
0 4 8 1 2 1 6 2 0 0.000 0.001 0.002 0.003 0.004 0.005 0.006 M o i s t u r e c o n c e n t r a t i o n ( g / c m 3) D is t a n c e f r o m e d g e ( m m ) H u m i d i t y F r e e z e _ 7 0 0 h H u m i d i t y F r e e z e _ 1 0 0 0 h Dam pHeat_700 h Dam pHeat_1000 h Figure 29: Comparison of moisture concentration versus distance from edge curves after 700 hours and 1000 hours between the damp heat and humidity freeze simulations. thickness of 0.1 mm. The thermal and mechanical properties of all the component materials are listed in Table 1. The size of this PV module is 48 cm, and the adopted mesh size of the finite element model is 1 mm. Since backsheet is permeable to water and thus moisture can penetrate from it and percolates along the edges of each solar cell, it is possible in the numerical simulation to directly impose moisture boundary condition at the edges of each solar cell embedded in the module, see Fig. 33. The temperature contour plots of the cross section of the minimodule at time point B and time point E are shown in Fig. 34. Owning to the different heat conductivity properties of different layers (glass, backsheet, and silicon solar cells), the temperature inside the PV minimodule presents significant spatial variation at both the cooling and heating stages. The temperature distribution inside the EVA layer of the module at six time points of the first cycle during the thermal cycling simulation is shown in Fig. 35. At time point A, heat has diffused inside the panel. From time point A to C, the thermal boundary condition is kept constant 63
0 200 400 600 800 1000 0.000 0.001 0.002 0.003 0.004 0.005 M o i s t u r e c o n c e n t r a t i o n ( g / c m 3) T i m e ( h ) H u m i d i t y F r e e z e _ 3 m m H u m i d i t y F r e e z e _ 5 m m H u m i d i t y F r e e z e _ 8 m m D a m p H e a t _ 3 m m D a m p H e a t _ 5 m m D a m p H e a t _ 8 m m Figure 30: 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. with the temperature of 85 ◦C, and as can be clearly seen, the EVA layer is almost fully heated at the end of this heating period. When the temperature boundary condition drops to -40 ◦C at time point D, this layer starts to cool down from its edges, but the temperature in the most area is still very high after this period, and the lowest temperature is up to around 50 ◦C in the time point F, see Fig. 35(f). This trend is quantified in Fig. 36 by plotting the temperature versus the distance from edge at six time points during the first cycle. Compared with the previous curves obtained from humidity freeze simulation, the platforms are relatively lower at the temperature level of around 0 ◦C because of the much shorter heating period with the constant temperature boundary condition 85 ◦C. The contour plots of moisture concentration inside the EVA layer above the silicon solar cell layer during the thermal cycling simulation are shown in Fig. 37. Notably, moisture diffusion can take place not only from the edges of module, but also along the interspaces between different solar cells as well as the crack indicated in Fig. 31. Thus the moisture 64
120th cycle 360th cycle 240th cycle 160th cycle 200th cycle 280th cycle 320th cycle 400th cycle 460th cycle Figure 31: EL images of the minimodule taken at different cyles during the thermal cycling test (adapted from [130]). 65
T e m p e r a t u r e ( °C ) T i m e ( h ) ABC D E F t*1 0t*2t*3t*4t*5 T*1 T*2 T*3 Figure 32: Temperature profile imposed inside the chamber during the thermal cycling test Si Cell Glass H2O H2O EVA EVA Backsheet H2O H2O H2OH2OH2O H2O H2O Figure 33: The sketch of the PV minimodule in the thermal cycling test. 85.0 C 3.01 C -3.26 C -40 C (a) Time point B (b) Time point E Figure 34: 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. 66
Chapter 3 Global-local phase field modeling of silicon cell cracking in photovoltaic laminates Nowadays the crystalline silicon photovoltaic (PV) modules have been massively deployed all over the world since the last century, and according to [132], the installed capacity of PV approximately increases from 40 GW to 715 GW in the recent 10 years. Among the major components of PV products, including tempered glass, ethylene-co-vinyl acetate (EVA), silicon cells, and backsheet, the silicon wafer accounts for more than 40 % of the manufacturing cost of crystalline modules [48]. It was proposed in the 2010 International Technology Roadmap for Photovoltaics that the thickness of silicon cell should be signifcantly decreased so as to save cost of cystalline PV modules. However, the decrease of silicon cell thickness could reduce its robustness under mechanical loading and thus lead to the formation of microcracking. Hence, the identification of crack initiation and propagation that may trigger performance degradation of PV modules has received a lot of attention in the recent decades. To accurately predict the failure and crack propagation, several for- 73
mulations on the basis of different numerical techniques have been developed in the past decades. Simplified phenomenological cohesive zone models specific for the crack modeling along internal boundaries have been proposed in [133, 134, 135, 136, 137], in which cracking events are triggered by the evaluation of traction-separation constitutive laws. This modeling technique has been widely employed to identify crack paths along the element edges by incorporating a characteristic length scale, see the discussions in [138, 85, 139, 140, 141], among many others. In contrast to the cohesive zone models, the extended finite element method has also been developed to model crack propagation, and it does not explicitly depend on the finite element discretization corresponding to the physical domain [142]. Other alternative computational techniques such as enhanced finite element method [143, 144] and generalized finite element method [145], can also be used to model crack events although possible operative difficulties to identify crack initiation and paths might arise. In order to overcome the disadvantages of the above-mentioned explicit modeling method for complex crack topologies, the phase field approach [146, 147] that is based on the Griffith’s theory has been proposed to address the crack modeling of quasi-brittle materials such as silicon solar cells [148, 149]. In this approach, the sharp crack is generally diffused through the definition of so-called phase field variable and the crack propagation is characterized by the evolution of the corresponding governing equations. According to the Γ-convergence discussions in [150], the crack discontinuities are regularized by a characteristic phase field length scale. Note that the phase field formulation shares several common aspects with the gradient enhanced damage formulations, which is one of the most appealing advantages of this crack modeling method [151, 152, 153]. In view of this feature, the phase field approach accounting for thermodynamic consistency [146, 154] has been extended to the modeling of ductile fracture [155], anisotropic fracture [156], multiphysics fracture [157], and polycrystalline materials [158, 159], among many others. Given the promising aspects of phase field approach, this technique has been employed to model the crack propagation of brittle silicon solar cells in this work. 74
As mentioned above, to reduce the production cost of silicon solar panels, the average thickness of solar cells has decreased from 300 µm to 150 µm during the recent decades [160, 161, 162], which increases the breakage rates and crack events. Based on the computational framework developed in [163], the solid shell formulation combined with anisotropic phase field model at finite deformation is proposed to model the preferential crack propagation in thin-film solar cells. In the literature, phase field models have been coupled with shell kinematics for the fracture modeling of thin-walled structures through many different ways, see [164, 165] and references therein. Previous attempts in [166] have been made to model geometrically linear problems according to the Reissner- Midlin theory, but the developments are limited to standard finite elements. Alternative investigation adopting the Local Maximum-Entropy approximation is proposed in [167], making it difficult to be implemented into the standard finite element method. The tension-compression split method by integrating the energy contribution through the shell element thickness has been proposed in [168, 169], and later extended to the isogeometric modeling of multipatch structures [170, 171]. Recently, this approach has been adopted in [172, 173] for phase field fracture modeling of the Reissner–Mindlin shells and plates. Besides, other models have been developed in combination of phase field theory with solid shell and Kirchhoff–Love formulations to model brittle and ductile fracture [174, 175, 176, 177]. Hence, the attempt to employ the phase field solid shell formulation, extended to account for the anisotropic fracture orientation, is a proper choice for the crack modeling of solar cells in thin-walled PV laminates . This chapter aims at presenting a global-local phase field approach with large-deformation solid shell element formulation so as to provide a very efficient modeling framework for the modeling of silicon solar cell crack propagation in the PV modules, and the main contents are taken from [163]. The two main innovative aspects of this chapter are: (i) the coupling of phase field approach for fracture modeling with solid shell formulation considering the EAS and ANS methods for the alleviation of locking effects, and its efficient and robust quasi-Newton monolithic 75
implementation into the commercial finite element code ABAQUS, which was comprehensively evaluated by comparison with the standard staggered Newton scheme through different benchmark examples, including single edge notched tension and shear, fracture of cylindrical structure under mixed bending and tension, and computationally demanding fatigue induced crack growth; (ii) a specific global-local approach in the 3D setting tailored for phase field modeling with solid shell elements is proposed here to save computational cost in fracture modeling of largescale thin-walled structures, and its capability was demonstrated by simulating the local crack growth of cylindrical structure under both static and fatigue cyclic loading. Finally, to model the preferential crack patterns observed from the experimental electroluminescence images of solar cells, the fracture anisotropy is taken into account in the phase field modeling, and it is showcased through the simulation of crack growth of silicon solar cells when the PV modules are subjected to the different loading, as depicted in Section 3.4.7. 3.1 Phase field solid shell formulation In this section, the phase field solid shell formulation is presented for the modeling of anisotropic crack propagataion in the thin-film brittle solar cells. Firstly, the basics of phase field approach are revisited in Section 3.1.1. For a more comprehensive description, the detailed derivation of this theory can be found in the significant work [146, 178, 179]. In Section 3.1.2, the solid shell kinetics at finite deformation are described in the convective curvilinear coordinate system. 3.1.1 Basics of phase field approach In the 3D setting, let B0⊂R3and Bt⊂R3denote the reference configuration and current configuration of an arbitrary solid body with existing crack represented by Γ, and Xand xstand for the position vectors in the corresponding configurations, respectively. Furthermore, let ∂B0 and ∂Btbe the exterior boundary of the solid body in the reference and 76
X 𝚪 𝟐𝒍 Y Z 𝛀 𝛀 𝚪𝒍 (a) (b) Figure 40: Schematic representation of solid body with (a) sharp crack topology, and (b) phase field approximation of diffusive crack current configuration, respectively. The motion of the material point X inside the body is denoted by φ(X, t) : B0×[0, t]→R3that maps into the corresponding position xin the current configuration during the time interval [0, t]. The deformation gradient Fuis defined as Fu:= ∂Xφ(X, t), where ∂Xrepresents the partial derivative with respect to the position X in the reference configuration, and its determinant Ju=det[Fu]denotes the Jacobian. The phase field approach is conceived as the regularization of sharp crack topology within a diffusive crack zone characterized by the scalar-valued function d:B0×[0, t]→[0,1] to model brittle fracture. The phase field variable dsmears out the sharp crack by a diffusive crack area of width l, as shown in Fig. 40. It is noted that the width of regularization area depends on the parameter l, which is called the phase field length scale that controls the transition between intact and damaged parts. In the modeling framework, dis regarded as a smooth function of (X, t), and d= 0 and d= 1 denote intact and cracked states, respectively. According to the variational theory of fracture [179], the total energy functional that governs the cracked body under external loadings takes the form of Π(u,Γ) = Πint (u,d,Γ) + Πext(u)(3.1) where udenotes the displacement field, and Πint and Πext represent the internal and external energy functionals, respectively. Note that the internal energy functional term Πint is defined as the sum of elastic energy stored in the solid body ΠBand the energy dissipated through crack propagation ΠΓ. Based on the Griffith’s theory, the internal energy func- 77
tional Πint is defined as Πint (u,d,Γ) = ΠB(u,d)+ΠΓ(Γ) = ZB\Γ Ψ(u,d) dΩ + ZΓGcdΓ (3.2) where Ψdenotes the specific elastic energy, and Gcis the critical energy release rate. The dissipated fracture energy during the crack propagation is evaluated through the Griffith theory. It should be pointed out that the competition between the elastic energy in the solid body and fracture energy is directly defined in the context of minimization problem. The evaluation of this competition is computationally challenging when the space discretization methods are used to track the crack propagation. To circumvent the use of tracking algorithms, the phase field approach that smears the crack over the whole domain of the body is employed in the same line with gradient damage formulations [150, 179, 180]. In the context of this approach, the dissipated surface energy can be approximated by a volume integral ZΓGcdΓ ≈ZB0Gcγ(d,∇Xd) dΩ (3.3) where γ(d,∇Xd)is the crack surface density per unit volume, and ∇Xd is the gradient of phase field variable. In this way, the sharp crack is regularized over the body inducing the diffusive representation, as shown in Fig. . The crack surface density γ(d,∇Xd)is defined according to the modified Ambrosio-Tortorelli functional [181] γ(d,∇Xd,ω) = 1 2ld2+l 2ω: (∇Xd⊗∇Xd)(3.4) where ωis the second order structural tensor characterising the material anisotropy. In order to make the energy release rate orientationdependent [158, 182], the structural tensor can be defined as ω=I2+αp(I2−Np⊗Np)(3.5) where I2is the second-order identity tensor, Nprepresents the unit vector normal to the preferential cleavage plane, and αpis the penalty parameter that is used to prevent damage to develop on planes not normal 78
to Np. According to the previous study [183], this parameter should be greater than 1.0, and in case of isotropic material, the value is equal to 0.0. The specific elastic energy Ψis influenced by the phase field variable d, which is motivated by the fact that this energy in the damage transition zone has to decrease so as to ensure the thermodynamic equilibrium. Furthermore, to avoid the asymmetric damage behaviour, the elastic energy density can be split into the active contribution Ψ+and the passive contribution Ψ−[178]. During the reversal loading process, the crack closing precludes the damage evolution, and also provokes the stiffness discovery [184]. Based on the aforementioned viewpoints, the elastic energy Ψcan be defined as Ψ = g(d)Ψ++ Ψ−(3.6) where g(d) = (1 −d)2+Kdenotes the monotonic degradation function, K ≈ 0is a positive parameter to ensure numerical stability in case of fully material degradation. 3.1.2 Solid shell kinematics In the concept of solid shell theory, the position vectors Xand xin the reference and current configurations, and the phase field variable dat any material point can be approximated by the corresponding vectors on the bottom and top surfaces of solid-like shell element, and can be defined as Xξ1, ξ2, ξ3=1 21 + ξ3Xtξ1, ξ2+1 21−ξ3Xbξ1, ξ2(3.7a) xξ1, ξ2, ξ3=1 21 + ξ3xtξ1, ξ2+1 21−ξ3xbξ1, ξ2(3.7b) dξ1, ξ2, ξ3=1 21 + ξ3dtξ1, ξ2+1 21−ξ3dbξ1, ξ2(3.7c) where the parametric space is identified as: A:= {ξ=ξ1, ξ2, ξ3∈R3 | −1≤ξi≤+1; i= 1,2,3}, the subscript b and t represent bottom and 79
G1 G2 G3 E1E2 E3Xt Xb e1e2 e3 xt xb Xx 𝚪𝐜 𝛏𝟏 𝛏𝟐 𝛏𝟑 Midsurface 𝑩𝟎𝑩𝐭 𝛗(𝚪𝐜, 𝐭) 𝟐𝒍 g2 g3 g1 𝟐𝒍 𝛗 𝐗, 𝐭 ∶ 𝐗 → 𝐱 𝐅u≔𝛛𝐱 𝛛𝐗 = 𝐠𝐢⊗ 𝐆𝐢 Figure 41: Definition of the cracking shell body in both the reference and current configurations with the solid shell concept and phase field approach for fracture modeling top surfaces, respectively, and ξ1, ξ2, ξ3represent the coordinates in the parametric space, see Fig. 41. The kinematics of solid shell element can be described by the use of the convective curvilinear system so that the ANS interpolation for the transverse normal and shear strain components can be implemented. The covariant tangent vectors Gi(ξ)and gi(ξ)in the reference and current configuration are defined as the partial derivatives of corresponding position vectors Xand xwith respect to the convective coordinates ξi Gi(ξ) := ∂X(ξ) ∂ξi,gi(ξ) := ∂x(ξ) ∂ξi, i = 1,2,3(3.8) The contravariant basis vectors can be determined in a standard manner by Gi·Gj=δj iand gi·gj=δj i, and metric tensors are defined as G=GijGi⊗Gj=GijGi⊗Gj,g=gijgi⊗gj=gijgi⊗gj. In curvilinear setting, the deformation gradient Fuis given by Fu=∂x ∂X=gi⊗Gi(3.9) where the Einstein summation convention on repeated indices is adopted here. Through the definition of metric tensor components Gij =Gi·Gj and gij =gi·gjin the reference and current configuration, the displace- 80
ment derived Green-Lagrange strain tensor is defined as Eu:= 1 2h(Fu)TFu−I2i=1 2[gij −Gij]Gi⊗Gj(3.10) The energetically conjugated second Piola-Kirchhoff stress tensor is defined as S=SijGi⊗Gj(3.11) where Sij represents the contravariant component. 3.2 Numerical implementation In this section, the numerical strategy rooted in the use of finite element method for the spatial approximation is described briefly. Since this work is restricted to quasi-static analysis, no temporal integration scheme is required, which leads to an equilibrium problem at each pseudotime step. Firstly, the variational formulation and finite element interpolations are derived in Section 3.2.1. Secondly, the solution schemes proposed for the phase field solid shell formulation are depicted in Section 3.2.2. 3.2.1 Finite element interpolation The mixed Hu-Washizu variational principle is adopted for the derivation of phase field solid shell formulation incorporating the Enhanced Assumed Strain (EAS) and Assumed Natural Strain (ANS) methods to alleviate the locking pathologies. It should be pointed out that the EAS method is employed to remedy volumetric and Poisson thickness locking, while the membrane and in-plane locking effects are tackled by the ANS method [111, 114]. Given the enhancement based on the EAS method at finite deformation, the Green-Lagrange strain consists of two parts following the approach proposed in [107], including the displacement derived compatible strain Euand incompatible strain ˜ E, and its complete form reads: E=Eu+˜ E. 81
The multi-field variational functional of the solid body takes the form of Π(S,˜ E,u,d) = ZB0 Gcl 2d2 l2+ω: (∇Xd⊗∇Xδd)dΩ +ZB0 ΨEu,˜ E,ddΩ −Πext (3.12) where Πext identifies the external energy functional. Note that the displacement u, the phase field variable d, and the imcompatible strain ˜ E are the independent variables. Given the orthogonality condition between the stress and strain fields, the stress field is ignored in line with the previous work [125]. The first variation of Eq. (3.12) with respect to the independent fields is given by Ru(u, δu,˜ E,d) = ZB0 ∂Ψ ∂E:∂Eu ∂uδudΩ −δΠext(u) = 0,∀δu∈Vu (3.13a) R˜ E(u,˜ E, δ ˜ E,d) = ZB0 ∂Ψ ∂E:δ˜ EdΩ = 0,∀δ˜ E∈V˜ E(3.13b) Rd(u,˜ E,o, δd) = ZB0Gcl1 l2dδd+ω: (∇Xd⊗∇Xδd)dΩ +ZB0−2(1 −d)δdΨ+dΩ = 0,∀δd∈Vd (3.13c) where Vu,V˜ Eand Vdare the admissible spaces of independent variables. By means of isoparametric concept, the approximations of position vectors Xand x, displacement vector uand its variation δuon the solid shell element level are given by X=N˜ X(3.14a) x=N˜ x(3.14b) u=Nd, δu=Nδd(3.14c) 82
Note that, though it is not fully coupled, the subsequent approximation of the stiffness matrix couples the phase field and displacement, see Eq. (3.37). The quasi-Newton stiffness will be reformed whenever the number of iterations exceed 8 without obtaining the convergent solution in this work. The computational procedure for the quasi-Newton monolithic implementation of phase field solid shell formulation is outlined in Algorithm 3. Algorithm 3: Numerical algorithm for the QN implementation of phase field solid shell formulation Data: dn,¯ dn,∆d(k) n+1,∆¯ d(k) n+1 Result: dn+1,¯ dn+1 Initialization of ζn,Rζ int,n,kζd,n,[kζζ,n]−1; while ||R∗ d|| >tolerance do Compute ∆ζ(k) n+1 =−[kζζ,n]−1hRζ int,n +kζd,n∆d(k) n+1i; Update the enhancing vector ζ(k) n+1 =ζn+ ∆ζ(k) n+1; for n←1to N integration points do Compute the curvilinear basis G(k) n+1 and g(k) n+1; Compute the B matrices B(k) n+1 and Bd(k) n+1 ; Modify the B(k) n+1 matrix according to the ANS method; Compute C(k) n+1,S(k) n+1,Ψ(k) n+1, and H(k) n+1; Compute the EAS operator M(k) n+1; end Compute the element stiffness matrices k(k) dd,n+1,k(k) dζ,n+1,k(k) ζd,n+1,k(k) ζζ,n+1, and k(k) dd,n+1; Compute the internal force vectors Rd(k) int,n+1,Rζ(k) int,n+1, and Rd(k) int,n+1; Perform the static condensation in the element level; Do the final assembly; end 3.3 The global-local fracture submodeling approach The global-local modeling approach has been widely employed in a variety of technical applications, such as J-integral calculation in fracture problems [188], three-dimensional stress analysis [189, 126] and computation of stress intensity factors [190]. Generally, a global finite ele- 89
Global element nodes Local element nodes Upper edges Bottom edges Figure 42: The sketch of solid shell elements from both global and local models ment model with coarse mesh is used to calculate the displacement field and determine the size and location of the critical region with low computational cost. Subsequently, the displacement output obtained from the global-scale calculation will be employed to interpolate the boundary condition of the local model, which is a fine scale representation of the critical area with more detailed mechanical degradation to achieve a greater accuracy [29]. However, very limited work has been focused on the global local approach for phase field fracture modeling, see [191, 192, 193], among others. To further reduce the computational cost in fracture modeling of largescale thin-walled structures, a specific global-local submodeling approach in the 3D setting tailored for phase field modeling using solid shell formulation is proposed. The submodeling technique is a computational strategy frequently employed in fiber-reinforced composites for modeling with two different length scales, such as macro scale and meso scale of the laminate. The global and local model are solved sequentially in the submodeling approach, in which the global model is solved firstly, and subsequently its results are taken as boundary conditions to drive the local model. As a result, both global and local models are analysed separately, and local model has no influence on the global model. In this work, this approach is introduced for phase field fracture modeling of shell structures to improve computational efficiency when the crack in the local region has minor effect on the global stiffness. For instance, it is pointed out in [194, 29] that the stiffness of the global photovoltaic panel 90
is hardly degraded by the cracks in the local solar cells, and in this case, it is viable to apply the submodeling technique uncoupling the global linear elastic finite element analysis of photovoltaic module from the local fracture modeling of solar cell. To simulate the realistic loading condition driving the crack propagation of brittle silicon solar cell, it is essential to model the whole photovoltaic laminate, and on the other hand, for the prediction of cracking phenomena at the cell level, phase field approach can be employed with very refined mesh for the local model. Algorithm 4: Numerical algorithm for displacement projection from global model to local model for t←1to T time increments do Impose boundary conditions on the global model; Solve the nodal displacements of global model using Newton-Rapson scheme; for i←1to Ssolid shells of the local model do for j←1to Nexterior facets of the local model do for k←1to Medges (M=2, upper and bottom edges) do Find the two closest nodes P1and P2belonging to the same edge and facet in the global model; Interpolate the nodal displacement uk=N1u1+N2u2, where u1and u2are the nodal displacement vector of P1and P2,N1and N2are linear shape functions. end end end end To efficiently model thin-walled global structures like photovoltaic panels, a solid-like continuum solid shell element is employed to simulate the structural response under complex loading conditions with relatively coarse mesh to save computational cost. Compared with Kirchhoff- Love shell elements and solid elements in thin-walled structure modeling [85], it is preferred since three dimensional constitutive laws are allowed and accuracy can be ensured with only one ply of element required through thickness of the global model. After obtaining the displacements of global coarse-scale model in the different directions by linear elastic finite element analysis, the boundary conditions applied to the nodes belonging to the edges of local fine-scale model can be determined by linear interpolation of global results. It is worth mentioning that the finite element discretization of the local model is not required to comply 91
with the one of the global finite element model, and a projection scheme can be used to map the displacement results from computation of global model onto the edge nodes of local model. In doing so, the mesh of local model can be more refined to perform the appropriate phase field fracture simulation. From the viewpoint of algorithm, since there is only one ply of solid shell element through the thickness direction, it is convenient to distinguish the nodes of exterior element facets belonging to the local model into upper and bottom edge nodes, see Fig. 42. For each node of local model on the exterior facets, the displacement boundary conditions are obtained by linear interpolation between the displacement values of two closest nodes of global model on the same upper or bottom edge. The projection scheme is illustrated in Algorithm 4. 3.4 Numerical examples In this section, to investigate the performance of quasi-Newton monolithic solution scheme in phase field modeling using solid shell formulation, four different numerical tests, including single edge notched tension, single edge notched shear, cylindrical structure under mixed tension and bending, and fatigue induced crack propagation, are performed. The performance in attaining convergence and computational efficiency is assessed by comparison with that of the widely adopted staggered scheme through the different types of boundary value problems. Then the capability of the proposed global-local approach to further reduce the computational cost in phase field modeling is demonstrated by the numerical tests of cylindrical structure subjected to both constant and cyclic fatigue loading conditions. With regard to the application of the computational framework for crack modeling of silicon solar cells in the photovoltaic modules, the feasibility to model anisotropic fracture using the current approach is firstly shown through a demonstration problem, and then the solar cell cracking modeling in the photovoltaic module under two different loading cases are presented. 92
Figure 43: (a) Sketch of the single edge notched tension and (b) finite element mesh of the specimen 3.4.1 Single edge notched tension The first numerical case is the mode I fracture of single edge notched specimen under monotonic loading, as indicated in Fig. 43(a). It is a very popular benchmark example for verification of phase field approach and challenging for the numerical algorithms used to solve the coupled governing equations [146]. Specimen dimensions and boundary conditions are shown in Fig. 43(a). The geometry is a square plate with unit lateral size of 1 mm and thickness of 0.001 mm. Note that there is an initial sharp crack with the length of 0.5 mm in the middle of the specimen. The degrees of freedom at the bottom edges of specimen are constrained, while the displacement boundary condition is imposed on the upper edges in the vertical direction, leading to a symmetric loading at the crack tip of specimen and subsequent mode I crack propagation. The Young’s modulus and Poisson’s ratio are equal to 210 GPa and 0.3, respectively, and regarding the fracture properties, the fracture toughness and phase field length scale are set to 2.7 N/mm and 0.024 mm, respectively. The model is discretised with solid shell elements, see Fig. 43(b), and notably, the characteristic elements along the potential crack path are refined to a size of 0.004 mm, being 6 times smaller than the length scale. The total number of elements is 7222. The contour plots of phase field values during the crack propagation are shown in Fig. 44. 93
0.0 1.0 0.0055 mm 0.0057 mm 0.006 mm Figure 44: The contour plots of phase field values at the loading displacements equal to 0.0055 mm, 0.0057 mm, and 0.006 mm during the crack propagation in the single edge notched tension The constrained force versus loading displacement curves are shown in Fig. 45. The numerical results obtained from the conventional staggered or AM scheme with different numbers of time increments are also presented here for comparison. It can be seen that the staggered scheme is very sensitive to the time increment size, and with the increase of the number of increments, the obtained force versus displacement curve will gradually coincide with that solved by the quasi-Newton monolithic scheme, which requires only 100 time increments. Both the quasi- Newton monolithic scheme and staggered scheme require a lot of iterations at the critical time increment when the crack starts to propagate. The total number of iterations versus loading displacement curves for the staggered schemes with different increment size and quasi-Newton monolithic scheme are shown in Fig. 46. Results indicate that reproducing the numerical result obtained by the quasi-Newton monolithic scheme through the staggered scheme requires at least around ten thousand iterations. It is also shown that the number of iterations by the staggered scheme with one thousand time increments is almost the same as that required by the quasi-Newton monolithic scheme, but its deviation from the convergent solution reaches around 20%. Note that the quasi-Newton monolithic scheme is roughly 15 times more efficient than the staggered scheme for the phase field modeling of this single edge 94
0.000 0.002 0.004 0.006 0.008 0 . 0 0 0 . 1 4 0 . 2 8 0 . 4 2 0 . 5 6 0 . 7 0 F o r c e ( N ) D is p l a c e m e n t ( m m ) A M 1 0 3 i n c . A M 2 X 1 0 3 i n c . A M 1 0 4 i n c . Q N 1 0 2 i n c . Figure 45: The force vs. displacement curves obtained from both the quasi- Newton monolithic and staggered schemes for the phase field modeling of single edge notched tension using solid shell formulation notched tension using enhanced assumed strain shell formulation. The significant differences in computational cost are mainly due to the greatly reduced number of iterations in the quasi-Newton monolithic solution. 3.4.2 Single edge notched shear In this section, the performance of the proposed quasi-Newton monolithic scheme in phase field modeling using solid shell element formulation is assessed in the context of mode II fracture of single edge notched specimen. The same material properties and specimen dimensions as the previous single edge notched tension are adopted but the shear loading is applied to the upper edges of specimen, see Fig. 47(a). The sheardominated crack conditions lead to crack propagation towards the right bottom part of the specimen, which is discretized with uniform solid shell elements. The contour plots of phase field values during the crack propagation in the single edge notched shear test are shown in Fig. 48, 95
0.000 0.002 0.004 0.006 0.008 1 1 0 100 1000 10000 T o t a l i t e r a t i o n s D is p l a c e m e n t ( m m ) A M 1 0 3 i n c . A M 2 X 1 0 3 i n c . A M 1 0 4 i n c . Q N 1 0 2 i n c . Figure 46: The total number of iterations vs. loading displacement curves obtained from both the quasi-Newton monolithic and staggered schemes for the phase field modeling of single edge notched tension using solid shell formulation Figure 47: (a) Sketch of the single edge notched shear and (b) finite element mesh of the specimen 96
0.0 1.0 0.0085 mm 0.010 mm 0.0135 mm Figure 48: The contour plots of phase field values at the loading displacements 0.0085 mm, 0.010 mm, and 0.0135 mm during the crack propagation in the single edge notched shear and it can be seen that the predicted crack trajectory agrees with results in the literature, see [146]. The constrained force versus displacement curves obtained from both the quasi-Newton and staggered schemes with different increment sizes are shown in Fig. 49. It can also be seen that reproducing the quasi- Newton monolithic result with staggered scheme requires the use of 104 number of increments. When using the less number of total increments, such as 103or 2×103, the resultant curve deviates from the predicted quasi-Newton one, see Fig. 49. The total number of iterations versus displacement curves obtained from both the quasi-Newton monolithic and staggered schemes for the phase field modeling of single edge notched shear test using solid shell element formulation are shown in Fig. 50. Difference of total number of iterations between the two solution schemes is smaller compared with the previous single edge notched tension, but remains significant in this case. To obtain the same result as the quasi- Newton monolithic solution, the required total number of time increments for the staggered scheme is still one order larger than the former scheme. It can be concluded that the quasi-Newton monolithic scheme in the single edge notched shear is also much more efficient as its computational cost is around ten times smaller than that required in the staggered solution with 104number of increments. 97
0.000 0.003 0.006 0.009 0.012 0.015 0 . 0 0 0 . 0 7 0 . 1 4 0 . 2 1 0 . 2 8 0 . 3 5 F o r c e ( N ) D is p l a c e m e n t ( m m ) A M 1 0 3 i n c . A M 2 X 1 0 3 i n c . A M 1 0 4 i n c . Q N 1 0 2 i n c . Figure 49: The force vs. displacement curves obtained from both the quasi- Newton monolithic and staggered schemes for the phase field modeling of single edge notched shear using solid shell formulation 3.4.3 Cylindrical structure under mixed tension and bending The performance of the quasi-Newton monolithic scheme for phase field modeling using enhanced assumed strain shell formulation is further assessed in capturing fracture events of thin-walled structure with curved geometry under mixed tension and bending. In this section, we consider a quarter of cylindrical structure with an initial crack notch in the middle of the specimen, see Fig. 51(a). The dimensions of the specimen are as follows: the length is equal to 240 mm, the internal radius is 100 mm, the external radius is 100.1 mm, the width of the initial notch is 3 mm and its length is 6 mm. To comply with the symmetric boundary conditions, the translational degrees of freedom at the X=0 side of the specimen is constrained in the X direction, and the Y=0 side is constrained in the translational vertical direction Y. Besides, the translational degree of freedom in the Z direction at the back of the specimen is constrained, 98
Figure 55: (a) Sketch of the single edge notched fatigue and (b) finite element mesh of specimen +U -U 0T 2T 3T 4T U = 0.002 mm T = 0.002 s Time (s) Displacement (mm) Figure 56: Cyclic fatigue loading with equal tension and compression for the single edge notched specimen 105
Fig. 55(b), and the total number of solid shell elements is 4878. The material and fracture properties are also the same as the previous study. The cyclic loading with equal tension and compression is applied to the upper surface of the specimen, and the amplitude and cycle are 0.002 mm and 0.002 s, respectively, see Fig. 56. The curves of crack length driven by fatigue loading versus number of cycles for quasi-Newton monolithic and staggered schemes with different time increment sizes are shown in Fig. 57. It is clear that extremely small time increment size is required for the staggered scheme so as to obtain the same result as that of the quasi-Newton monolithic solution, and convergence rate of the staggered scheme is quite low. The curves of total iterations versus number of cycles for the two solution schemes in the phase field fatigue simulation are shown in Fig. 58. Results show that the quasi-Newton solution requires almost the same number of total iterations as the staggered solution with 80 time increments per cycle. However, even the staggered scheme with 200 time increments per cycle cannot reproduce the quasi-Newton monolithic solution, and obviously more increments per cycle are required, which is computationally too expensive. Notably, the staggered solution with 200 increments per cycle requires 6 times more computation times than that in the quasi-Newton monolithic analysis. It can be concluded that quasi-Newton monolithic scheme shows the tremendous advantage over the staggered scheme in computational efficiency for fatigue problems. The contour plots of phase field values in fatigue induced fracture modeling with solid shell formulation solved by the quasi-Newton monolithic scheme are shown in Fig. 59. As can be clearly seen, the crack propagation starts from the 15th cycle of fatigue loading, and it proceeds with the increase of fatigue cycles. The single edge notched specimen for fatigue modeling is fully cracked at the 50th cycle. 106
0 2 0 4 0 6 0 8 0 0 . 0 0 . 1 0 . 2 0 . 3 0 . 4 0 . 5 C r a c k l e n g t h ( m m ) N u m b e r o f c y c l e s A M 1 6 i n c . / c y c l e A M 3 2 i n c . / c y c l e A M 8 0 i n c . / c y c l e A M 2 0 0 i n c . / c y c l e Q N 1 6 i n c . / c y c l e Figure 57: The curves of crack length vs. number of cycles obtained from both the quasi-Newton monolithic and staggered schemes for the phase field fracture modeling under fatigue loading 3.4.5 Global-local fracture modeling of thin-walled structure with curved geometry In this section, the global-local phase field approach using solid shells is tested through the fracture modeling of cylindrical structure with an initial notch in the center as shown in Fig. 60(a). The dimensions of this specimen in the global level are the same as that in Section 3.4.3. As shown in Fig. 60(a), the global specimen is separated into 9 different blocks in total and the initial notch lies within the central block, which is chosen as the region of local model. To save computational cost, the global model is discretized with relatively coarse mesh using the uniform size of 3 mm, while to appropriately capture the crack propagation in the phase field modeling of local region, very refined mesh is used to discretize the potential crack area around the notch of local model. Specifically, the mesh size around the edge of local model is 3 mm, which is the same as that of global model, and characteristic mesh size along the 107
0 2 0 4 0 6 0 8 0 1 1 0 100 1000 10000 T o t a l i t e r a t i o n s N u m b e r o f c y c l e s A M 1 6 i n c . / c y c l e A M 3 2 i n c . / c y c l e A M 8 0 i n c . / c y c l e A M 2 0 0 i n c . / c y c l e Q N 1 6 i n c . / c y c l e Figure 58: The curves of total iterations vs. number of cycles obtained from both the quasi-Newton monolithic and staggered schemes for the phase field fracture modeling under fatigue loading Cycle 15 Cycle 20 Cycle 25 Cycle 30 Cycle 40 Cycle 45 Cycle 50Cycle 35 0.0 1.0 Figure 59: The contour plots of phase field values in fatigue modeling using solid shell element solved by the quasi-Newton monolithic scheme 108
Figure 60: (a) Geometry and (b) finite element mesh of the global and local models in the phase field fracture submodeling of cylindrical structure potential crack path is 0.8 mm, see Fig. 60(b). The continuum solid shell is chosen for the elastic finite element analysis of global model, and phase field modeling for solid shell using quasi-Newton monolithic scheme is subsequently employed for the analysis of refined local model. The same material properties are adopted for the global and local model with the Young’s modulus equal to 210 GPa and the Poisson’s ratio equal to 0.3, and both the phase field length scale and fracture toughness of local model are the same as the previous study in Section 3.4.3, which are equal to 4 mm and 2.7 N/mm, respectively. The contour plots of phase field values of local model, displacement fields of both global and local models solved by quasi-Newton monolithic scheme at two different time points are shown in Fig. 61. The first time point is the time increment when the crack of local model starts to propagate under the drive of the boundary conditions interpolated from the displacement output of global model, while the second time point is the time increment when the local model with curved geometry is fully cracked. This global-local approach is very computationally efficient in fracture modeling of large thin-walled structures, since only the local re- 109
Type Time Time 1 Time 2 Displacement of global model Displacement of local model Phase field of local model Figure 61: The contour plot of global and local displacement as well as phase field variable in the fracture submodeling of cylindrical structure with solid shell element formulation gion of interest needs to be modeled by phase field approach. In the global model, the degree of freedom corresponding to phase field variable is ignored, and thus significant computational cost can be saved. Besides, mesh in the local fracture region of interest can be refined with more flexibility, which is very appealing for realistic application of phase field approach in fracture modeling of large-scale structures. To further assess the global-local phase field approach using solid shells, the above case study is extended to the analysis of fatigue induced crack growth by applying the cyclic loading to the frontal end of global model. The cyclic loading is characterized by the amplitude of 0.05 mm and the cycle of 0.01 s, as shown in Fig. 62. The time duration is 1 s, and thus there are 100 cycles in total. The phase field plots of local model in the fatigue submodeling solved by quasi-Newton monolithic scheme are shown in Fig. 63. It should be pointed out that the boundary conditions of local model are obtained by the interpolation of the displacement output from the global finite element analysis, which is subjected to cyclic loading shown in Fig. 62. In this fatigue submodeling case, the crack of local model starts to propagate from the 40th cycle and after 28 cycles of 110
+U -U 0T 2T 3T 4T U = 0.05 mm T = 0.01 s Time (s) Displacement (mm) Figure 62: Cyclic fatigue loading with equal tension and compression for the global model of the cylindrical structure fatigue loading, the local model is fully cracked along the expected crack path, which further demonstrates the capability of this global-local phase field approach for fracture modeling. 3.4.6 Benchmark example of anisotropic fracture In this example, the model will be assessed through a benchmark problem consisting in a square specimen that contains two parts and an initial crack, as shown in Fig. 64. The size of the specimen is 1 mm, and meshed with uniform solid shell elements. The elastic modulus and critical fracture energy release rate are set to 160 GPa and 4.32 N/m, respectively, which are the typical material properties of silicon solar cells in photovoltaic modules [148]. The preferential crack plane of the Grain 2 with an initial crack is oriented with a fixed angle of 45 degree, while the preferential crack direction of the Grain 1 varies with respect to the horizontal axis of the whole domain. Regarding the boundary condition, the bottom facet is constrained in the vertical direction, and a tensile loading is applied to the upper facet of the specimen. The predicted contour plots of phase field with different preferential crack orientations of Grain 1 are shown in Fig. 65. Different crack patterns can be observed with different preferential crack planes. The crack propagates from the initial crack following the fixed crack orientation 111
Cycle 40 Cycle 44 Cycle 48 Cycle 52 Cycle 60 Cycle 63 Cycle 68Cycle 56 0.0 1.0 Figure 63: The contour plots of phase field values of local model in fatigue submodeling solved by quasi-Newton monolithic scheme Grain boundary Crack Displacement loading Grain 2 Grain 1 45 degree Preferential direction Preferential direction Grain boundary Figure 64: Schematic diagram of the benchmark problem for anisotropic fracture modeling. 112
0.0 1.0 0 degree 22.5 degree 45 degree 67.5 degree Figure 65: The phase field contour plots of fully cracked specimen with the preferential crack orientation of Grain 1 equal to 0 degree, 22.5 degree, 45 degree, and 67.5 degree. 113
through the domain of Grain 2, and when approaching the boundary between the two grains, crack deviation occurs due to the different nature of Grain 1. As shown in Fig. 65, the crack path through the Grain 1 when the specimen is fully cracked exactly follows the expected orientation. The force versus displacement curves for the different crack orientations of Grain 1 are plotted in Fig. 66. It can be seen that the crack preferential orientation of Grain 1 has negligible influence on the mechanical response. The approximate linear behaviour can be observed before the specimen is fully cracked when the loading displacement reaches around 0.0005 mm. After the peak values, all the curves drops down to zero immediately except the curve with crack orientation angle of Grain 1 equal to 67.5 degree, which can be ascribed to the fact that part of specimen is still constrained even when reaching the fully cracked state in this case. It should be pointed out that the different crack patterns can create different electrically inactive area in the silicon solar cell, leading to different power loss of the photovoltaic modules. Hence, the fracture anisotropy must be taken into account for the modeling of cracking events [29]. 3.4.7 Silicon solar cell cracking in the photovoltaic module The modeling framework is further assessed in this section with application to the more realistic simulation of silicon solar cell cracking when the photovoltaic modules are subjected to the different external loading. According to the experimental investigations [194, 29], the cracking events of silicon solar cells have negligible influence on the global stiffness of photovoltaic modules, and as a result, the global-local approach proposed in [163] can be employed to reduce the computational cost by decoupling the elastic modeling of photovoltaic module with the fracture modeling of silicon solar cells. In this approach, the photovoltaic module is assumed to be purely elastic corresponding to the global model, and its solution at each time step is used to drive the local fracture modeling of silicon solar cell with the phase field approach. Hence, the global and local models are solved in a staggered manner, which features the one- 114
-76.55 mm 72.04 mm Loading Displacement 67.5 mm Loading Displacement 70.5 mm Loading Displacement 73.5 mm Loading Displacement 76.5 mm Figure 72: The predicted displacement contour plots of the whole photovoltaic module when subjected to three point bending at the loading displacement 67.5 mm, 70.5 mm, 73.5 mm, and 76.5 mm. 121
0.0 1.0 Loading Displacement 70.5 mm Loading Displacement 73.5 mm Loading Displacement 76.5 mm Loading Displacement 67.5 mm Figure 73: The predicted phase field contour plots of the local model with no fracture anisotropy at the loading displacement 67.5 mm, 70.5 mm, 73.5 mm, and 76.5 mm. 122
0.0 1.0 Loading Displacement 103.5 mm Loading Displacement 106.5 mm Loading Displacement 109.5 mm Loading Displacement 100.5 mm Figure 74: The predicted phase field contour plots of the local model with the preferential crack orientation angle equal to 45 degree at the loading displacement 100.5 mm, 103.5 mm, 106.5 mm, and 109.5 mm. 123
Chapter 4 Modeling and simulation of the recycling process for end-of-life photovoltaic laminates In the past decades, the development of environmentally friendly energy production techniques has resulted into the fast-growing increase of photovoltaic (PV) modules all over the world [27, 28, 196, 55, 197]. In view of the limited service lifetime (around 25 years) of PV modules [198, 85, 24, 29], a significant amount of these end-of-life products, which are defined as waste of electronic equipment [199], will continuously pose challenges to the sustainable energy production and environmental protection [200, 201, 202, 58, 1]. Hence, there is an urgent need to develop reliable techniques to dispose end-of-life PV modules [203, 204, 205, 206, 207, 208]. EVA is a kind of polymeric material consisting of 3D network of polymer molecules that are bonded to each other, and thus capable of recoverable large-deformation [209, 210, 211, 212]. During the recycling process using the solvent method, the solvent can penetrate into the EVA network when coming into contact, and the molecules can be absorbed by 124
the skeletal polymeric network of EVA, forming a swollen aggregate. The solvent molecules can interact with each other as well as with the polymer structure, and also migrate in and out of the swollen aggregate. The theory of mass transport in the continuum solids dates back to [213], and have been further developed for the description of coupled deformation and diffusion with application to the response of hydrogels, including swelling and deswelling induced by mechanical deformation and forced fluid permeation [214, 215, 216]. The purpose of this work is to develop a thermodynamically consistent large-deformation theory on the continuum level for the description of EVA polymer-solvent system in the PV recycling, so as to provide guidance to the design of experimental procedures for the nondestructive recovery of ultrathin and brittle silicon wafers. Similarly to the theories developed by [214], the polymer-solvent system is treated as homogenized continuum in which the solvent flux is allowed. In the development of constitutive theory, the Flory-Huggins theory for the free energy change due to the mixing of solvents with different solubility parameters and EVA polymer network is adopted, and the statistical-mechanical model is accounted for to describe the entropic change of EVA polymer chains due to the mechanical stretching. The adhesion strength of PV modules is greatly influenced by the operating environment, depending on different factors, including humidity, temperature, and irradiance [217, 218]. Moisture will weaken the ionic attraction of the bond and physical adsorption across the interface due to the high permittivity of water, which is highly relevant to the properties of high polarity adherend [219]. Both the ionic attraction force and van der Waals forces are inversely proportional to the relative permittivity, and a small amount of water will significantly reduce the adhesion strength [220]. Besides, the high surface tension resulting from the high polarity of water will reduce the thermodynamic work of adhesion and thus damage interface durability [221]. Temperature mechanically affects the performance of PV modules through differential thermal expansion, as can be seen in real life because of daily and seasonal thermal cycles. However, it is assumed that the loss of adhesion strength is primarily reduced by moisture ingress, and temperature acts as the 125
accelerator of the effects of humidity [222, 223]. This research aims at understanding the influence of moisture penetration accelerated by temperature on the adhesion durability of interfaces in the PV system. More specifically, the target is to develop a computational framework allowing for the accurate simulation of delamination in PV modules for different ageing conditions. In this regard, research is motivated by the fact that peeling techniques could be exploited to partially disassemble PV laminates at the end of their lifetime. This can be energetically and economically favourable over the actual recycling operations which consider a full crushing of the glass laminate and then separation of its recyclable components aftwerwards. During the past decades, cohesive zone models (CZM) have been widely incorporated into the nonlinear finite element method for engineering applications concerning the formation of free surfaces due to fracture. Compared with other methods such as the virtual crack closure technique (VCCT) [224], CZM can be easily implemented into research and commercial codes and is highly versatile to accommodate different physical phenomena. CZM is a nonlinear traction–relative displacement relationship between the interface tractions across the discontinuity and the sliding and opening displacements. Applications cover many fields, like quasi-static fracture in solids [225, 226], snap-back instabilities during crack propagation [227, 228], crack propagation in composites [126, 229], micromechanical and multiscale analyses [230, 135], and crack assessment in bio-inspired materials [136], among others. However, to the best knowledge of the authors, no previous research considers effect of environmental factors including both moisture and temperature on the interfacial strength in CZM for numerical prediction of fracture especially towards PV lifetime recycling. Peeling of layers of the laminate from the layer of solar cells would require large deformation simulations [231], since the deformed configuration cannot be approximated by the undeformed one. As a result, the interfacial gap calculated according to the initial undeformed configuration can result in errors that cannot be ignored when assessing the total energy required for peeling, which is the quantity of primary interest for 126
recycling. To accurately predict the large deformation framework, both the surface separation and relative rotation between two sides of the interface need to be taken into account. A first attempt has been proposed in [138], which adopted the middle surface of the cohesive element in the deformed configuration as the reference plane to calculate the normal and tangential directions of the interface. Nevertheless, this formulation stemming from differentiation of cohesive tractions with regard to normal unit vector led to non-symmetric geometric stiffness. Based on that research, a 3D interface finite element was used to simulate the standard fracture mechanics tests in thin aluminum panels, in which residual with rotation matrix is updated during the deformation process [232]. However, its linearization does not consider the dependence of differential operator with respect to the kinematic field. Another alternative formulation for large displacement interface element was proposed in [233] by introducing a co-rotational reference system coincident with one of the deformed surfaces of the interface element, but not easy to be implemented into the finite element method due to this co-rational description. In order to account for large deformation of the interface, the interface element formulation for both geometrical and material nonlinearity developed in [127] is herein adopted and extended to 3D application in this work. The first part of this section proposes a computational framework at finite strains for the modeling of EVA polymer-solvent system in PV recycling using the solvent method, which takes into account the coupled solvent diffusion, swelling and shrinking of EVA polymers, as well as effects of mechanical constraints and initial crosslinking of EVA during PV lamination. In Section 4.1.1 and 4.1.2, the basic continuum kinematics in the large-deformation framework and balance laws for force and mass balance are given according to the standard modern continuum mechanics. Section 4.1.3 presents the derivation of thermodynamics and constitutive theory, and the specific free energy definition is given in Section 4.1.4. The numerical implementation of this theory is derived in detail using the finite element method in Section 4.1.5, and then applied to the PV recycling through some examples in Section 4.1.6. In the second part 127
of this section, the polynomial CZM [234] coupled with strength degradation by the humidity dose model [222] is proposed for the peeling simulation in PV recycling, as shown in Section 4.2.1. Then it is incorporated into the 3D large deformation interface element extended from the 2D version proposed in [127], see Section 4.2.2. The numerical results obtained from this coupling strategy correlate with experimental results [222] very well for adhesion strength between backsheet and encapsulant layer in PV. In section 4.2.3, the difference between the specific energy required in peeling and traditional crushing for recycling of PV panels is compared and discussed. 4.1 Coupled modeling of polymer-solvent system in PV recycling 4.1.1 Continuum kinematics According to the standard concepts of continuum mechanics framework [235, 236], the EVA body in PV modules can be described as a macroscopically homogeneous continuum solid within the Euclidean space E3. It is assumed that kinematic fields in the continuum theory are applied at large length scales compared to that associated with the polymer molecule network and pore-structure of EVA [237, 238]. We identify the macroscopically homogeneous EVA with the space it occupies in the reference configuration B0, and describe the current configuration Bby a smooth one-to-one transformation mapping φ:B0→E3. The position vector of an arbitrary material point in the reference configuration B0is denoted by X, and the corresponding vector in the current configuration Bis represented by x=φ(X, t). The displacement field in the reference configuration is defined as u(X, t) = x−X=φ(X, t)−X, and the transformation mapping is described by the deformation gradient F, which reads F=∂φ(X, t) ∂X=∇Xφ(X, t)(4.1) where ∇Xdenotes the material gradient. The determinant of the defor- 128
mation gradient is defined as J=det F= 0, representing the volume change from the reference configuration to the current configuration. The velocity vand its gradient with respect to the position vector xin the current configuration Lare given by v=˙ φ(X, t)and L=∇xv=˙ F·F−1(4.2) where ∇xdenotes the spatial gradient. To describe the coupled mechanical behavior of polymer and solvent, the multiplicative decomposition of deformation gradient is given by F=Fe·Fs(4.3) where Feand Fsrepresent the elastic and swelling parts of deformation gradient F, respectively. The swelling distortion Fsis defined as Fs=λs1,with λs≥0(4.4) where λsis the swelling stretch, which accounts for the swelling of the polymeric body due to absorbed solvents that are pinned to the network structure in the intermediate local space at the material point X. For convenience, the determinant Jcan be written as J=JeJs,with Je=det Feand Js=det Fs= (λs)3(4.5) Remarkably, (Js−1) represents the volume change per unit reference volume of EVA resulting from the swelling behaviour due to the change in the solvent content. Analogous to the theory developed for hydrogels, the swelling constraint is defined as Js= 1 + υCRand ϕ= (1 + υCR)−1(4.6) where υdenotes the volume of a mole of solvent molecules, CRrepresents the number of moles of absorbed solvent molecules by the polymeric body per unit reference volume, and ϕis the polymer volume fraction. Following the standard decomposition of deformation gradient, we have F=R·U=V·Rand Fe=Re·Ue=Ve·Re(4.7) 129
where Rand Reare rotation tensors, Uand Ueare right stretching tensors, and Vand Veare left stretching tensors. Subsequently, the total right and left Cauchy-Green tensors Cand B, and the elastic corresponding tensors Ceand Beare defined as C=U2=FT·Fand B=V2=F·FT(4.8a) Ce=Ue2=FeT·Feand Be=Ve2=Fe·FeT(4.8b) Next, by substituting Eq. (4.3) into Eq. (4.2)2, the velocity tensor Lcan be expressed as L=Le+Fe·Ls·Fe−1(4.9) where the elastic and swelling velocity tensors Leand Lsare given by Le=˙ Fe·Fe−1and Ls=˙ Fs·Fs−1(4.10) 4.1.2 Balance laws In the current configuration, the local form of force balance in the absence of dynamic effects can be expressed as ∇x·σ+b=0(4.11) where bis the body force per unit deformed body, and σis the Cauchy stress tensor deriving from the balance of angular momentum. As is standard [239], the first Piola-Kirchhoff stress tensor Pand Kirchhoff stress tensor τcan be related to the symmetric Cauchy stress tensor σby P=Jσ·F−T(4.12a) τ=Jσ(4.12b) and thus σ=J−1P·FT(4.13) Let ∂Buand ∂Btbe the complementary surfaces of the polymeric body Bwith ∂B=∂Bu∪∂Btand ∂Bu∩∂Bt=∅. The boundary conditions imposed on the body in the current configuration can be given by u=¯ uon ∂Bu(4.14a) 130