Fundamental governing equations

CHAPTER 3

Fundamental governing equations Contents 3.1. Porous electrode theory and volume-averaging technique 3.1.1 Modeling domain 3.2. Governing equations 3.2.1 Conservation of mass and momentum 3.2.2 Conservation of species 3.2.3 Conservation of electrical charge 3.2.4 Energy balance 3.3. Volume-averaged governing equations 3.3.1 Conservation of mass 3.3.2 Conservation of momentum 3.3.3 Conservation of species 3.3.4 Conservation of electrical charge 3.3.5 Energy balance 3.4. Microscopic modeling 3.4.1 Specific active area 3.4.2 Species diffusion length 3.4.3 Microscopic Ohmic resistance 3.5. Summary 3.6. Problems

83 84 87 88 89 90 91 93 93 95 96 101 103 107 107 108 109 110 111

3.1 Porous electrode theory and volume-averaging technique Modern high-energy batteries are made of porous electrodes. Technically, porous electrodes contain more energy in comparison with solid flat ones since they have a more specific active area. Since electrochemical reactions take place at solid surfaces, increasing specific area (available surface area per volume) means a more intensive reaction in the same physical volume. From the simulation point of view, the theory of porous media should be applied for modeling and further simulating the electrodes. According to porous media modeling concepts, all the governing parameters should be averaged with volume as the main weight function. This means that the governing equation of battery dynamics should be volume-averaged. Simulation of Battery Systems https://doi.org/10.1016/B978-0-12-816212-5.00007-6

Copyright © 2020 Elsevier Inc. All rights reserved.

83

84

Simulation of Battery Systems

Figure 3.1 A typical lead–acid cell [modified from Ref. [27]]. (A) Battery section view. (B) Cell top view.

3.1.1 Modeling domain In Fig. 3.1, a typical lead–acid battery is shown. A cross section of a single cell can be seen in Fig. 3.1A. As can be seen, the cell consists of some positive and negative electrodes with proper separators in between. The top view of the cell shown in Fig. 3.1B clearly indicates that each positive electrode is surrounded by two negative electrodes. Consequently, each positive electrode makes two electrochemical cells by its two neighbors. The same argument is true for negative electrodes. In this point of view, we can imagine that each electrode can be divided into two half electrodes, each of which works against its own counter half-electrode. By this definition, for each cell set, there are many electrochemical cells that start from the center of the positive electrode to the center of the negative electrode. The above illustration implies that for each electrode, except for two boundary electrodes, a symmetric condition is valid. This assumption makes a negligible error in computation; as it is clear from the figures, there are normally many electrochemical cells inside each single cell set. The only exceptions are the two boundary electrodes for which the symmetric condition is not valid. Although this assumption introduces some errors in modeling, the amount of the error is negligible. As mentioned before, the symmetric condition suggests that the proper domain for modeling and simulation is the area starting from the center

Fundamental governing equations

85

Figure 3.2 A typical electrochemical cell model.

of the positive electrode and ends at the center of the negative electrode. Fig. 3.2 shows the domain of a typical cell model. The domain of interest, as discussed before, begins at the center of the positive electrode and ends at the center of the negative electrode. A separator is sandwiched between the electrodes. The separator may contain some ribs (as is typical in the lead–acid batteries shown in Fig. 3.1B) leaving a bulk volume of electrolyte next to one or both electrodes. In Fig. 3.2 the bulk electrolyte exists only near the positive electrode, but in some cases, the free electrolyte may also exist at the negative side. In the porous electrode model, each part of the cell is called a region. For example, the positive electrode, separator, and negative electrode in Fig. 3.2 are the regions of the cell. The regions are made of porous electrodes and separator that are filled with a specific electrolyte. In this figure, just a cross section of the cell is shown. Although it seems that the solid parts are not connected, they are actually connected in other sections. The current collector is located at the center of the positive electrode, and the other current collector is located at the center of the negative electrode. The positive electrode is made of positive active materials (PAM), and the negative electrode is made of negative active material (NAM). Concerning any porous media, the governing equations should be obtained according to the laws of porous media. In porous media theory, all governing equations should be volume-averaged inside the porous medium. A typical porous medium is shown in Fig. 3.3, which contains solid, gas, and electrolyte phases. This volume could be a representative electrode volume (REV) used for volume-averaging purposes. It is clear that each property (i.e., such as potential, concentration, etc.) is not uniform across

86

Simulation of Battery Systems

Figure 3.3 Illustration of a sample volume containing solid, gas, and electrolyte phase.

the volume. For example, the potential on solid phase differs from each solid particle in the sample REV. If the REV is considered to be small enough, then we can assign a unique solid-phase potential to the whole REV. This unique value is obtained from averaging the solid-phase potential all over the REV. However, since the averaging depends on the volume of the solid phase with respect to the volume of the REV, the averaging should be a volume-averaging. This fact is true for any characteristic parameter, which we denote by  . Assume that V◦ is the volume of the REV in Fig. 3.3 with solid phase (k = s), electrolyte phase (k = e), and even gas phase (k = g). General volume averaging theory for any property k in phase k is defined by the following equations: 

  ∂k ∂ k  1 − k w  k · nk dA = ∂t ∂t V◦ Ak

(3.1)

for temporal averaging and ∇k  = ∇ k  +

1 V◦

 k n k dA

(3.2)

Ak

for spatial averaging. In these equations, w k is the velocity of the boundary of the phase k, and nk is the normal unit vector pointing outward the phase boundaries. The integral in these equations is taken over all the phase surface areas or Ak that is in contact with other phases such as m; hence, the area is defined

Fundamental governing equations

as Ak =



Akm ,

m = k.

87

(3.3)

m

Finally, we define two intrinsic averaging operators, volume averaging and phase volume averaging, as 

1 Xk k d∀, V◦ V◦  1 k k = Xk k d∀. Vk V◦ k  =

