Three-Dimensional Lithium-Ion Battery Model

Understanding Spatial Variations in Battery Physics to Improve Cell 

Design, Operational Strategy, and Management

Gi-Heon Kim and Kandler Smith

4th International Symposium on Large Lithium Ion Battery Technology and Application

Tampa, Florida

May 12–14, 2008

gi_heon_kim@nrel.gov, kander_smith@nrel.gov

NREL/PR-540-43166
Multi-Scale Physics in Li-ion Battery

Li+

e

Li+

LixC6

Li+

e

Li+

LixC6

Li+

e

Li+

Liy-MetalOxide

Li+

e

Li+

Liy-MetalOxide

Solid-Phase Lithium Transport

Lithium Transport in Electrolyte

Electrochemical Kinetics

η = (φs −φe )−U

Charge Conservation/Transport

(Thermal) Energy Conservation

• Kinetics

• Phase transition

• Ion transport

• Energy dissipation

• Heat transfer 

Basic battery physics occurs in a 

wide range of length & time scales
Requirements & Resolutions 

10 [m] -10 10-8 10-6 10-4 10-2 100

“Requirements” are usually defined 

in a macroscale domain and terms.

Performance

Life

Cost

Safety

Design of Materials

Voltage

Capacity

Lattice stability

Kinetic barrier

Transport property

Design of Electrode 

Architecture

Li transport path (local)

Electrode surface area

Deformation & fatigue

Structural stability

Interface physics

Design of Electrodes 

Pairing and Lithium 

Transport

Electrodes selection

Li transport

Porosity, tortuosity

Layer thicknesses 

Load conditions

Design of 

Electronic Current 

& Heat Transport

Electric & thermal 

connections

Dimensions, form factor

Component shapes
0

50

100

150

200

0

50

100

150

3.08

3.09

3.1

3.11

3.12

3.13

X(mm)

Voltage [V]

Y(mm)

10-10 10-8 10-6 10-4 10-2 100 10-10 10-8 10-6 10-4 10-2 100

Li+

NREL’s Li-ion Battery Model Activities

focusing on different length scale physics

a) Quantum mechanical and molecular dynamic modeling

b) Numerical modeling resolving architecture of electrode materials

c) Electrode-scale performance model 

d) Cell-dimension 3D performance model 

[m]

a) b) c) d)
Why use a 3D model?

Enhanced understanding provides an opportunity for improving cell ...

• Design

• Operation Strategy

• Management

• Safety

Spatially Nonuniform Battery Physics

Dimension 

(Size effect, A/V )

Chemical Kinetics

Heat & Mass Transport

Charge Conservation

Spatial Distribution of ...

• Electrical Potentials

• Temperatures

Design of current and 

heat flow paths

Flux (gradient) 

Barrier (distance)

Spatial Difference

(gradient x distance)
Electrode-Scale Model

Chemical Kinetics

(k T ) q

t

T

c p = ∇⋅ ∇ + ′′′ ∂

∂

ρ

Heat & Mass Transport

Charge Conservation

This model captures relevant solid-state and 

electrolyte diffusion dynamics and predicts the 

current/voltage response of a battery. 

Composite electrodes are modeled using porous 

electrode theory, meaning that the solid and 

electrolyte phases are treated as superimposed 

continua without regard to microstructure. 

Li+

r

NREL’s Model

• Finite-Volume Method

• Matlab Environment

(Doyle, Fuller, and Newman, 1993)
Addressing the effects of:

 Nonuniform distributions 

 Thermal/electrical path design inside cells/batteries

 Localized phenomena 

 Geometries; shape and dimensions of cell component

3D Battery Dimension Model

Current Density

Temperature
Approach in the Present Study:

Multi-Scale Multi-Dimensional (MSMD) Modeling Current Collector (Cu)

Current Collector (Al)

p

Negative

Electrode

Separator

Positive

Electrode

Image source: www.dimec.unisa.it

Simulation

Domain

X

R

x

=

• Multi-scale physics from sub-micro-scale to battery-dimension-scales

• Difficulties in resolving microlayer structures in a computational grid

Macro Grid Micro Grid

(Grid for Sub-grid-scale Model)

+

