COMSOL 6.3 - Modeling Impedance in the Lithium-Ion Battery

Modeling Impedance in the Lithium-Ion Battery

Electrochemical impedance spectroscopy (EIS) is a valuable method to investigate electrochemical systems such as batteries. For a battery, where potential perturbations of varying frequency are applied on an electrode, the impedance response gives insight into several battery properties and processes. At high frequencies, short time-scale processes such as capacitance, electrochemical reactions, and local resistances affect the impedance. On the other hand, at low frequencies, diffusion in the electrolyte and active material particles (that is, large time-scale processes) contribute to the impedance.

This example simulates the impedance of a full lithium-ion battery cell using the Lithium-Ion Battery interface with an AC Impedance Stationary study. The model also reproduces the results by Abraham and others (Ref. 1) for sinusoidal potential perturbations between 10 mHz to 1 kHz after model fitting using the Parameter Estimation study step.

This tutorial model requires the Optimization Module.

Model Definition

The model consists of the following components:

•

Negative porous electrode: LTO (Li4Ti5O12) active material, electronic conductor, and filler.

•

Separator: Celgard 2325

•

Positive porous electrode: NCA (LiNi0.08Co0.15Al0.05O2) active material, electronic conductor, and filler.

•

Electrolyte: 1.2 M LiPF6 in EC:EMC (3:7 by weight)

A 1D model with three domains of different thickness is used: a 115 μm negative electrode, a 50 μm separator, and a 35 μm positive electrode. The model is set up in a similar fashion as the 1D Isothermal Lithium-Ion Battery, but with double-layer capacitances on the active electrode materials and on the electronic conductor in the positive electrode. The model approach using an additional double-layer at the electronic conductor has, for example, been used by Brown and others (Ref. 2) to describe experimental EIS spectra measured between 1 mHz and 1 kHz on NCA electrodes.

An AC Stationary Impedance study is used to solve the problem, implying that all variables are shifted from being time dependent to frequency dependent, as illustrated with the following expression:

where n is the variable and subscript 0 denotes the initial value around which the perturbation takes place and tilde the complex perturbation. Furthermore, i is the imaginary unit, f the frequency, and t the time.

The boundary of the positive electrode current-collector is set to a sinusoidal perturbation with a 10 mV amplitude and the boundary of the negative electrode current-collector fixed to 0 V, that is, ground.

The cell impedance, Z (SI unit: Ω·m2), is calculated at the boundary of the positive electrode current-collector as follows:

where

(SI unit: V) is the potential, n the boundary normal and Is (SI unit: Am−2) the current density in the solid. (Note that the lindev operator is used in the model file in order to use the perturbing values of the variables in the evaluation.)

Subtracting the reference electrode potential from the solid potential in the numerator enables the analysis of the impedance of the NCA positive electrode versus the reference electrode, thus separating the contribution of each electrode to the impedance. The LTO impedance versus the reference is consequently given by the difference between the cell impedance and NCA impedance versus the reference.

The model is computed for frequencies between 10 mHz and 1 kHz. This range is common in EIS to minimize any errors from experimental interference (for example, Ref. 2).

Optimization

The model is initially solved for the presented experimental set-up (Ref. 1) with the aid of parameters given in Ref. 2. However, since the simulated impedance deviates to a large extent from the experiments, an optimization is performed using the Parameter Estimation study step. Global Least-Squares Objective is set to minimize the least-squares error between the simulated and measured NCA impedance versus the reference electrode through variation of some control parameters. These parameters are given in Table 1.

Table 1: control parameters with bounds For a Cell voltage of 2.30V.
Control parameters	Initial value	Bounds
Exchange current density NCA	1 Am-2	Lower: 1 Am-2 Upper: 6 Am-2
Double-layer capacitance NCA	5·10-1 Fm-2	Lower: 1·10-1 Fm-2 Upper: 1 Fm-2
Film resistance NCA	1·10-4 Ωm2	Lower: 1·10-6 Ωm2 Upper: 5·10-3 Ωm2
Volumetric double-layer capacitance positive electronic conductor	1·105 Fm-3	Lower: 1·105 Fm-3 Upper: 3·105 Fm-3

