THERMOMECANICAL MODELING OF POLYMERIC BATTERY SEPARATORS By Royal Chibuzor Ihuaenyi A DISSERTATION Submitted to Michigan State University in partial fulfillment of the requirements for the degree of Mechanical Engineering - Doctor of Philosophy 2023 ABSTRACT One key challenge in the deployment of future e-mobility systems is to ensure the safe operating condition of high-energy density batteries. Therefore, understanding battery failure mechanisms and reducing safety risks are critical in the design of electrified systems. Although the response of battery materials and systems under various conditions has been extensively explored in recent years, there are still a lot of challenges with developing models for predicting failures. One such challenge is the development of accurate thermomechanical models to predict battery failure caused by combined thermal and mechanical loadings. Such thermomechanical models aim to identify the thermomechanical failure condition of batteries through battery materials such as the separator. The structural integrity of battery separators plays a critical role in battery safety. This is because the deformation and failure of the separator can lead to an internal short circuit which can cause thermal runaway. In thermal runaway scenarios, the separator first expands and then shrinks before reaching its melting temperature. Furthermore, this shrinkage induces tensile stresses in the separator. Hence, developing a thermomechanical model that can predict the response of separators in their entire range of deformation is necessary. Commonly used battery separators are dry-processed polymeric membranes with anisotropic microstructures and deformation modes that involve various physical processes that are difficult to quantify. These complexities introduce challenges in their characterization and modeling as their properties and structural integrity depend on multiple factors such as loading rate, loading direction, temperature, and the presence of an electrolyte. To predict the structural integrity of polymeric separators in abuse scenarios an understanding of the thermal and mechanical behavior of the separator is needed. Due to the multiple factors influencing the structural integrity of polymeric battery separators, developing models for the prediction of their thermomechanical response has always been challenging. Furthermore, computational models in the form of user-defined material models are needed to account for these factors since existing material models in commercial software do not have that capability. In this study, thermomechanical models are developed to predict the response of polymeric battery separators in thermal ramp scenarios. The time-dependent response of polymeric battery separators is taken into account and the material is modeled as viscoelastic in the deformation region before yielding and as viscoplastic under large deformations post-yield. As a first step, a linear thermoviscoelastic model developed on an orthotropic framework was extended to account for the temperature effect and the plasticization effect of electrolyte solutions to predict the thermomechanical response of separators within the linear range of its deformation. In the developed linear orthotropic thermomechanical model, the temperature effect was introduced through the time-temperature superposition principle (TTSP). To account for the plasticization effect of electrolyte solutions on the thermomechanical response of the separator, a time- temperature-solvent superposition method (TTSSM) was developed to model the behavior of the separator in electrolyte solutions based on the viscoelastic framework established in air. Furthermore, an orthotropic nonlinear thermoviscoelastic was developed to predict the material response under large deformations before the onset of yielding. The model was developed based on the Schapery nonlinear viscoelastic model and a discretization algorithm was employed to evaluate the nonlinear viscoelastic hereditary integral with a kernel of Prony series based on a generalized Maxwell model with nonlinear springs and dashpots. Temperature dependence was introduced into the model through the TTSP. Subsequently, the developed nonlinear viscoelastic model was coupled with a viscoplastic model developed on the basis of a rheological framework that considers the mechanisms involved in the initial yielding, change in viscosity, strain softening and strain hardening in the stress-strain response of polymeric battery separators. The coupled viscoelastic – viscoplastic model was developed to predict the thermomechanical response of separators in their entire range of deformation before the onset of failure. The material investigated in this work is Celgard®2400, a porous polypropylene (PP) separator. Experimental procedures were carried out under different loading and environmental conditions, using a dynamic mechanical analyzer (DMA), to characterize the material response, calibrate and validate the developed models. The developed thermomechanical models were implemented as user-defined subroutines in LS-DYNA® finite element (FE) package, which enables simulations with the thermal expansion/shrinkage behavior. Furthermore, analytical solutions were developed to verify the implementation and predictions of the viscoelastic models. The results from this study show that the model predictions of the material anisotropy, rate dependence, temperature dependence, and plasticization effect of electrolyte solutions agree reasonably well with the experimental data. The results also demonstrate that the non-isothermal simulations without considering the thermal expansion/shrinkage behavior of the separator resulted in large errors. ACKNOWLEDGEMENTS I would like to express my sincerest gratitude to my advisor, Professor Xinran Xiao, for her continuous support and guidance during my doctoral studies at the Composite Vehicle Research Center at MSU. I am so fortunate to be one of her last PhD students in her very successful academic career. In these four years, she worked very closely with me while giving me the freedom to run with my ideas. She has taught me how to carry out good scientific research, manage research projects, the basics of academic writing, communication skills and how to maintain a healthy work-life balance. I truly appreciate her input in my journey as these attributes I have been taught will become great treasures in my future academic career. I want to thank my Ph.D. thesis committee members, Professor Thomas Pence, Professor Andre Lee, and Professor Seungik Baek, for their constructive comments and advice during my doctoral studies. I want to thank Professor Thomas Pence, who instructed me on engineering analysis, viscoelasticity and thermoelasticity. Much of what I learnt from his classes was applied in my research endeavors and in the future, I hope to impact students in the way he has impacted me through his pedagogy. Additionally, special thanks to our collaborators Dr. Chulheung Bae and Dr. Jie Deng at Ford for their support, insightful comments and suggestions on our project during my doctoral studies. Working and collaborating these past four years has been a pleasure. I am also grateful to my labmates Dr. Shutian Yan and Sakib Iqbal for their collaboration, support and suggestions during my research endeavors. Many thanks to Shutian specifically as she helped me make a smooth transition into the lab and at the early stages of my research endeavors at MSU. I am also grateful to all my friends who helped me during my stay at MSU, especially Alessandro Bo and Lincoln Mtemeri. I want to thank Aaron Walworth at the MSU School of Packaging Research Labs for all his efforts in making it possible for me to complete my experiments. I also appreciate the training from Ed Drown and Dr. Per Askeland at the CMSC Labs. I am very thankful to the Mechanical Engineering Department, College of Engineering, and the Graduate School at MSU for their funding support. I want to express my deepest gratitude to my parents Sampson Ihuaenyi and Patience Ihuaenyi, for their unconditional love and support. They have imbibed the spirit of hard work, courage and steadfastness in me since I was a child. Although I am agnostic, I also want to thank them for their prayers (they seem to be working). I want to also thank my siblings for their encouragement and for always being there for me throughout this journey. I am also thankful to iv my in-laws, Igor Dolgun and Olga Dolgun and their entire family for their invaluable support and constant encouragement throughout this period. Finally, thanks are due to my dear wife Alexandra Dolgun for her understanding, support and belief in my dreams as much as hers. I am grateful to her for giving me the time and space for my academic work. She has been with me throughout this journey and I thank her for making my life a lot easier. This work is supported by the Ford/MSU alliance program. v TABLE OF CONTENTS Chapter 1 Introduction .................................................................................................................... 1 1.1 Lithium-ion Battery Separators............................................................................................. 1 1.2 Problem Definition................................................................................................................ 3 1.3 Research Objectives .............................................................................................................. 4 1.4 Thesis Outline ....................................................................................................................... 5 REFERENCES ........................................................................................................................... 7 Chapter 2 Literature review ............................................................................................................ 9 2.1 Linear Viscoelasticity ........................................................................................................... 9 2.1.1 Overview ........................................................................................................................ 9 2.1.2 Boltzmann Superposition Principle ............................................................................. 10 2.1.3 Rheological Models ..................................................................................................... 11 2.1.4 Elastic – Viscoelastic Correspondence Principle ......................................................... 14 2.2 Nonlinear Viscoelasticity .................................................................................................... 16 2.2.1 Overview ...................................................................................................................... 16 2.2.2 Schapery Nonlinear Viscoelastic Model ...................................................................... 16 2.2.3 Other Nonlinear Viscoelastic Models .......................................................................... 17 2.3 Environmental Effects and Long-Term Material Response ............................................... 19 2.3.1 Time-Temperature Superposition Principle (TTSP) .................................................... 19 2.3.2 Time-Temperature-Moisture Superposition Method (TTMSM) ................................. 21 2.4 Viscoplastic Modeling of Polymeric Materials .................................................................. 22 2.4.1 Overview ...................................................................................................................... 22 2.4.2 Viscoplastic Constitutive Models ................................................................................ 23 2.5 Thermomechanical Modeling of Polymeric Battery Separators ......................................... 26 2.5.1 Thermomechanical Behavior of Polymeric Battery Separators ................................... 26 2.5.2 Plasticization Effects of Electrolyte Solutions ............................................................. 28 2.5.3 Existing Models for Polymeric Battery Separators...................................................... 28 REFERENCES ......................................................................................................................... 30 Chapter 3 Orthotropic Linear Thermoviscoelastic Constitutive Modeling .................................. 39 3.1 Orthotropic Linear Viscoelastic Model Development ........................................................ 39 3.1.1 Constitutive Model Overview ...................................................................................... 39 3.1.2 Numerical Intergation of the Linear Viscoelastic Hereditary Integral with a Kernel of Prony Series .......................................................................................................................... 40 3.1.3 Plane Stress Stiffness Matrix ....................................................................................... 42 3.1.4 Three-Dimensional Stiffness Matrix with Transverse Orthotropy .............................. 43 3.2 Orthotropic Linear Thermoviscoelastic Material Characterization .................................... 44 3.2.1 Experimental Procedure ............................................................................................... 44 3.2.2 Determination of Linear Viscoelastic Strain Limit ...................................................... 46 3.2.3 Stress relaxation at different temperatures ................................................................... 49 3.2.4 Determination of In-Plane Shear Relaxation Modulus ................................................ 50 3.3 Temperature and Electrolyte Effects .................................................................................. 51 3.3.1 Implementation of the Time-temperature superposition principle (TTSP) ................. 51 3.3.2 Time-temperature-solvent superposition method (TTSSM) ........................................ 58 3.4 Analytical Solutions for Stress Relaxation ......................................................................... 62 3.4.1 Single Step Stress Relaxation at Constant Temperature .............................................. 62 vi 3.4.2 Step loading at constant temperature ........................................................................... 63 3.4.3 Step temperature history at constant strain level ......................................................... 64 3.5 Prony Series Fitting and Summary of TTSSM Parameters ................................................ 65 3.6 Simulation Details ............................................................................................................... 67 3.6.1 Shell Element Simulations ........................................................................................... 67 3.6.2 Hexahedral Solid Element Simulations ....................................................................... 68 3.7 Orthotropic Linear Thermoviscoelastic Model Verification .............................................. 69 3.8 Orthotropic Linear Thermoviscoelastic Model Validation ................................................. 77 3.9 Summary ............................................................................................................................. 82 REFERENCES ......................................................................................................................... 84 Chapter 4 Orthotropic Nonlinear Thermoviscoelastic Constitutive Modeling ............................. 86 4.1 Orthotropic Nonlinear Viscoelastic Model Development .................................................. 86 4.1.1 Constitutive Model Overview ...................................................................................... 86 4.1.2 Numerical Integration of the Nonlinear Viscoelastic Hereditary Integral with a Kernel of Prony Series ...................................................................................................................... 88 4.1.3 Modifications to the Schapery Nonlinear Viscoelastic Model .................................... 91 4.1.4 Plane Stress Stiffness Matrix ....................................................................................... 91 4.2 Orthotropic Nonlinear Thermoviscoelastic Material Characterization ............................... 92 4.2.1 Experimental Procedure ............................................................................................... 92 4.2.2 Stress Relaxation at Different Strain Levels ................................................................ 94 4.2.3 Determination of In-Plane Shear Relaxation Modulus ................................................ 95 4.3 Parameter Identification ...................................................................................................... 97 4.3.1 Overview ...................................................................................................................... 97 4.3.2 Determination of Temperature-Dependent Model Parameters .................................... 97 4.3.3 Prony Series Fitting .................................................................................................... 103 4.3.4 Determination of Strain-Dependent Nonlinear Parameters ....................................... 104 4.3.5 Calibration of Instrument Expansion ......................................................................... 107 4.4 Analytical Solutions for Stress Relaxation ....................................................................... 108 4.4.1 Single Step Stress Relaxation at Constant Temperature ............................................ 109 4.4.2 Step loading at constant temperature ......................................................................... 110 4.4.3 Step temperature history at constant strain level ....................................................... 110 4.5 Orthotropic Nonlinear Thermoviscoelastic Model Verification ....................................... 111 4.6 Orthotropic Nonlinear Thermoviscoelastic Model Validation ......................................... 115 4.7 Summary ........................................................................................................................... 120 REFERENCES ....................................................................................................................... 121 Chapter 5 Coupled Viscoelastic – Viscoplastic Modeling ......................................................... 124 5.1 Viscoplastic Model Overview........................................................................................... 124 5.2 Uniaxial Tensile Testing ................................................................................................... 126 5.3 Identification of Strain Rate-Dependent Model Parameters at Constant Temperature .... 129 5.4 Introducing Temperature Effect ........................................................................................ 134 5.5 Coupled Viscoelastic – Viscoplastic Model Implementation and Validation .................. 138 5.6 Summary ........................................................................................................................... 141 REFERENCES ....................................................................................................................... 143 Chapter 6 Conclusion and Future Work ..................................................................................... 144 6.1 Conclusion ........................................................................................................................ 144 vii 6.1.1 Orthotropic Linear Thermoviscoelastic Modeling ..................................................... 144 6.1.2 Orthotropic Nonlinear Thermoviscoelastic Modeling ............................................... 145 6.1.3 Coupled Viscoelastic – Viscoplastic Modeling ......................................................... 145 6.2 Future Work ...................................................................................................................... 146 6.2.1 Model Extension and Adaptation ............................................................................... 146 6.2.2 Model Implementation ............................................................................................... 147 REFERENCES ....................................................................................................................... 148 APPENDIX A. SUPPLEMENTARY EXPERIMENTAL RESULTS ....................................... 149 APPENDIX B. MATLAB SCRIPT FOR DETERMINATION OF PRONY SERIES PARAMETERS .......................................................................................................................... 152 APPENDIX C. PYTHON SCRIPT FOR NON-NEGATIVE LEAST SQUARE CURVE FITTING ..................................................................................................................................... 154 APPENDIX D. PUBLICATIONS AND PRESENTATIONS ................................................... 156 viii Chapter 1 Introduction This chapter presents an overview of battery separators, focusing on dry-processed polymeric microporous membranous separators used in lithium-ion batteries. The problem to be addressed in this thesis is defined and the research objectives and tasks are delineated. Finally, an outline, detailing the organization of this thesis is presented in this chapter. 1.1 Lithium-ion Battery Separators In lithium-ion batteries (LIBs), the separator is a porous membrane that is placed between electrodes of opposite polarity to provide electrical insulation while allowing ionic transport between the electrodes [1]. In the working condition of a LIB, the separator is immersed in an electrolyte through which ions flow. Hence, the separator must be chemically and electrochemically stable, and mechanically strong. Fig.1.1 shows a typical internal structure of a LIB cell comprising the anode, cathode and separator immersed in an electrolyte. According to composition and structure, separators are categorized into three major groups: (1) microporous polymer membranes [5], non-woven fabric mats [6], and inorganic composite separators [7]. However, polymeric membranes of 20-30 mm thickness are widely used in commercial LIB. This is due to their relatively low cost, thermal shutdown properties, and their provision of higher energy and power densities due to their small relative thickness [4]. The most commonly used polymeric membranes are semi-crystalline polyolefin materials such as polyethylene (PE), polypropylene (PP), their blend (PP-PE), and high-density polyethylene (HDPE). Dry-processed polymeric LIB separators show an orthotropic mechanical response, meaning that the material has three planes of symmetry and three corresponding mutually orthogonal axes [8]. These are the machine direction (MD), the transverse direction (TD), and the through-thickness direction (TTD). This anisotropy is introduced in the manufacturing process which involves annealing, stretching along the MD, and heat fixation [1,4]. 1 Figure 1.1 Typical internal schematic of a lithium-ion battery cell [9]. Furthermore, the surface microstructure of a typical polymeric separator is depicted in Fig.1.2(a), with its two in-plane material orientations (MD and TD). The 3D representation of a single-layer PP separator showing the three material orientations referred to as MD, TD and TTD is presented in Fig.1.2(b). The structural views of the separator provide a definition of its pore size, shape and structure. From Fig.1.2, the pores that allow for ionic transport are splits. These splits are a fraction of a micrometer long and ten to hundred nanometers apart from each other. The fiber-like structures between the splits are amorphous and the thick regions are semi-crystalline. 2 (a) (b) Figure 1.2 (a) Surface microstructure of Celgard®2400 [1] (b) 3D representation of Celgard®2500 [10]. 1.2 Problem Definition The structural integrity of the separator is critical in preventing an internal short circuit (ISC). Separators in LIBs may fail under abusive mechanical loadings such as impact and crash, be pierced by dendrites grown during excessive electrical over-charge or discharge, or melt due to local overheating initiated by manufacturing defects [2,3,11]. Failure of the separator can also be induced by combined thermomechanical load caused by local overheating associated with non- uniform heat generation and dissipation as the battery ages. Failure of the separator can cause the electrodes to come in contact with each other and lead to an ISC, which in turn leads to local temperature increases (i.e., thermal runaway), and can cause a thermal event. Furthermore, the relationship between separator failure and the onset of ISC has been studied and established using experimental techniques and finite element analysis (FEA) [11-13]. The evolution of ISC in LIBs can cause thermal runaway. To prevent ISC and enhance the safety of LIBs it is desirable to predict the response of separators under combined mechanical and thermal loadings. For this purpose, thermomechanical models that can accurately describe the orthotropic thermomechanical response of the polymeric battery separator are needed. Dry-processed polymeric separators are strongly anisotropic [14-20], 3 rate-dependent [15,16], and temperature-dependent [14,15] with a distinctive thermal expansion/shrinkage behavior [21,22]. Furthermore, polymeric separators soften when immersed in electrolyte solutions [15-17]. This plasticization effect of the electrolyte solution on the separator is a factor that needs to be taken into account. To accurately predict the thermomechanical response of polymeric battery separators to prevent ISC, it is vital to account for the above-mentioned constitutive behaviors in material models. Without taking these factors into account, developed material models will be inaccurate in predicting the mechanical response of separators in LIBs. 1.3 Research Objectives The objective of this research is to develop thermomechanical models for polymeric battery separators for simulations with combined thermal and mechanical loadings in crashworthiness analysis. Like typical polymeric materials, the response of polymeric battery separators to mechanical deformation is time and temperature dependent. Ideally, polymeric battery separators should be modeled as viscoelastic materials in the deformation region up to yielding and viscoplastic at the onset of permanent deformation and beyond. Hence, an orthotropic linear viscoelastic, orthotropic nonlinear viscoelastic and a coupled viscoelastic-viscoplastic thermomechanical model are developed in this work. The developed linear viscoelastic model is based on the framework of the well-established theory of viscoelasticity for isotropic materials [22-26] and counts for material anisotropy. In addition to considering the material anisotropy, the model will account for rate dependency, temperature effect, thermal expansion/shrinkage behavior of the polymeric separator, and the electrolyte effect. These are common factors that influence the mechanical and thermal response of polymeric separators in LIBs as discussed in section 1.2. To carry out this objective, the first task is to extend the orthotropic linear viscoelastic model developed by Yan et al [27], to account for temperature dependency and electrolyte effect. In this work, the temperature dependency will be introduced through the time-temperature superposition principle (TTSP) of linear viscoelasticity [22-25]. Furthermore, the electrolyte effect will be introduced through a time- temperature-solvent superposition method developed in the current work. The TTSSM uses the TTSP in air as a framework to predict the material behavior in electrolyte solutions at higher temperatures. This is in contrast to the time-temperature-moisture superposition method (TTMSM) [28-31] which superimposes the time-temperature master curves for different moisture 4 content (humidity) values to form a super time-temperature-moisture master curve to describe the polymer behavior at different temperatures and moisture content levels [28]. The second task is to develop an orthotropic nonlinear thermoviscoelastic model to account for the material response in the large deformation region up to yielding or the evolution of irrecoverable deformations. The temperature effect will also be introduced into the orthotropic nonlinear thermomechanical model through the TTSP. The model will be able to predict the anisotropic, rate-dependent, and temperature-dependent response of polymeric separators. This model will also take the thermal expansion behavior of the separator into account. Furthermore, the third task will involve coupling the orthotropic nonlinear thermoviscoelastic model with a viscoplastic model to account for the material response in the entire range of its deformation before failure. A yield criterion will be introduced into the model to establish the onset of permanent deformations and act as a coupling mechanism between the nonlinear viscoelastic and viscoplastic models. Finally, a viscoplastic formulation will be introduced to account for the rate and temperature-dependent stress-strain response of the polymeric separator once the yield criterion is satisfied. Experimental methods will be carried out for the determination of model parameters and validation of the developed model predictions. The model parameters will be determined for the response of a selected polymeric separator (Celgard®2400) in the MD, TD, and in-plane shear. Celgard®2400 is a semi-crystalline polypropylene (PP) separator of thickness 25𝜇𝑚. The developed model will be implemented in LS-DYNA® finite element (FE) package as a user- defined subroutine. The model predictions of the material anisotropy, rate dependence, temperature dependence, and electrolyte effect will be validated against experimental data. 1.4 Thesis Outline This document is organized as follows. Chapter 1 introduces the problem and presents an overview of Lithium-ion battery separators and pays close attention to polymeric battery separators. The research objectives for this work are also delineated in this chapter. Chapter 2 provides a literature review of linear and nonlinear viscoelasticity, viscoplasticity and constitutive modeling of polymeric materials. Also, a review is presented on the thermomechanical behavior of polymeric separators, the plasticization effect of electrolyte 5 solutions and available models for polymeric battery separators. This chapter also provides a review of the orthotropic modeling of viscoelastic materials. Chapter 3 presents the development of an orthotropic linear thermoviscoelastic model for polymeric battery separators. This chapter discusses the experimental methodology for material characterization and model parameter determination within the linear viscoelastic deformation limit. Also, the developed model is implemented, verified using analytical solutions and validated in this chapter. Chapter 4 presents the development of an orthotropic nonlinear thermoviscoelastic model for polymeric battery separators. A discretization algorithm of the nonlinear Schapery hereditary integral is presented. Methods for model parameter identification are discussed and the developed model is implemented, verified using analytical solutions and validated against experimental data. Chapter 5 presents the development, implementation and validation of a coupled viscoelastic-viscoplastic model for polymeric battery separators. Chapter 6 provides a summary and conclusion of this thesis, as well as future research work. 6 REFERENCES 1. P. Arora and Z. Zhang. Battery separators. Chemical reviews, 104(10), pp.4419-4462, 2004. 2. B. Liu, Y. Jia, C. Yuan, L. Wang, X. Gao, S. Yin, and J. Xu. Safety issues and mechanisms of lithium-ion battery cell upon mechanical abusive loading: A review. Energy Storage Materials, 24, pp.85-112, 2020. 3. J. Zhang, L. Lai, H. Wang, M. Chen, and Z. X. Shen. Energy storage mechanisms of anode materials for potassium ion batteries. Materials Today Energy, p.100747, 2021. 4. S. S. Zhang. A review on the separators of liquid electrolyte Li-ion batteries. Journal of power sources, 164(1), pp.351-364, 2007. 5. T. Ishihara, S. Miyaoka, K. Kono, D.J. Crowther, P. Brant, K. Yamada, in, Google Patents, 2017. 6. R. Butin, J. Keller, J. Harding, in, Google Patents, 1974. 7. S. S. Zhang, K. Xu, and T.R. Jow. An inorganic composite membrane as the separator of Li- ion batteries. Journal of Power Sources, 140(2), pp.361-364, 2005. 8. I.M. Daniel, O. Ishai, Engineering Mechanics of Composite Materials, Oxford University Press, 2006. 9. C. Liu, Z. G. Neale, and G. Cao. Understanding electrochemical potentials of cathode materials in rechargeable batteries. Materials Today, 19(2), pp.109-123, 2016. 10. T. Sarada, L. C. Sawyer, and M. I. Ostler. Three-dimensional structure of celgard® microporous membranes. Journal of Membrane Science, 15(1), pp.97-113, 1983. 11. B. Liu, X. Duan, C Yuan, L. Wang, J. Li, D. P. Finegan, B. Feng, and J. Xu. Quantifying and modeling of stress-driven short-circuits in lithium-ion batteries in electrified vehicles. Journal of Materials Chemistry A, 9(11), pp.7102-7113, 2021. 12. C. Yuan, L. Wang, S. Yin, and J. Xu. Generalized separator failure criteria for internal short circuit of lithium-ion battery. Journal of Power Sources, 467, p.228360, 2020. 13. D. C. Lee and C. W. Kim. Detailed Layered Nonlinear Finite Element Analysis for Lithium- Ion Battery Cells to Predict Internal Short Circuits Due to Separator Fractures under Hemisphere Indentation. Journal of The Electrochemical Society, 167(12), p.120511, 2020. 14. J. Zhu, T. Wierzbicki and W. Li. A review of safety-focused mechanical modeling of commercial lithium-ion batteries. Journal of Power Sources, 378, pp.153-168, 2018. 15. Avdeev, M. Martinsen, and A. Francis. Rate-and temperature-dependent material behavior of a multilayer polymer battery separator. Journal of materials engineering and performance, 23(1), pp.315-325, 2014. 16. J. Xu, L. Wang, J. Guan and S. Yin. Coupled effect of strain rate and solvent on dynamic mechanical behaviors of separators in lithium-ion batteries. Materials & Design, 95, pp.319- 328, 2016. 17. Sheidaei, X. Xiao, X. Huang and J. Hitt. Mechanical behavior of a battery separator in electrolyte solutions. Journal of Power Sources, 196(20), pp.8728-8734, 2011. 7 18. J. Cannarella, X. Liu, C. Z. Leng, P. D. Sinko, G. Y. Gor and C. B. Arnold. Mechanical properties of a battery separator under compression and tension. Journal of the Electrochemical Society, 161(11), p.F3117, 2014. 19. X. Zhang, E. Sahraei and K. Wang. Deformation and failure characteristics of four types of lithium-ion battery separators. Journal of Power Sources, 327, pp.693-701, 2016. 20. S. Yan, J. Deng, C. Bae, Y. He, A. Asta and X. Xiao. In-plane orthotropic property characterization of a polymeric battery separator. Polymer Testing, 72, pp.46-54, 2018. 21. C. T. Love. Thermomechanical analysis and durability of commercial micro-porous polymer Li-ion battery separators. Journal of Power Sources, 196(5), pp.2905-2912, 2011. 22. S. Yan, J. Deng, C. Bae, and X. Xiao. Thermal expansion/shrinkage measurement of battery separators using a dynamic mechanical analyzer. Polymer Testing, 71, pp.65-71, 2018. 23. J. D. Ferry, Viscoelastic Properties of Polymers, John Wiley & Sons, 1980. 24. I.M. Ward, J. Sweeney, An Introduction to the Mechanical Properties of Solid Polymers, John Wiley & Sons, 2004. 25. S. Wineman, K. R. Rajagopal, Mechanical response of polymers: an introduction. Cambridge university press, 2000. 26. H. F. Brinson and L. C. Brinson. Polymer engineering science and viscoelasticity. An introduction, Springer US, 2008. 27. R. M. Christensen. Theory of Viscoelasticity: An Introduction. Academic Press, New York, second edition, 1982. 28. S. Yan, J. Deng, C. Bae, S. Kalnaus and X. Xiao. Orthotropic viscoelastic modeling of polymeric battery separator. Journal of The Electrochemical Society, 167(9), p.090530, 2020. 29. P. C. Suarez-Martinez, P. Batys, M. Sammalkorpi and J. L. Lutkenhaus. Time–temperature and time–water superposition principles applied to poly (allylamine)/poly (acrylic acid) complexes. Macromolecules, 52(8), pp.3066-3074, 2019. 30. V. Fabre, G. Quandalle, N. Billon and S. Cantournet. Time-Temperature-Water content equivalence on dynamic mechanical response of polyamide 6, 6. Polymer, 137, pp.22-29, 2018. 31. A. Ishisaka, and M. Kawagoe. Examination of the time-water content superposition on the dynamic viscoelasticity of moistened polyamide 6 and epoxy. Journal of Appled Polymer Science, 93, 560–567, 2004. 32. S. M. Zhou, K. Tashiro, and T. Ii. Confirmation of universality of time–humidity superposition principle for various water‐absorbable polymers through dynamic viscoelastic measurements under controlled conditions of relative humidity and temperature. Journal of Polymer Science Part B: Polymer Physics, 39(14), pp.1638-1650, 2001. 8 Chapter 2 Literature review This chapter presents a review of the viscoelastic theory as well as constitutive models for polymeric materials based on viscoelastic models. The review on the viscoelastic theory pays attention to both the linear and nonlinear viscoelastic theories for the time-dependent response of polymeric materials. Furthermore, a review of the environmental effects and long-term response of polymeric materials is presented. Literature studies on the viscoplastic modeling of polymeric materials and the thermomechanical modeling of polymeric battery separators are also presented. 2.1 Linear Viscoelasticity 2.1.1 Overview Depending on the type of loading and environmental conditions, polymeric materials can display a wide range of mechanical behaviors, including those of an elastic solid and a viscous liquid. Polymers are typically referred to as viscoelastic, or time- and temperature-dependent, due to this mix of behaviors. As a result of this response, polymeric materials are susceptible to creep when a constant load is applied to the material. During the time at which the loading is kept constant, following an initial linear elastic response, strain continues to accumulate in the material. The typical creep response of viscoelastic materials is shown in Fig.2.1 below. Figure 2.1 Creep response for a constant applied stress. The creep response of viscoelastic materials for a constant applied loading can be expressed as: 𝜀(𝑡) = 𝐽(𝑡) ∙ 𝜎0 (2.1) 9 where 𝐽(𝑡) is the creep compliance function, 𝜀(𝑡) is the total creep strain, and is 𝜎0 the constant applied stress. Viscoelastic materials also exhibit another time-dependent response known as stress relaxation. When a constant strain is applied to a viscoelastic material, the stress in the material relaxes for the entire duration of the time during which the strain is kept constant. Fig.2.2 presents the typical stress relaxation response of viscoelastic materials. Figure 2.2 Stress relaxation response for a constant applied strain. The stress relaxation response of viscoelastic materials for a constant applied strain can be expressed as: 𝜎(𝑡) = 𝐺(𝑡) ∙ 𝜀0 (2.2) where 𝐺(𝑡) is the stress relaxation modulus, 𝜎(𝑡) is the total stress, and is 𝜀0 the constant applied strain. Furthermore, for the response of a material to be considered linear viscoelastic, the stress relaxation modulus or creep compliance function must be independent of the applied loading level [1]. Also, the creep and relaxation functions must be separable. 2.1.2 Boltzmann Superposition Principle The Boltzmann superposition principle is one of the fundamental mathematical representations of linear viscoelastic response. Boltzmann [2], proposed that the creep or stress relaxation response of a viscoelastic material is a function of its entire past loading history. This simply means that each loading step independently contributes to the final deformation. Hence, the total deformation can be determined by simply adding all of the individual contributions of the loading steps. 10 From the Boltzmann superposition principle, the total stress at time t for the stress relaxation response to a multistep loading program in which the incremental strains ∆𝜀1 , ∆𝜀2 , ∆𝜀3 , etc. are added at times 𝜏1 , 𝜏2 , 𝜏3 , etc. is given by: 𝜎(𝑡) = ∆𝜀1 ∙ 𝐺(𝑡 − 𝜏1 ) + ∆𝜀2 ∙ 𝐺(𝑡 − 𝜏2 ) + ∆𝜀3 ∙ 𝐺(𝑡 − 𝜏3 ) + ∙∙∙∙∙ (2.3) Based on the same principle, the total strain at time t for the creep response to a multistep loading program in which the incremental stresses ∆𝜎1 , ∆𝜎2 , ∆𝜎3 , etc. are added at times 𝜏1 , 𝜏2 , 𝜏3 , etc. is given by: 𝜀(𝑡) = ∆𝜎1 ∙ 𝐽(𝑡 − 𝜏1 ) + ∆𝜎2 ∙ 𝐽(𝑡 − 𝜏2 ) + ∆𝜎3 ∙ 𝐽(𝑡 − 𝜏3 ) + ∙∙∙∙∙ (2.4) Fig.2.3 below presents the stress relaxation response to a multistep loading program. Figure 2.3 Stress relaxation response of a linear viscoelastic material to multistep loading. Furthermore, the summations in the form of Eqns. 2.3 and 2.4 may be generalized into the hereditary integral representation of linear viscoelastic behavior. This form is also known as the Boltzmann superposition integral and can be expressed in the stiffness form as: 𝑡 𝑑𝜀 𝜎(𝑡) = ∫ 𝐺(𝑡 − 𝜏) 𝑑𝜏 (2.5) 0 𝑑𝜏 In terms of creep compliance, the hereditary integral takes the form: 𝑡 𝑑𝜎 𝜀 (𝑡 ) = ∫ 𝐽 (𝑡 − 𝜏 ) 𝑑𝜏 (2.6) 0 𝑑𝜏 Also, for a material to be considered as linearly viscoelastic, the Boltzmann superposition principle must be applicable. 2.1.3 Rheological Models Rheological models or mechanical analogs are used to establish differential equations that describe the stress – strain relationship in viscoelastic materials. These models combine springs 11 and dashpots to represent the dual nature of the material response. For linear viscoelastic response, combinations of Hookean springs and Newtonian dashpots are used to introduce elastic and viscous responses respectively. Also, for viscoelastic analysis and developing response equations, these mechanical analogs are assumed to be massless. The Maxwell Model The Maxwell model is one of the fundamental mechanical analogs. In this model, the spring and dashpot are connected in series as shown in Fig.2.4. Figure 2.4 Schematic representation of the Maxwell model. The elastic response of the spring is governed by Hooke’s law expressed in the form: 𝜎=𝐸∙𝜀 (2.7) where 𝐸 is the elastic modulus of the material. The viscous response portrayed by the dashpot is governed by Newton’s law of viscosity: 𝑑𝜀 𝜎=𝜇∙ (2.8) 𝑑𝑡 where 𝜇 is the viscosity of the material. Furthermore, the stress – strain relationship for the Maxwell model is governed by the differential equation expressed as: 𝑑𝜀 1 𝑑𝜎 𝜎 = ∙ + (2.9) 𝑑𝑡 𝐸 𝑑𝑡 𝜇 Consequently, the Maxwell model has certain deficiencies pertaining to the characterization of the complex response of viscoelastic materials. It fails to represent the complexity of creep behavior and a single exponential decay term from the Maxwell model is insufficient in representing stress relaxation behavior [3]. Kevin-Voigt Model 12 The Kevin-Voigt model is another basic mechanical analog in which a Hookean spring is connected to a Newtonian dashpot in parallel as shown in Fig. 2.5. Figure 2.5 Schematic representation of the Kevin-Voigt model. The stress – strain equation for this model takes the following form: 𝑑𝜀 𝜎 = 𝐸∙𝜀+𝜇∙ (2.10) 𝑑𝑡 While this model provides a good representation of the creep response of viscoelastic materials, it does not represent stress relaxation behavior quite well [3]. To overcome the limitations of the Maxwell and Kevin-Voigt models, different combinations of Maxwell and Kevin-Voigt elements can be used to better represent the response of viscoelastic materials. An example of such models is the standard linear solid model [3,4]. More realistic representations of viscoelastic response can also be achieved by combining different Maxwell elements in parallel or different Kevin-Voigt elements in series. Generalized Maxwell Model (N-Maxwell element in parallel) This is one of the standard generalizations of mechanical analogs used to represent the response of viscoelastic materials. In this model, Maxwell elements are arranged in parallel as shown in Fig.2.6. Figure 2.6 Generalized Maxwell Model. 13 The differential equation that expresses the stress – strain response described by the generalized Maxwell model is expressed as: (𝑝0 + 𝑝1 𝐷 + 𝑝2 𝐷2 + 𝑝3 𝐷3 + ∙ ∙ ∙ +𝑝𝑛 𝐷𝑛 )𝜎 = (𝑞0 + 𝑞1 𝐷 + 𝑞2 𝐷2 + 𝑞3 𝐷3 + ∙ ∙ ∙ +𝑞𝑛 𝐷𝑛 )𝜀. (2.11) where 𝑝𝑖 and 𝑞𝑖 are parameters dependent on the elastic modulus of the spring and viscosity of the dashpot, and 𝐷𝑖 is a differential operator for 𝑖 = 0,1,2, . . , 𝑛. 𝑛 is the total number of maxwell elements. Complete details on the mathematical evaluation of the generalized Maxwell model can be found in [4]. Furthermore, the stress relaxation response of the model is given as: 𝑡 𝑡 𝑡 𝑡 −𝜏 −𝜏 −𝜏 −𝜏 𝐺(𝑡) = 𝐺∞ + 𝐺1 ∙ 𝑒 𝑅1 + 𝐺2 ∙ 𝑒 𝑅2 + 𝐺3 ∙ 𝑒 𝑅3 + ∙ ∙ ∙ +𝐺𝑛 ∙ 𝑒 𝑅𝑛 (2.12) where 𝐺𝑖 are relaxation constants and 𝜏𝑅𝑖 is the relaxation time. The form of the stress relaxation function represented as Eqn.2.12 is also known as a Prony series. Another generalized model used to characterize viscoelastic response is the generalized Kevin-Voigt model and more details on this model can be found in [4]. In summary, these mechanical analogs are very useful for representing the basic response of viscoelastic materials. They can also be used as a foundation for introducing constitutive theories for modeling viscoelastic response. 2.1.4 Elastic – Viscoelastic Correspondence Principle The elastic-viscoelastic correspondence principle is a powerful tool used in simplifying the procedure for obtaining solutions in viscoelastic analysis. The correspondence principle applies to problems of a statistically determinate nature and involves viscoelastic bodies subjected to boundary conditions, applied initially and held constant. According to the correspondence principle, the stresses in such bodies can be obtained from the corresponding elastic solutions for the same bodies subjected to the same boundary conditions [4-6]. These elastic solutions can be converted to the appropriate viscoelastic solutions through integral transforms (i.e., Laplace transform). As a result, it can also be said that the equations of the related boundary value problem for an elastic body are defined by the Laplace transforms of the governing equations of the motion of a viscoelastic body. An illustration of an elastic body with its boundary conditions and its corresponding viscoelastic body is presented in Fig.2.7. 14 Figure 2.7 Correspondence of an elastic and a viscoelastic body in terms of time and transform variables. (𝑒) (𝑒) (𝑒) Let 𝑢𝑖 (𝒙), 𝜀𝑖𝑗 (𝒙), and 𝜎𝑖𝑗 (𝒙) which are the displacement, strain and stress field variables respectively represent the solution for the boundary value problem of the elastic body. The elastic-viscoelastic principle implies that the corresponding viscoelastic solutions are obtained by carrying out Laplace transforms of the elastic solutions. The viscoelastic solution is hence obtained by making the following substitutions [4]: 1. Material properties 𝐸𝑖𝑗 → 𝑎𝐺𝑖𝑗 (𝑎), 𝑣 → 𝑎𝑣(𝑎), 𝜇 → 𝑎𝜇𝑅 (𝑎) 2. Field variables (𝑒) (𝑒) (𝑒) 𝑢𝑖 (𝒙) → 𝑢𝑖 (𝒙, 𝑎), 𝜀𝑖𝑗 (𝒙) → 𝜀𝑖𝑗 (𝒙, 𝑎), 𝜎𝑖𝑗 (𝒙) → 𝜎𝑖𝑗 (𝒙, 𝑎) 3. Boundary conditions 𝑢𝑖∗ (𝒙) → 𝑢𝑖∗ (𝒙, 𝑎), 𝑇𝑖∗ (𝒙) → 𝑇𝑖∗ (𝒙, 𝑎) where 𝐸𝑖𝑗 is the elastic modulus, 𝐺𝑖𝑗 is the stress relaxation modulus, 𝑣 is the Poisson’s ratio, 𝜇 is viscosity, and 𝑇𝑖∗ is the surface traction vector. However, not all viscoelastic problems fit into the class of problems solvable by this principle. Boundary value problems for which the boundary conditions change with time at a point on the surface of the body cannot be solved using the elastic–viscoelastic correspondence principle. Furthermore, according to Wineman and Rajagopal [4], for the correspondence principle to be applicable, the following rules must apply: 1. A point 𝒙 on the body must be part of the body for all times 𝑡 ≥ 0. 15 2. The motion of the viscoelastic body must be quasi-static. 3. For a fixed point on the boundary, the boundary condition cannot be changed during the deformation. 2.2 Nonlinear Viscoelasticity 2.2.1 Overview As viscoelastic materials deform over time, at some point in their deformation history, the linear viscoelastic limit will be exceeded. When this limit is exceeded, the material response transitions from linear viscoelastic into nonlinear viscoelastic. The stress – strain response of the material in this deformation region becomes more complex. The emergence of a nonlinear viscoelastic response is marked by the dependence of the stress relaxation modulus (or creep compliance) function of the viscoelastic response on the applied strain (or stress) level. 𝜎(𝑡) 𝐺(𝑡, 𝜀) = (2.13) 𝜀0 𝜀(𝑡) 𝐽(𝑡, 𝜎) = (2.14) 𝜎0 Due to the modifications to the stress relaxation and creep compliance functions, the Boltzmann superposition principle can no longer be assumed to be applicable. To describe stepwise loading response in nonlinear viscoelasticity, modified superposition principles or methods will have to be employed. Furthermore, nonlinear viscoelastic models have been developed to predict the response of viscoelastic materials under deformations beyond the linear viscoelastic limit. Subsections 2.2.2 and 2.2.3 of this chapter present a summary of these models. 2.2.2 Schapery Nonlinear Viscoelastic Model The Schapery nonlinear viscoelastic model or the single integral model is one of the most widely used models for predicting the stress – strain response of viscoelastic materials beyond the linear viscoelastic limit. The constitutive theory was developed by Schapery using the principles of thermodynamics of irreversible processes [7-9]. The stiffness-based uniaxial formulation of the single integral model is expressed as: 𝑡 𝑑ℎ2 𝜀(𝜏) 𝜎(𝑡) = ℎ∞ ∙ 𝐺∞ ∙ 𝜀(𝑡) + ℎ1 ∙ ∫ ∆𝐺 [𝜌(𝑡) − 𝜌(𝜏)] 𝑑𝜏 (2.15) 𝑑𝜏 0 16 where 𝐺∞ and ∆𝐺(𝜌)are the equilibrium and transient linear viscoelastic relaxation modulus. ℎ∞ , ℎ1 , 𝑎𝑛𝑑 ℎ2 are strain-dependent nonlinearity parameters, and 𝜌 is the reduced-time defined as: 𝑡 𝑑𝑡 ′ 𝜌(𝑡) = ∫ (𝑎𝜀 > 0) (2.16) 0 𝑎𝜀 [𝜀(𝑡 ′ )] 𝜏 𝑑𝑡 ′ 𝜌(𝜏) = ∫ ′ (2.17) 0 𝑎𝜀 [𝜀(𝑡 )] According to Schapery [9], the nonlinearity parameters have specific thermodynamic significance. Variations in the first three nonlinearity parameters (ℎ∞ , ℎ1 , 𝑎𝑛𝑑 ℎ2 ) are due to third and higher-order strain effects in the Helmholtz free energy. Changes in 𝑎𝜺 result from strong strain influences in both entropy production and free energy. Irrespective of the thermodynamic significance of these parameters, their placement in the constitutive equation is indicative of how these parameters affect the physical response of the viscoelastic material. The ℎ∞ term represents the measure of nonlinearity in the equilibrium modulus. ℎ1 is the measure of the nonlinearity in the transient modulus and ℎ2 represents the measure of the nonlinearity in the strain rate effect. Lastly, the parameter 𝑎𝜺 is a time shift factor that can be strain and temperature dependent. Furthermore, when the input strain is within the linear viscoelastic strain limit, the values of the nonlinearity parameters ℎ∞ , ℎ1 , ℎ2 , and 𝑎𝜺 become equal to 1, and Eqn.2.15 reduces to the hereditary integral representation for linear viscoelastic materials (Eqn.2.5). The Schapery nonlinear viscoelastic model can be readily modified to include environmental effects such as temperature and can be tailored into numerical procedures. Furthermore, nonlinear viscoelastic models have been implemented in numeric simulations for predicting the response of polymeric materials and composites. Recurrent numerical algorithms [10,11] have been applied to the Schapery nonlinear viscoelastic model for numerical and finite element (FE) modeling of viscoelastic materials. Hence, the Schapery model has been successfully implemented over the past several decades to successfully model the nonlinear viscoelastic response of different materials [12-20]. Different parameter estimation methods as well as modifications to the Schapery model have been made to extend its applications. 2.2.3 Other Nonlinear Viscoelastic Models Over the past several decades, nonlinear viscoelastic models have been developed to predict the response of polymeric and composite materials. A multiple-integral formulation for the mechanics of nonlinear viscoelastic response was proposed by Green and Rivlin [21]. 17 Furthermore, Findley and Lai [22], developed a modified superposition principle for predicting the creep response in the nonlinear deformation range under time-dependent step loadings. Experimental results for polyvinyl chloride under tensile creep showed the superposition principle described the response of the polymer well. Pipkin and Rodgers [23] presented a finite strain integral model for nonlinear viscoelastic response to arbitrary stress or strain histories. Generalization to the three-dimensional form of the integral model, as well as the implications of isotropy was presented. Predictions of the single integral approximations agreed closely with experimental data from creep tests, for most cases. Furthermore, the concept of adaptive springs, originally introduced by Green and Tobolsky [24] was used by Drozdov et al [25] to model the chemical links between polymer molecules. The model was verified using creep data for polypropylene fibers and the predictions agreed well with experimental data. Drozdov continued to develop a nonlinear model based on a network of adaptive springs and was able to include thermal and aging effects [26-28]. Kolarik et al [29] proposed a Boltzmann-like superposition principle to analyze the multistep nonlinear creep response and demonstrated the method for three polypropylene materials. Knauss and Emri [30] used a single integral method to associate nonlinear parameters to free volume which allowed for the inclusion of stress-induced dilatation, moisture or other diffusion parameters in the theory. The free volume theory [30] was used by Bosi et al [31] to develop a nonlinear thermoviscoelastic model to characterize the response of orthotropic thin membranes up to yielding. The model accounted for the temperature and rate dependence under various mechanical loading conditions. Other nonlinear viscoelastic models such as the free volume-based model for a wide range of strain rates and temperatures [32], a fractional multiaxial nonlinear viscoelastic model [33], a nonlinear viscoelastic model based on the Ogden nonlinear elastic strain energy [34], and a differential nonlinear viscoelastic model for thermoset polymers under complex loadings [35] have also been successfully implemented in numerical modeling and FE analysis of viscoelastic materials. These constitutive models discussed above were successful in predicting the mechanical behavior of viscoelastic materials in the nonlinear deformation range and present an excellent basis for understanding nonlinear viscoelastic theory and the role constitutive modeling plays in understanding the mechanical response of viscoelastic materials under finite strains. 18 2.3 Environmental Effects and Long-Term Material Response The ability to predict material response over a long period or over its structural life span is a very important aspect of engineering design. Methods for predicting the long-term response of viscoelastic materials using short-term data provide a huge advantage in analyzing the response of viscoelastic materials. For polymeric materials, the time-temperature superposition principle (TTSP) presents a pathway for long-term material prediction. The TTSP has been extended in literature to introduce other environmental effects such as moisture in the form of time- temperature-moisture superposition principles. These superposition principles will be covered extensively in this section. 2.3.1 Time-Temperature Superposition Principle (TTSP) The theoretical foundation for the TTSP is the kinetic theory of polymers which has been extensively studied and extended for many applications [36-39]. According to the principle, there is a time-temperature equivalence in the response of viscoelastic materials. This simply implies that the viscoelastic response at a specific temperature is related to that at another by a change in the time scale only [3,4,39-41]. Hence, the time and temperature variation of the relaxation moduli (or compliances) of a viscoelastic material is often said to be equivalent. Using this principle, the response of a viscoelastic material for a long-time span can be determined using short-term tests for a range of temperatures. The shift in the time scale or the horizontal time shift (𝑎 𝑇 ) is expressed as the ratio of the relaxation time at one temperature (𝑇) to that at a reference temperature (𝑇0): 𝜏(𝑇) 𝑎𝑇 = (2.18) 𝜏(𝑇0 ) Typically, polymeric materials have many relaxation times. A polymeric material is considered thermorheologically simple if the same shift factor applies to all relaxation times. This means that all relaxation times of thermorheologically simple materials are affected by temperature in the same way. Hence, the TTSP is only applicable to thermorheologically simple materials. The shift factors are typically determined experimentally by building master curves. The horizontal time shift factor is usually fitted using the William-Landel-Ferry equation [40] or the Arrhenius equation [3]. The choice of the fitting equation depends on the considered temperature range and the glass transition temperature of the material (𝑇𝑔 ). 19 The TTSP was applied to a number of polymers by Williams, Landel and Ferry [40], and they found an empirical expression for the time shift factor: −𝐶1 ∙ (𝑇 − 𝑇0 ) 𝐿𝑜𝑔 (𝑎𝑇 ) = (2.19) 𝐶2 + (𝑇 − 𝑇0 ) where 𝐶1 and 𝐶2 are material constants. The WLF equation is applicable if the TTSP is implemented above 𝑇𝑔 and if considered temperature range is less than 𝑇𝑔 + 100°𝐶 or less than 𝑇0 + 50°𝐶 [40]. On the other hand, when the TTSP is implemented below 𝑇𝑔 , or at temperatures above 𝑇𝑔 + 100°𝐶 or 𝑇0 + 50°𝐶, it is appropriate to fit the time shift factor according to the Arhenius equation. The Arrhenius equation is expressed as: ∆𝐻 1 1 𝐿𝑜𝑔 (𝑎𝑇 ) = [ − ] (2.20) 2.303 ∙ 𝑅 𝑇 𝑇0 where ∆𝐻 is the activation energy and R is the universal gas constant. Furthermore, when implementing the TTSP it is common to introduce small vertical shifts (𝑏𝑇 ) along the relaxation modulus (or creep compliance) scale. It has been suggested that the inclusion of vertical shifts accounts for the change in the degree of crystallinity of the considered polymeric material (in the case of semi-crystalline polymers) [9] and improves the overall accuracy of the long-time prediction of the viscoelastic response [41]. To implement the TTSP, master curves are built by shifting individual relaxation (or creep) curves at different temperatures onto a curve at a reference temperature on a log-log plot. With respect to the vertical and horizontal shift factors, the stress relaxation modulus at a given temperature level can be expressed as: 𝑡 𝑏𝑇 ∙ 𝐺(𝑡, 𝑇) = 𝐺 [ , 𝑇0 ] (2.21) 𝑎𝑇 𝑡 where 𝐺 [𝑎 , 𝑇0 ] is the stress relaxation modulus function for the master curve built for the 𝑇 reference temperature 𝑇0. Graphically, the implementation of the TTSP is shown in Fig.2.8 below. 20 Figure 2.8 Implementation of the time-temperature superposition principle. The TTSP has been successfully applied in literature, to predict the long-time and temperature-dependent response of viscoelastic materials in the linear [42-46] and nonlinear deformation range [47-50]. 2.3.2 Time-Temperature-Moisture Superposition Method (TTMSM) The time-temperature-moisture superposition method (TTMSM) is an extension of the TTSP. This method was developed from the need to predict the long-time response of viscoelastic material in environmental or working conditions that cause a change in the moisture content of the material. The TTMSM was developed based on the assumption that the change in moisture content has the same effect on a viscoelastic material as the change in temperature. According to this principle, there is a time-moisture equivalence in the response of viscoelastic materials. This implies that the viscoelastic response of a material at one moisture content level can be related to that at another by changes in the time scale alone. This introduces a new shift factor known as the time-moisture shift factor. This method studies the plasticization effect of moisture content or/and temperature on the response of viscoelastic materials. The TTMSM has been implemented in literature for predicting the long-time viscoelastic response of polymeric materials from accelerated test data [51-58]. There are two main ways in which the TTMSM can be implemented. The first method involves carrying out stress relaxation, creep or frequency sweep tests at a constant temperature and different humidity levels. For this method, the test curves are superimposed onto a reference curve at a reference humidity level to form a time-moisture master curve. The second method 21 involves the effects of both moisture content and temperature. For this method, tests are carried out at different constant temperatures and different humidity levels. Time-temperature master curves a built for different constant humidity levels by shifting the curves onto that at a set reference temperature. Then, a super master curve is built by superimposing the time-temperature master curves onto that at a reference humidity level. This second method has been implemented in the work by Suarez-Martinez et al to study the dynamical and rheological behavior of polyelectrolyte complexes (PECs) [58]. Figure 2.9 Application of the time−water superposition principle for PECs. (a) Time−temperature master curves for relative humidity (RH) values of 50, 70, 80, 85, 90, and 95%. (b) Time−water super master curve made from time−temperature master curves in (a) with RHref = 80% and Tref = 40°C [58]. 2.4 Viscoplastic Modeling of Polymeric Materials 2.4.1 Overview The observed response of polymeric materials is generally time-dependent. Their stress response depends on the loading rate and the considered time scale. When the stress response beyond the limit where deformations are recoverable (yield point) is rate dependent as in Fig.2.10, the material response is said to be viscoplastic in that region. This rate-dependent response of polymeric materials is described by viscoplasticity (rate-dependent plasticity) models. This section covers a literature review on viscoelastic constitutive models for predicting the response of polymeric materials. 22 Figure 2.10 Stress-strain response of polymeric materials: (a) elastic-viscoplastic response and (b) viscoelastic-viscoplastic response. 2.4.2 Viscoplastic Constitutive Models In applied mechanics, constitutive models are used to predict how materials respond to large deformations. There are two main categories of these models outlined by Perić and Owen [59]. They are micromechanical and phenomenological models. Micromechanical models attempt to simulate the behavior of individual molecules in a material to understand its overall response, while phenomenological models are developed based on macroscopic observations of the material. Several micromechanical models have been proposed to model the large deformation response of polymeric materials. In the work by Parks and Ahzi [60], a viscoplastic model was developed that satisfied local kinematic constraints as well as global compatibility and applied to large deformation and texturing of orthorhombic polycrystalline materials. Lee et al [61], developed a micromechanically-based large deformation model for semi-crystalline polymers. The model was applied to predict the stress-strain response and texture evolution in high-density polyethylene under different straining modes. Furthermore, Drozdov and co-workers [62-65] have carried out extensive work on micromechanical modeling of polymeric materials under large deformations. Other physical micromechanical models have been developed and successfully implemented for the prediction of large deformation responses of polymeric materials [66,67]. However, a limitation of micromechanical models is that their predictive capabilities are restricted to specific conditions. Also, micromechanical models are not scalable to macroscale modeling of materials that exhibit different deformation characteristics at different scales. A typical example is semicrystalline polymers which are highly anisotropic on the microscale [68] but exhibit isotropic 23 deformation on the macroscale. This scaling discontinuity adds to the complexities of implementing micromechanical models on a practical scale. Macroscopic constitutive theories, also known as phenomenological models have been expressed in simple forms using rheological models. An elasto-viscoelastic-viscoplastic model termed “Burger’s model” was developed by Findlay et al [69] to model the stress-strain response of polymeric materials. The Burger’s model assumes that irrecoverable deformations occur at all stages of the material deformation and this assumption eliminates the yield surface which is typical of polymeric materials. A nonlinear elasto-viscoelastic-viscoplastic model was developed from thermodynamic principles by Schapery [70,71], to model the response of fiber-reinforced and unreinforced plastics, as well as semicrystalline polymers. Lai and Bakker [72] proposed an elasto- viscoelastic-viscoplastic integral model for polymeric materials. The developed model was validated using creep and recovery tests for high-density polypropylene. Furthermore, more coupled elasto-viscoelastic-viscoplastic models based on rheological frameworks have been developed and implemented for predicting the stress-strain response of polymeric materials in their entire deformation range [73-78]. Figure 2.11 One-dimensional elasto-viscoelastic-viscoplastic rheological models used for polymers by (a) Findley et al. [69], (b) Schapery [70], (c)Brusselle-Dupend et al. [73, 74] and (d) Kletschkowski et al. [75].[79]. Macroscopic viscoplastic constitutive models can be further classified into two groups. The first group is concerned with the consistency models and the other comprises the overstress models. The group of consistency viscoplasticity models was introduced by Wang et al. [80] and 24 was studied extensively by many authors [81-84]. The consistency viscoplastic models incorporate the time-dependent response in the large strain deformation region in a rate-dependent yield surface. These models are analogous to rate-independent plasticity in that they follow the Kuhn- Tucker optimality conditions. The second group of viscoplastic models is the overstress model whose main postulate is that the current stress state can be outside the yield surface and that the yield function can be greater than zero. These models were developed based on the ideas presented by Perzyna [85] and they disregard the Kuhn-Tucker optimality conditions. The Perzyna model [85] and the Duvaut-Lions model [86], are the most popular formulations in the group of overstress models. Both models have been widely implemented in small strain problems [87-95], as well as finite strain problems [96-99]. Furthermore, phenomenological constitutive models have been developed to predict the flow stress of polymeric materials. This group of phenomenological models are known as flow stress models. As long as the flow stress behavior of the polymeric material can be represented, these models can be easily calibrated by fitting experimental data. G’sell and Jonas [100] developed a constitutive relationship to predict the flow stress behavior of polymeric materials at different constant strain rates. This model has been implemented in its original and modified forms in different works to predict the stress-strain response of polymeric materials at different strain rates and temperatures [101-104]. The Johnson-Cook model was proposed by Johnson and Cook [105], to predict flow stress at different strain rates and temperatures. This model has been popularly applied to modeling the stress-strain response of metals. However, the Johnson-Cook model and its variations have been applied to describe the tensile and compressive response of PEEK [106,107]. A flow stress model based on the work by G’sell and Jonas [100], was proposed by Nasraoui et al [108] for strain rate and the temperature-dependent response of polymeric materials. The model was validated against uniaxial compression tests for PMMA (polymethylmethacrylate) under quasi-static and dynamic loadings. Furthermore, a uniform constitutive phenomenological model known as the DSGZ model was proposed by Duan et al [109] for glassy and semicrystalline polymers. The model was proposed based on several constitutive models [100,105,110] and has been used to predict the flow stress of PMMA [109] and PC [109,111]. Recently, Zhu et al [112], developed a new phenomenological model for predicting the flow stress of thermoplastics. In the model included a transition function was introduced to enable a smooth transition of the flow stress behavior under both small and large 25 deformations. The model was validated against uniaxial tensile and compression tests for PEEK and PC. The constitutive model adequately predicted the trends in the response of the tested materials. The constitutive models discussed in this section present different ways in which the rate- dependent response of materials in the large deformation region can be described and predicted. These works also show how viscoplastic models can be coupled with viscoelastic and/or elastic models to capture the stress-strain response of polymeric materials in their entire range of deformation. They also show how flow stress models can be employed to predict the viscoplastic response as well as the entire stress-strain response of polymeric materials. These works present a foundation for developing new constitutive theories that can be implemented for more accurate numerical modeling of polymeric materials. 2.5 Thermomechanical Modeling of Polymeric Battery Separators 2.5.1 Thermomechanical Behavior of Polymeric Battery Separators In order to develop an accurate material model for polymeric separators, an understanding of their behavior under various mechanical and thermal loading conditions is necessary. Considering their importance in battery safety, the thermomechanical response of polymeric battery separators has been studied extensively in literature. Dry-processed polymeric battery separators have been found to display strong material anisotropy [113-119]. This class of separators show higher tensile strength and elastic modulus in the MD than the TD due to the extra strain hardening given in the MD during their manufacturing process. The mechanical properties of polymeric separators are also rate-dependent [113-115]. Uniaxial tensile tests were carried out on Celgard®2400 polypropylene separators at different loading rates and the results were presented in the work by Yan et al [113]. From the results, there was a clear relationship between the loading rate and the measured modulus of elasticity. The measured elastic modulus was found to progressively increase with the increase in the loading rate in all three material orientations. The material orthotropy was also depicted in the results as shown in Fig 2.12. 26 (a) (b) (c) Figure 2.12 Stress-strain response of Celgard®2400 separator at different loading rates in the (a) MD, (b) TD, and (c) 45 off-axis direction [100]. The mechanical properties of the separators are temperature-dependent as commonly observed in polymeric materials. It has been shown that the material becomes significantly softer due to a reduction in its mechanical properties, with an increase in temperature [1114,116]. Avdeev et al [114] showed the rate and the temperature-dependent response of polymeric battery separators through the stress-strain behavior of a Celgard®C480 Tri-layered separator at different strain rates and temperatures. With an increase in temperature, the material stiffness was found to decrease in all orientations. Furthermore, investigating the temperature-dependent response of polymeric separators is of utmost importance from the point of safety in thermal runaway scenarios. Dry- processed polymeric separators also show a unique thermal expansion/shrinkage behavior [120,121]. Their thermal expansion/shrinkage response is unique in the sense that with increasing temperature, polymeric separators first expand and then shrink before failure. This response is not typical of polymeric materials and can be attributed to the unique microstructure of the material. 27 The amount of material shrinkage at elevated temperatures can be as high as 10%. In a constrained condition, such as in battery cells, the shrinkage will induce tensile stresses in separators and may play a role in causing an internal short circuit. Yan et al [120] measured the thermal expansion/shrinkage behavior of three commonly used polymeric separators. The coefficients of thermal expansion (CTE) as a function of temperature were determined for the separators. It was also determined that the fracture/necking temperatures of the TD samples for all three separators investigated were 15-35°C lower than the MD samples. This also shows evidence of the material’s anisotropic nature in its thermal expansion/shrinkage response. 2.5.2 Plasticization Effects of Electrolyte Solutions To understand the mechanical response of the separator in its working environment (battery cells), the electrolyte effect needs to be taken into account. Electrolyte solutions have been found to remarkably weaken the mechanical performance of LIB separators [114,116,119,122,123]. It has been observed that when being immersed in dimethyl carbonate (DMC) which is a common solvent for electrolyte solutions, the tensile modulus of Celgard ®2400 decreased to 48.5% and 87.7% in the MD and TD respectively in terms of ratio to the measured modulus in air [119]. It has also been observed that electrolyte solutions markedly weakened the thermal stability at elevated temperatures and mechanical performance, especially on the crack resistance of the polymeric separators [122]. Furthermore, Gor et al [123], investigated the effect of certain electrolyte solvents on the thickness and elastic modulus of Celgard®3501 polypropylene separator. Their results showed that electrolyte solvents such as dimethyl carbonate, diethyl carbonate, and ethyl acetate caused reasonable softening of the separator. The softening of the separator material was attributed to swelling. Hence, thermomechanical models developed based on the response of the polymeric separator in air (or dry condition) will lead to an overprediction in its mechanical response under various mechanical and thermal loading conditions. 2.5.3 Existing Models for Polymeric Battery Separators In literature, the mechanical response of separators has been modeled in numerical simulations with material models. Wu et al. [124] presented the development of a multiphysics model for stress analysis in polypropylene separators in lithium-ion battery cells that considered the mechanical loading, the Li-induced intercalation and the thermal expansion mismatch between the components of the battery. In addition to the considered stress-inducing factors, the time and 28 temperature-dependent behavior of the PP separator was also characterized experimentally by implementing the time-temperature superposition principle (TTSP). However, this model did not consider the anisotropic behavior of the separator. Yan et al [125] developed an orthotropic viscoelastic material model for polymeric separators implemented in LS-DYNA® as a user-defined model. The model was able to accurately predict the material anisotropy and rate dependence. Furthermore, X. Zhang et al. [126] implemented an anisotropic crushable foam model (MAT 126) from the LS-DYNA® material library to model the force vs displacement behavior of a polyethylene separator in the TD, MD and through-thickness direction in tension and compression. In [127], a large deformation elastic–plastic constitutive model of a PP separator was developed with the Rich–Hill large deformation elastoplastic constitutive theory. The constitutive model accurately captured the anisotropy behaviors and the elastic–plastic process considering the large deformation of the separator. Also, an implicit nonlinear dynamic method was implemented in modeling dry-processed PP separators as elastic-viscoplastic [128]. Xie et al [129] proposed a representative volumetric element modeling method in finite element simulation to predict the electrolyte-immersed, rate-dependent tensile properties of polypropylene (PP) separators. These works have offered some insights into the behavior of polymeric separators. However, these models are not adequate for simulations of the response of separators in thermal ramp conditions. Without the combined consideration of material anisotropy, temperature dependency, and the electrolyte effect, these models are inadequate in simulations with thermal ramp scenarios. Hence, for more accurate predictions of the response of polymeric battery separators, models that account for the constitutive behaviors of separators in their working conditions need to be developed. 29 REFERENCES 1. R. S. Lakes, Viscoelastic solids. CRC press, 2017. 2. L. Boltzmann, Pogg. Ann. Phys. Chem., 7, 624, 1876. 3. I. M. Ward and J. Sweeney. An Introduction to the Mechanical Properties of Solid Polymers (Wiley, New York, NY), 2004. 4. A. S. Wineman, K. R. Rajagopal. Mechanical Response of Polymers: An Introduction (Cambridge university press, Cambridge, UK), 2000. 5. R. M. Christensen. Theory of Viscoelasticity: An Introduction (Academic Press, New York, NY), 1982. 6. A. S. Argon. The physics of deformation and fracture of polymers. New York.: Cambridge, 2013. 7. R. A. Schapery. PTOC. 5th U.S. National Congress of Applied Mechanics, ASME, 511, 1966. 8. R. A. Schapery. Purdue University, Rept. No. 69-2, 1969. 9. R. A. Schapery. On the characterization of nonlinear viscoelastic materials. Polymer Engineering & Science, 9(4), pp.295-310, 1969. 10. R. M. Haj‐Ali and A. H. Muliana. Numerical finite element formulation of the Schapery non‐ linear viscoelastic material model. International Journal for Numerical Methods in Engineering, 59(1), pp.25-45, 2004. 11. J. Lai and A. Bakker. 3-D Schapery representation for non-linear viscoelasticity and finite element implementation. Computational Mechanics, 18(3), pp.182-191, 1996. 12. S. W. Park and R. A. Schapery. Methods of interconversion between linear viscoelastic material functions. Part 1 – a numerical method based on Prony series, International Journal of Solids and Structures, vol. 36, pp. 1653-1675, 1999. 13. E. Ségard, S. Benmedakhene, A. Laksimi, and D. Lai. Influence of the fibre-matrix interface on the behaviour of polypropylene reinforced by short glass fibres above the glass transition temperature, Composite Science and Technology, vol. 62, pp. 2029-2036, 2002. 14. J. Lai and A. Bakker. An integral constitutive equation for nonlinear plastoviscoelastic behaviour of high-density polyethylene, Polymer Engineering and Science, vol. 35, no. 17, pp. 1339-1347, 1995. 15. S. P. Zaoutsos, G. C. Papanicolaou, and A. H. Cardon. On the non-linear viscoelastic behaviour of polymer-matrix composites, Composites Science and Technology, vol. 58, pp. 883-889, 1998. 16. E. Marklund, J. Varna, and L. Wallström. Nonlinear viscoelasticity and viscoplasticity of flx/polypropylene composites, Journal of Engineering Materials and Technology, vol. 128, pp. 527-536, 2006. 17. P. Dasappa, P. Lee-Sullivan, X. Xiao, and H. P. Foss. Tensile creep of a long-fiber glass mat thermoplastic (GMT) composite Part II: Viscoelastic-Viscoplastic constitutive modeling, Polymer Composites, 2008. 30 18. A. Pasricha, M.E. Tuttle, A.F. Emery. Time-dependent response of IM7/5260 composites subjected to cyclic thermo-mechanical loading, Composites Science and Technology, vol. 55, pp. 49-56, 1995. 19. D. Peretz and Y. Weitsman. Nonlinear viscoelastic Characterization of FM-73 Adhesive, Journal of Rheology, vol. 26, no. 3, pp. 245-261, 1982. 20. D. J. Pooler and L.V. Smith. Nonlinear viscoelastic response of a wood-plastic composite including temperature effects, Journal of Thermoplastic Composite Materials, vol. 17, pp. 427- 445, 2004. 21. A. E. Green and R. S. Rivlin. The mechanics of non-linear materials with memory. Archive for rational mechanics and analysis, 1(1), pp.1-21, 1957. 22. W. N. Findley and J. S. Y. Lai. A modified superposition principle applied to creep of nonlinear viscoelastic material under abrupt changes in state of combined stress. Transactions of the Society of Rheology 11, no. 3, 361-380, 1967. 23. A. C. Pipkin and T. G. Rogers. A non-linear integral representation for viscoelastic behavior. Journal of the Mechanics and Physics of Solids, 16(1), pp.59-72, 1968. 24. M. S. Green and A. V. Tobolsky. A new approach to the theory of relaxing polymeric media. The Journal of chemical physics, 14(2), pp.80-92, 1946. 25. A. D. Drozdov and A. L. Kalamkarov. A constitutive model for nonlinear viscoelastic behavior of polymers, Polymer Engineering and Science, vol. 36, no. 14, pp. 1907-1919, 1997. 26. A. D. Drozdov. A model of adaptive links in finite viscoelasticity, Mechanics Research Communications, vol. 24, no. 2, pp 161-166, 1997. 27. A. D. Drozdov. A constitutive model in finite thermoviscoelasticity based on the concept of transient networks, Acta Mechanica, vol. 133, pp. 13-37, 1999. 28. A. D. Drozdov. A model for the non-linear viscoelastic response and physical aging in glassy polymers, Computations and Theoretical Polymer Science, vol. 9, pp. 73-87, 1999. 29. J. Kolařík and A. Pegoretti. Proposal of the Boltzmann-like superposition principle for nonlinear tensile creep of thermoplastics. Polymer Testing, 27(5), pp.596-606, 2008. 30. W. Knauss and I. Emri. Non-linear viscoelasticity based on free volume consideration. In Computational Methods in Nonlinear Structural and Solid Mechanics, pp. 123-128, Pergamon, 1981. 31. F. Bosi and S. Pellegrino. Nonlinear thermomechanical response and constitutive modeling of viscoelastic polyethylene membranes. Mechanics of Materials, 117, pp.9-21, 2018. 32. Gong, Y. Chen, T. Li, Z. Liu, Z. Zhuang, B. Guo, H. Wang, and L. Dai. Free volume based nonlinear viscoelastic model for polyurea over a wide range of strain rates and temperatures. Mechanics of Materials, 152, p.103650, 2021. 33. A. Muliana. A fractional model of nonlinear multiaxial viscoelastic behaviors. Mechanics of Time-Dependent Materials, pp.1-21, 2022. 34. Moghanian, A. Pazhouheshgar and A. Ghorbanoghli. Nonlinear viscoelastic modeling of synthesized silicate-based bioactive glass/polysulfone composite: theory and medical applications. Silicon, 14(2), pp.731-740, 2022. 31 35. F. Ellyin and Z. Xia. Nonlinear Viscoelastic Constitutive Model for Thermoset Polymers. ASME. Journal of Engineering Materials and Technology,128(4): 579–585, 2006. 36. L. E. Nielsen. The Mechanical Properties of Polymers, Reinhold, NY, 1965. 37. L. E. Nielsen and R. F. Landel. The Mechanical Properties of Polymers and Composites, 2nd ed. M. Decker, NY, 1994. 38. A. V. Tobolsky. Properties and Structure of Polymers, JW, NY, 1962. 39. J. D. Ferry. Viscoelastic Properties of Polymers, 3rd ed., JW, NY, 1980. 40. M. L. Williams, R. F. Landel, and J. D. Ferry. The Temperature Dependence of Relaxation Mechanisms in Amorphous Polymers and Other Glass Forming Liquids, Journal of the American Chemical Society, 77, 3701-3707, 1955. 41. H. F. Brinson, L. C. Brinson. Polymer engineering science and viscoelasticity. An introduction, 2008. 42. R. C. Ihuaenyi, S. Yan, J. Deng, C. Bae, I. Sakib, and X. Xiao. Orthotropic Thermo- Viscoelastic Model for Polymeric Battery Separators with Electrolyte Effect. Journal of The Electrochemical Society, 168(9), p.090536, 2021. 43. R. M. Guedes. A viscoelastic model for a biomedical ultra-high molecular weight polyethylene using the time–temperature superposition principle. Polymer testing, 30(3), pp.294-302, 2011. 44. J. Zhao, W. G. Knauss and G. Ravichandran, Applicability of the time–temperature superposition principle in modeling dynamic response of a polyurea. Mechanics of Time- Dependent Materials, 11(3), 289, 2007. 45. F. J. Wortmann and K.V. Schulz. Stress relaxation and time/temperature superposition of polypropylene fibres. Polymer, 36(2), pp.315-321, 1995. 46. J. Capodagli and R. S. Lakes. Isothermal viscoelastic properties of PMMA and LDPE over 11 decades of frequency and time: a test of time-temperature superposition. Rheologica Acta, 47:777–786, 2008. 47. E. Federico, J.L. Bouvard, C. Combeaud and N. Billon. Large strain/time dependent mechanical behaviour of PMMAs of different chain architectures. Application of time- temperature superposition principle. Polymer, 139, pp.177-187, 2018. 48. J. Diani, P. Gilormini and J.S. Arrieta. Direct experimental evidence of time-temperature superposition at finite strain for an amorphous polymer network. Polymer, 58, pp.107-112, 2015. 49. R. Nevière. An extension of the time–temperature superposition principle to non-linear viscoelastic solids. International journal of solids and structures, 43(17), pp.5295-5306, 2006. 50. P. A. O'connell and G.B. McKenna. Large deformation response of polycarbonate: Time‐ temperature, time‐aging time, and time‐strain superposition. Polymer Engineering & Science, 37(9), pp.1485-1495, 1997. 51. R. D, Maksimov, V. P. Mochalov, and Y. S. Urzhumtsev. Time—moisture superposition. Polymer Mechanics, 8(5), pp.685-689, 1972. 32 52. S. Onogi, K. Sasaguri, T. Adachi, and S. Ogihara. Time–humidity superposition in some crystalline polymers. Journal of Polymer Science, 58(166), pp.1-17, 1962. 53. F. Huber, H. Etschmaier, H. Walter, G, Urstöger, and P. Hadley. A time–temperature–moisture concentration superposition principle that describes the relaxation behavior of epoxide molding compounds for microelectronics packaging. International Journal of Polymer Analysis and Characterization, 25(6), pp.467-478, 2020. 54. A. Ishisaka, and M. Kawagoe. Examination of the time–water content superposition on the dynamic viscoelasticity of moistened polyamide 6 and epoxy. Journal of Applied Polymer Science, 93(2), pp.560-567, 2004. 55. V. Fabre, G. Quandalle, N. Billon, and S. Cantournet. Time-Temperature-Water content equivalence on dynamic mechanical response of polyamide 6, 6. Polymer, 137, pp.22-29, 2018. 56. J. Yao and G. Ziegmann. Equivalence of moisture and temperature in accelerated test method and its application in prediction of long-term properties of glass-fiber reinforced epoxy pipe specimen. Polymer testing, 25(2), pp.149-157, 2006. 57. S. M. Zhou, K. Tashiro, and T. Ii. Confirmation of universality of time–humidity superposition principle for various water‐absorbable polymers through dynamic viscoelastic measurements under controlled conditions of relative humidity and temperature. Journal of Polymer Science Part B: Polymer Physics, 39(14), pp.1638-1650, 2001. 58. P. C. Suarez-Martinez, P. Batys, M. Sammalkorpi, and J. L. Lutkenhaus. Time–temperature and time–water superposition principles applied to poly (allylamine)/poly (acrylic acid) complexes. Macromolecules, 52(8), pp.3066-3074, 2019. 59. Peric, and D. R. J. Owen. Finite-element applications to the nonlinear mechanics of solids. Reports on Progress in Physics, 61, 1495–1574, 1998. 60. M. Parks, and S. Ahzi. Polycrystalline plastic deformation and texture evolution for crystals lacking five independent slip systems. Journal of the Mechanics and Physics of Solids, 38(5), pp.701-724, 1990. 61. J. Lee, D. M. Parks, and S. Ahzi. Micromechanical modeling of large plastic deformation and texture evolution in semi-crystalline polymers. Journal of the Mechanics and Physics of Solids, 41(10), pp.1651-1687, 1993. 62. A. D. Drozdov and J. deC. Christiansen. Model for the viscoelastic and viscoplastic responses of semicrystalline polymers. Journal of Applied Polymer Science, 88,1438–1450, 2003. 63. A. D. Drozdov and J. deC. Christiansen. Constitutive equations for the viscoplastic response of isotactic polypropylene in cyclic tests: The effect of strain rate. Polymer Engineering and Science, 44, 548–556, 2004. 64. A. D. Drozdov and Q. Yuan. Effect of annealing on the viscoelastic and viscoplastic responses of low-density polyethylene. Journal of Polymer Science, Part B: Polymer Physics 41, 1638– 1655, 2003. 65. A. D. Drozdov. Modeling the viscoelastoplastic behavior of amorphous glassy polymers. Polymer Engineering and Science, 41, 1762–1770, 2001. 33 66. J. A. W.van Dommelin, D. M. Parks, M. C. Boyce, W. A. M. Brekelmans, and F. P. T. Baaijens. Micromechanical modeling of elasto-viscoplastic behavior of semi-crystalline polymers. Journal of the Mechanics and Physics of Solids, 51, 519-541, 2003. 67. A. Dahoun. Comportement plastique et textures de deformation des polymeres semi-cristallins en traction uniaxiale et en cisaillement simple. Ph.D. Dissertation, Institut National Polytechnique de Lorraine, France, 1992. 68. A. D. Drozdov, S. Agarwal, and R. K. Gupta. The effect of temperature on the viscoelastic behavior of linear low-density polyethylene. Archive Applied Mechanics 73, 591–614, 2004. 69. W. N. Findley, J. S. Lai, and K. Onaran. Creep and relaxation of nonlinear viscoelastic materials. New York: Dover Publications, Inc, 1976. 70. R. A. Schapery. Nonlinear viscoelastic and viscoplastic constitutive equations based on thermodynamics. Mechanics of Time-Dependent Materials, 1 pp.209–240, 1997. 71. R. A. Schapery. Nonlinear viscoelastic and viscoplastic constitutive equations with growing damage. International Journal of Fracture. 97, 33–66, 1999. 72. J. Lai and A. Bakker. An integral constitutive equation for nonlinear plasto-viscoelastic behavior of high-density polyethylene. Polymer Engineering and Science, 35,1339–1347, 1995. 73. N. Brusselle-Dupend, D. Lai, X. Feaugas, M. Guigon, and M. Clavel. Mechanical behavior of a semicrystalline polymer before necking. Part I: Characterization of uniaxial behavior. Polymer Engineering and Science, 41, 66–76, 2001. 74. N. Brusselle-Dupend, D. Lai, X. Feaugas, M. Guigon, and M. Clavel. Mechanical behavior of a semicrystalline polymer before necking. Part II: Modeling of uniaxial behavior. Polymer Engineering and Science, 43, 501–518, 2003. 75. T. Kletschkowski, U. Schomburg, and A. Bertram. An endochronic viscoplastic approach for materials with different behavior in tension and compression. Mechanics of Time-Dependent Materials, 8, 119–135, 2004. 76. C. Zhang and I. D. Moore. Nonlinear mechanical response of high-density polyethylene. Part I: Experimental investigation and model evaluation. Polymer Engineering and Science, 37, 404–413, 1997. 77. C. Zhang and I. D. Moore. Nonlinear mechanical response of high-density polyethylene. Part II: Uniaxial constitutive modeling. Polymer Engineering and Science, 37, 414–420, 1997. 78. D. Peric and W. Dettmer. A computational model for generalized inelastic materials at finite strains combining elastic, viscoelastic and plastic material behavior. Engineering Computations. 20, 768–787, 2003. 79. W. Holmes, J. G. Loughran, and H. Suehrcke. Constitutive model for large strain deformation of semicrystalline polymers. Mechanics of Time-Dependent Materials, 10(4), pp.281-313, 2006. 80. W. M. Wang, L. J. Sluys, and R. R. Borst. Viscoplasticity for instabilities due to strain softening and strain-rate softening, International Journal for Numerical Methods in Engineering 40 (20) 3839–3864, 2-6, (1997). 34 81. M. Ristinmaa and N. S. Ottosen. Consequences of dynamic yield surface in viscoplasticity, International Journal of Solids and Structures 37 (33) 4601–4622, 2000. 82. A. Carosio, K. Willam, and G. Etse. On the consistency of viscoplastic formulations, International Journal of Solids and Structures 37 (48) 7349–7369, 2000. 83. O. M. Heeres, A. S. J. Suiker, and R. De Borst. A comparison between the Perzyna viscoplastic model and the consistency viscoplastic model, European Journal of Mechanics, A/Solids 21 (1) 1–12, 2002. 84. R. Zaera and J. Fern´andez-S´aez. An implicit consistent algorithm for the integration of thermoviscoplastic constitutive equations in adiabatic conditions and finite deformations, International Journal of Solids and Structures 43 (6) 1594–1612, 2006. 85. P. Perzyna. Fundamental Problems in Viscoplasticity, Advances in Applied Mechanics 9 (C) 243–377, 1966. 86. Duvaut and J. Lions. Les Inequations en Mecanique et en Physique, Dunod, Paris, 1972. 87. O. C. Zienkiewicz and I. C. Cormeau. Visco-Plasticity-Plasticity and Creep in Elastic Solids - A unified numerical solution approach, International Journal for Numerical Methods in Engineering 8 (March) 821–845, 1974. 88. T. J. Hughes and R. L. Taylor. Unconditionally stable algorithms for quasi-static elasto/visco- plastic finite element analysis, Computers and Structures 8 (2) 169– 173, 1978. 89. I. Cormeau. Numerical stability in quasistatic elasto/viscoplasticity, International Journal for Numerical Methods in Engineering 9 (1) 109–127, 1975. 90. J. C. Simo, J. G. Kennedy, and S. Govindjee. Non‐smooth multisurface plasticity and viscoplasticity. Loading/unloading conditions and numerical algorithms. International Journal for Numerical Methods in Engineering, 26(10), pp.2161-218, 1988. 91. J. L. Chaboche. Constitutive equations for cyclic plasticity and cyclic viscoplasticity, International Journal of Plasticity, 247–302, 1989. 92. D. Peric. On a class of constitutive equations in viscoplasticity: Formulation and computational issues, International Journal for Numerical Methods in Engineering 36 (8), 1365–1393, 1993. 93. M. Ristinmaa and N. S. Ottosen. Viscoplasticity based on an additive split of the conjugated forces, European Journal of Mechanics, A/Solids 17 (2) 207–235, 1998. 94. K. Runesson, M. Ristinmaa, and L. Mahler. Comparison of viscoplasticity formats and algorithms, Mechanics of Cohesive-Frictional Materials 4 (1) 75–98, 1999. 95. A. Caggiano, E. Martinelli, D. Said Schicchi, and G. Etse. A modified Duvaut-Lions zero- thickness interface model for simulating the rate-dependent bond behavior of FRP-concrete joints. Composites Part B: Engineering, 149, pp.260-267, 2018. 96. A. Ibrahimbegovic and L. Chorfi. Viscoplasticity model at finite deformations with combined isotropic and kinematic hardening. Computers & Structures, 77(5), pp.509-525, 2000. 97. B. Nedjar. Frameworks for finite strain viscoelastic-plasticity based on multiplicative decompositions. Part I: Continuum formulations. Computer Methods in Applied Mechanics and Engineering, 191(15-16), pp.1541-1562, 2002. 35 98. V. Shutov and R. Kreißig. Finite strain viscoplasticity with nonlinear kinematic hardening: Phenomenological modeling and time integration. Computer Methods in Applied Mechanics and Engineering, 197(21-24), pp.2015-2029, 2008. 99. K. Kowalczyk-Gajewska, E. A. Pieczyska, K. Golasinski, M. Maj, S. Kuramoto, and T. Furutab. A finite strain elastic-viscoplastic model of Gum Metal. International Journal of Plasticity, 119, pp.85-101, 2019. 100. G'sell, and J. J. Jonas. Determination of the plastic behaviour of solid polymers at constant true strain rate. Journal of materials science, 14(3), pp.583-591, 1979. 101. Morin, G. Haugou, F. Lauro, B. Bennani and B. Bourel. Elasto-viscoplasticity behaviour of a structural adhesive under compression loadings at low, moderate and high strain rates. Journal of Dynamic Behavior of Materials, 1(2), pp.124-135, 2015. 102. J. I. Múgica, L. Aretxabaleta, I. Ulacia, and J. Aurrekoetxea. Rate-dependent phenomenological model for self-reinforced polymers. Composites Part A: Applied Science and Manufacturing, 84, pp.96-102, 2016. 103. M. Schoßig, C. Bierögel, W. Grellmann, and T. Mecklenburg. Mechanical behavior of glass-fiber reinforced thermoplastic materials under high strain rates. Polymer testing, 27(7), pp.893-900, 2008. 104. J. P. Torres, P. M. Frontini, and L. Aretxabaleta. Experimental characterization and computational simulations of the low‐velocity impact behaviour of polypropylene. Polymer international, 62(11), pp.1553-1559, 2013. 105. R. Johnson and W. H. Cook. A constitutive model and data for materials subjected to large strains, high strain rates, and high temperatures. Proceedings of the 7th International Symposium on Ballistics, pp.541-547, 1983. 106. Chen, H. Ou, B. Lu, and H. Long. A constitutive model of polyether-ether-ketone (PEEK). Journal of the mechanical behavior of biomedical materials, 53, pp.427-433, 2016. 107. Garcia-Gonzalez, A. Rusinek, T. Jankowiak, and A. Arias. Mechanical impact behavior of polyether–ether–ketone (PEEK). Composite Structures, 124, pp.88-99, 2015. 108. M. Nasraoui, P. Forquin, L. Siad, and A. Rusinek. Influence of strain rate, temperature and adiabatic heating on the mechanical behaviour of poly-methyl-methacrylate: experimental and modelling analyses. Materials & Design, 37, pp.500-509, 2012. 109. Y. Duan, A. Saigal, R. Greif, and M. A. Zimmerman. A uniform phenomenological constitutive model for glassy and semicrystalline polymers. Polymer Engineering & Science, 41(8), pp.1322-1328, 2001. 110. J. W. Brooks. Processing wrought nickel and titanium superalloys. Thermomechanical processing in theory, modelling and practice (TMP) 2, pp.52-77, 1997. 111. Wang, H. Zhou, Z. Huang, Y. Zhang, and X. Zhao. Constitutive modeling of polycarbonate over a wide range of strain rates and temperatures. Mechanics of Time-Dependent Materials, 21(1), pp.97-117, 2017. 112. Zhu, H. Ou, and A. Popov. A new phenomenological constitutive model for thermoplastics. Mechanics of Materials, 157, p.103817, 2021. 36 113. S. Yan, J. Deng, C. Bae, Y. He, A. Asta, and X. Xiao. In-plane orthotropic property characterization of a polymeric battery separator. Polymer Testing, 72, pp.46-54, 2018. 114. I. Avdeev, M. Martinsen, and A. Francis. Rate-and temperature-dependent material behavior of a multilayer polymer battery separator. Journal of materials engineering and performance, 23(1), pp.315-325, 2014. 115. Xu, L. Wang, J. Guan, and S. Yin. Coupled effect of strain rate and solvent on dynamic mechanical behaviors of separators in lithium-ion batteries. Materials & Design, 95, pp.319- 328, 2016. 116. Zhu, T. Wierzbicki, and W. Li. A review of safety-focused mechanical modeling of commercial lithium-ion batteries. Journal of Power Sources, 378, pp.153-168, 2018. 117. X. Zhang, E. Sahraei, and K. Wang. Li-ion battery separators, mechanical integrity and failure mechanisms leading to soft and hard internal shorts. Scientific reports, 6(1), pp.1-9, 2016. 118. J. Cannarella, X. Liu, C. Z. Leng, P. D. Sinko, G. Y. Gor, and C. B. Arnold. Mechanical properties of a battery separator under compression and tension. Journal of the Electrochemical Society, 161(11), p.F3117, 2014. 119. A. Sheidaei, X. Xiao, X. Huang, and J. Hitt. Mechanical behavior of a battery separator in electrolyte solutions. Journal of Power Sources, 196(20), pp.8728-8734, 2011. 120. S. Yan, J. Deng, C. Bae, and X. Xiao, Thermal expansion/shrinkage measurement of battery separators using a dynamic mechanical analyzer. Polymer Testing, 71, pp.65-71, 2018. 121. C.T. Love. Thermomechanical analysis and durability of commercial micro-porous polymer Li-ion battery separators. Journal of Power Sources, 196(5), pp.2905-2912, 2011. 122. J. Chen, H. Hu, S. Li, and Y. He. Evolution of mechanical properties of polypropylene separator in liquid electrolytes for lithium‐ion batteries. Journal of Applied Polymer Science, 135(27), p.46441, 2018. 123. Y. Gor, J. Cannarella, C. Z. Leng, A. Vishnyakov, and C. B. Arnold. Swelling and softening of lithium-ion battery separators in electrolyte solvents. Journal of Power Sources, 294, pp.167-172, 2015. 124. W. Wu, X. Xiao, X. Huang, and S. Yan. A multiphysics model for the in-situ stress analysis of the separator in a lithium-ion battery cell. Computational Materials Science, 83, pp.127-136, 2014. 125. S. Yan, J. Deng, C. Bae, S. Kalnaus, and X. Xiao, Orthotropic viscoelastic modeling of polymeric battery separator. Journal of The Electrochemical Society, 167(9), p.090530, 2020. 126. G. Zhao, H. Xi, and J. Yang. Transversely Isotropic Constitutive Model of the Polypropylene Separator Based on Rich–Hill Elastoplastic Constitutive Theory. Journal of Electrochemical Energy Conversion and Storage, 18(2), p.020911, 2021. 127. X. Zhang, E. Sahraei, and K. Wang. Deformation and failure characteristics of four types of lithium-ion battery separators. Journal of Power Sources, 327, pp.693-701, 2016. 37 128. C. Lee and C. W. Kim. Detailed Layered Nonlinear Finite Element Analysis for Lithium- Ion Battery Cells to Predict Internal Short Circuits Due to Separator Fractures under Hemisphere Indentation. Journal of The Electrochemical Society, 167(12), p.120511, 2020. 129. W. Xie, L. Wu, W. Liu, Y. Dang, A. Tang, and Y. Luo. Modelling electrolyte-immersed tensile property of polypropylene separator for lithium-ion battery. Mechanics of Materials, 152, p.103667, 2021. 38 Chapter 3 Orthotropic Linear Thermoviscoelastic Constitutive Modeling This chapter presents the development of an orthotropic linear viscoelastic material model for polymeric separators that accounts for the material anisotropy and rate dependence while accounting for temperature and electrolyte effects is presented. In this model, the temperature effect is introduced through the time-temperature superposition principle (TTSP). A time- temperature-solvent superposition method (TTSSM) is developed to model the behavior of the separator in electrolyte solutions based on the viscoelastic framework established in air. Analytical solutions for stress relaxation are developed and introduced for model verification. The model parameters are established for the selected separator material and presented in this chapter. The developed model is implemented in LS-DYNA® finite element (FE) package as a user-defined subroutine, which enables simulations with the thermal expansion/shrinkage behavior. The model predictions of the material anisotropy, temperature dependence, and solvent effect are validated against experimental data and the results are presented. The work presented in this chapter has been summarized and published [1]. 3.1 Orthotropic Linear Viscoelastic Model Development 3.1.1 Constitutive Model Overview The constitutive model developed in this work is based on the hereditary integral representation of the stress-strain response of linear viscoelastic materials. This hereditary integral is also known as the Boltzmann superposition integral [2-6]. The stiffness-based representation of the linear hereditary integral is expressed as: 𝑡 𝑑𝜀(𝜏) 𝜎(𝑡) = ∫ 𝐺(𝑡 − 𝜏) 𝑑𝜏 (3.1) 0 𝑑𝜏 where 𝜎 is the stress, 𝜀 is the strain and 𝐺(𝑡) is the stress relaxation stiffness matrix. To represent the components of the stress relaxation stiffness matrix, the stress-strain response of the separator is also expressed rheologically. In this work, this is done using the generalized Maxwell model which has been discussed in detail, in Chapter 2. The stress relaxation response of the generalized Maxwell model is expressed as a Prony series. Hence, each component of the stress relaxation stiffness matrix in Eqn.3.1 is usually expressed in terms of a Prony series expressed as: 𝑛 𝑡 𝐺(𝑡) = 𝐺∞ + ∑ 𝐺𝑖 ∙ 𝑒𝑥𝑝 [− ] (3.2) 𝜏𝑖 𝑖=1 39 where 𝐺∞ is the equilibrium or fully relaxed modulus, 𝐺𝑖 is the stress relaxation constant and 𝜏𝑖 is the relaxation time. 3.1.2 Numerical Intergation of the Linear Viscoelastic Hereditary Integral with a Kernel of Prony Series To implement the linear viscoelastic model for FE simulations involving arbitrary stress, strain and temperature histories, the hereditary integral has to be evaluated. A discretization algorithm for the evaluation of the linear viscoelastic hereditary integral with a kernel of a single exponential function was introduced by Puso and Weiss [7]. Using a similar approach, Yan et al [8] developed a discretization algorithm for the linear hereditary integral with a kernel of Prony series, implemented as an orthotropic linear viscoelastic model in FE packages. The simulations with this model compared well with analytical solutions with multi-step loading cases and with experimental results of non-isostress cases [8], proving that the algorithm is correct. This algorithm is used in this work as a framework for the orthotropic linear thermoviscoelastic model. The discretization algorithm starts from the analysis of Eqn.3.1 at a time 𝑡 + ∆𝑡: 𝑡+∆𝑡 𝑑𝜀 ′ 𝜎(𝑡 + ∆𝑡) = ∫ 𝐺(𝑡 + ∆𝑡 − 𝑡 ′ ) 𝑑𝑡 (3.3) 0 𝑑𝑡 ′ Evaluating Eqn.3.3 as a path integral allows the possibility of separating it into two terms that can be evaluated separately as follows: 𝑡 𝑡+∆𝑡 ′) 𝑑𝜀 ′ ′) 𝑑𝜀 ′ 𝜎(𝑡 + ∆𝑡) = ∫ 𝐺(𝑡 + ∆𝑡 − 𝑡 𝑑𝑡 + ∫ 𝐺(𝑡 + ∆𝑡 − 𝑡 𝑑𝑡 (3.4) 0 𝑑𝑡 ′ 𝑡 𝑑𝑡 ′ The second term in Eqn.3.4 can be solved approximately, using the Mean Value Theorem when ∆𝑡 is considered to be very small: 𝑡+∆𝑡 𝑑𝜀 ′ 𝜀(𝑡 + ∆𝑡) − 𝜀(𝑡) 𝑡+∆𝑡 Δ𝜀 ∫ 𝐺(𝑡 + ∆𝑡 − 𝑡 ′ ) ′ 𝑑𝑡 = ∫ 𝐺(𝑡 + ∆𝑡 − 𝑡 ′ ) 𝑑𝑡 ′ = ∙ 𝐺(∆𝑡) ∙ ∆𝑡 𝑡 𝑑𝑡 ∆𝑡 𝑡 ∆𝑡 = Δ𝜀 ∙ 𝐺(∆𝑡) (3.5) To evaluate the first term in Eqn.3.4, a specific kernel function was defined such that 𝐺(𝑡) was set initially to be a single term exponential function such that such that: 𝑡 𝐺(𝑡) = 𝐶1 ∙ 𝑒𝑥𝑝 (− ) (3.6) 𝜏 Hence, the first term in Eqn.3.4 becomes: 40 𝑡 𝑡 ′) 𝑑𝜀 ′ 𝑡 + ∆𝑡 − 𝑡 ′ 𝑑𝜀 ′ ∫ 𝐺(𝑡 + ∆𝑡 − 𝑡 𝑑𝑡 = ∫ 𝐶1 ∙ exp (− ) ′ 𝑑𝑡 0 𝑑𝑡 ′ 0 𝜏 𝑑𝑡 ∆𝑡 𝑡 𝑡 − 𝑡 ′ 𝑑𝜀 ′ = exp (− ) ∫ 𝐶1 ∙ exp (− ) ′ 𝑑𝑡 𝜏 0 𝜏 𝑑𝑡 ∆𝑡 = exp (− ) 𝜎(𝑡) (3.7) 𝜏 To represent the relaxation modulus over a broad time spectrum, using a single exponential function is insufficient. Hence, the algorithm was extended for the hereditary integral with a kernel of Prony series. This was based on the examination of the generalized Maxell model. From the model, the total stress is contributed to by the long-term equilibrium stress (𝜎∞ ) and the stresses in the individual Maxwell components (𝜎 𝑖 (𝑡)) expressed as: 𝑛 𝜎(𝑡) = 𝜎∞ + ∑ 𝜎 𝑖 (𝑡) (3.8) 𝑖=1 Applying the same mathematical operation as in Eqn.3.7, the first term of Eqn.3.4 is derived for a kernel of Prony series as: 𝑛 𝑡 𝑡 ′) 𝑑𝜀 ′ ∆𝑡 𝑡 − 𝑡 ′ 𝑑𝜀 ′ ∫ 𝐺(𝑡 + ∆𝑡 − 𝑡 𝑑𝑡 = 𝜎∞ + ∑ exp (− ) ∙ ∫ 𝐺𝑖 ∙ 𝑒𝑥𝑝 (− ) ′ 𝑑𝑡 0 𝑑𝑡 ′ 𝜏𝑖 0 𝜏𝑖 𝑑𝑡 𝑖=1 𝑛 ∆𝑡 = 𝜎∞ + ∑ exp (− ) ∙ 𝜎 𝑖 (𝑡) (3.9) 𝜏𝑖 𝑖=1 The flowchart for the implementation of this algorithm in LS-DYNA® user-defined subroutine is presented in Fig.3.1. At a step j, for time and strain increments ∆𝑡 and Δ𝜀 respectively, the stress increment is contributed to by the increment in the equilibrium relaxation stress and the sum of the incremental stresses in the individual Maxwell components. 𝑛 ∆𝜎 = ∆𝜎∞ + ∑ ∆𝜎 𝑖 (3.10) 𝑖=1 Furthermore, the stress increment in the individual Maxwell element is expressed as: ∆𝑡 ∆𝑡 ∆𝜎 𝑖 = 𝐺𝑖 exp (− ) ∙ ∆𝜀 + [exp (− ) ∙ 𝜎 𝑖 (𝑡) − 𝜎 𝑖 (𝑡)] (3.11) 𝜏𝑖 𝜏𝑖 In the user-defined subroutine, the stress in the individual Maxwell element and the total stress is updated and stored at every time step to be used in the next. 41 Figure 3.1 Flowchart for the implementation of evaluation of hereditary integral with Prony series [8]. 3.1.3 Plane Stress Stiffness Matrix To implement the discretization algorithm in section 3.1.2 for orthotropic materials, a stiffness matrix has to be introduced into the model. Polymeric battery separators are thin materials of thickness in the tens of micrometers and they are usually under a plane state of stress. For orthotropic linear elastic materials, the plane stress stiffness matrix is given as [9]: 𝐸1 𝑣12𝐸2 0 1 𝑬= ∙ [ 21 𝐸1 𝑣 𝐸2 0 ] (3.12) [1 − 𝑣12𝑣21] 0 0 [1 − 𝑣12 𝑣21 ]𝜏12 where 𝐸1 and 𝐸2 are the elastic moduli in the two principal in-plane directions. 𝜏12 is the in-plane shear modulus and 𝑣12 and 𝑣21 are the major and minor Poisons ratio respectively. The Poisons ratio and the elastic moduli for elastic materials are related by Betti’s reciprocal law such that [9]: 𝑣12 𝐸2 = 𝑣21 𝐸1 (3.13) Hence plane stress stiffness matrix in Eqn.3.12 above contains four independent engineering constants 𝐸1 , 𝐸2 , 𝑣12 , and 𝐺12. Furthermore, from the viscoelastic theory, the correspondence principle allows for the determination of solutions for viscoelastic problems, knowing its corresponding elastic solution 42 for the same loading, boundary conditions and geometry [2-6]. Hence, to obtain the plane stress stiffness matrix for an orthotropic linear viscoelastic material model, the elastic constants in Eqn.3.13 can be replaced with the corresponding viscoelastic functions according to the correspondence principle to give: 𝐺11 (𝑡) 𝑣12 𝐺22(𝑡) 0 1 𝑮= ∙ [𝑣12 𝐺22(𝑡) 𝐺22 (𝑡) 0 ] (3.14) [1 − 𝑣12 𝑣21 ] 0 0 [1 − 𝑣12 𝑣21]𝐺66(𝑡) Where 𝐺11 (𝑡), 𝐺22 (𝑡), and 𝐺66 (𝑡) are the relaxation moduli in the MD, TD and in-plane shear respectively, 𝑣12 and 𝑣21 are the major and minor Poison’s ratios. 3.1.4 Three-Dimensional Stiffness Matrix with Transverse Orthotropy In this work, a 3D stiffness matrix is also introduced for model implementation using solid elements. Here, transverse orthotropy is assumed in the 2-3 direction (i.e., 𝐸2 = 𝐸3 ) and the material response in the transverse direction (TD) is assumed to be the same in the through- thickness direction. Applying the elastic-viscoelastic correspondence principle for an elastic transversely orthotropic 3D stiffness matrix, the stiffness matrix in terms of viscoelastic functions becomes: 𝐶11 𝐶12 𝐶12 0 0 0 𝐶12 𝐶22 𝐶23 0 0 0 𝐶23 𝐶22 0 0 𝑮 = 𝐶12 0 (3.15) 0 0 0 𝐶44 0 0 0 0 0 0 𝐶66 0 [ 0 0 0 0 0 𝐶66 ] Where the components of the stiffness matrix 𝐶𝑖𝑗 are defined as: (1 − 𝑣232 ) ∙ 𝐺11 (𝑡) 𝐶11 = ∆ (𝑣12 + 𝑣12 ∙ 𝑣23 ) ∙ 𝐺22 (𝑡) 𝐶12 = ∆ (1 − 𝑣12 ∙ 𝑣21 ) ∙ 𝐺22 (𝑡) 𝐶22 = ∆ 𝐺22 (𝑡) 𝐶44 = 2 ∙ (1 + 𝑣23 ) (𝑣23 + 𝑣12 ∙ 𝑣21 ) ∙ 𝐺22 (𝑡) 𝐶23 = ∆ 43 𝐶66 = 𝐺66 (𝑡) 2 Also, ∆= [1 − 𝑣23 − 2 ∙ 𝑣12 ∙ 𝑣21 − 2 ∙ 𝑣12 ∙ 𝑣21 ∙ 𝑣23 ]. 𝑣23 is the poisons ratio in the 2-3 plane. Since transverse orthotropy is assumed here, it is also vital to note that 𝑣23 = 𝑣32. 3.2 Orthotropic Linear Thermoviscoelastic Material Characterization 3.2.1 Experimental Procedure The separator investigated is Celgard®2400, a porous PP material. Celgard®2400 is a single-layer film with a thickness of 25µm. The received separator was in the form of a roll of 255 mm in width. The in-plane viscoelastic responses of the separators were measured along the MD, TD, and in-plane shear using a 45° off-axis specimen as discussed in [10]. The nominal dimensions of each specimen were 45 mm in length and 5 mm in width. These specimens were cut using a razor blade. The sample gage length was ~17mm. The in-plane viscoelastic property of Celgard ®2400 was investigated in air and dimethyl carbonate (DMC). DMC is a common solvent for LIB electrolytes. It has been shown that the Young’s modulus of Celgard®2400 measured in DMC was about 95% of the value measured in an electrolyte solution of 1.1 M LiPF6 EC/DMC in a 1:1 volume ratio for both MD and TD [11]. Therefore, DMC is a good replacement for LIB electrolytes in the investigation of the electrolyte effect. (a) (b) Figure 3.2 (a) RSA-G2 dynamical mechanical analyzer testing set up, (b) external heating of DMC in water bath, (c) DMC transferred into solvent bath, and (d) closed solvent bath. 44 Figure 3.2 (Cont’d) (c) (d) Uniaxial stress relaxation tests were carried out using an RSA-G2 rheometer in tensile mode, as shown in Fig.3.2. For the test in air, the separator was tested over a temperature range of 20°C-110°C whereas, for the tests in DMC, the range was 20°C-40°C. The upper testing temperature in DMC was limited by the high volatility of the solvent. To saturate the sample during testing in DMC, a solvent bath was used as shown in Fig.3.2(a). The samples were clamped into the tensile fixture using a torque screwdriver to maintain the clamping force to ensure that there was no sample slippage during the tests. For the tests in DMC above ambient temperature, to reduce the evaporation during testing, the solvent was preheated. This was accomplished by heating the solvent in a water bath and monitoring its temperature with a thermocouple, as shown in Fig.3.2(b). When the solvent reached the desired testing temperature [30°C, 35°C, 40°C], it was poured into the testing fixture where the sample already resides as shown in Fig.1c. The heating chamber was closed and the temperature was equilibrated at the desired value. Regardless of whether it was in air or in DMC, the sample was kept at the specified temperature for 4min before testing. There was no observed difference in the relaxation behavior for soaking and heating times of 4min and 10min. A preload force of ~3g was applied to ensure the sample was straight and the gap was locked before starting each test. Each stress relaxation test was run for 20 mins at the given temperature level. 45 Isothermal tests were carried out in all three material orientations, in air, and in DMC according to the method described above. To validate the model, non-isothermal stress relaxation tests were performed. The non-isothermal stress relaxation tests were carried out in air using the samples in the MD with multi-step temperature histories. A fresh sample was used for each test and all tests were repeated two to three times for replicability and consistency of the results. The results presented here are the averaged results. Table 3.1 below provides a summary of the experimental procedures and their purpose. Table 3.1 Summary of the experimental procedures Experiment Specimen Test Conditions Purpose MD, TD, 0.2-1.2% strain levels at 20°C- Determination of Stress relaxation 45° off-axis 22°C and 80°C in air and at linear viscoelastic 20°C-22°C in DMC strain limit MD, TD, 0.2-0.3% strain levels at Establishment of the 45° off-axis temperatures ranging from TTSP and 20°C-110°C in air determination of Stress relaxation temperature- dependent model parameters. MD, TD, 0.4% strain level at Establishment of the 45° off-axis temperatures ranging from TTSSM and Stress relaxation 20°C-40°C in DMC determination of solvent-dependent model parameters. MD, TD 0.001, 0.01mm/s at 30°C, 40°C Validation under Tensile 50°C, 60°C in air iso-thermal conditions MD 30°C to 40°C, 60°C to 70°C, Validation under Stress relaxation and 30°C to 60°C with non iso-thermal T=10°C at 0.2% strain level condition MD, TD 0.067-3MPa at 3°C/min ramp Validation under Iso-stress rate non iso-thermal condition 3.2.2 Determination of Linear Viscoelastic Strain Limit To ensure that the stress relaxation experiments are conducted within the range of the linear viscoelastic response of the material, the linear viscoelastic strain limit needs to be determined. In linear viscoelasticity, the relaxation modulus is independent of strain levels. Hence, the linear viscoelastic limit was determined by measuring the stress relaxation responses at different strain 46 levels, as shown in Fig.3.3. From Fig.3.3, the strain limit for the TD is 0.3% in air and 0.4% in DMC. The strain limit for the linear stress relaxation response in air was determined at the lowest test temperature 20-22°C and at 80°C, and in DMC at 20-22°C. It has been observed that the strain limit increases as temperature crosses the glass transition temperature (𝑇𝑔 ) from below [12]. The reported 𝑇𝑔 of Celgard®2400 is -15°C in air [13]. In this work, all stress relaxation tests were conducted at temperatures higher than the 𝑇𝑔 . To observe the possible increase of the strain limit, the measurements were also determined at 80°C in air in the TD (Fig 3.4) and the strain limit was found indeed to be higher than the value at 20°C. When testing in DMC, molecules diffuse into the amorphous phase of the separator and act as plasticizers and lubricants[14]. It is expected that the strain limit would increase due to plasticization as reported in [15]. Table 3.2 summarizes the values of the strain limit for various conditions. (a) (b) Figure 3.3 Determination of the linear viscoelastic strain limit using stress relaxation experiments at 20-22°C. (a) TD in air (b) TD in DMC (c) MD in air (d) MD in DMC (e) 45° off- axis in air (f) 45° off-axis in DMC. 47 Figure 3.3 (Cont’d) (c) (d) (e) (f) Figure 3.4 Determination of the linear viscoelastic strain limit using stress relaxation experiments at 80°C for TD in air. 48 Table 3.2 Strain limit of the linear stress relaxation response in air and in DMC Strain limit of the linear stress relaxation response (%) Environment MD TD 45° Off-axis In air (20-22°C) 0.2 0.3 0.3 In air (80°C) - 0.4 - In DMC (20-22°C) 0.4 0.4 0.6 3.2.3 Stress relaxation at different temperatures Stress relaxation tests were conducted in air and in DMC at temperatures above the ambient temperature and within the strain levels listed in Table 3.2. The test results in air and in DMC are shown in Figs.3.5 and 3.6 respectively. The stress relaxation modulus is plotted vs log time and the results show the reduction in the relaxation with increasing temperature and the presence of the solvent. This is consistent with the viscoelastic theory. (a) (b) (c) Figure 3.5 Stress relaxation curves in air for samples cut along the (a) MD, (b) TD, and (c) off- axis 45°. 49 (a) (b) (c) Figure 3.6 Stress relaxation curves in DMC for samples cut along the (a) MD, (b) TD, and (c) off-axis 45° direction. 3.2.4 Determination of In-Plane Shear Relaxation Modulus Due to the difficulties involved in direct measurement of the in-plane shear properties of thin films, the time-dependent in-plane shear relaxation modulus of the separator in air and in DMC was determined using the 45° off-axis specimen and then computed analytically as described in [10]. This is based on a well-established method for determining the shear properties of unidirectional composites [9] : 1 𝑚2 (𝑚2 − 𝑛 2 𝑣12 ) 𝑛 2 (𝑛2 − 𝑚2 𝑣21 ) 𝑚2 𝑛 2 = + + (3.16) 𝐸𝑥 𝐸1 𝐸2 𝐺12 where m= cos θ, n= sin θ, and E1, E2, and G12 are the longitudinal, transverse, and shear elastic moduli. 50 According to the elastic-viscoelastic correspondence principle [2-6], the elastic constants in Eqn.3.16 can be substituted with the viscoelastic constants. For θ=45°, the time-dependent in- plane shear relaxation modulus was calculated using the relationship: −1 4 (1 − 𝑣12 ) (1 − 𝑣21 ) 𝐺66 (𝑡) = [ − − ] (3.17) 𝐺45° (𝑡) 𝐺11 (𝑡) 𝐺22 (𝑡) where G45° (t), G11(t), and G22(t) are the measured stress relaxation of the 45° off-axis, MD, and TD samples respectively. The major Poisson’s ratio 𝜈12 in DMC was assumed to be 0.17, the same value as in air. The minor poisons ratio 𝜈21 is dependent on v12, G22 , and G11 through Betti’s reciprocal law [9]. 𝜈21 was calculated analytically using the formula: 𝜈12 𝐺22 (𝑡) 𝜈21 = (3.18) 𝐺11 (𝑡) Furthermore, the stress relaxation curves for in-plane shear in air and in DMC are plotted in a logarithmic time scale and are presented in Fig.3.7. The relaxation curves followed the same trend, expressed by the reduction of the stress relaxation modulus with increasing temperature and presence of electrolyte solutions or solvent. (a) (b) Figure 3.7 In-plane shear relaxation modulus in (a) air and (b) DMC. 3.3 Temperature and Electrolyte Effects 3.3.1 Implementation of the Time-temperature superposition principle (TTSP) The TTSP can be used as a methodology to systematically characterize the time and temperature-dependent properties of viscoelastic materials and as a model to describe its thermo- viscoelastic behavior [2,16]. Following the TTSP, the stress relaxation of the material as a function 51 of time at different temperatures is expressed by a master curve and shifting factors. This master curve is built by shifting individual isothermal relaxation curves relative to the relaxation curve at a reference temperature. These isothermal curves are commonly shifted horizontally in a logarithmic time scale, as shown in Fig.3.8. Occasionally, vertical shifts in the logarithmic relaxation modulus scale in conjunction with the horizontal shifts, are used. This is done in order to produce the best master curve possible. For characterizing the temperature-dependent response of the separator in the linear range and for large times, only horizontal shifts were used in the implementation of the TTSP as they were sufficient to produce suitable master curves. In terms of the horizontal shift, the stress relaxation modulus as a function of time and temperature is expressed as [4]: 𝑡 𝐺(𝑡, 𝑇) = 𝐺 ( ,𝑇 ) (3.19) 𝑎 𝑇 (𝑇, 𝑇0 ) 0 Where 𝑎 𝑇 (𝑇, 𝑇0 ) is the horizontal shifting factor in the time scale which is dependent on the current (𝑇) and reference (𝑇0 ) temperatures. The TTSP was applied to the experimental data generated for all three material orientations with the lowest testing temperature as the reference temperatures. Fig.3.8 presents the master curves for MD, TD, and in-plane shear in air and in DMC. (a) (b) Figure 3.8 Application of the TTSP in air and in DMC: (a) Master curve for MD sample in air (b) Master curve for TD sample in air (c) Master curve for in-plane shear in air (d) Master curve for MD sample in DMC (e) Master curve for TD sample in DMC (f) Master curve for in-plane shear in DMC. 52 Figure 3.8 (Cont’d) (c) (d) (e) (f) The master curves together with the shiting factors describe the viscoelastic behavior of the material over a wide range of time and temperatures. In literature, the temperature dependency of the shifting factors has been explained with the William-Landel-Ferry (WLF) model and the Arrhenius equation [3]. The WLF model is based on the free volume theory for amorphous polymers whereas the Arrhenius equation is for the rate of thermally activated processes. Both models were examined for fitting the horizontal shifting factors. Firstly, the shifting factors were fitted using the WLF model using a least square regression and presented in Fig 3.8. The results show that the horizontal shifting factor (𝑎 𝑇 (𝑇, 𝑇0 )) data deviated from the WLF model at T > 70°𝐶 as detailed in Table 3.3. Furthermore, Williams et al [16], have reported that the WLF model is not expected to hold at temperatures 𝑇 > 𝑇𝑔 + 100°𝐶. In the current case, the 𝑇𝑔 of PP is -15°C in air and thus, the deviation from the WLF model is justified. 53 (a) (b) (c) Figure 3.9 Shifting factor fitting using the WLF model for master curves in air (a) MD, (b) TD, and (c) In-plane shear. Table 3.3 Predicted shifting factors using the WLF model compared with experimental shifting factors Machine Direction Transverse Direction In-Plane Shear (MD) (TD) T Experiment Predicted Experiment Predicted Experiment Predicted (ºC ) 20 0 0 0 0 0 0 30 -1.15 -1.18 -1.90 -1.89 -1.10 -1.17 40 -2.67 -2.42 -3.50 -3.67 -3.25 -2.45 50 -3.97 -3.73 -5.70 -5.32 -4.25 -3.86 60 -5.70 -5.09 -7.42 -6.86 -6.00 -5.42 70 - - -8.65 -8.31 - - 80 -7.60 -8.07 -9.70 -9.67 -7.75 -9.12 90 - - -10.50 -10.95 - - 100 -9.25 -11.38 -11.50 -12.16 -9.70 -13.85 110 -9.98 -13.19 -12.10 -13.31 -10.10 -16.74 54 Hence, for the model development and determination of the temperature-dependent parameters, the Arrhenius equation was used to fit the shifting factors. From Table 3.4 it is evident that the shifting factors predicted using the Arrhenius equation are comparable to the experimental values. The Arrhenius theory considers the temperature dependence of reaction rates. Assuming the stress relaxation process is thermally activated, it would conform to the kinetic rate theory. The Arrhenius equation for relating the shift factor to temperature is given by [17]: ∆𝐻 1 1 𝐿𝑜𝑔 (𝑎 𝑇_𝑎𝑖𝑟 ) = [ − ] (3.20) 2.303𝑅 𝑇 𝑇0 where ∆𝐻 is the activation energy, R is the universal gas constant, T and T0 are the current and reference temperatures, respectively. Here R=8.3145 J/(mol-K). The shifting factors for the master curves in air were fitted with the Arrhenius equation using the least square regression. The coefficients of determination R2 of these fittings were in the range of 0.996-0.998, indicating excellent fitting (Fig.3.10). Table 3.4 compares the experimental shifting factors and shifting factors calculated by the Arrhenius equation. From Table 3.4, it is evident that the temperature-dependent behavior of the separator material follows the Arrhenius equation over the entire temperature range. (a) (b) Figure 3.10 Shifting factor fitting using the Arrhenius equation for master curves in air (a) MD, (b) TD, and (c) In-plane shear. 55 Figure 3.10 (Cont’d) (c) Table 3.4 Experimental shifting factors in air and the predicted values using the Arrhenius equation Machine Direction (MD) Transverse Direction In-Plane Shear (TD) T Experiment Predicted Experiment Predicted Experiment Predicted (ºC ) T0 0 0 0 0 0 0 30 -1.15 -1.43 -1.90 -1.82 -1.10 -1.48 40 -2.67 -2.78 -3.50 -3.52 -3.25 -2.86 50 -3.97 -4.04 -5.70 -5.12 -4.25 -4.16 60 -5.70 -5.22 -7.42 -6.62 -6.00 -5.38 70 - - -8.65 -8.04 - - 80 -7.60 -7.39 -9.70 -9.37 -7.75 -7.61 90 - - -10.50 -10.63 - - 100 -9.25 -9.32 -11.50 -11.82 -9.50 -9.60 110 -9.98 -10.22 -12.10 -12.96 -10.10 -10.52 Assuming that for the same temperature range (20°C-110°C), the temperature-dependent behavior of the separator in DMC will also follow the same trend as in air, the Arrhenius equation was also used to fit the shifting factors for the master curve in DMC (aTDMC ). The fittings in all three orientations achieved R2>0.998 (Fig.3.11). Table 3.5 compares the experimental values and the prediction by the Arrhenius equation. The Arrhenius equation provided good fittings for the available experimental data obtained in DMC. Furthermore, Table 3.6 presents the values of the 56 activation energy in the Arrhenius equation determined by curve fitting for the shifting factors in air and in DMC. (a) (b) (c) Figure 3.11 Shifting factor fitting using the Arrhenius equation for master curves in DMC (a) MD, (b) TD, and (c) In-plane shear. Table 3.5 Experimental Shifting factors in DMC and the predicted values using the Arrhenius equation Machine Direction Transverse Direction In-Plane Shear (MD) (TD) T Experiment Predicted Experiment Predicted Experiment Predicted (ºC ) T0 0 0 0 0 0 0 30 -1.10 -1.06 -1.50 -1.46 -1.00 -0.98 35 -1.70 -1.62 -2.40 -2.33 -1.50 -1.45 40 -2.10 -2.17 -3.10 -3.17 -1.85 -1.89 57 Table 3.6 The activation energy in Arrhenius equation (Eqn.3.20) Parameters Machine Direction Transverse In-Plane Shear (MD) Direction (TD) ∆𝐻air (KJ) 243.88 309.30 251.11 T0 (°C) 20 20 20 ∆𝐻𝐷𝑀𝐶 (KJ) 201.49 311.60 166.78 T0 (°C) 21 22 20 3.3.2 Time-temperature-solvent superposition method (TTSSM) To understand the mechanical response of the separator in its working environment, the electrolyte effect needs to be taken into account. Electrolyte solutions have been found to remarkably weaken the mechanical performance [11] and thermal stability [18] of LIB separators. It has been observed that when being immersed in DMC, the tensile modulus of Celgard ®2400 decreased to 48.5% and 87.7% in the MD and TD respectively in terms of ratio to the measured modulus in air [11]. Hence, a model based on the behavior of the polymeric separator in air will lead to an overprediction in its mechanical response under stresses in batteries. Due to the difficulties in testing for the thermo-mechanical behavior in solvent at high temperatures, it is inherently important to develop a model to predict the mechanical response of the separator in electrolyte solutions. Krauklis et al [19] developed a time-temperature-plasticization superposition principle (TTPSP) to predict the long-term viscoelastic behavior of plasticized amorphous polymers below 𝑇𝑔 using short-term experimental data. However, in thermal ramp scenarios, polymeric battery separators experience temperatures well above 𝑇𝑔 , making the TTPSP not applicable. To overcome this limitation, a TTSSM is developed in this work to utilize the viscoelastic response of the separator in air as a framework for predicting its behavior in solvent at temperatures above 𝑇𝑔 . The TTSSM is an extension of the TTSP. It is based on the assumption that the temperature-dependent mechanical behavior of the separator in solvent at a temperature value (𝑇𝑠𝑜𝑙 ) is equivalent to the behavior in air at another temperature value (𝑇𝑎𝑖𝑟 ). In the current work, this principle was established, using the TTSP master curves in air and the stress relaxation curves in solvent. 58 Figure 3.12 Time-temperature-solvent superposition method (TTSSM). The process of applying the time-temperature-solvent superposition (TTSS) generates a time-temperature-solvent shift factor (aTTS ). aTTS describes the temperature-dependent behavior of the separator in solvent from its known behavior in air. In the current work, aTTS is decoupled into two shifting factors. The first is the shifting factor in DMC (𝑎 𝑇_𝐷𝑀𝐶 ) which accounts for the temperature effect on the material in the plasticized state. ATDMC is determined by horizontally shifting the experimental isothermal stress relaxation curves in DMC to create a TTSP master curve in DMC, as discussed in 3.3.1. The second shifting factor is the air-solvent shift (𝑎𝑎_𝑠 ) which is determined by shifting the TTSP master curve in DMC to the TTSP master curve in air, as illustrated in Fig.3.12. 𝑎𝑎_𝑠 accounts for the increase in the free volume, structural realignment, and decrease in the relaxation modulus of the separator due to the presence of the solvent. Hence, 𝑎𝑎_𝑠 introduces the plasticizing effect to the master curve in air, making it possible to use the master curve in air as a framework for predicting the electrolyte effect on the mechanical properties of the separator. The proposed decoupled aTTS is expressed as: 𝐿𝑜𝑔 (𝑎𝑇𝑇𝑆 ) = 𝐿𝑜𝑔 (𝑎𝑎_𝑠 ) + 𝐿𝑜𝑔 (𝑎 𝑇_𝐷𝑀𝐶 ) (3.21) 𝑎 𝑇𝑇𝑆 = 𝑎𝑎_𝑠 ∙ 𝑎 𝑇_𝐷𝑀𝐶 (3.22) Attempts to fit 𝑎𝑎_𝑠 data using the WLF and Arrhenius-like equations were made but were not successful. On the other hand, a linear log relationship was found to be sufficiently accurate in fitting the 𝑎𝑎_𝑠 (T, T0 ) function. 59 𝐿𝑜𝑔 (𝑎𝑎_𝑠 ) = 𝑄1 ∙ (𝑇 − 𝑇0 ) + 𝑄2 (3.23) In Eqn.3.23, T is the current temperature, T0 is the reference temperature of the super master curve (time-temperature-solvent master curve), Q1 as the slope which is dependent on the shape of the master curves representing the viscoelastic behavior in air and in solvent, and Q 2 is the air-solvent shift when the current temperature is equal to the reference temperature (T=T0 ). The air-solvent shift was fitted using a least square regression method according to Eqn. 3.23 and the results are presented in Fig.3.13. The fittings in all three material orientations achieved R2 > 0.99 showing excellent fittings. (a) (b) (c) Figure 3.13 Air-solvent shifting factor fitting in the (a) MD, (b) TD, and (c) In-plan shear. Implementing aTTS on the experimental data leads to the formation of an air-solvent super master curve. This super master curve, otherwise known as the time-temperature-solvent master curve will account for the solvent effect on the linear viscoelastic behavior of polymeric battery separators at different temperature values. Fig.3.14 presents the time-temperature-solvent 60 superposition (TTSS) process for the three material orientations. Figs.3.13(b, d, and f) show the super master curves for the MD, TD, and in-plane shear respectively. (a) (b) (c) (d) (e) (f) Figure 3.14 Implementation of the time-temperature-solvent shift 𝑎𝑎_𝑠 . (a-b) MD, (b-c) TD, and (e-f) in-plane shear. 61 3.4 Analytical Solutions for Stress Relaxation 3.4.1 Single Step Stress Relaxation at Constant Temperature The TTSP master curve for the relaxation modulus was fitted with the generalized Maxwell model in Prony series. The mathematical expression is given by [4]: 𝑛 𝑡 𝐺(𝑡) = 𝐺∞ + ∑ 𝐺𝑖 ∙ 𝑒𝑥𝑝 (− ) (3.24) 𝜏𝑖 𝑖=1 where 𝐺∞ and 𝐺I are the Prony series coefficients, t is the time, and 𝜏𝑖 are the characteristic relaxation times, respectively. To include the temperature effect, a reduced time t’ is introduced such that: 𝑡 𝑡′ = (3.25) 𝑎 𝑇_𝑎𝑖𝑟 Replacing t in Eq.3.24 with t’, we have the expression for the relaxation modulus with temperature dependence in air: 𝑛 𝑡 𝐺(𝑡) = 𝐺∞ + ∑ 𝐺𝑖 ∙ 𝑒𝑥𝑝 (− ) (3.26) 𝑎 𝑇𝑎𝑖𝑟 ∙ 𝜏𝑖 𝑖=1 Similarly, the reduced time t’ can include the time-temperature-solvent effect in the form: 𝑡 𝑡′ = (3.27) 𝑎𝑎_𝑠 ∙ 𝑎 𝑇_𝐷𝑀𝐶 Hence, the expression for the relaxation modulus with temperature dependence and solvent effect is expressed as: 𝐺(𝑡) = 𝐺∞ + ∑𝑛𝑖=1 𝐺𝑖 ∙ 𝑡 𝑒𝑥𝑝 (− 𝑎 ) (3.28) 𝑎_𝑠 ∙ 𝑎𝑇_𝐷𝑀𝐶 ∙ 𝜏𝑖 Equations 3.24, 3.26, and 3.28 are based on the same framework, i.e., the Prony series for the master curve in air obtained by TTSP. Furthermore, the relationship between the shifting factor in solvent and the shifting factor in air is assumed to be: 𝐿𝑜𝑔(𝑎 𝑇_𝐷𝑀𝐶 ) = Q 3 ∙ Log (𝑎 𝑇_𝑎𝑖𝑟 ) (3.29) Simplifying Eqn.3.29 by eliminating the logarithmic terms leads to: Q3 𝑎 𝑇_𝐷𝑀𝐶 = (𝑎 𝑇_𝑎𝑖𝑟 ) (3.30) where Q 3 , is a material constant. Since both 𝑎 𝑇_𝐷𝑀𝐶 and 𝑎 𝑇_𝑎𝑖𝑟 have been fitted with the Arhenius equation, we have: 62 ∆𝐻DMC Q3 = (3.31) ∆𝐻air From Eqn.3.31, Q 3 is the ratio of the activation energies for the stress relaxation process in solvent and in air. A value of Q 3 <1 indicates that the activation energy of the relaxation process is reduced in the presence of electrolyte solutions and vice versa. The modified equation for the relaxation modulus as a function of time in solvent can be expressed as: 𝑛 𝑡 𝐺(𝑡) = 𝐺∞ + ∑ 𝐺𝑖 ∙ 𝑒𝑥𝑝 (− Q3 ) (3.32) 𝑖=1 𝑎𝑎_𝑠 ∙ (𝑎 𝑇_𝑎𝑖𝑟 ) ∙ 𝜏𝑖 From Eqn. 3.32, we can simply revert the process to consider only the temperature effect in air (Eqn.3.26). In this case, there is no air-solvent shift (𝐿𝑜𝑔 (𝑎𝑎_𝑠 ) = 0, 𝑎𝑎_𝑠 = 1 ), and Q 3 = 1. Furthermore, the stress relaxation as a function of time can be calculated from the relaxation moduli expressed above by multiplying by the desired constant strain level, 𝜀0 . Therefore, the analytical solutions for the stress relaxation accounting for temperature dependence and temperature-solvent effect are given by Eqn.3.33 and Eq.3.34 respectively below: 𝑛 𝑡 𝜎(𝑡) = 𝜀0 ∙ [𝐺∞ + ∑ 𝐺𝑖 𝑒𝑥𝑝 (− )] (3.33) 𝑎 𝑇_𝑎𝑖𝑟 ∙ 𝜏𝑖 𝑖=1 𝑛 𝑡 𝜎(𝑡) = 𝜀0 ∙ [𝐺∞ + ∑ 𝐺𝑖 𝑒𝑥𝑝 (− 𝑄3 )] (3.34) 𝑖=1 𝑎𝑎_𝑠 ∙ (𝑎 𝑇_𝑎𝑖𝑟 ) ∙ 𝜏𝑖 3.4.2 Step loading at constant temperature According to the linear viscoelastic theory [3], the analytical expression for the total stress 𝜎(𝑡) at a given time t, when the constant strain for stress relaxation is added incrementally, is given by: 𝜎(𝑡) = ∆𝜀1 ∙ 𝐺(𝑡 − 𝑡1 ) + ∆𝜀2 ∙ 𝐺(𝑡 − 𝑡2 ) + ∙∙∙∙∙ +∆𝜀𝑛 ∙ 𝐺(𝑡 − 𝑡𝑛 ) (3.35) Eqn. 3.35 gives the basis the Boltzmann superposition principle, where 𝐺(𝑡 − 𝑡1 ) is the stress relaxation modulus, and ∆𝜀1, ∆𝜀2 and ∆𝜀𝑛 are the strain increments added at times 𝑡1 , 𝑡2 , and 𝑡𝑛 . The stress relaxation modulus can be represented by Eqn.3.26 for cases with temperature effect in air and Eqn.3.32 for cases in DMC. For multistep loading cases for stress relaxation at constant temperature, the analytical solution is given by Eqn.3.35 63 3.4.3 Step temperature history at constant strain level For the case where the temperature loading is applied stepwise and the temperature values 𝑇1, 𝑇2, up to 𝑇𝑛 for n steps are added at times 𝑡1 , 𝑡2 , up to 𝑡𝑛 . To develop the analytical solution for stress relaxation with a step temperature history, the reduced time has to be evaluated for the time intervals at which the temperature steps are held. For linear viscoelastic response with temperature dependence introduced via the TTSP, the reduced time is defined mathematically as: 𝑡 𝑑𝜏 𝜌(𝑡) = ∫ (3.36) 0 𝑎𝑇 To calculate the reduced time and vertical shifts, the total time is broken down into intervals: 𝑡 ∈ [𝑡1 , 𝑡2 ], 𝑡 ∈ [𝑡2 , 𝑡3 ], … , and 𝑡 > 𝑡𝑛 , where n is the total number of time segments. For each time interval, the reduced time is calculated as follows. For the interval 𝑡 ∈ [𝑡1 , 𝑡2 ], the reduced time becomes: 𝑡 𝑑𝜏 𝑡 − 𝑡1 𝜌(𝑡) = ∫ = (3.37) 𝑡1 𝑎 𝑇1 𝑎 𝑇1 For the interval 𝑡 ∈ [𝑡2 , 𝑡3 ], the reduced time is: 𝑡2 𝑡 𝑑𝜏 𝑑𝜏 𝑡2 − 𝑡1 𝑡 − 𝑡2 𝜌(𝑡) = ∫ +∫ = + (3.38) 𝑡1 𝑎 𝑇1 𝑡2 𝑎 𝑇2 𝑎 𝑇1 𝑎 𝑇2 In cases that involve more temperature steps, the reduced time can be evaluated for the intervals at which the temperature values are added following the procedure in Eqns. 3.37 and 3.38. The reduced time at the time interval 𝑡 > 𝑡𝑛 after the final temperature step has been added is given by: 𝑡2 𝑡3 𝑡 𝑑𝜏 𝑑𝜏 𝑑𝜏 𝑡2 − 𝑡1 𝑡3 − 𝑡2 𝑡 − 𝑡𝑛 𝜌(𝑡) = ∫ +∫ + ⋯+∫ = + + ⋯+ (3.39) 𝑡1 𝑎 𝑇1 𝑡2 𝑎 𝑇2 𝑡𝑛 𝑎 𝑇3 𝑎 𝑇1 𝑎 𝑇2 𝑎 𝑇𝑛 Hence, the analytical solution for the stress relaxation as a function of time and temperature with step history is given as: 𝑛 𝜌(𝑡) 𝜎(𝑡, 𝑇) = 𝜀0 ∙ {𝐺∞ + ∑ 𝐺𝑖 ∙ 𝑒𝑥𝑝 [− ]} (3.40) 𝜏𝑖 𝑖=1 With temperature change, the thermal strains due to thermal expansion have to be accounted for and introduced into the analytical solution. The thermal strain is expressed as: 𝜀𝑡ℎ = 𝛼 ∙ (𝑇 − 𝑇0 ) (3.41) 64 where 𝛼 is the coefficient of thermal expansion (CTE). Hence the analytical solution for stress relaxation accounting for thermal expansion effect takes the form: 𝑛 𝜌(𝑡) 𝜎(𝑡, 𝑇) = [𝜀0 − 𝜀𝑡ℎ ] ∙ {𝐺∞ + ∑ 𝐺𝑖 ∙ 𝑒𝑥𝑝 [− ]} (3.42) 𝜏𝑖 𝑖=1 The analytical solutions (Eqns.3.40 and 3.42) developed in this section are applicable for stress relaxation with step temperature increment in air (𝑎 𝑇_𝑎𝑖𝑟 ) and in DMC (𝑎 𝑇_𝐷𝑀𝐶 ). 3.5 Prony Series Fitting and Summary of TTSSM Parameters The stress relaxation modulus of the separator investigated in this work was expressed in terms of a Prony series. The relaxation constants were determined using the least squares method calculated using a MATLAB (MathWorks, Natick, MA, USA) code. The Prony series parameters for the separator material in MD, TD, and in-plane shear are presented in Table 3.7. To show the quality of the Prony series coefficients generated from the least square fitting, the curves for the fitted data and the master curves from the experimental data are compared in a log time scale for all three material orientations in Fig.3.15. Table 3.7 Prony series parameters Machine Direction (MD) Transverse Direction (TD) In-Plane Shear Relaxation Relaxation Relaxation Relaxation Relaxation Relaxation Coefficient time 𝝉𝒊 (s) Coefficient time 𝝉𝒊 (s) Coefficient time 𝝉𝒊 (s) 𝑮𝒊 (Pa) 𝑮𝒊 (Pa) 𝑮𝒊 (Pa) 2.08 × 107 - 9.18 × 107 - 9.18 × 106 - 8 7 7 3.64 × 10 1 7.61 × 10 1 2.24 × 10 1 8 4.22 × 10 102 8.47 × 10 7 102 3.16 × 10 7 102 8 4.02 × 10 104 9.16 × 10 7 104 2.79 × 10 7 104 3.01 × 108 106 9.95 × 107 106 2.67 × 107 106 6.47 × 107 107 1.02 × 108 108 6.92 × 106 107 8 1.89 × 10 108 1.04 × 10 8 1010 1.46 × 10 7 108 1.33 × 108 1010 8.15 × 107 1012 9.83 × 106 1010 8.71 × 107 1012 8.20 × 107 1014 3.85 × 106 1012 65 (a) (b) Figure 3.15 Prony series fitting for the stress relaxation master curves in log time scale for (a) MD and TD, and (b) In-plane shear. Furthermore, the TTSSM parameters are summarized in Table 3.8. It is worth noting that the Q3 values obtained by fitting agreed well with the values estimated by Eqn.3.31. Fig 3.16 below presents the least square regression fittings for the material parameter Q3. (a) (b) Figure 3.16 Curve fittings for the material parameter Q3 in the (a) MD, (b) TD, and (c) In-plane shear. 66 Figure 3.16 (Cont’d) (c) Table 3.8 TTSSM Parameters Parameters Machine Direction Transverse In-Plane Shear (MD) Direction (TD) Q1 0.016 -0.009 0.01 Q2 -4.11 -4.81 -3.70 Q3 0.83 1.01 0.66 ∆𝐻 0.83 1.01 0.66 Q3 = ∆𝐻DMC air 3.6 Simulation Details 3.6.1 Shell Element Simulations The developed model is validated in this work for uniaxial loading conditions. One-element simulations are used for model verification and validation. This is a standard practice among finite element (FE) code developers to test developed subroutines. A shell element is used for all simulations involving model validation. This is due to the large discrepancy in the thickness of considered thin film samples in comparison to the in-plane dimension. A square shell element of size 5mm with node numbers, constraints and coordinate system as shown in Fig.3.17 was used. The boundary conditions were defined as follows: node #1 was fixed in all translational degrees of freedom, while node #2 was constrained in the y-direction but allowed to translate in the x- direction, node #3 was constrained in the x-direction but allowed to move in the y-direction and node #4 has no translational constraints. Table 3.9 below summarizes the boundary conditions applied at the nodes with 1 meaning “constrained” and 0 meaning “allowed to translate”. For 67 simulations in the MD, the loadings are applied at nodes #2 and #4 (x-direction). For TD simulations, the loadings are applied at nodes #3 and #4 (y-direction). Furthermore, for 45° off- axis simulations, the loadings were applied in either the x-direction (nodes #2 and #4) or the y- direction (nodes #3 and #4) while the material coordinate was rotated in the anti-clockwise direction (𝛽 = 45°). For the 45° off-axis simulations, applying the loadings in the x or y-direction produced the same results. The appropriate temperature history was applied at all four nodes of the shell element. Explicit time integration was used in all simulations with the time, scaled so that 1ms will correspond to 1s in real-time. Figure 3.17 Shell element with boundary conditions. Table 3.9 Boundary conditions for shell element Node Number X - Direction Y - Direction 1 1 1 2 0 1 3 1 0 4 0 0 3.6.2 Hexahedral Solid Element Simulations The applicability of the developed model to different element geometries was explored by running simulations using a hexahedral solid element as shown in Fig.3.18. The element was assumed to be transversely isotropic with the properties in the TD assigned to the TTD. The tetrahedral element is of size 5mm x 5mm x 5mm and the node numbers, constraints and coordinate system are also shown in Fig 3.18. Furthermore, Table 3.10 below summarizes the boundary conditions applied at the nodes and their degrees of freedom. For simulations run with the solid 68 element, the three-dimensional stiffness matrix discussed in section 3.1.4 was implemented. Simulation results produced using the solid element will be compared to that produced using the shell element for the same loading conditions to validate the model’s applicability for solid elements. Figure 3.18 Hexahedral solid element with boundary conditions. Table 3.10 Boundary conditions for solid element Node Number X - Direction Y - Direction Z - Direction 1 1 1 1 2 0 1 1 3 0 1 0 4 1 1 0 5 1 0 1 6 0 0 1 7 0 0 0 8 1 0 0 3.7 Orthotropic Linear Thermoviscoelastic Model Verification The orthotropic thermomechanical material model based on TTSP and TTSSM framework is implemented in LS-DYNA® as a user material model for shell elements. The implemented model is verified using the analytical solutions introduced in section 3.4. The model predictions will also be compared with the isothermal stress relaxation experimental results. All simulations, except when indicated are performed with the shell element shown in Fig.3.16. For non-isothermal simulations *MAT_ADD_THERMAL_EXPANSION in the LS- DYNA® keyword [20] is introduced to include the thermal expansion/shrinkage effect in the simulations. In the work by Yan et. al. [21], the CTE as a function of temperature over the range 69 of 30º-160ºC has been determined for Celgard®2400 and two other separators. Table 3.11 presents the CTE in different temperature ranges as functions of temperature used in the current work. Table 3.11 CTE for Celgard®2400 [21] -6 Material Temperature (ºC) CTE (x10 /ºC) 2400 MD 30-65 44 65-175 −3.45 × 10−13 𝑇 7.6262 2400 TD 30-130 1.4T+26 130-140 0 140-160 -88T+12044 The analytical solutions are computed using Eqn.3.33 and Eqn.3.34 for cases with constant temperature history. Analytical solutions for stress relaxation with two-step loadings at constant temperatures were calculated using Eqn.3.35. For cases with step temperature histories, Eqn.3.40 and Eqn.3.42 are used. The TTSP parameters in Table 3.6, the Prony series parameters in Table 3.7, and the TTSSM parameters in Table 3.8 were the input parameters for the model. Fig.3.19 presents the stress relaxations in the MD in air at 20ºC, 50ºC, 80ºC, and 100ºC. The simulation results are compared with analytical solutions and experimental results for the selected temperature levels. As shown, the simulations and analytical solutions coincided at all four temperature levels, indicating the model has been implemented correctly in the user material model. Hence, verifying the numerical model. On the other hand, the simulations do not coincide with the experimental curves exactly. For 20ºC, 50ºC, and 100ºC, the simulation agreed with the experimental curves reasonably well. The largest discrepancy occurred at 80ºC. The predicted curve had a different shape from the experimental curve. Fig.3.20 and Fig.3.21 present the stress relaxations in air at 20ºC, 50ºC, 80ºC, and 100ºC in the TD and the 45° off-axis direction, respectively. Fig.3.22-Fig.3.24 present the stress relaxations in DMC for the MD, TD, and 45° off- axis direction, respectively. From the results, analytical solutions coincide with the simulation results showing that the electrolyte effect has been introduced into the model correctly. These results show that the model has captured the overall trend but may not be as accurate at certain temperatures. Furthermore, the model predictions agreed with experimental results in the TD better than the other two material orientations. 70 (a) (b) Figure 3.19 Comparison of the stress history predicted by simulations with analytical solutions and experimental data for stress relaxation in the MD at 0.2% strain in air at (a) 20ºC and 50ºC (b) 80ºC and 100ºC. (a) (b) Figure 3.20 Comparison of the stress history predicted by simulations with analytical solutions and experimental data for stress relaxation in the TD at 0.3% strain in air at (a) 20ºC, 50ºC (b) 80ºC,100ºC. 71 (a) (b) Figure 3.21 Comparison of the stress history predicted by simulations with experimental data for stress relaxation in the 45° off-axis direction at 0.3% strain in air at (a) 20ºC and 50ºC (b) 80ºC and100ºC. Figure 3.22 Comparison of the stress history predicted by simulations with analytical solutions and experimental data in the MD for stress relaxation at 0.4% strain in DMC at 21ºC, 30ºC, and 40ºC. 72 Figure 3.23 Comparison of the stress history predicted by simulation with analytical solutions and experimental data in the TD for stress relaxation at 0.4% strain in DMC at 22ºC, 30ºC, and 40ºC. Figure 3.24 Comparison of the stress history predicted by simulation with experimental data for stress relaxation in the 45° off-axis direction at 0.4% strain in DMC at 20ºC, 30ºC, and 40ºC. The difference in mechanical response in the MD and TD can be related to the formation of the microstructure of the separator. The porous structure of dry-processed separators was generated by stretching along the MD [22]. This microstructure is not in a state of equilibrium. With increasing temperature, the microstructure tends to return to its equilibrium state and this is manifested mainly in the MD. For Celgard®2400, the CTE changes from expansion to shrinkage 73 at 65ºC along the MD but at 130ºC along the TD [21]. With increasing temperature, the modulus decreases more rapidly in the MD than that in the TD. As shown in Figs.3.5a and 3.5b, at lower temperatures, the relaxation modulus is significantly higher in the MD than in the TD. At temperatures higher than 80ºC, however, the trend starts to reverse. This is also shown in Fig.3.14a. The MD and TD master curves cross over at long relaxation times. According to the TTSP, a longer time is equivalent to a higher temperature. Furthermore, the developed model was verified for stress relaxation with step loading at constant temperature and for stress relaxation with step temperature histories at constant strain levels. The model was verified according to the step histories presented in Fig.3.25 and Fig. 3.26 in the MD and TD. (a) (b) Figure 3.25 Strain loading for two-step stress relaxation at constant temperatures for (a) 0.1%- 0.2% strain levels (low to high), (b) 0.2%-0.1% strain levels (high to low) cases. (a) (b) Figure 3.26 Applied temperature history (a) 30ºC – 40ºC and (b) 30ºC – 40ºC – 50ºC – 60ºC. 74 (a) (b) (c) (d) Figure 3.27 Model verification with analytical solution for two-step (low to high) and two-step (high to low) loading cases respectively in the (a,b) MD and (c,d) TD. 75 (a) (b) (c) Figure 3.28 Model verification with the analytical solution for stress relaxation with step temperature history with and without considering thermal expansion in (a) MD 30ºC – 40ºC, (b) TD 30C – 40C and (c) MD 30ºC – 40ºC – 50ºC – 60ºC. The results presented in Fig.3.27 and Fig.3.28 show that the analytical solutions coincide with the simulation results showing that the model has been implemented correctly. Finally, the applicability of the developed model to solid elements is verified by comparing simulations carried out using the shell element to those produced using the hexahedral solid element. Fig.3.29 compares single and multi-step stress relaxation simulations in the MD and TD. The results show that the predictions from simulations using both elements coincide for the different loading conditions and for the different material orientations. 76 (a) (b) (c) (d) Figure 3.29 Model verification by comparing shell and solid element simulations for stress relaxation cases in air for (a) MD 20ºC – 0.2% strain, (b) TD 20C – 0.3% strain (c) MD 30ºC – 40ºC – 50ºC – 60ºC – 0.2% strain and (d) TD 30ºC – 40ºC – 50ºC – 60ºC – 0.3% strain. 3.8 Orthotropic Linear Thermoviscoelastic Model Validation In this section, the developed model is validated by comparing simulations with experimental data carried out under various combined thermomechanical loadings separate from those used to determine the model parameters. The experimental data used to validate the model were generated from uniaxial tensile tests at constant temperatures, stress relaxation tests under non-isothermal conditions and iso-stress temperature ramp tests. Firstly, the model is validated against uniaxial tensile tests carried out on samples cut along the MD and TD at temperatures ranging from 30ºC to 60ºC. The tests were also carried out at displacement rates of 0.001mm/s and 0.01mm/s. The comparison between the experimental data 77 and simulation results is shown in Fig.3.30. From the results, the model predictions agree well with the experimental data as they capture the trends in the material anisotropy, rate and temperature dependence. Furthermore, the modulus of elasticity (Young’s modulus) was determined from the experimental and simulated stress-strain curves. The Young’s modulus was determined by measuring the slope of the stress-strain curves in the deformation region from 0.1% to 0.5% strain. The results are presented in Table 3.12 and they show that the predicted modulus is comparable to the experimental modulus for the same loading conditions. (a) (b) (c) (d) Figure 3.30 Comparison between simulations and experimental data for uniaxial tensile tests at different temperatures and displacement rates in (a) MD 30C, (b) MD 40C, (c) MD 60C (d)TD 30C (e) TD 40C (f) TD 50C (g) TD 60C. 78 Figure 3.30 (Cont’d) (e) (f) (g) Table 3.12 Experimental and predicted young’s modulus of Celgard® 2400 T Machine Direction (MD) Transverse Direction (TD) (ºC) Experiment Simulation Experiment Simulation (MPa) (MPa) (MPa) (MPa) 0.001 0.01 0.001 0.01 0.001 0.01 0.001 0.01 mm/s mm/s mm/s mm/s mm/s mm/s mm/s mm/s 30 1103.00 1306.3 1143.6 1304.9 609.46 655.19 627.75 655.73 40 808.81 1091.1 809.68 1089.8 525.57 580.18 542.87 577.27 50 - - - - 414.06 467.05 454.09 513.26 60 373.63 495.93 439.34 615.87 340.91 383.19 371.16 440.75 Furthermore, the model is validated by comparing simulations with experimental data under non-isothermal conditions. Stress relaxation experiments were performed in the MD in air 79 under non-isothermal conditions following a multi-step temperature history, as shown in Fig.3.31. The experiments were performed for three cases: two with two-step 30ºC-40ºC and 60ºC-70ºC, and one with four-step 30ºC-40ºC-50ºC-60ºC. The measurements were repeated at least twice for consistency and the averaged results are presented here. The recorded temperature histories are also shown in Fig.3.31. (a) (b) (c) (d) Figure 3.31 Comparison of the stress history predicted by LS-DYNA® simulation using the orthotropic thermo-mechanical material model with experimental data in the MD for non- isothermal stress relaxation at temperature steps of (a) 30ºC-40ºC, (b) 60ºC-70ºC, and (c) 30ºC- 40ºC-50ºC-60ºC. (d) Thermal expansion/shrinkage of Celgard®2400 MD in the temperature range of 30º-80ºC [21]. The simulations were carried out for stress relaxation with the recorded temperature history. To examine the effect of CTE, simulations were performed with and without considering 80 CTE. The two-step results showed that the CTE had a significant effect at a 30ºC-40ºC increment but little effect at a 60ºC-70ºC increment. To understand this, the CTE data of Celgard ®2400 MD was inspected. The data in the range of 30º-80ºC is presented in Fig.3.31d. Celgard®2400 MD sample displays a nearly linear thermal expansion up to 65ºC and then starts to shrink with increasing temperature. As a result, the thermal strains at 60ºC and 70ºC are almost the same and hence the effect of CTE is negligible for the 60ºC-70ºC step. For the four-step case, the simulation with CTE agrees with experimental results quite well. Without the consideration of CTE, the simulation overestimated the stress value by over 100%. These cases demonstrate the importance of thermomechanical simulations. Finally, the model is validated against iso-stress temperature ramp experimental results obtained in [21]. As shown in Fig.3.32, the experiments were performed in the MD at three stress levels: 0.067, 0.13, and 3MPa, and in the TD at 0.067, and 0.13 MPa over the temperature range of 25ºC -110ºC. Simulations were performed with and without the effect of thermal expansion. The experimental results showed that the thermal expansion/shrinkage effect exceeded that of creep at lower stress levels (0.067 and 0.13MPa) in both the MD and TD. Without considering the CTE, the predicted strains were the accumulated creep strains only, which were far from experimental curves. Simulations with the CTE effect correctly predicted the trend of strain evolution over this broad temperature range for each case. From the experimental results for the MD under iso-stress of 3MPa, the creep strain accumulation exceeded the shrinkage effect and the strain increased with temperature. The simulation with the CTE effect captured this trend. The predicted curve was closer to the experimental curve than the simulation without the CTE. 81 (a) (b) (c) (d) Figure 3.32 Comparison between simulations and experimental data for iso-stress temperature ramp at stress levels of (a) MD - 0.067MPa, (b) MD – 0.13MPa and 3MPa, (c) TD – 0.06MPa (d)TD – 0.13MPa. 3.9 Summary In this chapter, the development of an orthotropic thermo-mechanical model has been presented. The model is built upon a linear viscoelastic framework with temperature dependence based on the time-temperature superposition principle (TTSP). To consider the effect of electrolyte solvents, a time-temperature-solvent superposition principle (TTSSM) has been proposed. Orthotropic viscoelastic characterizations have been performed for a PP separator, Celgard ® 2400 in air and in DMC. The model parameters for Celgard® 2400 have been established. 82 The developed model was implemented as a user material model in LS-DYNA®. Using one element FE model, the user material model was tested in simulations for a wide range of experimental cases, including uniaxial stress relaxation under isothermal and non-isothermal conditions in the MD, TD, and 45º off-axis direction, uniaxial tensile tests, and iso-stress experiments with a temperature ramp. Analytical solutions were developed and implemented to verify the model. The model was validated against experimental data and although the predictions do not coincide exactly with the experimental results for each particular case, they predicted the overall trend in material anisotropy, temperature dependence, and solvent effect. The results also demonstrated that the simulations without the thermal expansion/shrinkage behavior of the separator resulted in large errors. This current model is limited to predictions of the separator response under small strains (linear viscoelastic region), the extension of the model to capture and predict the separator response under large deformations will be addressed in the subsequent chapters. 83 REFERENCES 1. R. C. Ihuaenyi, S. Yan, J. Deng, C. Bae, I. Sakib, and X. Xiao, Orthotropic Thermo- Viscoelastic Model for Polymeric Battery Separators with Electrolyte Effect. Journal of The Electrochemical Society, 168(9), p.090536, 2021. 2. J. D. Ferry, Viscoelastic Properties of Polymers, John Wiley & Sons, 1980. 3. I. M. Ward, J. Sweeney, An Introduction to the Mechanical Properties of Solid Polymers, John Wiley & Sons, 2004. 4. A. S. Wineman, K. R. Rajagopal, Mechanical response of polymers: an introduction. Cambridge university press, 2000. 5. H. F. Brinson, L. C. Brinson, Polymer engineering science and viscoelasticity. An introduction 2008. 6. R. M. Christensen. Theory of Viscoelasticity: An Introduction. Academic Press, New York, second edition, 1982. 7. M.A. Puso and J.A. Weiss, Finite element implementation of anisotropic quasilinear viscoelasticity using a discrete spectrum approximation. Journal of Biomechanical Engineering, 120:62-70, 1998. 8. S. Yan, J. Deng, C. Bae, S. Kalnaus, and X. Xiao, Orthotropic viscoelastic modeling of polymeric battery separator. Journal of The Electrochemical Society, 167(9), p.090530, 2020. 9. I. M. Daniel, O. Ishai, Engineering Mechanics of Composite Materials, Oxford University Press, 2006. 10. S. Yan, J. Deng, C. Bae, Y. He, A. Asta, and X. Xiao, In-plane orthotropic property characterization of a polymeric battery separator. Polymer Testing, 72, pp.46-54, 2018. 11. A. Sheidaei, X. Xiao, X. Huang, and J. Hitt, Mechanical behavior of a battery separator in electrolyte solutions. Journal of Power Sources, 196(20), pp.8728-8734, 2011. 12. I.V. Yannas, Transition from linear to nonlinear viscoelastic behavior. III. Linearity below and above Tg. Rubber Chemistry and Technology, 45(6), pp.1623-1629, 1972. 13. J. Xu, L. Wang, J.Guan, and S. Yin, Coupled effect of strain rate and solvent on dynamic mechanical behaviors of separators in lithium-ion batteries. Materials & Design, 95, pp.319- 328, 2016. 14. S. Yan, X. Xiao, X. Huang, X. Li, and Y. Qi, Unveiling the environment-dependent mechanical properties of porous polypropylene separators. Polymer, 55(24), pp.6282-6292, 2014. 15. O. Starkova and A. Aniskevich, Limits of linear viscoelastic behavior of polymers. Mechanics of Time-Dependent Materials, 11(2), pp.111-126, 2007. 16. M. L. Williams, R. F. Landel, J.D. Ferry, Journal of the American Chemical society, 77(14), 3701, 1955. 17. Y. Miyano, M. Nakada, Durability of Fiber-Reinforced Polymers. Wiley-VCH Verlag, 2018. 18. J. Chen, H. Hu, S. Li, and Y. He. Evolution of mechanical properties of polypropylene separator in liquid electrolytes for lithium‐ion batteries. Journal of Applied Polymer Science, 135(27), p.46441, 2018. 84 19. A. E. Krauklis, A. G. Akulichev, A. I. Gagani, and A. T. Echtermeyer. Time–temperature– plasticization superposition principle: predicting creep of a plasticized epoxy. Polymers, 11(11), p.1848, 2019. 20. LS-DYNA Keyword Manual, Livermore Software Technology Corporation (LSTC), 2014. 21. S. Yan, J. Deng, C. Bae, and X. Xiao. Thermal expansion/shrinkage measurement of battery separators using a dynamic mechanical analyzer. Polymer Testing, 71, pp.65-71, 2018. 22. P. Arora, and Z. Zhang. Battery separators. Chemical reviews, 104(10), pp.4419-4462, 2004. 85 Chapter 4 Orthotropic Nonlinear Thermoviscoelastic Constitutive Modeling In this chapter, the development of an orthotropic nonlinear thermoviscoelastic material model for polymeric separators that accounts for the material anisotropy and rate dependence is presented. In the previous chapter (Ch.3), we developed, verified, and validated an orthotropic linear viscoelastic material model for polymeric separators. The results showed that the model was capable of accurately predicting the thermomechanical response of the separator under loadings within the linear viscoelastic limit. This present study extends the model to predict the thermomechanical response of the material in the nonlinear deformation region up to the onset of yielding or evolution of irrecoverable deformations. This model is developed based on the Schapery nonlinear viscoelastic model and the temperature dependence is introduced through the time-temperature superposition principle (TTSP). The model is implemented in LS-DYNA® finite element (FE) package as a user-defined subroutine. The model parameters are determined for a polypropylene (PP) separator. Analytical solutions for step loading and step temperature history conditions are introduced and compared to the model predictions to verify the model implementation. The predicted material responses under large deformations in isothermal and non- isothermal temperature conditions for stress relaxation, creep and tensile loadings at different rates are compared with the experimental data to validate the model predictions. The work presented in this chapter has been summarized and published [1]. 4.1 Orthotropic Nonlinear Viscoelastic Model Development 4.1.1 Constitutive Model Overview Nonlinearities arise in the stress-strain response of viscoelastic materials when the applied loading is larger than the linear viscoelastic strain or stress limit. Within the linear viscoelastic limit, the relaxation and creep response of a viscoelastic material are functions of time alone. However, above this limit, creep and relaxation functions become functions of both time and stress or strain. The viscoelastic response of the material becomes nonlinear above the linear viscoelastic limit. This nonlinearity in the material response evolves before the onset of irrecoverable deformations or yielding. Fig.4.1 below shows the stress relaxation response of the TD of Celgard®2400 PP separator carried out at different strain levels. From Fig.4.1 it is clear that above 0.3% the stress relaxation response of the material becomes dependent on the strain level. Hence, 86 when developing models for viscoelastic materials under large deformations, it is vital to consider the nonlinearity in the stress-strain response of the material before yielding. Figure 4.1 Stress relaxation response of a PP separator in the TD at 20-22°C. The nonlinearities in the stress-strain response of the separator material arising due to large deformations before yielding are accounted for through a modified version of the stiffness-based formulation of the Schapery single integral viscoelastic model. For uniaxial loading, the Schapery single integral model in its stiffness form is written as: 𝑡 𝑑ℎ2 (𝜀)𝜀(𝜏) 𝜎(𝑡) = ℎ∞ (𝜀) ∙ 𝐺∞ ∙ 𝜀(𝑡) + ℎ1 (𝜀) ∙ ∫ ∆𝐺[𝜌(𝑡) − 𝜌(𝜏)] 𝑑𝜏 (4.1) 𝑑𝜏 0 where 𝜎 is the stress, 𝜀 is the strain, 𝐺∞ and ∆𝐺(𝑡)are the equilibrium and linear viscoelastic transient relaxation modulus. ℎ∞ , ℎ1 , ℎ2 , are strain-dependent nonlinearity parameters, and 𝜌 is the reduced time defined as: 𝑡 𝑑𝑡 ′ 𝜌(𝑡) = ∫ (4.2) 0 𝑎𝜀 [𝜀(𝑡 ′ )] ∙ 𝑎 𝑇 [𝑇, 𝑇0 ] 𝜏 𝑑𝑡 ′ 𝜌(𝜏) = ∫ ′ (4.3) 0 𝑎𝜀 [𝜀(𝑡 )] ∙ 𝑎 𝑇 [𝑇, 𝑇0 ] In Eqns.4.2 and 4.3, 𝑎𝜺 is the time-strain shifting factor and 𝑎 𝑇 is the time-temperature shifting factor, where T and T0 are current and reference temperatures respectively. According to TTSP, the effect of temperature is equivalent to that of time. When T > T0 , the time needed to produce the equivalent viscoelastic response will be reduced. In Eqn.4.1, the effects of temperature and strain are not coupled. The other strain-dependent nonlinearity parameters ( ℎ∞ , ℎ1 , and ℎ2 ) 87 are considered temperature-independent. In this way, the strain-dependent and temperature- dependent parameters are decoupled in the current model. When the input strain is below the linear viscoelastic strain limit, the values of the nonlinear parameters ℎ∞ , ℎ1 , ℎ2 , and 𝑎𝜺 become equal to 1, and Eqn.4.1 reduces to the hereditary integral representation for linear viscoelastic materials. 4.1.2 Numerical Integration of the Nonlinear Viscoelastic Hereditary Integral with a Kernel of Prony Series To evaluate the nonlinear viscoelastic hereditary integral for loading with arbitrary stress, strain, and temperature histories, numerical solutions have to be developed. The discretization algorithm implemented in this section is an extension of the approach introduced by Yan et al [2] based on the work by Puso and Weiss [3] for the linear viscoelastic hereditary integral. In the present work, this algorithm is extended for nonlinear viscoelastic modeling with the Schapery hereditary integral. To evaluate the Schapery hereditary integral, the generalized Maxwell model is assumed to be composed of nonlinear springs and dashpots. From the generalized Maxwell model, if the stress in each component is known, the total stress is given as the sum of the stresses in the individual components in the form: 𝑛 𝜎(𝑡) = 𝜎∞ + ∑ 𝜎 𝑖 (𝑡) (4.4) 𝑖=1 where 𝜎∞ is the equilibrium stress and from Eqn.4.1, it can be represented as: 𝜎∞ (𝑡) = ℎ∞ ∙ 𝐺∞ ∙ 𝜀(𝑡) (4.5) At a time 𝑡 + ∆𝑡 , the equilibrium stress increment becomes: ∆𝜎∞ = 𝜎∞ (𝑡 + ∆𝑡) − 𝜎∞ (𝑡) = ℎ∞ ∙ 𝐺∞ ∙ [𝜀(𝑡 + ∆𝑡) − 𝜀(𝑡)] = ℎ∞ ∙ 𝐺∞ ∙ ∆𝜀 (4.6) Furthermore, the stress in each individual Maxwell component is given by: 𝑡 𝑑ℎ2 𝜀(𝜏) 𝜎 𝑖 (𝑡) = ℎ1 ∙ ∫ ∆𝐺[𝜌(𝑡) − 𝜌(𝜏)] 𝑑𝜏 (4.7) 𝑑𝜏 0 At a time 𝑡 + ∆𝑡 , the stress increment is given by: 𝑡+∆𝑡 𝑑ℎ2 𝜀(𝜏) 𝜎 𝑖 (𝑡 + ∆𝑡) = ℎ1 ∙ ∫ ∆𝐺[𝜌(𝑡 + ∆𝑡) − 𝜌(𝜏)] 𝑑𝜏 (4.8) 𝑑𝜏 0 88 Eqn. 4.8 is evaluated by separating it into two parts that are evaluated separately: 𝑡 𝑖 (𝑡 𝑑ℎ2 𝜀(𝜏) 𝜎 + ∆𝑡) = ℎ1 ∙ ∫ ∆𝐺[𝜌(𝑡 + ∆𝑡) − 𝜌(𝜏)] 𝑑𝜏 𝑑𝜏 0 𝑡+∆𝑡 𝑑ℎ2 𝜀(𝜏) + ℎ1 ∫ ∆𝐺[𝜌(𝑡 + ∆𝑡) − 𝜌(𝜏)] 𝑑𝜏 (4.9) 𝑑𝜏 𝑡 Prior to the evaluation of both parts of Eqn.4.8, we should note that the reduced time increment is defined as: Δ𝜌 = 𝜌(𝑡 + ∆𝑡) − 𝜌(𝑡) (4.10) The second term in Eqn.4.8 can be solved approximately by using the mean value theorem when ∆𝑡 is small and hence: 𝑡+∆𝑡 𝑑ℎ2 𝜀(𝜏) ℎ1 ∙ ∫ ∆𝐺[𝜌(𝑡 + ∆𝑡) − 𝜌(𝜏)] 𝑑𝜏 𝑑𝜏 𝑡 𝑡+∆𝑡 𝜀(𝑡 + ∆𝑡) − 𝜀(𝑡) = ℎ1 ∙ ℎ2 { } ∫ ∆𝐺[𝜌(𝑡 + ∆𝑡) − 𝜌(𝜏)]𝑑𝜏 ∆𝑡 𝑡 ∆𝜀 = ℎ1 ∙ ℎ2 { } ∆𝐺(∆𝜌)∆t = ℎ1 ∙ ℎ2 ∙ ∆𝜀 ∙ ∆𝐺(∆𝜌) (4.11) ∆𝑡 To evaluate the first term in Eqn.4.8, a kernel function is defined by setting ∆𝐺 to be a multiple-term exponential function which is the Prony series expansion of the generalized Maxwell model for the transient relaxation modulus represented as: 𝑛 𝑡 ∆𝐺(t) = ∑ 𝐺𝑖 ∙ 𝑒𝑥𝑝 [− ] (4.12) 𝜏𝑖 𝑖=1 This shows that the transient portion of the stress relaxation modulus of the material is contributed by n individual Maxwell components, each with its characteristic relaxation time 𝜏𝑖 (𝑖 = 1,2, … 𝑛). Hence the first term of Eqn.4.8 becomes: 𝑡 𝑑ℎ2 𝜀(𝜏) ℎ1 ∙ ∫ ∆𝐺[𝜌(𝑡 + ∆𝑡) − 𝜌(𝜏)] 𝑑𝜏 𝑑𝜏 0 𝑛 𝑡 𝜌(𝑡 + ∆𝑡) − 𝜌(𝜏) 𝑑ℎ2 𝜀(𝜏) = ℎ1 ∙ ∑ ∫ 𝐺𝑖 ∙ 𝑒𝑥𝑝 (− ) 𝑑𝜏 𝜏𝑖 𝑑𝜏 𝑖=1 0 89 𝑛 𝑡 ∆𝜌 𝜌(𝑡) − 𝜌(𝜏) 𝑑ℎ2 𝜀(𝜏) = ∑ 𝑒𝑥𝑝 [− ] ∙ ℎ1 ∙ ∫ 𝐺𝑖 ∙ 𝑒𝑥𝑝 (− ) 𝑑𝜏 𝜏𝑖 𝜏𝑖 𝑑𝜏 𝑖=1 0 𝑛 ∆𝜌 = ∑ 𝑒𝑥𝑝 [− ] ∙ 𝜎 𝑖 (𝑡) (4.13) 𝜏𝑖 𝑖=1 At a time 𝑡 + ∆𝑡 , the stress increment for each i-th Maxwell component is given by: ∆𝜎 𝑖 = 𝜎 𝑖 (𝑡 + ∆𝑡) − 𝜎 𝑖 (𝑡) (4.14) ∆𝜌 ∆𝜌 ∆𝜎 𝑖 = ℎ1 ∙ ℎ2 ∙ ∆𝜀 ∙ 𝐺𝑖 ∙ 𝑒𝑥𝑝 [− ] + 𝑒𝑥𝑝 [− ] ∙ 𝜎 𝑖 (𝑡) − 𝜎 𝑖 (𝑡) (4.15) 𝜏𝑖 𝜏𝑖 From Eqn. 4.4, we know that the total stress increment is the sum of the equilibrium stress increment and the stress increment for each i-th Maxwell component. The total stress increment can now be expressed as: 𝑛 ∆𝜌 ∆𝜌 ∆𝜎 = ℎ∞ ∙ 𝐺∞ ∙ ∆𝜀 + ∑ {ℎ1 ∙ ℎ2 ∙ ∆𝜀 ∙ 𝐺𝑖 ∙ 𝑒𝑥𝑝 [− ] + 𝑒𝑥𝑝 [− ] ∙ 𝜎 𝑖 (𝑡) − 𝜎 𝑖 (𝑡)} (4.16) 𝜏𝑖 𝜏𝑖 𝑖=1 For the stiffness-based formulation of the Schapery nonlinear model, the nonlinearity parameters ℎ1 𝑎𝑛𝑑 ℎ2 are often determined as a product (ℎ1 ∙ ℎ2) using single-step relaxation tests. Separating the two parameters requires a two-step stress relaxation test [4]. Furthermore, the formulation in Eqn. 4.16 takes away the need to separate the parameters ℎ1 𝑎𝑛𝑑 ℎ2. The product (ℎ1 ∙ ℎ2 ) can be represented as a single term ℎ𝑘 and hence Eqn. 4.16 becomes: 𝑛 ∆𝜌 ∆𝜌 ∆𝜎 = ℎ∞ ∙ 𝐺∞ ∙ ∆𝜀 + ∑ {ℎ𝑘 ∙ ∆𝜀 ∙ 𝐺𝑖 ∙ 𝑒𝑥𝑝 [− ] + 𝑒𝑥𝑝 [− ] ∙ 𝜎 𝑖 (𝑡) − 𝜎 𝑖 (𝑡)} (4.17) 𝜏𝑖 𝜏𝑖 𝑖=1 Figure 4.2 below presents the flow chart for the implementation of this algorithm in LS- DYNA® as a user-defined subroutine. 90 Figure 4.2 Flowchart for the implementation of evaluation of the nonlinear hereditary integral with Prony series. 4.1.3 Modifications to the Schapery Nonlinear Viscoelastic Model For the numerical implementation of the evaluation of the hereditary integral (Eqn.4.17), to allow for small time steps and stability in explicit simulations, the strain-time shift factor (𝑎𝜺 ) is set to equal to 1. Furthermore, the individual contributions of ℎ1 and ℎ2 may be difficult to separate. Also, as shown above, the final result of the implementation of the discretization algorithm (Eqn.4.17) takes away the necessity to separate the nonlinearity parameters ℎ1 and ℎ2. Furthermore, in literature, the nonlinear viscoelastic response of polymeric materials has been modeled with the assumption that some but not all of the Schapery nonlinear parameters are equal to 1 [5-8]. Hence, the nonlinearity in this model is considered through only two strain-dependent nonlinearity parameters ℎ∞ and ℎ𝑘 , with ℎ𝑘 = ℎ1 ∙ ℎ2. 4.1.4 Plane Stress Stiffness Matrix Using the elastic-viscoelastic correspondence principle, the stiffness matrix for the plane- stress condition of a viscoelastic material is given by: 𝐺11 (𝑡) 𝑣12 𝐺22(𝑡) 0 1 𝑮= ∙ [𝑣21 𝐺11(𝑡) 𝐺22 (𝑡) 0 ] (4.18) [1 − 𝑣12 𝑣21 ] 0 0 [1 − 𝑣12 𝑣21 ]𝐺66 (𝑡) 91 where 𝐺11 (𝑡), 𝐺22 (𝑡), and 𝐺66 (𝑡) are the relaxation moduli in the MD, TD and in-plane shear respectively, 𝑣12 and 𝑣21 are the major and minor Poison’s ratios. For orthotropic viscoelastic materials, 𝑣12 and 𝑣21 are not independent but related by Betti’s reciprocal law [9] such that: 𝑣12 𝐺22 = 𝑣21 𝐺11 (4.19) In nonlinear viscoelastic models for PE thin films [10,11], the above relation was used resulting in symmetric in-plane material compliance matrices. Furthermore, the bi-axial strain fields have been measured using digital image correlation (DIC) in tensile experiments for Celgard®2400 in the strain range of 0.1% - 5% in air by Yan et al [12]. It was found that 𝑣12 had a constant value of 0.17 in the entire strain range whereas the value of 𝑣21 increased with 𝜀1 , the strain in the MD direction. 𝑣21 was fitted with a linear relationship expressed as: 𝑣21 = 0.0354 ∙ 𝜀1 + 0.0801 (4.20) where the unit of 𝜀1 is in percentage strain. When 𝜀1 is greater than 1%, the value of 𝑣21 determined using Eqn.4.19 starts to deviate from that calculated using Eqn.4.20 [12]. Therefore, Eqn.4.20 was used to calculate 𝑣21 in the current model. In other words, the model has five independent parameters in the plane stress stiffness matrix. This behavior is suspected to be caused by the unique anisotropic porous microstructure of the material. 4.2 Orthotropic Nonlinear Thermoviscoelastic Material Characterization 4.2.1 Experimental Procedure The polymeric separator used in this study is Celgard®2400, a porous polypropylene (PP) material. Celgard®2400 is a single-layer film with a thickness of 25µm. For the experimental procedures, the specimens were cut along the MD, TD and 45° off-axis direction using a razor blade as described in a previous work by Yan et al [12]. The nominal dimensions of the specimen were 45 mm in length and 5 mm in width. The experiments were carried out using the RSA-G2 dynamical mechanical analyzer (DMA) as shown in Fig. 4.3. For measurements below ambient temperature conditions (20°C), liquified nitrogen (LN2) was used to provide the needed cooling. For tests carried out under isothermal conditions, the sample was kept at the specified temperature for 4 minutes to allow for temperature equilibration before loading. 92 Figure 4.3 RSA-G2 dynamical mechanical analyzer testing set up with a mounted sample. Stress relaxation, creep, iso-stress, and tensile experiments were performed at various loading and temperature conditions. Table 4.1 presents a summary of the experiments and the purpose of each set of experiments. Table 4.1 Experimental summary Experiment Specimen Test Conditions Purpose MD, TD, 0.2-5% strain levels at 20°C Determination of Stress relaxation 45° off-axis nonlinearity parameters AL thin foil 0.02MPa at 3°C/min ramp rate Instrument expansion Iso-stress calibration Creep MD, TD 3-10MPa, 20°C MD, TD 0.7-3% strain, 30°C, 40°C, Stress relaxation Validation under iso- 60°C thermal condition MD, TD 0.001, 0.01mm/s at 30°C, 40°C Tensile 50°C MD, TD 30°C to 60°C with T=10°C Validation under non Stress relaxation iso-thermal condition To characterize the nonlinear viscoelastic response of the separator, stress relaxation tests were carried out at 20°C in the MD, TD and 45° off-axis direction. The tests were performed at constant strain levels ranging from 0.2% to 5%. To provide data for model validation under creep loading, creep tests were carried out at 20°C for the MD and TD directions. The linear viscoelastic stress limits of Celgard ®2400 have 93 been determined at ambient temperature conditions as 5MPa for the MD and 2MPa for TD [12]. For nonlinear model validation, creep tests were performed at 4, 5, 8 and 10MPa, for MD, and 3 and 5MPa for TD. The creep tests were run for 20 mins. For validation at elevated temperatures up to 60°C, stress relaxation tests were carried out at several strain levels in the nonlinear region. For MD, stress relaxations were performed at 1% strain at 30°C and 40°C, and at 0.7% strain at 60°C. For TD, 1% and 2% strain at 30°C, and 40°C. To validate the thermomechanical model for non-isothermal conditions, stress relaxation experiments were performed with a multi-step temperature history of 30°C-40°C-50°C-60°C. At the first step of the test, each sample was kept isothermal for 4 minutes at the desired temperature before loading and data collection. At the subsequent temperature steps, data were collected simultaneously with the temperature ramp. The stress relaxation tests were run for 20 mins. In the non-isothermal stress relaxation procedure, the total applied strain is contributed to, by the mechanical strain that induces relaxation and the thermal strain. Furthermore, in DMA systems, the strain is determined by measuring the grip displacement with a sensor located on the drive shaft. The thermal expansion of the grip and the drive shaft also contribute to the measured thermal strain, in addition to that of the sample [13]. The contribution of the instrument to thermal strain can be determined by calibration, using a sample with known CTE. To accurately determine the contribution of the instrument to the thermal strain, iso-stress temperature ramp tests were performed using an Al thin foil with a thickness of about 25 µm and an initial gage length of 15mm. During the test, a constant small stress level of 0.02MPa was applied to keep the sample straight and the temperature was ramped at a rate of 3℃/𝑚𝑖𝑛 from ambient to 150℃. To validate the model predictions for the stress-strain response, uniaxial tensile tests were run at displacement rates of 0.001mm/s and 0.01mm/s. The tests were carried out at constant temperatures of 30°C, 40°C, and 50°C in the TD, and at 30°C and 40°C in the MD. For every test, each sample was mounted and aligned as shown in Fig. 4.3. A preload force of ~3g was applied to ensure the sample was straight before starting each test. All tests were repeated two to three times for consistency and the average results are presented. 4.2.2 Stress Relaxation at Different Strain Levels Stress relaxation tests were carried out on Celgard ®2400 samples cut along the MD, TD and 45° off-axis direction at different strain levels ranging from 0.2% - 5% as described in the 94 experimental procedure. These tests were carried out at the reference temperature of the model calibration (20°C) to determine the strain-dependent model parameters. The results are presented in Fig.4.4 and they show that the stress relaxation curves follow the trend described by the viscoelastic theory. With an increase in the applied constant strain, the stress relaxation modulus decreases and vice versa. (a) (b) (c) Figure 4.4 Stress relaxation results at different strain levels and at 20°C for tests carried out on samples cut along the (a) MD, (b) TD, and (c) 45° off-axis direction. 4.2.3 Determination of In-Plane Shear Relaxation Modulus To implement the in-plane stiffness matrix (Eqn.4.18) in the thermoviscoelastic model, the in-plane shear modulus needs to be determined. However, it is difficult to experimentally measure the in-plane shear property of thin films. In Chapter 3, the linear viscoelastic in-plane shear relaxation modulus was determined from the modulus measured from samples cut along the MD, 95 TD, and 45° off-axis through a transformation equation based on the elastic-viscoelastic correspondence principle [14,15]. In literature, the elastic-viscoelastic correspondence principle has been applied to analyze the response of viscoelastic materials under small deformations [2,16,17-20] and by Zhang et al [21] to consider the viscoelastic recovery response in the nonlinear deformation region. To determine the in-plane shear relaxation modulus in the nonlinear viscoelastic region, the transformation equation [16] is modified through the introduction of strain-dependent stiffness functions in the form: −1 4 (1 − 𝑣12 ) (1 − 𝑣21 ) 𝐺66 (𝑡, 𝜀) = [ − − ] (4.21) 𝐺45° (𝑡, 𝜀) 𝐺11 (𝑡, 𝜀) 𝐺22 (𝑡, 𝜀) where 𝐺11 (𝑡, 𝜀), 𝐺22 (𝑡, 𝜀), and 𝐺45° (𝑡, 𝜀) are the experimentally measured stress relaxation moduli of the samples cut along the MD, TD, and 45° off-axis direction respectively. Eqn.4.21 above was utilized to analytically determine the in-plane shear relaxation modulus for different constant strain levels at 20°C for the investigated material and the results are presented in Fig.4.5 below. The stress relaxation curves produced also followed the trend marked by the negative effect of the increasing constant strain level on the relaxation modulus. This is consistent with the viscoelastic theory. Figure 4.5 In-plane shear relaxation modulus for different constant strain levels at 20°C. 96 4.3 Parameter Identification 4.3.1 Overview This section presents the methodologies involved in the determination of the parameters required to calibrate the developed model. The parameters for the Schapery nonlinear viscoelastic model have been calibrated in literature [22-26], and these methods used present a basis for the methodology used in this work. The challenge has always been the determination of both the strain and temperature-dependent parameters. In literature, these parameters have been assumed to be decoupled [10,22,24-27] and coupled [28]. For decoupled formulations, the Schapery nonlinear parameters have been assumed to be both strain and temperature dependent in additive [24] and multiplicative forms [25,26]. The additive form in terms of stiffness is expressed as: ℎ𝑖 (𝜀, 𝑇) = ℎ𝑖 (𝜀) + ℎ𝑖 (𝑇) 𝑓𝑜𝑟 𝑖 = ∞, 1, 𝑜𝑟 2 (4.22) 𝑎(𝜀, 𝑇) = 𝑎(𝜀) + 𝑎(𝑇) (4.23) The multiplicative form is expressed as: ℎ𝑖 (𝜀, 𝑇) = ℎ𝑖 (𝜀) ∙ ℎ𝑖 (𝑇) 𝑓𝑜𝑟 𝑖 = ∞, 1, 𝑜𝑟 2 (4.24) 𝑎(𝜀, 𝑇) = 𝑎(𝜀) ∙ 𝑎(𝑇) (4.25) In this work, the strain and temperature-dependent parameters are assumed to be decoupled. The nonlinearity parameters carry the strain dependence (ℎ∞ (𝜀), ℎ𝑘 (𝜀)) and the temperature dependence is introduced through the TTSP using shifting factors. The temperature dependence is introduced through the transient linear viscoelastic modulus. Hence, the temperature-dependent parameters are determined in the linear viscoelastic strain limit for temperatures ranging from the reference temperature (20°C) to 60°C. Finally, the strain-dependent parameters are determined at the reference temperature for strain values ranging from the linear viscoelastic strain limit (0.2%-0.3%) to 5%. 4.3.2 Determination of Temperature-Dependent Model Parameters The temperature dependence is modeled through master curves based on TTSP [14,15,29- 31]. To build a master curve, viscoelastic properties such as the relaxation moduli or creep compliance measured at different temperatures are plotted vs time in a log-log plot. Figure 4.6 presents the relaxation moduli for MD, TD and in-plane shear directions measured within the linear viscoelastic strain limit and at constant temperatures varying from 20°C-60°C [16]. The master 97 curve is constructed by shifting individual curves horizontally and vertically to form a composite curve. (a) (b) (c) Figure 4.6 Stress relaxation curves produced from tests carried out under small strains for samples cut along the (a) MD, (b) TD, and (c) In-plane shear [16]. In Chapter 3, the master curves were built for calibrating the orthotropic linear viscoelastic model with relaxation moduli measured from 20°C to 110°C and by using horizontal shifts only as they were sufficient at the time for predictions in the linear regions. As a result, the master curves covered a wide period from 10-1 to 1013 seconds. To accurately fit the curves was challenging. With limited terms, some local regions were not accurately described. Furthermore, eliminating vertical shifts resulted in errors in some regions. For a semi-crystalline polymer like Celgard®2400, it has been suggested that the horizontal shift accounts for the change in the temperature-dependent response in the amorphous region and the vertical shift accounts for the change in the degree of crystallinity [4]. Hence, in the current work, the master curves were 98 reconstructed using the relaxation moduli curves measured from 20°C-60°C as shown in Fig.4.6. Both horizontal shifts in the log time scale and small vertical shifts in the log modulus scale were used to build the master curves. It should be mentioned that both the horizontal and vertical shifts were determined according to the need to construct smooth master curves rather than being linked to specific physical parameters such as the change in the degree of crystallinity. Fig.4.7 presents the master curves constructed from the data in Fig.4.6. (a) (b) (c) Figure 4.7 Time-temperature superposition master curves in the (a) MD, (b) TD, and (c) In- plane shear. In terms of the horizontal and vertical shifting factors, the stress relaxation modulus can be expressed as: 1 𝑡 𝐺(𝑡, 𝑇) = ∙ 𝐺 [ , 𝑇0 ] (4.26) 𝑏𝑇 𝑎𝑇 99 The viscoelastic behavior of a material over a wide range of time and temperatures can be described by the master curves and the shifting factors. In literature, the temperature dependency of the horizontal shifting factor has been modeled with the William-Landel-Ferry (WLF) model and the Arrhenius equation [14, 29-31]. In this chapter, the horizontal time-temperature shift factor (aT ) is modeled using the WLF model since the considered temperature range is less than 𝑇𝑔 + 100°𝐶 or less than 𝑇0 + 50°𝐶. The reported 𝑇𝑔 of Celgard®2400 is -15°C in air [32] and the reference temperature (𝑇0 ) is 20°C. The WLF model is expressed mathematically as: −𝐶1 ∙ (𝑇 − 𝑇0 ) 𝐿𝑜𝑔 (𝑎𝑇 ) = (4.27) 𝐶2 + (𝑇 − 𝑇0 ) where C1 and C2 are material parameters. The horizontal shifting factors for the master curves were fitted with the WLF model using the least square regression. The coefficients of determination R2 of these fittings were in the range of 0.995-0.9998, indicating excellent fitting (Fig.4.8). Table 4.2 compares the experimental horizontal shifting factors and the horizontal shifting factors calculated by the WLF model. (a) (b) Figure 4.8 Horizontal shifting factor fitting using the WLF model for master curves in the (a) MD, (b) TD, and (c) In-plane shear. 100 Figure 4.8 (Cont’d) (c) Table 4.2 Predicted horizontal shifting factors (𝐿𝑜𝑔(𝑎 𝑇 )) using the WLF model compared with experimental horizontal shifting factors Machine Direction Transverse Direction In-Plane Shear (MD) (TD) T Experiment Predicted Experiment Predicted Experiment Predicted (ºC ) 20 0.00 0.00 0.00 0.00 0.00 0.00 30 -0.60 -0.60 -0.60 -0.60 -0.80 -0.81 40 -1.40 -1.38 -1.30 -1.27 -1.90 -1.75 50 -2.40 -2.42 -2.00 -2.02 -2.65 -2.86 60 -3.80 -3.88 -2.80 -2.86 -4.10 -4.19 Furthermore, it has been argued that the vertical shift factor can be expressed empirically [33]. Hence, the vertical shifting factor (bT ) was fitted using a simple linear relationship that expresses the shifting factor as a function of temperature in the form: 𝐿𝑜𝑔(𝑏𝑇 ) = 𝑘1 ∙ 𝑇 + 𝑘2 (4.28) where 𝑘1 and 𝑘2 are material constants. The vertical shifting factors for the master curves were also fitted using a least square regression. The vertical shifting factors for the master curves were fitted according to the linear relationship (Eqn.4.28) using the least square regression. The coefficients of determination R2 of these fittings were in the range of 0.9688-0.9976, indicating excellent fitting (Fig.4.9). Table 4.3 compares the experimental vertical shifting factors and the vertical shifting factors calculated by the WLF model. Furthermore, the values of the constants C1 , C2 , 𝑘1 , and 𝑘2 are summarized in table 4.4. 101 (a) (b) (c) Figure 4.9 Vertical shifting factor fitting using a linear relationship for master curves in (a) MD, (b) TD, and (c) In-plane shear. Table 4.3 Predicted vertical shifting factors (𝐿𝑜𝑔(𝑏𝑇 )) using the linear relationship compared with experimental vertical shifting factors Machine Direction Transverse Direction In-Plane Shear (MD) (TD) T Experiment Predicted Experiment Predicted Experiment Predicted (ºC ) 20 0.00 0.00 0.00 0.00 0.00 0.00 30 0.03 0.033 0.037 0.028 0.016 0.018 40 0.07 0.068 0.065 0.072 0.08 0.06 50 0.1 0.10 0.13 0.12 0.10 0.09 60 0.14 0.138 0.175 0.16 0.14 0.13 102 Table 4.4 Model parameters for the horizontal and vertical shifting factors Parameters Machine Direction Transverse Direction In-Plane Shear (MD) (TD) C1 -4.74 -11.40 -10.56 C2 -88.93 -199.43 -140.73 𝑘1 0.0035 0.0044 0.0036 𝑘2 -1.0275 -1.3052 -1.0721 4.3.3 Prony Series Fitting The stress relaxation modulus of the separator in the temperature range of 20°C-60°C was expressed in terms of a Prony series. The parameters were determined by using the least squares curve fitting method implemented in MATLAB. The values of the Prony series parameters for the separator material in MD, TD, and in-plane shear are summarized in Table 4.5. To show how well the Prony series coefficients generated from the least square fitting compare with the experimental data, the curves for the fitted data and the master curves from the experimental data are compared in a log – log plot for all three material orientations and presented in Fig.4.10. Table 4.5 Prony series parameters Machine Direction (MD) Transverse Direction In-Plane Shear (TD) Relaxation Relaxation Relaxation Relaxation Relaxation Relaxation Coefficients time 𝝉𝒊 (s) Coefficient time 𝝉𝒊 (s) Coefficient time 𝝉𝒊 (s) Gi (Pa) Gi (Pa) Gi (Pa) 3.75 × 108 - 3.54 × 108 - 2.31 × 107 - 2.69 × 108 10−1 4.79 × 107 10−1 4.08 × 107 10−1 1.96 × 108 100 3.65 × 107 100 2.11 × 107 101 2.09 × 108 101 3.31 × 107 101 1.12 × 107 102 2.18 × 108 102 4.80 × 10 7 102 2.04 × 107 103 2.72 × 108 103 5.96 × 107 103 1.70 × 107 104 2.21 × 108 104 7.56 × 107 104 1.96 × 107 105 2.10 × 108 105 9.50 × 10 7 105 1.22 × 107 106 2.00 × 108 106 9.45 × 107 106 1.64 × 107 107 103 (a) (b) (c) Figure 4.10 Prony series fitting for the stress relaxation master curves in log – log scale for (a) MD, (b) TD, and (c) In-plane shear. 4.3.4 Determination of Strain-Dependent Nonlinear Parameters The current model assumes that the strain-dependent and temperature-dependent parameters are decoupled. This assumption has been used in literature for modeling the nonlinear thermomechanical response of viscoelastic materials [10,22,24-27]. Additionally, this assumption greatly simplified the experimental procedure for the determination of the model parameters. The nonlinearity in the current model is introduced through two strain-dependent nonlinear parameters that introduce nonlinearity in the equilibrium modulus (ℎ∞ ), and the transient stiffness (ℎ𝑘 ). As ℎ∞ and ℎ𝑘 are independent of temperature, at the reference temperature, the relaxation modulus is expressed as: 𝑛 𝑡 𝐺(𝑡, 𝜀) = ℎ∞ (𝜀) ∙ 𝐺∞ + ℎ𝑘 (𝜀) ∙ ∑ 𝐺𝑖 ∙ 𝑒𝑥𝑝 [− ] (4.29) 𝜏𝑖 𝑖=1 where 𝐺𝑖 , 𝐺∞ and 𝜏𝑖 , are the determined Prony series parameters. 104 The nonlinear parameters ℎ∞ and ℎ𝑘 were determined by least square curve fittings of the stress relaxation curves measured at the reference temperature. As these two parameters should not be negative, a non-negative criterion was enforced in the least square fitting procedure. The fittings for the determination of ℎ∞ and ℎ𝑘 in the MD, TD, and in-plane shear are shown in Fig.4.11. From Fig.4.11 it is evident that the curves from the fitting and the experimental results overlap nicely. (a) (b) (c) Figure 4.11 Curve fitting for stress relaxation curves at 20°C in the (a) MD, (b) TD, and (c) In- plane shear. The nonlinear parameters ℎ∞ and ℎ𝑘 were expressed by polynomial functions of strain, given in Eqn. 4.30 through Eqn. 4.32. Fig. 4.12 shows that the values computed by the polynomial functions are in good agreement with the values determined from the curve fitting. Furthermore, comparisons between the fitted nonlinear parameters and the parameters predicted by the polynomial functions are presented in tables 4.6 and 4.7 for ℎ∞ and ℎ𝑘 respectively. 105 ℎ = 1.25 × 106 𝜀 4−1.59 × 105𝜀 3 +7.41 × 103𝜀 2 − 1.46 × 102 𝜀. +1.25 𝑀𝐷 ∶ { ∞ (4.30) ℎ𝑘 = −1.25 × 104 𝜀 3 +1.18 × 103 𝜀 2 − 41.06𝜀. +1.08 ℎ = 7.95 × 105 𝜀 4 −1.15 × 105 𝜀 3 +6.10 × 103 𝜀 2 − 1.44 × 102 𝜀 + 1.38 𝑇𝐷 ∶ { ∞ (4.31) ℎ𝑘 = −9.61 × 105 𝜀 4 +1.27 × 105 𝜀 3 −5.81 × 103 𝜀 2 + 88.15𝜀 + 0.79 ℎ∞ = −1.05 × 106 𝜀 4 +1.16 × 105 𝜀 3−3.71 × 103 𝜀 2 + 15.93𝜀. +0.64 𝑆ℎ𝑒𝑎𝑟 : { (4.32) ℎ𝑘 = −6.39 × 105 𝜀 4+7.51 × 104𝜀 3 −2.80 × 103𝜀 2 + 22.62ε + 1.04 (a) (b) (c) Figure 4.12 Polynomial fitting for strain-dependent nonlinearity parameters ℎ∞ and ℎ𝑘 in the (a) MD, (b) TD, and (c) In-plane shear. 106 Table 4.6 Fitted values of the nonlinearity in equilibrium modulus (ℎ∞ ) compared with predicted values using the polynomial equations. Machine Direction Transverse Direction In-Plane Shear (MD) (TD) Strain Fitted Predicted Fitted Predicted Fitted Predicted (% ) value value value value value value 0.2-0.3 1.00 0.98 1.00 1.00 0.65 0.65 0.7 0.49 0.53 0.63 0.63 - - 1 0.38 0.38 0.45 0.45 0.54 0.53 1.5 0.29 0.24 0.26 0.25 - - 2 0.17 0.20 0.14 0.15 0.21 0.23 2.5 0.22 0.20 0.10 0.11 - - 3 0.18 0.21 0.12 0.10 0.06 0.04 4 0.22 0.20 0.06 0.07 0.03 0.04 5 0.26 0.26 0.04 0.04 0.02 0.02 Table 4.7. Fitted values of the combined nonlinearities in the transient stiffness and strain rate effect (ℎ𝑘 ) compared with predicted values using the polynomial equations. Machine Direction Transverse Direction In-Plane Shear (MD) (TD) Strain Fitted Predicted Fitted Predicted Fitted Predicted (% ) value value value value value value 0.2-0.3 1.00 1.00 1.00 1.01 1.08 1.08 0.7 0.84 0.84 1.18 1.16 - - 1 0.76 0.77 1.22 1.21 1.07 1.05 1.5 0.69 0.68 1.16 1.19 - - 2 0.64 0.63 1.11 1.10 0.84 0.87 2.5 0.59 0.59 0.98 0.98 - - 3 0.57 0.57 0.87 0.87 0.74 0.71 4 0.53 0.53 0.72 0.72 0.62 0.64 5 0.42 0.43 0.59 0.59 0.52 0.57 4.3.5 Calibration of Instrument Expansion The contribution of the instrument to the thermal strain during non-isothermal tests was determined by calibrating the instrument using an aluminium thin foil sample of known CTE as described in the experimental section. Using a sample with known CTE, it becomes possible to identify the thermal strains due to the instrument expansion from the total measured strain during the iso-stress test. The contribution of the instrument is determined by subtracting the predicted 107 dimensional change from the dimensional change measured by the DMA. The predicted dimensional change of the calibration sample was computed as: ∆𝑙𝑝𝑟𝑒𝑑𝑖𝑐𝑡𝑒𝑑 = 𝛼 ∙ 𝑙0 ∙ (𝑇 − 𝑇0 ) (4.33) where 𝑙0 is the initial gage length and 𝛼 is the CTE of pure aluminium (Al) determined to be 24.5 × 10−6/℃ over the temperature range of 20-200°C [34]. Fig.4.13 presents the displacement- temperature curves for the measured, predicted and calibrated displacements. Figure 4.13 Displacement versus temperature curves for the measured, predicted and calibrated displacements. The calibrated displacement vs temperature curve was fitted using a linear least square regression relationship to determine the rate at which the instrument grip and the drive shaft displace for every degree rise in temperature. From the result, the instrument displaces at a rate of 0.0013mm/°C and for an initial gage length of 15mm, the thermal strain contributed by the instrument for every 1°C increment is 0.0087%. The contribution of the instrument to the thermal strain will be accounted for in simulations for model validation in non-isothermal conditions. 4.4 Analytical Solutions for Stress Relaxation For step-loading cases, the hereditary integral in Eqn.4.1 may be evaluated analytically. In this sub-section, analytical solutions for stress relaxation are derived for multi-step strain loading at a constant temperature, and for multi-step temperature histories under a constant strain. These solutions are useful for verifying the implementation of the numerical model. 108 4.4.1 Single Step Stress Relaxation at Constant Temperature The TTSP master curves were fitted according to the generalized Maxwell model in Prony series. The mathematical expression for the stress relaxation response of the generalized maxwell model for a step strain loading is given by [14]: 𝑛 𝑡 𝐺(𝑡) = 𝐺∞ + ∑ 𝐺𝑖 ∙ exp (− ) (4.34) 𝜏𝑖 𝑖=1 where 𝐺∞ and 𝐺I are the Prony series coefficients, t is the time, and 𝜏𝑖 are the characteristic relaxation times, respectively. With the introduction of nonlinear parameters, the stress relaxation modulus as a function of time and strain at the reference temperature is given by: 𝑛 𝑡 𝐺(𝑡, 𝜀) = ℎ∞ (𝜀) ∙ 𝐺∞ + ℎ𝑘 (𝜀) ∙ ∑ 𝐺𝑖 ∙ 𝑒𝑥𝑝 [− ] (4.35) 𝜏𝑖 𝑖=1 A reduced time t’ is used to account for the temperature dependence through the horizontal shifting factor, such that: 𝑡 𝑡′ = 𝑎 (4.36) 𝑇 Replacing t in Eqn.4.35 with t’, and introducing the vertical shifting factor (bT ) to Eqn.4.35, we have the expression for the relaxation modulus as a function of time, strain and temperature: 𝑁 1 𝑡 𝐺(𝑡, 𝜀, 𝑇) = ∙ {ℎ∞ (𝜀) ∙ 𝐺∞ + ℎ𝑘 (𝜀) ∙ ∑ 𝐺𝑖 ∙ 𝑒𝑥𝑝 [− ]} (4.37) 𝑏𝑇 (𝑇) 𝑎 𝑇 ∙ 𝜏𝑖 𝑖=1 Furthermore, the stress relaxation as a function of time, strain and constant temperature can be calculated from the relaxation modulus expressed above by multiplying by the desired constant strain level, 𝜀0 . Therefore, the analytical solutions for the stress relaxation accounting for temperature dependence for single step loading is given below: 𝑁 1 𝑡 𝜎(𝑡) = 𝜀0 ∙ [ ∙ {ℎ∞ (𝜀) ∙ 𝐺∞ + ℎ𝑘 (𝜀) ∙ ∑ 𝐺𝑖 ∙ 𝑒𝑥𝑝 [− ]}] (4.38) 𝑏𝑇 (𝑇) 𝑎 𝑇 ∙ 𝜏𝑖 𝑖=1 109 4.4.2 Step loading at constant temperature From the linear viscoelastic theory [23], the total stress 𝜎(𝑡) at a given time t, is given by: 𝜎(𝑡) = ∆𝜀1 ∙ 𝐺(𝑡 − 𝑡1 ) + ∆𝜀2 ∙ 𝐺(𝑡 − 𝑡2 ) + ∙∙∙∙∙ +∆𝜀𝑛 ∙ 𝐺(𝑡 − 𝑡𝑛 ) (4.39) Eqn. 4.39 gives the basis the Boltzmann superposition principle, where 𝐺(𝑡 − 𝑡1 ) is the stress relaxation modulus, and ∆𝜀1, ∆𝜀2 and ∆𝜀𝑛 are the strain increments added at times 𝑡1 , 𝑡2 , and 𝑡𝑛 . For loading cases with strains larger than the linear viscoelastic limit, the relaxation modulus becomes a nonlinear function which is dependent on strain and temperature. Extending the linear superposition to the nonlinear deformation range, for the one-step loading case at a constant temperature T, the total stress is given by: 𝜎(𝑡) = ∆𝜀1 ∙ 𝐺(𝑡 − 𝑡1 , 𝜀1 , 𝑇) (4.40) For the multi-step loading case, the total stress at time 𝑡 is given by: 𝜎(𝑡) = ∆𝜀1 ∙ 𝐺(𝑡 − 𝑡1 , 𝜀1 , 𝑇) + ∆𝜀2 ∙ 𝐺(𝑡 − 𝑡2 , 𝜀2 , 𝑇) + ∙∙∙∙∙ +∆𝜀𝑛 ∙ 𝐺(𝑡 − 𝑡𝑛 , 𝜀𝑛 , 𝑇) (4.41) Where the stress relaxation modulus at a given time interval, as a function of strain and temperature (𝐺(𝑡 − 𝑡𝑖 , 𝜀𝑖 , 𝑇)) is expressed by Eqn.4.37. 4.4.3 Step temperature history at constant strain level In section 3.4.3 of chapter 3, the case where the temperature loading is applied stepwise and the temperature values 𝑇1, 𝑇2, up to 𝑇𝑛 for n steps are added at times 𝑡1 , 𝑡2 , up to 𝑡𝑛 was considered. It was concluded that to develop the analytical solution for stress relaxation with a step temperature history, the reduced time (𝜌(𝑡)) had to be evaluated for the time intervals at which the temperature steps are held. The total time was broken down into intervals: 𝑡 ∈ [𝑡1 , 𝑡2 ], 𝑡 ∈ [𝑡2 , 𝑡3 ], … , and 𝑡 > 𝑡𝑛 , where n is the total number of time segments. For each time interval, the reduced time was calculated and presented as Eqns. 3.36 through 3.39 in chapter 3. However, the temperature dependence introduced in the previous chapter did not consider the vertical shifting factor. Hence, to complete the development of the analytical solution for stress relaxation with step temperature history, the vertical shifting factor needs to be considered. From the evaluation of the reduced time at the different time segments (Ch.3), it is evident that the reduced time for temperature steps other than the first step carries the fading memory effect of the temperature response at previous steps. In other words, the part of the reduced time due to 110 𝑡−𝑡𝑖 a recent temperature increment [ 𝑎 ] will contribute more to the stress relaxation than parts of the 𝑇𝑖 𝑡𝑖+1 −𝑡𝑖 reduced time due to temperature increments in the distant past [ ]. 𝑎𝑇 𝑖 Furthermore, the vertical shift needed to match the reduced time with fading memory effect at temperature steps other than the first step is the change in the vertical shift from the previous temperature step to the current. This can be represented mathematically as: 𝐿𝑜𝑔(𝑏𝑇𝑖+1 ) = 𝐿𝑜𝑔(𝑏𝑇𝑖 ) ± 𝐿𝑜𝑔(∆𝑏𝑇𝑖→𝑖+1 ) (4.42) where 𝐿𝑜𝑔(∆𝑏𝑇𝑖→𝑖+1 ) is the change in the vertical shift required to go from the stress relaxation response at one temperature step to another and i=0,1, 2...n. Reducing Eqn. 4.42 by eliminating the logarithmic expression, the change in the vertical shift can be represented as: 𝑏𝑇𝑖+1 ∆𝑏𝑇𝑖→𝑖+1 = (4.43) 𝑏𝑇𝑖 Hence, the analytical solution for the stress relaxation as a function of time, strain and temperature with step histories is given as: 𝑛 1 𝜌(𝑡) 𝜎(𝑡, 𝜀, 𝑇) = 𝜀0 ∙ ∙ {ℎ∞ (𝜀) ∙ 𝐺∞ + ℎ𝑘 (𝜀) ∙ ∑ 𝐺𝑖 ∙ 𝑒𝑥𝑝 [− ]} (4.44) ∆𝑏𝑇 (𝑇) 𝜏𝑖 𝑖=1 With temperature change, the thermal strains due to thermal expansion have to be accounted for and introduced into the analytical solution. The thermal strain is expressed as: 𝜀𝑡ℎ = 𝛼 ∙ (𝑇 − 𝑇0 ) (4.45) where 𝛼 is the coefficient of thermal expansion (CTE). Hence the analytical solution for stress relaxation accounting for thermal expansion effect takes the form: 𝑛 1 𝜌(𝑡) 𝜎(𝑡, 𝜀, 𝑇) = [𝜀0 − 𝜀𝑡ℎ ] ∙ ∙ {ℎ∞ (𝜀) ∙ 𝐺∞ + ℎ𝑘 (𝜀) ∙ ∑ 𝐺𝑖 ∙ 𝑒𝑥𝑝 [− ]} (4.46) ∆𝑏𝑇 (𝑇) 𝜏𝑖 𝑖=1 4.5 Orthotropic Nonlinear Thermoviscoelastic Model Verification The orthotropic nonlinear thermoviscoelastic material model was written as a user material model for shell element simulations in LS-DYNA®. An executable file was then compiled with the addition of the user material model to run the simulations. In this section, the implemented model is verified using the analytical solutions for stress relaxation for one and two-step loading cases at constant temperatures of 20°C and 30°C, and for non-isothermal cases with a step temperature history of 30°C-40°C-50°C-60°C at a constant strain level. 111 All simulations were performed with a one-element model, with boundary conditions as described in Chapter 3. Also, in these simulations, the required temperature history was applied at all four nodes of the shell element. Explicit time integration was used in all simulations with the time, scaled so that 1ms will correspond to 1s in real-time. For non-isothermal simulations, the thermal expansion was also introduced into the simulation through *MAT_ADD_THERMAL_EXPANSION in the LS-DYNA® keyword [35] as in the orthotropic linear thermoviscoelastic model. The CTE of Celgard®2400 for samples cut along the MD and TD as functions of temperature have been determined and are summarized in the work by Yan et. al. [13] and are shown in Table 3.11. Furthermore, the contribution of the testing instrument to the total thermal strain was also introduced in the non-isothermal simulations for model validation. Figure 4.15 presents the model verification in the MD and TD at 20ºC and 30ºC with analytical solutions computed using Eqns.4.40 and 4.41 for step loadings at constant temperatures. For the one-step loading case, the stress relaxation was monitored at a constant strain level of 1% and for the two-step loading case, the strain loading varied according to Fig.4.14. As shown, the simulations and analytical solutions coincided for all cases, indicating the model has been implemented correctly. (a) (b) Figure 4.14 Strain loading for two-step stress relaxation for (a) low to high and (b) high to low cases. 112 (a) (b) (c) (d) (e) (f) Figure 4.15 Model verification with analytical solutions for one-step, two-step (low to high), and two-step (high to low) loading cases respectively in the (a-c) MD and (d-f) TD. 113 The developed nonlinear thermomechanical user material model was also verified using analytical solutions computed using Eqns.4.44 and 4.46 for stress relaxation with step temperature histories with the applied strain level kept constant at 1%. The applied temperature history is shown in Fig.4.16. The analytical solutions were computed for cases with and without considering the thermal expansion effect (i.e., thermal stains due to sample expansion alone). From Fig.4.17, the predicted stress histories from the simulations match with the stress histories calculated from the analytical solutions. The results also show that the numerical model was correctly implemented. Figure 4.16 Applied temperature history. (a) (b) Figure 4.17 Model verification with the analytical solutions for stress relaxation with step temperature history with and without considering thermal expansion in the material response in the (a) MD and (b) TD. 114 4.6 Orthotropic Nonlinear Thermoviscoelastic Model Validation This section covers the model validation by comparing simulations with experimental results from stress relaxation tests at different strain levels and temperature histories, creep tests carried out at different loading levels, and tensile tests at different displacement rates, as summarized in Table 4.1. Stress relaxation tests were carried out on samples cut along the MD, TD, and 45° off-axis direction at different strain levels above the linear viscoelastic limit and constant temperature values ranging from 20ºC-60ºC. The stress histories predicted by the model compare well with the experimental results, as shown in Fig. 4.18. Furthermore, the excellent agreement between the experimental results and the model predictions in Fig.4.18(d) verifies the applicability of the transformation equation (Eqn.4.21) based on the elastic-viscoelastic correspondence principle in the nonlinear deformation region considered in this work. (a) (b) Figure 4.18 Comparison of the stress history predicted by LS-DYNA® simulation using the orthotropic nonlinear thermoviscoelastic material model with experimental data for stress relaxation tests in the (a) MD at 30ºC and 1% strain: 40ºC and 1% strain; 60ºC and 0.7% strain, (b) TD at 30ºC and 1%, 2%, and 3% strain levels, (c) TD at 40ºC and at 1%, and 2% strain levels, and (d) 45º off-axis direction at 20ºC and 0.3%-5% strain levels. 115 Figure 4.18 (Cont’d) (c) (d) Furthermore, the implemented model is validated against experimental results from stress relaxation tests performed under non-isothermal conditions following a multi-step temperature history, as shown in Fig.4.16. The experiments were performed at constant strain levels of 1% and 0.7% in the MD and TD. Simulations were carried out for stress relaxation with the recorded temperature history. To examine the thermal expansion effects, the simulations were performed with and without considering the CTE. Without the consideration of CTE, the simulation results overpredicted the stress histories from the experimental procedures. The results (Fig.4.19) show that with the consideration of the thermal expansion effects, the model predictions of stress histories in non-isothermal conditions agree well with experimental data. 116 (a) (b) (c) (d) Figure 4.19 Comparison of the stress histories predicted with and without considering CTE by LS-DYNA® simulation using the orthotropic nonlinear thermoviscoelastic material model with experimental data for non-isothermal stress relaxation with a multi-step temperature history of 30ºC-40ºC-50ºC-60ºC in the MD at (a) 0.7% strain level, (b) 1% strain level, and in the TD at (c) 0.7% strain level and (d) 1% strain level. The model is also validated against creep tests carried out at different stress levels in the nonlinear range at 20ºC. Figure 4.20 shows the comparison between predicted and experimental strain histories under constant stress levels of 4, 5, 8 and 10MPa in the MD, and 3 and 5MPa in the TD. The results show a good agreement between simulation results and the experimental data. 117 (a) (b) Figure 4.20 Comparison of the strain history predicted by LS-DYNA® simulation using the orthotropic nonlinear thermoviscoelastic material model with experimental data for creep tests at 20ºC in the (a) MD and (b) TD. Uniaxial tensile tests were carried out on specimens at different displacement rates and constant temperatures as summarized in Table 4.1. Figure 4.21 presents a comparison between the experimental and the simulation results for the uniaxial stress-strain response in the MD and TD. The experimental results showed that with increasing the displacement rate, the material showed higher stiffness in its stress-strain response. With increasing temperature, a lower stiffness in the material response was the case. Simulation results using the orthotropic linear viscoelastic model developed in the previous chapter were introduced in Fig 4.21 to show the discrepancy in accuracy of the predictions of the linear and nonlinear viscoelastic models in the large deformation region for samples cut along the MD and TD respectively. From the results, the nonlinear viscoelastic deformation range of the material was observed to be up to 1.5% strain in the TD. With an increase in temperature, the material softened and the onset of permanent deformations and yielding began to occur at lower strains for samples cut along the TD. However, the stress-strain behavior of the samples cut along the MD showed no distinctive yield point. The MD samples continuously deformed with no signs of yielding. The results show that the stress-strain curves predicted by the orthotropic nonlinear thermoviscoelastic model agree well with the experimental results and capture the trend in the rate and temperature dependence of the stress-strain response of the material in the nonlinear viscoelastic deformation range. 118 (a) (b) (c) (d) (e) Figure 4.21 Comparison of the stress-strain response predicted by LS-DYNA® simulation using the orthotropic nonlinear thermoviscoelastic material model with experimental data from uniaxial tensile tests carried out at displacement rates of 0.001mm/s and 0.01mm/s rates for samples cut along the MD at (a) 30ºC, (b) 40ºC, and the TD at (c) 30ºC, (d) 40ºC, and (e) 50ºC. 119 4.7 Summary The development of an orthotropic nonlinear thermoviscoelastic model for predicting the thermomechanical response of polymeric battery separators in thermal ramp scenarios has been presented in this work. A discretization algorithm was employed to evaluate the nonlinear viscoelastic hereditary integral with a kernel of Prony series based on a generalized Maxwell model with nonlinear springs and dashpots. Temperature dependence was introduced into the model through the time-temperature superposition principle (TTSP). The model parameters were determined for Celgard®2400, which is a porous polypropylene (PP) separator. The model was implemented in LS-DYNA® finite element (FE) package as a user material model. Analytical solutions for stress relaxation with step loading and step temperature history were formulated based on the viscoelastic theory to verify the implementation of the model. The developed model was validated against stress relaxation tests, creep tests, uniaxial tensile tests at constant temperatures, and stress relaxation at non-isothermal conditions accounting for the thermal expansion of the material. The results show that the model predictions of the material anisotropy, rate dependence, and temperature dependence of the separator under large deformations in the nonlinear viscoelastic region agree well with the experimental data. Furthermore, it was observed that the nonlinear viscoelastic deformation region of polymeric materials like the separator investigated in this work is considerably larger compared to the linear viscoelastic region, especially for samples cut along the MD. Hence, the nonlinear recoverable response of the polymeric material could be mistakenly considered irrecoverable if the linear viscoelastic theory is solely used. The validation results also confirm that assuming the strain and temperature-dependent parameters are decoupled was sufficient to describe the wide range of thermomechanical behaviors examined in the current work. This model is limited to predictions of the thermomechanical response of the separator under large deformations before yielding, the coupling of the nonlinear thermoviscoelastic model with a viscoplastic model to capture and predict the separator response under large irrecoverable deformations will be addressed in the subsequent chapter. 120 REFERENCES 1. R. C. Ihuaenyi, J. Deng, C. Bae, and X. Xiao. An orthotropic nonlinear thermoviscoelastic model for polymeric battery separators. Journal of the Electrochemical Society, 170, 010520, 2023. 2. S. Yan, J. Deng, C. Bae, S. Kalnaus, and X. Xiao, X. Orthotropic viscoelastic modeling of polymeric battery separator. Journal of The Electrochemical Society, 167(9), p.090530, 2020. 3. M. A Puso and J. A Weiss. Finite element implementation of anisotropic quasilinear viscoelasticity using a discrete spectrum approximation. Journal of Biomechanical Engineering, 120:62-70, 1998. 4. R. A. Schapery. On the characterization of nonlinear viscoelastic materials. Polymer Engineering & Science, 9(4), pp.295-310, 1969. 5. S. Y. Zhang amd X.Y. Xiang. Creep characterization of a fiber reinforced plastic material. Journal of reinforced plastics and composites, 11(10), pp.1187-1194, 1992. 6. M. Katouzian, O.S. Bruller and A. Horoschenkoff. On the effect of temperature on the creep behavior of neat and carbon fiber reinforced PEEK and epoxy resin. Journal of composite materials, 29(3), pp.372-387, 1995. 7. A. Horoschenkoff. Characterization of the creep compliances J22 and J66 of orthotropic composites with PEEK and epoxy matrices using the nonlinear viscoelastic response of the neat resins. Journal of composite Materials, 24(8), pp.879-891, 1990. 8. J. Henderson, G. Calderon, and J. Rand. A nonlinear viscoelastic constitutive model for balloon films. In 32nd Aerospace Sciences Meeting and Exhibit, p. 638, 1994. 9. I.M. Daniel and O. Ishai, Engineering Mechanics of Composite Materials (Oxford University Press, Oxford, UK), 2006. 10. F. Bosi and S. Pellegrino. Nonlinear thermomechanical response and constitutive modeling of viscoelastic polyethylene membranes. Mechanics of Materials, 117, pp.9-21, 2018. 11. J. Rand, D. Grant and T. Strganac. The nonlinear biaxial characterization of balloon film. In 34th Aerospace Sciences Meeting and Exhibit, p. 574, 1996. 12. S. Yan, J Deng, C. Bae, Y. He, A. Asta, and X. Xiao. In-plane orthotropic property characterization of a polymeric battery separator. Polymer Testing, 72, pp.46-54, 2018. 13. S. Yan, J. Deng, C. Bae, and X. Xiao. Thermal expansion/shrinkage measurement of battery separators using a dynamic mechanical analyzer. Polymer Testing, 71, pp.65-71, 2018. 14. A. S. Wineman, K. R. Rajagopal, Mechanical Response of Polymers: An Introduction (Cambridge university press, Cambridge, UK), 2000. 15. R. M. Christensen, Theory of Viscoelasticity: An Introduction (Academic Press, New York, NY) 2nd ed., 1982. 16. R. C. Ihuaenyi, S. Yan, J. Deng, C. Bae, C., I. Sakib, and X. Xiao. Orthotropic Thermo- Viscoelastic Model for Polymeric Battery Separators with Electrolyte Effect. Journal of The Electrochemical Society, 168(9), p.090536, 2021. 121 17. G. A. C. Graham. The correspondence principle of linear viscoelasticity theory for mixed boundary value problems involving time-dependent boundary regions. Quarterly of Applied Mathematics, 26(2), pp.167-174, 1968. 18. S. Mukherjee and G.H. Paulino. The elastic-viscoelastic correspondence principle for functionally graded materials, revisited. J. Appl. Mech., 70(3), pp.359-363, 2003. 19. L. Khazanovich. The elastic–viscoelastic correspondence principle for non-homogeneous materials with time translation non-invariant properties. International Journal of Solids and Structures, 45(17), pp.4739-4747, 2008. 20. E. V. Dave, G. H. Paulino and W. G. Buttlar. Viscoelastic functionally graded finite-element method using correspondence principle. Journal of Materials in Civil Engineering, 23(1), pp.39-48, 2011. 21. C. Zhang and W. Zhang. Elasticity recovery correspondence principles for physically nonlinear viscoelastic problems for a class of materials. International journal of solids and structures, 38(46-47), pp.8359-8373, 2001. 22. S. Sawant and A. Muliana. A thermo-mechanical viscoelastic analysis of orthotropic materials. Composite Structures, 83(1), pp.61-72, 2008. 23. R. M. Haj-Ali and A. H. Muliana. A multi-scale constitutive formulation for the nonlinear viscoelastic analysis of laminated composite materials and structures. International Journal of Solids and Structures, 41(13), pp.3461-3490, 2004. 24. M. E. Tuttle, A. Pasricha, A. F. Emery.The nonlinear viscoelastic-viscoplastic behaviour of IM7/5260 composites subjected to cyclic loading. Journal of Composite Materials, vol. 29, pp. 2025-2046, 1995. 25. D. Peretz and Y. Weitsman. Nonlinear viscoelastic Characterization of FM-73 Adhesive,” Journal of Rheology, vol. 26, no. 3, pp. 245-261, 1982. 26. D. Peretz and Y. Weitsman. The Nonlinear Thermoviscoelastic Characterizations of FM-73 Adhesives,” Journal of Rheology, vol. 27, iss. 2, pp. 97-114, 1983. 27. M. Jafaripour and F. Taheri-Behrooz. Creep behavior modeling of polymeric composites using Schapery model based on micro-macromechanical approaches. European journal of mechanics-A/Solids, 81, p.103963, 2020. 28. S. E. Boyd, J. J. Lesko, and S. W. Case, The Thermo-Viscoelastic, Viscoplastic Characterization of Vetrotex 324/Derakane 510A-40 through Tg. Journal of Engineering Materials and Technology, Transactions of the ASME. 128(4): p. 586-594, 2006. 29. J. D. Ferry, Viscoelastic Properties of Polymers (Wiley, New York, NY), 1980. 30. I. M. Ward and J. Sweeney, An Introduction to the Mechanical Properties of Solid Polymers (Wiley, New York, NY), 2004. 31. H. F. Brinson, L. C. Brinson, Polymer Engineering Science and Viscoelasticity: An Introduction (Springer, Berlin), 2008. 32. J. Xu, L. Wang, J. Guan, and S. Yin. Coupled effect of strain rate and solvent on dynamic mechanical behaviors of separators in lithium-ion batteries. Materials & Design, 95, pp.319- 328, 2016. 122 33. R. F. Landel and L. E. Nielsen, Mechanical properties of polymers and composites (CRC press), 1993. 34. P. Hidnert and H.S. Krider. Thermal expansion of aluminum and some aluminum. J Res Natl Bur Stand, 48, p.209, 1952. 35. LS-DYNA Keyword Manual, Livermore Software Technology Corporation (LSTC). Livermore, California, 2018. 123 Chapter 5 Coupled Viscoelastic – Viscoplastic Modeling This chapter presents the development of a phenomenological viscoplastic constitutive model for polymeric battery separators. The developed viscoplastic model is coupled with the nonlinear viscoelastic model developed in the previous chapter (Ch.4) to predict the stress-strain response of the considered separator material in the entire range of its deformation before the onset of failure. The viscoplastic model is developed based on a rheological framework that considers the mechanisms involved in the initial yielding, change in viscosity, strain softening and strain hardening of polymeric separators. The developed model accounts for the material anisotropy, rate and temperature dependence. The model parameters are determined for a polypropylene (PP) separator and are implemented in LS-DYNA® finite element (FE) package as a user-defined subroutine. Simulations of uniaxial tensile stress-strain response at different strain rates and temperatures are compared with experimental data to validate the model predictions. 5.1 Viscoplastic Model Overview The mechanical response of polymeric materials under large deformations beyond the onset of yielding has been widely studied and shown to be rate-dependent. Hence to accurately predict the response of polymeric materials, viscoplastic models have to be considered. Phenomenological viscoplastic constitutive models have been applied widely in literature to model the temperature and rate-dependent large strain response of polymeric materials [1-5]. Hence, in this chapter, a phenomenological constitutive model is developed to predict the strain rate and temperature-dependent post-yield response of polymeric battery separators. To accomplish this, a rheological framework is introduced as the basis of the phenomenological viscoplastic constitutive model. The rheological framework is proposed with the approach that the stress evolution in the material beyond the yield threshold follows a governing relationship that can be simplified using a combination of a sliding frictional element, a dashpot and two springs as shown in Fig.5.1 below. 124 Figure 5.1 Rheological framework for viscoplastic model. From Fig. 5.1, the rheological model consists of a sliding frictional element to represent the yield threshold, a dashpot to account for the changes in the viscosity of the material due to the accumulation of plastic strains and changes in the strain rate, and two springs to account for the strain softening and strain hardening mechanisms. The elements in the rheological model are arranged in such a way that the total stress can be determined by a summation of the stresses in the individual elements. Hence, after the yield threshold is exceeded, the stress as a function of the plastic strain and the constant strain rate is expressed as: 𝜎(𝜀𝑝 , 𝜀̇) = 𝜎𝑦0 (𝜀̇) + 𝜇𝑣𝑝 (𝜀̇) ∙ {1 − 𝑒𝑥𝑝(−𝑘(𝜀̇) ∙ 𝜀𝑝 )} + 𝐻1 (𝜀̇) ∙ 𝜀𝑝 + 𝐻2 (𝜀̇) ∙ 𝜀𝑝𝑚 (5.1) where 𝜎𝑦0 is the initial value of the stress at the onset of permanent deformations or the initial yield stress. 𝜎𝑦0 is strain rate and temperature dependent. With an increase in strain rate and a decrease in temperature, the initial yield stress increases in value and vice versa. 𝜇𝑣𝑝 (𝜀̇) and 𝑘(𝜀̇) are the first and second viscosity parameters. 𝜇𝑣𝑝 (𝜀̇) accounts for the rate of the change in the viscosity of the material in the large deformation region on the flow stress due to changes in strain rate and temperature while 𝑘(𝜀̇) regulates the rate at which softening and or hardening occurs after the initial yield stress has been exceeded. 𝐻1 (𝜀̇) is the rate-dependent softening coefficient that accounts for the decrease in flow stress with an increase in strain. 𝐻2 (𝜀̇) and 𝑚 are the hardening coefficients that combine to depict the increase in flow stress under a further increase in strain. To implement the coupled viscoelastic – viscoplastic model, the calibration of the viscoplastic model parameters has to be completed. Section 5.2 presents the experimental procedure for uniaxial tensile testing of the selected separator material carried out at different strain rates and temperatures. Furthermore, the procedures for the determination of the strain rate and 125 temperature-dependent model parameters using uniaxial tensile tests are presented in sections 5.3 and 5.4 respectively. The model implementation procedure and validation are presented in section 5.5 and the summary of the findings is given in section 5.6. 5.2 Uniaxial Tensile Testing Uniaxial tensile tests were carried out on the porous PP Celgard ®2400 samples at strain rates of 0.1/s, 0.01/s, 0.001/s, and 0.0001/s. The nominal dimensions of the specimen were as described in Chapter 3. The tensile tests were also carried out at temperatures ranging from 20°C to 50°C. For the tests at 20°C using samples cut along the MD and TD of the polymeric separator. A liquified nitrogen (LN2) tank was connected to the RSA G2 rheometer as shown in Fig.5.2(b), to cool each tested sample to the desired temperature. The experimentally determined stress-strain curves for the samples cut along the MD and TD are presented in Figs.5.3 and 5.4 respectively. Furthermore, the typical post-test geometry of the samples cut along the MD and TD are presented in Fig.5.5 to show the deformation mechanisms the samples undergo under uniaxial tension and large strains. (a) (b) Figure 5.2 Experimental set up for uniaxial tensile tests: (a) sample fixture (b) liquified nitrogen tank connected to RSA-G2. 126 (a) (b) (c) (d) Figure 5.3 Uniaxial tensile stress-strain curves for samples cut along the MD at different strain rates and temperatures of: (a) 20°C, (b) 30°C, (c) 40°C, and (d) 50°C. (a) (b) Figure 5.4 Uniaxial tensile stress-strain curves for samples cut along the TD at different strain rates and temperatures of: (a) 20°C, (b) 30°C, (c) 40°C, and (d) 50°C. 127 Figure 5.4 (Cont’d) (c) (d) (a) (b) Figure 5.5 Typical post-test geometry of samples under tension in the (a) MD and (b) TD. From the results, the effects of strain rate and temperature are manifested in the stress- strain response of the separator in the MD and TD. With increasing strain rate, there is an observed increase in the stiffness of the material and vice versa. Also, with increasing temperature, the material softened and there was an observed reduction in stiffness. Furthermore, the results also show that the stress-strain curves for the samples cut along the TD are that of a typical semicrystalline polymer with five major phases in its deformation process. The first phase consists of fully recoverable linear viscoelastic deformation. In phase two, the material transitions into a non-linear viscoelastic deformation region. Phase three characterizes the yield mechanism after which the deformation becomes irrecoverable. Phases four and five are strain softening and strain 128 hardening respectively. For semicrystalline polymers, strain softening is attributed to broken crystallization and strain hardening related to recrystallization at large strains [5]. Also, a neck formation is observed in the TD sample deformation. The necking region is characterized by a rapid decrease in the cross-sectional area at a particular point along the gage length of the specimen. This necking phenomenon is expressed in the posttest sample geometry in Fig.5.5(b). However, the material response in the MD is different from that of a typical semicrystalline material in the sense that there is no clear yielding mechanism following phases one and two which constitutes linear and non-linear viscoelastic response. Furthermore, in the MD, there is an observed continuous hardening in the material response under very large strains. These distinct differences in the uniaxial stress-strain response of the material in the MD and TD introduce complexities in the constitutive modelling of polymeric separators under large strains. These complexities have to be taken into account. Furthermore, the uniaxial tensile test results give us an understanding of the total stress-strain behavior of the separator material in its different orientations and present an excellent case for model validation. 5.3 Identification of Strain Rate-Dependent Model Parameters at Constant Temperature The viscoplastic model parameters at a constant temperature are calibrated using uniaxial tensile test results carried out at different strain rates at 20°C. This temperature value is selected intentionally as it is the reference temperature for the viscoelastic model that will be coupled with the viscoplastic model to predict the stress-strain response of the investigated polymeric separator material. The parameters to be determined are the initial yield stress, viscosity coefficients, strain hardening and strain softening coefficients. The initial yield stress was determined based on the need to identify it within the boundaries of the range of prediction of the nonlinear viscoelastic model (0 - 5% strain). Also, identifying the initial yield point of polymeric materials is not as straightforward as it is with metals. Hence, the offset method was employed to identify the initial yield stresses in the MD and TD of the polymeric separator at different strain rates. From the test results for samples cut along the MD, the yield stresses were identified by taking a 1% offset equivalent to a value of 2.5% strain at yield. While the yield stresses for samples cut along the TD were determined by taking a 0.3% offset equivalent to 1.5% strain at yield. The rationale behind taking a larger offset or strain value at yield in the MD than in the TD is supported by the superiority in the mechanical response of the material in the MD compared to that in the TD. Hence, it is expected that the material will yield at lower strain values in the TD than in the MD. 129 To determine the other model parameters (𝜇𝑣𝑝 (𝜀̇) , 𝑘(𝜀̇), 𝐻1 (𝜀̇), and 𝐻2 (𝜀̇)), an optimization fitting method was scripted in python to minimize the error between the fitted results and experimental data and generate the model parameters at each constant strain rate value. The fitting was carried out and the determined model parameters for the MD and TD responses are summarized in tables 5.1 and 5.2 respectively. Comparisons between the experimental data and the fitted curves are shown in Fig.5.6 to show the quality of the implemented fitting method. Table 5.1 Viscoplastic model parameters for MD response Strain rate 𝝈𝒚𝟎 𝝁𝒗𝒑 𝒌 𝑯𝟏 𝑯𝟐 𝒎 (𝒔−𝟏 ) (𝐌𝐏𝐚) (𝐌𝐏𝐚) (𝐌𝐏𝐚) (𝐌𝐏𝐚) 0.0001 19.33 22.19 25.21 0 228.77 1 0.001 21.20 18.11 34.44 0 281.24 1 0.01 25.61 16.84 41.29 0 307.71 1 0.1 34.45 9.70 115.36 0 343.93 1 Table 5.2 Viscoplastic model parameters for TD response Strain rate 𝝈𝒚𝟎 𝝁𝒗𝒑 𝒌 𝑯𝟏 𝑯𝟐 𝒎 (𝒔−𝟏 ) (𝐌𝐏𝐚) (𝐌𝐏𝐚) (𝐌𝐏𝐚) (𝐌𝐏𝐚) 0.0001 7.39 6.06 51.25 -4.79 2.76 2 0.001 8.61 6.13 52.74 -5.98 4.24 2 0.01 9.83 6.74 63.94 -8.06 7.85 2 0.1 10.33 7.35 77.86 -8.74 8.49 2 (a) (b) Figure 5.6 Curve fitting for identifying viscoplastic model parameters at 20°C in the (a) MD and (b) TD. 130 The curve fitting results from Fig.5.6 shows how nicely the fitted curves overlap with the experimental curves. This verifies the accuracy in the determination of the model parameters given in Tables 5.1 and 5.2. However, from Tables 5.1 and 5.1, we have unique sets of values for each viscoplastic model parameter at different strain rate values. To implement the viscoplastic model in predictive modeling, given arbitrary strain rate values within the considered range, equations relating the viscoplastic parameters as functions of strain rate have to be defined. From the current data generated, we have to define 5 equations for the TD and 4 equations for the MD model parameters respectively. To reduce the number of relationships generated for the model parameters and the number of input parameters for the model implementation overall, a parameter reduction method is applied. To reduce the number of parameters required, a strain rate scaling factor (𝑎𝜀̇ (𝜀̇)) is introduced to modify Eqn.5.1 such that the stress evolution after the yield threshold is exceeded at the reference temperature is expressed as: 0 𝜎(𝜀𝑝 , 𝜀̇) = 𝜎𝑦0 (𝜀̇) + {𝜇𝑣𝑝 ⋅ {1 − 𝑒𝑥𝑝(−𝑘(𝜀̇) ∙ 𝜀𝑝 )} + 𝐻10 ∙ 𝜀𝑝 + 𝐻20 ∙ 𝜀𝑝𝑚 } ⋅ 𝑎𝜀̇ (𝜀̇) (5.2) 0 where 𝜇𝑣𝑝 , 𝐻10 , and 𝐻20 are the rate of change in viscosity, strain softening and strain hardening coefficients at the reference strain rate (0.0001/s). Furthermore, the reduced viscoplastic model parameters for the MD and TD responses are summarized in Tables 5.3 and 5.4 respectively. Figure 5.7 below also shows how well the fitted curves produced using the reduced model parameters compare with the experimental data. Table 5.3 Reduced viscoplastic model parameters for MD response Strain 𝝈𝒚𝟎 𝝁𝒗𝒑 𝒌 𝑯𝟏 𝑯𝟐 𝒎 𝒂𝜺̇ rate (𝐌𝐏𝐚) (𝐌𝐏𝐚) (𝐌𝐏𝐚) (𝐌𝐏𝐚) (𝒔−𝟏 ) 0.0001 19.33 22.19 25.21 0 228.77 1 1 0.001 21.20 22.19 25.21 0 228.77 1 1.09 0.01 25.61 22.19 25.21 0 228.77 1 1.16 0.1 34.45 22.19 25.21 0 228.77 1 1.12 Table 5.4 Reduced viscoplastic model parameters for TD response Strain 𝝈𝒚𝟎 𝝁𝒗𝒑 𝒌 𝑯𝟏 𝑯𝟐 𝒎 𝒂𝜺̇ rate (𝐌𝐏𝐚) (𝐌𝐏𝐚) (𝐌𝐏𝐚) (𝐌𝐏𝐚) (𝒔−𝟏 ) 0.0001 7.39 6.06 51.25 -4.79 2.76 2 1 0.001 8.61 6.06 52.74 -4.79 2.76 2 0.98 0.01 9.83 6.06 63.94 -4.79 2.76 2 1.06 0.1 10.33 6.06 77.86 -4.79 2.76 2 1.16 131 (a) (b) Figure 5.7 Curve fitting for reduced viscoplastic model parameters at 20°C in the (a) MD and (b) TD. From the implementation of the parameter reduction methodology, it is clear that a combination of the model parameters at the reference temperature along with the strain rate scaling factor is enough to fit the experimental curves with good agreement. Furthermore, it is important to note that the viscosity parameter (𝑘) is reduced and has a constant value for the material response in the MD. However, it is dependent on the strain rate for the TD response. This is mainly due to the difference in the stress-strain responses in the MD and TD. From the TD stress-strain curves, the rates at which softening and hardening occur after the initial yield stress has been exceeded cannot be constant for every strain rate. Keeping this value constant will also cause the fitted curves to diverge from the experimental data. However, from the MD curves, this value can be kept constant as evident in Fig.5.7. The strain rate dependent parameters are plotted vs the logarithm of strain rate and fitted using the least-squares regression method to generate simple linear and polynomial relationships that can be introduced into the constitutive model to span out the values of the parameters within the range of strain rates. These fittings as well as their mathematical expressions are presented in Figs 5.8 and 5.9 for the MD and TD parameters respectively. 132 (a) (b) Figure 5.8 Least-squares regression fittings for viscoplastic parameters (a) initial yield stress and (b) strain rate scaling factor for MD response. (a) (b) (c) Figure 5.9 Least-squares regression fittings for viscoplastic parameters (a) initial yield stress (b) strain rate scaling factor, and (c) viscosity coefficient for TD response. 133 5.4 Introducing Temperature Effect To predict the temperature-dependent response of the investigated polymeric material under large deformations, a new term ℎ(𝑇), modified from the G’sell-Jonas model [1] is introduced into Eqn.5.2 such that: 1 1 ℎ(𝑇) = 𝑒𝑥𝑝 [𝑎 ∙ ( − )] (5.3) 𝑇 𝑇0 0 𝜎(𝜀𝑝 , 𝜀̇, 𝑇) = {𝜎𝑦0(𝜀̇) + {𝜇𝑣𝑝 ⋅ {1 − 𝑒𝑥𝑝(−𝑘(𝜀̇) ∙ 𝜀𝑝 )} + 𝐻10 ∙ 𝜀𝑝 + 𝐻20 ∙ 𝜀𝑝𝑚 } ⋅ 𝑎𝜀̇ (𝜀̇) } ∙ ℎ(𝑇) (5.4) where 𝑎 is a material parameter that accounts for the effect of temperature on the flow stress of the material. 𝑇 and 𝑇0 are the current and reference temperatures respectively. The modified term in Eqn.5.3 is different from the original term used in the G’sell-Jonas model [1] in the sense that it introduces the reference temperature to give a clearer meaning to the ℎ(𝑇) term. At the reference temperature ℎ(𝑇) becomes equal to 1 and Eqn.5.4 reduces to Eqn.5.3 which accounts for only the strain rate effect at the reference temperature. To implement the fully developed model, the values of the parameter 𝑎, have to be determined for the material response in the MD and TD. At the same strain value (𝜀), the value of the flow stress (𝜎) is related to the strain rate and current temperature (𝑇). To determine parameter 𝑎, at a constant large strain value and strain rate, some points are selected as: (𝜀, 𝜎1 ), (𝜀, 𝜎2 ), 𝑎𝑛𝑑 (𝜀, 𝜎𝑖 ) at different temperatures 𝑇1 , 𝑇2 , and 𝑇𝑖 respectively. Where 𝑖 is the maximum number of constant temperature values considered. At the reference temperature, the reference flow stress is obtained by Eqn.5.2. Therefore, the stress at the current temperature can be expressed as: 𝜎(𝜀𝑝 , 𝜀̇, 𝑇) = 𝜎(𝜀𝑝 , 𝜀̇, 𝑇0 ) ∙ ℎ(𝑇) (5.5) In this case, the reference temperature was selected as 20°C (293K) and reference stress is already known. Hence, parameter 𝑎 can be determined by fitting the relationship between the flow stress and temperature with Eqn.5.5. A similar methodology for determining the parameter 𝑎 was applied by Zhu et al [5] in constitutive modeling of thermoplastics. Furthermore, least square curve fittings were carried out according to Eqn.5.5 at different strain rates in the MD and TD as shown in Figs.5.10 through 5.13. The analysis for the determination of parameter 𝑎 was carried out at different large strain values. For MD samples, the analysis was carried out at 20% constant strain in the MD and the determined parameter 𝑎 values were different at different strain rates. The 134 determined values of parameter 𝑎 showed a linear relationship with the strain rate as shown in Fig.5.14. However, for the TD samples, the analysis was carried out at 20%, 30% and 40% constant strain values and the values of parameter 𝑎, determined were similar for the different strain values and the different strain rates. Hence, the average value of 1100.64 was selected as parameter 𝑎 to calibrate the material response in the TD. The determined values of parameter 𝑎 for the material response in the MD and TD are summarized in Table 5.5. (a) (b) (c) (d) Figure 5.10 Least-squares regression analysis at 20% strain for the determination of parameter 𝑎 for MD response at (a) 0.0001/s, (b) 0.001/s, (c) 0.01/s, and (d) 0.1/s strain rates. 135 (a) (b) (c) (d) Figure 5.11 Least-squares regression analysis at 20% strain for the determination of parameter 𝑎 for TD response at (a) 0.0001/s, (b) 0.001/s, (c) 0.01/s, and (d) 0.1/s strain rates. (a) (b) Figure 5.12 Least-squares regression analysis at 30% strain for the determination of parameter 𝑎 for TD response at (a) 0.0001/s, (b) 0.001/s, (c) 0.01/s, and (d) 0.1/s strain rates. 136 Figure 5.12 (Cont’d) (c) (d) (a) (b) (c) (d) Figure 5.13 Least-squares regression analysis at 40% strain for the determination of parameter 𝑎 for TD response at (a) 0.0001/s, (b) 0.001/s, (c) 0.01/s, and (d) 0.1/s strain rates. 137 Figure 5.14 Linear least-squares regression fitting of parameter 𝑎 for MD response. Table 5.5 Determined parameter 𝑎 values Strain rate 𝑴𝑫 𝟐𝟎% 𝑻𝑫 𝟐𝟎% 𝑻𝑫 𝟑𝟎% 𝑻𝑫 𝟒𝟎% TD Average (𝒔−𝟏 ) 0.0001 1581.70 1298.60 1293.60 1305.6 1299.27 0.001 1156.40 918.87 910.48 908.12 912.49 0.01 763.02 1144.5 1144.50 1154.30 1149.47 0.1 632.76 965.24 1063.4 1095.4 1041.35 5.5 Coupled Viscoelastic – Viscoplastic Model Implementation and Validation The developed coupled viscoelastic – viscoplastic model was implemented as a user- defined material model for shell element simulations in LS-DYNA®. To couple the orthotropic nonlinear viscoelastic model with the viscoplastic model, a uniaxial Von-Misses yield criterion was introduced to mark the onset of the yielding mechanism. The yield function introduced is expressed as: 𝑓 = 𝜎 𝑇𝑟𝑖𝑎𝑙 − 𝜎𝑦0 (5.6) where 𝜎 𝑇𝑟𝑖𝑎𝑙 is the viscoelastic trial stress. The onset of yielding is marked when the yield function 𝑓 is equal to or greater than zero. Once the yield threshold is exceeded, the flow stress evolution becomes governed by the viscoplastic mechanism. The coupled viscoelastic – viscoplastic model is implemented according to the following steps: Step 1: Implement nonlinear viscoelastic predictor 𝑛 ∆𝜌 ∆𝜌 ∆𝜎 = ℎ∞ ∙ 𝐺∞ ∙ ∆𝜀 + ∑ {ℎ𝑘 ∙ ∆𝜀 ∙ 𝐺𝑖 ∙ 𝑒𝑥𝑝 [− ] + 𝑒𝑥𝑝 [− ] ∙ 𝜎 𝑖 (𝑡) − 𝜎 𝑖 (𝑡)} (5.7) 𝜏𝑖 𝜏𝑖 𝑖=1 138 𝜎 𝑇𝑟𝑖𝑎𝑙 = 𝜎 + ∆𝜎 (5.8) Step 2: Check for plastic flow – uniaxial Von-Misses criterion If the yield function is less than zero (𝑓 < 0), then the total stress is the same as the viscoelastic trial stress (𝜎 = 𝜎 𝑇𝑟𝑖𝑎𝑙 ) and there is no plastic strain evolution (𝜀𝑝 = 0). However, the yield threshold is exceeded when the yield function becomes greater than or equal to zero (𝑓 ≥ 0). Step 3: Implement the flow stress predictor This step is only activated when the yield function becomes greater than or equal to zero. Once this is the case, the total stress is determined as the flow stress given by Eqn.5.4 and also there is plastic strain accumulation in the form: 𝜀𝑝 = 𝜀𝑝 + ∆𝜀 (5.9) The model predictions were validated by comparing simulation results with uniaxial tensile test results carried out at different strain rates and temperatures for samples cut along the MD and TD. All simulations were carried out using a one-element model, with boundary conditions as described in Chapter 3 and the required temperature history was applied at all four nodes of the shell element. Figs.5.15 and 5.16 presents the comparison between the experimental and the simulation results for the uniaxial stress-strain response in the MD and TD respectively. The results show that the stress-strain curves predicted by the coupled viscoelastic - viscoplastic model are in good agreement with the experimental results and capture the trend in the rate and temperature dependence of the stress-strain response of the polymeric separator. 139 (a) (b) (c) (d) Figure 5.15 Comparison of the stress-strain response predicted by LS-DYNA® simulation using the coupled viscoelastic - viscoplastic material model with experimental data from uniaxial tensile tests carried out at different strain rates for samples cut along the MD at (a) 20ºC, (b) 30ºC, (c) 40ºC, and (d) 50ºC. 140 (a) (b) (c) (d) Figure 5.16 Comparison of the stress-strain response predicted by LS-DYNA® simulation using the coupled viscoelastic - viscoplastic material model with experimental data from uniaxial tensile tests carried out at different strain rates for samples cut along the TD at (a) 20ºC, (b) 30ºC, (c) 40ºC, and (d) 50ºC. 5.6 Summary This chapter presents the development of a coupled viscoelastic-viscoplastic model for predicting the thermomechanical response of polymeric battery separators. The theoretical framework for the viscoplastic model is given on the basis of a rheological model consisting of a sliding frictional element, viscoplastic dashpot and two springs representing the strain softening and hardening mechanisms. The formulation for the stress-strain response of the polymeric material after the yield threshold has been exceeded takes the form of a modified G’sell – Jonas model. 141 The strain rate and temperature-dependent model parameters were determined for a selected PP Celgard®2400 separator. A parameter reduction method was introduced to decrease the number of strain rate-dependent parameters needed to calibrate the model. The model was implemented in LS-DYNA® finite element (FE) package as a user material model. The developed model was validated against uniaxial tensile tests carried out at constant strain rates and temperatures. The results show that the model predictions of the material anisotropy, rate dependence, and temperature dependence of the separator in its range of deformation before failure agree well with the experimental data. Furthermore, the validation results also confirm the validity of the offset method used to determine the yield stresses in the MD and TD as it was sufficient to describe the stress-strain response of the separator material before failure. The developed model is not only limited to the prediction of uniaxial tensile response as it can be calibrated to predict uniaxial compression. It can also be extended to account for biaxial responses of the separator. The developed coupled viscoelastic – viscoplastic model also has applications in thermomechanical modeling of a wide range of polymers and biomaterials. 142 REFERENCES 1. C. G'sell, and J. J. Jonas. Determination of the plastic behaviour of solid polymers at constant true strain rate. Journal of materials science, 14(3), pp.583-591, 1979. 2. M. Nasraoui, P. Forquin, L. Siad, and A. Rusinek. Influence of strain rate, temperature and adiabatic heating on the mechanical behaviour of poly-methyl-methacrylate: experimental and modelling analyses. Materials & Design, 37, pp.500-509, 2012. 3. Y. Duan, A. Saigal, R. Greif, and M. A. Zimmerman. A uniform phenomenological constitutive model for glassy and semicrystalline polymers. Polymer Engineering & Science, 41(8), pp.1322-1328, 2001. 4. H. Wang, H. Zhou, Z. Huang, Y. Zhang, and X. Zhao. Constitutive modeling of polycarbonate over a wide range of strain rates and temperatures. Mechanics of Time-Dependent Materials, 21(1), pp.97-117, 2017. 5. H. Zhu, H. Ou, and A. Popov. A new phenomenological constitutive model for thermoplastics. Mechanics of Materials, 157, p.103817, 2021. 143 Chapter 6 Conclusion and Future Work 6.1 Conclusion In this study, thermomechanical models have been developed to predict the response of polymeric battery separators in thermal ramp scenarios. Experimental techniques have been employed to characterize the mechanical and thermal responses of a selected PP separator. The experimental results were also employed to calibrate and validate the developed models under different combined thermal and mechanical loadings. The developed models can be implemented in commercial finite element analysis (FEA) software for simulating the thermomechanical response of polymeric battery separators. The developed models can also be incorporated into multiscale and multiphysics models for safety and health assessment of LIBs as well as vehicle crash simulations. It is important to remark that the models developed in this study can also be applied for modeling the thermomechanical response of other viscoelastic and viscoplastic materials such as polymers, polymer films, polymer matrix composites, and biomaterials. The major conclusions and findings from this study are as follows: 6.1.1 Orthotropic Linear Thermoviscoelastic Modeling An orthotropic linear thermoviscoelastic model has been developed in this work to predict the thermomechanical response of polymeric battery separators in thermal ramp scenarios. The model is built upon a linear viscoelastic framework and the temperature effect is introduced through the TTSP. Stress relaxation tests were carried out at small constant strain values to identify the linear viscoelastic limit of the material response. It was observed that an increase in temperature and the presence of a solvent increased the strain limit of the linear viscoelastic response of the material. A new TTSSM was also proposed to account for the plasticization effect of electrolyte solutions and predict the response of the polymeric separator in electrolyte solutions using its response in air as a framework. The developed model was calibrated using stress relaxation tests performed for a PP separator, Celgard®2400, within the linear viscoelastic limit at temperatures ranging from 20C to 110C in air and at temperatures of 20C to 40C in DMC. The model was implemented as a user material model in LS-DYNA®. Analytical solutions for stress relaxation were developed to verify the model predictions and the model predictions were validated against experimental results from stress relaxation, iso-strain temperature ramp tests and uniaxial tensile tests. The model predictions were found to accurately replicate the material anisotropy, rate 144 dependence, temperature dependence and electrolyte effect on the material response within the linear viscoelastic domain of the material deformation. The results from the model validation in non-isothermal conditions show that the simulations without considering the thermal expansion/shrinkage behavior of the separator resulted in large errors. 6.1.2 Orthotropic Nonlinear Thermoviscoelastic Modeling To predict the thermomechanical response of polymeric battery separators in thermal ramp scenarios under large deformation before the onset of irrecoverable deformations, an orthotropic nonlinear thermoviscoelastic model has been developed. A discretization algorithm was employed to evaluate the nonlinear viscoelastic hereditary integral with a kernel of Prony series based on a generalized Maxwell model assumed to consist of nonlinear springs and dashpots. Temperature dependence was introduced into the model through the transient linear relaxation modulus by implementing the TTSP. In this model, the nonlinearity and temperature effect are decoupled as the nonlinear parameters are only dependent on the strain and not temperature. The strain-dependent model parameters were calibrated using stress relaxation tests carried out at different constant strain values ranging from within the linear viscoelastic strain limit (~0.2%) up to 5% strain at the reference temperature (20°C). The TTSP was implemented using stress relaxation tests carried out within the linear viscoelastic limit at temperatures ranging from 20°C to 60°C. The model was implemented in LS-DYNA® finite element (FE) package as a user material model. Analytical solutions for stress relaxation were formulated based on the viscoelastic theory to verify the implementation of the model. The developed model was validated against stress relaxation tests, creep tests, uniaxial tensile tests at constant temperatures, and stress relaxation at non-isothermal conditions accounting for the thermal expansion of the material. The results show that the model predictions of the material anisotropy, rate dependence, and temperature dependence of the separator under large deformations in the nonlinear viscoelastic region agree well with the experimental data. 6.1.3 Coupled Viscoelastic – Viscoplastic Modeling The developed orthotropic nonlinear thermoviscoelastic model has been coupled with a viscoplastic model to predict the thermomechanical response of polymeric battery separators under large strains. A theoretical framework for the viscoplastic model was presented on the basis of a rheological model consisting of a sliding frictional element, viscoplastic dashpot and two springs 145 representing the strain softening and hardening mechanisms. The rheological framework captures the deformation mechanisms of polymeric materials and these mechanisms were reflected in the stress-strain formulation. The strain rate and temperature-dependent model parameters were calibrated for the selected PP Celgard®2400 separator. A parameter reduction method was introduced and employed to reduce the number of parameters needed for the calibration of the strain rate–dependent model parameters. The temperature dependence of the flow stress at constant large strain values was exploited to determine the temperature-dependent parameter necessary to introduce temperature dependence into the viscoplastic model. The model was implemented in LS-DYNA® finite element (FE) package as a user material model. The developed model was validated against uniaxial tensile tests carried out at constant strain rates and temperatures. The results show that the model predictions of the material anisotropy, rate dependence, and temperature dependence of the separator in its range of deformation before failure agree well with the experimental data. The developed coupled viscoelastic–viscoplastic model is capable of predicting the thermomechanical response of polymeric battery separators in every stage of their deformation before the onset of failure. 6.2 Future Work The developed thermomechanical models have been implemented in LS-DYNA® as user- defined material models, verified and validated against experimental results produced by tests carried out under combined thermal and mechanical loadings. Although great strides have been made in this study, there are a few more efforts that need to be taken for the thermomechanical modeling of battery separators to be complete. The future work necessary to complete the thermomechanical modeling of battery separators is divided into two major groups. The first group deals with model extension and adaptations while the second group deals with model implementation. 6.2.1 Model Extension and Adaptation To complete the thermomechanical model development and to make sure the constitutive responses and the deformation patterns of separators in their working condition are taken into account, the following works are required: 146 Validation for biaxial loading cases: In literature [1-3], the biaxial punch test has been identified as a representative loading case for separators in most real-world mechanical abuse scenarios. The developed orthotropic linear thermoviscoelastic model has been validated against biaxial punch tests carried out at ambient conditions [3]. The predicted strain patterns agreed well with the experimental results. However, more work needs to be done by adopting this same approach to validate the orthotropic nonlinear thermoviscoelastic model against biaxial punch tests carried out at different punch rates and temperatures. Furthermore, the viscoplastic model has to be extended for predictions of the separator response in biaxial loading cases. Plasticization effect of electrolyte solution: As detailed in Chapter 3, the separator in LIBs is immersed in an electrolyte solution that acts as a plasticizer and reduces its mechanical response. This effect has been accounted for in orthotropic linear thermoviscoelastic modeling through the time-temperature-solvent superposition method. However, further research has to be carried out to account for the solvent effect in the nonlinear thermoviscoelastic and viscoplastic modeling of battery separators. Multiaxial failure criteria: Future work will also involve introducing a multiaxial failure criterion to the coupled viscoelastic – viscoplastic model to mark the onset of failure of the separator under large deformations in multiaxial loading scenarios. 6.2.2 Model Implementation The developed thermomechanical models have been implemented for a single-layer PP separator. The developed model is also applicable for modeling tri-layered dry-processed polymeric separators as well as wet-processed separators with more isotropic microstructures, and mechanical and thermal behaviors. Future work will involve calibrating the developed models for application to these different types of battery separators. Furthermore, the extended model will be incorporated into multiscale and multiphysics models for coupled thermal and mechanical analysis of LIBs in thermal ramp scenarios. Aside from integration into models for thermal ramp scenarios in vehicle crash analysis, the developed models can be introduced into multiphysics models for the safety and health assessment of LIBs. 147 REFERENCES 1. X. Zhang, E. Sahraei, and K. Wang. Li-ion battery separators, mechanical integrity and failure mechanisms leading to soft and hard internal shorts. Scientific reports, 6(1), pp.1-9, 2016. 2. S. Kalnaus, A. Kumar, Y. Wang, J. Li, S. Simunovic, J. A. Turner, and P. Gorney. Strain distribution and failure mode of polymer separators for Li-ion batteries under biaxial loading. Journal of Power Sources, 378, pp.139-145, 2018. 3. S. Yan, J. Deng, C. Bae, S. Kalnaus, and X. Xiao. Orthotropic viscoelastic modeling of polymeric battery separator. Journal of The Electrochemical Society, 167(9), p.090530, 2020. 148 APPENDIX A. SUPPLEMENTARY EXPERIMENTAL RESULTS (a) (b) (c) Figure A.1 Uniaxial stress strain curves at different displacement rates for samples cut along the MD direction at (a) 30C, (b) 40C, and (c) 60C. 149 (a) (b) (c) (d) Figure A.2 Uniaxial stress strain curves at different displacement rates for samples cut along the TD direction at (a) 30C, (b) 40C, (c) 50C and (d) 60C. Figure A.3 Uniaxial stress strain curves for samples cut along the off-axis 45 direction at 20C at different strain rates. 150 Figure A.4 Stress relaxation test in air at 0.3% constant strain level with step temperature history (30C- 40C) for samples cut along the TD direction. (a) (b) Figure A.5 Iso-stress (iso-force) temperature ramp tests at different stress (force) levels in the (a) MD and (b) TD. 151 APPENDIX B. MATLAB SCRIPT FOR DETERMINATION OF PRONY SERIES PARAMETERS # MATLAB Script for the determination of Prony series parameters from stress relaxation tests # Royal Chibuzor Ihuaenyi, Michigan State University 2023 # Date of last modification: 05/1/2022 # Code Notes # The code is split into two parts. # 1) Subroutine for the stress relaxation function # 2) Subroutine for least square curve fitting, parameter and plot generation # Subroutine for the stress relaxation function function j=relaxation(x,t) ! Input relaxation times T1=10^(-1);T2=10^(0); T3=10^(1); T4=10^(2);T5=10^(3); T6=10^(4); T7=10^(5); T8=10^(6); ! Execute stress relaxation function j=x(1)*x(1)+x(2)*x(2).*exp(-t/T1)+x(3)*x(3).*exp(- t/T2)+x(4)*x(4).*exp(-t/T3)+x(5)*x(5).*exp(- t/T4)+x(6)*x(6).*exp(-t/T5)+x(7)*x(7).*exp(- t/T6)+x(8)*x(8).*exp(-t/T7)+x(9)*x(9).*exp(-t/T8); # Subroutine for least square curve fitting, parameter and plot generation clc clear x0 = [1;1;1;1;1;1;1;1;1]; xdata=dlmread('xdata.txt'); ydata=dlmread('ydata.txt'); [xe,resnorm] = lsqcurvefit(@relaxation,x0,xdata,ydata); T1=10^(- 1);T2=10^(0);T3=10^(1);T4=10^(2);T5=10^(3);T6=10^(4);T7=10^(5);T 8=10^(6); j=xe(1)*xe(1)+xe(2)*xe(2).*exp(-xdata/T1)+xe(3)*xe(3).*exp(- xdata/T2)+xe(4)*xe(4).*exp(-xdata/T3)+xe(5)*xe(5).*exp(- xdata/T4)+xe(6)*xe(6).*exp(-xdata/T5)+xe(7)*xe(7).*exp(- xdata/T6)+xe(8)*xe(8).*exp(-xdata/T7)+xe(9)*xe(9).*exp(- xdata/T8); subplot(2,1,1); subplot(2,1,1); plot(xdata,j,'r+') hold plot(xdata,ydata,'b*') legend('Fitted Curve','Experimental Data'); subplot(2,1,2) 152 semilogx(xdata,j,'r+') hold semilogx(xdata,ydata,'b*') legend('Fitted Curve','Experimental Data'); 153 APPENDIX C. PYTHON SCRIPT FOR NON-NEGATIVE LEAST SQUARE CURVE FITTING # Python Script for non-negative least square curve fitting for the determination of strain dependent nonlinear viscoelastic parameters # Royal Chibuzor Ihuaenyi, Michigan State University 2023 # Date of last modification: 10/10/2022 #Import necessary modules and curvefit functions import numpy as np import math import matplotlib.pyplot as plt %matplotlib inline from scipy.optimize import curve_fit from lmfit import minimize, Parameters, Parameter, report_fit #Read stress relaxation data in form of a csv file xdata,ydata = np.loadtxt("data.csv", usecols = (0,1), unpack= True, delimiter=',') #Input Prony series parameters for TTSP master curve G=[] T=[] #Execute curve fitting, parameter identification and plotting bounds=[] i=0 while i<2: i=i+0.001 bounds.append(i) print(bounds) def func2(params, x, data): h1 = params['h1'].value hk = params['hk'].value model = G0*h1 + (hk*(G1*np.exp(-x/T1)+G2*np.exp(- x/T2)+G3*np.exp(-x/T3)+G4*np.exp(-x/T4)+G5*np.exp(- x/T5)+G6*np.exp(-x/T6)+G7*np.exp(-x/T7)+G8*np.exp(-x/T8))) return model - data #This is what is minimized #Create a set of Parameters params = Parameters() for i in bounds: 154 params.add('h1', value = 1, min = i, max = i+2) #value is the initial condition params.add('hk', value = 1, max = 2) result = minimize(func2, params, args =(xdata, ydata)) final = ydata + result.residual rss = (result.residual**2).sum() # same as result.chisqr tss = sum(np.power(ydata - np.mean(ydata), 2)) r_squared=1 - rss/tss if r_squared>=0.95: #value can be regulated depending on your requirement on accuracy print(result.params) print(f"R² = {1 - rss/tss:.3f}") try: import pylab pylab.plot(xdata, ydata, 'k+') pylab.plot(xdata, final, 'r') pylab.show() except: pass else: pass 155 APPENDIX D. PUBLICATIONS AND PRESENTATIONS Refereed journal publications related to this thesis [1] R.C. Ihuaenyi, J. Deng, C. Bae, and X. Xiao, A coupled nonlinear viscoelastic-viscoplastic model for the thermomechanical response of polymeric battery separators (In Preparation). [2] R.C. Ihuaenyi, J. Deng, C. Bae, and X. Xiao, An orthotropic nonlinear thermoviscoelastic model for polymeric battery separators. Journal of the Electrochemical Society. DOI: 10.1149/1945-7111/acb178, 2023. [3] R.C. Ihuaenyi, S. Yan, J. Deng, C. Bae, S. Iqbal, and X. Xiao, Orthotropic Thermo- Viscoelastic Model for Polymeric Battery Separators with Electrolyte Effect. Journal of The Electrochemical Society. DOI: 10.1149/1945-7111/ ac24b6, 2021. Refereed journal and conference publications outside the scope of this thesis [4] S. Iqbal, B. Li, K. Sonta, R.C. Ihuaenyi, and X. Xiao. Probabilistic finite element analysis of sheet molding compound composites with an extended strength distribution model. Finite Elements in Analysis and Design, 214, p.103865, 2023. [5] S. Iqbal, B. Li, K. Sonta, R.C Ihuaenyi, and X. Xiao. Simulations of sheet molding compound composite structures with a randomly distributed tensile strength, Proceedings of CAMX 2021, Dallas, TX, October 2021. Presentations [6] R.C. Ihuaenyi, J. Deng, C. Bae, and X. Xiao, A nonlinear thermoviscoelastic model for polymeric battery separators, 19th U.S National Congress on Theoretical and Applied Mechanics, Austin, TX, June 2022. [7] R.C. Ihuaenyi, J. Deng, C. Bae, and X. Xiao, nonlinear thermoviscoelastic modeling of polymeric battery separators, MSU Engineering Graduate Research Symposium, MI, April 2022. 156