To Address ...
T(X , R,t) (X , R, x,t) s φ

(X , R, x,t) e φ

c (X , R, x,r,t) s

c (X , R, x,t) eCurrent Collector (Cu)

Current Collector (Al)

p

Negative

Electrode

Separator

Positive

Electrode

Image source: www.dimec.unisa.it

Detailed 

Structure

X

R

x

≈ Cell Dimension 

Transport Model

Electrode Scale 

Submodel (1D) +

V (X , R,t)

Q (X , R, x,t) i

j (X , R, x,t) Li

SOC(X , R,t)

i(X , R,t)

V

Adx Q X R t Q

x

= ∫ i ( , , )

NOTE:

Selection of the “sub-grid 

electrochemical model” is 

independent of the “macrogrid model” selection.

Solution Variables
T∞=35oC

h= 25W/m2K 163mm

3mm

27mm

kX= 29W/mK

kR= 0.8W/mK

X

R

time step size: 0.5 sec, time step #: 360

Computation time: 98 minutes (Windows/PC)

Simulation time: 3 minutes

Axisymmetric FVM Model for Macro-Domain Model 

+ 1D FVM Model for Electrochemistry Submodel

Model Combination 

150A Constant Current Discharge

X R x r

Grid# 12 8 14 4

0 0.5 1 1.5 2 2.5

2.5

3

3.5

4

Time(min)

Voltage (V)

V

o

V

ocv
0 0.5 1 1.5 2 2.5

2.5

3

3.5

4

Time(min)

Voltage (V)

Vo

V

Simulation Results ocv Snapshot 2 minutes after start of discharge

X(mm)

R(mm)

Normalized ANODE SURFACE CONCENTRATION

0 20 40 60 80 100 120 140 160

5

10

15

20

25

0.99

0.995

1

1.005

1.01

X(mm)

R(mm)

Normalized CHATHODE SURFACE CONCENTRATION

0 20 40 60 80 100 120 140 160

5

10

15

20

25

0.99

0.995

1

1.005

1.01

X(mm)

R(mm)

TEMPERATURE [o

C]

0 20 40 60 80 100 120 140 160

5

10

15

20

25

48

49

50

X(mm)

R(mm)

CURRENT PRODUCTION [A/m2]

0 20 40 60 80 100 120 140 160

5

10

15

20

25

162

164

166

168

X

R
For details about the State Variable Model:

See the Poster Presentation by Kandler Smith (NREL) titled,

“Fast Running Electrochemistry-Based Models for Battery Design, Integration, and Control”

SVM is preferred because of its fast execution

SVM, developed by Kandler Smith (NREL), quickly solves “Newman type” 

governing equations using numerical schemes for calculating load reduction.

Dropping very fast battery responses (approx. 60 Hz or more) is one of the main 

calculation order reduction methods used in the model.

SVM is promising for use in on-board BMS reference model because of its 

fast execution and capability to provide nonmeasurable electrochemical 

parameters and current and voltage responses with potentially better accuracy.

MSMD model incorporating SVM Submodel runs ~1.75 faster than real time.

Axisymmetric FVM Model for Macro-Domain Model 

+ State Variable Model (SVM) for Submodel

Another Combination Choice 
Analysis

Temperature Variation in a Cylindrical Cell

Uniform Potential Assumption

Impact of Aspect Ratio

Impact of Cell Size

Temperature & Potential Variation in a Prismatic Cell

Impact of Tab Location and Size
Analysis

Temperature Variation in a Cylindrical Cell

Uniform Potential Assumption

Impact of Aspect Ratio

Impact of Cell Size

Temperature & Potential Variation in a Prismatic Cell

Impact of Tab Location and Size
Considerations for Addressing Thermal 

Issues in PHEV-type Cells

High energy and high power requirements

Large format may be preferred to small cells

Fewer number of components 

Fewer interconnects

Less monitoring & balancing circuitry

Less expensive

Less weight

Significant heating may be possible, depending on power 

profile

Internal temperature imbalance can lead to unexpected 

performance and life degradation
For a fixed capacity (electrode volume), surface area for heat rejection 

can be increased via:

• Reducing D/H ratio