The first three parameters have bounds selected as anticipated by the information given by Abraham and others (Ref. 1). The last parameter depend on the size of the electronic conductor particles. It is also common to run parametric sweeps before selecting control parameters and to narrow down the bounds. Normally the optimization time is reduced and the optimization accuracy is improved when few control parameters and narrow bound intervals are selected.

Study settings

The problem is solved with an AC Impedance, Initial Values study. The optimization solver is set up by adding a Parameter Estimation study step to the study and selecting the IPOPT method with the numeric gradient method.

Results and Discussion

Simulated and experimental Nyquist spectra from the LTO/NCA cell at 2.30 V are displayed in Figure 1. The simulation utilizes the parameters from Ref. 2. Characteristic semi-circles and low-frequency tails are seen. The largest contribution to the cell impedance is in the positive NCA electrode. The difference between the experimental and simulated impedance spectra of the NCA electrode versus reference is evident: The simulation shows one large semi-circle and an indication of another semi-circle at high frequencies. The experimental spectrum has a smaller depressed (elongated) semi-circle.

Figure 1: Simulated and experimental Nyquist spectra for 10 mHz to 1 kHz. Parameters taken from Ref. 2.

Since the difference is most prominent in the semi-circle frequency region (that is, mid-high frequencies) mainly the parameters describing short time-scale processes do not fit the experiments. Four of these parameters are investigated with a parametric sweep to determine their impact on the impedance. The sweep is shown in Figure 2.

Figure 2: Nyquist plots when testing the impact of five parameters for 10 mHz to 1 kHz.

As a comparison, the sweep also includes variation of the radius of the positive electrode material particles, mainly affecting the large time-scale process of diffusion in the particles. Its impact is therefore more evident at lower frequencies. This is seen in Figure 2, where both the angle and shape of the low-frequency tail change. The results also show that the parameters describing large time-scale processes seem to be quite accurate in the model.

These findings are consistent with the observations by Abraham and others (Ref. 1) indicating that the electrochemical reaction is less limited in their electrodes. With this input Table 1 is set up and the optimization is run. The resulting spectra are shown in Figure 3.

Figure 3: Nyquist spectra with optimized parameters for 10 mHz to 1 kHz. Squares denote the data points used in the optimization.

The simulated and experimental spectra for NCA versus reference are now seen to overlap well. The optimized control parameters are tabulated in Table 2.

Table 2: Change of control parameters After optimization.
Control Parameters	Optimized Value
Exchange current density NCA	2.6 Am-2
Double-layer capacitance NCA	6.0·10-1 Fm-2
Film resistance NCA	1.0·10-6 Ωm2
Volumetric double-layer capacitance positive electronic conductor	1.9·105 Fm-3

References

1. D.P. Abraham, S. Kawauchi, and D.W. Dees, “Modeling the impedance versus voltage characteristics of LiNi0.8Co0.15Al0.05O2”, Electrochim. Acta, vol. 53, pp. 2121–2129, 2008.

2. S. Brown, N. Mellgren, M. Vynnycky, and G. Lindbergh, “Impedance as a Tool for Investigating Aging in Lithium-Ion Porous Electrodes. II. Positive Electrode Examination,” J. Electrochem. Soc, vol. 155, p. A320, 2008.

Battery_Design_Module/Lithium-Ion_Batteries,_Performance/li_battery_impedance

Modeling Instructions

From the menu, choose .

New

In the window, click

Model Wizard

In the window, click

In the tree, select > > .

Click .

Click

In the tree, select > .

Click

Global Definitions

Load parameters from a text file.

Parameters 1

In the window, under click .

In the window for , locate the section.

Click

Browse to the model’s Application Libraries folder and double-click the file li_battery_impedance_parameters.txt.

Definitions

The impedance is calculated using variables defined in a text file. It is defined at the rightmost boundary versus ground (at the leftmost boundary) or a reference electrode. The lindev operator explicitly evaluates the deviation from the stationary solution in the perturbed solution.

Variables 1

In the window, under right-click and choose .

In the window for , locate the section.

Click

Browse to the model’s Application Libraries folder and double-click the file li_battery_impedance_variables.txt.

Geometry 1

The model geometry includes two porous electrodes with a separator in between. A reference electrode is positioned in the middle of the separator.

Interval - Negative Electrode

In the window, under right-click and choose .

In the window for , type Interval - Negative Electrode in the text field.

Locate the section. In the table, enter the following settings:


Coordinates (m)
0
L_neg

Locate the section. Select the checkbox.

Interval - Separator

In the window, right-click and choose .

In the window for , type Interval - Separator in the text field.

Locate the section. From the list, choose .

In the text field, type L_neg.

In the table, enter the following settings:


Lengths (m)
L_sep

Locate the section. Select the checkbox.

Interval - Positive Electrode

Right-click and choose .

In the window for , type Interval - Positive Electrode in the text field.

Locate the section. From the list, choose .

In the text field, type L_neg+L_sep.

In the table, enter the following settings:


Lengths (m)
L_pos

Locate the section. Select the checkbox.

Point - Reference Electrode

Right-click and choose .

In the window for , type Point - Reference Electrode in the text field.

Locate the section. In the text field, type L_neg+L_sep/2.

Locate the section. Select the checkbox.

Click

Click the

button in the toolbar.

Add Material

The cell materials are available in the Battery Material Library.

In the toolbar, click

to open the window.

Go to the window.

In the tree, select > > .

Click the button in the window toolbar.

In the tree, select > > .

Click the button in the window toolbar.

In the tree, select > > .

Click the button in the window toolbar.

In the toolbar, click

to close the window.

Materials

LTO, Li4Ti5O12 (Negative, Li-ion Battery) (mat2)

In the window for , locate the section.

From the list, choose .

NCA, LiNi0.8Co0.15Al0.05O2 (Positive, Li-ion Battery) (mat3)

In the window, click .

In the window for , locate the section.

From the list, choose .

Definitions

Integration 1 (intop1)

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

Select Boundary 5 only.

Lithium-Ion Battery (liion)

Double layers exist on both active materials.

Separator 1

In the window, under > click .

In the window for , locate the section.

In the εl text field, type epsl_sep.

Locate the section. From the list, choose . In the fl text field, type epsl_sep^brugl_sep.

From the list, choose . In the fDl text field, type epsl_sep^brugl_sep.

Porous Electrode 1

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

Locate the section. From the list, choose .

Locate the section. In the σs text field, type sigmas_neg.

Click to expand the section. From the list, choose .

In the Rfilm text field, type Rfilm_neg.

Locate the section. In the εs text field, type epss_neg.

In the εl text field, type epsl_neg.

Locate the section. From the list, choose . In the fl text field, type epsl_neg^brugl_neg.

From the list, choose . In the fDl text field, type epsl_neg^brugl_neg.

Particle Intercalation 1

In the window, click .

In the window for , locate the section.

From the list, choose .

Locate the section. In the cs,init text field, type mat2.def.csmax*soc0_neg.

Locate the section. In the rp text field, type rp_neg.

Porous Electrode Reaction 1

In the window, click .

In the window for , locate the section.

From the list, choose .

Locate the section. In the i0,ref(T) text field, type i0ref_neg.

Porous Electrode 1

In the window, click .

Porous Matrix Double Layer Capacitance 1

In the toolbar, click

and choose .

In the window for , locate the section.

In the Cdl text field, type cdl_neg.

Porous Electrode 2

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

Locate the section. From the list, choose .

Locate the section. In the σs text field, type sigmas_pos.

Locate the section. From the list, choose .

In the Rfilm text field, type Rfilm_pos.

Locate the section. In the εs text field, type epss_pos.

In the εl text field, type epsl_pos.

Locate the section. From the list, choose . In the fl text field, type epsl_pos^brugl_pos.

From the list, choose . In the fDl text field, type epsl_pos^brugl_pos.

Particle Intercalation 1

In the window, click .

In the window for , locate the section.

From the list, choose .

Locate the section. In the cs,init text field, type mat3.def.csmax*soc0_pos.

Locate the section. In the rp text field, type rp_pos.

Porous Electrode Reaction 1

In the window, click .

In the window for , locate the section.

From the list, choose .

Locate the section. From the list, choose .

In the i0 text field, type i0_pos.

Locate the section. From the list, choose . In the av text field, type as_pos.

Porous Electrode 2

In the window, click .

Porous Matrix Double Layer Capacitance 1

In the toolbar, click

and choose .

In the window for , locate the section.

In the Cdl text field, type cdl_pos.

From the list, choose . In the av,dl text field, type as_pos.

Initial Values 1

In the window, under > click .

In the window for , locate the section.