(3.4) (3.5)

In these operators, Xk is a phase operator equal to unity inside phase k and zero otherwise. Moreover, Vk is defined as the volume of phase k inside the volume of REV or V◦ . By this definition the volume averaging and phase volume averaging operators are related to each other by the relation k  = εk k k ,

(3.6)

where εk is the volume fraction of phase k in V◦ . In other words, εk is the fraction of volume occupied by phase k and is called the phase porosity: εk =

Vk . V◦

(3.7)

Eqs. (3.1) and (3.2) provide necessary tools for obtaining volumeaveraged governing equations out of normal equations. In this chapter, we first introduce the governing equations of battery dynamics. Then using the mentioned volume-averaging technique, all the governing equations are averaged, which gives us the proper form of the equations, which is useful for numerical simulation.

3.2 Governing equations Electrochemical reactions occur at the electrode/electrolyte interface. For each electrode, there can be more than one electrochemical reaction, either primary or side reactions. In general, we can describe each reaction by the equation 

species p=1

z

sp Mp p = nj e− ,

(3.8)

88

Simulation of Battery Systems

where the summation takes place over active species in reaction j for which sp is the stoichiometric coefficient of the active species Mp , and nj is the charge number or number of electrons that participates in the reaction. According to the kinetic rate, an electrical charge is produced or consumed by electrochemical reactions with rate rse , which can be obtained from the equation rse = −

  sj j

nj F



(3.9)

ij ,

where sj is the stoichiometric coefficient for the reaction, F is Faraday’s constant, and ij is the reaction rate described by the Butler–Volmer equation      αaj F αcj F ij = i◦j exp ηj − exp − ηj .

RT

RT

(3.10)

The subscript se in Eq. (3.9) indicates that the source term should be considered at electrode/electrolyte interface. It should be noted that bulk values differ from their corresponding values at the reaction site, due to the available sources or sinks. In the Butler–Volmer equation (3.10), anodic and cathodic charge transfer coefficients are αa and αc , respectively, which are correlated to each other by the relation αaj + αcj = n.

(3.11)

In Eq. (3.10) the overpotential ηj of reaction j is defined as ηj = φs − φe − Uj ,

(3.12)

which indicates the driven force by which the electrochemical reactions take place.

3.2.1 Conservation of mass and momentum The Electrolyte and active materials of many battery technologies are liquid or gas, whose movement can be described by fluid mechanics laws. These laws are the conservation of mass or continuity equation for each phase k, which is shown in the equation ∂ρk + ∇ · (ρk vk ) = 0, ∂t

(3.13)

Fundamental governing equations

89

and the equation of conservation of linear momentum, ∂ (ρk vk ) + ∇ · (ρk vk vk ) = −∇ pk + ∇ · τk + Bk . ∂t

(3.14)

In these equations, the subscript k refers to phase k, and τ and B are the stress tensor and body force, respectively. Full discussions about fluid motion can be obtained from various classical fluid mechanics textbooks.

3.2.2 Conservation of species In an electrochemical cell, various species exist, some of which are neutral, and some others are reactive. The concentration of these species is crucial for battery simulation. According to Fick’s law, the following equation is valid for each species in phase k: ∂ ck  k, = −∇ · N ∂t

(3.15)

where c is the molar concentration of that species, and N describes the species flux vector. It is worth mentioning that in an interface the species may jump from one level to another level. In an interface like k − m, we have the following relation between the concentration of each species:  k − ck w  m − cm w (N  k ) · n k + (N  m ) · nm = −rkm .

(3.16)

This equation clearly indicates that across an interface, species may jump  and/or interface movement w.  To use due to generation rkm , species flux N,  and interthe equation, we need to pay attention to the velocity vector w, face normal vector n must be considered outward from the phase boundary. Finally, we should emphasize that the generation rate rkm not only depends on electrochemical reactions at k − m interface, but also on physical phenomena such as evaporation, solidification, and so on contribute to its value.  k . In an electroAn important part of Eq. (3.15) is flux of species or N chemical cell, in general, species transport depends on diffusion, migration, and convection and can be written as  k = − D k ∇ ck + N

tk  ik + ck vk , zF

(3.17)

where Dk and tk are the diffusion coefficient and transference number of that species in phase k, respectively. It is clear that the transference number

90

Simulation of Battery Systems

is calculated with respect to a specific velocity. Although different references can be defined for the velocity, it is better to take the averaged mass velocity into account. By this manner the velocity field obtained from Eq. (3.14) is used in all the equations and vastly simplifies numerical calculations. In solid electrolytes the velocity field is zero, and mass transfer occurs only via diffusion and migration. Finally, it is worth mentioning that mass transfer by migration is almost zero for solid species. In other words, ts ∼ 0.

(3.18)

3.2.3 Conservation of electrical charge Electrochemical reactions occur at active material and electrolyte interface; consequently, no electrical charge is produced or consumed inside any phase k. This indicates that the net electrical charge production inside a single phase is zero and can be mathematically expressed as ∇ · ik = 0.

(3.19)

This equation is called electroneutrality, which is valid inside any phase. However, electrochemical reactions contribute to the production or consumption of electrons governed by the equation ik · nk = −



inj .

(3.20)

j

In these equations, ik is the current density in phase k. For the solid phase, the electrical charge is transferred by electrons whose transport is governed by Ohm’s law: is = −σ ∇φs ,

(3.21)

where is and φs are the current density and potential in the solid phase, respectively, and σ is the electrical conductivity of the phase. On the other hand, in the electrolyte phase, no matter if it is solid or liquid, electrical charge is transferred by ions. Ion transfer takes place by two different mechanisms: 1. Diffusion of ions from more or less concentrated locations governed by Fick’s law. 2. Migration of ions due to the electrical field. Consequently, Ohm’s law in the electrolyte phase should be modified to govern both effects. This form is sometimes is called the modified Ohm’s

