A MULTIPHYSICS MODEL FOR THE STRESS ANALYSIS OF THE
SEPARATOR IN A LITHIUM-ION BATTERY
By
Wei Wu

A DISSERTATION
Submitted to
Michigan State University
in partial fulﬁllment of the requirements
for the degree of
Mechanical Engineering - Doctor of Philosophy

2013

ABSTRACT
A MULTIPHYSICS MODEL FOR THE STRESS ANALYSIS OF THE SEPARATOR IN
A LITHIUM-ION BATTERY
By
Wei Wu
The separator is a critical component to the durability and safety of lithium-ion (Liion) batteries. It prevents the physical contact between the positive and negative electrodes
while enabling ionic transport. The integrity of the separator is vital to the performance
and reliability of a battery.
In a Li-ion battery, stresses arise as the results of mechanical loading and constraints,
intercalation-induced deformation in the active materials in the electrodes, and thermal
expansion mismatch between the battery components. Currently, there are no methods
available to evaluate the in − situ stresses of the separator.
This work appertains to the development of a modeling tool for the evaluation of stresses
in the separator in a Li-ion battery. Towards this end, a multiphysics model for a basic
Li-ion battery cell had been developed and implemented in the multiphysics ﬁnite element
package COMSOL® . The model considered charge/species transport through diﬀusion
and migration, charge balance, electrochemical reaction kinetics, heat generation and heat
transfer, as well as the stress and deformation due to mechanical loading, Li intercalation
and temperature variation in the battery.
The multiphysics model was based on a validated one-dimensional (1D) battery model developed by Newman’s group and its implementation in COMSOL® for a Lix C6 |LiPF6 |Liy Mn2 O4
cell. For stress analysis, a 1D “Thermal”sub-model and a 2D meso-scale representative volume element (RVE) “Stress”sub-model were subsequently added. These sub-models were

fully coupled with each other. The temperature dependence of the transport properties and
reaction rate parameters was also included. The separator was modeled with a linear viscoelastic model based on experimental data. Its temperature dependence was established
through the time-temperature superposition principle.
The model was used to investigate the stress in the separator with battery cycles. The
simulation results revealed that the stress in the separator varies in phase with the battery cycles. The eﬀects of Li interaction and temperature rise in the battery could not be
superimposed, requiring both to be considered concurrently. The stress state and magnitude depended upon the mechanical properties of separator, active particle size and packing.
Besides addressing separator stresses, the model was also used to evaluate the inﬂuence
of the active particle size and component thickness on the heat generation rate and battery performance. The simulations revealed the evolution of the rate limiting mechanisms
with charge/discharge rates and a set of complex relationships between the battery design
parameters and battery utilization.
The 1D battery model was limited by the assumption that the electrodes were constructed
with uniformly distributed spherical particles of equal size. To consider the inﬂuence of real
microstrucutural eﬀects, a 2D microstructure resolved model has also been developed. This
model was used to investigate the stresses in the active particles and conductive binder.
This work presents a feasible approach to understanding the relationships among diﬀerent
physical phenomena in a Li-ion battery, evaluating the stress state and magnitude in battery
components such as the separator, and battery design optimization.

ACKNOWLEDGMENTS

I would foremost like to thank my supervisor, Dr. Xinran Xiao. She guided me onto the
road of research. What I learned from her instruction is not only the speciﬁc scientiﬁc knowledge and skill, but also the methodology to do research on the scientiﬁc problems. Without
her guidance, it would have been impossible for me to ﬁnish my Ph.D study at Michigan
State University. Her integrity, work ethic, and standards of conduct set an excellent example
for me to follow. I sincerely thank everything she did for me.
I would also deeply thank people who kindly helps me throughout my research. First of
all, I would like to thank Shutian Yan, for her beautiful work on the separator creep test
under diﬀerent temperatures. I really appreciate that she is willing to admit me to use the
measured property data at my thesis. Secondly, I would like to thank Danghe Shi for his
time to discuss on problems during my project. His valuable suggestion often renews my
ideas and provide new crew to solve the problem. Finally, I want to give special thanks to
Dr. Xiaosong Huang for the technique support and helpful discussion.
Deep appreciation is to thesis committee professors, Professor Andre Benard, Professor
Alfred Loos and Prefessor Wei Lai, for their helpful suggestions and comments to my research.
Their interest and response helps me to ﬁnish this research.
I also want to express my appreciation to my friends at Michigan State University: Wei
Zhang, Marie Wang, Fang Hou, Xuefei Chen, Chris Cater and Yalla Abushawashi. Thank
them to give me good memory of my Ph.D life. I am also indebted to the my family that have
helped me at every step. I would like to thank my parents, for their persistently support.
Most importantly, I heartily appreciate my husband Jack Shen. Without his patience and
encouragement, my thesis would never have been accomplished.
iv

This research is supported by NSF CMMI 1030821 and the start-up funding to Dr. Xiao
from Michigan State University.

v

TABLE OF CONTENTS
LIST OF TABLES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

ix

LIST OF FIGURES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

xi

KEY TO SYMBOLS AND ABBREVIATIONS . . . . . . . . . . . . . . . . . . . . .

xx

Chapter 1 Introduction . . . . . . . . .
1.1 Lithium-Ion Battery . . . . . . . .
1.2 Separator . . . . . . . . . . . . . .
1.3 Problem Deﬁnition . . . . . . . .
1.4 Literature Review . . . . . . . . .
1.4.1 Mathematical Modeling of
1.4.2 Stress in a Li-Ion Battery
1.4.3 Thermal Analysis . . . . .
1.5 Objective of Research . . . . . . .
1.6 Scope of the Thesis . . . . . . . .

. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
Phenomenon
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .

1
1
4
5
7
7
8
11
13
16

. . . . . . . . . . . . .
. . . . . . . . . . . . .
. . . . . . . . . . . . .
. . . . . . . . . . . . .
. . . . . . . . . . . . .
the Electrochemical
. . . . . . . . . . . . .
. . . . . . . . . . . . .
. . . . . . . . . . . . .
. . . . . . . . . . . . .

Chapter 2 Multiphysics Phenomena . . . . .
2.1 Electrochemical Phenomenon . . . . . .
2.1.1 Ionic Transport . . . . . . . . . . .
2.1.2 Ionic Charge Balance . . . . . . .
2.1.3 Mass Transport . . . . . . . . . . .
2.1.4 Electronic Charge Balance . . . .
2.1.5 Reaction Kinetics . . . . . . . . .
2.2 Thermal Phenomenon . . . . . . . . . . .
2.2.1 Heat Generation: Global Model
2.2.2 Heat Generation: Local Model .
2.3 Mechanical Phenomenon . . . . . . . . .
2.4 Multiphysics Coupling . . . . . . . . . .

.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.

17
17
18
19
21
22
23
24
25
26
29
32

Chapter 3 Multiphysics Model . . . . . . . .
3.1 Model depiction . . . . . . . . . . . . .
3.1.1 1D Battery Model . . . . . . . .
3.1.1.1 “Battery”Sub-Model .
3.1.1.2 “Electrode”Sub-Model
3.1.1.3 Boundary Conditions .
3.1.1.4 Model Implementation

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

35
35
37
39
40
41
43

vi

.
.
.
.
.
.
.

3.1.2 “Thermal”Sub-Model . . . . . . . . . . . . . . . . . . . . .
3.1.3 “Stress”Sub-Model . . . . . . . . . . . . . . . . . . . . . .
3.2 Material Properties and Simulation Condition . . . . . . . . .
3.2.1 Temperature Dependent Variables . . . . . . . . . . . . .
3.2.1.1 Variables Related to the Electrodes . . . . . . .
3.2.1.2 Variables Related to the Li+ Transport in the
trolyte . . . . . . . . . . . . . . . . . . . . . . . . .
3.2.2 Viscoelastic Property of Separator . . . . . . . . . . . . .
3.2.2.1 Time Dependent Viscoelastic Behavior . . . . .
3.2.2.2 Temperature Eﬀect . . . . . . . . . . . . . . . . .
3.2.3 Battery Parameters and Other Properties . . . . . . . .
3.3 Simulation Condition . . . . . . . . . . . . . . . . . . . . . . . . .

. . . .
. . . .
. . . .
. . . .
. . . .
Elec. . . .
. . . .
. . . .
. . . .
. . . .
. . . .

43
44
46
46
48
49
53
54
56
63
68

Chapter 4 Heat Generation in a Li-Ion Battery . . . . . . . . . . . . .
4.1 Model Validation . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.2 Thermal Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.3 The Eﬀect of Temperature Dependence of Physical Properties
4.4 Heat Generation . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.5 Eﬀects of Battery Design Parameters . . . . . . . . . . . . . . . .
4.5.1 Particle Size . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.5.2 Component Thickness . . . . . . . . . . . . . . . . . . . . . .
4.5.3 Component Thickness and Particle Size Eﬀects . . . . . .
4.5.4 Tortuosity . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.6 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.

71
71
72
73
78
80
83
88
88
90
95

Chapter 5 Stress in the Separator . . . . . . . . . . . .
5.1 Numerical Experiments . . . . . . . . . . . . . . .
5.2 Numerical Investigations . . . . . . . . . . . . . .
5.2.1 The Eﬀect of Particle Packing . . . . . .
5.2.2 Concentration Change and Temperature
5.2.3 Thickness Variation . . . . . . . . . . . . .
5.2.3.1 Particle Indentation . . . . . . . .
5.2.3.2 Cell Thickness Variation . . . . .
5.2.4 Strain and Stress in the Separator . . . .
5.2.5 The Eﬀect of Discharge Rate . . . . . . .
5.2.6 Anisotropy . . . . . . . . . . . . . . . . . . .
5.3 Summary . . . . . . . . . . . . . . . . . . . . . . . .

. . .
. . .
. . .
. . .
Rise
. . .
. . .
. . .
. . .
. . .
. . .
. . .

Chapter 6 Model Extension . . . . . . . . . . . . . . . . .
6.1 The Needs of a Microstructural Resolved Model
6.2 Model Depiction . . . . . . . . . . . . . . . . . . . .
6.2.1 “Electrochemical”Sub-model . . . . . . . .
6.2.1.1 Active Electrode Particles . . . . .
6.2.1.2 Electrolyte . . . . . . . . . . . . . .
6.2.1.3 Separator . . . . . . . . . . . . . . .
vii

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.

97
97
100
100
102
104
104
109
114
120
125
126

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

128
128
131
131
131
134
136

6.3

6.4
6.5
6.6

6.7

6.2.1.4 Binder . . . . . . . . . . . . . . .
6.2.1.5 Current Collector . . . . . . . .
6.2.1.6 Model Implementation . . . . .
6.2.2 “Stress”Sub-model . . . . . . . . . . . . .
Material Properties and Battery Parameters
6.3.1 Conductive Binder . . . . . . . . . . . . .
6.3.2 Other Battery Components . . . . . . .
Numerical Experiment . . . . . . . . . . . . . .
Model Validation . . . . . . . . . . . . . . . . . .
Results and Discussion . . . . . . . . . . . . . .
6.6.1 Concentration Change . . . . . . . . . .
6.6.2 Stress in the Electrode . . . . . . . . . .
6.6.2.1 Stress in the Active Particles .
6.6.2.2 Stress in the Binder . . . . . . .
6.6.3 The eﬀect of Discharge Rate . . . . . .
Summary . . . . . . . . . . . . . . . . . . . . . . .

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

136
137
137
137
138
138
141
142
143
145
146
148
148
154
155
162

Chapter 7 Conclusion and Future Work . . . . . . . . . . . . . . . . . . . . . 166
7.1 Summary and Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . 166
7.2 Outlook of Future Work . . . . . . . . . . . . . . . . . . . . . . . . . . . 169
APPENDICES . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . 172

Appendix A Non-Linear Least-Square Subroutine . . . . . . . . . . . . . . . . 173
Appendix B Laplace Transformation Subroutine . . . . . . . . . . . . . . . . . 174
BIBLIOGRAPHY

. . . . . . . . . . . . .

viii

. . . . . . . . . . . . . . . . . 176

LIST OF TABLES

Table 3.1

Temperature dependency of six sets of physical properties. . . . . .

Table 3.2

Solid phase diﬀusion coeﬃcient D1 and reaction rate constant k0

.

48

Table 3.3

Prony series parameters for the stress relaxation modulus in 1.1M
LiPF6 in EC/DMC (1:1 by volume). . . . . . . . . . . . . . . . . .

57

Table 3.4

Critical battery parameters used in the battery model. . . . . . . . .

64

Table 3.5

Physical properties of the battery components in “Thermal”sub-model. 65

Table 3.6

Mechanical properties of diﬀerent battery components used in the
“Stress”sub-model. . . . . . . . . . . . . . . . . . . . . . . . . . . .

66

Table 3.7

Mechanical properties of electrolyte LiPF6 in EC/DMC (1:1 by volume) 67

Table 3.8

Element information in four sub-models of the multiphysics model. .

68

Table 4.1

Comparison of the predicted minimum and maximum heat generation
rates with the available experimental data in literature. . . . . . . .

82

Models with diﬀerent particle sizes R and component thicknesses L
investigated in the design parameter study . . . . . . . . . . . . . .

83

The description of exponent p, tortuosity τt and eﬀective electrolyte
transport properties for battery I and battery II respectively. . . . .

92

Parameters in the 1D “Battery”sub-model with the original battery
size and in the multi-scale model with a quarter-sized or a half-sized
1D RVE “Stress”sub-model. . . . . . . . . . . . . . . . . . . . . . .

98

Table 4.2

Table 4.3

Table 5.1

Table 5.2

The comparison of Li concentration range in “Electrode”and “Stress”submodels at the end of discharge (1000 s) and the end of charge (2500
s) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104

ix

Table 5.3

The volume fractions of active materials and liquid phase in Battery
1 and 2. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 110

Table 5.4

The range of Li concentration in sub-models at DOD = 0.3 at diﬀerent discharge rates. . . . . . . . . . . . . . . . . . . . . . . . . . . . 120

Table 6.1

Prony series parameters for the stress relaxation modulus of PVDFHFP (FC1278). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 140

Table 6.2

Critical battery parameters used in the half-sized 2D microstructure
resolved model. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 143

x

LIST OF FIGURES

Figure 1.1

Figure 1.2

Figure 2.1

Figure 3.1

Figure 3.2

Schematic of a Li-ion battery being charged [1]. The lower pictures
in the ﬁgure are scanning electron microscope (SEM) images of different battery component materials. From left to right, they are a
current collector, which is a metal foil to collect current from outside
circuit, a negative electrode using graphite as active particle materials [1], a polymeric separator (Celgard) [2], the positive particles
using Liy Mn2 O4 active particles [3]. For interpretation of the references to color in this and all other ﬁgures, the reader is referred to
the electronic version of this dissertation. . . . . . . . . . . . . . . .

2

The fully coupled multiphysics model consists of four sub-models.
The schematic shows the physical phenomena in each sub-model and
the coupling relationships among four sub-models. . . . . . . . . .

15

The coupling relationships among thermal, electrochemical and mechanical phenomena. . . . . . . . . . . . . . . . . . . . . . . . . . .

33

The meso-scale representative volume element (RVE) represents the
stress state at the center of the battery cell. . . . . . . . . . . . . . .

37

COMSOL® implementation of the fully coupled thermal-electrochemicalstress model. The electrochemical performance of the battery is evaluated with an existing 1D Li-ion battery model in COMSOL® which
consists of two sub-models; (a) “Battery”and (b) “Electrode”. The
thermal analysis is performed using (c) “Thermal”sub-model. The
stress analysis for the components in a cell is performed concurrently
using (d) a meso-scale RVE “Stress”sub-model. . . . . . . . . . . . . 38

Figure 3.3

The boundary conditions for stress analysis.

Figure 3.4

Temperature dependence of (a) the diﬀusion coeﬃcient D1 and (b)
reaction rate constant k0 for Li in the negative (graphite) and positive
(Liy Mn2 O4 ) electrodes. The values of D1 at 25 ◦ C are 3.9 × 10−14
m2 s−1 for graphite and 1 × 10−13 m2 s−1 for Liy Mn2 O4 . The
value of k0 at 25 ◦ C is 2 × 10−6 A m2.5 mol−1.5 for both electrodes. 49

xi

. . . . . . . . . . . . .

45

The entropy change as a function of SOC at 25 ◦ C for (a) graphite
[104] and (b) Liy Mn2 O4 [105]. . . . . . . . . . . . . . . . . . . . . .

50

OCP of (a) graphite and (b) Liy Mn2 O4 vs. lithium metal as a function of SOC and temperature. The curves at the reference temperature 25 ◦ C were determined by experiments [32, 104]. . . . . . . . .

50

The ionic conductivity k2 as a function of the salt concentration and
temperature in 1M LiPF6 in EC/DMC (2:1 by volume), calculated
with Equation 3.16. . . . . . . . . . . . . . . . . . . . . . . . . . . .

52

The Li+ diﬀusion coeﬃcient D2 as a function of the salt concentration and temperature in LiPF6 in EC/DMC (2:1 by volume), calculated with Equation 3.14. . . . . . . . . . . . . . . . . . . . . . . .

53

The product of the thermodynamic factor and the anion transference
number v as a function of the salt concentration and temperature in
LiPF6 in EC/DMC (2:1 by volume), calculated with Equation 3.15.

54

SEM image of the surface of a single layer Celgard 2400 separator [9].
The arrow indicates the transverse direction (TD). . . . . . . . . . .

55

The stress relaxation function for the separator Celgard 2400 in 1.1
M LiPF6 in EC/DMC (1:1 by volume) [110]. . . . . . . . . . . . . .

57

The images depicting the sample installation and isothermal test establishment. (a) sample installation; (b) & (c) the establishment of
the isothermal test by closing the furnace. . . . . . . . . . . . . . . .

59

Figure 3.13

The creep compliance of Celgard 2400 at diﬀerent stress levels. . . .

60

Figure 3.14

The creep compliance
◦ C, 40 ◦ C, 45 ◦ C, 50
MPa. . . . . . . . . .

of Celgard 2400 at diﬀerent temperatures: 25
◦ C, 55 ◦ C and 60 ◦ C. The stress level is 1.5
. . . . . . . . . . . . . . . . . . . . . . . . . .

61

Construction of a master curve for PP using shear stress relaxation
experimental data. The stress plots are at 25 ◦ C, 40 ◦ C, 45 ◦ C, 50 ◦ C,
55 ◦ C and 60 ◦ C. The above stress relaxation plots are obtained from
creep compliance measured curves at the corresponding temperatures
by Laplace transformation. . . . . . . . . . . . . . . . . . . . . . . .

61

The relationships of aT and temperature for PP. . . . . . . . . . . .

62

Figure 3.5

Figure 3.6

Figure 3.7

Figure 3.8

Figure 3.9

Figure 3.10

Figure 3.11

Figure 3.12

Figure 3.15

Figure 3.16

xii

Figure 3.17

Figure 4.1

Figure 4.2

Figure 4.3

Figure 4.4

Figure 4.5

Figure 4.6

Figure 4.7

Figure 4.8

Finite element meshes for (a) the “Battery”; (b) the “Electrode”; (c)
the “Thermal”and (d) RVE “Stress”sub-models. . . . . . . . . . . .

69

Comparison of predicted temperature rises with the experimental results of Ye et al. [41]. . . . . . . . . . . . . . . . . . . . . . . . . . .

73

Inﬂuence of the temperature dependent parameters on the predicted
temperature rises during a galvanostatic discharge at 1C (a-c), 4C (df) and 8C (g-i). Simulations were performed with one temperature
dependent parameter at a time. The corresponding temperature rises
from coupled and decoupled model were also plotted in each ﬁgure. .

75

(a) The temperature rise and (b) cell utilization as a function of
discharge rate. The cut oﬀ voltage is 3.0 V. . . . . . . . . . . . . .

78

Partial heats due to diﬀerent sources contributing to the total heat
during 1C discharge. . . . . . . . . . . . . . . . . . . . . . . . . . . .

79

Heat generation rate versus DOD for discharge rates from 0.1C to
6C: (a) 0.1C-2C; (b) 0.1C-6C. . . . . . . . . . . . . . . . . . . . . .

81

Comparison of the predicted heat generation rates with the upper and
lower bounds of experimental data in literature (a) the maximum heat
generation rate; (b) the minimum heat generation rate. . . . . . . .

82

(a) The cell temperature rise and (b) utilization at the end of the
discharge in the 0.1R/1L, 1R/1L and 1.5R/1L models. . . . . . . . .

84

Comparison of the Li+ concentration proﬁles in the 0.1R/1L, 1R/1L
and 1.5R/1L models under three discharge rates. (a, b) at DOD=0.91,
1C; (c, d) at DOD=0.55, 3C; (e, f) 5C; when one of the cells reached
the cut-oﬀ voltage. (a, c, e) Li+ concentration in the electrolyte, (b,
d, f) Li concentration at the surface of the active particles. . . . . .

86

Figure 4.9

The evolution of Li+ concentration proﬁles in electrolyte in the 0.1R/1L
and 1R/1L models at 1C discharge rate; (a) at DOD=0.026 (100 s),
(b) DOD=0.28 (1000 s), (c) DOD=0.53 (2000 s), (d) DOD=0.91
(3400 s). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87

Figure 4.10

The temperature rise and cell utilization at the end of the discharge.
(a) temperature rise, 1R/0.5L and 1R/1L cells; (b) temperature rise,
0.1R/0.5L and 0.1R/1L cells; (c) cell utilization, 1R/0.5L and 1R/1L
cells; (d) cell utilization, 0.1R/0.5L and 0.1R/1L cells. . . . . . . . .
xiii

89

Figure 4.11

Figure 4.12

Figure 4.13

Figure 4.14

Figure 4.15

Figure 5.1

Figure 5.2

The temperature rise and cell utilization at the end of the discharge
of 0.1R/0.5L and 1R/0.5L. (a) temperature rise; (b) cell utilization.

90

Comparison of the discharge proﬁles between cells with a thickness
of 0.5L and 1L. Discharge rate: (a) 1C, (b) 4C, and (c) 10C. . . . .

91

The temperature rise and cell utilization at the end of the discharge
from the battery models with diﬀerent particle size: (a) temperature
rise; (b) cell utilization. The models are simulated with p values 3.3
and 2, respectively. . . . . . . . . . . . . . . . . . . . . . . . . . . .

92

The temperature rise and cell utilization from the battery models with
two battery thicknesses: (a) temperature rise; (b) cell utilization. The
models are simulated with a Bruggeman exponent of p=3.3 and p=2.

93

Comparison of the Li+ concentration proﬁles of the 1R/0.5L and
1R/1L cells under a discharge rate of 3C for the p=2 battery and the
p=3.3 battery respectively. (a) Li+ concentration in the electrolyte
for the p=3.3 battery; (b) Li concentration at the surface of the active particles for the p=3.3 battery; (c) Li+ concentration in the
electrolyte for the p=2 battery; (d) Li concentration at the surface
of the active particles for the p=2 battery. (a) and (b) were made at
DOD=0.706, when 1R/1L ﬁrstly reaches the cut-oﬀ voltage 3V; and
(c) and (d) were made at DOD=0.908, when 1R/0.5L ﬁrstly reaches
the cut-oﬀ voltage 3.0 V. . . . . . . . . . . . . . . . . . . . . . . . .

94

(a) A battery cycle for a 1R/0.5L unit cell at a 2C discharge and
charge rate. The cycle consists of discharge (0 s-1000 s), open circuit period (1000 s-1500 s), charge(1500 s-2500 s), and open circuit
period(2500 s-3000 s). (b) The cell voltage history for the dischargecharge cycle. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

99

The predicted through thickness strain distribution at the end of
discharge (1000 s) under a ﬁxed-free (10 psi) boundary condition
using a quarter-sized 2D “Stress”sub-model. The areas with a void
like appearance have a strain value beyond the range and hence are
removed from the display. Particles 2 and 4 have a close-packed
pattern whereas particles 1, 3 and 5 have a loose-packed pattern.
The maximum strain in the separator is found near a loose-packed
particle of larger diameter. . . . . . . . . . . . . . . . . . . . . . . . 100

xiv

Figure 5.3

The Li concentration distribution in the negative and positive electrode at diﬀerent times during discharge-charge cycle: (a) 0 s, (b)
500 s, (c) 1000 s, (d) 1500 s, (e) 2000 s, and (f) 2500 s. . . . . . . . 103

Figure 5.4

The temperature rise history in the Li-ion battery during dischargecharge cycle: discharge (0 s-1000 s), open circuit period (1000 s-1500
s), charge(1500 s-2500 s), and open circuit period(2500 s-3000 s). . . 105

Figure 5.5

The through thickness displacement ﬁeld at t = 0 s predicted by
a half-sized 2D RVE “Stress”sub-model with close-packed particles.
The black lines in the separator represent the position where the
through thickness direction at the end of discharge was checked in
Figure 5.6. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106

Figure 5.6

The indentations into the separators made by negative and positive
electrode particles. (a) caused by both Li intercalation and thermal
expansion; (b) caused by Li intercalation; (c) caused by thermal expansion. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 108

Figure 5.7

The displacement-history plots at the right hand end of the cell. Battery 1 (ε1,pos = 0.362) and Battery 2 (ε1,pos = 0.301) were investigated. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 110

Figure 5.8

The displacement-history plots at the right hand end of the cell. (a)
the thickness caused by thermal expansion (b) the thickness caused
by Li intercalation. . . . . . . . . . . . . . . . . . . . . . . . . . . . 111

Figure 5.9

The through thickness displacement and corresponding Li particle
surface concentration at diﬀerent discharge times: (a) through thickness displacement, 360 s; (b) Li concentration, 360 s; (c) through
thickness displacement, 720 s; (d) Li concentration, 720 s; (e) through
thickness displacement, 1000 s; (f) Li concentration, 1000 s. The data
was from the simulation that intercalation stress is considered alone. 112

Figure 5.10

The through thickness stress and strain across the RVE ”Stress” submodel at the end of discharge. . . . . . . . . . . . . . . . . . . . . . 115

Figure 5.11

Stress components in the separator in the cell at the end of discharge.
(a) σx ; (b) σy ; (c) σxy . . . . . . . . . . . . . . . . . . . . . . . . . . 116

xv

Figure 5.12

The history plots of the principal stresses and the Von Mises stress
in the separator near the negative and positive electrode particles.
(a) (c) (e) and (g): near negative electrode; (b) (d) (f) and (h): near
positive electrode. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 118

Figure 5.13

The history plots of the principal stresses and the Von Mises stress
in the separator near the negative and positive electrode particles.
(a) 1st principal stress; (b) 2nd principal stress; (c) 3rd principal
stress; (d) Von Mises stress. Three diﬀerent mechanical properties of
separator are considered: constant linear elasticity (E(const)); time
dependent viscoelasticity (E(t)); and time and temperature dependent viscoelsaticity (E(t,T)). . . . . . . . . . . . . . . . . . . . . . . 121

Figure 5.14

The temperature rise Vs. DOD under discharge rate 2C and 4C. The
DOD range is 0 to 0.7. . . . . . . . . . . . . . . . . . . . . . . . . . 123

Figure 5.15

The comparison of the thickness variations for the half-sized cell under
discharge rate 2C and 4C. The DOD range is 0 to 0.7. . . . . . . . . 123

Figure 5.16

The evolution of the stress components in the separator near the
negative particles under diﬀerent discharge rate: 2C and 4C: (a) the
ﬁrst principal stress; (b) second principal stress; (a) third principal
stress and (d): Von Mises stress. The DOD range is from 0 to 0.7. . 124

Figure 5.17

The cell voltage for the half-sized cell under discharge rate 2C and
4C. The DOD range is 0 to 0.7. . . . . . . . . . . . . . . . . . . . . 125

Figure 6.1

(a) The cross-section image of a typical LiC6 -LiCoO2 rechargeable
Li-ion battery [37]. The upper, middle and bottom layer correspond
to the graphite negative electrode, polymeric separator and LiCoO2
positive electrode, respectively. (b) An individual slice from a 3D
X-ray tomography of graphite electrode [138]. . . . . . . . . . . . . 130

Figure 6.2

The schematic of the 2D microstructure resolved model. . . . . . . . 132

Figure 6.3

A schematic diagram of part of the Li-ion battery (one current collector, one electrode and the separator) with the governing equations. 132

Figure 6.4

The boundary equations for the Electrochemical sub-model. . . . . . 133

Figure 6.5

The boundary conditions for stress analysis.

xvi

. . . . . . . . . . . . . 138

Figure 6.6

The strain-stress curves for P(VDF-HFP) materials [150–152]. . . . 139

Figure 6.7

The comparison of the cyclic stress strain curves from simulation and
experiment for PVDF-HFP (FC1278) [152]. . . . . . . . . . . . . . . 141

Figure 6.8

The stress strain curve of the PVDF-HFP−LiClO4 [150] . . . . . . . 142

Figure 6.9

The mesh for the geometry of the 2D microstructure resolved model.
Sections (a) and (b) correspond to the zoomed-in regions. (a): mesh
in the negative electrode; (b): mesh in the positive electrode. . . . . 144

Figure 6.10

The comparison of the predicted cell voltage with the experimental
results [32]: (a) 0.5C; (b) 1C. . . . . . . . . . . . . . . . . . . . . . . 145

Figure 6.11

The Li concentration distribution in the negative and positive electrode at diﬀerent times during discharge-charge cycle: (a) 0 s; (b)
500 s; (c) 1000 s; (d) 1500 s; (e) 2000 s; and (f) 2500 s. The unit for
concentration in the graph is molm−3 . . . . . . . . . . . . . . . . . 147

Figure 6.12

The ﬁrst principal stress distribution in the
diﬀerent times during discharge-charge cycle:
1000 s; (d) 1500 s; (e) 2000 s; and (f) 2500 s.
the graph is Pa. . . . . . . . . . . . . . . . .

Figure 6.13

The ﬁrst principal stress distribution in the positive electrode at different times during discharge-charge cycle: (a) 0 s; (b) 500 s; (c) 1000
s; (d) 1500 s; (e) 2000 s; and (f) 2500 s. The unit for stress in the
graph is Pa. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 150

Figure 6.14

The evolution of the stress components in the active particles. The
stress curves were plotted from diﬀerent locations: (a) negative electrode particle center, location (2) in Figure 6.12; (b) negative electrode particle boundary, location (1) in Figure 6.12; (c) positive electrode particle center, location (3) in Figure 6.13; (d): positive electrode particle boundary, location (4) in Figure 6.13. . . . . . . . . . 152

Figure 6.15

The evolution of the stress components at the normal point of particle
boundary. The stress curves were plotted from diﬀerent locations:
(a) negative electrode particle boundary with small radius, location
(5) in Figure 6.12; (b) negative electrode particle boundary with large
radius, location (6) in Figure 6.12; (c) positive electrode particle with
small radius, location (7) in Figure 6.13; (d) positive electrode particle
boundary with large radius, location (8) in Figure 6.13. . . . . . . . 153
xvii

negative electrode at
(a) 0 s; (b) 500 s; (c)
The unit for stress in
. . . . . . . . . . . . . 149

Figure 6.16

