Mathematical modeling of a lithium ion battery with thermal effects in COMSOL Inc. Multiphysics (MP) software

Mathematical modeling of a lithium ion battery with thermal effects in COMSOL Inc. Multiphysics (MP) software

Journal of Power Sources 196 (2011) 5985–5989 Contents lists available at ScienceDirect Journal of Power Sources journal homepage: www.elsevier.com/...

1MB Sizes 3 Downloads 153 Views

Journal of Power Sources 196 (2011) 5985–5989

Contents lists available at ScienceDirect

Journal of Power Sources journal homepage: www.elsevier.com/locate/jpowsour

Short communication

Mathematical modeling of a lithium ion battery with thermal effects in COMSOL Inc. Multiphysics (MP) software Long Cai, Ralph E. White ∗ Department of Chemical Engineering, University of South Carolina, 301 Main Street, Columbia, SC 29208, United States

a r t i c l e

i n f o

a b s t r a c t

Article history: Received 2 February 2011 Received in revised form 7 March 2011 Accepted 8 March 2011 Available online 15 March 2011

The existing lithium ion battery model in Multiphysics (MP) software (COMSOL Inc., Palo Alto, CA) is extended to include the thermal effects. The thermal behavior of a lithium ion battery is studied during the galvanostatic discharge process with and without a pulse. © 2011 Elsevier B.V. All rights reserved.

Keywords: Lithium ion battery Thermal model Pulse discharge Temperature

1. Introduction

where i = p for the positive electrode and i = n for the negative electrode. At the center of the particle, there is no flux, that is:

The existing lithium ion battery model in COSMOL Inc. Multiphysics 3.5a is extended here by adding an energy balance and the temperature dependence of properties of the battery. This thermal model is developed based on the pseudo twodimensional (P2D) model which was described in [1,2] and a thermal, electrochemistry coupled model. The diffusion coefficient of Li ions in the solid phase and electrolyte, the reaction rate constants of the electrochemical reactions, the open circuit potentials and the thermal conductivity of the binary electrolyte depend on the temperature in the model presented here. 2. Mathematical model A schematic of a lithium ion battery is shown in Fig. 1. Generally, a lithium ion battery consists of the current collector, the positive electrode, the separator and the negative electrode. A lithiated organic solution fills the porous components and serves as the electrolyte. The material balance for the Li ions in an active solid material particle is governed by Fick’s second law in spherical coordinates: ∂cs,i ∂t

= Ds,i

1 ∂ r 2 ∂r

 r2

∂cs,i



∂r

(1)

at r = 0 :

− Ds,i

∂cs,i

=0

∂r

On the surface of the particle, the flux is equal to the consuming/producing rate of Li ions due to the electrochemical reaction occurring at the solid/liquid interface: at r = Rs,i :

− Ds,i

∂cs,i ∂r

= Ji

where Ji is the flux of lithium ions away from the surface of the spherical particles. The material balance for the binary electrolyte in the liquid phase is given by: εi

∂ci ∂ = ∂t ∂x



Deff,i

∂ci ∂x



0 + (1 − t+ )ai Ji

(2)

where i = p, s, and n and ai is the electrode surface area per unit volume of the electrode. In the separator the pore wall flux Js is equal to zero. At the two ends of the cell in the x-direction, there is no mass flux, that is:



−Deff,p

∂cp  ∂x 

x=0 ∂cn  −Deff,n ∂x 

=0 =0

x=Lp +Ls +Ln

∗ Corresponding author. Tel.: +1 803 777 3270; fax: +1 803 777 0973. E-mail address: [email protected] (R.E. White). 0378-7753/$ – see front matter © 2011 Elsevier B.V. All rights reserved. doi:10.1016/j.jpowsour.2011.03.017

At the interfaces between the positive electrode/separator and separator/negative electrode, the concentration of the binary

5986

L. Cai, R.E. White / Journal of Power Sources 196 (2011) 5985–5989

Nomenclature as c Cp cs De Ds Ecell F h Iapp J k L R Rs t T t+ U Ui,ref x ε εf  i  1 2

