COMSOL 6.4 - Lithium-Sulfur Battery

Lithium-sulfur (Li-S) batteries are used in niche applications with high demands for specific energy densities, which may be as high as 500–600 Wh/kg.

The chemistry is complex, with multiple polysulfide species participating in several charge transfer reactions. The chemistry also involves precipitation and dissolution of multiple solid species.

This example models discharge of a Li-S cell at two different discharge rates. The electrolyte charge and mass transport of a lithium salt and 6 polysulfides is included, as well as the precipitation-dissolution of solid octasulfur (S8) and lithium sulfide (Li2S) in the separator and positive electrode.

The model is based on a paper by Zhang and others (Ref. 1).

Model Definition

The 1D model geometry is shown in Figure 1. The separator and the positive (cathode during discharge) porous electrode are defined as domains in the geometry. The leftmost boundary represents the negative lithium metal electrode, and the rightmost boundary represents the positive metal current collector.

Figure 1: Model geometry.

The electrolyte mass and charge transport are modeled using the Nernst-Planck equations with electroneutrality.

The following five charge transfer reactions and six participating polysulfide species are considered in the positive electrode domain:

(1)

(2)

(3)

(4)

(5)

The charge transfer reactions are defined using the Nernst equation in combination with the Butler–Volmer equation, assuming the law of mass action.

In addition to the above six polysulfide species, also the species Li+ and A- are included in the electrolyte, where A- is the counter anion of the lithium salt. The same species are also present in the separator domain.

Two precipitation-dissolution (nonfaradaic) reactions involving the solid species S8 and Li2S are included in the model (in both the domains) according to

(6)

(7)

The model is solved for two different discharge rates: 0.2C and 1C.

Results and Discussion

Figure 1 shows the cell voltages during a 0.2C and a 1C discharge. The voltage and capacity gets significantly lowered for the higher discharge rate.

Figure 2: Voltage discharge curves at 0.2C and 1C.

Figure 3 and Figure 4 show the concentration profiles in the cell at end of discharge at 0.2C and 1C, respectively. The amount of remaining polysulfide species correspond to the lowered capacity at 1C compared 0.2C that was seen in Figure 2.

Figure 3: Concentration of electrolyte species at end of discharge at 0.2C.

Figure 4: Concentration of electrolyte species at end of discharge at 1C.

Finally, Figure 5 and Figure 6 show the volume fraction profiles of Li2S(s) during discharge at 0.2C and 1C, respectively. The 1C profiles are generally less uniform than the profiles at 0.2C, indicating a less uniform current distribution for the higher discharge rate.

Figure 5: Volume fraction profiles of Li2S(s) during discharge at 0.2C.

Figure 6: Volume fraction profiles of Li2S(s) during discharge at 1C.

Reference

1. T. Zhang, M. Marinescu, S. Walus, and G. Offer, “Modeling transport-limited discharge capacity of lithium-sulfur cells,” Electrochimica Acta, vol. 219, pp. 502–508, 2016.

Battery_Design_Module/Batteries,_General/lithium_sulfur

Modeling Instructions

From the menu, choose .

New

In the window, click

Model Wizard

In the window, click

In the tree, select > > .

Click .

In the text field, type 8.

Type in the names of the electrolyte concentration variables in the table as follows. (The solid S8(s) and Li2S(s) species will be added later.)

In the table, enter the following settings:


S8
S8_2m
S6_2m
S4_2m
S2_2m
S_2m
Li_1p
A_1m

Click

In the tree, select > .

Click

Tertiary Current Distribution, Nernst–Planck (tcd)

The Nernst-Planck model assumes electroneutrality, and eliminates one of the concentration-dependent variables based on the electroneutrality condition. For this model it is suitable to select the A_1m (the electrolyte salt anion) as to be taken from electroneutrality since its concentration is relatively high in relation to the other species, and since it does not participate in any electrode reactions.

In the window, under click .

In the window for , locate the section.

From the list, choose .

Global Definitions

Parameters 1

Load the model parameters from a text file.

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 lithium_sulfur_parameters.txt.

Geometry 1

Interval 1 (i1)