The contour of the ﬁrst principal stress in the binder corresponding
to diﬀerent binder materials at the end of discharge (1000 s): (a)
PVDF-HFP(FC1278), E=1.1 MPa; (b) PVDF-HFP-LiClO4 , E=184
MPa. The unit for stress in the graph is Pa. . . . . . . . . . . . . . 154

Figure 6.17

Stress history in binder: (a) binder with low modulus at a location
in the negative electrode shown in Figure 6.16 (a); (b) binder with
low modulus at a location in the positive electrode in Figure 6.16 (a);
(c) binder with high modulus at a location in the negative electrode
in Figure 6.16 (b); (d) binder with high modulus at a location in the
positive electrode in Figure 6.16 (b). . . . . . . . . . . . . . . . . . . 156

Figure 6.18

Strain history in binder: (a) binder with low modulus at a location
in the negative electrode shown in Figure 6.16 (a); (b) binder with
low modulus at a location in the positive electrode in Figure 6.16 (a);
(c) binder with high modulus at a location in the negative electrode
in Figure 6.16 (b); (d) binder with high modulus at a location in the
positive electrode in Figure 6.16 (b). . . . . . . . . . . . . . . . . . . 157

Figure 6.19

The evolution of the ﬁrst principal stress in the active particles under diﬀerent discharge rate: (a) negative electrode particle center,
location (2) in Figure 6.12; (b) negative electrode particle boundary,
location (1) in Figure 6.12; (c) positive electrode particle center, location (3) in Figure 6.13; (d) positive electrode particle boundary,
location (4) in Figure 6.13. . . . . . . . . . . . . . . . . . . . . . . . 159

Figure 6.20

The evolution of the Von Mises stress in the active particles under different discharge rate: (a) negative electrode particle center, location
(2) in Figure 6.12; (b) negative electrode particle boundary, location
(1) in Figure 6.12; (c) positive electrode particle center, location (3)
in Figure 6.13; (d) positive electrode particle boundary, location (4)
in Figure 6.13. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 160

Figure 6.21

The evolution of the ﬁrst principal stress in the polymeric binder under diﬀerent discharge rate. The stress curves were plotted from high
stress location with diﬀerent binder material: (a) PVDF-HFP(FC1278)
(low modulus), location (1) in Figure 6.16(a); (b) PVDF-HFP(FC1278)
(low modulus), location (2) in Figure 6.16(a); (c) PVDF-HFP-LiClO4
(high modulus), location (3) in Figure 6.16(b); (d) PVDF-HFP-LiClO4
(high modulus), location (4) in Figure 6.16(b). . . . . . . . . . . . . 163

xviii

Figure 6.22

The evolution of the Von Mises stress in the polymeric binder under diﬀerent discharge rate. The stress curves were plotted from high
stress location with diﬀerent binder material: (a) PVDF-HFP(FC1278)
(low modulus), location (1) in Figure 6.16(a); (b) PVDF-HFP(FC1278)
(low modulus), location (2) in Figure 6.16(a); (c) PVDF-HFP-LiClO4
(high modulus), location (3) in Figure 6.16(b); (d) PVDF-HFP-LiClO4
(high modulus), location (4) in Figure 6.16(b). . . . . . . . . . . . . 164

xix

KEY TO SYMBOLS AND ABBREVIATIONS

a

ion number

Acell

cross section area of the battery (cm2 )

c

Li concentration (mol m−3 )

Cp

speciﬁc heat capacity (J (kg K)−1 )

D1

diﬀusion coeﬃcient of Li+ in the solid electrode particles (m2 s− 1)

D2

electrolyte salt diﬀusivity (m2 s−1 )

f±

average molar activity coeﬃcient

E

Young’s modulus (MPa)

Ea,D

the activation energies associated with D1 (KJ mol−1 )

Ea,R

the activation energies associated with k0 (KJ mol−1 )

F

Faradays constant, 96487 (C mol−1 )

G

stress relaxation modulus (Pa)

h

lumped heat transfer coeﬃcient (W m−2 K−1 )

¯
Hij

partial molar enthalpy (J mol−1 )

∆Hp

enthalpy change in the chemical reaction p (J mol−1 )

i

current density (A m−2 )

j0

exchange current density (A m−2 )

jn

local current density (A m−2 )

J0

instantaneous compliance (N m−2 )

J

compliance (N m−2 )

k0

reaction rate constant (A m2.5 mol−1.5 )

k1

solid phase electronic conductivity (S m−1 )
xx

k2

liquid phase ionic electrocity (S m−1 )

Ktherm

thermal conductivity (W (m K)−1 )

li,j

constants

L

thickness of each battery component i (µm)

n

charge number

M

mobility of solute (m2 (V s)−1 )

N0

Li+ ﬂux (mol m−2 s−1 )

p

the exponent in empirical relation for eﬀective transport properties

pc

constant

q

heat generation rate (W L−1 )

r

radial distance from the center of electrode active particle (µm)

R

electrode active particle radius (µm)

Rc

gas constant, 8.314 (J (mol K)−1 )

Rcol

electronic resistivity of current collector (Ω m)

Sa

speciﬁc surface area (m−1 )

St

maximum tensile stress (Pa)

t

time (s)

t+

transport number of Li+

t−

transport number of anion

T

absolute temperature (K)

U

open circuit potential of the electrode (V)

v

the product of thermodynamic factor and the anion transference number

v0

solvent velocity (m s−1 )

xxi

V

cell potential (V)

z

charge number

Greek letters
α

reaction coeﬃcient

αij

thermal expansion coeﬃcient (K−1 )

αT

horizontal time-temperature shift factor

γp

reaction rate of a side chemical reaction p (mol m−1 s−1 )

Γi

active particle/electrolyte interface

Γn

negative electrode/currrent collector interface

Γp

positive electrode/currrent collector interface

δij

Dirac delta function

∆Si

entropy change (J (mol K)−1 )

ε

porosity of the battery electrode

εij

strain component

η

local surface overpotential (V)

µ

chemical potential (V)

ν

Poissons ratio

Π

interface heat generation due to the entropy change in the electrode (W L−1 )

ρ

material density (kg m−3 )

σ

stress (Pa)

σb

battery domain

σe

electrolyte domain

σh

hydrostatic component of stress (Pa)
xxii

σs

solid, active particle domain

σY

yield stress (Pa)

τ

retardation time (s)

τt

tortuosity

φ

electrical potential (V)

Subscripts and superscripts
1

solid phase or ﬁrst principal

2

liquid phase or second principal

3

third principal

c

current collector

i

cell component or constant

n

negative electrode or reaction charge number

p

positive electrode

s

separator

ave

average value

ef f

eﬀective value

ei, T

eigen strain due to thermal expansion

ei, c

eigen strain due to intercalation

ini

initial value

max

maximum value

me

mechanical

ref

reference value

surf

surface value
xxiii

VM

Von Mises

∞

inﬁnity

+

positive

−

negative

xxiv

Chapter 1
Introduction
After the ﬁrst commercial product of lithium-ion(Li-ion) batteries by Sony Corporation in
1991, the market has grown rapidly. Currently, Li-ion batteries are widely used in consumer
electronics such as laptops, digital cameras, and cell phones. The applications for Li-ion
batteries have recently expanded to include electric or hybrid-electric vehicles by automakers
as well as aerospace missions conducted by NASA [1] . Compared to other battery systems,
the Li-ion battery oﬀers numerous advantages such as light weight, high energy density and
ease of manufacture.

1.1

Lithium-Ion Battery

A basic Li-ion battery cell comprises a current collector, a negative electrode (anode), a separator, a positive electrode (cathode), and another current collector, as depicted in Figure 1.1
[1]. The material and microstructures of electrodes and separators vary widely. Figure 1.1
presents the materials used in today’s Li-ion batteries: a graphite negative electrode, a
Liy Mn2 O4 positive electrode, and a polymeric separator. The entire battery is immersed in
the electrolyte.
The electrodes are made of active particles and conductive additives. They are mixed with
the polymeric binders to form a paste and coated onto a metal foil (current collector). The
function of the binder is to maintain physical integrity and to help maintain an electrical

1

pathway between active material particles and the current collector. The size of active
particles (red and green balls), with which Li+ ions react and into which lithium inserts,
ranges from the submicrometer level to 10 µm.

Figure 1.1: Schematic of a Li-ion battery being charged [1]. The lower pictures in the ﬁgure
are scanning electron microscope (SEM) images of diﬀerent battery component materials.
From left to right, they are a current collector, which is a metal foil to collect current
from outside circuit, a negative electrode using graphite as active particle materials [1], a
polymeric separator (Celgard) [2], the positive particles using Liy Mn2 O4 active particles [3].
For interpretation of the references to color in this and all other ﬁgures, the reader is referred
to the electronic version of this dissertation.

For current Li-ion batteries, graphite is one of the most common anode material, because
2

of its relatively high speciﬁc capacity, small irreversible capacity, and good cycle ability.
Various kinds of carbonaceous materials, from crystalline to strongly disordered carbon, have
been examined as anodes in Li-ion battery [5]. Besides graphite, amorphous carbon, silicon,
nitrides, tin oxides and tin-based alloys have also been investigated as the anode material
[4]. The common cathode material in Li-ion battery could be spinel Liy Mn2 O4 , and other
layered oxides such as LiNi1−y Coy O2 , LiMnO2 , LiNi1−y Mny O2 , LiNi1−y−z Mny Coz O2
and Li1+x M1−x O2 [6]. The crystal structures of these cathode materials are similar to the
α-NaFeO2 type structure. Because of the compact crystal lattice, the α-NaFeO2 type of
structure has an inherent advantage in energy stored per unit volume.
An electrolyte contains ions and releases cations and anions by electrolysis. The electrolyte acts as a good ionic conductor. The cations and anions move towards cathode and
anode respectively thus conducting current. Generally, there are two types of electrolyte
in Li-ion batteries: non-aqueous liquid electrolyte, which is a mixture of lithium salt and
organic solvent; and gel-type polymer electrolyte, which is composed of polymer matrices
and organic solvents or plasticizers [7].
The electrochemical process in a Li-ion battery depends on a Li intercalation reaction.
The active material in the positive electrode is a compound in which Li is thermodynamically
stable. It may be any material that can store Li but at a lower thermodynamic potential
than the negative electrode material. The movement of ions is driven by the thermodynamic
force. When an external electrical potential is imposed, Li+ is forced to move from the
positive electrode to the negative electrode. At the surface of the negative electrode, Li+
reacts with an electron and either reduces to Li or forms an ion-electron pair, which then
moves into the electrode material by diﬀusion. When the external power supply is removed,
the reverse reaction takes place. A charge balance is maintained in the battery through the
3

movement of ions and electrons. At the graphite negative electrode, the reaction is

discharge
−− −− −
LiC6 − −− −− − xLi+ + Li1−x C6 + e−
− −− −−
charge

(1.1)

At the Lix Mn2 O4 positive electrode, the reaction is
discharge
−− −− −
xLi+ + xe− + Li1−x M n2 O4 − −− −− − LiM n2 O4
− −− −−
charge

(1.2)

For a Lix C6 |Polyproplyene(PP)|Li1−x Mn2 O4 cell, the reaction is
discharge
−− −− −
LiC6 + M n2 O4 − −− −− − Lix C6 + Li(1−x) M n2 O4
− −− −−
charge

(1.3)

To distinguish Li+ in the active particles from that in the electrolyte, it is designated as
Li in the following text.

1.2

Separator

The separator is a porous membrane that prevents physical contact between the positive and
negative electrodes while enabling ionic transport [2, 8, 9]. Three common types of separators
are polymeric membranes, nonwoven mats, and ceramic enhanced membranes. Currently,
polymeric membranes are used predominantly due to their low cost and relative thinness. A
thin separator facilitates the ionic transport and provide higher energy and power densities

4

[2]. The commercial PP membranes, such as Celgard series, with several dozen micrometer
thickness and 30-50% porosity are typical separator materials. The integrity of the separator
is vital to the performance and stability of Li-ion batteries. A short circuit due to separator
failure can lead to a thermal event [10, 11]. One of the major remaining challenges for largescale Li-ion batteries used in high energy applications, such as electric or hybrid-electric
vehicles is to further improve the performance of the separator.
It is well recognized that stress can have a signiﬁcant impact on the performance of
the separator. Currently, there is no method to evaluate the stress inside a battery. Consequently, there are no clear guidelines for the mechanical property requirements for the
separators. To improve the performance of battery separators, the stress in the separator
must be known. This includes the origin of the stress, its state, amplitude, and evolution.
This investigation is crucial to future separator development for improved battery abuse
tolerance and durability. It will contribute to the establishment of the relationships between
the stress and other physical phenomena involved in Li-ion batteries. This knowledge is
needed for the optimization of the separator materials, battery design, and charge-discharge
strategies for enhanced battery durability.

1.3

Problem Deﬁnition

In addition to the mechanical stress, two types of non-mechanical stresses can arise in a
Li-ion battery. The ﬁrst type is the intercalation stress, a unique problem for rechargeable
Li-ion batteries [12–15]. The second type of non-mechanical stresses is the thermal stress
due to self-heating of the battery during its operation or due to environmental temperature
variation.

5

The electrochemical process in a rechargeable Li-ion battery depends on a Li intercalation
and deintercalation reaction. One of the consequences of the reaction is the volume expansion
and contraction of the electrode materials. The intercalation stress arises due to the Li
concentration gradient inside an active electrode particle, which is the same as the diﬀusion
induced stress referred to in physics [16]. In the literature, the studies on intercalation stress
often focus on the stresses inside an electrode particle and the fracture of the particle itself
[12–15, 17–19]. In a battery, the eﬀect of intercalation stress is not limited to the individual
particles. Due to the constraint of the surrounding components, the expanding/contracting
particle will induce stresses in the surrounding materials.
When subjected to a temperature change, thermal stress arises in structures composed
of constituents with diﬀerent coeﬃcients of thermal expansion. A battery consists of a wide
variety of materials including metals, ceramics, and polymers. These components are stacked
up and packed into a given volume under pressure. Thermal stress is inevitable in such
structures. In the literature, the sources of heat generation [20, 22–24] and heat transfer
[21, 25–30] have been investigated for Li-ion batteries. However, there are no published
results concerning thermal stress.
In a Li-ion battery, a membrane separator is sandwiched between two electrodes. During
the charge and discharge cycles, the active particles expand in one electrode and contract in
another. To evaluate the stress in a separator in a battery during operation, both electrodes
must be considered.
At the present time, besides measuring the response of a battery pack in its through
thickness direction, there is no method to evaluate the stress in a separator inside a battery
with ongoing electrochemical reactions. Furthermore, the stress measured on a battery
pack provides an average global stress in the through thickness direction. At the local
6

level, the stress state and its amplitude are unknown. The electrodes used in state-of-theart Li-ion batteries are composites of discrete particles, binders, and other additives. The
surfaces of the electrodes are rough at the microscopic scale. Thus, the local stresses may
depart signiﬁcantly from the global stress measured at the surface. A rechargeable battery
is expected to last for a designated number of discharge-charge cycles. Thus the stress in
the separator may be cyclic. What its amplitude is and how to predict its fatigue life are
open questions.
Because the separator is a polymeric material, the temperature dependence of the separator’s mechanical properties is of great concern. A temperature variation of 20-30 ◦ C is not
uncommon for a battery during operation. For a polymer, this temperature variation may
cause a signiﬁcant acceleration in stress relaxation which in turn will inﬂuence the stress in
the material. Without the knowledge of the thermal history of the battery, it is diﬃcult to
estimate this inﬂuence.

1.4

Literature Review

In literature, some aspects of the above mentioned problems have been investigated. This
section provides an overview on current understandings relevant to this problem.

1.4.1

Mathematical Modeling of the Electrochemical Phenomenon

The most well-established model describing the electrochemical phenomenon in the porous
electrode was developed by Newman, s group [31, 32]. This model adopted much of the theory
developed previously from other battery systems, such as lead acid, nickel hydride or molten
salt LiAl-FeS [33]. Both the porous electrode theory and the concentrated solution theory

7

were employed in Newman’s model. The porous electrode theory treated the electrode homogeneously at the macro level without going into the exact geometric microstructure detail.
Most electrolyte solutions in Li-ion batteries exhibit the concentrated solution behavior as a
non-ideal solution. Therefore the concentrated solution theory was used to describe the Li+
transport in the liquid electrolyte. It includes all species present in the solution.
Comparisons of the simulation and experimental results [31, 32, 40, 41] indicated that
Newman’s model is suﬃcient to describe the electrochemical phenomenon in Li-ion batteries.
Thus various researchers incorporated this model in their work to investigate electrochemicalthermal [34, 39, 40] and electrochemical-mechanical [13–15, 35] phenomena.
Recently, volume averaging of generalized Poisson-Nernst-Planck equations on a de Levie
straight pore model was performed to simulate the electrochemical phenomenon of porous
electrodes by Lai et al.[38]. Compared to Newman’s model, these equations are more thermodynamically consistent, as they consider the innate correlation between the thermodynamic
and kinetic material properties. For example, the chemical diﬀusivity in Lai’s work is dependent on the other thermodynamic and kinetic properties of the battery. While it is treated
as a constant in Newman’s work.

1.4.2

Stress in a Li-Ion Battery

Non-mechanical stress arises in diﬀerent battery systems during the battery operation process. In a lead-acid battery system, the internal stress originates from two issues: volume
change due to the formation of sulfate and thermal expansion [42] due to battery self-heating.
In nickel-metal hydride batteries [43, 44], the non-mechanical stress generates by the lattice
expansion and contraction by hydrogen intercalation and deintercalation. Furthermore, the
numerical model including reaction kinetics, heat generation and mechanical analysis ad8

dressed there are another two types of non-mechanical stress: the chemical stress induced
by the oxygen potential gradient [45, 46] and the thermal stress induced by the mismatch
of thermal expansion of diﬀerent cell components [47]. Overall, the non-mechanical stress
generated in battery systems is from two aspects: material volume change by speciﬁc battery
electrochemical mechanism and thermal mismatch by battery temperature variation.
In a Li-ion battery, Li ions are in motion during charge and discharge. The ions either
intercalate into or deintercalate from the active electrode particles. This creates a Li concentration gradient from the center of the particle to its surface and therefore a non-uniform
strain ﬁeld in the particle. Similarly as the diﬀusion-induced stress referred to in physics
[16], intercalation stress arises in a particle with the Li concentration gradient.
The analysis of the intercalation stress requires the knowledge of the Li concentration
ﬁeld in the particles which is determined by battery kinetics. Christensen and Newman
[13, 14] presented a mathematical framework for the intercalation stress inside a spherical
particle. The diﬀusion in a spherical particle is one-dimensional (1D). The stress, diﬀusion
and battery kinetics were described by a set of 20 equations. These coupled equations were
solved by the ﬁnite diﬀerence method. For a more general shape, this approach becomes
unmanageable.
Another approach to this problem is the multiscale method. A common treatment is to
divide the diﬀerential equations into two sets. The ﬁrst set consists of the equations associated with the electrochemical kinetics of the battery; these include the expressions for Li+
transport and migration in the electrolyte (ionic charge balance, mass conservation), Ohm’s
law, and the interface reaction kinetic equation. The second set consists of the diﬀusion
equation and constitutive equations between stress and strain with the inclusion of eigenstrain (non-mechanical strain). This approach allows diﬀerent physical phenomena to be
9

represented by models at diﬀerent scales. Based on the porous electrode and concentrated
solution theories [48], a 1D macroscopic model has been well established for the electrochemical kinetics of Li-ion batteries [31, 32]. Utilizing this 1D battery model signiﬁcantly
reduces the complexity of the multiphysics problem. The Li+ ﬂux from the 1D model was
mapped to a three-dimensional (3D) model of a particle in a deformable matrix according
to its depth in the electrode. The stresses in the composite electrode were then evaluated
using a micromechanics model [49] based on the Mori-Tanaka eﬀective-ﬁeld theory [50].
Alternatively, the analysis of intercalation stress in particles can also be conducted by applying a prescribed Li+ ﬂux at the particle boundary without a battery model [12, 15, 18, 51–
53]. Zhang et al. [15] investigated the intercalation stress for a single active particle using
a detailed 3D ﬁnite element (FE) model using this approach. The study focused on the
stress in a particle with a ﬁxed boundary condition. Cheng et al. investigated 1D analytical solutions for intercalation stress under galvanostatic and potentiostatic conditions for
spherical particles [18] and cylinder nanowire electrodes [52]. These 1D analytical solutions
in nano-sized particles with the consideration of surface tension and surface modulus were
also provided [12]. Furthermore, these solutions for nano-sized particles were applied to investigate the cyclic stress inside a particle under a predeﬁned periodic potential boundary
condition [51]. On the other hand, Golmon et al. [53] investigated the intercalation stress
in a silicon spherical particle by applying the linear decrease overpotential to the particle
surface.
Two dimensional (2D) and 3D battery models have also been developed [19, 36]. Garcia
et al. [19] developed a 2D battery model for a Lix C6 —Liy Mn2 O4 cell and implemented it
in the object-oriented ﬁnite element program. This model considered the intercalation stress
in the composite cathode consisting of Liy Mn2 O4 particles in a matrix of an electrolyte.
10

The electrolyte was assumed to be a compliant gel with negligible resistance to deformation.
With a low volume fraction of active particles in the model, this treatment limited the
eﬀect of intercalation stress to individual particles. Wang and Sastry [36] developed a 3D
mesoscale battery model for the prediction of the electrochemical performance of the Li-ion
battery utilizing the Nernst-Planck equation in COMSOL® , a multi-physics FE package.
The stress analysis was not included in the mesoscale model.
So far, the analysis of intercalation stress has only been considered within individual
particles.
An experimental investigation of the intercalation stress in the thin battery electrode ﬁlm
was developed. A common technique is the laser beam deﬂection method (LBDM), initially
suggested by Stoney [54]. The stress in a thin ﬁlm deposited on a thick plate substrate will
cause the bilayer to bend. Stoney derived an equation relating the stress to the curvature of
the ﬁlm.
In Li-ion batteries, the curvature change of the electrode ﬁlm during battery charge
and discharge can be monitored by the LBDM method [55–58]. Using Stoney’s equation,
the in − situ stress of the electrode ﬁlm was calculated by the experimental value of ﬁlm
curvature. Several groups adopted this method [57] or a modiﬁed form such as the multibeam optical sensor(MOS) technique [59, 60] to study the stress in a large volume change
electrode material such as Si or Sn.
In open literature, there are no reports about the thermal stress in a Li-ion battery.

1.4.3

Thermal Analysis

The temperature variation during the battery working process not only causes the thermal
stress in the Li-ion battery, but also signiﬁcantly aﬀects the viscoelastic mechanical property
11

of polymeric separator. Therefore, understanding the heat generation is a necessary step in
assessing the stress within the separators of Li-ion batteries.
Newman and co-workers developed a number of heat generation models [20, 22, 23]
based on a framework derived by Bernardi et al. [20] from the view of energy balance in the
battery system. For Li-ion batteries, four heat generation sources were identiﬁed, namely:
irreversible resistance heating, reversible entropy heating, heat of mixing and heat of phase
change. The irreversible resistance heating includes joule heating and the energy dissipated
due to electrode overpotential. The reversible entropy heating is related to the change of
the Gibbs function of the cell contents. During electrochemical reactions, the composition
of individual species in each phase departs from their initial composition and hence the
total entropy of the cell varies. The heat of mixing represents the heat associated with
the generation or relaxation of concentration proﬁles. Thomas and Newman [23] found
that for a battery well designed to mitigate concentration overpotential, the heat of mixing
was negligible. Zhang et al. [24] implemented the heat generation model by Thomas and
Newman [23] into COMSOL® and examined the heat generation in a single Liy Mn2 O4
particle using a detailed 3D model. Their results conﬁrmed that the heat of mixing was
negligible compared to the irreversible resistance heating and reversible heating. Without
considering the side chemical reactions within a Li-ion battery, there is only one reaction
happening in each electrode, and no phase change occurs. Therefore, heat of phase change
was negligible [20, 61].
The heat generation expression proposed by Bernardi et al. [20] in its simpliﬁed form is
cited in the global heat generation model [25, 26, 29, 30, 62, 63]. The global model has
been applied in thermal analysis of a single cell [29] and cell stacks [30]. In this model, the
electrode potential is obtained based on the average composition.
12

Inside a battery cell, the local current and potential vary from point to point. To consider
this eﬀect, Rao and Newman [22] derived a local heat generation model in the subsequent
study. The heat associated with electrochemical reactions was considered at each electrochemical interface, which includes the irreversible resistance heat due to local overpotential
and reversible entropy contributions. The remaining irreversible resistance heat was further
divided into electronic and ionic types in which the electronic ohmic heat occurs in all solid
phases and the ionic ohmic heat occurs in the electrolyte phase. The local heat generation model has been incorporated into the electrochemical model to investigate the thermalelectrochemical behavior of a Lix C6 |LiPF6 |Liy Mn2 O4 [34, 39] and a Lix C6 |LiPF6 |Liy CoO2
cell [40, 64].
On the other hand, the heat generation rate of Li-ion batteries was also measured by
several groups by three primary experimental methods: accelerated-rate calorimetry [65–67],
isothermal heat conduction calorimetry [68–78, 103] and radiative calorimeter [80]. Most of
the experiments were focused on discharge or charge rates lower than 1C during discharge
[65–78, 80, 103] or charge [66–68, 70–77, 80, 103]. Due to the variations in battery materials
and the uncertainties in battery geometries [61], the reported volumetric heat generation
rate had a wide range. The peak values changed from 0.31 W L−1 [71] to 84.5 W L−1 [69]
under 1C discharge rate. Furthermore, the reversible entropy heating can be relatively large
under low rates [61, 68, 70, 77], and the consideration of this heat source is essential.

1.5

Objective of Research

The objective of this research is to develop a modeling tool for the evaluation of stresses in
the separator in a Li-ion battery.

13

In addition to the mechanical stress, the model must consider the intercalation stress and
thermal stress as well as the physics associated with these two stresses. Overall, the model
should include the battery kinetics, mass transport, thermal behavior, and stress analysis at
the level of a basic cell with all battery components. Such models have not been reported in
open literature thus far.
To accomplish this goal within the time frame of a Ph.D thesis, the model will be built
upon a validated 1D battery model that has been implemented in the multiphysics code
COMSOL® [32, 83]. It has two sub-models, “Battery”and “Electrode”. To consider the
stress, two additional sub-models, “Thermal”and “Stress”, will be developed. All coupling
relationships among diﬀerent sub-models will be considered. This will create a multiphysics,
multiscale model that allows one to assess the stress in the separator of a Li-ion battery.
The model development consists of the following tasks:
• “Stress”sub-model
• “Thermal”sub-model
• Coupling among the “Thermal”, “Battery”, and “Electrode”sub-models
• Coupling between the “Stress”and “Battery”sub-models
• Coupling between the “Thermal”and “Stress”sub-models
the physical phenomena considered in each sub-model and the coupling relationship between them are depicted in Figure 1.2.
Based on this multiphysics model, parameter studies will be conducted to investigate
the possibility of optimization of battery design parameters such as the thickness of the

14

Figure 1.2: The fully coupled multiphysics model consists of four sub-models. The schematic
shows the physical phenomena in each sub-model and the coupling relationships among four
sub-models.

components and the size of the active particles on the battery heat generation, as well as
the performance of the battery.
Several contributions to the ﬁeld of a Li-ion battery research are expected to be achieved
through this model:
• Improvement of the fundamental understanding of the stress in the separator;
• Enabling of the investigation of interrelated mechanisms and provide insights into the
interactions between diﬀerent physical phenomena in a Li-ion battery;
• Enabling of battery design optimization for better performance and longer cycle life;
• Estiblishment of useful information in the evaluations of battery control strategies.

15

1.6

Scope of the Thesis

The purpose of this work is developing a multiphysics model to investigate the stress in the
separator of a Li-ion Battery.
Chapter 1 provides the background information about this research, describes the objective and scope of the work, and presents the organization of this thesis.
Chapter 2 introduces the main phenomena occurring inside the Li-ion battery, and provides the corresponding governing equations and coupling relationships among various physical phenomena.
Chapter 3 describes the multiphysics model developed in this research and provides a
summary of the input parameters in the models.
Chapter 4 presents the simulation results of heat generation and temperature rise in a
battery. The temperature rise not only causes the thermal stress in the separator, but also
aﬀects the viscoelastic behavior of the polymeric separator. The relationship between the
temeprature and electrochemical reaction as well as cell utilization are discussed.
Chapter 5 presents the simulation results for the stress and strain in the separator. The
thermal eﬀect on the stress in the separator is examined. The eﬀect of discharge rate on the
stress in the separator is also involved in this chapter.
Chapter 6 describes the microstructure resolved model as the extension of the multiphysics model. The stress in the active particles and binders was examined based on the
realistic electrode microstructure geometry.
Chapter 7 summarizes the major ﬁndings in this work and suggests possible directions
for the future work.

16

Chapter 2
Multiphysics Phenomena
In order to evaluate the stress in the separator during the working process of the battery,
the electrochemical, thermal and mechanical phenomena in the battery must be considered.
The speciﬁc description of each phenomenon and the relationships among multiphysics phenomena are introduced in this chapter.

2.1

Electrochemical Phenomenon

In this work, Newman’s model [31, 32] is employed to describe the electrochemical phenomenon in a Li-ion battery. The porous electrode theory and the concentrated solution
theory are adopted in Newman’s model. The model considers an electrochemical system consisting of two porous electrodes and a separator saturated with the electrolyte. The porous
solid together with the electrolyte is modeled as a homogenized composite medium. The
volume fraction of the electrolyte equals the porosity. On the other hand, the concentrated
solution theory is used to describe the mass transport in the electrolyte. The electrochemical reaction that occurs at the active material surfaces is represented by the Bulter-Volmer
equation. The solid phase is denoted as phase 1 and the electrolyte phase is denoted as
phase 2. The variables and parameters that appear in both phases are distinguished by a
subscript 1 or 2.

17

2.1.1

Ionic Transport

In the electrolyte phase, Li+ is transported through two mechanisms: diﬀusion and migration. Based on the concentrated solution theory, the governing equation for material
conservation is:

ε2

dc2
=
dt

ef
· D2 f c2 −

1
i2 ·
F

t+ +

Sa jn
(1 − t+ )
F

in Ωe , the electrolyte domain

(2.1)

ef f
where c2 is the Li+ concentration, t is the time, D2 is the eﬀective diﬀusion coeﬃcient
in the electrolyte, ε2 is the porosity (volume fraction of the electrolyte), F is Faraday’s
constant, t+ is the cation transport number representing the percentage of the current in
the electrolyte carried by Li+ , the vector i2 is the current density in the electrolyte, jn is
the charge transfer current density at the interface, and Sa is the speciﬁc surface area of the
electrode, which is deﬁned as the surface area of the active materials per unit volume of the
total electrode. For the spherical particle, Sa is expressed as:

Sa =

3ε1
R

18

(2.2)

where ε1 is the volume fraction of the active materials, and R is the active particle radius.
The relationship between i2 and jn obeys the following equation:

· i2 − Sa jn = 0

(2.3)

Equation 2.1 states that the change in Li+ concentration is the net result of Li+ diﬀusion
(ﬁrst term), Li+ migration in an electrical ﬁeld (second term), and the concentration change
induced by the Li+ that leaves or enters the electrolyte through the interface with the active
particles (third term). The third term vanishes for a volume that is not adjacent to an active
particle.