Fundamental governing equations

91

law and is frequently used for simulation of current density for electrolyte: ie = −k∇φe − kD ∇(ln cei ),

i = + or − .

(3.22)

In this equation, k and kD are the effective electrolyte conductivity and diffusion conductivity, respectively. The diffusion conductivity accounts for the share of current density due to the electrolyte concentration gradient.

3.2.4 Energy balance Energy is described by the general equation  ρk cpk

  ∂ Tk

k ∇ · Jk , + vk ∇ Tk = −∇ · qk + H ∂t species

(3.23)

where ρ and cp are the density and specific heat at constant pressure respectively, T is the temperature, v is the velocity vector, q is the heat flux

vector, J is the molar flux of species due to diffusion and migration, and H is the partial molar enthalpy of species. In all these parameters, the subscript k refers to phase k. By these definitions the second term at the right-hand side of the equation determines the heat transfer due to the diffusion and migration of species. This is why the summation takes place over all species. In general, the heat flux vector contains the effect of Fourier conductive heat transfer, heat flux due to species diffusion, and Dufour heat flux, which is the energy flux due to a mass concentration gradient occurring as a coupled effect of irreversible processes. Usually, the Dufour heat flux can be neglected in battery systems; hence we can describe the heat flux vector q as qk = −λk ∇ Tk +



k Jk , H

(3.24)

species

where λk is the heat transfer coefficient in phase k. Combining Eq. (3.24) and continuity in phase k, that is, ∇ · vk = 0,

(3.25)

  ∂ Tk

k . + ∇ · (vk Tk ) = ∇ · (λk ∇ Tk ) − Jk · ∇ H ∂t species

(3.26)

we obtain:  ρk cpk

92

Simulation of Battery Systems

From the thermodynamic relation

k = μk − T H

∂μk ∂T

(3.27)

and the definition of chemical potential, μk = μ0k + RT ln

ak

ak,ref

+ zF φk ,

(3.28)

where μ0k and ak respectively are the standard chemical potential and species activity in phase k, and ak,ref is species activity at standard conditions, we obtain ⎡



ak

⎤

⎢ ∂ ln ak,ref ⎢

k = −R ∇ ⎢ T 2 ∇H ⎢ ∂T ⎣

⎥ ⎥ ⎥ + zF ∇(φk − T ∂φk ). ⎥ ∂T ⎦

(3.29)

Eq. (3.29) is valid without any assumption, but in practical cases, some appropriate assumptions simplify the equation. For instance, the first term at the right-hand side is mixing enthalpy and is neglected in practical applications. Furthermore, we can ignore the dependency of phase potential on temperature. Consequently, Eq. (3.29) reduces to

k = zF ∇φk . ∇H

(3.30)

Substituting Eq. (3.30) into Eq. (3.26) yields  ρk cpk

  ∂ Tk + ∇ · (vk Tk ) = ∇ · (λk ∇ Tk ) − zF Jk · ∇φk . ∂t species

(3.31)

By the electro-neutrality assumption the current density in phase k takes place by diffusion and migration. In other words, ik =



zF Jk .

(3.32)

species

Therefore Eq. (3.31) can be rewritten as  ρk cpk

 ∂ Tk + ∇ · (vk Tk ) = ∇ · (λk ∇ Tk ) − ik · ∇φk . ∂t

(3.33)

Fundamental governing equations

93

The second term on the right-hand side is called the Joule heating and acts as a heating source term in electrical systems such as batteries. At any electrode/electrolyte interface the balance of energy becomes λe ∇ Te · ne + λs ∇ Ts · n s = ¯in + ¯in η.

(3.34)

This is because electrochemical reactions occur at the interface. In this equation, n is the phase normal vector pointing outward the phase, and subscripts s and e indicate the solid and electrolyte phase. Moreover, in is the current density due to electrode reactions. The right-hand side of the equation accounts for the generated heat due to electrode reactions. The generated heat is divided into the reversible and irreversible parts. The first term on the right-hand side indicates the reversible part and changes sign when the battery is charged or discharged, whereas the second term accounts for the irreversible part and is always positive. Summing up, Eq. (3.31) with interface balance shown by Eq. (3.34) constitutes a system of equations for the balance of energy. However, this form of energy balance is not suitable for mathematical simulation and like other governing equations, should be volume-averaged.

3.3 Volume-averaged governing equations As discussed before, the governing equations obtained do not apply to porous media. Since almost all electrochemical batteries are made of porous electrodes and separators, we have to apply volume-averaging technique to these governing equations to obtain their volume-averaged forms. The volume-averaged form applies to porous electrodes and will be used in other chapters to simulate different batteries. In the following sections of the chapter, for generality, we discuss all governing equations. For each specific battery technology, these equations should be simplified according to its own specific physical phenomena. For example, for solid electrolyte batteries (such as lithium-ion batteries), the conservation of momentum is not used; for low-current batteries such as zinc–silveroxide button batteries, we may drop energy equation since the battery does not get warm during operation.

3.3.1 Conservation of mass Applying Eqs. (3.1) and (3.2) to Eq. (3.13) results in

94

Simulation of Battery Systems

 ∂(εk ρk ) + ∇ · (εk ρk vk k ) =

km , ∂t m

where 1

km = V◦

m = k,

(3.35)

 ρk (w  k − vk ) · n k dA.

(3.36)

Akm

In these equations, km is the rate of phase change from phase m to k at electrode/electrolyte interface. Using the mean theorem of integration, we can replace the integral of Eq. (3.36) with a production of a mean flux and an average surface area. Therefore Eq. (3.36) becomes

km = akm ρk w¯ nkm ,

(3.37)

where akm is called the specific interface area at the interface between phase m-k and is obtained by the equation akm = Akm /V◦ .

(3.38)