The geometry consists of two domains: the positive porous electrode and the separator. The negative electrode is modeled as a boundary condition and is not added as a domain at this point.

In the window, under right-click and choose .

In the window for , locate the section.

From the list, choose .

In the table, enter the following settings:


Lengths (m)
L_sep
L_pos

Click

Tertiary Current Distribution, Nernst–Planck (tcd)

In the window, under click .

In the window for , locate the section.

In the Ac text field, type A_cell.

Start setting up the physics for charge and mass transport in the separator.

Species Charges 1

In the window, under > click .

In the window for , locate the section.

In the zS82m text field, type -2.

In the zS62m text field, type -2.

In the zS42m text field, type -2.

In the zS22m text field, type -2.

In the zS2m text field, type -2.

In the zLi1p text field, type 1.

In the zA1m text field, type -1.

Separator 1

In the toolbar, click

and choose .

Select Domain 1 only.

In the window for , locate the section.

In the DS8 text field, type D_S8.

In the DS82m text field, type D_S8_2m.

In the DS62m text field, type D_S6_2m.

In the DS42m text field, type D_S4_2m.

In the DS22m text field, type D_S2_2m.

In the DS2m text field, type D_S_2m.

In the DLi1p text field, type D_Li_1p.

In the DA1m text field, type D_A_1m.

Locate the section. In the εl text field, type epsl_sep_0.

Note that, by default, the diffusion coefficients will be corrected by a Bruggeman relation to account for the effect of the pore network.

Porous Electrode 1

In the toolbar, click

and choose .

Select Domain 2 only.

In the window for , locate the section.

In the DS8 text field, type D_S8.

In the DS82m text field, type D_S8_2m.

In the DS62m text field, type D_S6_2m.

In the DS42m text field, type D_S4_2m.

In the DS22m text field, type D_S2_2m.

In the DS2m text field, type D_S_2m.

In the DLi1p text field, type D_Li_1p.

In the DA1m text field, type D_A_1m.

Locate the section. From the σs list, choose . In the associated text field, type sigma_s.

Locate the section. In the εs text field, type 1-epsl_pos_0.

In the εl text field, type epsl_pos_0.

For the electrode phase assume that the conductivity value entered is the effective conductivity and hence already accounts for the porous structure:

Locate the section. From the list, choose .

Porous Electrode Reaction 1

Now start setting up the electrode kinetics. In total, there are five active electrode reactions in the porous electrode.

In the window, click .

In the window for , locate the section.

In the νS8 text field, type -1/2.

In the νS82m text field, type 1/2.

Locate the section. In the Eeq,ref(T) text field, type Eeq_1_ref.

Click to expand the section. In the table, enter the following settings:


Electrolyte species	Reference concentrations (mol/m^3)
S8	c_S8_ref
S8_2m	c_S8_2m_ref

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

Locate the section. In the av text field, type Av_pos.

Porous Electrode 1

In the window, click .

Porous Electrode Reaction 2

In the toolbar, click

and choose .

In the window for , locate the section.

In the νS82m text field, type -3/2.

In the νS62m text field, type 2.

Locate the section. In the Eeq,ref(T) text field, type Eeq_2_ref.

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


Electrolyte species	Reference concentrations (mol/m^3)
S8_2m	c_S8_2m_ref
S6_2m	c_S6_2m_ref

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

Locate the section. In the av text field, type Av_pos.

Porous Electrode 1

In the window, click .

Porous Electrode Reaction 3

In the toolbar, click

and choose .

In the window for , locate the section.

In the νS62m text field, type -1.

In the νS42m text field, type 3/2.

Locate the section. In the Eeq,ref(T) text field, type Eeq_3_ref.

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


Electrolyte species	Reference concentrations (mol/m^3)
S6_2m	c_S6_2m_ref
S4_2m	c_S4_2m_ref

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

Locate the section. In the av text field, type Av_pos.

Porous Electrode 1

In the window, click .

Porous Electrode Reaction 4

In the toolbar, click

and choose .

In the window for , locate the section.

In the νS42m text field, type -1/2.

In the νS22m text field, type 1.

Locate the section. In the Eeq,ref(T) text field, type Eeq_4_ref.

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