2.1.2

Ionic Charge Balance

For ionic charge balance in the electrolyte phase, the governing equation is given as
ef f
2Rc T k2
ef f
· −k2
φ2 +
F





1+

∂ ln f±
∂ ln c2

1 − t+

(ln c2 ) = Sa jn

in Ωe , the electrolyte domain

(2.4)

ef f
Here Rc denotes the gas constant, T denotes the temperature, kk2
is the eﬀective
ionic conductivity, φ2 is the ionic potential, and f± is the mean molar activity coeﬃcient.
Equation 2.4 is derived for an electrolyte of a binary salt using concentrated solution theory.

19

It can be considered as a modiﬁed Ohm’s law. The ﬁrst term in the bracket is equivalent
to the expression describing the electron movement in a solid conductor. The second term
accounts for the eﬀect of Li+ concentration on the ionic current [48].
In order to consider the eﬀect of material porosity and tortuosity on the ionic diﬀusion
coeﬃcient and conductivity, the eﬀective values of these parameters are used in equation
2.4. Both the ionic diﬀusion coeﬃcient and conductivity depending on the porosity and
tortuosity follow the relationship:

ε D
ef f
D2 = ε2 p D2 = 2 2
τt

(2.5)

ε k
ef f
k2 = ε2 p k2 = 2 2
τt

(2.6)

where D2 and k2 are the ionic diﬀusion coeﬃcient and conductivity measured from the pure
electrolyte, p is a constant, which is the exponent in the above empirical relations for eﬀective
transport properties, and τt is the tortuosity of the battery material.
For a system made of a continuous conductive phase, for example the dilute solution with
the porosity near zero, the value of p is 1.5 based on the asymmetrical Bruggeman eﬀective
model [38, 48]. This value sometimes is used as a default value in battery models. However,
the porosity of the normal porous battery material is around 0.4 or even higher. Thus the

20

default value 1.5 is unable to correctly predict the mass transport resistance [38, 48, 84, 85].
One way to obtain the value of the tortuosity or exponent p of the electrode or separator
is using discharge curves to ﬁt as the adjustable value [32, 41, 88, 89, 133]. The mathematical
model of porous networks in relation to battery electrodes and separators is another way to
get the value [85]. Recently, the experimental measurement of the tortuosity of the porous
electrode and separators was developed [89, 90]. All of three methods depict that the value
of p is higher than the default value 1.5.

2.1.3

Mass Transport

The active materials in the electrodes are considered to be spherical particles. The Li
diﬀusion in the active particles follows Fick’s second law with the spherical coordinate:

∂c1
1 ∂
+
∂t
r2 ∂r

−r2 D1

∂c1
∂r

=0

in Ωs , the solid, active particle domain

(2.7)

where c1 is the Li concentration, D1 is the solid phase diﬀusion coeﬃcient, and r is the
radial position in the spherical particle.
At the particle surface, the Li+ ﬂux is related to the local current density by:

N0 = −jn /F
21

(2.8)

where N0 is the Li+ ﬂux at the particle surface.

2.1.4

Electronic Charge Balance

The charge balance in the solid phase of the electrode obeys Ohm’s law:

· (−k1 φ1 ) = i1

in Ωs , the solid, active particle domain

(2.9)

where φ1 is the electronic potential of the solid phase and k1 denotes the electronic conductivity of the porous electrode. It is found by measuring the resistance to the passage of
current through the porous electrode. The vector i1 is the current density in the solid phase
and

· i1 + Sa jn = 0

22

(2.10)

2.1.5

Reaction Kinetics

The electrochemical reactions occur at the active particle surface. The reaction kinetics is
described by Butler-Volmer equation:

jn = j0 exp

ηαn F
RT

− exp

(−η)αp F
RT

on Γi , the active particle/electrolyte interface

(2.11)

where αn and αp are the transfer coeﬃcient of electrochemical reaction in the negative and
positive electrode respectively. η is the local surface overpotential for the electrode, which is
deﬁned as

η = φ1 − φ2 − Ui (c1 ,surf )

(2.12)

and j0 denotes the exchange current density. This is determined by

j0 = k0

c2 (c1,max − c1,surf )c1,surf

23

(2.13)

In equations 2.12 and 2.13, k0 is the reaction rate constant, c1,surf is the Li concentration
at the surface of the particle, c1,max is the maximum achievable Li concentration in the
particle, and Ui is the electrode’s open circuit potential (OCP) or equilibrium potential,
which is the potential that no current ﬂows. The OCP can be treated as a function of the
intrinsic nature of the species and is inﬂuenced by concentration or temperature.

2.2

Thermal Phenomenon

The heat transfer in a 1D isotropic material is described by:

ρCp

∂T
= Ktherm 2 T + q
∂t

in Ωb , the battery domain

(2.14)

where Ktherm is the thermal conductivity, ρ is the density, Cp is the speciﬁc heat capacity
of the material, and q is the heat generation rate.
The heat generation accompanying the electrochemical phenomenon is the critical part
of the thermal phenomenon in a Li-ion battery. Currently, both the global model and the
local model are used to simulate the heat generation in a Li-ion battery [61].

24

2.2.1

Heat Generation: Global Model

The general heat generation expression for Li-ion batteries was ﬁrst derived by Bernardi et
al. [20] from the energy balance of the battery system. The general expression for the heat
generation rate is [20]

q = I(V − U ) + IT

dU
ave
+
∆Hp γp +
dT
p

j

i

∂cij
¯
¯ ave
(Hij − Hij )
dv
∂t

(2.15)

where U is the OCP of the full Li-ion cell, I is the cell current, V is the cell potential,
¯
Hij is the partial molar enthalpy of species i in phase j, γp is the reaction rate of a side
ave
chemical reaction p and ∆Hp is the average enthalpy change in the chemical reaction p.
In Equation 2.15, the four terms correspond to the heat generation by irreversible resistance
heating, reversible entropy heating, heat of phase change from side chemical reactions, and
heat of mixing, respectively.
The global heat generation model has been applied in thermal analysis of a single cell
[29] and cell stacks [25–27, 30, 62]. Due to the lack of data, the heat of mixing and phase
change were neglected in these studies and Equation 2.15 is further reduced to

q = I(V − U ) + IT

25

dU
dT

(2.16)

2.2.2

Heat Generation: Local Model

In order to consider the local current and potential variation from point to point in a Li-ion
battery, Rao and Newman [22] derived a local heat generation model. The local model is
adopted in this work.
In the local model, the total heat generation is obtained by adding up the contributions
at electrochemical interfaces and in bulk phases of the cell [22], thus

q = qp + qs + qn

(2.17)

where subscripts p, s and n correspond to the positive electrode, the separator, and the
negative electrode, respectively.
For a porous electrode, the heat generation at an electrochemical interface includes a
reversible part due to the entropy change in the electrode, Π, and an irreversible part due
to the local surface overpotential, η, which is expressed by Equation 2.12 for both negative
and positive electrodes. The reversible part Π is described as:

∆Si
Π=T
F
i

∆Si
∂Ui
F is proportional to ∂T [48, 64], such that

26

(2.18)

∆Si
dUi
=
nF
dT

(2.19)

where n is the charge number contributing to the reaction. For a Li-ion battery, n=1 and,
therefore,

∆Si
dUi
=
Π=T
F
dT
i

(2.20)

where ∆Si is the electrode entropy change (i=p or n). It is typically measured in the
reduction reaction in the working electrode, which is the discharge reaction for a cathode in
a half cell, where the anode typically is Li metal foil. Considering the convention, the sign
of heat generation can be either negative or positive for the exothermic reaction. When the
heat generation for the exothermic reaction is positive, the anodic current is positive, and
the cathodic current is negative. In this work, the positive sign for heat generation in an
exothermic reaction is adopted.
The heat generation of a porous electrode in a local area is given by

qi = Sai in (ηi + Πi ) + i1 ·

27

φ1 + i2 ·

φ2

(2.21)

where i1 and i2 are the current densities in the solid phase and the electrolyte phase, respectively.
In the separator, the current is conducted through ion transport. The heat generation
due to the ionic ohmic heat in the porous separator is given by

qs = i 2 ·

φ2

(2.22)

The heat generation in each current collector is calculated by

qc = i2 Rcol
col

(2.23)

where icol is the current density in the current collector and Rcol is the resistivity of the
current collector.
Thomas and Newman [23] have developed a heat generation model which considered
the heat of mixing and compared the simulation results with measured heat of mixing for
a Li|LiPF6 (EC/DMC)|LiAl0.2 Mn1.8 O4−δ F0.2 cell. It suggests that the heat of mixing
during high rate pulses may be negligible considering that the concentration relaxation is
relatively small at these rates. Therefore, the heat of mixing is not considered in this work.
In addition, the heat from the side reaction is also not considered.

28

2.3

Mechanical Phenomenon

Besides mechanical stress, the intercalation stress and thermal stress arise during the working
process of the battery.
ei,c
The diﬀusion-induced eigenstrain εij (non-mechanical strain) in an isotropic electrode
material is considered through a diﬀusion/thermal analogy such that [12, 16]:

ei,c 1
εij = ∆c1 Ωδij
3

(2.24)

where ∆c1 is the change in the concentration of the diﬀusion species, Ω is the partial molar
volume of the solute in the host material, and δij is the Dirac delta function.
The partial molar volume can be calculated using the following Equation [15]

Ω=

ei,c
3εij
∆ycs,max

(2.25)

ei,c
where εij is the strain caused by intercalation expansion along the major axis of the
material, ∆y refers to the change in the number of lithium atoms that occurs during the
charge and discharge cycles, and cs,max is the stoichiometric maximum concentration in the
electrode.
On the other hand, stress will aﬀect the diﬀusion process. The eﬀect of stress on diﬀusion
has been considered through the chemical potential of a solute in a solid solution [91] as:

29

µ = µ0 + RT ln c1 − σh Ω

(2.26)

where µ0 is the chemical potential at the reference condition, and σh is the hydrostatic
component of the stress, which is deﬁned as:

σh =

σxx + σyy + σzz
3

(2.27)

where σxx , σyy and σzz are the normal stresses along x, y, and z directions respectively.
The driving force for diﬀusion is a gradient in chemical potential. In this sense, the
diﬀusion ﬂux is expressed by the following equation:

N0 = −M c1 µ = M c1 RT

c1
− Ω σh
c1

(2.28)

where M is the mobility of the solute. Noting that the diﬀusion coeﬃcient D1 = M RT ,
Equation 2.28 changes to:

N0 = −D1

c1 −
30

Ωc1
σh
RT

(2.29)

The law of mass conservation is:

· N0 = −

∂c1
∂t

(2.30)

Substituting Equation 2.29 into Equation 2.30, the diﬀusion equation with stress eﬀects
is derived as

D1

2c − Ω
1 RT

c1 ·

σh −

Ωc1 2
σh
RT

=

∂c1
∂t

(2.31)

In Equation 2.31, D1 is assumed to be independent of stress.
The mechanical property evolution as a function of Li concentration is a very active
research topic [92–94].The ﬁrst-principles density functional theory is employed to reveal
the concentration dependence of Young’s modulus of electrode materials such as graphite or
silicon [92, 94]. However, at the present time, there is no one generally accepted model. In
this study, the Li concentration dependence of the mechanical property is not considered.
ei,T
The other eigenstrain εij caused by thermal expansion or shrinkage is calculated by:

ei,T
εij = αij ∆T δij

31

(2.32)

where αij is the coeﬃcient of thermal expansion.
For an isotropic elastic material, the constitutive relationship for the mechanical strain
εme due to the mechanical stress is given as [95]
ij

εme =
ij

1
(1 + ν)σij − νσkk δij
E

(2.33)

where E is Young’s modulus and ν is Poisson’s ratio of the material. In this work, E and
ν were treated as constants. The variation of the elastic constants with temperature and Li
concentration was not considered.
The total strain is the summation of two types of eigenstrains and the mechanical strain
in the elastic battery material as:

ei,T
ei,c
εij = εme + εij + εij
ij
1
1
= E (1 + ν)σij − νσkk δij ) + αij ∆T δij + 3 ∆cΩδij

2.4

(2.34)

Multiphysics Coupling

The electrochemical, thermal, and mechanical phenomena contemporarily exist in the battery when it is in the operation process. They are highly related with one another, which is
clearly demonstrated in Figure 2.1.
32

Figure 2.1: The coupling relationships among thermal, electrochemical and mechanical
phenomena.

According to the description of the thermal phenomenon in Section 2.2, the electric and
ionic transport in the bulk materials and the interface reaction at the active particle surfaces
induce the heat generation in the Li-ion battery.
On the other hand, mass transport and reaction rate strongly depend on temperature.
The inﬂuence of temperature on the mass transport and reaction rate is considered for six
variables. They are the solid phase diﬀusion coeﬃcient D1 , reaction rate constant k0 , OCP,
electrolyte ionic conductivity k2 , electrolyte diﬀusion coeﬃcient D2 , and the product of the
thermodynamic factor and the anion transference number v. D1 , OCP, and k0 are related
to the electrodes; k2 , D2 and v are related to the Li+ transport in the electrolyte system.
The concentration gradient causes intercalation stress in the active electrode particles in
a Li-ion battery. In addition, the temperature variation due to the heat generation causes

33

the thermal stress in a Li-ion battery.
The inﬂuence of temperature on the viscoelastic behavior of the polymeric separator is
considered by using the time-temperature superposition principle [96]. A thermo-viscoelastic
model based on the experimental data will be implemented to the multiphysics model.
The inﬂuence of Li concentration on the elastic properties of the active particles as well
as the eﬀect of stress on Li diﬀusion in solid particle will not be considered in the current
work.

34

Chapter 3
Multiphysics Model
A multiphysics model for the stress analysis of a basic Li-ion battery cell has been developed.
In this model, the multi-physics equations were solved concurrently using four sub-models
at diﬀerent scales. The battery electrochemical kinetics, current balance, ionic transport
and mass transport were solved by a 1D macroscopic “Battery”sub-model and a pseudo 2C
“Electrode”sub-model based on an existing 1D battery model in COMSOL® . The thermal
behavior was considered by the 1D “Thermal”sub-model. The stress analysis at the level
of a battery cell was conducted with a meso-scale representative volume element (RVE)
“Stress”sub-model of the battery cell. Based on this approach, a multi-scale model with
the inclusion of intercalation stress was developed for a Lix C6 |PP|Liy Mn2 O4 cell, with the
electrolyte of 1M LiPF6 in a solvent of ethylene carbonate (EC)/dimethylene carbonate
(DMC) (2:1 by volume).

3.1

Model depiction

A major diﬃculty in modeling the stress due to Li intercalation at the battery cell level is
the multi-scale, multi-physics nature of the problem. A separator is sandwiched between
two electrodes in a basic cell. In the thickness direction, a basic cell ranges from 150 to 300
µm. The intercalation stress occurs at the level of individual active particles which have an
average size of 5-20 µm. On the other hand, in a planer (pouch) cells used in large-scale

35

battery, the in-plane dimension of a single cell is in the range of 20-30 cm. The aspect ratio
for a single planer cell is at the order of 103 .
To overcome the diﬃculty of the large aspect ratio of a planer cell, a RVE approach was
employed together with a multi-scale approach. RVE by deﬁnition is the smallest material
volume element of a heterogeneous material for which the eﬀective constitutive representation
is a suﬃciently accurate model to represent the mean constitutive response of the material
[97]. In this work, the concept of RVE is applied in a slightly diﬀerent sense. RVE is the
smallest volume of a battery cell that represents the mechanical and physical response at
the center region of a planar cell. Ideally, a RVE model should be suﬃcient to reveal the
response of individual components in a battery. The RVE used in this work is illustrated
in Figure 3.1. The deformation of the RVE is restricted in the plane of the pouch cell and
hence may be considered under a plane strain condition and analyzed with a 2D RVE.
Figure 3.2 presents the schematics of the individual sub-model implemented in COMSOL® .
In this model, the equations of the multi-physics problem are divided into four sets and solved
by four sub-models at three diﬀerent scales. The battery electrochemical kinetics is solved
by a 1D macroscopic “Battery”sub-model based on that of Doyle et al. [31, 32]. The Li
diﬀusion in the active particles and the interaction between the particle diﬀusion and battery kinetics is solved using a microscopic “Electrode”sub-model for spherical particles in the
electrode. The thermal behavior is solved by the 1D macroscopic “Thermal”sub-model. The
stress analysis was carried out at the level of a basic battery cell using a RVE. Compared
to the other sub-models, the RVE “Stress”sub-model is at a meso-scale. It allows the investigation of mechanical responses of individual components in a battery environment. The
ﬁrst two sub-models have been reported in literature [32] and available in COMSOL® for
Lix C6 |PP|Liy Mn2 O4 battery cell. The current multiphysics model is built upon the 1D
36

Figure 3.1: The meso-scale representative volume element (RVE) represents the stress state
at the center of the battery cell.

battery model in COMSOL® . A description of the ﬁrst two sub-models and their coupling
are provided in the following text, in addition to the discussion of “Thermal”sub-model and
“Stress”sub-model.

3.1.1

1D Battery Model

The Newman’s 1D Battery model, including the 1D “Battery”sub-model and pseudo 2D
“Electrode”sub-model, was used to describe the electrochemical phenomenon in the Li-ion
battery.

37

Figure 3.2: COMSOL® implementation of the fully coupled thermal-electrochemical-stress
model. The electrochemical performance of the battery is evaluated with an existing 1D
Li-ion battery model in COMSOL® which consists of two sub-models; (a) “Battery”and (b)
“Electrode”. The thermal analysis is performed using (c) “Thermal”sub-model. The stress
analysis for the components in a cell is performed concurrently using (d) a meso-scale RVE
“Stress”sub-model.

38

3.1.1.1

“Battery”Sub-Model

The ionic charge transport, charge balance, and reaction kinetics in the negative electrode,
separator and positive electrode were described in the 1D “Battery ”sub-model. The current
collectors were not involved. The two ends of the 1D geometry in the “Battery”sub-model
as shown in Figure 3.2(a) represented the electrode/current collector boundaries.
Substituting Equation 2.3 into Equation 2.1, the ionic charge transport in the electrolyte
was derived as:

ε2
=

dc2
=
dt

1
eﬀ
· D2 c2 − F i2 ·
i2
eﬀ
D2 c2 + F (1 − t+ )

·

·i
t+ + F 2 (1 − t+ )

(3.1)

The dimensionless length method was adopted in the 1D “Battery”sub-model. Consequently, Equation 3.1 was changed to:

dc
Li ε2 2 −
dt

·

eﬀ
D2
Li

c2

=

· i2

Li
F

(1 − t+ ) =

Li Sa jn
(1 − t+ )
F

(3.2)

In the same manner, Equations 2.4 and 2.9, which describe the ionic balance in the
electrolyte and electronic balance in the solid phase, respectively, were changed to:



 k ef f

∂ ln f±
2RT
2
·
− φ2 +
(1 +
)(1 − t+ ) (ln c2 )
= Li S a jn
 Li

F
∂ ln c2
39

(3.3)

and:

ef f
k
− 1
·
L2
i


Li


φ1  = −Li Sa jn

(3.4)

Here, the local charge transfer current density jn was described by the Bulti-Volmer
equation, as shown in Equation 2.11.

3.1.1.2

“Electrode”Sub-Model

In the pseudo 2D “Electrode”sub-model as shown in Figure 3.2(b), the dimensionless radius
y, y = r/R, where R was the particle radius, was introduced to Li solid state diﬀusion.
Equation 2.7 was rearranged as:

y2R

∂c1
∂
+
∂t
∂y

D ∂c
−y 2 1 1
R ∂y

=0

(3.5)

In Figure 3.2(b), there are two squares in the pseudo 2D “Electrode”sub-model representing negative and positive electrodes, respectively. The x-axis (horizontal-axis) corresponds
to the position along the width of the electrodes and the y-axis (vertical-axis) represents
the particle radius. The Li concentration at the surface of the active particles, y = 1, was
coupled to the 1D “Battery”sub-model for the battery kinetics calculation.
The Li diﬀusion was considered within the individual spherical particles in this model.
40

The mass transport among diﬀerent spherical particles was neglected. Therefore, the anisotropic
diﬀusion is employed in the pseduo 2D “Electrode”model. The negligible diﬀusion coeﬃcient
D1
was applied in the x-axis direction, while a scaled diﬀusion coeﬃcient ( rp ) was applied in
the y-axis direction.

3.1.1.3

Boundary Conditions

The model was constructed by a galvanostatic charge/discharge battery operation mode.
By the electronic balance, the applied current density iapp was speciﬁed at the positive
electrode/current collector boundary (a2 in Figure 3.2(a)). It is:

i1 = iapp

on Γp = 3, the positive electrode/current collector interface

(3.6)

The cell voltage was calculated relative to the negative electrode. Therefore, the electronic
potential of the solid phase at the interface of the negative electrode and current collector
(a1 in Figure 3.2(a)) is set to zero:

φ1 = 0

on Γn = 0, the negative electrode/current collector interface
41

(3.7)

For the ionic charge balance in the electrolyte, the electrode/current collector boundaries
were insulating:

i2 = 0

on Γn = 0,

Γp = 3

(3.8)

For the ionic transport, ions do not diﬀuse to the current collectors. Therefore, the
insulating condition was applied at the electrode/current collector boundaries:

c2 = 0

on Γn = 0,

Γp = 3

(3.9)

The lithium ion ﬂux N0 calculated from the 1D “Battery”sub-model was mapped to the
surface of the particles, which is the boundary of y = 1 in the “Electrode”sub-model (see
Figure 3.2(b) ). All the other boundaries for the pseudo 2D model were insulated.

42

3.1.1.4

Model Implementation

For the 1D “Battery”sub-model, the Equations 3.3 and 3.4 describing ionic balance in the
electrolyte and electronic balance in the solid phase were implemented into PDE Interfaces
mode. The ionic transport in the electrolyte was modeled using Diﬀusion and Convection
Application mode in COMSOL® . This mode is also used to describe the mass transport in
active particles in the pseudo 2D “Electrode”sub-model.

3.1.2

“Thermal”Sub-Model

The 1D “Thermal”sub-model considered the heat generation and thermal analysis of a basic
Li-ion battery consisting of a negative electrode, a separator, a positive electrode and two
current collectors. The heat generation was considered by the local heat generation model
as described in Section 2.2.2. As the same as the 1D “Battery”sub-model, The electrodes
and the separator are considered as composites consisting of solid and liquid phases. The
current collectors contribute to the thermal energy balance only.
The 1D “Thermal”sub-model has a symmetric boundary condition (insulated) at the end
of Cu current collector(b1 in Figure 3.2(c)) and a lumped heat transfer coeﬃcient of h=0.2
W m−2 K−1 assigned at the end of Al current collector (b2 in Figure 3.2(c)).The lumped
heat transfer coeﬃcient represented the combined eﬀect of heat conduction, convection heat
transfer and radiative heat transfer.
The “Thermal”was implemented into the Heat Transfer mode in COMSOL® .

43

3.1.3

“Stress”Sub-Model

To carry out the stress and deformation analysis, the model should be at least in 2D. The
“Stress”sub-model presented in Figure 3.2(d) was a prototype quarter-sized 2D model. From
left to right, the components included were current collector, negative electrode, separator,
positive electrode, and current collector. The circles in the electrodes represent the active
particles. In its thickness direction, the electrode in the quarter-sized sub-model had a
dimension that was equivalent to a quarter of the original electrode thickness.
The multi-physics phenomena considered in the “Stress”sub-model included the diﬀusion
of Li inside the electrode particles and the stress due to Li intercalation and thermal expansion, as well as mechanical loading. To reduce the modeling eﬀort and allow active particles
of diﬀerent sizes and shapes to be considered, the diﬀusion process in the “Stress”sub-model
had one-way coupling with the 1D “Battery”sub-model. The Li ﬂux in the 1D “Battery”submodel was mapped to the active particle surface in the “Stress”sub-model. In this way, even
if the volume fraction of the active particles of a certain shape and packing pattern cannot
be realized accurately in the “Stress”sub-model, it would not aﬀect the battery kinetics in
the 1D “Battery”sub-model. The temperature rise obtained from the “Thermal”sub-model
was also mapped to the “Stress”sub-model. The deformation from the mechanical stress,
volume changes due to Li intercalation and thermal expansion, as shown in Equation 2.34,
are considered in this analysis. At this stage, the eﬀects of the stress on the Li diﬀusion as
expressed by Equation 2.31 is not considered.
The boundary condition of the stress analysis is presented in Figure 3.3. A periodical
boundary was assigned at the lower and upper boundaries in the in-plane direction. A ﬁxedfree boundary condition was prescribed in the through-thickness direction where the left end

44

was ﬁxed while the right end was subjected to a pressure normal to the boundary. The
pressure level investigated in this work was 10 psi.
The Li solid state diﬀusion and stress analysis was implemented into the Diﬀusion and
Convection Application mode and the Solid Mechanics (plane strain) mode in COMSOL®
, respectively.

Figure 3.3: The boundary conditions for stress analysis.

45

3.2

Material Properties and Simulation Condition

The battery considered in this work is a Lix C6 |LiPF6 (EC/DMC)|Liy Mn2 O4 cell. The
electrolyte is the solution electrolyte of 1M(1000 mol m−3 ) LiPF6 dissolved in a mixture of
EC and DMC with a volume ratio of 2:1.

3.2.1

Temperature Dependent Variables

The temperature dependence is considered for six sets of physical properties. These are
(1) k0 , the reaction rate constant; (2) D1 , the diﬀusion coeﬃcient of Li in the solid active
particles; (3) OCP of the electrodes; (4) D2 , the diﬀusion coeﬃcient of Li+ in the electrolyte;
(5) k2 , the electrolyte ionic conductivity; and (6) v= 1 +

∂ ln f±
∂ ln c2

1−t+ , the product of

the thermodynamic factor and the anion transference number. Table 3.1 provides a summary
for the expressions used to describe their temperature dependency, Equations 3.10 through
3.15. The decoupled model is the referred as the model that the temperature dependence of
the physical properties is not considered.

46

Table 3.1: Temperature dependency of six sets of physical properties.

Physical Property

Expressions

Equation

Ref.

3.10

[64]

3.11

[64]

3.12

[64]

3.13

[106]

3.14

[106]

3.15

[32, 106]

Active material
(1) Li Ddiﬀusion coeﬃcient
(2) Reaction rate constant
(3) Open circuit potentials

Ea
1
D1 = D1,ref exp[ RD ( T 1 − T )]
ref
EaR 1
1
k0 = k0,ref exp[ R ( T
− T )]
ref
∂Ui
U = U (SOCsurf , Tref ) +
(T − Tref )
∂T SOC
,Tref
surf
∆Si
= nF
(T − Tref )
SOCsurf ,Tref

Electrolyte
(4) Ionic conductivity
(5) Li diﬀusion coeﬃcient
(6) v= 1 +

∂ ln f±
∂ ln c2

1−t+

m n
k2 (c2 ,T )
li,j ci T j
=
c2
2
i=0 j=0
54
− 0.22c2
log10 D2 = −4.43 −
T −(229+5c2 )
v = 1−0.363 0.601 − 0.24c2 0.5 0.982 [1 − 0.0052(T − 294)] c2 1.5
1−0.399

47

3.2.1.1

Variables Related to the Electrodes

The variables related to the electrodes include k0 , D1 and OCP. For Li in graphite and
Liy Mn2 O4 , the temperature dependences of these variables has been reported in open literature [32, 39, 64, 99–105].
Solid state diﬀusion coeﬃcient D1 and the reaction rate constant k0
The temperature dependences of D1 and k0 follow the Arrhenius relationship [64], as
shown in Equations 3.10 and 3.11. The constants at Tref =25 ◦ C for the battery system
considered in this work are listed in Table 3.2. Figures 3.4(a) and (b) present D1 and k0 as
a function of temperature, respectively. The state of charge (SOC) dependence of D1 and
k0 was not considered at this stage of the work due to lack of data.

Table 3.2: Solid phase diﬀusion coeﬃcient D1 and reaction rate constant k0
Property
Solid phase diﬀusion coeﬃcient at 25 ◦ C
D1,ref (m2 s−1 )
Diﬀusion coeﬃcient activation energy
EaD (KJ s−1 )
Reaction rate constant at 25 ◦ C
k0,ref (A m2.5 mol−1.5 )
Reaction activation energy
EaR (KJ mol−1 )

OCP

48

Negative electrode
(Lix C6 )

Positive electrode
(Liy Mn2 O4 )

3.9 × 10−14 [32]

1 × 10−13 [32]

35 [99]

20 [39]

2 × 10−6 [32]

2 × 10−6 [32]

20 [100]

53.07 [101]

(a)

(b)

Figure 3.4: Temperature dependence of (a) the diﬀusion coeﬃcient D1 and (b) reaction
rate constant k0 for Li in the negative (graphite) and positive (Liy Mn2 O4 ) electrodes. The
values of D1 at 25 ◦ C are 3.9 × 10−14 m2 s−1 for graphite and 1 × 10−13 m2 s−1 for
Liy Mn2 O4 . The value of k0 at 25 ◦ C is 2 × 10−6 A m2.5 mol−1.5 for both electrodes.

The temperature dependence of the OCP can be expressed with a Taylors expansion
[64], Equation 3.12. As shown, the OCP is a function of the local SOC at the surface of
the active particle , the entropy change ∆S, and temperature. For graphite and Liy Mn2 O4
electrodes, these two functions at 25 ◦ C have been obtained by experiments [32, 103–105].
Figures 3.5(a) and (b) represent the SOC dependent ∆S at 25 ◦ C for the two electrodes.
Figures 3.6(a) and (b) represent OCP as a function of SOC and temperature calculated by
Equation 3.12 for graphite and Liy Mn2 O4 , respectively. As shown, the eﬀect of temperature
is negligible on OCP.

3.2.1.2

Variables Related to the Li+ Transport in the Electrolyte

Due to the lack of experimental data, the temperature dependence of Li+ transport in LiPF6
in EC/DMC (2:1 by volume) was assumed to follow the trends of similar electrolyte systems.
49

(a)

(b)

Figure 3.5: The entropy change as a function of SOC at 25 ◦ C for (a) graphite [104] and
(b) Liy Mn2 O4 [105].

(a)

(b)

Figure 3.6: OCP of (a) graphite and (b) Liy Mn2 O4 vs. lithium metal as a function of
SOC and temperature. The curves at the reference temperature 25 ◦ C were determined by
experiments [32, 104].

50

Ionic conductivity k2
The ionic conductivity k2 for LiPF6 in PC/EC/DMC (10:27:63 by volume) has been
investigated experimentally [106]. The concentration and temperature dependence of k2
was found to follow an empirical relationship given by Equation 3.13 and the PC/EC/DMC
system (10:27:63 by volume) exhibited a trend such that k2 → 0 as c2 → ∞. Combining
the parameters for LiPF6 in EC/DMC (2:1 by volume) at 25 ◦ C [32] with the trend of
the temperature dependence for the PC/EC/DMC system [106], for m=n=2, the following
expression was obtained:

k2 (c2 ,T )
= 1.12 × (−8.2488 + 0.053248T − 0.000029871T 2 + 0.26235c2
c2
−0.0093063c2 T + 0.000008069c2 T 2 + 0.22002c2 2 − 0.0001765c2 2 T )

(3.16)

where the unit of c2 is mol L−1 , the unit of T is Kelvin, and the unit of k2 is mS cm−1 .
Figure 3.7 plots k2 as a function of Li+ concentration at temperatures from -20 to 60 ◦ C.
The value of k2 used in the decoupled model is from Ref. [32].
Diﬀusion coeﬃcient D2
The diﬀusion coeﬃcient D2 in LiPF6 in PC/EC/DMC (10:27:63 by volume) was studied
in [106]. The experimental results of the temperature and concentration dependence of D2
were ﬁtted to Equation 3.14. The magnitude of D2 given in Ref. [106] is similar to what
was reported for LiPF6 in EC/EMC for a set of volume ratios [107] and LiN(CF3 SO2 )2 in
EC/DEC (2:3 by volume) with a copolymer PVdF-HFP as a host [108]. The above results
indicate that the transport property of LiPF6 in carbonate solvents is strong functions of
c2 and T but a much weak function of the formulation of the carbonate mixture. Based on
51

Figure 3.7: The ionic conductivity k2 as a function of the salt concentration and temperature
in 1M LiPF6 in EC/DMC (2:1 by volume), calculated with Equation 3.16.

this observation, Equation 3.14 was employed to describe the temperature and concentration
dependence of D2 in the current work, even though the solvent mixture was diﬀerent.
Figure 3.8 presents D2 as a function of the salt concentration at diﬀerent temperatures.
As seen, using a constant D2 can lead to signiﬁcant errors when the local environment
deviates from the standard condition of 1M concentration and 25 ◦ C.
Product of the thermodynamic factor and the anion transference number v
Product of the thermodynamic factor and the anion transference number is given as
v= 1 +

∂ ln f±
∂ ln c2

1−t+ where t+ represents the measured conductivity associated with

Li+ mobility. For LiPF6 in PC/EC/DMC (10:27:63 by volume), it was found that t+ =0.399
and this value appeared to be independent of Li+ concentration [106].

For LiPF6 in

EC/DMC (2:1 by volume), the reported value was t+ = 0.363 [32]. Combining the em-

52

Figure 3.8: The Li+ diﬀusion coeﬃcient D2 as a function of the salt concentration and
temperature in LiPF6 in EC/DMC (2:1 by volume), calculated with Equation 3.14.

pirical result of the concentration and temperature dependence of from Ref. [106] and
t+ = 0.363 for LiPF6 in EC/DMC, Equation 3.15 was derived.
Figure 3.9 presents v as a function of Li+ concentration under diﬀerent temperatures. As
shown, v increases with the salt concentration but decreases with the temperature. Among
the two factors, the salt concentration has a much stronger eﬀect on v.

3.2.2

Viscoelastic Property of Separator

The viscoelastic behavior reﬂects the combined viscous and elastic response of the materials
under mechanical stress [96]. The stress-strain relation of a viscoelastic material is time and
temperature dependent. As a polymeric membrane, the separator exhibits the viscoelastic
property. The stress relaxation and thermal softening of the polymeric separator are consid-

53

Figure 3.9: The product of the thermodynamic factor and the anion transference number
v as a function of the salt concentration and temperature in LiPF6 in EC/DMC (2:1 by
volume), calculated with Equation 3.15.

ered in this work. At this stage, the linear viscoelastic material model is adopted to describe
the viscoelastic behavior of separator.
The polymeric separator has the anisotropic microstructure. Figure 3.10 presents a typical image of Celgard 2400 observed under SEM. As seen, the membrane has two distinct
material directions, which are commonly referred as machine direction (MD) and transcerse
direction (TD). In this work, the viscoelastic behavior of the separator was tested along MD.

3.2.2.1

Time Dependent Viscoelastic Behavior

The viscoelastic property of the polymeric separator can be characterized by the creep compliance and stress relaxation function. A creep test measures the strain with time under a
constant load. There was a series of tensile creep tests performed with a commercially avail-

54

Figure 3.10: SEM image of the surface of a single layer Celgard 2400 separator [9]. The
arrow indicates the transverse direction (TD).

able PP membrane, Celgard 2400, in the electrolyte solution of 1.1M LiPF6 in EC/DMC( 1:1
by volume) using a dynamic mechanical analyzer (DMA)[110]. The creep tests was applied
along machine direction (MD).
The measured creep compliance in the linear viscoelastic range was ﬁtted by a model
with Kelvin-Voigt elements by Sheidaei et al [110]. The Kelvin-Voigt model can be described
mathematically by a Prony series [118]

n
J(t) = J0 +

t
−τ
Ji e i

(3.17)

i=1

where J0 is the instantaneous compliance, Ji is the creep constant and τi is the retardation
55

times associated with Kelvin-Voigt element i.
The stress relaxation function E(t) is related to the creep compliance function, J(t), in
the Laplace domain by the following relation [112, 117, 118, 118]:

E(s)J(s) =

1
s2

(3.18)

Therefore, the Prony series for stress relaxation function is:

n
E(t) = E∞ +

t
Ei e− τ

(3.19)

i=1

The measured creep compliance was ﬁtted to Prony series using a non-linear least-square
subroutine in Matlab 2006, which was presented in Appendix A. The parameters in the stress
relaxation function expressed in Prony series were obtained by Laplace transformation using
Mathematica 8, which was presented in Appendix B. The results were presented in Table 3.3
and Figure 3.11. As the samples were tested under the tensile mode, the data represented
Young’s modulus. The stress relaxation function was implemented into the multiphysics
model.

3.2.2.2

Temperature Eﬀect

The temperature dependence of the separator’s mechanical properties is of great concern.
The time-temperature dependence of viscoelastic solids is modeled based on the framework
56

Figure 3.11: The stress relaxation function for the separator Celgard 2400 in 1.1 M LiPF6
in EC/DMC (1:1 by volume) [110].

Table 3.3: Prony series parameters for the stress relaxation modulus in 1.1M LiPF6 in
EC/DMC (1:1 by volume).
i

Stress relaxation constant (Pa)

Relaxation time (s)

E∞ = 1.23 × 108
1

E1 = 1.69 × 108

τ = 66.33

2

E2 = 1.22 × 107

τ = 968.07

3

E3 = 1.92 × 108

τ = 3860.3

4

E4 = 2.81 × 105

τ = 99756

57

of the time-temperature superposition principle (TTSP) of the linear viscoelasticity theory
[112, 118]. TTSP states that the relaxation or deformation behavior of a linear viscoelastic
material at an elevated temperature over a time period is equivalent to that occurs at a
lower temperature over a longer period time. Accordingly, the eﬀect of temperature can
be represented in a viscoelastic function through the use of a master curve and a reduced
time t = t/αT where αT is the so-called horizontal time-temperature shift factor. αT
is obtained through constructing a master curve at a reference temperature from a set of
curves measured under an isothermal condition over a range of temperatures. Considering
time-temperature superposition, the Prony series is adjusted as [114]:

4
E(t) = E∞ +

−a t τ
Ei e T

(3.20)

i=1

This dependence was measured for Celgard 2400 separator by conducting creep test at
diﬀerent temperatures in air [115]. It was then assumed that the temperature dependence
obtained in air is equivalent to that in electrolyte solutions. Electrolytes tend to evaporate.
During room temperature tests, the amount of electrolyte can be monitored and hence
maintained. This option is not possible for testing at elevated temperatures.
The Celgard 2400 in this test had a nominal thickness of 25 µm. The clamp had a gage
length ﬁxed at 15 mm. The samples in the form of long strips were cut from the separator
ﬁlm by a razor blade. The sample width varied from 4.88 mm to 5.20 mm.
The creep tests were carried out under a tensile mode using a DMA Q800 with a ﬁlm
clamp as shown in Figure 3.12(a). The isothermal tests were established by accommodating

58

(a)

(b)

(c)

Figure 3.12: The images depicting the sample installation and isothermal test establishment.
(a) sample installation; (b) & (c) the establishment of the isothermal test by closing the
furnace.

the clamp and ﬁlm in a closed furnace as shown in Figure 3.12(b) and (c). The creep
experiments were carried out at 25 ◦ C, 40 ◦ C, 45 ◦ C, 50 ◦ C, 55 ◦ C and 60 ◦ C.
The temperature range from 25 ◦ C to 60 ◦ C was the common range for a battery in
the working situation [116]. The heat generated from the electrical board below clamp
would dissipate into the furnace. The temperature in the furnace would gradually increase
if external cooling by connecting liquid nitrogen tank to the furnace was not provided.
The creep test was conducted by three steps. At the beginning, the temperature was
increased from the room temperature (23.5 ◦ C to 25 ◦ C) to the speciﬁed temperature by the
ramp rate 4C min1 . A preload force 0.001N was applied to maintain the sample at straight
and tight status. Then the sample was hold at the isothermal state for 5 min. At the last
step, a series of creep test were performed with samples in MD at diﬀerent load levels. A
new sample was used in each creep test. The samples were subjected to a constant load for

59

Figure 3.13: The creep compliance of Celgard 2400 at diﬀerent stress levels.

20 minutes.
For a linear viscoelastic material, the creep compliance is independent of the stress level
[110]. At the frist step, the stress range of the linear creep response needs to be determined.
Figure 3.13 presents the creep compliances measured at 1, 1.5, 1.75 and 2 MPa at 60 ◦ C. As
shown, the creep compliances at the 1 and 1.5 MPa were close to each other while the creep
compliances at the higher stress levels were obviously higher. The results indicate that the
creep response is linear at least up to 1.5 MPa. At 1.75 and 2 MPa, the creep response is
nonlinear. In further investigation, the creep tests were conducted under 1.5 MPa at diﬀerent
temperatures.
Figure 3.14 plots the creep compliance under 1.5 MPa at diﬀerent temperatures. The
creep compliance curves were converted to the stress relaxation curves by Laplace transformation [117, 118] and plotted versus logarithmic time in Figure 3.15.

60

Figure 3.14: The creep compliance of Celgard 2400 at diﬀerent temperatures: 25 ◦ C, 40
◦ C, 45 ◦ C, 50 ◦ C, 55 ◦ C and 60 ◦ C. The stress level is 1.5 MPa.

Figure 3.15: Construction of a master curve for PP using shear stress relaxation experimental
data. The stress plots are at 25 ◦ C, 40 ◦ C, 45 ◦ C, 50 ◦ C, 55 ◦ C and 60 ◦ C. The above stress
relaxation plots are obtained from creep compliance measured curves at the corresponding
temperatures by Laplace transformation.

61

Figure 3.16: The relationships of aT and temperature for PP.

The process of constructing a master curve is illustrated in Figure 3.15. First, a reference
temperature Tref has to be chosen. For the current case, Tref = 25 ◦ C. Next, the
curves measured at other temperatures are shifted horizontally to join the reference curve in
order to form a composite curve, or the so-called master curve. The amount of horizontal
shift required for these curves to form the master curve yields the αT function, as shown in
Figure 3.16. When plotted versus 1/T, aT is found to follow a Williams-Landel-Ferry (WLF)
relation for amorphous polymers, and Arrhenius relation for semi-crystalline polymers [117].
In Figure 3.16, in the range of 25 to 60 ◦ C, the aT function displayed a bilinear shape, i.e.
Arrhenius type of response. A change in slope corresponds to a change in the activation
energy which corresponds to a change in the dominant relaxation mechanism [113].
By ﬁtting shift factors at diﬀerent temperatures in Figure 3.16, the temperature dependent aT can be expressed as:
62

25 ◦ C - 50 ◦ C:

1
LogaT = 14230( ) − 47.76
T

(3.21)

1
LogaT = 7526.9( ) − 27
T

(3.22)

50 ◦ C - 60 ◦ C:

3.2.3

Battery Parameters and Other Properties

Table 3.4 provides a summary of the battery parameters used in this work, which include
the battery design parameters and the electrochemical related variables. The cross-section
of the cell is 24 cm2 and the total volume of the battery cell is 0.830 cm3 . The thickness of
each collector is 20 µm.
The eﬀect of the tortuosity of the composite electrode and the separator on the ionic
conductivity and diﬀusion coeﬃcient in the electrolyte was considered using the Bruggeman
relation by Equations 2.5 and 2.6, respectively. For the composite electrodes investigated in
this study, p = 3.3 [32]. For the separator, p = 3.3 is estimated with τt = εp−1 based on
the measured tortuosity value of τt = 6.0 [84, 86] and a porosity of ε = 0.46 [133].
The physical properties of the battery component materials used in the local heat generation model are listed in Table 3.5. The property parameters of the negative electrode,
the positive electrode, the separator, the Cu foil and the Al foil are transposed from Ref.

63

Table 3.4: Critical battery parameters used in the battery model.
Quantity
Thickness
L (µm)
Particle radius
R (µm)
Initial concentration
cini (mol m−3 )
Maximum concentration
cmax (mol m−3 )
Volume fraction of
solid phase ε1
Volume fraction of
liquid phase ε2
Electronic conductivity
of solid phase k1 (S m−1 )
Exponent p in
Bruggeman relation
Cross sectional area
Acell (cm2 )

Negative
electrode

Positive
electrode

Separator

Electrolyte

Cu
foil

Al
foil

100a

174a

52a

/

10b

10b

12.5a

8.5a

/

/

/

/

23751b
(SOC=0.9)

4572b
(SOC=0.2)

/

1000b

/

/

26390a

22860a

/

/

/

/

0.471a

0.362b

/

/

/

/

0.357a

0.444a

0.46c

/

/

/

100a

3.8a

/

/

/

/

3.3a

3.3a

3.3a

/

/

/

24a

24a

24a

/

24a

24a

a: Ref. [32]
b: estimated or assumed value.
c: Ref. [133]

64

Table 3.5: Physical properties of the battery components in “Thermal”sub-model.

Materials

Density
ρ
kg m−3

Heat capacity
Cp
(J (kg K)−1 )

Thermal
conductivity
Ktherm
(w (m K)−1 )

Electrical
resistivity
rc × 10−8
(Ω m)

Negative electrode

2223

641

6.5

/

Positive electrode

4202

672

6.2

/

Separator

900

1883

2.16

/

Electrolyte

1210

1518

0.099

/

Cu foil

8700

396

398

1.7

Al foil

2700

897

237

1.8

1861

845

4.21

/

2874

830.2

3.49

/

1043

1688

1.21

/

Electrolyte ﬁlled
negative electrode
Electrolyte ﬁlled
positive electrode
Electrolyte ﬁlled
separator

[119], with the exception of the electrolyte heat capacity. Both the experimental data and
calculated result of the heat capacity of EC and DMC were reported in Ref. [120]. Therefore
the electrolyte heat capacity in this work was calculated by the data from Ref. [121] based on
the rule of mixtures. The porous electrode can be treated as a composite of active particles
embedded in a matrix of the electrolyte [121, 122]. The densities, thermal conductivities,
and heat capacities of the two electrodes and separator with ﬁlled electrolyte were calculated
based on the rule of mixtures for composites as described in Ref. [121].
The mechanical properties of the diﬀerent battery components are listed in Table 3.6 and

65

Table 3.6: Mechanical properties of diﬀerent battery components used in the “Stress”submodel.
Parameter
Thermal expansion
coeﬃcient
α × 10−6 (K−1 )
Partial molar
volume
−6 (m3 mol−1 )
Ω × 10
Young’s modulus
E (GPa)
Poisson’s ratio
ν

Negative
electrode

Positive
electrode

Cu
foil

Al
foil

Separator

4.06a

8.62b

17b

23.6b

130b

3.5c

4.17c

/

/

/

12b

10b

117c

70c

section 3.2.2

0.3b

0.3b

0.35c

0.34c

a: Ref. [123]
b: Ref. [124]
c: Ref. [125]

Table 3.7 .
The elastic modulus of Liy C6 calculated in theoretical studies is in the range of 28-108
GPa [92, 126]. These values are high as compared with those reported in the data sheet
of graphite particles, which are in the range of 9-15 GPa [127]. This discrepancy is due
to the imperfection and porosity in commercial grade graphite particles. The present work
concerned a battery with commercially available graphite particles. A value of 12 GPa was
therefore selected.
The value of the Young’s modulus of graphite particles was reported in the range of 9-15
MPa [127]. Concerning a battery with the commercially available graphite particles, the
value of 12 GPa was selected in this work.
The partial molar volume for Lix C6 was calculated using Equation 2.25. The stoichiometric maximum concentration for Lix C6 was given as 2.64 × 104 mol m−1 [32] and 3.06 × 104
66

Table 3.7: Mechanical properties of electrolyte LiPF6 in EC/DMC (1:1 by volume)
Parameter

Value

Thermal expansion coeﬃcient α10×−6 (K−1 )

1000

Bulk modulus K (MPa)

4

Shear modulus G (MPa)

0.3

Branch Gi (i=1) (MPa)

0.1 (300s)

mol m−1 [128]. From C to Lix C6 , ∆y = 1. Assuming a maximum volume change of 11%
and the stoichiometric maximum concentration of 2.64 × 104 mol m−3 , the partial molar
volume of Lix C6 was estimated as:

Ω=

ei,c
3εij
∆ycs,max

≈

3 × 0.111/3
= 4.17 × 10−6 m3 mol−1
4
1 × 2.64 × 10

(3.23)

The partial volume of Liy Mn2 O4 was already estimated using the same method in Ref.
[15].
The viscoelastic behavior of the electrolyte was represented using the Prony series [118]
with an instantaneous shear modulus of 0.2 MPa. To avoid the convergence problem, the
value used in this work is about an order higher than that reported for the compliant gels
[129]. The assumed value of 0.2 MPa was more than two orders lower than the modulus of
the separator and, therefore, was judged as acceptable.

67

3.3

Simulation Condition

The multiphysics models for the lithium ion battery were constructed and analyzed using
software COMSOL® 3.5a. The sub-models of “Battery”and “Electrode”had a dimension
of a unit length whereas the “Thermal”and “Stress”sub-models had a real dimension. The
type, number and the maximum size of the elements in each sub-model is summarized in
Table 3.8.

Table 3.8: Element information in four sub-models of the multiphysics model.
Type of element

Dimensionless

Number of element

Max. size

Battery

1D linear

140

0.048

Electrode

2D triangle

11776

0.025

Thermal

1D linear

x

64

6.5µm

Stress

2D triangle

x

180522

1µm

The displacement ﬁeld across diﬀerent components was assumed to be continuous and
modeled with a continuous mesh. It should be noted that in a battery, the active particles
are in contact with, instead of bonded to, the separator. The relative movement between
the separator and electrodes may happen at a local area at microscopic level. To take into
consideration of this eﬀect requires a contact deﬁnition between the particles and separator.

68

Figure 3.17: Finite element meshes for (a) the “Battery”; (b) the “Electrode”; (c) the
“Thermal”and (d) RVE “Stress”sub-models.

It was possible to establish a contact relationship between the particles and separator in a
stand alone model. When incorporating the RVE “Stress”sub-models with contact deﬁnition
in the multi-scale model, convergence problems arose. The problem became worsen when a
compliant electrolyte phase was added. A continuous mesh, therefore, was adopted in the
RVE sub-models in this work. The stress revealed by this model, at best, may represent the
condition in the separator where the particle roughness created a nonslip condition in the
local area.
The models were solved using the “Direct (UMFPACK)”time dependent solver. The
relative tolerance was set as 0.0001 and the absolute tolerance was set as 0.001. UMPACK
is a robust, highly eﬃcient linear solver for non-symmetric systems. The average solution

69

times was within 2 hours using a single CPU.

70

Chapter 4
Heat Generation in a Li-Ion Battery
The preceding chapters have described the development of a multiphysics model to simulate
a basic Lix C6 |LiPF6 (EC/DMC)|Liy Mn2 O4 cell. The “Thermal”sub-model is capable of
simulating the heat generation and temperature rise in a Li-ion battery. Understanding
the heat generation is a necessary step in assessing the thermal stress in a battery cell.
Additionally, since the separator is a polymeric material, the temperature dependence of the
separator’s mechanical properties is of great concern. In this chapter, the “Thermal”submodel was validated against literature data. The predicted heat generation rates were within
the range of available experimental results. The model was used to evaluate the inﬂuence
of active particle size and component thickness on the heat generation rate and battery
performance.

4.1

Model Validation

The “Thermal”sub-model is validated by the comparison of the simulated data with the
experimental data from the similar battery system. Ye et al. [41] measured the temperature
rise of a Lix C6 |1.2M LiPF6 (PC/EC/DMC)|Liy Mn2 O4 prismatic battery using a thermocouple at the surface of the battery under a thermally insulated condition. The lumped heat
transfer coeﬃcient determined under this condition was h=0.38 W m−2 K−1 . Experiments
were performed under galvanostatic charge and discharges at 0.2C, 1C and 2C rates.

71

Ye’s experiments were simulated using the current model with the same geometrical
battery parameters for a similar Lix C6 |1.2M LiPF6 (EC/DMC)|Liy Mn2 O4 cell. Simulations
were performed considering the lumped heat transfer coeﬃcient at the surface of both two
current collectors. With h=0.38 W m−2 K−1 , the predicted temperature rises were lower
than Ye et al. [41]. Nevertheless, good agreements were obtained at all three rates with
h=0.2 W m−2 K−1 , as shown in Figure 4.1. This result provided a certain justiﬁcation for
the current model approach. It should be noted that the battery system investigated here
is not exactly the same as in Ye et al. [41]. Furthermore, the parameters in the current
model were based on literature data or selected from similar systems as described in section
3.2.3. There was no additional adjustment or ﬁtting as the model in Ref. [41]. Under
this condition, the agreement between the model predictions and the experimental results is
rather good. Therefore, in the following text, h=0.2 W m−2 K−1 was applied at the surface
of both two current collectors.

4.2

Thermal Analysis

The validated model was used in the investigation of the inﬂuence of individual physical and
battery design parameters in the heat generation of the battery. In simulations, the thermal
analysis was performed with an initial temperature of 25 ◦ C. The battery was operated
under galvanostatic conditions. The applied current for a 1C rate was 28 A m−2 . The
discharge cut-oﬀ voltage was 3.0 V. The discharge was stopped when the battery cell voltage
reached the cut-oﬀ voltage or when the reaction stopped because of the shortage of Li+
supply at the reaction interface [32].
Because the heat is generated in all battery components and the thin thickness of the

72

Figure 4.1: Comparison of predicted temperature rises with the experimental results of Ye
et al. [41].

cell, the temperature distribution is uniform through the thickness of a basic cell. This
was observed even at the highest charge rate. Thus only the temporal variation of the
temperature will be discussed.

4.3

The Eﬀect of Temperature Dependence of Physical
Properties

The inﬂuence of temperature dependence of individual properties on battery performance
was investigated by keeping one property as temperature independent at a time. Figure
4.2 presents the inﬂuence of each property on the predicted temperature rise under three
discharge rates: 1C, 4C and 8C. The results were compared with those from coupled and de-

73

coupled models, when all six sets of properties were with or without temperature dependence,
respectively.
At the same depths of discharge (DOD) level, the decoupled model predicted a higher
temperature rise than the other models. At 4C and 8C, the decoupled model also predicted
much lower achievable DOD. At 1C, only the temperature dependence of the Li+ diﬀusion
coeﬃcient D2 in electrolyte and reaction rate constant k0 for the positive electrode had a
noticeable eﬀect on the temperature rise. At 4C, the inﬂuence of the ionic conductivity k2 ,
and the Li diﬀusion coeﬃcient D1 and reaction rate constant k0 of the negative electrode
also became distinguishable. The inﬂuence of OCP was omitted in Figure 4.2 since its eﬀect
was negligible.
Kumaresan et al. [88] investigated the eﬀects of the temperature dependence of D2 ,
k2 , and ν for the electrolyte and D1 for the negative electrode on the cell capacity for
a Lix C6 |LiPF6 (PC/EC/DMC) |Liy CoO2 cell at 1C. Judging by the discharge proﬁles,
their results indicated that, with a starting temperature of 25 ◦ C, the eﬀect of temperature
dependence of k2 was negligible. These are in agreement with our results at 1C. However,
the current study indicated that the inﬂuences of temperature dependence of k2 became important with increasing rates. The same was true for some other properties. For properties
with a temperature dependence described by the Arrhenius equation, the one with a larger
value of activation energy exhibited a stronger inﬂuence, such as the k0 for the positive
electrode and D1 of the negative electrode.
The diﬀerence between the decoupled and coupled models under diﬀerent discharge rates
was then investigated. Figures 4.3(a) and 4.3(b) present the temperature rise and cell utilization at the end of discharge obtained for rates from 0.1C to 4C. The cell utilization is
deﬁned as the ratio of the actual achieved capacity to the maximum capacity. As shown, the
74

(a) 1C, electrolyte

(b) 1C, negative electrode

(c) 1C, positive electrode
Figure 4.2: Inﬂuence of the temperature dependent parameters on the predicted temperature rises during a galvanostatic discharge at 1C (a-c), 4C (d-f) and 8C (g-i). Simulations
were performed with one temperature dependent parameter at a time. The corresponding
temperature rises from coupled and decoupled model were also plotted in each ﬁgure.
75

Figure 4.2: (cont’d)

(d) 4C, electrolyte

(e) 4C, negative electrode

(f) 4C, positive electrode

76

Figure 4.2: (cont’d)

(g) 8C, electrolyte

(h) 8C, negative electrode

(i) 8C, positive electrode

77

(a) Temperature rise

(b) Cell utilization

Figure 4.3: (a) The temperature rise and (b) cell utilization as a function of discharge rate.
The cut oﬀ voltage is 3.0 V.

decoupled model predicted a much lower utilization at higher rates. The lower temperature
rise predicted by the decoupled model at rates greater than 4C was due to earlier termination
of the discharge process. In reality, a moderate temperature rise improves the ion transport
and reduces concentration polarization which leads to a better cell utilization. To consider
these eﬀects, it is necessary to use coupled thermal-electrochemical battery models. The
results presented thereafter are produced with the coupled model.

4.4

Heat Generation

In numerical simulations, the contributions of individual heat generation mechanisms can be
computed separately. Figure 4.4 plots the heat generation rates vs. DOD under 1C discharge
rate by individual mechanisms as well as their summation. The values from overpotential and
ohmic heats were all positive and among which the overpotential heat at the electrochemical

78

Figure 4.4: Partial heats due to diﬀerent sources contributing to the total heat during 1C
discharge.

interface and the ohmic heat in the electrolyte phase appear to be dominant. The amplitude
of the entropic heat was high, but it exhibited an S shape with DOD. Its shape resembles a
mirror image of the entropy change versus SOC curve of the Liy Mn2 O4 cathode in Figure
3.5(b). Due to the entropic heat, the total heat generation rate also exhibited an S shape. Sshaped heat generation curves have been observed in experimental measurements [23, 61, 68–
71] and reported in model predictions [23, 61, 68, 69, 71] of Li-ion batteries with diﬀerent
electrode materials.
Figure 4.5 plots the total heat generation rate vs. DOD as a function of the discharge
rate. Both the magnitude of the heat generation rate and its undulation increased with the
discharge rate. At 6C, the heat generation rate had a range from 227 W L−1 to 400 W
L−1 . For example, at 6C, the maximum achievable DOD was about 20%. The battery

79

performance at higher rates will be further discussed in section .
Several groups have measured the heat generation rates for Li-ion batteries [9,16-30]. As
discussed in the review paper of Bandhauer et al. [61], due to the diﬃculty in measuring
transient high rate heat generation, most of the experiments were conducted at lower rates
(1C) during discharge [16-28] or charge [16,18,19,21-30].
Table 4.1 compares the predicted minimum and maximum generation rates with the
available experimental data in the literature under discharge rates of 0.1C, 0.5C, and 1C.
The same comparison is provided in Figure 4.6. Because of the variations in battery materials
and uncertainties in battery geometries [61], the reported values varied widely. The purpose
of comparison in Table 4.1 and Figure 4.6 is to examine whether the prediction of this study
is within the same order of magnitude as the literature data. The results show that both the
maximum and the minimum heat generation rates at all three rates determined in this work
were within the range of the literature data. The predicted maximum heat generation rates
were within the range but appeared to be high compared to most of the experimental data
point. This might be caused by the range of DOD in the experiments. Figure 4.4 shows that
the maximum heat generation occurs at the end of discharge. Therefore, experiments with
a smaller DOD will yield a lower maximum heat generation rate.

4.5

Eﬀects of Battery Design Parameters

The battery cells with diﬀerent design parameters were designated as the following. The
baseline model with the set of parameters listed in Table 3.4 was referred to as 1R/1L,
where R and L stand for the active particle size and component thickness, respectively. In
other models, the radii of the positive and negative particles and/or the thicknesses of the

80

(a) 0.1C-2C

(b) 0.1C-6C
Figure 4.5: Heat generation rate versus DOD for discharge rates from 0.1C to 6C: (a) 0.1C2C; (b) 0.1C-6C.

81

Table 4.1: Comparison of the predicted minimum and maximum heat generation rates with
the available experimental data in literature.

Heat Generation Rate
q (W L−1 )

Discharge Rate
0.1C

0.5C

1C

Minimum values from
the experimental data

-0.45 to 0.79

-0.32 to 4.12

-0.05 to 13.3

Predicted minimum value

-0.4

1.273

8.67

Maximum values from
the experimental data

0.03 to 11

1.77 to 44

0.31 to 84.5

Predicted maximum value

1.79

13.9

37.0

References

[68–73, 78]

[65, 66, 69, 70, 80]

[65, 66, 68–71, 73, 76]

(a) 0.1C-2C

(b) 0.1C-6C

Figure 4.6: Comparison of the predicted heat generation rates with the upper and lower
bounds of experimental data in literature (a) the maximum heat generation rate; (b) the
minimum heat generation rate.

82

battery components were changed proportionally while the porosity and thickness of the
current collectors remained the same. A list of the models with diﬀerent particle size and
component thickness used in the design parameter study is given in Table 4.2.

Table 4.2: Models with diﬀerent particle sizes R and component thicknesses L investigated
in the design parameter study
Model
Parameter

Negative
electrode
L
R
(µm) (µm)

Positive
electrode
L
R
(µm) (µm)

L
(µm)

Cu
foil
L
(µm)

Al
foil
L
(µm)

Separator

0.1R/1L

1.25

174

0.85

52

10

10

1R/1L

100

12.5

174

8.5

52

10

10

1.5R/1L

100

18.75

174

12.75

52

10

10

1R/0.5L

50

12.5

87

8.5

26

10

10

0.1R/0.5L

4.5.1

100

50

1.25

87

0.85

26

10

10

Particle Size

Lu et al. [131] reported that reducing the particle size of the Liy Mn2 O4 positive electrode
improved the battery capacity. Botte et al. [132] used a mathematical model to investigate
the eﬀect of particle size in the negative electrode on temperature rise in a Lix C6 |LiPF6
(EC/DMC)|Liy NiO2 cell. The eﬀect of particle size has been recognized in the intercalation

83

(a)

(b)

Figure 4.7: (a) The cell temperature rise and (b) utilization at the end of the discharge in
the 0.1R/1L, 1R/1L and 1.5R/1L models.