specific interfacial area, m2 m−3 concentration of the binary electrolyte, mol m−3 specific heat, J kg−1 K−1 concentration of lithium ion in solid, mol m−3 salt diffusion coefficient, m2 s−1 diffusion coefficient of lithium ion in solid electrode particles, m2 s−1 cell voltage, V Faraday’s constant, 96,487 C eq−1 effective heat transfer coefficient, W m−2 K−1 applied current density, A m−2 pore wall flux of lithium ions, mol m−2 s−1 electrochemical reaction rate constant, m2.5 mol−0.5 s−1 thickness of battery component, m gas constant, 8.3145 J mol−1 K−1 radius of electrode particle, m time, s temperature, K transference number of lithium ion open-circuit potential, V open circuit potential of electrode i under the reference temperature, V spatial coordinate porosity volume fraction of fillers density, kg m−3 electronic conductivity of solid matrix, S m−1 thermal conductivity, W m−1 K−1 potential in the solid phase, V potential in the electrolyte, V

Subscripts eff effective max maximum n negative electrode p positive electrode s separator

electrolyte and its flux are continuous cp |x=L− = cs |x=L+ p

p

cs |x=(Lp +Ls )− = cn |x=(Lp +Ls )+





∂cp  −Deff,p ∂x 



 x=Lp

∂cs  −Deff,s ∂x 

∂cs  = −Deff,s ∂x 

x=(Lp +Ls )−

x=Lp+



∂cn  = −Deff,n ∂x 

x=(Lp +Ls )+

The effective diffusion coefficient of the binary electrolyte in the liquid phase is corrected by porosity: Deff,i = De,i εi bruggi where the diffusion coefficient of the binary electrolyte (De,i ) is given in Ref. [3] as follows: −3

De,i = 10−4 × 10−4.43−(54/(T −229−5.0×10

ci ))−0.22×10−3 ci

(3)

The charge balance in the solid phase is governed by Ohm’s law: eff,i

∂2 1,i ∂x2

= ai FJi

(4)

where i = p and n. The effective conductivity is determined by  eff,i =  i (1 − εi − εf,i ). The specific area can be calculated by ai = (3/Rs,i )(1 − εi − εf,i ). At the interface of the current collector and the positive electrode, the charge flux is equal to the current density applied to the cell:

  = Iapp ∂x  x=0

∂1,p 

−eff,p

where Iapp is the applied current density. At the interfaces of the positive electrode/separator and the separator/negative electrode, there is no charge flux:

  =0 ∂x  x=Lp ∂1,n  −eff,n =0 ∂x 

−eff,p

∂1,p 

x=Lp +Ls

The potential of the solid phase at the right end of the cell is set to zero 1,n |x=Lp +Ls +Ln = 0, and the potential of the solid phase at x = 0, 1,p |x=0 is equal to Ecell (the cell voltage). The charge balance in the liquid phase is based on Ohm’s law and is given by: −

∂ ∂x



eff,i

∂2,i ∂x



+