Electrolyte species	Reference concentrations (mol/m^3)
S4_2m	c_S4_2m_ref
S2_2m	c_S2_2m_ref

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

Locate the section. In the av text field, type Av_pos.

Porous Electrode 1

In the window, click .

Porous Electrode Reaction 5

In the toolbar, click

and choose .

In the window for , locate the section.

In the νS22m text field, type -1/2.

In the νS2m text field, type 1.

Locate the section. In the Eeq,ref(T) text field, type Eeq_5_ref.

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


Electrolyte species	Reference concentrations (mol/m^3)
S2_2m	c_S2_2m_ref
S_2m	c_S_2m_ref

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

Locate the section. In the av text field, type Av_pos.

Separator 1

There are also two solid species in the model. Define these in the separator as follows:

In the window, under > click .

In the window for , click to expand the section.

Click

In the table, enter the following settings:


Species	Density (kg/m^3)	Molar mass (kg/mol)
Li2S_s	rho_Li2S_s	M_Li2S_s

Add the settings for S8(s) in a second row in the table as follows:

Click

In the table, enter the following settings:


Species	Density (kg/m^3)	Molar mass (kg/mol)
S8_s	rho_S8_s	M_S8_s

Porous Electrode 1

The same solid species are also present in the electrode.

In the window, click .

In the window for , click to expand the section.

Click

In the table, enter the following settings:


Species	Density (kg/m^3)	Molar mass (kg/mol)
Li2S_s	rho_Li2S_s	M_Li2S_s

Also for the porous electrode, add the settings for S8(s) in a second row in the table:

Click

In the table, enter the following settings:


Species	Density (kg/m^3)	Molar mass (kg/mol)
S8_s	rho_S8_s	M_S8_s

Definitions

Now define the dissolution-precipitation rates of the solid species.

Variables - Separator

Import some variable expressions for the separator that define the volume fractions of the solid species. These expressions are based on the automatically defined concentration variables for each species set up in the section of the separator node.

In the window, under right-click and choose .

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

Locate the section. From the list, choose .

Select Domain 1 only.

Locate the section. Click

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

Variables - Positive Electrode

In the window, right-click and choose .

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

Locate the section. From the list, choose .

Select Domain 2 only.

Locate the section. Click

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

Note that an expression for the specific surface area, which depends on the local electrolyte volume fraction, is also present in the variables you just imported.

Variables - All Domains

Right-click and choose .

In the window for , type Variables - All Domains in the text field.

Locate the section. Click

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

The next step is to add these rate expressions as in each domain, together with initial values for the solid concentration variables.

Tertiary Current Distribution, Nernst–Planck (tcd)

Non-Faradaic Reactions - Li2S(s)

In the window, under > > click .

In the window for , type Non-Faradaic Reactions - Li2S(s) in the text field.

Use the R_Li2S_s rate variable to set the rate for S_2m, Li_1p, and Li2S_s as follows:

Locate the section. In the RS2m text field, type -R_Li2S_s/tcd.epsl.

In the RLi1p text field, type -2*R_Li2S_s/tcd.epsl.

In the table, enter the following settings:


Species	Reaction rate (mol/(m^3*s))
Li2S_s	R_Li2S_s

Separator 1

In the window, click .

Non-Faradaic Reactions - S8(s)

In the toolbar, click

and choose .

In the window for , type Non-Faradaic Reactions - S8(s) in the text field.

Here, use the R_S8_s rate variable to set the rate for S8 and S8_s:

Locate the section. In the RS8 text field, type -R_S8_s/tcd.epsl.

In the table, enter the following settings:


Species	Reaction rate (mol/(m^3*s))
S8_s	R_S8_s

Separator 1

In the window, click .

Initial Values for Dissolving–Depositing Species 1

In the toolbar, click

and choose .

In the window for , locate the section.

In the table, enter the following settings:


Species	Concentration (mol/m^3)
Li2S_s	eps_Li2S_s_sep_0/Vm_Li2S_s
S8_s	eps_S8_s_sep_0/Vm_S8_s

Porous Electrode 1

In the window, under > click .

Non-Faradaic Reactions - Li2S(s)

In the toolbar, click

and choose .