stress model [12].
The eﬀect of the electrode particle size was investigated with the 0.1R/1L, 1R/1L, and
1.5R/1L models. Figure 4.7(a) presents the temperature rise at the end of discharge as a
function of the discharge rate. As shown, the inﬂuence of particle size on temperature rise
was small at discharge rates below 2C. Above 2C, the temperature rise apparently decreased
with reducing particle size. A further inspection revealed that the lower temperature rise
was related to the low battery utilization. As shown in Figure 4.7(b), the 0.1R/1L model
exhibited the highest utilization at the lower and higher rates. In the range of 2C to 4C,
however, the trend reversed. In other words, while a smaller particle size in general reduced
the temperature rise and improved the cell utilization, the trend was not monotonic across
the rates.
To understand this phenomenon, the concentration proﬁles of diﬀerent particle sizes were
examined. Figures 4.8(a) through 4.8(e) present the concentration proﬁles in the electrolyte
84

and at the surface of the active particles at three discharge rates: 1C, 3C and 5C. The comparison was made for DOD=0.91 (1C), DOD=0.55 (3C), and DOD=0.24 (5C) respectively,
when the cell voltage of the model with the lowest value reached the cut-oﬀ value of 3.0 V.
At 1C, the rate limiting mechanism was Li diﬀusion in the negative particles. As shown
in Figure 4.8(b), with a large particle size of 1.5R, Li at the surface of the particles was
depleted across the negative electrode. With the same volume fraction, reducing the particle
size increases the surface area and reduces the local current density [48]. A smaller particle
size also reduces the distance of Li diﬀusion from the particle interior to the surface and
sustains the supply of Li at the surface.
At 3C, the rate limiting mechanism switched to Li+ transport in the electrolyte. As
shown in Figure 4.8(c), Li+ depletion in the electrolyte in the positive electrode started
from the end next to the current collector. Because of a greater concentration gradient
earlier in the 0.1R/1L cell, Li+ depletion occurred earlier in this cell.
At 5C, Li+ transport in the electrolyte remained to be the rate limiting mechanism,
Figure 4.8(e). However, Li+ depletion occurred ﬁrst in the 1.5R/1L cell at DOD=0.24. At
this DOD level, 1.5R/1L exhibited the highest Li+ concentration gradient in the electrolyte.
Apparently, the eﬀect of particle size on Li+ concentration gradient varied during the
discharge. This phenomenon was further investigated. Figure 4.9 compares the development
of the Li+ concentration proﬁle during a discharge at 1C between the 1R/1L and 0.1R/1L
models at four DOD levels: 0.026, 0.28, 0.53, and 0.91. As shown, at the two lower DOD
levels, 0.1R/1L model had a lower Li+ gradient. The trend changed with increasing DOD.
At the two higher DOD levels, 0.1R/1L model had a greater Li+ gradient. This explains the
non-monotonic eﬀect of the particle size. The maximum DOD that a cell can reach decreases
with the discharge rate. At lower discharge rates, the cell can reach a higher maximum DOD
85

(a)

(b)

(c)

(d)

(e)

(f)

Figure 4.8: Comparison of the Li+ concentration proﬁles in the 0.1R/1L, 1R/1L and
1.5R/1L models under three discharge rates. (a, b) at DOD=0.91, 1C; (c, d) at DOD=0.55,
3C; (e, f) 5C; when one of the cells reached the cut-oﬀ voltage. (a, c, e) Li+ concentration
in the electrolyte, (b, d, f) Li concentration at the surface of the active particles.
86

(a)

(b)

(c)

(d)

Figure 4.9: The evolution of Li+ concentration proﬁles in electrolyte in the 0.1R/1L and
1R/1L models at 1C discharge rate; (a) at DOD=0.026 (100 s), (b) DOD=0.28 (1000 s), (c)
DOD=0.53 (2000 s), (d) DOD=0.91 (3400 s).

when cells with a smaller particle size exhibit a greater Li+ gradient. At higher discharge
rates, the maximum DOD of the cell is much lower. The cells with a smaller particle size
have a smaller Li+ gradient at lower DODs and hence exhibit a higher utilization than the
cells with a larger particle size.

87

4.5.2

Component Thickness

The eﬀect of the component thickness was examined using the 0.1R/0.5L, 0.1R/1L, 1R/0.5L
and 1R/1L models. Figure 4.10 compares the temperature rises and utilizations of these four
cells. With the same particle size, cells with thinner components resulted in a much lower
temperature rise. Reducing the component thickness reduced the concentration polarization.
The apparently higher temperature rise of 0.1R/0.5L and 1R/0.5L cells at discharge rates
greater than 4C was due to their much higher utilization rate.
It should be noted that the result presented here is at the cell level. At the battery
level, for the same total power, a battery constructed with thinner basic cells would require
more cells. The added weight and volume of the additional current collector layers must be
considered for the overall optimization of the battery.

4.5.3

Component Thickness and Particle Size Eﬀects

As discussed in Section 4.5.2, the inﬂuence of particle size on battery utilization is not
monotonic across the rates. The same trend was observed in models with a thin component.
Figures 4.11(a) and 4.11(b) compare the temperature rise and cell utilization for 0.1R/0.5L
and 1R/0.5L models up to a discharge rate of 18C. The 0.1R/0.5L model had a higher
utilization up to a discharge rate of 8C. It was overtaken by 1R/0.5L in the range of 8C to
15C. The trend was reversed back at above 15C.
Figures 4.12(a) through 4.12(c) compare the discharge proﬁles of 0.1R/0.5L, 0.1R/1L,
1R/0.5L and 1R/1L cells at 1C, 4C and 10C. The results clearly indicated that, with the
same particle size, reducing the component thickness is a rather eﬃcient method to improve
the cell utilization. On the other hand, the eﬀect of particle size is not monotonic.

88

(a)

(b)

(c)

(d)

Figure 4.10: The temperature rise and cell utilization at the end of the discharge. (a)
temperature rise, 1R/0.5L and 1R/1L cells; (b) temperature rise, 0.1R/0.5L and 0.1R/1L
cells; (c) cell utilization, 1R/0.5L and 1R/1L cells; (d) cell utilization, 0.1R/0.5L and 0.1R/1L
cells.

89

(a)

(b)

Figure 4.11: The temperature rise and cell utilization at the end of the discharge of 0.1R/0.5L
and 1R/0.5L. (a) temperature rise; (b) cell utilization.

4.5.4

Tortuosity

The tortuosity of a porous material may be aﬀected by material design and manufacturing
processes [9, 135]. To investigate the eﬀect of tortuosity, batteries with Bruggeman’s exponents of p=3.3 and p=2 were compared. Table 4.3 presents the eﬀective ionic conductivities
and diﬀusion coeﬃcients in all porous components in these two batteries. The porosity of
each component was the same as in Table 3.4.
The eﬀect of active particle size on temperature rise and cell utilization was examined
in Figure 4.13 for the batteries with p=3.3 and p=2. The impact of tortuosity was much
more signiﬁcant as compared to that of the particle size. Reducing tortuosity reduced the
temperature rise and improved the cell utilization. With p=2, the temperature rises at 2C
discharge rate were about 51% and 33% of those for p=3.3 batteries for the 1R/1L and
0.1R/1L models, respectively. The cell utilization was about 91% at 6C even for the 1R/1L
model when p=2.
90

(a)

(b)

(c)
Figure 4.12: Comparison of the discharge proﬁles between cells with a thickness of 0.5L and
1L. Discharge rate: (a) 1C, (b) 4C, and (c) 10C.

91

Table 4.3: The description of exponent p, tortuosity τt and eﬀective electrolyte transport
properties for battery I and battery II respectively.
Eﬀective Li+ diﬀusion
coeﬃcient in electrolyte
Def f

Exponent
p

Tortuosity
τt

Eﬀective ionic
conductivity
kef f

Battery I

3.3

10.7
6.0
6.5

0.033k2 (Lix C6 )
0.077k2 (PP)
0.069k2 (Liy Mn2 O4 )

0.033D2 (Lix C6 )
0.077D2 (PP)
0.069D2 (Liy Mn2 O4 )

Battery II

2

2.8
2.2
2.3

0.13k2 (Lix C6 )
0.21k2 (PP)
0.2k2 (Liy Mn2 O4 )

0.13D2 (Lix C6 )
0.21D2 (PP)
0.2D2 (Liy Mn2 O4 )

(a)

(b)

Figure 4.13: The temperature rise and cell utilization at the end of the discharge from the
battery models with diﬀerent particle size: (a) temperature rise; (b) cell utilization. The
models are simulated with p values 3.3 and 2, respectively.

92

(a)

(b)

Figure 4.14: The temperature rise and cell utilization from the battery models with two
battery thicknesses: (a) temperature rise; (b) cell utilization. The models are simulated
with a Bruggeman exponent of p=3.3 and p=2.

Figure 4.14 compares the eﬀect of the component thickness. With p=2, the battery with
1L thickness actually had a better utilization at high discharge rates compared to the 0.5L
battery. The improved utilization of p=2 batteries is attributed to reduced concentration
polarization, as shown in Figure 4.15. The gradient of Li+ in electrolyte in Figure 4.15(c)
was much smaller as compared to that in the p=3.3 battery in Figure 4.15(a). The rate
limiting factor for p=2 batteries was the Li depletion at the surface of the negative particles,
as opposed to the Li+ depletion at the positive electrode/current collector interface for p=3.3
batteries.
The above case studies demonstrate the complexity in battery design. There are competing rate limiting mechanisms. Due to the close couplings among thermal, transport,
and electrochemical processes, the controlling mechanism can switch from one to another
depending upon the battery operating condition and battery design. The inﬂuence of one de-

93

(a)

(b)

(c)

(d)

Figure 4.15: Comparison of the Li+ concentration proﬁles of the 1R/0.5L and 1R/1L cells
under a discharge rate of 3C for the p=2 battery and the p=3.3 battery respectively. (a) Li+
concentration in the electrolyte for the p=3.3 battery; (b) Li concentration at the surface of
the active particles for the p=3.3 battery; (c) Li+ concentration in the electrolyte for the
p=2 battery; (d) Li concentration at the surface of the active particles for the p=2 battery.
(a) and (b) were made at DOD=0.706, when 1R/1L ﬁrstly reaches the cut-oﬀ voltage 3V;
and (c) and (d) were made at DOD=0.908, when 1R/0.5L ﬁrstly reaches the cut-oﬀ voltage
3.0 V.

94

sign parameter, therefore, may be opposite under diﬀerent conditions. This presents a great
challenge to an experiment based design optimization. In this sense, mathematical and numerical models are advantageous. They are much more eﬃcient in examining a wide variety
of conditions and combinations and in identifying the underlying relationships and mechanisms. In Chapter 5, the “Thermal”sub-model was incorporated with the “Stress”sub-model
[125] to form a framework for the stress investigation in the separator.

4.6

Summary

The “Thermal”sub-model in the multiphysics model computes the heat generation and temperature rise for a Lix C6 |LiPF6 (EC/DMC)|Liy Mn2 O4 cell. The model prediction on temperature rise was in good agreement with experimental results of a similar battery at 1C
and 2C. The predicted heat generation rates were in good agreements with the experimental
data available in literature in the range of 0.1 to 1C.
The inﬂuence of the temperature dependence of the variables on temperature rises varied
with rates. Their importance was given in the following sequence: Li+ diﬀusion coeﬃcient
in electrolyte, reaction rate constant in the positive particles, Li diﬀusion coeﬃcient in the
negative particles, ionic conductivity, reaction rate constant in the negative particles.
Among various heating sources, the overpotential at the electrochemical interface and
the ohmic heat in electrolyte were the dominant sources. The magnitude of the heat generation rate and its undulation increased with the discharge rate. At discharge 6C, the heat
generation rate had a range from 227 W L−1 to 400 W L−1 .
The primary rate limiting mechanism was Li+ transport in the electrolyte. The thickness
of battery components had a great impact on concentration polarization. The batteries with

95

a thin component thickness had a lower temperature rise and better battery utilization. The
second rate limiting mechanism was Li diﬀusion in the negative particles. Reducing the
active particle size was most eﬃcient in improving the cell utilization when this mechanism
became dominant. The relationship between the particle size and cell utilization was not
monotonic across the discharge rates. The cell with a smaller particle size in general exhibited
a higher utilization but the cells with larger particle sizes had better utilization at certain
intermedium rates. This behavior was related to the evolution of the Li+ concentration
gradient with DOD.

96

Chapter 5
Stress in the Separator
In Li-ion batteries, stresses in the separator arise as the results of mechanical loading and
constrains, intercalation induced deformation in the active materials in the electrodes, and
thermal expansion mismatch between the battery components. This chapter focuses on the
simulation results about the stress and strain in the separator by the multiphysics model
including all three types of deformation. The Li-ion battery with the parameter 1R/0.5L was
chosen to investigate. The simulation results revealed that the stress in the separator varied
in phase with the battery cycles. Its state and magnitude depended upon the mechanical
property of the separator, electrode particle size and packing of the cell. Moreover, in
the separator, the thermal stress was comparable with the stress caused by intercalation
deformation of the electrode particle. Considering thermal eﬀect modiﬁed the stress curves
in the separator but the results remained in the same order.

5.1

Numerical Experiments

To begin with, the prototype quarter-sized and half-sized 2D RVE “Stress”sub-models were
built for the 1D battery LiC6 |PP|Liy Mn2 O4 battery in COMSOL® [32]. The original particle size and the component thickness are listed in Table 3.4. In its thickness direction, the
electrode in the quarter-sized sub-model had a dimension that was equivalent to a quarter of
the original electrode thickness. The thickness of the separator was about a half of that in

97

Table 3.4. Keeping the same particle size, the resulted quarter-sized “Stress”sub-model had
one layer of particles in the negative electrode and three layers of particles in the positive
electrode. The construction of the half-sized RVE followed the same rule.
In the multi-scale model, some of the parameters in the 1D battery sub-models need to
be modiﬁed with the change in battery dimension. Table 5.1 presents the battery component
thickness of the multi-scale model with a quarter-sized and a half-sized 2D RVE “Stress”submodels. The initial DOD in the negative electrode and positive electrode is 0.75 and 0.2,
respectively. All other parameters in the multi-scale model were kept the same as in Table
3.4 unless otherwise speciﬁed in the following text.

Table 5.1: Parameters in the 1D “Battery”sub-model with the original battery size and in
the multi-scale model with a quarter-sized or a half-sized 1D RVE “Stress”sub-model.

Parameter

Unit

“Stress”sub-model
Quarter-sized

Half-sized

Negative active particle radius

µm

12.5

12.5

Positive active particle radius

µm

8.5

8.5

Negative electrode thickness

µm

25

50

Positive electrode thickness

µm

51

85

Separator thickness

µm

26

26

Applied current density
(2C rate)

A m−2

11.7

23.4

98

The numerical experiments for the stress in separator were conducted with a dischargecharge cycle with a discharge rate 2C. The cycle consists of a discharge (0 s - 1000 s),
open circuit (1000 s - 1500 s), charge(1500 s - 2500 s), and another open circuit period
(2500 s - 3000 s). Figure 5.1(a) plots the current density history for a 2C discharge-charge
cycle for a half-sized cell. The history of cell voltage obtained by simulations using the 1D
“Battery”sub-model is presented in Figure 5.1(b).

(a)

(b)

Figure 5.1: (a) A battery cycle for a 1R/0.5L unit cell at a 2C discharge and charge rate. The
cycle consists of discharge (0 s-1000 s), open circuit period (1000 s-1500 s), charge(1500 s-2500
s), and open circuit period(2500 s-3000 s). (b) The cell voltage history for the dischargecharge cycle.

99

5.2
5.2.1

Numerical Investigations
The Eﬀect of Particle Packing

The initial numerical investigations were conducted with a prototype quarter-sized RVE
sub-model. The particles in this RVE had two packing patterns. As shown in Figure 5.2,
particles 2 and 4 are in contact with their neighbors. This packing pattern is referred as a
close packed pattern thereafter. On the other hand, particles 1, 3 and 5 are separated from
their neighbors by a space ﬁlled with electrolyte, which is referred as a loose-packed pattern
in the following text.

Figure 5.2: The predicted through thickness strain distribution at the end of discharge (1000
s) under a ﬁxed-free (10 psi) boundary condition using a quarter-sized 2D “Stress”sub-model.
The areas with a void like appearance have a strain value beyond the range and hence are
removed from the display. Particles 2 and 4 have a close-packed pattern whereas particles
1, 3 and 5 have a loose-packed pattern. The maximum strain in the separator is found near
a loose-packed particle of larger diameter.

100

Numerical simulations revealed that the stress and strain ﬁelds in the separator were
not uniform and varied in-phase with the discharge-charge cycle. For the discharge-charge
cycle studied in this work, a stress and strain maxima occurred at the end of the discharge
period at the surface in contact with the negative particles. At the side in contact with the
positive particles, it happened at the end of charge. When other conditions were equal, the
maximum strains and stresses were found near particles with a larger diameter.
Figure 5.2 displays the distribution of strain component along thickness direction at the
end of discharge. In the “Stress”sub-model, the strain in the electrolyte phase was much
higher as compared to that in other components. To ease the analysis of strain in the
separator, in Figure 5.2, the range of the display scale is set to ±0.05 and therefore, the
areas with a higher strain have a void like appearance.
Figure 5.2 reveals that the strain in the separator is higher in the area adjacent to a
particle with a loose-packed pattern. This is evident by comparing the strain distribution in
the vicinity of particle 1 with that of particle 2. The same trend was observed for strain along
in plane direction. It was suspected that this high stress might be due to the limitation of the
prototype RVE model that a continuous mesh is employed, instead of a contact deﬁnition,
between the particles and the separator.
FE models with a continuous mesh between components may be employed to study the
stress behavior of a “Stress”sub-model with close-packed particles. With a suﬃcient number
of particles, the central area of the “Stress”sub-model may reveal information about the stress
and strain history in a battery cell. For this purpose, half-sized 2D “Stress”sub-models were
investigated. The ﬁnal “Stress”sub-model consisted of a total 76 particles, comprised of 16
negative particles and 60 positive particles. Further investigations were carried out with this
RVE “Stress”sub-model.
101

The multiphysics model including the half sized 2D RVE “Stress”sub-model was employed
to do further investigation. In order to compare the importance of intercalation stress with
thermal stress in the separator, another two half-sized models, one considering the intercalation stress alone (labeled as intercalation only) and one considering the thermal stress alone
(labeled as thermal only), were also involved in this work. All the material properties and
boundary conditions were kept the same unless otherwise speciﬁed in the following text.

5.2.2

Concentration Change and Temperature Rise

In order to consider the strain caused by Li insertion/removal in a Li-ion battery, the Li
concentration distribution in the electrode is calculated in the “Stress”sub-model. Figure
5.3 depicts the Li concentration distribution in the electrode particles. As shown, Li concentration varied along the through thickness direction in the electrodes. The regions near
separator behaved larger concentration changes than those far away from the separator. It
demonstrates that the reaction rate is initially highest at this location, as the solid phase
electronic conductivity is higher than that of the electrolyte [48, 133]. In addition, there
was larger concentration gradient in the negative electrode than in the positive electrode, as
diﬀusion coeﬃcient D1 of the negative electrode is lower than that of the positive electrode.
Table 5.2 compares the Li concentration range from “Stress”and “Electrode”sub-models
at diﬀerent times. It demonstrated that the recalculated Li concentration from “Stress”submodel is at a similar level of that in “Electrode”sub-model. This recalculation is capable to
obtain the actual Li concentration distribution in the active particles under the current 2C
discharge-charge cycle .
The temperature rise in the Li-ion battery under the discharge-charge cycle was examined
by the “Thermal”sub-model. The temperature rise in the cell during the discharge-charge
102

(a)

(b)

(c)

(d)

(e)

(f)

Figure 5.3: The Li concentration distribution in the negative and positive electrode at
diﬀerent times during discharge-charge cycle: (a) 0 s, (b) 500 s, (c) 1000 s, (d) 1500 s, (e)
2000 s, and (f) 2500 s.

103

Table 5.2: The comparison of Li concentration range in “Electrode”and “Stress”sub-models
at the end of discharge (1000 s) and the end of charge (2500 s)
Electrode (mol m−3 )

Stress (mol m−3 )

Lix C6 at 1000 s

5659-13628

4540-14029

Liy Mn2 O4 at 1000 s

4540-13757

9440-14045

Lix C6 at 2500 s

16957-23610

16364-25354

Liy Mn2 O4 at 2500 s

4244-5526

3759-5912

cycle at 2C is plotted in Figure 5.4. As shown, the temperature rise during discharge is
about 8 ◦ C and during charge is about 7 ◦ C. The temperature rise decreased during the
open circuit periods (1000 s-1500 s and 2500 s-3000 s), as there was no heat generated. For
the entire discharge-charge cycle, the temperature rise is about 13 ◦ C. This temperature
proﬁle obtained from the “Thermsl ”sub-model was mapped to the “Stress”sub-model to
calculate the thermal stress in the Li-ion battery.

5.2.3

Thickness Variation

The dimensional change in the thickness direction of the battery is a measurable parameter.
To validate the model, the displacement in the thickness direction was investigated.

5.2.3.1

Particle Indentation

It was observed, upon applying a surface pressure, there was an immediate reduction in the
cell thickness. This dimensional change was absorbed by the deformation of the separator.

104

Figure 5.4: The temperature rise history in the Li-ion battery during discharge-charge cycle:
discharge (0 s-1000 s), open circuit period (1000 s-1500 s), charge(1500 s-2500 s), and open
circuit period(2500 s-3000 s).

The amount of deformation depended upon the elastic modulus of the separator and the
pressure at the surface of the cell. Figure 5.5 presents the displacement ﬁeld of the RVE
“Stress”sub-model at t = 0 s. The simulation was conducted with a separator under a
pressures of 10 psi. At t = 0 s, the volume of the active particles remained the same.
The deformation was, therefore, entirely due to the pressure. As seen in Figure 5.5, at
the microscopic level, the deformation was not uniform in the separator. The areas in
contact with the particles exhibited a larger deformation. Further evidence is provided by
the displacement proﬁles obtained at the surfaces of the separator, as displayed in Figure
5.6.
The Li intercalation/deintercalation in the active particles and battery temperature rise
result in the dimensional change of the battery cell and its components. By numerical
105

Figure 5.5: The through thickness displacement ﬁeld at t = 0 s predicted by a half-sized
2D RVE “Stress”sub-model with close-packed particles. The black lines in the separator
represent the position where the through thickness direction at the end of discharge was
checked in Figure 5.6.

106

simulations, the information on displacement, stress and strain components for individual
battery components at a speciﬁc location can be examined. This provides additional insights
into the fundamental mechanisms of stress in the separator in a Li-ion battery.
Figure 5.6(a), (b) and (c) present the through thickness displacement in the separator
at the end of discharge. The curves were plotted along a line which is in close proximity
and parallel to the interface with the negative electrode and positive electrode, respectively.
This clearly showed the indentation of particles into the separator. The periodical pattern in
displacement curve corresponded to the periodical packing pattern of the electrode particles.
The maximum displacement value appeared at the point where the particle is in contact
with the separator.
In the separator, as there was no active particles for Li diﬀusion, the value of the displacement caused by Li intercalation from the place near negative electrode and positive
electrode were similar in Figure 5.6(b). However, as shown in Figure 5.6(c), the value of the
displacement caused by thermal expansion near the positive electrode was 0.04 µm higher
than that near the negative electrode, because the thermal expansion coeﬃcient of the separator was 130 K−1 , which was much higher than the value of the other battery components.
As a result, the diﬀerence of the displacement involving two factors near negative electrode
and that near positive electrode, as shown in Figure5.6(a), was mainly caused by thermal
expansion of the separator.
The local indentation of particles at the surface of the separator had further implications.
The particles in the current model have a spherical shape. If irregular shaped particles are
employed, they may be pressed into the separator. Similar situations are when an irregular
shaped foreign particle is landed at the surface of the separator during manufacturing or
a dendrite forms at a particle surface near the separator. In these cases, a non-slipping
107

(a)

(b)

(c)
Figure 5.6: The indentations into the separators made by negative and positive electrode
particles. (a) caused by both Li intercalation and thermal expansion; (b) caused by Li
intercalation; (c) caused by thermal expansion.

108

interface may be created and a local tensile stress may arise.

5.2.3.2

Cell Thickness Variation

The variation of the cell thickness during the discharge-charge cycle was measured by the
through thickness displacement at the midpoint of the right side of the cell. As shown, for the
simulated case, the deformations due to intercalation only or thermal (thermal expansion)
only were in the same order. As in the previous investigation [125], the intercalation induced
deformation was cyclic. The maximum value of the deformation due to intercalation only
was about 0.04 µm in this work, about an order lower as compared to the previous results
[125] of 0.4 µm obtained under 1C rate. The deformation due to thermal expansion had the
same shape as the temperature history curve. The total maximum deformation was 0.1 µm
over a cell thickness of 181 µm.
The smaller thickness variation in the current simulation may be due to the viscoelastic
relaxation of the separator and/or a higher particle volume fraction (lower porosity) in
the positive electrode in the current simulation. To examine this eﬀect, batteries with
two diﬀerent volume fractions in the positive electrode were compared. Table 5.3 lists the
electrode porosity of two batteries. Battery 1 has the original electrode porosity as shown in
Table 5.1, while Battery 2 presents higher porosity in positive electrode. The other material
property and battery geometry were kept the same in both batteries. The results in Figure
5.8 showed that the diﬀerence between the two batteries was relatively small in the thermal
deformation but large in the interaction deformation. An increase in particle volume fraction
from 0.301 to 0.362 resulted in around 0.21 µm thickness variation increment.
A further inspection revealed that the greater thickness change is from the positive electrode thickness change. Figure 5.9(a), (c) and (e) present the through thickness change at
109

Figure 5.7: The displacement-history plots at the right hand end of the cell. Battery 1
(ε1,pos = 0.362) and Battery 2 (ε1,pos = 0.301) were investigated.

Table 5.3: The volume fractions of active materials and liquid phase in Battery 1 and 2.
Material property
volume fraction of active materials
ε1
volume fraction of liquid phase
(porosity) ε2

Negative electrode
Battery 1 Battery 2

Positive electrode
Battery 1 Battery 2

0.471

0.471

0.362

0.301

0.357

0.357

0.444

0.505

110

(a)

(b)

Figure 5.8: The displacement-history plots at the right hand end of the cell. (a) the thickness
caused by thermal expansion (b) the thickness caused by Li intercalation.

diﬀerent times. The diﬀerence of thickness change in the positive electrode between Battery 1 and 2 increased with discharge time. The diﬀerence between the two batteries was
negligible in the negative electrode.
The corresponding Li concentrations at the surface of the active particles at diﬀerent
times were also plotted in Figure 5.9(b), (d) and (f). In the Battery 2, the value at the positive
electrode was higher than that of Battery 1. This concentration diﬀerence was increased
with time. The higher surface concentration in the positive electrode would represent higher
average Li concentration for the electrode with same properties and active particle size. As
a result, the Li concentration diﬀerence ∆c1 in Battery 2 was larger than that in Battery
1. According to Equation 2.24, the large Li concentration diﬀerence induced the large strain
caused by Li intercalation, and ﬁnally caused the large thickness change in the positive
electrode of Battery 2.
The predicted maximum thickness variations for Battery 1 and Battery 2 was about
111

(a)

(b)

(c)

(d)

(e)

(f)

Figure 5.9: The through thickness displacement and corresponding Li particle surface concentration at diﬀerent discharge times: (a) through thickness displacement, 360 s; (b) Li
concentration, 360 s; (c) through thickness displacement, 720 s; (d) Li concentration, 720 s;
(e) through thickness displacement, 1000 s; (f) Li concentration, 1000 s. The data was from
the simulation that intercalation stress is considered alone.
112

0.060% and 0.18%, respectively. This value seems to be low as compared to the values obtained in an experimental study on LiC6 |LiPF6 |LiCoO2 battery [81]. The reported thickness
variation was 7.69 µm when averaged for a cell using natural graphite and 2.29 µm for using
mesophase carbon micro beads (MCMB). The cell thickness was 171 µm. This translates
into a linear strain of 4.49% for nature graphite and 1.34% for MCMB respectively.
It should be noted that the thickness variation predicted by the current model is due to
Li insertion/removal and thermal expansion. The volume expansion due to gas generation
in the Li-ion batteries [82] will also result in volume change. The volume variation due to Li
insertion/removal in cyclic, and that due to thermal expansion is mimic to the battery temperature rise curve; whereas the volume expansion due to gas generation increases gradually
over time.
The smaller dimensional change of the battery cell in the present study can also be
partially attributed to a better cell design. The maximum expansion of LiCoO2 positive
electrode is in the range of 1% [81] whereas the maximum volume change of Liy Mn2 O4
positive electrode is 6.5%. Because the volume variation of the two electrodes is opposite
in sign, it is the diﬀerence in volume change between the two electrodes, instead of the
absolute expansion of one electrode, determines the volume variation of the battery cell.
From this perspective, a LiC6 |Liy Mn2 O4 battery would have a smaller net volume variation
as compared to a LiC6 |LiCoO2 battery. Furthermore, in the present “Stress”sub-model,
Battery 1 with low porosity demonstrated a much lower thickness change than Battery 2.
This result suggests that the dimensional change of the battery may be reduced by optimizing
the electrode porosity in accordance with their volume expansion ratio. The stress analysis
in the following text was based on Battery 1.

113

5.2.4

Strain and Stress in the Separator

The stress state and its magnitude were found to vary greatly across the cell in both temporal
and spatial senses. As an example, Figure 5.10 presents the proﬁles of stress and strain in
the thickness direction at the end of discharge (t = 1000 s) across the cell thickness alone
a line here the separator is in contact with electrode particles at both side. As seen, there
are periodical hills and valleys in the proﬁles corresponding to the periodical location of the
particles. At 1000 s, the discharge had just stopped, and a Li concentration gradient still
existed in the particles. For the negative particles, the surface Li concentration was lower
than at the center. This resulted in a strain gradient in the particle. The stress gradient
was opposite in sign with the strain gradient. A maximum compressive intercalation stress
of 44 MPa was found at the center of the negative particles. The trend was opposite in the
positive particles. The center of the positive particle had a lower Li concentration than its
surface and the intercalation stress was therefore intensive.
The value of the strain and stress caused by Li diﬀusion and that caused by thermal
extension in diﬀerent battery components was also examined in Figure 5.10. In the electrode
particles, both the stress and strain caused by Li diﬀusion were about two orders higher than
the stress or strain caused by thermal extension. However, in the separator, the two kinds
of stress were comparable. They were around several mega pascals. The strain caused by
thermal extension was around 0.2%, which was higher than that caused by Li diﬀusion.
Figure 5.11 presents the distributions of the normal stresses and the shear stress in the
separator at t = 1000 s under a pressure of 10 psi. High compressive stresses can be seen in
the area in contact with the particles. The stress was higher at the side with larger diameter
of particles. As above-mentioned, due to diﬃculty with the contact modeling, a continuous

114

Figure 5.10: The through thickness stress and strain across the RVE ”Stress” sub-model at
the end of discharge.