Also, V◦ is the volume of REV with unit of cm2 cm−3 , and w¯ nkm is the normal velocity of the interface of phases m − k with respect to phase k pointing outward phase k. Concerning porous electrodes, w¯ nkx or interface velocity is the consequence of phase or volume change due to all electrochemical reactions. Just for an example, we can refer to lead–acid batteries in which porosity of the electrode changes during discharge because porous electrodes (Pb in negative and PbO2 in positive electrodes) convert to lead sulfate, which occupies more space. Similarly in Ni−Cd batteries, cadmium converts to Cd(OH)2 , causing volume change resulting in surface velocity. Another example is hydration and dehydration of metals in Ni−Cd batteries. At the electrode/electrolyte interface, the interface velocity is mostly due to the structure change of electrodes, which is a direct result of species reactions. The reactions cause the species conversion with different molecular weight and density. Therefore, a different amount of space is occupied during charge and discharge resulting in interface movement, and hence w¯ nse =

   sj j

species

nj F

 ¯inj V¯ s .

(3.39)

Fundamental governing equations

95

In Eq. (3.39), ¯inj is the exchange current density of reaction j obtained from the Butler–Volmer equation (3.10) with overpotential defined by η¯ j = φ¯ se − φ¯ es − Uj .

(3.40)

Also, the partial molar volume of a species in the solid phase is denoted by V¯ s . Combining Eq. (3.39) with Eq. (3.37) results in

se = ase ρs

   sj j

species

nj F

 inj V¯ s .

(3.41)

Applying Eq. (3.35) in solid phase (k = s) and knowing that the solid surface velocity is equal to zero result in     sj ∂(εs ρs ) inj V¯ s = aes ρs ∂t nj F j species

(3.42)

with the mean partial molar volume V¯ s = MW /ρ,

(3.43)

where MW and ρ are the molecular weight and density of species, respectively. Since the porosity of an electrode ε is related to the solid phase porosity εs by the relation ε = (1 − εs ),

(3.44)

we can use Eq. (3.42) for calculation of electrode porosity change.

3.3.2 Conservation of momentum Volume-averaging of momentum equation results in     k    ∂  εk ρk vk k + ∇ · εk ρk vk k vk k = −εk ∇ pk + ∇ · τk  + τkt ∂t   d

Mkm + Mkm (3.45) , + εk Bk k + m

where   t  τk = (vk − vk k )(vk − vk k ) ,  1 d = τk · n k dA, Mkm

V◦

Akm

(3.46) (3.47)

96

Simulation of Battery Systems

1 Mkm = V◦

 ρk vk (w  k − vk ) · n k dA.

(3.48)

Akm

 

d

and Mkm In these equations, τkt is the dispersion shear stress tensor, Mkm respectively are the momentum transfer rates due to viscous drag and interface movement. The macroscopic shear and dispersion stress in Eq. (3.45) can be simplified by using an effective mean viscosity:

  τk  + τkt = μ∗k εk ∇ vk k ,

(3.49)

where μ∗ is the total transport property containing both effects of microscopic property and the effect of microstructures in the porous medium. Since the velocity inside the porous electrodes is very small, for the first estimation, we may choose μ∗ equal to its microscopic value, neglecting the effect of velocity dispersion due to the microstructures of porous media. In porous media problems, it is customary to use the modified Darcy law for accounting the effect of interface drag, which includes both effects d

of Mkm and Mkm . The modified Darcy law in this case is 



d

Mkm + Mkm = −εk2

m

μk

Kkrk

· vk k .

(3.50)

In this equation, K is the absolute permeability, and krk is the relative permeability of phase k, including the effect of decreasing of the area due to fluid vortexes that occupy the pores. The relative permeability is a function of the volumetric phase ratio and is obtained from experimental data [28]. Substituting Eqs. (3.49) and (3.50) into Eq. (3.45), the macroscopic equation of momentum can be obtained as       k ∂  εk ρk vk k + ∇ · εk ρk vk k vk k = −εk ∇ pk + ∇ · μ∗k εk ∇ vk k ∂t μk · vk k . (3.51) + εk Bk k − εk2

Kkrk

3.3.3 Conservation of species Applying volume-averaging technique to Eq. (3.15) results in   ∂ ck  1 k − = −∇ N ∂t V◦

 Ak

 k − ck w (N  k ) · nk dA.

(3.52)

Fundamental governing equations

97

In comparison   with Eq. (3.17), the volume-averaging of species flux in  k can be determined by the equation phase k or N 



 k = − Dk ∇ ck  + N

tk   ik + ck vk  , zF

(3.53)

where the coefficient tk /zF is the result of volume averaging. This value is almost constant in volume V◦ . From the physical point of view, the first term in the right-hand side of Eq. (3.53) determines the macroscopic diffusion of species, which normally is obtained using the effective mean diffusion coefficient: Dk ∇ ck  = Dkeff ∇ ck k ,

(3.54)

where Dkeff is the effective diffusion coefficient. Since the electrolyte is confined in porous electrodes and separator, this coefficient is modified so that the effect of porosity and tortuosity of the porous media is considered. The convective term in Eq. (3.53) (the third term in the right-hand side) normally is replaced by the product of the mean concentration and mean velocity plus an additional term, which accounts the effect of hydrodynamic dispersion causing from electrolyte movement inside the microscopic pores. In mathematical form the convective term becomes ck vk  = εk ck k vk k − Da ∇ ck k ,

(3.55)

where Da is called the dispersion coefficient, which includes the axial dispersion caused by electrolyte movement near the porous electrode walls. The dispersion coefficient is not a fundamental property of fluid and depends on fluid movement. When the electrolyte is stagnant (e.g., in solid electrolytes), this term is zero. Substituting Eqs. (3.54) and (3.55) into Eq. (3.53), we obtain 







 k = − Dkeff + Da ∇ ck k + N

tk   ik + ε ck k vk k . zF

(3.56)