• Increasing number of cells in parallel (#P) 1P

2P

3P

40Ah

D

H

*Surface area includes side, top & bottom of can. All cells assumed to 

have inactive inner mandrel with 8mm diameter.

Analysis Parameters

20Ah

13.3 Ah
The two cases explored in this presentation:

2

Tamb= 35 ̊C

SOC = 50%

• Steady-State, Charge Sustaining

• Liquid Cooling (150 W/m2K)

2 200A Geometric Cycle

Tamb= 35 ̊C

SOC = 50%

Tamb= 35 ̊C

SOC = 50%

• Steady-State, Charge Sustaining

• Liquid Cooling (150 W/m2K)

200A Geometric Cycle

Moderate Thermal Condition Severe Thermal Condition

Two Usage Profiles

1

• Transient, Charge Depleting

• Natural Convection (15 W/m2K)

Tinit = 35 ̊C

Tamb= 35 ̊C

SOCinit = 100%

1 150A Max. Cont. Discharge

• Transient, Charge Depleting

• Natural Convection (15 W/m2K)

Tinit = 35 ̊C

Tamb= 35 ̊C

SOCinit = 100%

Tinit = 35 ̊C

Tamb= 35 ̊C

SOCinit = 100%

150A Max. Cont. Discharge

Air

1

• Transient, Charge Depleting

• Natural Convection (15 W/m2K)

Tinit = 35 ̊C

Tamb= 35 ̊C

SOCinit = 100%

1 150A Max. Cont. Discharge

• Transient, Charge Depleting

• Natural Convection (15 W/m2K)

Tinit = 35 ̊C

Tamb= 35 ̊C

SOCinit = 100%

Tinit = 35 ̊C

Tamb= 35 ̊C

SOCinit = 100%

150A Max. Cont. Discharge

Air
T

Q

i

D/H = 1⁄4

h = 15 W/m2K

Tamb = 35°C

Transient Results

Results:

150 A Continuous 

Discharge

After 500 seconds of discharge:

• Cell center is slightly warmer 

than exterior

• Preferential reaction current at 

cell center

Near the end of discharge:

• Cell center depleted/saturated

• Preferential reaction current at 

cell exterior
D/H = 1⁄4

h = 15 W/m2K

Tamb = 35 C

1P

2P

3P

Moderate usage + air convection = small internal gradients

2.7 ̊C

150 A Single Discharge (at End)

40Ah

D

H

20Ah

13.3 Ah
200 A Geometric Cycling

∆T= ~17 oC 

~16% difference in

local current production

At Steady State

( )A

I

i

T −Tavg

D/H = 1⁄4

h = 150 W/m2K

Tamb = 35°C

1P size cell (40 Ah)
1P

D/H = 1⁄4

h = 150 W/m2K

Tamb = 35 C

2P

3P

Severe usage + liquid cooling = large internal gradients

18 ̊C

200 A Geometric Cycle (Steady-State)

40Ah

D

H

20Ah

13.3 Ah
h = 150 W/m2K

Tamb = 35°C

• Under severe usage, low D/H 

and/or >1P designs significantly 

reduce thermal stress

• Larger diameter leads to higher 

internal gradient

• Multidimensional electrochemical 

cell model quantified the impacts 

of D/H aspect ratio and cell size 

on the internal temperature 

difference.

Internal Temperature Difference

D/H ratio ~4.5 ̊C ↕

2P ~6.0 ̊C, 3P ~9.0 ̊C ↓ Tavg

Tamb

ΔTint

-1

200 A Geometric Cycle (Steady-State)
Analysis

Temperature Variation in a Cylindrical Cell

Uniform Potential Assumption

Impact of Aspect Ratio

Impact of Cell Size

Temperature & Potential Variation in a Prismatic Cell

Impact of Tab Location and Size
Impact of Tab Location & Size

200 mm

150 mm

• Thickness: 12 mm

• 40 Ah

• 2-minute discharge, 200 A

• 200A geometric cycle

30 mm

60 mm
200A Discharge for 2 minutes

0 20 40 60 80 100 120

0

50

100

150

200

250

300

t (s)

I(t) (A)
Current Field – 2-min 200 A discharge

X(mm)

Y(mm)

Current Density in positive foil

0 20 40 60 80 100 120 140 160 180 200

0

50

100

150

X(mm)

Y(mm)

Current Density in negative foil

0 20 40 60 80 100 120 140 160 180 200

0

50

100

150

X(mm)

Y(mm)

Current Density in positive foil

0 20 40 60 80 100 120 140 160 180 200

0

50

100

150

X(mm)

Y(mm)

Current Density in negative foil

0 20 40 60 80 100 120 140 160 180 200

0

50

100

150
Voltage across Current Collector Foils 

– 2-min 200 A discharge

Vmax-Vmin= 0.0364 V Vmax-Vmin = 0.0154 V

0

50

100

150

200

0

50

100

150

3.08

3.09

3.1

3.11

3.12

3.13

X(mm)

Voltage [V]

Y(mm) 0

50

100

150

200

0

50

100

150

3.08

3.09

3.1

3.11

3.12

3.13

X(mm)

Voltage [V]

Y(mm)
Temperature – 2-min 200 A discharge

0

50

100

150

200

0

50

100

150

46

47

48

49

50

51

52

X(mm)

Temperature [°C]

Y(mm) 0

50

100

150

200

0

50

100

150

46

47

48

49

50

51

52

X(mm)

Temperature [°C]

Y(mm)

Tmax-Tmin = 5.03oC Tmax-Tmin = 1.35oC

0 20 40 60 80 100 120

35

40

45

50

t (s)

T(t) (°C)

avarage temperature

maximum temperature

minimum temperature

0 20 40 60 80 100 120

35

40

45

50

t (s)

T(t) (°C)

avarage temperature

maximum temperature

minimum temperature
imax-imin = 4.54 A/m2

0

50

100

150

200

0

50

100

150

155

160

165

170

175

X(mm)

Current Production [A/m2]

Y(mm) 0

50

100

150

200

0

50

100

150

155

160

165

170

175

X(mm)

Current Production [A/m2]

Y(mm)

imax-imin = 13.2 A/m2

Current Production – 2-min 200 A discharge
SOC – 2-min 200 A discharge

SOCmax-SOCmin = 1.91%

0

50

100

150

200

0

50

100

150

59.5

60

60.5

61

61.5

62

X(mm)

soc [%]

Y(mm) 0

50

100

150

200

0

50

100

150

59.5

60

60.5

61

61.5

62

X(mm)

soc [%]

Y(mm)

SOCmax-SOCmin = 0.76%
200A Geometric Cycling

0 200 400 600 800

-300

-200

-100

0

100

200

300

t (s)

I(t) (A)
Temperature Variation

0 200 400 600 800

35

40

45

50

55

t (s)

T(t) (°C)

avarage temperature

maximum temperature

minimum temperature

0

50

100

150

200

0

50

100

150

-3

-2

-1

0

1

2

3

X(mm)

T - avg(T)

Y(mm)
0

50

100

150

200

0

50

100

150

42.5

42.6

42.7

42.8

42.9

43

X(mm)

soc [%]

Y(mm)

SOC swing

0

50

100

150

200

0

50

100

150

56.5

57

57.5

58

X(mm)

 []

Y(mm)

SOC [%]

0

100

200

0

50

100

150

14

14.5

15

15.5

X(mm)

soc swing [%]

Y(mm)
Summary

Nonuniform battery physics, which is more probable in large-format 

cells, can cause unexpected performance and life degradations in 

lithium-ion batteries. 

A three-dimensional cell performance model was developed by 

integrating an electrode-scale submodel using a multiscale modeling

scheme. 

The developed tool will be used to provide better understanding and 

help answer engineering questions about improving cell design, cell 

operational strategy, cell management, and cell safety.

Engineering Questions to be addressed in future works include ...

What is the optimum form-factor and size of a cell?

Where are good locations for tabs or current collectors?

How different are measured parameters from their non-measurable internal values? 

Where is the effective place for cooling? What should the heat-rejection rate be?

How does the design of thermal and electrical paths impact under current-related 

safety events, such as internal/external short and overcharge? 
Acknowledgments

DOE and FreedomCAR-Fuel Partnership Support

• Tien Duong

• Dave Howell

NREL Energy Storage Task

• Ahmad Pesaran