mesh has been used between the particles and separator. The results only represent the
condition where the particle roughness created a nonslip condition in the local area.
When subjected to load, a material may yield due to shear stress or rupture due to
high tensile stress depending upon the stress state in the material, the available deformation
mechanisms in the material and other factors. Under 3D stress, the commonly used criteria
are the maximum distortion energy theory (i.e. Von Mises theory, Equation 5.1) and the
maximum principal stress theory (i.e. the theory of maximum tensile stress, Equation 5.2,
the Mohr and modiﬁed Mohrs theories etc.), respectively [134].

σV M =

(σ1 − σ2 )2 + (σ2 − σ3 )2 + (σ3 − σ1 )2 )
= σY
2

115

(5.1)

Figure 5.11: Stress components in the separator in the cell at the end of discharge. (a) σx ;
(b) σy ; (c) σxy .

Maximum(σ1 , σ2 , σ3 ) = St

(5.2)

where σ1 , σ2 , σ3 represent ﬁrst, second and third principal stress, respectively; σV M is
the Von Mises stress; σY is the yield stress; and St is the maximum tensile stress.
The variation of the stress during the discharge-charge cycles at a location of high stress
near a negative active particle and a positive active particle was inspected. Figure 5.12
presents a summary of stress history plots for the principal stresses and the Von Mises
stress. The left and right columns in Figure5.12 present the stresses from the place near
116

the negative electrode particles ans positive electrode particles, respectively. Figure 5.12
shows that the three principal stresses in the separator were in the negative range, i.e. under
compression, regardless it was near a negative or a positive particle. This result is consistent
with previous study [125]. The maximum Von Mises was 3.6 MPa, found at the side of
positive electrode. This value was lower than the measured strength of a value greater than
10 MPa for the separator at 25 ◦ C with the presence of electrolyte solutions [110]. The above
results indicate that under a normal operation condition with a closely packed particle bed,
the stress level in the separator is far from in danger of immediately yield or rupture in the
center of the cell.
Figure 5.12 also compares the stresses computed with intercalation only, thermal only,
and when both factors were in eﬀect. Unlike the cell thickness variation where the two eﬀects
were cumulative (Figure 5.7), the eﬀect of two factors on stress was not a simple summation
of the two individual cases. Considering thermal eﬀect modiﬁed the previously observed
interaction induced stress but the results remain in the same order.
The results in Figure 5.12 were obtained with the consideration of thermal-viscoelastic
relaxation of the separator. Neglecting this eﬀect, the estimated stress would be much higher
[125]. The eﬀect of the constitutive relationship on the stress was investigated in Figure 5.13.
Simulations were carried out with three types of constitutive behaviors for the separator:
linear elastic(E(const), Young’s modulus of E = 496 MPa, which is the stress relaxation
modulus at the beginning time), viscoelastic without temperature dependence (E(t), Figure
3.11), and thermal-viscoelastic (E(t, T), Equations 3.21 and 3.22). The thermal eﬀect was
considered in all three cases in the other three sub-models. It was found that the stress values
were at least twice higher when the separator was modeled as linear elastic as compared to the
two viscoelastic cases. Moreover, the thermal eﬀect on the stress relaxation barely reduced
117

(a)

(b)

(c)

(d)

Figure 5.12: The history plots of the principal stresses and the Von Mises stress in the
separator near the negative and positive electrode particles. (a) (c) (e) and (g): near negative
electrode; (b) (d) (f) and (h): near positive electrode.

118

Figure 5.12: (cont’d)

(e)

(f)

(g)

(h)

119

the amplitude of the cyclic stress.

5.2.5

The Eﬀect of Discharge Rate

For individual particles, increasing the charge/discharge rate would result in a greater Li
concentration gradient in the particle and therefore higher stresses. The eﬀect of charge
rate on stresses in the separator is not as unambiguous. The stress in the separator in a
Li-ion battery cell is the results of mechanical loading and constrains, intercalation induced
deformation in the active materials in the electrodes, and thermal expansion mismatch between the battery components. The value also depends upon the viscoelastic relaxation of
the separator itself.

Table 5.4: The range of Li concentration in sub-models at DOD = 0.3 at diﬀerent discharge
rates.
Discharge Rate

Electrode (mol m−3 )

Stress (mol m−3 )

2C

7710 − 17721

6604 − 18134

4C

5903 − 19084

2142 − 19325

6C

4636 − 19568

0 − 19683

To investigate the eﬀect of charge/discharge rate on stress in the separator, simulations
were carried out for discharge at 4C and 6C rates. As mentioned before, the Li concentration

120

(a)

(b)

(c)

(d)

Figure 5.13: The history plots of the principal stresses and the Von Mises stress in the
separator near the negative and positive electrode particles. (a) 1st principal stress; (b) 2nd
principal stress; (c) 3rd principal stress; (d) Von Mises stress. Three diﬀerent mechanical
properties of separator are considered: constant linear elasticity (E(const)); time dependent
viscoelasticity (E(t)); and time and temperature dependent viscoelsaticity (E(t,T)).

121

distribution in the active particles in “Stress”sub-model should be at a similar level of that in
the “Electrode”sub-model. However, it was observed that the discrepancy between the two
sub-models increased at higher rates. An example is shown for the depth of discharge (DOD)
= 0.3 in Table 5.4. The “Stress”sub-model tends to have a large concentration gradient as
compared to that in the “Electrode”sub-model. Moreover, the discrepancy also increased
with the DOD under the discharge/charge rate. This is caused by the diﬀerence in the volume
fraction of the active particles in the two sub-models. To have a consistent representation of
the micro structure in “Stress”sub-model requires a micro-structural resolved model, which
was described in Chapter 6. A qualitative discussion is given at here based on the observation
of the trend with the results of 2C and 4C from DOD = 0 to DOD = 0.7. When DOD was
higher than 0.7, the concentration in “Stress”sub-model was not at the same level of that in
“Electrode”sub-model.
Figures 5.14 and 5.15 compare the temperature rise and the cell thickness variation for
discharge at 2C and 4C. A higher discharge rate resulted in a higher temperature rise in the
cell. At DOD = 0.7, the temperature rise was 9.6 ◦ C and 18 ◦ C for 2C and 4C, respectively.
The higher cell temperature resulted in a larger change in the cell thickness. As shown in
Figure 5.15, the value increased about 50% at the rate of 4C.
The changes in stress values, however, were relatively small as shown in Figure 5.16.
This may be attributed to the combined eﬀect of a higher rate. At one hand, it reduces the
time for relaxation. On the other hand, a higher temperature rise at higher rates accelerates
the relaxation process. In summary, a higher discharge rate did not lead to a substantial
increase in stress in the separator.
The stress in the separator with the diﬀerent cell voltage range was also examined by
Figure 5.17. When DOD = 0.7, the cell voltage decreased to 3.64 V and 3.54 V under
122

Figure 5.14: The temperature rise Vs. DOD under discharge rate 2C and 4C. The DOD
range is 0 to 0.7.

Figure 5.15: The comparison of the thickness variations for the half-sized cell under discharge
rate 2C and 4C. The DOD range is 0 to 0.7.

123

(a)

(b)

(c)

(d)

Figure 5.16: The evolution of the stress components in the separator near the negative
particles under diﬀerent discharge rate: 2C and 4C: (a) the ﬁrst principal stress; (b) second
principal stress; (a) third principal stress and (d): Von Mises stress. The DOD range is from
0 to 0.7.

124

discharge rate 2C and 4C, respectively, as shown in Figure 5.17. There was no obvious stress
value change shown in Figure5.17 when the cell voltage boundary was decreased to 3.54V
under 4C.

Figure 5.17: The cell voltage for the half-sized cell under discharge rate 2C and 4C. The
DOD range is 0 to 0.7.

5.2.6

Anisotropy

The polymeric separator Celgard 2400 considered in this work has anisotropic mechanical
properties. A typical SEM image of Celgard 2400 is presented at Figure 3.10 in Chapter 3
[9]. The material displays a higher stiﬀness and strength along MD than TD. The separator

125

in the “Stress”sub-model is in a plane consisting of MD and the through thickness direction. Limited by the isotropic viscoelastic material model in COMSOL® the separator was
modeled in this work as an isotropic material with properties measured with MD specimens
tested in DMC.
When tested in air, separator exhibits strong anisotropy between MD and TD. Nevertheless, the diﬀerence in mechanical properties between the two directions diminishes when
tested in electrolyte solution [110]. The ratio of the elastic modulus between TD and MD in
1.1M LiPF6 EC/DMC is about 0.88. The ratio is 0.51 when tested in air. The mechanical
properties in the through thickness direction is not available currently. A close inspection of
Celgard 2400 suggests that the microstructure in through thickness direction has some similarity as in TD. Therefore, it is reasonable to assume Celgard 2400 as a transverse isotropic
solid, i.e., the property in the through thickness is the same as TD. Based on the above
information, the error induced by using an isotropic model is limited.

5.3

Summary

A multi-scale modeling approach for the stress analysis of the separator in a basic Li-ion
battery cell has been proposed. The stress analysis at the level of a battery cell was conducted
with a meso-scale RVE “Stress”sub-model of a basic cell. Currently, the model parameters
are available for a Lix C6 |LiPF6 |Liy Mn2 O4 cell. It is applicable for other insertion batteries.
The Simulation results provided new insights into the deformation of battery components
caused by Li insertion/removal in a charge and discharge cycle. It was observed that the net
cyclic dimensional variation of the battery cell was determined by the diﬀerence of the volume
change of the two electrodes and this amount may be regulated by optimizing the porosity

126

of the two electrodes. For the LiC6 |LiPF6 |Liy Mn2 O4 battery cell analyzed, the maximum
thickness variation caused by Li insertion/removal with diﬀerent electrode porosities was
about 0.2% and 0.05% respectively under a pressure of 10 psi.
The multiphysics model was used to analyze the stress in the separator. Simulations
were also conducted to investigate the inﬂuences of interaction and thermal phenomena and
their coupling, and the types of constitutive relationship of the separators on the stress. The
results showed that eﬀects of interaction and thermal are not simple summation and hence
must be considered concurrently. With an elastic behavior for the separator, the stresses in
the separator were about twice of those with a viscoelastic constitutive behavior.
The results indicate that under normal operation conditions the stress level in the separator is far from in danger of immediately yield or rupture in the center of the cell. It should
be noted that due to the diﬃculty in solving contact problems in COMSOL® , a continuous
mesh has been used between the particles and the separator. The results only represented
the condition where the particle roughness had created a nonslip condition in the local area
in a battery with a closely packed particle bed.

127

Chapter 6
Model Extension
The multiphysics model described in previous chapters is based on a 1D battery model. It is
suﬃcient when the eﬀect of the microstructure of the electrodes on the battery performance
is negligible. To consider the inﬂuence of microstructure, and particularly its evolution with
battery cycles, a 2D microstructure resolved model is desired. This chapter presents the
development and validation of a 2D microstructure resolved multiphysic model for a Li-ion
battery cell. The model is then used in the investigation of stresses in the active particles
and binders.

6.1

The Needs of a Microstructural Resolved Model

In the previous chapters, a multiphysics model for a Li-ion battery cell has been developed.
The model considered species and charge transport in electrodes, battery kinetics, heat
generation and heat transfer, intercalation and thermal induced stress as well as mechanical
stress. These phenomena were described with four sub-models that were coupled together.
The species and charge transport, battery kinetics, and Li diﬀusion in active particles were
described using the 1D battery model developed by Newman’s group [32, 48, 137], which
considers an averaged active particle size and porosity in the electrode. The intercalation,
thermal and mechanical stresses within and between battery components were analyzed using
a 2D RVE model with idealized microstructure. Within a certain range, the Li concentration

128

proﬁles were similar between the 1D “Battery”sub-model and “Stress”sub-model. However,
it was diﬃcult to maintain the consistency between the two sub-models. Furthermore, the
model cannot consider the eﬀect of real microstructure and their evolution during battery
cycle on the battery performance.
Figure 6.1 shows the cross-section of a graphite electrode of a commercial battery. As
illustrated, the microstructure of the electrodes is far from uniform. Particles are diﬀerent
shapes and sizes, as are the pores. The variation in local microstructures will inﬂuence the
species transport and therefore the local current density. The volume change of the active
particles with charge and discharge may also result in a change in local tortuosity. This eﬀect
would be substantial in electrodes with high capacity materials. To consider the inﬂuence
of microstructure, particularly its evolution with battery cycles, it is desirable to have a
microstructural resolved model.
The need for a microstructure resolved model has been discussed by diﬀerent groups.
Wang and Sastry [36] have developed the 3D microstructure resolved model constructed on
the porous cathode in Li|PEO2 − LiClO4 |LiMn2 O4 half cell. The modiﬁed Nernst-Planck
equation [48, 139] was involved to describe the mass transfer and ionic charge balance in the
electrolyte. Simulations with regular and random particle arrays and diﬀerent particle size
revealed that the electrode microstructure was critical to battery performance. Graham et
al. [140] employed the 3D particle resolved model of Li-ion batteries with regular particle
arrays to evaluate the approximations and to assist in the estimation of empirical Bruggeman
exponents in Newman’s 1D model. However, in the above works, the active materials were
treated as spherical particles. The irregular particle geometry was ignored.
In order to consider a more realistic electrode geometry, the 2D microstructure resolved
model developed by Smith et al. [37] was constructed based on the cross-section SEM image
129

Figure 6.1: (a) The cross-section image of a typical LiC6 -LiCoO2 rechargeable Li-ion battery
[37]. The upper, middle and bottom layer correspond to the graphite negative electrode,
polymeric separator and LiCoO2 positive electrode, respectively. (b) An individual slice
from a 3D X-ray tomography of graphite electrode [138].

of LiC6 |LiCoO2 . The investigation revealed that the Li diﬀusion in the active particles is
aﬀected by particle size, morphology, packing and surface roughness. Moreover, Less et al.
[141] built the 3D microstructure resolved model based on the FIB-SEM micrography of
the electrode. The result showed the battery performance can be optimized by varying the
geometric parameters of the electrodes.
To date, the microstructure resolved model has been employed to investigate the eﬀect
of microstructure on the battery performance. However, the microstructure eﬀect on the
stress in battery components has not been examined. The stress state and magnitude at
the local level may change largely depending on the electrode microstructure. In this work,
the 2D microstructure resolved model was developed to study the mechanical phenomenon

130

of the active particles and binders in the electrodes of a basic Lix C6 |PP|Liy Mn2 O4 cell.
This model was constructed based on the realistic electrode geometry. By using this model,
the eﬀect of complex electrode geometry on the mechanical phenomenon of the battery was
investigated.

6.2

Model Depiction

In this work, a 2D microstructure resolved model is built within the framework COMSOL® multiphysics. Figure 6.2 depicts the structure of the model. The entire microstructure resolved
model consists of two sub-models: a “Electrochemical”sub-model and a “Stress”sub-model.
The “Electrochemical”sub-model computes the spices and charge transport, as well as electrochemical reaction. The “Stress”sub-model computes the intercalation induced strain, the
mechanical strain and the stress. At this stage, the heat generation and thermal eﬀect is not
considered.

6.2.1

“Electrochemical”Sub-model

The battery kinetics,mass transfer and electronic balance are described in the “Electrochemcial” sub-model. Figures 6.3 and 6.4 illustrate the governing equations and boundary conditions, respectively.

6.2.1.1

Active Electrode Particles

The governing equations for Li diﬀusion and charge balance in the active particles in the
“Electrochemical”sub-model are the same as those used in Newman’s 1D “Battery ”model.
The Li diﬀusion in the solid particles follows Fick’s second law. The charge balance in the

131

Figure 6.2: The schematic of the 2D microstructure resolved model.

Figure 6.3: A schematic diagram of part of the Li-ion battery (one current collector, one
electrode and the separator) with the governing equations.

132

Figure 6.4: The boundary equations for the Electrochemical sub-model.

solid particles obeys Ohm’s law, as shown in Equation 2.9.
The Li ion ﬂux N0 and local current density jn are applied at the particle boundary. The
Bulter-Volmer equation is used to describe the relationship among the local electrochemical
reaction, concentration, and potential:

jn =k0

c2 c1,max − c1,surf c1,surf exp

ηαn F
RT

− exp

(−η)αp F
RT

(6.1)

The Li ion ﬂux has a linear relationship with the local current density jn , as shown in
Equation 2.8.

133

6.2.1.2

Electrolyte

The Nernst-Planck equations are used to describe the Li+ diﬀusion and charge balance in
the liquid phase[36, 137].
Based on the concentrated solution theory [137], the ﬂux of binary ionic species is described by:

d ln c0
N+ = − 1 −
d ln c

D2

i t
c + 2 + + c+ v0 = −D2 c +
z+ F

i2 t+
+ c+ v0
z+ F

(6.2)

d ln c0
N− = − 1 −
d ln c

D2

i t
c + 2 − + c− v0 = −D2 c +
z− F

i2 t−
+ c− v0
z− F

(6.3)

where N+ and N− are the negative and positive ion molar ﬂuxes, respectively; t+ and t−
are the transport number of cations and anions into which a molecule of electrolyte LiPF6
dissociates, respectively; z+ and z− are the charge number of the positive and negative ions,
respectively; c+ and c− are positive and negative ion concentration, respectively; and v0
is the solvent velocity. The magnitude of the dlnc0 is zero if the initial Li+ concentration
is constant. The value of v0 usually is assumed as 0. Therefore, equations 6.2 and 6.3 are
reduced to:

i t
N+ = −D2 c + 2 +
z+ F
134

(6.4)

i t
N− = −D2 c + 2 −
z− F

(6.5)

The electroneutrlity requires:

zi ci = 0

(6.6)

i

For the lithium salt, the charge number z+ = 1 and z− = −1. Therefore, Equation 6.6
is changed to:

c+ = c−

(6.7)

The current density can be expressed as:

2RT
i2 = −k 2 φ2 + k 2
cF

1+

135

d ln f±
d ln c

1 − t+

c

(6.8)

Substituting i2 with Equation 6.8, equations 6.4 and 6.5 become:

N+ = (−D2 +

2k2 t+ t− RT
d ln f±
) c+ −
)(1 +
d ln c
c+ z+ F 2

N− = (−D2 +

2k2 t− 2RT
d ln f±
)(1 +
) c− −
d ln c
c− z− F 2

k2 t+
φ2
z+ F

k2 t−
φ2
z− F

(6.9)

(6.10)

Because the charged mobile species in the electrolyte cannot diﬀuse into the current
collectors, the ﬂuxes of the two charged mobile species are set as zero at the interface of the
electrolyte and the current collectors, as shown in Figure 6.4. Finally, the periodic boundary
conditions are applied at the upper and lower boundaries.

6.2.1.3

Separator

The separator is treated as a homogeneous matrix, as the microstructure of the separator is
not explicitly modeled. Its porosity is considered through eﬀective properties for diﬀusion
coeﬃcient D2 and ionic conductivity k2 in the Nernst-Planck equations 6.9 and 6.10. The
relationships of D2,ef f , k2,ef f and D2 , k2 are explained in equations 2.5 and 2.6.
6.2.1.4

Binder

In the model, the conductive binders are depicted as a continuous thin layer over the particle
surface. The electrical charge balance in the binder obeys Ohm’s law as Equation 2.9. On
the other hand, the Nernst-Planck equation is extended over the domain occupied by binder.
136

In other words, binder is modeled as being conductive for electronic and ionic charge.

6.2.1.5

Current Collector

The electrons are the only mobile charge carrier in the current collectors. Equation 2.9 is
employed to determine the rate of electron transfer in the current collectors.
The voltage is set to zero at the left side of the left current collector in the battery. The
constant current density is applied at the right side of the right current collector. The model
simulates a galvanostatic condition.

6.2.1.6

Model Implementation

The current balance in the active particles, binders and current collectors was solved by
the Electric Currents Interface in COMSOL 4.2® ; the diﬀusion in the active particles was
solved by the Transport of Diluted Species Interface; the current balance and mass transport
in the electrolyte were solved by the Nernst-Planck Equations Interface.

6.2.2

“Stress”Sub-model

The “Stress”sub-model computes the stress and deformation in the active particles and
conductive binders. As described in Equation 2.34, the total strain is the summation of
diﬀusion-induced eigenstrain, thermal strain and the mechanical strain in the elastic battery
material. The Li concentration computed by the Electrochemical sub-model is used as input
to solve for diﬀusion-induced strain. The thermal strain is neglected at this stage.
The boundary conditions are depicted in Figure 3.3. They are the same as those for the
“Stress”sub-model in the multiphysics model in Section 3.1.3.

137

Figure 6.5: The boundary conditions for stress analysis.

The stress and deformation was solved by the Solid Mechanics (plane strain) interface
from COMSOL® 4.2.

6.3
6.3.1

Material Properties and Battery Parameters
Conductive Binder

The electric conductive performance of the binder in the electrode is controlled by binderparticle interaction [143], the ratio of conductive additive, and even the battery temperature.
The electrical conductivity value of 400 S m−1 is typical for PVDF material [143], which is
adopted in this work.
The mechanical properties of the diﬀerent binder types vary widely with temperature,
additives, or the preparation method [149–153]. Figure 6.6 presents the stress strain curves
138

of PVDF-HFP materials reproduced from literature data [150–152]. As seen, the response
of the materials varies greatly from very ﬂexible to have a moderate stiﬀness. In general, the
material is capable of large deformation with the exception of the grade with 5% Chitin.

Figure 6.6: The strain-stress curves for P(VDF-HFP) materials [150–152].

In order to comprehensively investigate the stress in the binder, two types of the PVDFHFP were selected: PVDF-HFP(FC1278) [152], a lower modulus material with E = 1.1
MPa; and PVDF-HFP−LiClO4 [150], a higher modulus material with a Young’s modulus
of E = 184 MPa.
The Prony series was employed to describe the viscoelastic behavior of PVDF-HFP
(FC1278) [152], which was labeled in Figure 6.6 by green color. By ﬁtting the stress re-

139

laxation curves calculated from the four element linear model for PVDF-HFP (FC1278)
[152], the Prony series parameters were determined, as shown in Table 6.1.

Table 6.1:
(FC1278).

Prony series parameters for the stress relaxation modulus of PVDF-HFP

i

Stress relaxation constant (Pa)

Relaxation time (s)

E∞ = 0.52 × 106
1

E1 = 0.045 × 106

τ = 100

2

E2 = 0.89 × 106

τ = 700

3

E3 = 0.0025 × 106

τ = 1500

4

E4 = 0.088 × 106

τ = 3000

5

E5 = 0.000207 × 106

τ = 8000

Figure 6.7 depicts the comparison of the simulated stress strain curve of the binder by
the stress relaxation modulus in Table 6.1 and the experimental cyclic stress strain curve
[152]. The cyclic stress strain curve was measured by a strain rate of 1% per minute. The
simulated curve showed good agreement with the original stress strain curve.
Figure 6.8 presents the stress strain curve of the PVDF-HFP-LiClO4 . The material has
an elastic modulus of 184MPa, a yield stress of 9.2 MPa and an ultimate tensile strength
of 13.4 MPa. PVDF-HFP-LiClO4 binder was described using a piece-wise linear plasticity

140

Figure 6.7: The comparison of the cyclic stress strain curves from simulation and experiment
for PVDF-HFP (FC1278) [152].

law, as shown in Figure 6.8.

6.3.2

Other Battery Components

The mechanical properties of the other battery components are identical to those presented
in the Chapter 3, as shown in tables 3.5, 3.6 and 3.7. In particular, the time dependent viscoelastic behavior of the separator is considered as expressed in Table 3.3. The temperature
eﬀect on the viscoelastic behavior is not involved.
Table 6.2 collects the other battery parameter used in the microstructure resolved model.
For consistency, most battery parameters, such as the battery thickness and the active
particle electrical conductivity are identical to the values in Table 3.4. The volume fraction
of the conductive binder is 0.080 and 0.071 in the negative electrode and positive electrode,

141

Figure 6.8: The stress strain curve of the PVDF-HFP−LiClO4 [150]
.

respectively. Because the microstructure resolved geometry shown in the Figure 6.2 already
reﬂects the electrode porosity and tortuosity, the Bruggeman relationship and exponent p
are no longer needed in the electrode domain.
The temperature dependence of the physical properties is not considered in the current
model. The values of OCP of the graphite and LiMn2 O4 , the ionic conductivity k2 and the
diﬀusion coeﬃcient D2 of Li in the electrolyte, etc., are those at 25 ◦ C.

6.4

Numerical Experiment

Figure 6.9 shows the FE mesh for the 2D microstructure resolved model. It includes 299996
triangular elements, and 2426258 degrees of freedom.
The solution was obtained using the “Direct (UMFPACK)”time dependent solver with

142

Table 6.2: Critical battery parameters used in the half-sized 2D microstructure resolved
model.
Quantity
Thickness
1 L (µ m)
2
Initial concentration
cini (mol m−3 )
Maximum concentration
cmax (mol m−3 )
Volume fraction of
active particles εa
Volume fraction of
conductive binder εc
Electrical conductivity
of solid phase k1 (S m−1 )

Negative
electrode

Positive
electrode

Separator

Electrolyte

Cu
foil

Al
foil

50

87

26

/

10b

10b

14870b

3900b

/

1000b

/

/

26390a

22860a

/

/

/

/

0.520b

0.430b

/

/

/

/

0.080b

0.071b

0.46c

/

/

/

100a

3.8a

/

/

/

/

a: Ref. [32]
b: estimated or assumed value.
c: Ref. [133]

a relative tolerance of 0.0001 and an absolute tolerance of 0.001 for all variables.

6.5

Model Validation

The “Electrochemical”sub-model in 2D microstructure resolved model was validated by comparing the simulated data with the experimental data from the same battery system. For
comparison, simulations were also carried out using Newman’s 1D battery model. Doyle et
al. [32] measured the cell voltage of the a Lix C6 |PP|Liy Mn2 O4 battery system, with the
polymeric electrolyte 2M LiPF6 in EC/DMC (the volume ratio 2:1) with the polymeric matrix PVDF-HFP. This experiment was simulated by the 2D microstructure resolved model
with the same battery system and geometry parameters. Figures 6.10 (a) and (b) show the
143

Figure 6.9: The mesh for the geometry of the 2D microstructure resolved model. Sections
(a) and (b) correspond to the zoomed-in regions. (a): mesh in the negative electrode; (b):
mesh in the positive electrode.

cell voltage comparison at discharge rate 0.5C and 1C, respectively. Good agreement was
obtained at both 0.5C and 1C.
The result from 2D microstructure resolved model had better agreement with the experiment result than that from 1D battery model at 1C. Whereas the 1D battery utilized
uniform spherical shapes for the active particles, the 2D microstructure model employed realistic microstructure geometry, which has various particle shapes and sizes. A recent study
revealed that the battery performance is aﬀected by the battery microstructure [37]. During
the discharge process, the graphite particles of smaller than average size, and of complex
morphology induce the localized lithium depletion, which leads to an increase in the polarization losses of the anode [37]. As a result, the 2D microstructured resolved model predicts
a lower cell voltage than the 1D Battery model, in which the active material is assumed as
144

(a)

(b)
Figure 6.10: The comparison of the predicted cell voltage with the experimental results [32]:
(a) 0.5C; (b) 1C.

spherical particles of the same size.

6.6

Results and Discussion

The validated model was used to investigate the stress in the binder and active particles in
the Li-ion battery. The numerical experiments for the stress in the binder were conducted

145

with a discharge-charge cycle with a rate 2C, which equates to an input current density 19.2
A m−2 . The cycle consists of a discharge (0 s - 1000 s), open circuit (1000 s - 1500 s), charge
(1500 s - 2500 s) and another open circuit period (2500 s - 3000 s). Two types of PVDF
were selected for investigation: PVDF-HFP(FC1278) [152], a lower modulus material with
E = 1.1 MPa; and PVDF-HFP-LiClO4 [150], a higher modulus material with a Young’s
modulus of E = 184 MPa.

6.6.1

Concentration Change

Figure 6.11 summarizes the Li distribution in negative and positive electrodes at multiple
points in time during the discharge and charge periods. Speciﬁcally, the points in time
are 0, 500, and 1000 seconds (corresponding to the discharge) and 1500, 2000, and 2500
seconds (corresponding to the charge). The simulation result showed that the concentration
distribution in the negative electrode was larger than that in the positive electrode. This
was attributed to the smaller Li diﬀusion coeﬃcient and larger particle size in the negative
electrode.
As shown in Figure 6.11, the mass transport in the smaller particles (i.e., region β)
was faster than that in the larger particles (i.e., region α). Furthermore, the particles with
complex morphology (i.e., region γ) also showed faster Li diﬀusion. This eﬀect of particle
size and morphology on the Li diﬀusion was also reported by Smith et al. [37]. In the same
electrode, the particles with small size or complex morphology result in a large surface area
per volume, enabling faster mass transfer rates [37]. The local Li diﬀusion was enhanced in
those edges and corners.

146

(a)

(b)

(c)

(d)

(e)

(f)

Figure 6.11: The Li concentration distribution in the negative and positive electrode at
diﬀerent times during discharge-charge cycle: (a) 0 s; (b) 500 s; (c) 1000 s; (d) 1500 s; (e)
2000 s; and (f) 2500 s. The unit for concentration in the graph is molm−3 .

147

6.6.2

Stress in the Electrode

6.6.2.1

Stress in the Active Particles

The stress distributions in the active particles of negative and positive electrodes at diﬀerent
time frames are plotted in Figure 6.12 and Figure 6.13, respectively. The particles were
surrounded by polymeric binders with low modulus. As seen, in the negative electrode,
the stress value at the particle boundary was higher than that at the particle center during
discharge, while conversely the stress at the boundary was lower than that at the particle
center during charge. The trend was opposite in the positive particles. This observation was
consistent with the result in Chapter 5. Furthermore, the high tensile stress was predicted
at the particle boundary and particle center in both electrodes. The stress value of a sharp
corner at the particle boundary may be even higher than that at the particle center.
The variation of the stress during discharge-charge cycles at the particle center (location
(2) in Figure 6.12 and location (3) in Figure 6.13) and the sharp corner with the high stress
(location (1) in Figure 6.12 and location (4) in Figure 6.13 ) was inspected in Figure 6.14. As
shown, the maximum stress value in the negative electrode was much higher than that in the
positive electrode. There are several factors resulting in the lower stress value in the positive
electrode active particles. According to the analytical solution [18], the radial or tangential
stress was proportional to the particle mechanical property, Li diﬀusion coeﬃcient, partial
volume ratio, particle size and local current density :

σ(r,orθ) ∝

j0 REΩ
3F D(1 − ν)

148

(6.11)

(a)

(b)

(c)

(d)

(e)

(f)

Figure 6.12: The ﬁrst principal stress distribution in the negative electrode at diﬀerent
times during discharge-charge cycle: (a) 0 s; (b) 500 s; (c) 1000 s; (d) 1500 s; (e) 2000 s; and
(f) 2500 s. The unit for stress in the graph is Pa.

149

(a)

(b)

(c)

(d)

(e)

(f)

Figure 6.13: The ﬁrst principal stress distribution in the positive electrode at diﬀerent times
during discharge-charge cycle: (a) 0 s; (b) 500 s; (c) 1000 s; (d) 1500 s; (e) 2000 s; and (f)
2500 s. The unit for stress in the graph is Pa.

150

In this work, the average positive electrode particle size was 8.26 µm, which was 0.6 time of
the negative particle size. Moreover, the negative electrode surface area was about one sixth
of the positive electrode surface area. As a result, the average local current density in the
negative electrode was about 6 times of that in the positive electrode. This large average
local current density was one of the critical reasons for the higher stress in the negative
particle.
The maximum tensile stress was found at the sharp corner of the negative particle boundary at the end of discharge. The value of the ﬁrst principal stress at location (1) was 470
MPa at 1000 s (the end of discharge), which was 4 times of the maximum stress value found
in the negative particle center at 2500 s (the end of charge). In the positive active particles,
the maximum stress value at the particle boundary was close to that at the particle center.
In order to reﬂect the eﬀect of the sharp corner on the stress value, the stress history
of a smooth surface point at the particle boundary was plotted at Figure 6.15. For a given
particle, the points located at the small radius boundary and the large radius boundary
were examined. The maximum ﬁrst principal stress value of a smooth surface point at the
negative particle boundary was approximately 200 MPa, which was less than one half of that
at the sharp corner.
The maximum value of Von Mises stress in the smooth surface boundary points of positive
electrode particles (location (7) and (8)) was around 85 MPa. It is higher than that in the
sharp corner (location (4)). In contrast to the shape of the particle at location (4), the
particles located at location (7) and (8) were long and narrow. This indicates that the Von
Mises stress at the positive electrode particle was not highly aﬀected by the sharp corner,
but by the particle shape itself. The sharp corner eﬀect was more obvious in the large active
151

(a)

(b)

(c)

(d)

Figure 6.14: The evolution of the stress components in the active particles. The stress
curves were plotted from diﬀerent locations: (a) negative electrode particle center, location
(2) in Figure 6.12; (b) negative electrode particle boundary, location (1) in Figure 6.12; (c)
positive electrode particle center, location (3) in Figure 6.13; (d): positive electrode particle
boundary, location (4) in Figure 6.13.

152

(a)

(b)

(c)

(d)

Figure 6.15: The evolution of the stress components at the normal point of particle boundary. The stress curves were plotted from diﬀerent locations: (a) negative electrode particle
boundary with small radius, location (5) in Figure 6.12; (b) negative electrode particle
boundary with large radius, location (6) in Figure 6.12; (c) positive electrode particle with
small radius, location (7) in Figure 6.13; (d) positive electrode particle boundary with large
radius, location (8) in Figure 6.13.

153

particles, for example, the negative electrode particles.

6.6.2.2

Stress in the Binder

Figure 6.16 presents the distribution of the ﬁrst principal stress in the two types of binders
at the end of discharge. Because of the thinness of the binder layer relative to the model,
one has to zoom-in to see the distribution. The zoomed-in area corresponds to the location
of maximum stress. The higher stresses were found at regions near sharp corners of the
particles, between the particles, or between the particles and current collector.

(a)

(b)

Figure 6.16: The contour of the ﬁrst principal stress in the binder corresponding to diﬀerent
binder materials at the end of discharge (1000 s): (a) PVDF-HFP(FC1278), E=1.1 MPa;
(b) PVDF-HFP-LiClO4 , E=184 MPa. The unit for stress in the graph is Pa.

The histories of the principal stresses and Von Mises stress in the binders PVDF-HFP(FC1278)

154

(low modulus) and PVDF-HFP-LiClO4 (high modulus) were plotted in Figure 6.17. The
stress varied with the cycle, reaching the maximum at the end of discharge, relaxing somewhat during the open circuit period and reducing to near zero at the end of charge. For the
battery system in this study, the binder stress was always higher in the negative electrode.
This may be due to the larger particle size in the negative electrode. The high modulus
binder had a higher stress magnitude. The maximum Von Mises stress was about 10 MPa
in the high modulus binder as compared to 0.18 MPa in the low modulus binder. The ratio
of the maximum stress between the two types of binders was about 56. The maximum value
of the ﬁrst principal stresses in PVDF-HFP-LiClO4 at the negative electrode was 15.2 MPa,
which was higher than the tensile strength 13 MPa shown in Figure 6.8. This result indicated
that in this battery system, PVDF-HFP-LiCl4 binder may fail during discharge.
Figure 6.18 plots the strain histories. The strain levels were about 20% and 7% for the
low modulus and high modulus binders, respectively. Since the material’s elongation tends
to reduce with improving in strength, binders with lower modulus appear to have advantages
over the higher modulus one in terms of reducing binder related failure, particularly for active
particles with large deformation. This result is consistent with the experimental observation
that ﬂexible binders improve the cycle life of silicon anode [154].

6.6.3

The eﬀect of Discharge Rate

In order to examine the eﬀect of discharge rate on the stress at the diﬀerent locations of the
active particles, the ﬁrst principal stress and Von Mises stress at both the particle center
and boundary under discharge rate 2C, 4C and 6C were compared in Figures 6.19 and 6.20,
respectively. As the ﬁgures indicate, the battery stopped discharging earlier under the higher
discharge rate as discussed in Chapter 4. These results also indicate that under the discharge
155

(a)

(b)

(c)

(d)

Figure 6.17: Stress history in binder: (a) binder with low modulus at a location in the
negative electrode shown in Figure 6.16 (a); (b) binder with low modulus at a location in
the positive electrode in Figure 6.16 (a); (c) binder with high modulus at a location in the
negative electrode in Figure 6.16 (b); (d) binder with high modulus at a location in the
positive electrode in Figure 6.16 (b).

156

(a)

(b)

(c)

(d)

Figure 6.18: Strain history in binder: (a) binder with low modulus at a location in the
negative electrode shown in Figure 6.16 (a); (b) binder with low modulus at a location in
the positive electrode in Figure 6.16 (a); (c) binder with high modulus at a location in the
negative electrode in Figure 6.16 (b); (d) binder with high modulus at a location in the
positive electrode in Figure 6.16 (b).

157

rate of 6C, the smallest stress magnitude can be attributed to the earliest determination of
the discharge process.
Under the same DOD, the ﬁrst principal stress magnitude was increasing with the increase
of the discharge rate. The maximum stress value appeared at the sharp corner of the negative
particle boundary, and reached 810 MPa at the end of discharge (DOD = 0.55) at 4C.
The eﬀect of discharge rate on Von Mises stress was most pronounced at the positive
particle boundary. The maximum Von Mises stress value at the positive particle boundary
was 630 MPa at the end of discharge of 4C, compared to 380 MPa at the end of discharge
of 2C. At the center of the negative and positive particles, the Von Mises stress value was
not obviously aﬀected by the discharge rate.
Several developed models [15, 24, 35, 157] for the stress in the active particles and the
measured results [155, 156] for the stress in the electrode have been established. The results
of the previous simulation indicated that the radial and tangential components of stress
are increasing with the discharge/charge rate but at the low rate range. After reaching a
peak of stress at a certain discharge/charge rate, the stress started decreasing at higher
rates due to the earlier battery termination [157]. In this work, the ﬁrst principal stress
under diﬀerent discharge rate showed the same trend. A wide range of the stress values
existed among previously simulated and experimental results even for the same electrode
material. Speciﬁcally, the maximum stress value ranged from 10 MPa to 400 MPa for the
stress in the graphite electrode [155, 156]. It should be noted that there were several factors
that may be responsible for the wide range of stress values observed. The ﬁrst factor was
attributed to the lack of experimentally measured property data for various constituents of
the composite electrode [155]. The second factor was the neglected complex geometry of the
active particle. For all the previously available models, the active particles were assumed to
158

(a)

(b)

(c)

(d)

Figure 6.19: The evolution of the ﬁrst principal stress in the active particles under diﬀerent
discharge rate: (a) negative electrode particle center, location (2) in Figure 6.12; (b) negative
electrode particle boundary, location (1) in Figure 6.12; (c) positive electrode particle center,
location (3) in Figure 6.13; (d) positive electrode particle boundary, location (4) in Figure
6.13.

159

(a)

(b)

(c)

(d)

Figure 6.20: The evolution of the Von Mises stress in the active particles under diﬀerent
discharge rate: (a) negative electrode particle center, location (2) in Figure 6.12; (b) negative
electrode particle boundary, location (1) in Figure 6.12; (c) positive electrode particle center,
location (3) in Figure 6.13; (d) positive electrode particle boundary, location (4) in Figure
6.13.

160

be with regular geometry, such as sphere or cylinder. In addition, the previous models only
considered the stress in the individual particles; while the eﬀect of the surrounding materials
on the active particle stress was largely ignored. In this work, the maximum Von Mises stress
value and the ﬁrst principal stress value at the negative electrode particle center were 95
MPa and 52 MPa at 6C, respectively, which were within the stress value range from previous
results. Furthermore, the maximum stress value at the sharp corner at the particle boundary
exceeded the predicted value by the previously developed models substantially. This reﬂects
the eﬀect of the localized particle’s complex geometry on the stress in the active particles.
The latter was impossible to examine if the previous models or experimental measurements
were employed.
The stress ﬁeld in an electrode may lead to the electrode morphological change [158].
The graphite anode suﬀers from particle cracking and damage resulting in surface structure
disordering upon prolonged cycling [159, 160]. As seen in Raman spectra, the graphite
electrode shows an increased intensity of the carbon D-band (1350 cm−1 ) with respect to
G-band (1580 cm−1 ) in a cyclic Li-ion battery [155, 161]. Such structure damage has a direct
eﬀect on the thickness and composition of the solid-electrolyte-interphase (SEI) layer [162,
163]. It is particularly true for the active particles under high intercalation/deintercalation
rate [159]. This work provides a quantiﬁable method for the characterization of the stress ﬁeld
and evolution in the electrode. It adds to the current state of knowledge on the relationship
of electrode mechanical damage and morphological change at microstructure level, especially
for the instances with high intercalation/deintercalation rate.
The eﬀect of discharge rate on the stress in the binder was examined next. Figures 6.21
and 6.22 compared the stress value at the locations showing high stress value at both PVDFHFP(FC1278) and PVDF-HFP-LiClO4 , respectively. The stress change in the binder under
161

diﬀerent discharge rates was not so obvious as in the active particles. Under the same DOD,
both the ﬁrst principal stress and Von Mises stress magnitudes increased slightly with the
discharge rate. The battery operation rate was not the dominant factor for the stress in the
binder.

6.7

Summary

A 2D microstructure resolved multiphysics model including electrochemical and mechanical
analysis has been developed. The realistic microstructure was represented in the model. It
was used to study the stress in the active particle and the polymeric binder of a Lix C6 |PP|Liy Mn2 O4
cell.
The stress in the negative electrode particle was higher than that in the positive electrode
particle. The larger particle size and higher average local current density were the critical
reasons. The eﬀect of the irregular particle shape on the stress in the electrode particles
was examined. The maximum stress value was at the sharp corner of the negative electrode
particle boundary. The value was much higher than that at the particle center or the smooth
surfacer point at the particle boundary.
In the polymeric binders, the maximum stress was found in regions near sharp corners
of the particles, between the particles, or between the particles and current collector. The
binder stress was higher in the negative electrode where the particle size was larger. The
maximum stress was observed at the end of discharge. The binder with a higher modulus
exhibited higher stresses. The stress in the binder was proportional to its modulus while the
strain remained at the order of 10-20%. Flexible binders appear to have advantages over the
more rigid binders in terms of reducing binder related failure. The eﬀect of the discharge

162

(a)

(b)

(c)

(d)

Figure 6.21: The evolution of the ﬁrst principal stress in the polymeric binder under diﬀerent
discharge rate. The stress curves were plotted from high stress location with diﬀerent binder
material: (a) PVDF-HFP(FC1278) (low modulus), location (1) in Figure 6.16(a); (b) PVDFHFP(FC1278) (low modulus), location (2) in Figure 6.16(a); (c) PVDF-HFP-LiClO4 (high
modulus), location (3) in Figure 6.16(b); (d) PVDF-HFP-LiClO4 (high modulus), location
(4) in Figure 6.16(b).

163

(a)

(b)

(c)

(d)

Figure 6.22: The evolution of the Von Mises stress in the polymeric binder under diﬀerent
discharge rate. The stress curves were plotted from high stress location with diﬀerent binder
material: (a) PVDF-HFP(FC1278) (low modulus), location (1) in Figure 6.16(a); (b) PVDFHFP(FC1278) (low modulus), location (2) in Figure 6.16(a); (c) PVDF-HFP-LiClO4 (high
modulus), location (3) in Figure 6.16(b); (d) PVDF-HFP-LiClO4 (high modulus), location
(4) in Figure 6.16(b).

164

rate on the stress in the binder was limited.

165

Chapter 7
Conclusion and Future Work

7.1

Summary and Conclusion

A multiphysics model for the stress analysis of the separator in a basic Li-ion battery cell
has been developed. Three groups of phenomena that occurred inside the battery cell were
considered in order to perform the stress analysis. They are electrochemical, thermal and
mechanical phenomena. Besides mechanical stress, the intercalation stress and thermal stress
in a separator were investigated in this study. This work is the ﬁrst work that computes the
stress in battery components within a multiphysics framework with ongoing electrochemical
reaction.
In this model, the multi-physics equations were solved concurrently using four coupled
sub-models at diﬀerent scales. The battery electrochemical kinetics, current balance, ionic
transport and mass transport were solved by a 1D macroscopic “Battery”sub-model and a
pseudo 2D “Electrode”sub-model based on an existing 1D battery model in COMSOL® ,
where the porous electrode theory and the concentrated solution theory were employed. The
thermal behavior was considered by the 1D “Thermal”sub-model. The local model was used
to simulate the heat generation in a Li-ion battery. The stress analysis at the level of a
battery cell was conducted with a 2D idealized meso-scale RVE “Stress”sub-model of the
battery cell. In order to consider the intercalation stress in the separator, the Li ﬂux in
the 1D “Battery”sub-model was mapped to the active particle surface to calculate the Li
166

concentration in the sub-model. The battery temperature obtained from the “Thermal”submodel was also mapped to the “Stress”for thermal stress calculation. The multiphysics
model was implemented in COMSOL® , a commercially available multiphysics code.
The major limitation of the model is the inability to get a converged solution for the
multiphysics problem when a contact condition between materials with a large diﬀerence
in stiﬀness was introduced, as in the case between the separator and active particles. To
overcome this obstacle, a continuum mesh had to be used at the separator and electrode
interfaces. The stress revealed by this model, at best, may represent the condition in the
separator where the particle roughness created a nonslip condition in the local area.
The developed multiphysics model was used to investigate the stress in polymeric separator PP in a Lix C6 |PP|Liy Mn2 O4 cell, with the liquid electrolyte 1M LiPF6 in a solvent
of EC/DMC (2:1 by volume). The measured thermo-viscoelastic property of the separator
was used in the work.
The “Thermal”sub-model was validated against available experimental data in terms
of cell temperature rise and heat generation rates and then used in the investigation of
relationships among the battery temperature rise, cell utilization and the design parameters.
The thickness of battery components had a great impact on Li+ concentration polarization.
The batteries with a thin component thickness had a lower temperature rise and better
battery utilization. Reducing the active particle size was most eﬃcient in improving the cell
utilization when Li diﬀusion in the negative particles became dominant.
The simulation results from “Stress”sub-model provided new insights into the battery
components deformation in a discharge/charge cycle. It was observed that the net cyclic
dimensional variation of the battery cell was determined by the diﬀerence of the volume
change of the two electrodes and this amount may be regulated by optimizing the porosity
167

of the two electrodes.
The electrode particle made small indention into the separator under pressure. The maximum stress in the separator was found in indented areas. The stress state and its magnitude
in the separator depended upon the particle size, particle packing, and the boundary constraint of the battery cell. The eﬀects of interaction and thermal are not simple summation
and hence must be considered concurrently. The results indicate that under normal operation conditions the stress level in the separator is far from in danger of immediately yield or
rupture in the center of the cell.
In summary, the multiphysics model developed in this study can be used to investigate
the stress in the separator, the heat generation in the battery, and the relationship among
electrochemical, thermal and mechanical phenomena. Due to the lack of the complete set
of model parameters and experimental data, the model was qualitatively rather than quantitatively validated. However, as the homogenized active particle size and porosity were
considered in the 1D “Battery”sub-model and pseudo 2D “Electrode”sub-model, the multiphysics cannot consider the eﬀect of electrode microstructure and its evolution on the battery
performance.
In order to involve the real microstructure eﬀect, a 2D microstructure resolved model was
developed in Chapter 6 for the same Lix C6 |PP|Liy Mn2 O4 battery system. This 2D model
was constructed based on the real electrode geometry. There were two sub-models sharing
the same geometry in the model: “Electrochemical”sub-model and “Stress”sub-model. The
“Electrochemical”sub-model computes the spices and charge transport, and electrochemical
reaction. The “Stress”sub-model computes the intercalation induced strain, the mechanical
strain and the stress. At this stage, the heat generation and thermal eﬀect is not considered.
The 2D microstructure model was used to investigate the stress in the active particles and
168

conductive binders.
By using the 2D microstructure resolved model, The eﬀect of the irregular particle shape
on the stress in the electrode particles was examined. The maximum stress value was at the
sharp corner of the negative electrode particle boundary. The value was much higher than
that at the particle center or the normal particle boundary point. Furthermore, the stress in
the negative electrode particle was larger than that in the positive electrode particle. The
larger particles size and average local current density were the critical reasons.
In the polymeric binders, the maximum stress was found in regions near sharp corners
of the particles, between the particles, or between the particles and current collector. The
binder with a higher modulus exhibited higher stresses. The stress in the binder was proportional to its modulus while the strain remained at the order of 10-20%. Flexible binders
appear to have advantages over the more rigid binders in terms of reducing binder related
failure.

7.2

Outlook of Future Work

This work developed a ﬁnite-element based multi-scale, multiphysics approach to investigate
the stress in a Li-ion cell. It is validated and used to predict the stress in the separator
in LiC6 |Li2 MnO4 battery system, as it is one of the systems with a complete set of model
parameters and experimental data.
The application of the current model can be expanded to battery systems with high
capacity electrode materials, which is one of the most active topics in Li-ion battery research.
Compared to graphite, Sn, Ge and Si have 10 times better volumetric capacity and 210 times better speciﬁc capacity. Employing these materials will signiﬁcantly improve the

169

energy density of the batteries. However, the cyclic stability of the electrodes with high
capacity materials in general is rather poor. A problem accompanied by high Li storage
capacity is the large volume changes during charging and discharging cycles. To improve the
cycle life, electrodes with various sizes, geometries and constructions have been investigated
experimentally. A model that can assess the battery performance and stress development in
systems with these electrodes will be invaluable. It will provide insight into the relationships
between diﬀerent parameters. It can also be used in design optimization.
A main challenge in multiphysics model development is to obtain the material properties
and parameters required for the multiphysics model for these systems.
Among all high capacity anode materials, Si has received most attention. Crystalline
Si was found to go through a series of phase changes and eventually becomes amorphous
Si when reacting with Li+ [165]. Shenoy et al. [164] and Zhao [158] employed the ﬁrst
principal model to calculate the SOC dependent elastic properties of Li-Si alloys. The
results indicated that the struture change of the lithiation of crystalline Si electrode renders
the Si anode based system the capability of large plastic deformation. Zhao et al. [166]
proposed a phenomenological model to describe the yield surface of Si with Li interaction.
The model prediction was in good agreement with the stress measurements of Si ﬁlm during
cyclic intercalation/deintercalation [167]. This model may be incorporated into the present
multiphysics model to investigate the battery with Si anode. Similar work is required for
other high capacity material systems.
Besides inﬂuencing the deformation behavior of the host, Li concentration [168] and
stress state [15] also inﬂuence the Li diﬀusion in the host material. Further experiments and
numerical investigations are required to consider these eﬀects in the multiphysics model for
high capacity electrodes.
170

Another challenge in modeling Li-Ion battery with high capacity electrodes is the numerical technique itself. As stated in Chapter 4, models built upon COMSOL® had convergence
problems when modeling contact between materials with large diﬀerence in stiﬀness and a
continuous mesh had to be used. With electrodes capable of large volume variation, the
issue with contact modeling needs to be resolved.
The current microstructure resolved model is limited to 2D. The local level porosity and
tortuosity in 2D electrode geometry are not the same as those in the realistic 3D geometry.
The current imaging techniques can reveal 3D microstructures to certain accuracy. The mesh
generation tools are available to automatically generate the mesh for 3D microstructure
geometry. However, the computer resource required to solve 3D models far exceeds the
current capability.

171

APPENDICES

172

Appendix A
Non-Linear Least-Square Subroutine
The following program is the non-linear least-square subroutine for creep compliance of Celgard 2400 at 40 ◦ C. The subroutine for creep compliance at diﬀerent times is the similar
but using diﬀerent original creep compliance data.
clear
x0 = [0.1;20;100;20;49;10];

% initial point, active set algorithm only

xdata=[2.04,4.04,10.04,101.05,308.06,640.08,921.99,1199.01];

% diﬀerent time

ydata=[1082.967,1120.411,1175.19,1386.174,1563.476,1734.761,1839.002,1919.932];
%creep value corresponding to diﬀerent time value
[xe, resnorm] = lsqcurveﬁt(@creep, x0, xdata, ydata);
T1=0.1;

% nonlinear curve-ﬁtting least-square

% retardation time

T2=100; T3=1000; T4=10000; T5=100000; j=xe(1)*xe(1)+xe (2)*xe (2).*(1-exp (-xdata/T1))
*100000+xe (3)*xe (3).*(1-exp (-xdata/T2))*150000+xe (4)*xe (4).*(1-exp (-xdata/T3))
+xe (5)*xe (5).*(1-exp (-xdata/T4))*14+xe (6)*xe (6).*(1-exp (-xdata/T5));
% Prony series expression for creep curve
plot (xdata,j,’+’)

% plot original creep curve

hold
plot (xdata,ydata,’*’)

% plot nonlinear curve-ﬁtting curve

legend (’Fitted Curve’,’DMA Data’);

% add legend

173

Appendix B
Laplace Transformation Subroutine
The following program is the Laplace transformation subroutine to obtain stress relaxation
function for Celgard 2400 at 40 ◦ C. The subroutine for stress relaxation function at diﬀerent
temperature is similar but using diﬀerent creep constant data.This program is written by
Mathematica 8.

J0 = 1027

* creep constant*

J1 = 87.85
J2 = 365.21
J3 = 0.821
J4 = 3850
J5 = 178
τ1 = 0.1

* retardation time*

τ2 = 100
τ3 = 1000
τ4 = 10000
τ5 = 100000

−t
−t
−t
−t
−t
J = J0 + J1 ∗ (1 − e τ1 ) + J2 ∗ (1 − e τ2 ) + J3 ∗ (1 − e τ3 ) + J4 ∗ (1 − E τ4 ) + J5 ∗ (1 − E τ5 )
* Prony series of creep compliance*
JS = LaplaceTransform[J, t, s]

* Laplace transformation of creep compliance*
174

GS = (1/s2 )/JS

* stress relaxation function in Laplace domain*

InverseLaplaceTransform[GS, s, t]

* inverse Laplace transformation*

175

BIBLIOGRAPHY

176

BIBLIOGRAPHY

[1] S. J. Harris, A. Timmons, D. R. Baker, and C. Monroe, “Direct in−situ Measurements
of Li Transport in Li-Ion Battery Negative Electrodes”, Chemical Physics Letters, 485,
2010, 265-274.
[2] S. S. Zhang, “A Review on the Separators of Liquid Electrolyte Li-Ion Batteries”,
Journal of Power Sources, 164(1), 2007, 351-364.
[3] www.Limn2o4.com
[4] V. Etacheri, R. Marom, R. Elazari, G. Salitra and D. Aurbach, “Challenges in the Development of Advanced Li-Ion Batteries: a Review”, Energy & Environmental Science,
4, 2011, 3243-3262.
[5] Z. Ogumi and M. Inaba, Advances in Lithium-Ion Batteries, Chapter 2: Carbon Anodes, Kluwer Academic/Plenum Publishers, United States, 2002.
[6] M.S. Whittingham, “Lithium Batteries and Cathode Materials”, Chemical Review,
104, 2004, 4271-4301.
[7] Y. Nishi, “Lithium-Ion Secondary Batteries with Gelled Polymer Electrolytes”,
Advances in Lithium-Ion Batteries, Chapter 7, Kluwer Academic /Plenum Publisher,2002.
[8] R. Baldwin, “A Review of State-of-the-Art Separator Materials for Advanced LithiumBased Batteries for Future Aerospace Missions”, NASA/TM-2009-215590, 2009.
[9] P. Arora and Z.M. Zhang, “Battery Separators”, Chemical Reviews, 104(10), 2004,
4419-4462.
[10] M.D. Farrington, “Proposed Amendments to UNST/SG/AC.10/11: Transport of Dangerous Goods-Lithium Batteries”, Journal of Power Sources, 80(1-2), 1999, 278-285.
[11] M.D. Farrington, “Safety of Lithium Batteries in Transportation”, Journal of Power
Sources, 96(1), 2001, 260-265.

177

[12] Y. T. Cheng and M.W. Verbrugge, “The Inﬂuence of Surface Mechanics on Diﬀusion
Induced Stresses within Spherical Nanoparticles”, Journal of Applied Physics, 104(8),
2008, 083521(1-6).
[13] J. Christensen and J. Newman, “Stress Generation and Fracture in Lithium Insertion
Materials”, Journal of Solid State Electrochemistry, 10(5), 2006, 293-319.
[14] J. Christensen and J. Newman, “A Mathematical Model of Stress Generation and
Fracture in Lithium Manganese Oxide”, Journal of he Electrochemical Society, 153(6),
2006, A1019-A1030.
[15] X. C. Zhang, W. Shyy, and A.M. Sastry, “Numerical Simulation of IntercalationInduced Stress in Li-Ion Battery Electrode Particles”, Journal of The Electrochemical
Society, 154(10), 2007, A910-A916.
[16] S. Prussin, “Generation and Distribution of Dislocations by Solute Diﬀusion”, Journal
of Applied Physics, 32(10), 1961, 1876-1881.
[17] K. E. Aifantis and J. P. Dempsey, “Stable Crack Growth in Nanostructured LiBatteries”, Journal of Power Sources, 143(1-2) 2005, 203-211.
[18] Y. T. Cheng, and M. W. Verbrugge, “Evolution of Stress within a Spherical Insertion Electrode Particle under Potentiostatic and Galvanostatic Operation”, Journal of
Power Sources, 190(2), 2009, 453-460 and 196, 2011, 2430-2431.
[19] R. E. Garcia, Y. M. Chiang, W. C. Carter, P. Limthongkul and C. M. Bishop, “Microstructural Modeling and Design of Rechargeable Lithium-Ion Batteries”, Journal
of The Electrochemical Society, 152(1), 2005, A255-A263.
[20] D. Bernardi, E. Pawlikowski and J. Newman, “A General Energy-Balance for Battery
Systems”, Journal of The Electrochemical Society, 132(1), 1985, 5-12.
[21] T. Ohshima, M. Nakayama, K. Fukuda, T. Araki and K. Onda, “Thermal Behavior
of Small Lithium-Ion Secondary Battery during Rapid Charge and Discharge Cycles”,
Electrical Engineering in Japan, 157(3), 2006, p. 17-25.
[22] L. Rao and J. Newman, “Heat-Generation Rate and General Energy Balance for Insertion Battery Systems”, Journal of The Electrochemical Society, 144(8), 1997, 26972704.
[23] K. E. Thomas and J. Newman, “Thermal Modeling of Porous Insertion Electrodes”,
Journal of The Electrochemical Society, 150(2), 2003, A176-A192.
178

[24] X. C. Zhang, A. M. Sastry and W. Shyy, “Intercalation-Induced Stress and Heat
Generation within Single Lithium-Ion Battery Cathode Particles”, Journal of The
Electrochemical Society, 155(7), 2008, A542-A552.
[25] S. C. Chen, C. C. Wan, and Y. Y. Wang, “Thermal Analysis of Lithium-Ion Batteries”,
Journal of Power Sources, 140(1), 2005, 111-124.
[26] Y. F. Chen and J. W. Evans, “Heat-Transfer Phenomena in Lithium PolymerElectrolyte Batteries for Electric Vehicle Application”, Journal of The Electrochemical
Society, 140(7), 1993, 1833-1838.
[27] Y. F. Chen and J. W. Evans, “Thermal-Analysis of Lithium Polymer Electrolyte Batteries by a 2-Dimensional Model-Thermal Behavior and Design Optimization”, Electrochimica Acta, 39(4), 1994, 517-526.
[28] Y. F. Chen and J. W. Evans, “Thermal Analysis of Lithium-Ion Batteries”, Journal
of The Electrochemical Society, 143(9), 1996, 2708-2712.
[29] C. R. Pals and J. Newman, “Thermal Modeling of the Lithium/Polymer Battery .1.
Discharge Behavior of a Single-Cell”, Journal of The Electrochemical Society, 142(10),
1995, 3274-3281.
[30] C. R. Pals and J. Newman, “Thermal Modeling of the Lithium/Polymer Battery .2.
Temperature Proﬁles in a Cell Stack”, Journal of The Electrochemical Society, 142(10),
1995, 3282-3288.
[31] M. Doyle, T. F. Fuller, and J. Newman, “Modeling of Galvanostatic Charge and Discharge of the Lithium Polymer Insertion Cell”, Journal of The Electrochemical Society,
140(6), 1993, 1526-1533.
[32] M. Doyle, J. Newman, A. S. Gozdz, C. N. Schmutz and J. M. Tarascon, “Comparison
of Modeling Predictions with Experimental Data from Plastic Lithium Ion Cells”,
Journal of The Electrochemical Society, 143(6), 1996, 1890-1903.
[33] K. E. Thomas, J. Newman, and R. M. Darling, Advances in Lithium-Ion Batteries,
Chapter 12: Mathematical Modeling of Lithium Batteries, Kluwer Academic/Plenum
Publishers, United States, 2002.
[34] W. Gu, “Thermal-Electrochemical Coupled Modeling of Dual-Insertion Batteries”,
Ph.D Thesis, The Pennsylvania State University, 2000.

179

[35] S. Golmon, K. Maute, and M. L. Dunn, “Numerical Modeling of ElectrochemicalMechanical Interactions in Lithium Polymer Batteries”, Computers and Structures,
87(23-24), 2009, 1567-1579.
[36] C. W. Wang and A.M. Sastry, “Mesoscale Modeling of a Li-Ion Polymer Cell”, Journal
of The Electrochemical Society, 154(11), 2007, A1035-A1047.
[37] M. Smith, R. E. Garcia and Q. C. Horn, “Characteration of the 3-Dimensional Microstructure of a Graphite Negative Electrode from a Li-Ion Battery”, Journal of The
Electrochemical Society, 156(11), 2009, A896-A904.
[38] W. Lai and F. Ciucci, “Mathematical Modeling of Porous Battery Electrodes- Revisit
of Newmans Model”, Electrochimica Acta, 56, 2011, 4369-4377.
[39] V Srinivasan and C. Y. Wang, “Analysis of Electrochemical and Thermal Behavior of
Li-Ion Cells”, Journal of The Electrochemical Society, 150(1), 2003, A98-A106.
[40] L. Cai and R. E. White. “An Eﬃcient Electrochemical-Thermal Model for a LithiumIon Cell by Using the Proper Orthogonal Decomposition Method”, Journal of The
electrochemical Society, 157(11), 2010, A1188-A1195.
[41] Y. Ye, Y. Shi, N. Cai, J. Lee and X. He, “Electro-thermal modeling and experimental
validation for lithium ion battery”, Journal of Power Sources, 199, 2012, 227-238.
[42] J. Alzieu, N. Koechlin and J. Robert, “Internal Stress Variations in Lead-Acid Batteries
during Cycling”, Journal of The Electrochemical Society, 134(8), 1987, 1881-1884.
[43] P. Bauerlein, C. Antonius, J. Loﬄer and J. Kumpers, “Progress in High-Power NickelMetal Hydride Batteries”, Journal of Power Sources, 176, 2008, 547-554.
[44] A. Etiemble, H. Idrissi, S. Meille and L. Roue, “In Situ Investigation of the Volume Change and Pulverization of Hydride Materials for Ni-MH Batteries by Concomitant Generated Force and Acoustic Emission Measurements”, Journal of Power
Sources,205, 2012, 500-505.
[45] H. Yakabe and I. Yasuda, “Model Analysis of the Expansion Behavior of LaCrO3
Intercalation under Solid Oxide Fuel Cell Operation”, Journal of The Electrochemical
Society, 150(1), 2003, A35-A45.
[46] H. Yakabe, T. Ogiwara, M. Hishinuma and I. Yasuda, “3-D model calculation for
planar SOFC”, Journal of Power Sources, 102, 2011, 144-154.
180

[47] A. Nakajo, Z. Wuillemin, J. V. Herle and Daniel Favrat, “Simulation of Thermal
Stresses in Anode-Supported Solid Oxide Fuel Cell Stacks. Part I: Probility of failure
of the cells”, Journal of Power Sources, 193, 2009, 203-215.
[48] J. Newman, Electrochemical Systems, 3th ed., United States: Wiley-Interscience, 2004.
[49] Y. Benveniste, “A New Approach to the Application of Mori-Tanaka Theory in
Composite-Materials”, Mechanics of Materials, 6(2), 1987, 147-157.
[50] T. Mori and K. Tanaka, “Average Stress in Matrix and Average Elastic Energy of
Materials with Misﬁtting Inclusions”, Acta Metallurgica, 21(5), 1973, 571-574.
[51] M. W. Verbrugge and Y. T. Cheng, “Stress and Strain-Energy Distributions within
Diﬀusion-Controlled Insertion-Electrode Particles Subjected to Periodic Potential Excitations”, Journal of The Electrochemical Society, 156(11), 2009, A927-A937.
[52] R. Deshpande, Y. T. Cheng and M. W. Verbrugge, “Modeling Diﬀusion-Induced Stress
in Nanowire Electrode Structures”, Journal of Power Sources, 195, 2010, 5081-5088.
[53] S. Golmon, K. Maute, S. H. Lee and M. L. Dunn, “Stress Generation in Silicon Particles
during Lithium Insertion”, Applied Physics Letters, 97, 2010, 03311(1-3).
[54] G. G. Stoney, “The Tension of Metallic Films Deposited by Electrolysis”, Proceedings
of the Royal Society of London”, 82, 1909, 172-175.
[55] J. Scarminio, S. N. Sahu, F. Decker, “In Situ Measurements of the Stress Changes in
Thin-Film Electrodes”,Journal of Physics E: Scientiﬁc Instruments, 22, 1989, 755-757.
[56] K. Y. Chung and K. B. Kim, “Investigation of Structural Fatie in Spinel Electrodes
Using in situ Laser Probe Beam Deﬂection Technique”, Journal of The Electrochemical
Society, 149 (1), 2002, A79-A85.
[57] H. Mukaibo, T. Momma, Y. Shacham-Diamand, T. Osaka and M. Kodaira, “In
situ Stress Transition Observations of Electrodeposited Sn-Based Anode Materials
for Lithium-Ion Secondary Batteries”, Electrochemical and Solid-State Letters, 10(3),
2007, A70-A73.
[58] Su-I1 Pyun, J. Y. Go, T. S. Jang, “An Investigation of Intercalation-Induced Stresses
Generated during Lithium Transport through LiCoO2 ﬁlm electrode using a laser beam
deﬂection method”, Electrochimica Acta, 49, 2004, 4477-4486.

181

[59] V. A. Sethuraman, V. Srinivasan, A. F. Bower and P. R. Guduru, “In Situ Measurement of Stress-Potential Coupling in Lithiated Silicon”, Journal of The Electrochemical
Society, 157(11), 2010, A1253-A1261.
[60] V. A. Sethuraman, M. J. Chon, M. Shimshak, V. Srinivasan and P. R. Guduru, “In
Situ Measurements of Stress Evolution in Silicon Thin Films during Electrochemical
Lithiation and Delithiation”, Journal of Power Sources, 195, 2010, 5062-5066.
[61] T. M. Bandhauer, S. Garimella and T. F. Fuller, “A Critical Review of Thermal Issues
in Lithium-Ion Batteries”, Journal of The Electrochemical Society, 158(3), 2011, R1R5.
[62] L. Song and J. W. Evans, “Electrochemical-Thermal Model of Lithium Polymer Batteries”, Journal of The Electrochemical Society, 147(6), 2000, 2086-2095.
[63] Y. F. Chen, and J. W. Evans, “Three-Dimensional Thermal Modeling of Lithium Polymer Batteries under a Galvanostatic Discharge and Dynamic Power Proﬁle”, Journal
of The Electrochemical Society, 141(11), 1994, 2947-2955.
[64] M. Guo, G. Sikha and R. E. White, “Single-Particle Model for a Lithium-Ion Cell:
Thermal Behavior”, Journal of The Electrochemical Society, 158(2), 2011, A122-A132.
[65] S. A. Hallaj, J. Prakash and J. R. Selman, “Cheracterization of Commercial Li-Ion Batteries Using Electrochemical-Calorimetric Measurements”, Journal of Power Sources,
87, 2000, 186-194.
[66] J. S. Hong, H. Malki, S. Al Hallaj, L. Redey and J. R. Selman, “ElectrochemicalCalorimetric Studies of Lithium-Ion Cell”, Journal of The Electrochemical Society,
145(5), 1998, 1489-1501.
[67] S. A. Hallaj, R. Venkatachalapathy, J. Prakash and J. R. Selman, “Entropy Changes
Due to Structural Transformation in the Graphite Anode and Phase Change of the
LiCoO2 Cathode”, Journal of The Electrochemical Society, 147(7), 2000, 2432-2436.
[68] H. Bang, H. Yang, Y. K. Sun and J. Prakash, “In Situ Studies of Lix Mn2 O4 and
Lix Li0.7 Mn1.83 O3.97 S0.03 Cathode by IMC”, Journal of The Electrochemical Society,
152(2), 2005, A421-A428.
[69] K. Onda, H. Kameyama, T. Hanamoto and K. Ito, “Experimental Study on Heat
Generation Behavior of Small Lithium-Ion Secondary Batteries”, Journal of The Electrochemical Society, 150(3), 2003, A285-A291.
182

[70] J. S. Kim, J. Prakash and J. R. Selman, “Thermal Characteristics of Lix Mn2 O4
Spinel”, Electrochemical and Solid-State Letters, 4(9), 2001, A141-A144.
[71] W. Lu and J. Prakash, “In Situ Measurements of Heat Generation in a Li/Mesocarbon
Microbead Half-Cell”, Journal of The Electrochemical Society, 150(3), 2003, A262A266.
[72] W. Lu, I. Belharouak, S. H. Park, Y. K. Sun and K. Amine, “Isothermal Calorimetry
Investigation of Li1+x Mn2−y Alz O4 Spinel”, Electrochimica Acta, 52, 2007, 5837-5842.
[73] H. Yang and J. Prakash, “Determination of the Reversible and Irreversible Heats of
a LiNi0.8 Co0.15 Al0.05 O2 /Natural Graphite Cell Using Electrochemical-Calorimetric
Technique”, Journal of The Electrochemical Society, 151(8), 2004, A1222-A1229.
[74] Y. Kobayashi, N. Kihira, K. Takei, H. Miyashiro, K. Kumai, N. Terada and R. Ishikawa,
“Electrochemical and Calorimetric Approach to Spinel Lithium Manganese Oxide”,
Journal of Power Sources, 81-82, 1999, 463-466.
[75] Y. Kobayashi, H. Miyashiro, K. Kumai, K. Takei, T. Iwahori and I. Uchida, “Precise
Electrochemical Calorimetry of LiCoO2 /Graphite Lithium-Ion Cell”, Journal of The
Electrochemical Society, 149(8), 2002, A978-A982.
[76] W. Lu, I. Belharouak, D. Vissers and K. Amine, “In Situ Thermal Study of
Li1+x [Ni1/3 Co1/3 Mn1/3 ]1−x O2 Using Isothermal Micro-clorimetric Techniques”,
Journal of The Electrochemical Society, 153(11), 2006, A2147-A2151.
[77] W. Lu, H. Yang and J. Prakash, “Determination of the Reversible and Irreversible
Heats of LiNi0.8 C00.2 O2 /Mesocarbon Microbead Li-Ion Cell Reactions Using Isothermal Microcalorimetery”,Electrochimica Acta, 51, 2006, 1322-1329.
[78] Y. Saito, K. Kanari and K. Takano, “Thermal Studies of a Lithium-Ion Battery”,
Journal of Power Sources, 68, 1997, 451-454.
[79] Y. Saito, K. Takano, K. Kanari, A. Negishi, K. Nozaki and K. Kato, “Comparative
Study of Thermal Bahaviors of Various Lithium-Ion Cells”, Journal of Power Sources,
97-98, 2001, 688-692.
[80] H. Vaidyanathan, W. H. Kelly and G. Rao, “Heat Dissipation in a Lithium Ion Cell”,
Journal of Power Sources, 93, 2001, 112-122.

183

[81] M. Majima, T. Tada, S. Ujiie, E. Yagasaki and S. Inazawa, “Design and Characteristics
of Large-Scale Lithium Ion Battery”, Journal of Power Sources, 81-82, 1999, 877-881.
[82] J. Shin, C. Han, U. Jung, S. Lee, H. Kim and K. Kim, “Eﬀect of LiCO3 Additive on
Gas Generation in Lithium-Ion Batteries”, Journal of Power Sources, 109 (1), 2002,
47-52.
[83] COMSOL 3.5a, “Rechargeable Lithium-Ion Battery”, Chemical Engineering Module:
Model Library, 2009.
[84] I. V. Thorat, D. E. Stephenson, N. A. Zacharias, K. Zaghib, J. N. Hard and D. R.
Wheeler, “Quantifying Tortuosity in Porous Li-ion Battery Materials”, Journal of
Power Sources, 188, 2009, 592-600.
[85] K. K. Patel, J. M. Paulsen and J. Desilvestro, “Numerical Simulation of Porous Networks in Relation to Battery Electrodes and Separators”, Journal of Power Sources,
122(2), 2003, 144-152.
[86] D. Djian, F. Allion, S. Martinet, H. Lignier, J.Y. Sanchez, “Lithium-Ion Batteries with
High Charge Rate Capacity: Inﬂuence of the Porous Separator”, Journal of Power
Sources, 172, 2007, 416-421.
[87] D. M. Bernadi and J. Y. Go, “Analysis of Pulse and Relaxation Behavior in LithiumIon Batteries”, Journal of Power Sources, 196(1), 2011, 412-427.
[88] K. Kumaresan, G. Sikha and R. E. White, “Thermal Model for a Li-Ion Cell”, Journal
of The Electrochemical Society, 155(2), 2008, A162-A171.
[89] D. E. Stephenson, E. M. Hartman, J. N. Harb, and D. R. Wheeler, “Modeling of
Particle-Particle Interactions in Porous Cathodes for Lithium-Ion Batteries”, Journal
of The Electrochemical Society, 154(12), 2007, A1146-A1155.
[90] D. E. Stephenson, B. C. Walker, C. B. Skelton, E. P. Gorzkowski, D. J. Rowenhorst
and D. R. Wheeler, “Modeling 3D Microstructure and Ion Transport in Porous Li-Ion
Battery Electrodes”, Journal of The Electrochemical Society, 158(7), 2011, A781-A789.
[91] F. Q. Yang, “Interaction between Diﬀusion and Chemical Stresses”, Materials Science
and Engineering A, 409, 2005, 153-159.
[92] Y. Qi, H. Guo, L. G. J. Hector and A. Timmons, “Threefold Increase in the Young’s
Modulus of Graphite Negative Electrode during Lithium Intercalation”, Journal of
The Electrochemical Society, 157(5), 2010, A558-566.
184

[93] R. Deshpane, Y. Qi and Y. T. Cheng, “Eﬀect of Concentration-Dependent Elastic
Modulus on Diﬀusion induced Stresses for Battery Applications”, Journal of The Electrochemical Society, 157, 2010, A967-A971.
[94] V. B. Shenoy, P. Johari and Y. Qi, “Elastic Softening of Amorphous and Crystalline
Li-Si Phases with Increasing Li Concentration: a First-Principle Study ”, Journal of
Power Sources, 195(19), 2010, 6825-6830.
[95] S. Timoshenko, Theory of Elasticity, 3rd edition, McGraw-Hill, United States, 1980.
[96] M. T. Shaw and W. J. Macknight, Introduction to Polymer Viscoelasticity, 3th ed.,
United States: Wiley-Interscience, 2005.
[97] T. Kanit, S. Forest, I. Galliet, V. Mounoury and D. Jeulin, “Determination of the
Size of the Representative Volume Element for Random Composites: Stratistical ans
Numerical Approach”, International Journal of Solids and Structures, 40(13-14), 2003,
3647-3679.
[98] A. D. Drozdov, S. Agarwal and R. K. Gupta, “Linear Thermo-Viscoelasticity of Isotactic Polypropylene”, Computational Materials, Science, 29(2), 2004, 195-213.
[99] T. L. Kulova, A. M. Skundin, E. A. Nizhnikovskii and A. V. Fesenko, “Temperature Effect on the Lithium Diﬀusion Rate in Graphite ”, Russian Journal of Electrochemistry,
42(3), 2006, 259-262.
[100] O. Yu. Egorkina and A. M. Skundin, “The Eﬀect of Temperature on Lithium Interaction into Carbon Materials”, Journal of the Solid State Electrochemistry, 2, 1998,
216-220.
[101] W. Tao, Q. Zhuang, C. Wu, Y. Cui, L. Fang and S. Sun, “Eﬀects of Temperature on
the Intercalation-Deintercalation Process of Lithium Ion in the Spinel LiMnO4 ”, Acta
Chimica Sinica, 68(15), 2010, 1481-1486.
[102] S. Kobayashi and Y. Uchimoto, “Lithium Ion Phase-Transfer Reaction at the Interface
between the Lithium Manganese Oxide Electrode and the Nonaqueous Electrolyte”,
The Journal of Physical Chemistry B, 109, 2005, 13322-13326.
[103] V. A. Sethuraman, L. J. Hardwick, V. Srinivasan and R. Kostecki, “Surface Structural Disordering in Graphite upon Lithium Intercalation/Deinterclation ”, Journal of
Power Sources,195, 2010, 3655-3660.

185

[104] Y. Reynier, R. Yazami and B. Fultz, “The Entropy and Enthalpy of Lithium Intercalation into Graphite”, Journal of Power Sources, 119-121, 2003, 850-855.
[105] K. E. Thomas, C. Bogatu and J. Newman, “Measurement of the Entropy of Reaction
as a Function of State of Charge in Doped and Undoped Lithium Manganese Oxide ”,
Journal of The Electrochemical Society, 148(6), 2001, A570-A575.
[106] L. O. Valoen and J. N. Reimers, “Transport Properties of LiPF6 -Based Li-Ion Battery
Electrolytes ”, Journal of The Electrochemical Society, 152(5), 2005, A882-A891.
[107] C. Capiglia, Y. Saito, H. Kagayama, P. Mustarelli, T. Iwamoto, T. Tabuchi and H.
Tukamoto, “Li and F Diﬀusion Coeﬃcients and Thermal Properties of Non-Aqueous
Electrolyte Solutions for Rechargeable Lithium Batteries ”, Journal of Power Sources,
81-82, 1999, 859-862.
[108] C. Capiglia, Y. Saito, H. Yamamoto, H. Kageyama and P. Mustarelli, “Transport
Properties and Microstructure of Gel Polymer Electrolytes”, Electrochimica Acta, 45,
2000, 1341-1345.
[109] P. Arora, M. Doyle, A. S. Gozdz, R. E. White and J. Newman, “Comparison between
Computer Simulations and Experimental Data for High-Rate Discharges of Plastic
Lithium-ion Batteries ”, Journal of Power Sources, 88(2), 2000, 219-231.
[110] A. Sheidaei, X.R. Xiao, X. S. Huang and J. Hitt, “Mechanical Behavior of a Battery
Separator in Electrolyte Solutions ”, Journal of Power Sources, 196(20), 2011, 87288734.
[111] I. M. Ward and J. Sweeney, An Introduction to the Mechanical Properties of Solid
Polymers, 2nd ed., United States: John Wiley & Sons, 2004.
[112] A. S. Wineman, K. R. Rajagopal, “Mechanical Response of Polymers: An Introduction”, Cambridge University, 2000.
[113] A. D. Drozdov, S. Agarwal and R. K. Gupta, “Linear Thermo-viscoelasticity of Isotactic Polypropylene ”, Computational Materials Science, 29, 2004, 195-213.
[114] G. A. Huber and D. S. Decker, Engineering Properties of Asphalt Mixtures and the
Relationship to Their Performance, ASTM special technical publication; 1265, 1995.
[115] Shutian Yan, Private Communication.

186

[116] W. Wu, X. Xiao and X. Huang, “The Eﬀect of Battery Design Parameters on Heat
Generation and Utilization in a Li-Ion Cell”, Electrochimica Acta, 83, 2012, 227-240.
[117] J. D. Ferry, Viscoelastic Properties of Polymers, Wiley, United States, 1980.
[118] I. M. Ward and J. Sweeney, An Introduction to the Mechanical Properties of Solid
Polymers, 2nd ed., John Wiley & Sons, West Sussex, England, 2004.
[119] W. Wu, X. Xiao, X. Huang, “Modeling Heat Generation In a Lithium Ion Battery”,
ASME 5th international conference on energy sustainability & 9th fuel cell science,
Washington DC, 2011, ESFuelCell2011-54785.
[120] Z. K. Liu, “Thermodynamic Modeling of Organic Carbonates for Lithium Batteries ”,
Journal of the Electrochemical Society, 150(3), 2003, A359-A365.
[121] D. Shi, X. Xiao, X. Huang and H. Kia, “Modeling Stresses in the Separators in a Pouch
Lithium-Ion Cell ”, Journal of Power Sources, 196, 2011, 8129-8139.
[122] J. Kim, T. V. Nguyen and R. E. White, “Thermal Mathematical Modeling of a Multi
Cell Common Pressure Vessel Nickel-Hydrogen Battery ”, Journal of The Electrochemical Society, 139(10), 1992, 2781-2786.
[123] D. K. L. Tsang, B. J. Marsden, S. L. Fok and G. Hall, “Graphite Thermal Expansion
Relationship for Diﬀerent Temperature Ranges ”, Carbon, 43, 2005, 2902-2906.
[124] J. F. Shackelford and W. Alexande, Materials Science and Engineering Handbook, 3rd
ed., CRC, United States, 2001.
[125] X. Xiao, W. Wu and X. Huang, “A Multi-Scale Approach for the Stress Analysis of
Polymeric Separators in a Lithium-Ion Battery ”, Journal of Power Sources, 195, 2010,
7649-7660.
[126] K.R. Kganyago, P.E. Ngoepe, “Structural and Electronic Properties of Lithium Intercalated Graphite LiC6 ”, Physical Review B, 68, 2003, 205111(1-16).
[127] http://www.cevp.co.uk/generalgraphite.htm.
[128] R. P. Ramasamy, J. W. Lee and B. N. Popov, “Simulation of Capacity Loss in Carbon
Electrode for Lithium-Ion Cells during Storage ”, Journal of Power Sources, 166, 2007,
266-272.

187

[129] X. L. Peng, J. Y. Huang, L. Qin, C. Y. Xiong and J. Fang, “A Method to determine
Young’s Modulus of Soft Gels for Cell Adhesion “, Acta Mechanica Sinica, 25(4), 2009,
565-570.
[130] E. Romano, J. L. Trenzado, E. Gonzalez, J. S. Matos, L. Segade and E. Jimenez,
“Thermophysical Properties of Four Binary Dimethyl Carbonate+ 1-Alcohol Systems
as 288.15-313.15K ”, Fluid Phase Equilibria, 211, 2003, 219-240.
[131] C.H. Lu and S. W. Lin, “Inﬂuence of the Particle Size on the Electrochemical Properties
of Lithium Manganese Oxide ”, Journal of Power Sources, 97-98, 2001, 458-460.
[132] G. G. Botte, B. A. Johnson, and R. E. White, “Inﬂuence of Some Design Variables on
the Thermal Behavior of a Lithium-Ion Cell ”, Journal of The Electrochemical Society,
146, 1999, 914-923.
[133] D. M. Bernardi and J. Y. Go, “Analysis of Pulse and Relaxation Behavior in LithiumIon Batteries”, Journal of Power Sources, 196, 2011, 412-427.
[134] J. R. Barber, Intermediate MEchanics of Materials, Solid Mechanics and Its Applications, 2nd edition, Springer, United States, 2011.
[135] D. Kehrwald, P. R. Shearing,N. P. Brandon, P. K. Sinha, S. J. Harris, “Local Tortuosity Inhomogeneities in a Lithium battery Composite Electrode ”, Journal of The
Electrochemical Society, 158, 2011, A1393-A1399.
[136] Y. H. Chen, C. W. Wang, G. Liu, X. Y. Song, V. S. Battaglia and A. M. Sastry,
“Selection of Conductive Additives in Li-Ion Battery Cathodes: A Nemerical Study”,
Journal of The Electrochemical Society, 154(10), 2007, A978-A986.
[137] C. M. Doyle, “Design and Simulation of Lithium Rechargeable Batteries”, Ph.D Thesis,
University of California, Berkeley, 1995.
[138] P. R. Shearing, L. E. Howard, P. S. Jorgensen, N. P. Brandon and S. J. Harris, “Characteration of the 3-Dimensional Microstructure of a Graphite Negative Electrode from
a Li-Ion Battery”, Electrochemistry Communications, 12, 2010, 374-377.
[139] M. C. Lonergan, D. F. Shriver, and M. A. Ratner, “Polymer Electrolytes-the Importance of Ion-Ion Interactions In-Diﬀusion Dominted Behvior”, Electrochimica Acta,
13-14, 1995, 2041-2048.

188

[140] G. M. Goldin, A. M. Colclasure, A. H. Wiedemann and R. J. Kee, “Three-Dimensional
Particle-Resolved Models of Li-Ion Batteries to Assist the Evaluation of Empirical
Parameters in One-Dimensional Models”, Electrochimica Acta, 64, 2012, 118-129.
[141] G. B. Less, J. H. Seo, S. Ha, A. M. Sastry, J. Zausch, A. Latz, S. Schmidt, C. ieser,
D. Kehrwald and S. Fell, “Micro-Scale Modeling of Li-Ion Batteries: Parameterization
and Validation”, Journal of The Electrochemical Society, 159(6), 2012, A697-A704.
[142] M. Park, X. Zhang, M. Chung, G. B. Less nd A. M. Sastry, “A Review of Conduction
Phenomena in Li-Ion Batteries”, Journal of Power Sources, 195, 2010, 7904-7929.
[143] G. Liu, H. Zheng, X. Song and V. S. Battaglia, “Particles and Polymer Binder Interaction: A Controlling Factor in Lithium-Ion Electrode Performance”, Journal of The
Electrochemical Society, 159(3), 2012, A214-A221.
[144] R. B. Wright, C. G. Motloch, J. R. Belt, J. P. Christophersen, C. D. Ho, R. A.
Richardson, I. Bloom, S. A. Jones, V. S. Battaglia, G. L. Henriksen, T. Ukelhaeuser, D.
Ingersoll, H. L. Case, S. A. Rogers and R. A. Sutula, “Calendar- and Cycle- Life Studies
of Advanced Technology Development Program Generation 1 Lithium-Ion Batteries”,
Journal of Power Sources, 110, 2002, 445-470.
[145] H. Zheng, R. Yang, G. Liu, X. Song and V. S. Battaglia, “Cooperation between Active
Material, Polymeric Binder and Conductive Carbon Additive in Lithium Ion Battery
Cathode”, The Journal of Physical Chemistry C, 116, 2012, 4875-4882.
[146] I. Kovalenko, B. Zdyrko, A. Magasinski, B. Hertzberg, Z. Milicev, R. Burtovyy, I. Luzinov and G. Yushin, “A Major Constituent of Brown Algae for Use in High-Capacity
Li-Ion Batteries”, Science, 334, 2011, 75-79.
[147] L. A. Riley, A. S. Cavanagh, S. M. George, S-H Lee and A. C. Dillon, “Improved Mechanical Intergrity of ALD-Coated Composite Electrodes for Li-Ion Batteries”, Electrochemical and Solid-State Letters, 14(3), 2011, A29-A31.
[148] N. Yoshizawa, O. Tamaike, H. Hatori, K. Yoshizawa, A. Kondo and T. Abe, “TEM
and Electron Tomography Studies of Carbon Nanospheres for Lithium Secondary Batteries”, Carbon, 44(12), 2006, 2558-2564.
[149] Y. Ding, W. Di, Y. Jiang, F. Xu, Z. Long, F. Ren and P. Zhang, “The Morphological
Evolution, Mechanical Properties and Ionic Conductivities of Electrospinning P(VDFHFP) Membranes at Various Temperatures ”, Ionics, 15, 2009, 731-734.

189

[150] N. Angulakshmi, S. Thomas, K. S. Nahm, A. M. Stephan and R. N. Elizabeth, “Electrochemical and Mechanical Properties of Nanochitin-Incorporated PVDF-HFP-Based
Polymer Electrolytes for Lithium Batteries”, Ionics, 17, 2011, 407-414.
[151] V. Djokovic, D. Kostoski and M. D. Dramicanin, “Viscoelastic Behavior of Semicrystalline Polymer at Elevated Temperatures on the Basis of a Two-Process Model for
Stress Relaxation”, Journl of Polymer Science, 38, 2000, 3239-3246.
[152] Z. Chen, L. Christensen and J. R. Dahn, “A Study of the Mechanical and Electrical
Properties of a Polymer/Carbon Black Binder System Used in Battery Electrodes”,
Journal of Applied Polymer Science, 90, 2003, 1891-1899.
[153] Z. Chen, L. Christensen and J. R. Dahn, “Mechanical and Electrical Properties of
Poly(vinylidene ﬂuoride-tetraﬂuoroethylene-propylene)/Super-S Carbon Black Swelled
in Liquid Solvent as an Electrode Binder for Lithium-Ion Batteries”, Journal of Applied
Polymer Science, 91, 2004, 2958-2965.
[154] G. Liu, S. Xun, N. Vukmirovic, X. Song, P. Olalde-Velasco, H. Zheng, V. S. Battaglia,
L. Wang and W. Yang, “Polymers with Tailored Electronic Structure for High Capacity
Lithium Battery Electrodes”, Advanced Materials, 23, 2011, 46794683.
[155] V. A. Sethuraman, N. V. Winkle, D, P. Abraham, A. F. Bower and P. R. Guduru,
“Real-Time Stress Measurements in Lithium-Ion Battery Negative-Electrodes”, Journal of Power Sources, 206, 2012, 334-342.
[156] A. Mukhopadhyay, A. Tokranov, K. Sena, X. Xiao and B. W. Sheldon, “Thin Film
Graphite Electrodes with Low Stress Generation during Li-Intercalation”, Carbon, 49,
2011, 2742-2749.
[157] J. Christensen, “Modeling Diﬀusion-Induced Stress in Li-Ion Cells with Porous Electrodes”, Journal of The Electrochemical Society, 157, 2010, A366-A380.
[158] K. Zhao, “Mechanics of Electrodes in Lithium-Ion Batteries”, Ph. D Thesis, Harvard
University, 2012.
[159] R. Kostecki and F. McLarnon, “Microprobe Study of the Eﬀect of Li Intercalation on
the Structure of Graphite”, Journal of Power Sources, 119-121, 2003, 550-554.
[160] L. J. Hardwick, H. Buqa, M. Holzapfel, W. Esheifele, F. Krumeich and P. Novak, “Behaviour of Highly Crystalline Graphitic Materials in Lithium-Ion Cells with Propylene
Carbonate Containing Electrolytes: An In Situ Raman and SEM Study”, Electrochimica Acta, 52, 2007, 4884-4891.
190

[161] Y. Wang, D. C. Alsmeyer and R. L. McCreery, “Raman Spectroscopy of Carbon Materials: Structural Basis of Observed Spectra”, Chemistry of Materials, 2(5), 1990,
557-563.
[162] K. Suzuki, T. Hamada and T. Sugiura, “Eﬀect of Graphite Surface Structure on Initial Irreersible Reaction in Graphite Anodes”, Journal of The Electrochemical Society,
146(3), 1999, 890-897.
[163] K. A. Striebel, J. Shim, E. J. Carins, R. Kostecki, Y. J. Lee, J. Reimer, T. J. Richardson, P. N. Ross, X. Song and G. V. Zhuang, “Diagnostic Analysis of Electrodes from
High-Power Lithium-Ion Cells Cycled under Diﬀerent Conditions”, Journal of The
Electrochemical Society, 151(6), 2004, A857-A866.
[164] V. B. Shenoy, P. Johari and Y. Qi, “Elastic Softening of Amorphous and Crystalline
Li-Si Phases with Increasing Li Concentration: A First-Principles Study”, Journal of
Power Sources, 195, 2010, 6825-6830.
[165] P. Limthongkul, Y.-I. Jang, N. J. Dudney, and Y.-M. Chiang, “ElectrochemicallyDriven Solid-State Amorphization In Lithium-Silicon Alloys and Implications for
˚
Lithium Storage’, Acta Material, 51, 2003, 1103.
[166] Kejie Zhao, G. A. Trisaris, M. Pharr, W. Wang, O. Okeke, Z. Suo, J. J. Viassak, E.
Kaxiras, “Reactive Flow in Silicon Electrodes Assisted by the Insertion of Lithium”,
Nano Letters, 12, 2012, 43974403.
[167] V. A. Sethuraman, M. J. Chon, M. Shimshak, V. Srinivasan and P. R. Guduru, “In
Situ Measurements of Stress Evolution in Silicon Thin Films during Electrochemical
Lithiation and Delithiation”, Journal of Power Sources, 195(15), 2010, 5062-5266.
[168] D. Dini, S. Passerini, B. Scrosati and F. Decker, “Stress Changes in Electrochromic
Thin Film Electrodes: Laser Beam Deﬂection Method (LBDM) as a Tool for the
Analysis of Intercalation Processes”, Solar Energy Materials & Solar Cells, 56, 1999,
213-221.

191