Hence Eq. (3.52) can be rewritten as     ∂(εk ck k ) + ∇ · εk ck k vk k = ∇ · (Dkeff + Da )∇ ck k ∂ t      tk   1 1 tk  −∇ · ik + Dk ∇ ck · nk dA − ik · nk dA +

1 V◦



zF

V◦

Ak

ck (w k − vk ) · nk dA. Ak

V◦

Ak

zF

(3.57)

98

Simulation of Battery Systems

By applying the averaging technique to Eq. (3.19), Eq. (3.57) can be further simplified. The result is    1 ik · nk dA = 0.  ∇ · ik +

V◦

(3.58)

Ak

Examining this equation reveals that a part of right-hand side of Eq. (3.57) is canceled out by the fourth term. Hence this equation reduces to ∂(εk ck k ) + ∇ · (εk ck k vk k ) = ∇ · [(Dkeff + Da )∇ ck k ] ∂t      tk d

+ (Jkm + Jkm ) − ik · ∇ ,

zF

m

(3.59)

where 

1 = Dk ∇ ck · nk dA, V◦ Akm  1

= ck (w k − vk ) · nk dA. Jkm V◦ Akm

d Jkm

(3.60) (3.61)

d

and Jkm represent the surface transport of a species in The terms Jkm phase k due to microscopic diffusion and surface movement, respectively. Eq. (3.59) determines the conservation of species in solid and liquid electrolytes. If transport equation of surface species in Eqs. (3.60) and (3.61) are accurately modeled, then we can use this equation to determine the macroscopic equations from microscopic ones. For this purpose, by considering the mean value theorem in integration the transport of species term due to the surface movement can be replaced by the product of the mean surface concentration and mass flux at the interface:

= c¯km km . Jkm

(3.62)

Similarly, the integral of Eq. (3.60) can be replaced by the product of an interfacial specific area and a mean diffusion flux at the interface. From physical point of view, the mass transfer term at interface represents the diffusion phenomena by which concentration gradient occurs. The diffusion flux is proportional to its own driving force, which is the concentration

Fundamental governing equations

99

Figure 3.4 Illustration of species diffusion length and Ohmic drop at an interface [29].

gradient. On the other hand, the flux is inversely proportional to the diffusion length l. The diffusion length determines the resistance to diffusion; hence we have d Jkm

 ∂ ck  = akm Dk  ∂ n k 

= akm Dk km

¯ckm − ck k

lkm

.

(3.63)

Hence mathematically, the diffusion length of species in phase k can be written as lkm =

¯ckm − ck k  . ∂ ck  −  ∂ n k 

(3.64)

km

In these equations, Dk is the diffusion coefficient of species in phase k, and ¯ckm is the mean concentration of species at m − k interface. Fig. 3.4 schematically illustrates the microscopic species distribution at solid and electrolyte interface. The interface shown in this figure is just a tiny portion of what is shown in Fig. 3.3. In this figure the definition of diffusion length is illustrated. The diffusion length defined in Eq. (3.64) is a very complex function of microscopic phenomena, where its determination requires specific atten-

100

Simulation of Battery Systems

tion. In the following sections, we discuss this characteristic length in more detail. The species transport equations are correlated to each other by means of surfaced balancing. To carry out the balance, Eq. (3.16) is integrated at the interface of phases k and m inside volume V◦ , which results in 

1 tk  (−Dk ∇ ck + ik + ck vk − ck w k ) · nk dA+ V◦ Akm zF  1 tm  (−Dm ∇ cm + im + cm vm − cm w m ) · nm dA = −akm r¯km , V◦ Akm zF

(3.65)

where rrm is the mean reaction rate at Akm interface. By Eqs. (3.60) and (3.61) and the relations given by Eq. (3.21), Eq. (3.65) can be rewritten as 





d

d

Jkm + Jkm + Jmk + Jmk



⎛

⎞

tk + tm  ¯ ⎠ inj . = akm ⎝r¯km − zF j

(3.66)

In this equation, two important case are to be considered: 1. When the m − k interface is electrochemically active at which ions exist only at one side (e.g., in phase k). Therefore Eq. (3.66) can be simplified as 

d

Jkm + Jkm



⎡ ⎤    tk s j ¯inj ⎦ . + = −akm ⎣ j

zF

nj F

(3.67)

To obtain this equation, Eq. (3.9) is averaged over the surface area. This equation not only can be used directly, but also can be used to obtain surface species concentration ¯ckm . To do this, we should use the constitute equations defined by Eqs. (3.62) and (3.63) in Eq. (3.67), which results in ⎡ ⎤      tk s akm Dk akm Dk j ¯inj ⎦ , (3.68) ck k − akm ⎣ ¯ckm km + + =

lkm

lkm

j

zF

nj F

where the phase change rate km is obtained from Eq. (3.41). Note that Eq. (3.68) is coupled with the Butler–Volmer equation so that ¯inj can be obtained. This is because the exchange current density and reversible potential are strong functions of species surface concentration c¯km . 2. The second case is where the m − k interface is electrochemically inactive. For an example, we can refer to soluble hydrogen or oxygen at

Fundamental governing equations

101

liquid or gas phases. In this situation the balance of species is reduced to 







d

d

Jkm + Jkm + Jmk + Jmk = 0.

(3.69)

Substituting Eqs. (3.62) and (3.63) into this equation, we have akm Dk amk Dm (¯ckm − ck k ) + (¯cmk − cm m ) lkm lmk

km = , ¯cmk − c¯km

(3.70)

where the mass balance mk = − km at interface is used. From Eq. (3.70) the phase change rates ¯ckm , c¯mk at the interface can be determined.

3.3.4 Conservation of electrical charge Macroscopic conservation of electrical charge is obtained by applying volume-averaging to Eq. (3.19):    ∇ · ik − Ikm = 0,

m = k.

(3.71)

m

