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. KolarÌŒiÌ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