Faster Lead-Acid Battery Simulations from Porous-Electrode Theory: I. Physical Model

An isothermal porous-electrode model of a discharging lead-acid battery is presented, which includes an extension of concentrated-solution theory that accounts for excluded-volume effects, local pressure variation, and a detailed microscopic water balance. The approach accounts for three typically neglected physical phenomena: convection, pressure diffusion, and variation of liquid volume with state of charge. Rescaling of the governing equations uncovers a set of fundamental dimensionless parameters that control the battery's response. Total volume change during discharge and nonuniform pressure prove to be higher-order effects in cells where variations occur in just one spatial dimension. A numerical solution is developed and exploited to predict transient cell voltages and internal concentration profiles in response to a range of C-rates. The dependence of discharge capacity on C-rate deviates substantially from Peukert's simple power law: charge capacity is concentration-limited at low C-rates, and voltage-limited at high C-rates. The model is fit to experimental data, showing good agreement.


Introduction
Lead-acid batteries are the most widely used electrochemical storage technology, with applications including car batteries and off-grid energy supply. Models can improve battery management-for example, by minimising overcharge to extend cycle life.
The most rigorous mechanistic approach to battery-cell modelling begins with a detailed microscopic description, wherein the electrolyte and electrodes occupy discrete spatial domains; volume averaging is then performed to produce a macroscopic model [1][2][3]. The details are beyond the present scope, but such a homogenization underpins the model presented below.
Our development of a detailed macrohomogeneous model of a typical lead-acid battery augments standard approaches by explicitly considering the balance of water, the variation of acid density with molarity, and the distribution of pressure. As well as ensuring thermodynamic consistency and retaining model closure when systems span multiple spatial dimensions [4], this reveals novel convective phenomena that may occur at high discharge currents. Nondimensionalization of the model helps to assess the relative importances of different multiphysical phenomena within it. This facilitates numerical implementation because it greatly improves the conditioning of the system, and hence improves the speed of computations, by making most dependent variables close to unity. This paper focuses on the development of a detailed nonisobaric physical model of battery discharge. A novel nondimensionalization is presented, which helps to identify the key dimensionless parameters that control the battery's transient response. Finally, a numerical procedure is developed to solve the nonlinear system of partial differential equations comprising the model. The results are fit against experimental discharge data and used to show how discharge capacity depends on C-rate. The nonlinear model shows that the battery response does not satisfy Peukert's law. Instead, the capacity follows one power law at low C-rates, where average acid concentration controls the response, and a different power law at at high C-rates, where overpotentials control the response.
We show below that some of the excluded-volume and pressure effects can be neglected in the detailed physical model, on the basis that certain dimensionless factors are small for the particular materials used in the battery considered here. In cases where dimensionless parameters in a model are small or large, asymptotic analysis can be employed, producing simplified, more computationally efficient models that achieve the same level of physical accuracy. The full nonlinear model presented here serves as the baseline for testing a hierarchy of computationally efficient reduced-order models, which will be developed by perturbation expansion in part II.