In this equation the interfacial current with the unit of A cm−3 is calculated by the relation  1 ik · n k dA. (3.72) Ikm = − V◦ Akm Ohm’s law can be used to relate the electrical potential and current density in the solid phase. This law in proper vector form is ∇ · (σ eff ∇ φe e ) +



Ism = 0,

m = s.

(3.73)

m

Similarly, for the electrolyte phase, Ohm’s law indicates:   e e ∇ · (keff ∇ φe e ) + ∇ · (keff Iem = 0, D ∇ ln ci ) +

m = e.

(3.74)

m

We can use the Taylor expansion rule to obtain 

e  e ln cei = ln cei ,

(3.75)

and substituting into Eq. (3.74) yields   i e ∇ · (keff ∇ φe e ) + ∇ · (keff Iem = 0, D ∇ ln ce ) + m

m = e.

(3.76)

102

Simulation of Battery Systems

The current density Ikm at electrode/electrolyte interface is Ise = −Ies = ase



¯inj .

(3.77)

j

Similarly, for any interface between phases m and k, the interface current density can be modeled as Ikm = akm

φ¯ km − φk k

Rkm

,

(3.78)

where Rkm with units of  cm−2 is defined as Rkm =

φ¯ km − φk k  . ∂φk  −σ  ∂ n k 

(3.79)

km

This is the Ohmic resistance between phases k and interface m − k. In Eq. (3.78), it is assumed that phase k is a solid phase in which Ohm’s law is valid. For liquid phases at the interface, the effect of migration should also be considered. By this assumption Ohm’s law becomes Ikm = akm

φ¯ km − φk k

Rkm

+ akm

 kD  c¯km − ck k , ¯ckm lkm

(3.80)

where phase k is liquid. By combining Eq. (3.78) or Eq. (3.80) with Eq. (3.77) we can obtain the surface potential φ¯ se or φ¯ es . These potentials are used to calculate the surface overpotential, which is necessary for calculation of the overpotential used in the Butler–Volmer equations. In practical applications the electrode conductivity of electrodes σ is usually very high. This implies that φ¯ se must be equal to its averaged value φs s . However, this argument is not true for liquid electrolytes. In other words, the conductivity of electrolytes is relatively low; hence φ¯ se and its phase mean value φe e are not equal especially in high-current applications. This means that we cannot use φe e instead of φ¯ se to calculate the overpotential used in the Butler–Volmer equation. The relation between surface potential and volume-averaged potential is illustrated in Fig. 3.4.

Fundamental governing equations

103

3.3.5 Energy balance Volume averaging from energy equation (3.23) results in # ρk cpk

$   ∂εk Tk k k k k   + ∇ · (εk Tk  vk  ) = ∇ · (λeff + λ )∇ T , a ,k k k ∂T   Joule d

Qkm + Qkm (3.81) − ik  · ∇ φk k + Qk , m

where 

1 λk ∇ Tk · nk dA, V◦ Akm  1

= ρk cpk Tk (w  k − vk ) · nk dA, Qkm V◦ Akm

d Qkm =

(3.82) (3.83)

and QkJoule

 k 1  = − ik ·

1 − V◦



V◦

(φk − φk k )nk dA

Akm

  ik − ik k · ∇(φk − φt k )d∀.

(3.84)

Vk

In these equations, λeff k is the effective thermal conductive coefficient in phase k, and λa,k is called the dispersion coefficient in phase k. Like other properties, λeff k should be modified so that the effect of porosity and its tortuosity is included. Normally, these effects are included using the Bruggman relation [30] 1.5 λeff k = λk ε k .

(3.85)

In contrast to λk , which is an intrinsic property of any material, λa,k accounts for hydrodynamic dispersion effects, which obviously depend on the velocity and temperature of the fluid. In other words, λa,k is a flow property and is zero when the electrolyte is stagnant. The second term in the right-hand side of Eq. (3.81) is the total sum

of heat transfer at the interface; Qkm indicates conductive heat transfer,

whereas Qkm is the heat transfer due to interface movement. Using the

mean value theorem in integrals, we can replace Qkm by the product of the mean velocity value and the mean temperature of the interface:

Qkm = cpk T¯ km km ,

(3.86)

104

Simulation of Battery Systems

where T¯ km is the mean surface temperature at m − k interface. The last term in the right-hand side of Eq. (3.81), QkJoule , comes from volume averaging from Eq. (3.33). This term vanishes when phase k is electrically at equilibrium or φk = φk k . This means that the electrical potential is uniform all over phase k. When the electrical conductivity is not very high or the current density is too high, the equilibrium condition is not valid, and QkJoule should be considered. Such a case happens when semiconductor electrodes are in use like NiOOH electrode with electrical conductivity about 10−5 S/cm. Summation of Eq. (3.81) over all contributing phases (i.e., solid, liquid, and gas phases) yields %

#