0) 2RT (1 − t+ ∂ F ∂x



eff,i

∂(ln ci ) ∂x



= ai FJi

(5)

where i = p, s and n and the specific conductivity of the binary electrolyte is a function of temperature and the concentration of the binary electrolyte in the liquid phase: eff,i = i εi bruggi where the expression of the ionic conductivity for the binary electrolyte is given in Ref. [3] as follows: i = 10−4 × ci (−10.5 + 0.668 × 10−3 ci + 0.494 × 10−6 ci2 + 0.074T −1.78 × 10−5 ci T − 8.86 × 10−10 ci2 T − 6.96 × 10−5 T 2 +2.80 × 10−8 ci T 2 )

2

(6)

At the two ends of the cell, there is no charge flux in the liquid phase: −eff,p

    = 0 and −eff,n ∂2,n   ∂x  ∂x x=0 x=L

∂2,p 

=0

p +Ls +Ln

Fig. 1. Schematic of a lithium ion battery.

The potential in the solution phase and its flux are continuous at the interfaces of the electrodes and the separator.

L. Cai, R.E. White / Journal of Power Sources 196 (2011) 5985–5989

5987

In the above equations the pore wall flux, Ji , is determined by the Bulter–Volmer equation: 0.5 ci0.5 Ji = ki (cs,i,max − cs,i,surf )0.5 cs,i,surf



× exp

 0.5F  RT

i

 0.5F 

− exp −

RT

i

(7)

and the overpotential of the intercalation reaction is given by: i = 1,i − 2,i − Ui The energy balance is given by [4,5]: Cp

dT ∂2 T =  2 + Qrxn + Qrev + Qohm dt ∂x

(8)

with the boundary conditions determined by Newton’s cooling law:



∂T  −  ∂x

x=0 ∂T  −  ∂x

= h (T∞ − T ) = h (T − T∞ )

x=Lp +Ls +Ln

where h is the heat transfer coefficient, T∞ is the environmental temperature, Qrxn is the total reaction heat generation rate, Qrev is the total reversible heat generation rate, Qohm is the total ohmic heat generation rate. The heat fluxes are defined by: Qrxn = FaJ(1 − 2 − U) Qrev = FaJT

∂U ∂T

 Qohm = eff

∂1 ∂x

2

 + eff

∂2 ∂x

2 +

2eff RT 0 ∂(ln c) ∂2 ) (1 − t+ F ∂x ∂x

The temperature dependent open circuit potential of electrode i is approximated by Taylor’s first order expansion around a reference temperature:

 dU 

Ui = Ui,ref + (T − Tref )

dT

i

where Ui,ref is the open circuit potential under the reference temperature Tref .

Fig. 2. Geometries and variables coupling between the geometries using COSMOL Multiphysics 3.5a. (The top 1D geometry consists of three segments which denote the positive electrode, the separator and the negative electrode. The bottom two rectangles represent the solid phases in the positive electrode and the negative electrode, respectively. The pore wall flux is extracted from the 1D geometry and projected to the top boundary of the 2D geometry. The concentration of Li ions on the top boundary of the 2D geometry is projected to the 1D domain as the surface concentration of the solid particles.)

The thermal behavior of the Li ion battery during pulse discharge is also simulated in COSMOL Multiphysics. The battery is discharged for 3000 s at C/2 rate first and then discharged at 3 C rate until the cell voltage drops to 2.5 V. The change of the applied current density is implemented by using the smoothed Heaviside function “flsmhs” and is shown in Fig. 3. The model parameters can be found in Ref. [5].Simulation results Fig. 4 shows the temperature on the cell surface at 1 C discharge process under three different cooling conditions where the heat transfer coefficient is 10.0, 1.0 and 0.1 W m−2 K−1 , respectively, and two limiting conditions: the isothermal condition and the adiabatic condition. The thermal effect on the cell voltage shown in Fig. 5. The cell provides more discharge capacity when it is placed in a better heat isolation environment (i.e. adiabatic condition). In a better isolated environment, the cell temperautre increases faster during the 1 C discharge process which results in the higher diffusion coefficient for the binary electrolyte and reduces the diffusion limitations. The reduction of the diffusion limit in the binary electrolyte can be verified by comparing the concentration profile of the electrolyte under different cooling conditions. Fig. 6 shows the concentration profiles of the electrolyte at the end of 1 C discharge process under the two limiting conditions. The concentration profile under the adiabatic condition is flatter than that in the isothermal case, which indicates a better diffusion property in the electrolyte under the adiabatic condition than under the isothermal condition.

3. The COMSOL Multiphysics (MP) 3.5a model The mathematical model described in Section 2 is a multi-scale model. We developed several geometries using this software: a 1D geometry which consists of three sequentially connected lines to represent the positive electrode, the separator and the negative electrode, respectively, a 2D geometry which consists of two rectangles to denote the solid phase in the electrodes. The geometries are shown in Fig. 2. The vertical coordinate in the 2D geometry indicates the radial direction of the solid particles. Since we ignore the diffusion of Li ions in the x-direction in the particle, the corresponding diffusion coefficient is set to zero in this direction. The concentration of the binary electrolyte, the potential in the electrolyte, the potential in the solid phase and the pore wall flux are solved in the 1D geometry. The concentration of Li ions in the solid phase is solved in the 2D geometry. The pore wall flux is extruded from the 1D domain and projected to the top boundary of the 2D geometry by using “subdomain extrusion coupling variables” in COMSOL Multiphysics. The concentration of Li ions on the top boundary in the 2D geometry is projected to the 1D domain by using “boundary extrusion coupling variables”.

Fig. 3. Current density profile in the discharge process including a 3 C pulse.

5988

L. Cai, R.E. White / Journal of Power Sources 196 (2011) 5985–5989

Fig. 4. Temperature on the cell surface during 1 C discharge process under different cooling conditions.

Fig. 5. Cell voltage for 1 C discharge process under different cooling conditions.

Fig. 7 shows the cell temperature during the 1 C discharge process at different current rates as the heat transfer coefficient is 1.0 W m−2 K−1 . As expected, the cell gets hotter as the dishcarge current rate increases. It is also noticed that the wave part which appears in beginning of the temperature curve at low current rate (less than 2 C) does not exist in the high current rate cases. The wave part on the temperature curve is characterised by the reversible heat generation during discharging. Under low current

rate discharging, the reversible heat is roughly equivelant to the ohmic heat, but becomes unimportant as the discharge current rate increases. The P2D model mentioned in Section 2 is also useful for simulating the discharge process with pulse. Fig. 8 shows the cell voltage during the C/2 discharge for 3000 s followed by a 3 C pulse discharge until the cell voltage drops to 2.5 V. The corresponding temperature

Fig. 6. Comparison of the concentration profiles of the binary electrolyte at the end of the 1 C discharge process under the isothermal condition and the adiabatic condition.

Fig. 7. Temperature on the cell surface during discharge process under different current rates while the heat transfer coefficient h = 1.0 W m−2 K−1 where the DOD is defined as: DOD = time * C rate/3600.

L. Cai, R.E. White / Journal of Power Sources 196 (2011) 5985–5989

Fig. 8. Cell voltage at C/2 discharge for 3000 s followed by a 3 C pulse discharge.

5989

Fig. 10. Concentration of the binary electrolyte at the two ends of the cell (the bottom line is at x = 0, the top line is at x = 1.65 × 10−4 m) in the C/2 discharge process with a 3 C pulse. Mathematical modeling of a lithium ion battery with thermal effects in COMSOL.

the beginning of the pulse, the concentration of the electrolyte changes extremely, after that it relaxes and tend to a stable value. 5. Conclusions The thermal behavior of a lithium ion battery during galvanostatic discharge w/o pulse can be predicted by using COMSOL Multiphysics (MP) 3.5a. Though better thermal isolation environment improves the discharge capacity, the higher cell temperature raises the risk of thermal runaway and more rapid degradation of the cell. References

Fig. 9. Temperature on the cell surface in the discharge process with 3 C pulse.

on the surface of the cell is also plotted in Fig. 9. The surface temperature at the end of the 3 C pulse is slightly less than that in the pure 3 C discharge process. Fig. 10 shows the concentration of the binary electrolyte at the two ends of the cell during the pulse discharge process. At

[1] M. Doyle, T.F. Fuller, J. Newman, Journal of Electrochemical Society 140 (1993) 1526–1533. [2] T.F. Fuller, M. Doyle, J. Newman, Journal of Electrochemical Society 141 (1994) 1–10. [3] N. Lars Ole Valoen, Journal of Electrochemical Society 152 (2005) A882– A891. [4] W.B. Gu, C.Y. Wang, Journal of Electrochemical Society 147 (2000) 2910– 2922. [5] K. Kumaresan, G. Sikha, R.E. White, Journal of Electrochemical Society 155 (2008) A164–A171.