Unscaled governing system
Battery configuration and chemistry. Typically, a valve regulated lead-acid battery comprises six 2 V cells wired in series. Figure 1 depicts one such cell, which consists of five lead (Pb) electrodes and four lead dioxide (PbO 2 ) electrodes, sandwiched alternatingly around a porous, electrically insulating separator to produce eight electrode pairs, wired in parallel at the top edge of the electrode pile. The pile has height H, depth W , and cross-sectional area A cs = HW . The negative (Pb) and positive (PbO 2 ) electrodes have half-widths l n and l p respectively, and the separator has width l sep , giving each internal electrode pair total width L = l n + l sep + l p . The layers repeat periodically, allowing the whole pile to be modelled via analysis of one electrode pair. Some thermal effects have been documented experimentally [5,6] and considered theoretically [7,8]; we assume an isothermal system for simplicity. The electrodes are porous, permeated by an aqueous sulfuric acid (H 2 SO 4 ) electrolyte that also permeates the separator. The negative and positive half-reactions are  [7,[9][10][11] have observed that gas evolution can become important when recharging. To avoid the need for considering such side reactions, the present analysis is limited to discharges.
On any scale much larger than the Debye length, local electroneutrality holds in the liquid phase. Both ion concentrations then relate to the H 2 SO 4 molarityĉ througĥ In an isothermal state, thermodynamic consistency of the liquid volume requires thatĉ w depends onĉ alone if the bulk modulus of the liquid is very large [15], througĥ whereV e =V + +V − is the partial molar volume of the acid. Consequently the total liquid molarityĉ T = 2ĉ +ĉ w also depends only onĉ. Acid molarity further determines the mass densityρ of the liquid, becausê in which M e = M + +M − stands for the acid's molar mass. Experiments show thatρ varies almost linearly withĉ [16], so the partial molar volumes can be assumed constant. The mechanical state of the liquid is described by the external pressurep, while its mixing free energy is parameterized by a thermodynamic Darken factorχ (ĉ).
Balances. The positionx within each electrode pair traverses three subdomains: the negative electrode (0 <x < l n ), separator (l n <x < L − l p ), and positive electrode (L − l p <x < L). The model structure is identical in each subdomain. Parameters describing the electrolyte are similar everywhere, while those describing pore geometry generally differ among subdomains; quantities that parameterise reactions differ in the negative-and positiveelectrode subdomains, and vanish in the separator.
Liquid volume fractionε and reactive solid area per volume A characterize the homogenized geometry within the electrode pair. In electrode subdomains, deposition (removal) of solid PbSO 4 on pore surfaces is accompanied by removal (deposition) of solid Pb or PbO 2 , so these geometric parameters can generally vary locally. Solid PbSO 4 does not tend to deposit as a compact thin film, so mechanistic lead-acid battery models usually let A vary with state of charge [9,[11][12][13][14]. Since the functionality of this variation is disputed, we instead let the area be constant, following the approach of Liu et al., who showed this is sufficient to model Li/air-battery discharge [17]. From this perspective A describes an immobile Gibbs dividing surface that partitions the layers of Pb or PbO 2 that contribute to the battery's charge state from the currentcollecting Pb layer beneath them, which does not. Thus the time change in pore volume relates to the molar volumes of the solid species in scheme (2.1) through in whichV k is the partial molar volume of species k andŜ k is the rate at which k is generated by interfacial reactions per unit of pore area. After homogenization, the local mass balance of species k ∈ {w, −, +} in the pore-filling liquid phase implies that in whichv =V wNw +V −N− +V +N+ signifies the volumeaveraged liquid velocity. Under Faraday's law, ion balances (2.6) also combine to demonstrate charge continuity in the liquid, ∇ ·î = Aĵ. (2.8) Letting F be Faraday's constant, the liquid-phase current density isî = F (N + −N − ) and the current density associated with interfacial charge exchange isĵ = F (Ŝ + −Ŝ − ).
(Positiveĵ flows into pores.) Note that any current leaving the liquid at a given location enters the solid there. Thuŝ whereî s is the solid-phase current density. In a general isothermal setting, liquid convection is determined by a momentum balance, such as Cauchy's equation [4,15]. This governs the distribution of momentum densityρv, naturally expressed in terms of the massaveraged velocityv. The kinematic relation specifies howv must relate tov .
Flux constitutive laws. Two Onsager-Stefan-Maxwell flux laws govern transport in the liquid phase. The law for the thermodynamic force on water [15] can be inverted [18] to giveN where R is the gas constant, T is the absolute temperature, t w + is the cation transference number relative to the water velocity, andD eff is the effective diffusivity of acid in water; the pressure-diffusion factorψ depends onĉ througĥ (2.12) (To put equation (2.11) in the form given, one must have [19].) The law for the thermodynamic force on cations can be linearly transformed to produce a current/voltage relation: Hereκ eff is the effective conductivity;Φ is the potential measured by a reversible hydrogen electrode atp [20]. Effective properties appear in equations (2.11) and (2.13) because pore connectivity affects apparent transport rates. We letD eff =D(ĉ)ε 3/2 andκ eff =κ(ĉ)ε 3/2 , (2.14) following Bruggeman's tortuosity correlation. The solid phase is electronically but not ionically conductive. Thus the current density there relates to the solidphase potential,Φ s , through Ohm's law, where σ eff is the effective conductivity, Here σ eff is constant becauseε s , the volume fraction of nonreacting Pb bounded within electron-exchange surface A, does not vary. Note thatε s does not count the volume fraction occupied by solid reactants, soε s = 1 −ε. A constitutive law for liquid stress closes the model [4]. Below, terms involving momentum will be shown to be negligible in the first approximation. To analyse the scale of these terms, it will suffice to assume that current density induces a flow with low Reynolds number, in which case the homogenization of Cauchy's equation produces Darcy's law,v = −K µ∇p , (2.17) in whichμ is the liquid viscosity; pore geometry determines the Kozeny-Carman permeabilityK.
Interfacial constitutive laws. Electrons are the only solidphase charge carrier in scheme (2.1), making interfacial electronic current a proxy for the half-reaction rates. No interfacial reactions occur in the separator domain, soĵ = 0 uniformly there and the reactive area A does not need to be defined there.
In an electrode subdomain where a single half-reaction involving n e electrons occurs, the molar flux of reactants isĵ/n e F and hence everyŜ k is given bŷ where s k is the signed stoichiometric coefficient of species k in the half-reaction. (For a reduction half-reaction, s k is positive for a product and negative for a reactant.) Relationship (2.18) permits simpler expressions to be used in place of the general material balances (2.6). First, equations (2.5), (2.6), (2.7) and (2.18) combine to givê This liquid-phase volume balance introduces the volume of reaction ∆V , related to half-reaction stoichiometry by The generation term here includes a single additional parameter, Three aspects of the acid balance (2.21) are new. First, the convection term, and second, the volume of reaction ∆V , appear because state equation (2.3) imposes constraints on balances (2.6). Finally, a pressure-diffusion term appears because the flux laws are based on thermodynamic forces. To complete the model, the current densityĵ across the surface A is governed by a chemical-kinetic constitutive law. Generally such laws involve the voltage difference between the liquid and solid, the equilibrium potential of the half-reaction, and the chemical activities of the reactants. We assume the half-reactions in scheme (2.1) are elementary, following Butler-Volmer kinetics. Butler-Volmer laws naturally include a symmetry factor, which we take to be one half [13], yieldinĝ whereĵ 0 is the concentration-dependent exchange-current density,ĵ andη is the surface overpotential HereÛ stands for the half-reaction's open-circuit potential (OCP) relative to a particular reference electrode. (The reference must be the same for all half-reactions.) Terms involving interfacial capacitance C dl help to smooth out numerics, but have a negligible effect on model predictions because the capacitive time-scale is very short [21].
Kinematic relation (2.10) shows thatv · n = 0 at these boundaries too. Flux laws (2.11), (2.13), and (2.17) further require the surface-normal gradients ∂ĉ/∂n, ∂p/∂n and ∂Φ/∂n to vanish at these boundaries. Above the electrodes, there is a region of free electrolyte, with heightĥ(t). At the top surface of this region, we impose a known external pressurep ext , and the absence of flux relative to the surface, which moves with velocitŷ v head = ∂ĥ/∂t e z : Note the final condition in (2.27) determines the a priori unknown heightĥ(t).
Hereafter, subscripts n and p denote property values in the negative-and positive-electrode subdomains. We choose the negative electrode to be the ground state, and define the cell voltage to be the potential at the positive electrode tab:Φ One can either control the voltage, or consider a currentcontrolled discharge where the voltage is determined by − tabnî s · n dS = tabpî s · n dS =Î circuit (t)/8, (2.29) whereÎ circuit is the current drawn from the battery, which is positive for a discharge; the factor of 8 appears because the cell comprises eight electrode pairs in parallel. This paper focuses on experiments under 'galvanic control', following condition (2.29), which allowÎ circuit t to be any function of time. Since six cells are connected in series, the voltage in the external circuit isV circuit = 6V .
Relationships between subdomains. The liquid phase permeates all three subdomains. Therefore scalar variableŝ c,p, andΦ, as well as the normal components of all flux vectors, are continuous across electrode/separator boundaries.
There is no solid-phase current at either edge of the separator subdomain or pore-surface charge exchange within it, soî s vanishes uniformly there. Since the separator subdomain electronically isolates the positive and negative electrodes,Φ s is not continuous across it.
Integrating the interfacial current distributions in equation (2.9) and applying the divergence theorem, boundary conditions (2.29) and the fact thatî s vanishes at the electrode/separator interfaces leads to integral constraints, In short, these say that the total current leaving the negative electrode domain must enter the positive electrode domain. Expressions (2.30) will help to analyze the scales of pressure and velocity in (2.42a).
Initial conditions. The initial electrolyte concentration and electrode porosities are spatially uniform, but depend on the state of charge, q, which we define as (2.31) Let c max , ε max n and ε max p be the values of electrolyte concentration, negative electrode porosity, and positive electrode porosity, respectively, at full state of charge. In a lead-acid battery, the state of charge is closely linked to the concentration of the electrolyte. Hence q is chosen to be unity when the concentration of the electrolyte is at its maximum value, c max , and zero when the concentration of the electrolyte is zero, so that the initial conditions arê (2.32) Parameters q max and ε ∆ are chosen to make (2.31) and (2.32) consistent with (2.21) and (2.5) (see Appendix A for details): Finally, because interfacial capacitance effects are included, initial conditions are needed for the potentials; we choose spatially homogeneous values such thatĵ = 0 att = 0: (2.35b)

Dimensionless model
Nondimensionalization of the system presented in the Unscaled governing system section helps to determine the dominant effects in the system. If ∂ĥ/∂t is sufficiently small and the electrode conductivity is sufficiently high, one can assume uniformity in the plane normal tox, reducing the problem to one spatial dimension. In this case, the tabs cover the whole electrodes atx = 0 andx = L, and so the boundary condition (2.29) becomeŝ i s · e x =î cell (t) =Î circuit /8A cs atx = 0, L, (2.36) which introduces the variableî cell to stand for the total current density at the cell level. LetN + ,î,î s ,v andv represent thex-components of vectorsN + ,î,î s ,v and v. A first set of dimensionless variables is formed by the natural rescalings Reduce the number of parameters by identifying the scale of interfacial current density,ī/AL, and discharge time-scale, c max F L/ī. Hence nondimensionalize interfacial current density, exchange-current density, and time as (j n , j 0,n ) = (ĵ n ,ĵ 0,n ) ī (2.39b) QuantitiesD (and henceD eff ) andĉ w are appropriately scaled with their values at c = c max : where D max =D(c max ) and c max w =ĉ w (c max ). The conductivity and Darken thermodynamic factor rescale as where p ref is a reference pressure (e.g. atmospheric). This scaling transforms the equations governing v , v and p to ∂v ∂x = βj,  Table 1, the dimensionless parameters β and π os are small. We assume a limit where both of these parameters are zero (and hence ∂ĥ/∂t = 0), in which case v vanishes everywhere, and v and p decouple from the other variables. The full model can then be solved to find the voltage without needing the velocity and pressure. 1 Having decoupled flow velocity and pressure from the rest of the model, the following dimensionless system for c, j, ε, i, Φ, i s and Φ s results: and initial conditions Integral condition (2.30) nondimensionalizes to Composition dependences of the properties D eff , χ, κ eff , j 0 , U Pb , and U PbO2 are established through the functions listed in table A.2. The four dominant dimensionless parameters, C d , ι s , β surf and γ dl , are defined in Table 1, which also states physical interpretations and typical values. In particular, the diffusional C-rate can be written as where 8A cs c max F L is the volumetric capacity of the battery (in Ah), Q is the capacity of the battery (in Ah), and C = 8A csī /Q is the C-rate of operation. Alternatively, one can identify C d to be the ratio between the applied current scale,ī, and the scale of the liquid-phase limiting current, i L = c max D max F/L. In the Results section, q 0 will be taken to be unity (the battery starts from a fully charged state) unless explicitly stated.

Numerical solution
The system of equations (2.45) was solved numerically. Code used to solve the model and generate the results below is available publicly on GitHub [22].
To facilitate numerical solution, the model was first manipulated into a form suitable for application of the Finite Volume Method. Letting φ = Φ s − Φ and noting that i s = i cell − i, one can replace equations (2.45d) to (2.45f) with We also eliminate ∂Φ/∂x from equation (3.1) to find i as a functional of c and φ: The result is a closed system of PDEs for c, ε and φ: where j = ∂i/∂x is given by equation (3.2a) in each electrode and vanishes in the separator, with boundary conditions and initial conditions derived from (2.45j). Equation system (3.2) is solved by discretising the spatial domain using Finite Volumes, choosing the spatial discretisation to be uniform within each subdomain and as uniform as possible across subdomains The results are robust to the total number of points chosen.
Having discretised in space, the resulting system of transient ordinary differential equations is solved using  (Table A.1), for a range of C-rates.
scipy.integrate in Python [23]. Finally, Φ is calculated as where the final term comes from the fact that Φ = −φ at x = 0. The voltage drop across the electrode pair is computed with Crucially, the system of partial differential equations (3.2) is much easier to solve than the differential-algebraic system (2.45). Figure 2 shows the (dimensional) voltage against capacity used for a range of C-rates. As is to be expected, the total capacity available decreases as the C-rate increases. This dependence is further elucidated by exploring the total capacity of the battery for constant-current discharges across a higher range of C-rates, summarized on Figure 3. By Peukert's Law, one would expect the graph to be linear on this log-log axis. However, in this case deviations from linearity occur because there are two distinct capacitylimiting mechanisms. At low C-rates (below 1C), the battery capacity is concentration-limited, with full discharge occurring when the electrolyte concentration reaches zero. At high C-rates, the battery is voltage-limited, because the cut-off voltage of 10.5V is achieved before the bulk concentration falls to zero.

Results
The model also allows internal variables to be explored, particularly the local water concentration. Electrolyte and viscous pressure scale osmotic pressure scale 3.6 × 10 −5 C Table 1: Dimensionless parameters, relative to the C-rate, C = 8A csī /Q. C d is the diffusional C-rate.   (Table A.1).
water concentrations are depicted in Figure 4. At a very low C-rate of 0.1C (Figure 4a), both concentrations remain almost uniform throughout the discharge; electrolyte concentration decreases, and water concentration increases proportionally. At a higher C-rate of 0.5C ( Figure 4b) and 2C (Figure 4c), the concentrations become spatially inhomogeneous, which leads to concentration overpotentials that limit the accessible capacity. A Derivative-Free Gauss-Newton method [24] was also used to fit the model to data from a series of constantcurrent discharges of a 17 Ah BBOXX Solar Home battery at intervals of 0.5 A from 3 A to 0.5 A ( Figure 5). Each constant-current discharge is followed by a two-hour rest period during which the current is zero. The fit is good, but this approach is slow since it requires solving the full PDE system at each iteration. A faster approach will be developed in part II.

Conclusions
Three novel phenomena were included in a porouselectrode model for lead-acid batteries. First, the massaveraged and volume-averaged velocities of the electrolyte were both considered, the former associated with momentum transport, and the latter with kinematics. Due to density variation in the liquid, and volume changes associated with the electrode reactions, neither of these velocities   remains solenoidal after homogenisation. Second, an extra convective term, associated with the volume-averaged velocity, drives cation transport. Third, a pressure-diffusion term drives cation transport, and appears also in the MacInnes equation (modified Ohm's law) describing the liquid. Although these terms are small in magnitude for lead-acid batteries, they could be important for other chemistries where large volume changes occur during charge/discharge, such as lithium-ion batteries with silicon anodes [25][26][27][28][29][30], or to understand the impedance signature of electromechanical/transport coupling [31]. Nondimensionalisation of this model allows us to identify key parameter groupings, and could also easily be extended to other models and chemistries. Two distinct mechanisms determine the total cell capacity: concentration limitation at low C-rates, and voltage limitation at high C-rates. In addition to capacitylimitation and voltage-limitation, there may be additional capacity limitations from additional physical effects, such as pore occlusion at very high C-rates [17].
The model developed in this paper provides physical detail about the electrochemical processes occurring in a lead-acid battery during discharge, but is ultimately too computationally intensive to be used for advanced battery management systems. In part II, asymptotic analysis is used to derive three simplified models valid at low-tomoderate discharge rates. These can be solved much faster than the detailed model developed here, while giving additional physical insights.

Acknowledgements
This publication is based on work supported by the EP-SRC center For Doctoral Training in Industrially Focused Mathematical Modelling (EP/L015803/1) in collaboration with BBOXX. JC, CP, DH and CM acknowledge funding from the Faraday Institution (EP/S003053/1).

List of symbols
Variables The dimensional parameters are given in Table A.1, and the concentration dependences of coefficients are laid out in Table A.2. Formulae for open-circuit potentials were obtained empirically by Bode [32]. Note that these could be written aŝ to reflect their relation to the Nernst equation (e.g. [33]), with U Pb = −0.295 and U PbO2 = 1.628.
Consistent initial conditions. For consistent initial conditions, we take 'starting at x% SOC' to mean an internally equilibrated (i.e. spatially homogeneous) initial state at uniform concentration corresponding to this charge. The following determines initial conditions forĉ andε consistent with this choice. Considering a one-dimensional model and integrating      [14], citing [36], and agree with data of [16]. (*) Our fit to data of given reference(s). ( ‡) m(ĉ) =ĉV w /[(1 −ĉV e )M w ].