$& ∂εk Tk k k k ρk cpk + ∇ · (εk Tk  vk  ) ∂t k ( '  k = ∇ · (λeff k + λa,k )∇ Tk 



k

+

 1  k=m m

V◦

     ¯ ik · ∇ φk k . λk ∇ Tk · nk dA + cpk Tkm km −

Akm

k

(3.87) The second sum in the right-hand side of Eq. (3.87) comes from the electrochemical reactions at solid/electrolyte interface. From the energy point of view, to obtain the sum, we need to write heat balance at the interface, that is, λe ∇ Te · ne + λs ∇ Ts · ns = ¯in η + ¯in ,

(3.88)

where n is the normal surface unit vector pointing outward from the phase, and the subscripts e and s as usual indicate the electrolyte and solid phases, respectively. Moreover, in is the local transfer current density, which is the result of electrochemical reactions. The right-hand side of Eq. (3.88) is the overall generated heat due to electrochemical reactions that can be divided into irreversible (the first term) and reversible (the second term) parts. Torabi and Esfahanian [31,32] presented a full detailed discussion about reversible and irreversible parts of the generated heat and introduced the concept of general Joule heating because the irreversible part of the generated heat is of Joule heating type. While the irreversible part is always positive, the reversible part changes the sign in charge and discharge conditions, the fact that is very obvious

Fundamental governing equations

105

from their definition. In other words, the irreversible part is the product of the local current density in and the local overpotential η, whose values are always positive since both parameters simultaneously change the sign in charge and discharge. However, the reversible part is the product of the local current density and local Peltier coefficient j . The Peltier coefficient was first introduced by Newmann [33], where by neglecting Duffour energy flux it is related to the entropy of reaction by the relation j =

T Sj . nj F

(3.89)

Since the entropy Sj for each reaction is a known value, during charge and discharge, it is fixed, and its value does not change. Therefore, the product of the Peltier coefficient and local current density changes sign in charge and discharge resembles a reversible process. It should be strongly pointed out that the represented heat balance is carried out only for heat of electrochemical reactions and does not include the effect of any physical phenomena such as evaporation, condensation, crystallization, and so on. If we need to include these effects, an additional source therm should be added to the right-hand side of Eq. (3.88). Hence for a general formulation, we can express Eq. (3.88) as  λk ∇ Tk · nk + λm ∇ Tm · n m =

asj¯inj (ηj + j ) + (hk − hm ) km

j

akm

,

(3.90)

where km is the phase change rate from phase m to phase k, asj is the specific active surface for electrochemical reactions, and akm is the specific surface area at interface between phases k and m. Finally, h indicates the enthalpy with respect to phases k or m according to its subscript. Applying Eq. (3.90) to Eq. (3.87) results in %

#

j

  − ( ik · ∇ φk k ).

$&   ∂εk Tk k k k k ρk cpk + ∇ · (εk Tk  vk  ) = ∇ · (λeff k ∇ Tk  ) ∂t k k  ) * ¯ (hk − hm ) km + (cpk − cpm )T¯ km km + asj inj (η¯ j + j ) +



k=m m

k

(3.91)

106

Simulation of Battery Systems

When local thermal equilibrium exists, that is, Tk k = Tm m = T¯ km = T¯ mk = T ,

(3.92)

we can drop averaging symbols; hence Eq. (3.91) is simplified to ∂(ρ cp T ) + ∇ · (vT ) = ∇ · λ∇ T + q, ∂t

(3.93)

where ρ cp =



εk ρk cpk ,

(3.94)

εk ρk cpk vk k ,

(3.95)

k

v =

 k

 λ= (λeff k + λa,k ),

(3.96)

k

and the heat source is defined as q=



asj¯inj (η¯ j + j ) +



     ik · ∇ φk k h∗ km −

k=m m

j

(3.97)

k=e,s

with h∗ = (hk − hm ) + (cpk − cpm )T ,

(3.98)

the mean volumetric current density at electrolyte phase    i ik = −keff ∇ φe e − keff D ∇ ln ce i,

and at solid phase

  is = −σ eff φs s .

i = + or −,

(3.99)

(3.100)

In these equations, kD is ionic conductivity at phase k, and σ is the solid conductivity. The superscript eff like before indicates that these values are corrected to account the effect of the porous media. Moreover, note that no current flows through the gas phase (k = g), meaning that the gas acts as an insulator. Eq. (3.97) shows that heat sources due to Joule heating, electrochemical reactions, and phase change exist both in solid and electrolyte phases. Since this equation is described in vector form, we can make use of local current density to obtain the temperature distribution for systems that are not

Fundamental governing equations

107

in thermal equilibrium. Finally, Eq. (3.97) is a temporal equation, which means that by its solution the temperature variation concerning time can be obtained. The heat sources and sinks in any battery system play a vital role in determining its thermal behavior. According to the importance of source term, that is, q in Eq. (3.97), an entire chapter is dedicated for a full discussion about the issue. Hence a detailed information is available in Chapter 4.

3.4 Microscopic modeling The battery model discussed here is a macroscopic model. This means that all fundamental parameters are averaged; hence the microscopic distribution of these parameters are rendered. However, three important microscopic parameters are crucial for the evaluation of accurate surface properties in the macroscopic model. In other words, the effect of microscopic surface phenomena that plays an important role in determining battery behavior is to be included in the model. These effects can be added to the model by accurate modeling of the following three parameters: 1. Specific active area akm . 2. Diffusion length lkm . 3. Ohmic resistance Rkm . Accurate evaluation and implementation of these parameters extremely affect the accuracy of the model. These parameters are usually very hard to obtain, and their modeling requires complex analysis of microscopic tortuosity of the porous media and analysis of the diffusion of species and charge transfer at the solid interface. Here an acceptable model for considering the effect of these parameters is presented. It will be shown that by the present modeling the effect of microscopic phenomena such as solid diffusion and ohmic resistance in semiconductors can be modeled and added to the overall macroscopic model.

3.4.1 Specific active area As it was previously shown, the specific active area or electrode/electrolyte interface area is defined as ase = Ase /V◦

(3.101)

and has a crucial role in modeling interfacial properties. In practice the interface area contains much important information about the electrochem-

108

Simulation of Battery Systems

ical reactions but unfortunately are canceled in volume-averaging process. The canceled information is really important; thus their effect should be added to the model by constitute equations. The specific active area is normally obtained by adsorbing methods (like BET test) or by double-layer charging method. However, here an estimation method for calculation of specific active area is presented. The specific active area for electrodes that are made of spherical particles (like porous electrodes made of powders) having radius rs can be calculated by a◦se =

3εs 3(1 − ε◦ ) = , rs rs

(3.102)

where ε◦ is the porosity of the electrodes. Note that the solid phase porosity and electrode porosity are related such that εs = 1 − ε◦ .