In the phil text field, type -mat2.def.Eeq(soc0_neg).

In the cl text field, type cl_init.

Initial Values 2

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

Locate the section. In the phil text field, type -mat2.def.Eeq(soc0_neg).

In the cl text field, type cl_init.

In the phis text field, type mat3.def.Eeq(soc0_pos)-mat2.def.Eeq(soc0_neg).

Electric Ground 1

Select ground on the current collector of the negative LTO electrode. Set a 10 mV perturbation around the cell voltage of 2.30 V.

In the toolbar, click

and choose .

Select Boundary 1 only.

Electric Potential 1

In the toolbar, click

and choose .

Select Boundary 5 only.

In the window for , locate the section.

In the ϕs,bnd text field, type E_cell_init.

Locate the section. Select the checkbox.

In the Rc text field, type R_curr.

Harmonic Perturbation 1

In the toolbar, click

and choose .

In the window for , locate the section.

In the Δϕs text field, type E_pert.

Additional Porous Electrode Material 1

Add an Additional Porous Electrode Material feature to model the electronic conductor in NCA.

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

Locate the section. From the list, choose .

Porous Electrode Reaction 1

No electrochemical reaction takes place on the electronic conductor, but a double layer capacitance exists.

In the window, right-click and choose .

Additional Porous Electrode Material 1

In the window, click .

Porous Matrix Double Layer Capacitance 1

In the toolbar, click

and choose .

In the window for , locate the section.

In the Cdl text field, type 1[F/m^2].

From the list, choose . In the av,dl text field, type cdlvol_cs_pos/1[F/m^2].

Reference Electrode 1

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

Mesh 1

In the window, under click .

In the window for , locate the section.

From the list, choose .

Size

In the window, under > click .

In the window for , locate the section.

Click the button.

Locate the section. In the text field, type 1e-5.

Size 1

In the window, click .

Select Domains 2 and 3 only.

In the window for , locate the section.

Click the button.

Locate the section.

Select the checkbox. In the associated text field, type 1e-6.

Size 2

In the window, click .

Select Boundaries 2–4 only.

In the window for , locate the section.

Click the button.

Locate the section.

Select the checkbox. In the associated text field, type 5E-7.

Study 1

Solve the model for frequencies between 10 mHz and 1 kHz.

Step 1: Frequency-Domain Perturbation

In the window, under click .

In the window for , locate the section.

In the text field, type 10^{range(-2,0.2,3)}.

In the toolbar, click

Solution 1 (sol1)

Save the solution of the non-optimized impedance results.

In the window, expand the > node.

Right-click and choose > .

Not Optimized

In the window, under > click .

In the window for , type Not Optimized in the text field.

Results

Follow these steps to generate Figure 1:

Start by loading experimental data of real and imaginary impedance into a table for plotting purposes. Organize the tables by naming them.

Experimental Impedance Table

In the toolbar, click

In the window for , type Experimental Impedance Table in the text field.

Locate the section. Click

Browse to the model’s Application Libraries folder and double-click the file li_battery_impedance_Zexp.csv.

Nyquist Plots

In the window, under click .

In the window for , type Nyquist Plots in the text field.

Click to expand the section. Locate the section. From the list, choose .

Locate the section. From the list, choose .

Locate the section. Select the checkbox.

Cell Impedance vs. Ground

In the window, expand the node, then click .

In the window for , type Cell Impedance vs. Ground in the text field.

Click to expand the section. Select the checkbox.

Find the subsection. Select the checkbox.

Clear the checkbox.

Nyquist Plots

Right-click > > and choose .

Impedance NCA vs. Reference

In the window for , type Impedance NCA vs. Reference in the text field.

Select Boundary 5 only.

Locate the section. In the text field, type -imag(Z_ref_NCA).

Click to expand the section. From the list, choose .

Locate the section. From the list, choose .

In the text field, type real(Z_ref_NCA).

Click to expand the section. Select the checkbox.

Find the subsection. Select the checkbox.

Clear the checkbox.

Impedance LTO vs. Reference

Right-click and choose .

In the window for , type Impedance LTO vs. Reference in the text field.

Locate the section. In the text field, type -imag(Z_ref_LTO).

Locate the section. In the text field, type real(Z_ref_LTO).

Experimental Impedance NCA vs. Reference

