University of South Carolina University of South Carolina
Scholar Commons Scholar Commons
Faculty Publications Chemical Engineering, Department of
1987
A Mathematical Model of a Lead-Acid Cell: Discharge, Rest, and A Mathematical Model of a Lead-Acid Cell: Discharge, Rest, and
Charge Charge
Hiram Gu
T. V. Nguyen
Ralph E. White
University of South Carolina - Columbia
, white@cec.sc.edu
Follow this and additional works at: https://scholarcommons.sc.edu/eche_facpub
Part of the Chemical Engineering Commons
Publication Info Publication Info
Journal of the Electrochemical Society
, 1987, pages 2953-2960.
© The Electrochemical Society, Inc. 1987. All rights reserved. Except as provided under U.S. copyright law,
this work may not be reproduced, resold, distributed, or modiBed without the express permission of The
Electrochemical Society (ECS). The archival version of this work was published in the Journal of the
Electrochemical Society.
http://www.electrochem.org/
DOI: 10.1149/1.2100322
http://dx.doi.org/10.1149/1.2100322
This Article is brought to you by the Chemical Engineering, Department of at Scholar Commons. It has been
accepted for inclusion in Faculty Publications by an authorized administrator of Scholar Commons. For more
information, please contact digres@mailbox.sc.edu.
Vol. 134, No. 12
EFFICIENT VANADIUM REDOX FLOW CELL 2953
The independent sizing of the redox system's power
and storage capacity makes it particularly attractive for
long term electricity storage in remote areas where solar
arrays or wind generators have been installed. An at-
tractive feature of the vanadium redox battery is that
since there is no solution contamination due to cross-
mixing, after the original capital investment, there
would be negligible running and maintenance costs.
Since the cell can be completely discharged without any
deterioration in performance, it would be ideally suited
for large scale energy storage in load leveling applica-
tions.
A larger scale five-cell battery unit is currently under
construction so that the system performance can be fur-
ther tested before scaling up to a 1 kW unit.
Acknowledgments
Support for this project was provided under the Na-
tional Energy Research Development and Demonstra-
tion Programme which is administered by the Australian
Commonwealth Department of Resources and Energy.
The authors are grateful to Dr. Miron Rychcik for useful
discussions.
Manuscript submitted Sept. 2, 1986; revised manu-
script received Dec. 16, 1986.
The University of New South Wales assisted in meeting
the publication costs of this article.
REFERENCES
1. L. H. Thaller, NASA TMX-71540, National Aeronautics
and Space Administration, U.S. Department of En-
ergy, 1974; and U.S. Pat. 3,996,064 (1974).
2. L. H. Thaller, NASA TMM-79143, National Aeronautics
and Space Administration, U.S. Department of En-
ergy, 1979.
3. Redox Flow Cell Development and Demonstration
Project, NASA TM-79067, National Aeronautics and
Space Administration, U.S. Dept. of Energy, 1979.
4. H. H. Hagedorn and L. H. Thaller, NASA TM-81464,
National Aeronautics and Space Administration,
U.S. Department of Energy, 1980.
5. D. A. Johnson and M. A. Reid, NASA TM-82913, Na-
tional Aeronautics and Space Administration, U.S.
Department of Energy, 1982.
6. K. Nozaki, H. Kaneko, A. Negishi, and T. Ozawa, in
"Proceedings of Meeting" pp. 1641-1646, 18th IECEC
(1983).
7. R. F. Gahn, N. H. Hagedorn, and J. A. Johnson, NASA
TM87034, National Aeronautics and Space Adminis-
tration, U.S. Department of Energy, 1985.
8. M. Skyllas-Kazacos, M. Rychcik, and R. Robins,
Patent Applications filed.
9. M. Skyllas-Kazacos, M. Rychcik, R. Robins, A. Fane,
and M. Green, This Journal, 133, 1057 (1986).
10. R. Rychcik and M. Skyllas-Kazacos, J. Power Sources,
19, 45 (1986).
11. D.-G. Oei, J. Appl. Electrochem., 15, 231 (1985).
A Mathematical Model of a Lead-Acid Cell
Discharge, Rest, and Charge
Hiram Gu* and T. V. Nguyen, **'1
Electrochemistry Department, General Motors Research Laboratories, Warren, Michigan 48090-9055
R. E. White*
Department of Chemical Engineering, Texas A&M University, College Station, Texas 77843-3122
ABSTRACT
A mathematical model of a lead-acid cell is presented which includes the modeling of porous electrodes and various
physical phenomena in detail. The model is used to study the dynamic behavior of the acid concentration, the porosity of
the electrodes, and the state of charge of the cell during discharge, rest, and charge. The dependence of the performance of
the cell on electrode thicknesses and operating temperature is also investigated.
The lead-acid system is used in the largest number of
secondary batteries manufactured in the world. The
most important market remains the car battery for start-
ing, lighting, and ignition, with approximately 50 • 106
units sold per year in the U.S.A. (1). Other applications
are in emergency power supplies, load-leveling, and
more recently for instruments, radio, and other electrical
apparatus. The design and improvement of these
batteries are mostly done by trial-and-error.
This traditional approach, which consists of experi-
mental cell build-ups and extensive testing, is costly and
time consuming. Furthermore, results from such tests
provide only global information and do not provide in-
sight into the governing phenomena. It is advantageous
to develop a mathematical model of the cell which would
allow one to gain a better understanding of the cause and
effect relationships and the phenomena involved, and
suggest directions for improvements.
Complementing experimental testing with mathemat-
ical modeling is a cost effective approach to the develop-
ment and design of batteries. Testing is still needed to
verify predictions of the model and to uncover physical
phenomena that may not have been included in the
model. But with the help of this mathematical tool, ex-
*Electrochemical Society Active Member.
**:Electrochemical Society Student Member.
~Present address: Department of Chemical Engineering, Texas
A&M University, College Station, Texas 77843-3122.
tensive experimental testing is no longer needed. Great
savings in material, labor, and time can be realized in the
development of a new battery.
The first use of mathematical models to describe the
behavior of the lead-acid system was applied to the po-
rous positive electrode (PbO2) by Stein (2, 3) and Euler
(4), with further improvements made by Simonsson (5-7),
Micka and Rousar (8), Gidaspow and Baker (9), and oth-
ers. A good review of the development in the theory of
flooded porous electrodes prior to 1975 has been pro-
vided by Newman and Tiedemann (10). Recently, Tiede-
mann and Newman (11) and Sunu (12) applied Newman's
theory to the development of a.complete cell model de-
scribing the discharge behavior of the lead-acid battery
system. However, a model for predicting the cell behav-
ior during charge and rest, as well as the effects of cy-
cling, is not available. To assist designers and engineers
in the further development of the lead-acid batteries
with improved performance and cycle life, a detailed
mathematical model of a lead-acid cell is presented that
can be used to predict the dynamic behavior of the cell
not only during discharge, but also during charge, rest,
and cycling.
Model Development
A schematic for the lead-acid cell is shown in Fig. 1.
The cell consists of the following boundaries and re-
Downloaded 13 Jun 2011 to 129.252.106.227. Redistribution subject to ECS license or copyright; see http://www.ecsdl.org/terms_use.jsp
2954
J. Electrochem. Soc.:
ELECTROCHEMICAL SCIENCE AND TECHNOLOGY
December 1987
Region 1, positive electrode.--
porosity variation--
Oe
1 [ MWebso4
Ot
2F PPbSO4
Ohm's law in solution--
MWebo2 I Oi~
PPb02 ] 0X
0 [6]
is Od)~ RT
0 in (cf)
--+ -- (1- 2t~ -- -0 [7]
eexlK 0x F 0x
Ohm's law in solid--
i2
exml
a~l
_
I
=
0 [8]
-- {[ O'Pb02 0X
material balance--
OC
Fig. 1. One-dimensional macro-homogeneous model of a lead-acid cell 6
--
gions: a lead grid current collector at x = 0, which is at
the center of the positive electrode, the positive (PbO2)
electrode (region 1), the positive electrode/reservoir in-
terface, the reservoir (region 2), the reservoir/separator
interface, the separator (region 3), the separator/negative
electrode interface, the negative (Pb) electrode (region
4), and the center of the negative electrode where an-
other grid is located. Details of the geometry are ignored
and the whole cell is regarded as a homogeneous macro-
scopic entity with distributed quantities in the direction
perpendicular to the grid. An extensive discussion of av-
erage quantities used in the devel.opment of the model
has been given by Dunning (13) and Trainham (14). Addi-
tionally, isothermal conditions are assumed here. The
electrolyte is concentrated H2SO4 which is considered to
be a binary electrolyte that dissociates into H + and
HSO4- in H~O.
The electrode reactions during discharge
are
PbO2(s) + HSO4 + 3H + + 2e-
and
Pbls) + HSO(
discharge
) PbSO4(~) + 2H20
(positive electrode)
discharge
) PbSO4(~) + 2e- + H +
(negative electrode)
There are five explicit unknowns in the model: acid
concentration (c), electrode porosity (6), superficial cur-
rent density in the electrolyte (i2), potential in the solid
phase (qbl), and potential in the electrolyte ((b2). The gov-
erning equations and boundary conditions are presented
next.
Center of the positive electrode, x = O.--
Oc
- 0 [1]
ox
o6
- 0 [2]
Ox
i2 = 0 [3]
(bl =
0
[4]
O62
= 0 [5]
OX
Equations [1], [2], and [5] are in accordance with the as-
sumption of symmetry. Equation [3] states that, at the
center of the positive electrode, all the current is in the
current collector and none is in the electrolyte. In Eq. [4],
as a convenient choice, the solid phase potential, d)l, is
designated to be 0 V at this boundary. Without this refer-
ence potential, a particular solution cannot be obtained.
1F( wPbso4 WPbo2)
PPbS04
0t 2F PPb02
(3
2t~ + 2V ] i2 ac O ( Oc )
-- -- o --
E exl D
Ox Ox Ox
- . ~ 9 (1-cYo)+~ ox
electrode kinetics--for discharge and rest
~ -- 6ol
0i2
ama• ~ C~-e~ erna• Co1
0X
[ RT (~)1 - ~D2 --
P-l} o
exp [~- (q), - r - =
and for charge
0i2 * (C tYI( ~--~o1 )~1 e .... --E
--OX -- amaxll~ \ C-~ef ] 6maxl -- 6ol ~maxl -- ~ol
where
[9]
[10a]
exp [ c~IF Uebo~)] -
t RT (~1- ~b2-
[ -CtelF "'
Upb02)] } 0
exp [~ (r - r - =
[lOb]
Upb02 ~ Upbo20- Upb 8 [11]
Equation [6] describes the change in porosity with time
due to the conversion of the active material in the elec-
trode reaction. Equation [7] is a modified Ohm's law for
the electrolyte which states that the current in the elec-
trolyte is driven by the electric potential and chemical
potential gradients. Equation [8] is Ohm's law applied to
the solid matrix. Equation [9] states that the electrolyte
concentration at any point in space changes with time
because of the electrode reaction, diffusion, and migra-
tion. Equations [10a] and [10b] are kinetic expressions
for the electrode reaction. Equation [10b] includes a fac-
tor to account for the depletion of lead sulfate as lead di-
oxide is being formed. For convenience, the dependence
of the overpotential on the acid concentration is neg-
lected, and a concentration independent lead electrode
is used as a reference electrode in Eq. [10a], [10b], and
[11].
Interface between region 1 and 2.--
~C
region
OC
region
e .... [12]
~X ~ aX 2
0e 1 ( MWebso4
MWpboz)OiZ=o[13]
Ot
2F
PPbS04 '. PPbO2 0X
Downloaded 13 Jun 2011 to 129.252.106.227. Redistribution subject to ECS license or copyright; see http://www.ecsdl.org/terms_use.jsp
Vol. 134, No. 12
MATHEMATICAL MODEL OF A LEAD-ACID CELL 2955
i2 = I [14] e = %p [25]
~r i2 = I [26]
o~ = 0 [15]
Ox
~ $~ = 0 [27]
~ext 0520x ~g~o~ I - 0520x ~egio, 2 [16]
Equations [12] and [16] satisfy the requirement that the
flux of a species i, Ni, and the superficial current density,
i2, are continuous across the interface. In other words,
Ni,region 1 =
Ni,region 2 and
i2,regio n 1 = i2,region
2; with i = +, -, or
o. But since [Ref. (15)]
V |= ~--~ CiviY i
[17]
i
and
N i = civ i
[18]
with Vi assumed constant, the requirement that the flux
of species i is continuous implies that the volume aver-
age velocity is also continuous,
i.e., v~" = vs'.
Equation [13] describes the porosity variation with
time and Eq. [14] states that all the current at the inter-
face is in the electrolyte phase. Finally, Eq. [15] states
that the electrode solid phase potential gradient is equal
to zero at the interface because the electrode solid phase
ends there.
Region 2, reservoir.-
porosity-
Ohm's law in solution--
i~ 0r RT
+--
K
0x F
solid phase potential--
= 1 [19]
In (c3O
-- (1 - 2t~ - 0 [20]
0x
~, = 0 [21]
material balance--
0C 1 [( MWpbso 4 MWpbo2 )
Ot
2F PPbSO4 PPbO2
2t~
0c _ D 9-cc =0
(3
o J OX OX 2
current in solution--
[22]
i~ = I [23]
Equation [19] states that the reservoir is filled with
electrolyte. Equation [20] is Ohm's law applied to the
electrolyte written with respect to Pb reference elec-
trode. The electrode solid phase potential is zero (Eq.
[21]) because there is no electrode solid phase in the res-
ervoir. Equation [22] is the material balance for the elec-
trolyte in the reservoir. Finally, Eq. [23] states that all the
applied current is flowing through the solution phase
since there is no electrode solid phase in the reservoir.
Interface between region 2 and region 3.--
0c
region .C~u.~. region
D
-- CV2 g = D~sep exa
-- CV3 I
OX z 3
where
[24]
V2 m
1 1( MWpbso4 MWpbo2
2F PPDSO4 PPbO2
- (3 - 2t~ + 2VolI
1 [( MWpbso4
v3" = 2F PPbSO4
MWpb ) +(1- PPb 2t~ I
[0
(~2
RT
0 In (cf) ]
Ox
F (1 - 2t~ 0x r~g~on
=[e~opex30~)2
RT
01n(cf) ] [28]
0--x-- - e~Ptx~ F (1 - 2t~ 0~ ~gio~ 3
Equation [24] satisfies the condition that the fluxes of the
species are continuous. As shown by the expressions for
v2" and v3" below Eq. [24], v2" is not equal to v3". v2" is
equal to Vl" as stated earlier and v3" is equal to v4" which
is the volume average velocity due to migration and
changes in the structure of the lead electrode
(i.e.,
region
4). The implicit requirement that v2" = v3" is not met here.
This inconsistency would be removed in a two-dimen-
sional model which included the fact that the level of the
acid changes during charge and discharge. Previous
workers (11, 16, 17) have avoided this problem by assum-
ing that the reservoir is well mixed and, consequently,
separates the convective flows out of (or into) the porous
electrodes. Since the magnitudes of the volume average
velocities are typically small, this inconsistency is ig-
nored in the model presented here.
Equation [25] sets the porosity at this interface to be
that of the separator porosity, which is a constant value.
Equation [26] indicates that the superficial current den-
sity in the solution is equal to the applied current density
since no charge transfer reaction occurs in either the res-
ervoir or the separator. Equation [27] sets the electrode
solid phase potential to zero since there is no electrode
here. Equation [28] satisfies the condition that the cur-
rent is continuous across the interface. Unlike Eq. [16]
where the terms with the concentration gradients cancel
according to Eq. [12], the concentration gradient terms in
Eq. [28] do not cancel because they are not equal accord-
ing to Eq. [24].
Region 3, separator.--
porosity--
e = esep [29]
Ohm's lawin solution--
i2 a$2 RT
--+ -- -- (1-- 2t~ --
~sep ex3
K 0X F
0
in (cf)
- 0 [30]
Ox
solid phase potential--
material balance--
~, = 0 [31]
1 [( wpbso4 Mwpb)
oc +--
esep O~ 2F PPbSO4 {OPt)
-2t~ 1
I ac ~a D
02c
+
(1
0x -~' ~=0
current in solution--
[32]
i~ = I [33]
Equation [32] is the material balance for the electrolyte
in the separator. The second term in the equation is v3"
(Oc/Ox),
where v3" is equal to that at the interface of re-
gions 3 and 4. The porosity of the separator is fixed (Eq.
[29]), and the current density in the solution is equal to
the applied current density (Eq. [33]). The electrode solid
phase potential again is set to zero because there is no
electrodehere (Eq. [31]). Finally, Eq. [30] is Ohm's law
applied to the solution phase in the separator.
Interface between region 3 and region 4.--
~!isepeX 3
OC
~ ~ex4 C~
OX region 3 OZ region 4
[34]
Downloaded 13 Jun 2011 to 129.252.106.227. Redistribution subject to ECS license or copyright; see http://www.ecsdl.org/terms_use.jsp
2956
J. Electrochem. Soc.:
ELECTROCHEMICAL SCIENCE AND TECHNOLOGY
December 1987
O e + 1 ( MWpb) 0i2 0 [35 ]_ =
Ot
2F PPbSO4 PPb c]X
i2 = I [36]
00010x
~g~on 4 =
0
[37]
esepeX3
0(~2 eex4 0~)2 region
OX
reg[o, 3 :
OX
4
[38]
The equations for this interface were established with
the same reasoning that was used for the interface be-
tween region 1 and region 2.
Region 4, negative electrode.-
porosity variation--
0: + 1 ( . MWpb ) 0i2 0 [39
Ot
2F pvbso4 pw 0x
Ohm's law in solution--
i2 or RT
0 In (cJ)
--+-- (1-2t~ [40]
eex4K
0X F 0X
Ohm's law in solid--
04,
ig
- I~exmgo'pb- -- I = 0 [41]
Ox
material balance--
1 [( wPbso4 M pb)
oc +
Ot 2F-
PPuso4 PPb
+(1-2t~
OxO (e
ex4D
oxOC )
t - 2t~ )
Oig = 0
- (1 -
C~re)-oX
[421
~f
electrode kinetics--
for discharge and rest
0~- -- amax4?'~ \
ere f ] emax4
-- ~io4 /
{ exp [ r
RT (~b,-
r exp kl[ -aOgFRT
(~)l
-- (~2)]} = 0
[43a]
and for charge
(ct, ( ;(
0X amaxg?'~ ........
\
ref
/ ~max4 -- ~o4 ~max4 -- ~o4
[43b]
Equations [39]
through [43] are counterparts of the equa-
tions established for the positive electrode. The concep-
tual reference electrode used to measure the potential
difference is a lead electrode, the same kind as the nega-
tive electrode.
Center of the negative electrode,
x =/.--Equations [1],
[2], [3], [51, and [43a1 or [43b] apply. Since the electrode
solid phase potential was set to zero at the other bound-
ary (x = 0), a kinetic expression is used to calculate 4, at
this boundary.
Numerical Procedure
The governing equations were put into finite differ-
ence forms and solved using the numerical technique of
Newman (18, 19). The cell was divided into
NJ
node
points with J = 1 designated to be the center of the posi-
tive electrode and
J = NJ
to be the center of the negative
electrode. Each region was evenly divided, but node
point spacings were different between regions.
The finite difference approximations of the derivatives
for an internal mesh point can be written as
02C~(J) Ck(J +
1) + Ck(J - 1) - 2Ck(J)
[44]
Ox g ( ax) g
aCk(J) Ck(J +
1) -
C~(J -
1)
ox
2(5x)
and for a boundary node
0Ck(J) --Ck(J
+
2) +
4Ck(J +
1) - 3Ck(J)
[45]
[46]
ox
2(AX)
OCk(J) Ck(J -
2) -
4Ck(J -
1) + 3Ck(J)
[47]
Ox
2(hx)
where hx denotes the distance between node points and
Ck(J)
represents the k th unknown at node J.
The accuracy of the finite difference approximations of
the above derivatives is (hx) 2. Equation [46], in the for-
ward difference form, is used at J = 1; and Eq. [47], in the
backward difference form, is used at
J = NJ.
for the inter-
nal boundaries, Eq. [46] is used on the higher region
number side and Eq. [47] is used on the lower region
number side. For example, Eq. [24] written in the finite
difference form is
CI(J - 2) -
4C,(J -
1) + 3CI(J)
D - C,(J)v2"
2(hx2)
-C,(J + 2) + 4C~(J + i) - 3C,(J)
= Desep ex3 - CI(J)v3"
2(Axe)
[48]
where the electrolyte concentration c is written as "un-
known 1". The subscript on Ax refers to the region num-
ber. We note that five node points are needed to describe
the continuity of the flux at an internal boundary to
maintain
(5x) 2
accuracy. Newman's BAND(J) can be
used for only three node points (18, 19). Consequently,
since up to five node points are used at an internal
boundary, a modified version of Newman's subroutine,
called Pentadiagonal BAND(J) (20, 21), was used in-
stead. As the name implies, this subroutine allows up to
five node points to be used at any position.
For the time derivative, implicit stepping was used
OCk(J) Ck(J) -
CKk(J)
[49]
ot At
where CKk(J) refers to the value at the previous time
step,
t - At.
Initial distributions of variables can be deter-
mined by taking a small time step (10-4s,
e.g.).
Results and Discussion
In this section, we examine some simulated results
and their implications. The parameters used in the calcu-
lations are given in Table I. Parameters not referenced
are arbitrary but reasonable quantities. The same ex-
change current densities were assigned to the positive
and negative electrodes. These values were chosen so
that the calculated cell voltages agree well with mea-
sured ones.
Effect of temperature on discharge.--The
simulated
discharge behavior at 25 ~ and -18~ under a constant
current density of 340 mA/cm ~ is presented first. The dis-
charge cutoff was chosen to be 1.55V. The temperature
dependence of the kinetic expression is explicitly ex-
pressed in the exponential terms. Implicitly, the ex-
change current density in the kinetic expression is also a
strong function of temperature. As mentioned earlier,
the exchange current density was adjusted for the two
temperatures so that the calculated cell voltages agreed
relatively well with experimental observations.
Downloaded 13 Jun 2011 to 129.252.106.227. Redistribution subject to ECS license or copyright; see http://www.ecsdl.org/terms_use.jsp