3.4.2 Species diffusion length In practical cases, battery active materials are made of planar, cylindrical, or spherical particles. Each geometry has its own relation for calculation of species diffusion length. In cylindrical geometries the active materials are supposed to be uniformly pasted on a substrate. Since in making the electrodes a very thin layer of active material is used, the effect of electrode curvature is negligible, and a parabolic equation is enough for modeling of species diffusion length: cs = a◦ + a1 r + a2 r 2 ,

(3.103)

where r is the radial distance. Boundary conditions depict ∂ cs =0 ∂r cs = c¯se

at

r = r◦ ,

(3.104)

at

r = rs ,

(3.105)

where r◦ is the radius of the core, and rs is the radius of electrode/electrolyte interface. The local mean species concentration cs s is defined as cs s =

1 Vs



cs rdr = Vs



2 rs2

− r2 ◦

rs

cs rdr . r◦

(3.106)

Fundamental governing equations

109

Using these three equations, we can obtain the coefficients a◦ , a1 , and a2 in Eq. (3.103). Substituting the coefficients into Eq. (3.103) results in (rs − r )(rs − 2r◦ + r ) ¯cse − cs = . c¯se − cs s rs2 − r◦2 2 4r◦3 − r s r◦ + 2 3 3(rs + r◦ )

(3.107)

Consequently, using Eq. (3.64), we easily calculate the diffusion length: lse =

r s + r◦ rs r◦ 3r◦3 − + . 4 3(rs − r◦ ) 3(rs2 − r◦2 )

(3.108)

Using this equation when r◦ = 0, we simplify Eq. (3.108) rs lse = . 4

(3.109)

This equation shows that the species diffusion length is proportional to the species particle size and has almost the same order of magnitude. Similarly, we can calculate the species diffusion length for Cartesian and spherical cases. For planar electrodes with a layer of active materials with a half thickness of rs , rs (3.110) lse = , 3 and for spherical particles with radius of rs , rs lse = . 5

(3.111)

3.4.3 Microscopic Ohmic resistance By the definition given in Eq. (3.79), the Ohmic resistance from electrode/electrolyte and electrode/substrate interfaces to the bulk of active materials are respectively defined as Rse =

φ¯ se − φs s  , ∂φs  −σs  ∂r 

(3.112)

r =rs

φ¯ b − φs s  , Rsb = ∂φs  −σ◦  ∂r  r =r◦

(3.113)

110

Simulation of Battery Systems

where the subscript b refers to the substrate, and φs is the solid phase potential in active materials. Note that the electrical conductivity σ is a function of SoC; consequently, it is variable across the electrode surface. The symbols σ◦ and σs show the electrical conductivity of electrode/substrate and electrode/electrolyte interfaces, respectively. It is assumed in this equation that the electrical current flows radially and the assumption of electroneutrality is valid. Integration of the equation over the active material layer results in     ∂φs ∂φs = σ ·r . σ ·r ∂ r r =r◦ ∂ r r =rs

(3.114)

Two other boundary conditions are: % φs =

at at

φb φ¯ se

r =r◦ , r =rs .

(3.115)

Assuming a quadratic profile for potential distribution in the solid phase of the form φs = a◦ + a1 (r − r◦ ) + a2 (r − r◦ )2 ,

(3.116)

the coefficients a◦ , a1 , and a2 are to be obtained from Eqs. (3.114) and (3.115). After obtaining the coefficients and substituting into Eqs. (3.112) and (3.113), we calculate the Ohmic resistances: 





r s r s − r◦ rs + 3r◦ 3rs + 5r◦ + , 12 rs + r◦ ) σ◦ r◦ σs rs    r ◦ r s − r◦ 5rs + 3r◦ 3rs + r◦ + Rsb = . 12 rs + r◦ σ◦ r◦ σs rs

Rse =

(3.117) (3.118)

It is obvious that these resistances are proportional to species diffusion length and are inversely proportional to electrical conductivity, σ .

3.5 Summary In summary, Eqs. (3.42), (3.45), (3.59), (3.73) or (3.74) and (3.93) constitute a complete set of governing equations for obtaining macroscopic behavior of batteries. The solution of the system of equations gives a variation of the following unknowns: 1. Electrode porosity εk . 2. Velocity vector of electrolyte vk k .

Fundamental governing equations

Table 3.1 Governing equations at interfaces. Transport terms at interface balance at interface

km = akm ρk w¯ nkm

km + mk =0 d = a Dk (¯c − c k ) Jkm km k lkm km

= c J Jkm km km

Ism = asm Iem = aem

zF

φ¯ sm − φs s

Rsm

φ¯ em − φe

Rem

d + J ) + (J d + J ) (Jkm mk km mk ⎛ ⎞  t + t m k ¯inj ⎠ = akm ⎝r¯km −

e

111

related equation

mass species

j

Ism + Ims = 0

electrode charge

Iem + Ime = 0

electrolyte charge

3. Concentration of chemical species ck k . 4. Potential of solid and electrolyte phases φk k . 5. Temperature Tk k . Note that Eq. (3.59) should be repeated for each active species; hence the number of concentration equations and related unknowns are as many as the number of species under consideration. Besides the governing equations, boundary conditions at any interface are necessary. These conditions are tabulated in Table 3.1, as well as the balance of parameters at the same interface. Using these equations, we can obtain the necessary information used for obtaining surface properties and, consequently, the rate of reaction at the interfaces.

3.6 Problems 1. Using the balance of mass, obtain a proper relation for porosity change during charge and discharge for NiCd batteries. 2. Using the balance of mass, obtain a proper relation for porosity change during charge and discharge for Li-ion batteries. 3. In a lead–acid factory, the size of lead particles used for making the electrodes is rs ≈ 2.5 × 10−4 cm. Calculate the approximate active area versus porosity and plot the result. 4. Discuss on mass transport mechanisms in a Li-ion battery. Which mechanisms are dominant and which we can ignore. 5. In a gelled lead–acid battery the rate of heat generation is assumed to be dependent only on Joule heating. Using Eq. (3.93), give a proper relation for energy equation.