In the window, right-click and choose .

The Table Graph plots the experimental data that was imported into the Experimental Impedance Table.

In the window for , type Experimental Impedance NCA vs. Reference in the text field.

Locate the section. From the list, choose .

From the list, choose .

In the list, select .

Click to expand the section. Select the checkbox.

Find the subsection. Select the checkbox.

Clear the checkbox.

Nyquist Plots

In the window, click .

In the window for , locate the section.

From the list, choose .

Locate the section. Select the checkbox.

In the toolbar, click

Click the

button in the toolbar.

The Bode Plots are redundant in this study.

Impedance with Respect to Ground, Imaginary Part (liion), Impedance with Respect to Ground, Real Part (liion)

In the window, under , Ctrl-click to select and .

Right-click and choose .

Study 1

Add a parametric sweep to investigate the impact four key parameters have on the impedance of NCA versus the reference and how it relates to the experimental measurements.

Parametric Sweep

In the toolbar, click

In the window for , locate the section.

Click

Browse to the model’s Application Libraries folder and double-click the file li_battery_impedance_parametric_sweep.txt.

In the window, click .

In the window for , locate the section.

Clear the checkbox.

In the toolbar, click

Rename the solution for the parametric sweep.

Parametric Sweep

In the window, under > click .

In the window for , type Parametric Sweep in the text field.

Follow these steps to generate Figure 2:

Results

Nyquist Plots, Parametric Sweep

In the window, right-click and choose .

In the window for , type Nyquist Plots, Parametric Sweep in the text field.

Cell Impedance vs. Ground

In the window, expand the node.

Right-click and choose .

Impedance LTO vs. Reference

In the window, right-click and choose .

Impedance NCA vs. Reference

In the window, under > click .

In the window for , locate the section.

From the list, choose .

Locate the section. From the list, choose .

In the table, enter the following settings:


Legends
Impedance NCA vs. Reference, Original Settings
i0_pos=5 Am<sup>-2</sup>
cdl_pos=0.50 Fm<sup>-2</sup>
Rfilm_pos=1\<sup\>.\</sup\>10\<sup\>-5\</sup\> Ωm\<sup\>2\</sup\>
cdlvol_cs_pos=1<sup>.</sup>10<sup>4</sup> Fm<sup>-3</sup>
rp_pos=4<sup>.</sup>10<sup>-7</sup> m

Nyquist Plots, Parametric Sweep

In the window, click .

In the toolbar, click

Click the

button in the toolbar.

Add Study

In the toolbar, click

to open the window.

Go to the window.

Find the subsection. In the tree, select > .

Click the button in the window toolbar.

In the toolbar, click

to close the window.

Study 2

In the window for , locate the section.

Clear the checkbox.

Parameter Estimation

Select i0_pos, Rfilm_pos, cdl_pos, and cslvol_cs_pos as the four parameters to be optimized.

In the toolbar, click

and choose .

The experimental data is already available as a result table.

In the window for , locate the section.

From the list, choose .

Locate the section. In the table, click to select the cell at row number 1 and column number 2.

From the drop-down list, choose .

The second column of the experimental data contains the real-impedance global value.

In the table, click to select the cell at row number 2 and column number 3.

In the text field, type comp1.intop1(real(Z_ref_NCA)).

In the text field, type Real_impedance.

From the list, choose .

In the text field, type 1e-3.

The third column of the experimental data contains the imaginary-impedance global value.

In the table, click to select the cell at row number 3 and column number 3.

In the text field, type -comp1.intop1(imag(Z_ref_NCA)).

In the text field, type Imaginary_impedance.

From the list, choose .

In the text field, type 1e-3.

Locate the section. Click

four times.

Row by row, select the parameter name in the first column, then set the corresponding initial value, scale, and bounds as follows:


Parameter name	Initial value	Scale	Lower bound	Upper bound
i0_pos (Exchange current density positive electrode)	1[A/m^2]	1	1	5
Rfilm_pos (Film resistance positive electrode)	1e-4[m^2/S]	1e-3	1e-6	5e-3
cdl_pos (Double layer capacitance positive electrode material)	0.5[F/m^2]	1	0.10	1
cdlvol_cs_pos (Volumetric capacitance of electronic conductor in positive electrode)	1e5[F/m^3]	1e6	1e5	3e5