In the window for , type Non-Faradaic Reactions - Li2S(s) in the text field.

Locate the section. In the RS2m text field, type -R_Li2S_s/tcd.epsl.

In the RLi1p text field, type -2*R_Li2S_s/tcd.epsl.

In the table, enter the following settings:


Species	Reaction rate (mol/(m^3*s))
Li2S_s	R_Li2S_s

Porous Electrode 1

In the window, click .

Non-Faradaic Reactions - S8(s)

In the toolbar, click

and choose .

In the window for , type Non-Faradaic Reactions - S8(s) in the text field.

Locate the section. In the RS8 text field, type -R_S8_s/tcd.epsl.

In the table, enter the following settings:


Species	Reaction rate (mol/(m^3*s))
S8_s	R_S8_s

Porous Electrode 1

In the window, click .

Initial Values for Dissolving–Depositing Species 1

In the toolbar, click

and choose .

In the window for , locate the section.

In the table, enter the following settings:


Species	Concentration (mol/m^3)
Li2S_s	eps_Li2S_s_pos_0/Vm_Li2S_s
S8_s	eps_S8_s_pos_0/Vm_S8_s

Electrode Surface 1

Specify the negative electrode and electrode reaction (Li/Li+) as follows:

In the toolbar, click

and choose .

Select Boundary 1 only.

This boundary is to be grounded, so the default voltage setting of 0 V does not need to be changed.

Electrode Reaction 1

In the window, click .

In the window for , locate the section.

In the νLi1p text field, type -1.

Locate the section. In the Eeq,ref(T) text field, type Eeq_Li_ref.

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

Load Cycle 1

In the toolbar, click

and choose .

Select Boundary 3 only.

In the window for , locate the section.

From the list, choose .

Locate the section. From the list, choose .

In the Emin text field, type 1.7[V].

Locate the section. Select the checkbox.

The probe will plot the cell voltage while solving.

Definitions

Load Cycle Probe (tcd_lc1_volt)

In the window, under > click .

In the window for , locate the section.

Select the checkbox. In the associated text field, type Cell voltage.

Tertiary Current Distribution, Nernst–Planck (tcd)

Load Cycle 1

In the window, under > click .

Current 1

In the toolbar, click

and choose .

Select Boundary 3 only.

In the window for , locate the section.

In the Iset text field, type -I_1C*C.

You will vary the parameter C later when solving the model to simulate different discharge rates.

Initial Values 1

Now set the initial concentrations for the electrolyte species.

In the window, under > click .

In the window for , locate the section.

In the S8 text field, type c_S8_ref.

In the S82m text field, type c_S8_2m_ref.

In the S62m text field, type c_S6_2m_ref.

In the S42m text field, type c_S4_2m_ref.

In the S22m text field, type c_S2_2m_ref.

In the S2m text field, type c_S_2m_ref.

In the Li1p text field, type c_Li_1p_ref.

In the phis text field, type Eeq_1_ref.

Mesh 1

All the necessary physics settings are now completed. Due to steep gradients of the Li2S(s) species, this model needs a very well resolved mesh close to the separator-positive electrode boundary. Modify the default mesh as follows:

In the window, under right-click and choose .

Size 1

In the window, right-click and choose .

In the window for , locate the section.

From the list, choose .

Select Boundary 2 only.

Locate the section. Click the button.

Locate the section.

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

Edge 1

Right-click and choose .

Study 1

Next, set up the solver sequence. Add a parametric sweep to solve for a range of different C-rates.

Parametric Sweep

In the toolbar, click

In the window for , locate the section.

Click

In the table, enter the following settings:


Parameter name	Parameter value list	Parameter unit
C (C-rate)	0.2 1

The above settings mean that the study will solve for two different current-density boundary conditions as defined in the node you added before.

Step 2: Time Dependent

Now set the times of the time-dependent solver. Use the C-parameter to solve for a range of time steps corresponding to a range from 100% to (minimum) 0% nominal state of charge.

In the window, click .

In the window for , locate the section.

From the list, choose .

In the text field, type range(0,0.01/C,1/C).

The model is now ready for solving. It will take about a minute to solve.

In the toolbar, click