A zero dimensional model of lithium–sulfur batteries during charge and discharge

Abstract

Lithium–sulfur cells present an attractive alternative to Li-ion batteries due to their large energy density, safety, and possible low cost. Their successful commercialisation is dependent on improving their performance, but also on acquiring sufficient understanding of the underlying mechanisms to allow for the development of predictive models for operational cells. To address the latter, we present a zero dimensional model that predicts many of the features observed in the behaviour of a lithium–sulfur cell during charge and discharge. The model accounts for two electrochemical reactions via the Nernst formulation, power limitations through Butler–Volmer kinetics, and precipitation/dissolution of one species, including nucleation. It is shown that the flat shape of the low voltage plateau typical of the lithium–sulfur cell discharge is caused by precipitation. During charge, it is predicted that the dissolution can act as a bottleneck, because for large enough currents the amount that dissolves becomes limited. This results in reduced charge capacity and an earlier onset of the high plateau reaction, such that the two voltage plateaus merge. By including these effects, the model improves on the existing zero dimensional models, while requiring considerably fewer input parameters and computational resources than one dimensional models. The model also predicts that, due to precipitation, the customary way of experimentally obtaining the open circuit voltage from a low rate discharge might not be suitable for lithium–sulfur. This model can provide the basis for mechanistic studies, identification of dominant effects in a real cell, predictions of operational behaviour under realistic loads, and control algorithms for applications.

---

1 Introduction

Lithium sulfur batteries (Li–S) are seen as a candidate for the next generation of batteries, as they are potentially low cost and have a theoretical energy density of ∼2500 W h kg⁻¹ and a practical one estimated at ∼600 W h kg⁻¹.¹ These figures of merit represent a step change improvement compared to the energy density of current Li-ion technology of 250 W h kg⁻¹ (theoretical maximum 400 W h kg⁻¹)¹ and are possible due to the fact that the Li–S mechanism is different from that of Li-ion cells. Energy is stored and released through phase conversion rather than intercalation, and each sulfur atom can, in principle, contribute sixteen electrons to the circuit. The phase conversion mechanism, however, brings about many new challenges that must be addressed before Li–S can be considered competitive for large scale deployment. Currently, relatively high energy densities are achieved at cell level (300–400 W h kg⁻¹), but with the limited cycle life of a few hundred cycles, as commercialised by OXIS Energy and inferred from cathode development data.² Conversely, thousands of cycles are achieved for significantly lower energy densities, similar to those exhibited by lithium ion chemistries, as commercialised by OXIS Energy. Long cycle life can be achieved either by using excess electrolyte and Li in the cell,³ or by limiting the operating voltage window.⁴ The use of Li–S in applications is currently limited by low charge efficiency and relatively high self-discharge rates, but their inherent safety is an advantage in some niche applications.⁵

In order to improve the practical energy densities and cycle life of Li–S, most studies are currently aimed at decreasing the polysulfide shuttle and increasing sulfur utilisation, by improving material properties,^6,7 and by revealing the underlying mechanisms through characterisation experiments,⁸ amongst others. The necessary step of linking characterisation with design has not been explored sufficiently in the literature. Mechanistic models based on characterisation studies have been proposed, as discussed below, but none have been used yet to help design better cells or direct materials research. This state is similar to that of the early days of incumbent lithium ion chemistries.⁹

The first documented Li–S model was developed by Mikhaylik and Akridge¹⁰ and focused on studying the interplay between polysulfide shuttle and charging current. Their two-step reaction zero dimensional model includes heat generation as a result of the shuttle phenomenon, based on conclusions from their earlier study.¹¹ An expansion on the Nernst description of a zero dimensional Li–S cathode is explored by Moy et al.,¹² by adding intermediate reaction steps in the chain of polysulfide reduction, which they relate to the typical regions of the discharge curve. A one dimensional steady state model of polysulfide diffusion through the separator is added, with the aim of interpreting the results of their proposed method for shuttle rate measurement and of predicting capacity fade caused by shuttle. While providing a good fit to experimental data, their model is only applicable to situations in which the cell is maintained at constant voltage. Capacity fade predictions based on a similarly simple reaction chain were obtained by Risse et al.¹³ from a Markov chain model. None of these models account for diffusion limitations, precipitation/dissolution of insulating polysulfides and kinetic limitations, ignoring activation overpotentials. These phenomena can contribute to the cell voltage performance and capacity fade. Because of this and despite being useful for understanding some of the cell mechanisms, these models are of limited use for predictions of cell performance under operational conditions.

A more comprehensive one dimensional model proposed by Kumaresan et al.¹⁴ includes activation overpotentials in the form of Butler–Volmer kinetics, diffusion limitations of multicomponent transport in a dilute electrolyte, and precipitation of species via a rate constant, including nucleation. Polysulfide shuttle, however, is not included. The model allows for a detailed analysis of the interplay between the mechanisms in the system during cell discharge, providing direct access to the time evolution and spatial distribution of polysulfide species and to the various contributions to the cell potential. While it reproduces some of the main features of the discharge curve, such as the presence of two plateaus with a voltage dip in-between, a detailed sensitivity analysis by Ghaznavi and Chen^15–17 shows that their model as published is not suitable to predict behaviour during charge and capacity fade.

The framework of a similar model was developed by Neidhardt et al.,¹⁸ and used to reproduce the cell response to constant current charge, discharge and electrochemical impedance spectroscopy experiments.¹⁹ The model includes the presence of electrical double layers at the surface of either electrode. Polysulfide shuttle via an added reaction at the anode with resulting irreversible precipitation is included by Hofmann et al.,²⁰ allowing for the exploration of further Li–S features, such as low Coulombic efficiency and capacity fade.

While no existing model is able to reproduce all known features of the Li–S behaviour, one dimensional mechanistic models are the most promising in their ability to direct research towards improved cell performance. However, there are two major drawbacks to the use of such models for prediction in applications. Firstly, they require a large number of physical and chemical parameters whose values are difficult to obtain, such that their quantitative predictions are usually unreliable. More information helpful in determining the values of rate constants, equilibrium concentrations, or reaction pathways should become available as a result of advances in characterisation studies.⁸ A second drawback is the models' need for significant computational power, making them unsuitable as a basis for identification and control algorithms. Attempts to produce reduced order models derived directly from physical models of Li–S cells are expected to be extremely useful for control and application engineers.

For the purpose of Li–S cell simulation for use in applications, a suitable model needs to include the following functionality:

(i) retrieve the main features of a typical constant current performance, shown in Fig. 1; during discharge: the existence of two plateaus, with a dip in-between; during charge: the initial sharp increase in voltage, and a less pronounced difference between the two plateaus,

Fig. 1 Typical discharge and charge performance of a Li–S cell.

(ii) predict the cell response to dynamic loads, by allowing charge, discharge, and switching between the two,

(iii) include the effect of current on the voltage response, and thus account for power limitation and shuttle effects, and

(iv) provide information on the the amount of stored energy throughout operation, offering the possibility of including cell degradation and capacity fade.

We present a zero-dimensional model that aims to fulfill these requirements with relatively modest computational expense. The model is based on the simple two-step electrochemical reaction chain proposed by Mikhaylik and Akridge.¹⁰ The choice of a short reduction chain is also supported by density functional theory calculations,²¹ where the S²⁻₄ → S²⁻₂ and S²⁻₂ → S²⁻ reactions were found to have relatively similar standard potentials (2.22 V vs. 2.18 V). Reaction kinetics limitations are introduced via the Butler–Volmer relation. The shuttling of high order polysulfides is modelled via a constant shuttle rate, as in Mikhaylik and Akridge,¹⁰ such that it can take place both during charge and discharge and is only indirectly dependent on voltage. While precipitation can be a determining mechanism for the performance of Li–S cells, the debate regarding which species are precipitating alongside Li₂S is ongoing. Experimental characterisation studies usually infer the presence of Li₂S₂,²² while first principles studies predict various levels of stability for the Li₂S₂ compound.^23–25 In the present model, only the precipitation of the most reduced sulfur species is allowed, modelled via a constant precipitation rate including nucleation and the effect of a saturation concentration.

As a zero dimensional model, transport limitations cannot be retrieved. Despite this simplification, charge and discharge predictions are similar to those obtained from the more computationally intensive one dimensional models. This observation is supported by the conclusion of Ghaznavi and Chen¹⁷ that species transport becomes limiting in the one-dimensional model of Kumaresan et al.¹⁴ only in the limit of much lower diffusion coefficients than the values used.

Finally, it is worth mentioning that the concept of an open circuit voltage (OCV), which lies at the core of operational models for Li-ion intercalation batteries, is not well defined in the case of Li–S. Moreover, obtaining the OCV curve by using the same experimental procedure as for Li-ion intercalation batteries is not necessarily valid. Mikhaylik and Akridge¹⁰ used the common interpretation of the OCV as the discharge voltage curve at low enough constant current, in order to parametrise the standard potentials for the two redox reactions. This approach is not strictly valid in the upper plateau, if the phenomenon of shuttle is present during discharge. Precipitation was also shown to have a determining effect on the flatness of the lower voltage plateau.^14,26 Shuttle and precipitation/dissolution affect the experimentally obtained OCV for both charge and discharge in ways that need to be understood and quantified further. For Li–S cells, discharge/charge at decreasing current rates does not necessarily lead to an equilibrium state, or a state that is an appropriate input for operational models. As a result, it might not be possible to directly employ the same framework as that used for intercalation Li-ion for Li–S operational models, be they equivalent circuit models or mechanistic models. This situation makes the approach of equivalent circuit modelling, the tool of choice for engineers, especially unreliable when modelling Li–S cells. Indeed, all Li–S equivalent circuit models currently available in the literature are used for interpreting impedance spectroscopy data,^27,28 are developed for a quasi-static state and thus would not be appropriate for the simulation of a cell during operation.

2 Model

2.1 Equilibrium

In the present model, the same two-step reaction chain as in Mikhaylik and Akridge¹⁰ is used, under the assumption that a single electrochemical reaction dominates each of the two discharge regions:

S⁰₈ + 4e⁻ ↔ 2S²⁻₄, (1a)

S²⁻₄ + 4e⁻ ↔ 2S²⁻↓ + S²⁻₂. (1b)

While experimental evidence strongly suggests that the reaction pathways for charge and discharge might be different,²⁹ the exact mechanisms are not yet agreed upon. In the present model, the same reactions govern both charge and discharge. The form of sulfur in the fully charged cell depends on the chosen upper limit voltage or cell type, i.e. cathode vs. catholyte. We assume that the entire sulfur mass is in the form of dissolved S⁰₈. Alternative assumptions include solid S⁰₈ (ref. 14) and dissolved S²⁻₈,¹² the latter leading to reduced discharge capacity.

The equilibrium potential for the two reactions is given by the Nernst equations

where E⁰_H and E⁰_L are the standard potentials for eqn (1a) and (1b) respectively, R is the gas constant, F the Faraday constant, T the temperature, and S²⁻, S²⁻₂, S²⁻₄ and S⁰₈ the amounts of sulfur dissolved in the electrolyte in the respective forms in grams. The constants f_H and f_L in eqn (2a) and (2b) arise from the desire to calculate amounts of sulfur rather than concentrations. From the authors' understanding, this was not considered in Mikhaylik and Akridge,¹⁰ but the missing factor could prove essential when comparing model predictions to experimental data. The two constants take the formwhere v is the volume of electrolyte in the system, M_S8 the molar mass of sulfur, and the n-numbers represent the number of atoms per S⁰₈ molecule and S²⁻₄, S²⁻₂ and S²⁻ ions. Alternatively, f_H and f_L can be omitted from the Nernst equations, dictating that polysulfide quantities are calculated as concentrations, not masses in the system. If this is the case, the equations presented in Section 2.4 describing the time evolution of species must include an additional term accounting for the volume of electrolyte in the cell.

The voltage of the cell was approximated by Mikhaylik and Akridge¹⁰ as E_H = E_L = V. The first equality is valid for two reactions occurring simultaneously at the same electrode. The second equality is strictly valid only for the case of zero net current, when the forward and backward reactions are in equilibrium. This is not the case for a cell under charge/discharge.

The Nernst expressions in eqn (2) warrant a short discussion, as the inherent assumptions are sometimes overlooked. The Nernst equation as written in eqn (2a) and (2b) is strictly valid for a half cell reduction reaction, at the cathode. A battery, however, contains two half cells, since the electrochemical reactions take place at the two electrodes simultaneously. The expression of the anode potential as a result of Li oxidation is omitted in Mikhaylik and Akridge,¹⁰ under the assumption that its value is comparatively small. The Nernst potential of the reaction at the lithium anode is included in Kumaresan et al.,¹⁴ where its value is calculated to vary by one order of magnitude less than the variation of reaction potentials at the cathode. The small variation is a result of their model predicting a relatively constant average Li-ion concentration throughout cell operation. The prediction itself has not yet been validated. The charge transfer overpotential at the anode is not included. In the light of their estimations, the anode potential is ignored in the present model. The validity of this assumption, however, remains to be verified experimentally.

A second observation is that the Nernst potential in its original form contains the ratio of reactant to product activities rather than concentrations. For example, for the reaction in eqn (1a) the Nernst potential should read

with the powers of the terms under the logarithm corresponding to the stoichiometric coefficients. The activity of a species is proportional to its concentration as a = γc/c⁰, where γ is the dimensionless activity coefficient and c⁰ is the reference concentration, chosen to correspond to the standard half cell potential E⁰. The only way to obtain eqn (2) from eqn (4) is to make the following assumptions for all species considered:

(1) γ = 1, an assumption strictly valid at very low concentrations, and

(2) c⁰ = 1 mol L⁻¹, corresponding also to the concentration at which E⁰_L and E⁰_H should be measured.

While neither condition above is expected to be met for the Li–S system, these assumptions continue to be made here, as in previous models. Various models to approximate the activity coefficient at moderate and strong concentrations, such as extensions to the Debye–Hückel approximation³⁰ or the Pitzer's equations,³¹ have been derived mainly for aqueous systems. Their predictions have not been validated for the concentrated multispecies non-aqueous electrolyte solutions typical of the Li–S system, nor are appropriate parameter values available. While obtaining reliable E⁰ values might possible via computational simulations,²¹ many other parameter values are unlikely to be within experimental reach in Li–S cells, such that leaving them as fitting parameters may be a necessary lasting compromise.

2.2 Reaction kinetics

The presence of current in the outside circuit corresponds to a non-equilibrium state: electrochemical reactions with a net flow of electrons occur at the electrode. It is assumed that the currents related to the two electrochemical reactions are described by the Butler–Volmer approximation:Here n_e is the number of electrons transferred in each reaction, in this case in both reactions n_e = 4, i_H,0, i_L,0 are the exchange current densities, η_H, η_L the surface overpotentials for the two reactions, and a_r the active surface area available for the reaction, assumed constant. The exchange current densities depend, in general, on local properties at the interface, such as temperature, structure of the electrode surface, and composition of the solution. For a single electron reaction, this latter dependence can be obtained experimentally as a function of product and reactant concentrations. As such data is not available for the Li–S system, the exchange current density is assumed constant. The expressions in eqn (5) are written for the case in which the cathodic and anodic reactions are promoted equally, such that the symmetry factor is 0.5, in the absence of any evidence of the contrary. Finally, no double layer effects at the electrolyte/electrode interface are accounted for in the present formulation.

A non-zero surface overpotential is a measure of the driving force for a reaction to occur and is given by the difference between the apparent (measured) voltage of the cell, here of the cathode, and the reaction Nernst potential:

η_H = V − E_H (6a)

η_L = V − E_L. (6b)

The sign of the surface overpotential establishes the direction of the reaction.

Charge conservation dictates that the measured cell current I is given by the combined contribution of the two reactions

I = i_H + i_L. (7)

In this model, a positive i_H or i_L value corresponds to reduction, such that positive I corresponds to cell discharge.

2.3 Shuttle and precipitation

The shuttling of high order polysulfides and the precipitation of low order polysulfides are characteristics of Li–S cells that have been observed in various electrolytes. In the present model, the two effects are described in a conceptually similar manner. A simple but effective model for the shuttle with qualitatively good estimates of the resulting decrease in charge efficiency was proposed by Mikhaylik and Akridge.¹⁰ The effect is modelled by a shuttle rate k_s acting to decrease S⁰₈, in the same direction as the reduction reaction S⁰₈ → S²⁻₄, but without contributing to the reaction current.

Similar to shuttle, the precipitation reaction removes from electrolyte an amount of S²⁻ proportional to a precipitation rate k_p. There are three direct effects related to precipitation, regarding electrolyte conductivity, nucleation and cathode active surface area. Firstly, the fact that sulfur ions are removed from electrolyte causes a lowering of ionic concentration and results in a variation of the electrolyte conductivity.²⁶ The discharge capacity of the cell is reduced if the species precipitating could otherwise participate in further reduction reactions.¹⁴ In this model the sole precipitating species is the last one in the reacting chain, such that precipitation does not lead to reduced capacity during discharge. A second effect related to precipitation is the nucleation phenomenon: the presence of already precipitated material offers a nucleation surface for further precipitation to occur, as long as the concentration in the electrolyte is above a given saturation concentration. Here a simple model for the effect of nucleation is considered, which assumes that the precipitating amount is proportional to the amount of material already precipitated, as described by eqn (8e). Thirdly, the precipitated sulfur is electrically insulating; when it covers an active cathode/electrolyte interface, rather than agglomerate in solution or over already inactive cathode area, the active surface area is decreased. As further reactions are hindered from taking place, the reaction overpotential should increase, affecting, in turn, the exchange current term in the Butler–Volmer relation in eqn (5b). This last effect is ignored here, due to lack of experimental data. Unutilised sulfur due to pore blocking is a capacity fade phenomenon beyond the scope of this study.

2.4 Time evolution of species

The two electrochemical reactions allowed, together with the assumptions on the shuttle and precipitation phenomena described above, lead to the following time-evolution relations for the various sulfur species in the system:where S_p is the mass of precipitated sulfur, ρ_S its density and the saturation mass of S²⁻, assumed to be constant.

2.5 Computational implementation and initial conditions

Eqn (2) and (5)–(8) form a differential algebraic system that can be solved for the twelve unknowns: the Nernst potentials E_H, E_L, the contributions to the total current from the two reactions i_H, i_L, the cell voltage V, here equal to the cathode voltage, the overpotentials η_H, η_L, and the mass of the five forms of sulfur S⁰₈, S²⁻₄, S²⁻₂, S²⁻ and S_p. The Jacobian for the system is calculated analytically and the system is solved in Matlab using a second order solver. For the Matlab code please see the ESI.†

The initial conditions for all variables are calculated self consistently from chosen values of V, S⁰₈ and S_p for discharge, and V, S_p and S²⁻ for charge. If the initial split between i_H and i_L in eqn (7) is known, the initial conditions for the other parameters can be calculated as follows. For discharge, η_H is obtained from eqn (5a), E_H, from eqn (6a) and S²⁻₄ from eqn (2a). A similar sequence is followed for the low plateau reaction, with the term (S²⁻)²S²⁻₂ obtained from eqn (2b), and the assumption of S²⁻₂ = S²⁻ + S_p allowing for the calculation of S²⁻ and S²⁻₂. For charge, a similar procedure is used to obtain values for E_H and E_L; S²⁻₂ is calculated as above, and used in eqn (2b) to obtain S²⁻₄, which in turn is used in eqn (2a) to solve for S⁰₈. In this way, the initial conditions are calculated without the need of solving the non-linear system of equations. As a drawback, the total sulfur mass in the system is not specifically constrained such that the three initial values must be chosen to reflect a system with the desired total sulfur mass.

For an initial state far enough from the boundary between plateaus, the split of I can be assumed as I = i_H, i_L = 0 for discharge and I = i_L, i_H = 0 for charge. In the following analysis, only discharge from a fully charged state and charge from a fully discharged state are considered, making this assumption valid.

2.6 Model use

A fully operational model can be built on the basis of the model developed here. To this purpose, contributions from all Ohmic losses in the system should be subtracted from the voltage output shown in Fig. 2. These include losses caused by the contact resistance of the studied cell, as well as by the resistivity of the electrolyte as a function of the varying polysulfide concentration, shown in Zhang et al.²⁶ to be required for correct estimates of cell impedance.

Fig. 2 (a) Simulated discharge for constant current rates of 1.7 A (C/2, △) and 6.8 A (2C, ●) and (b) simulated charge for constant current rates of 1.7 A (C/2, △) and 3.4 A (1C, ●), for three levels of model complexity. During discharge, the model predicts power limitations and a flat lower plateau only by including a description of current overpotential and precipitation, respectively. During charge, slow dissolution of the precipitate can eliminate the distinction between the two plateaus and allow for a dramatic effect of charging current on cell capacity. Symbols distinguish between data sets rather than denote computed values.

While the two step mechanism assumed here predicts many known performance features, experimental studies indicate that Li–S charge/discharge occurs through more complex mechanisms.³² Further reactions and species can be included in the following manner, consistent to the current model. Electrochemical reactions are described by an equilibrium potential as given by their Nernst equation, and by a current–overpotential relation given by the Butler–Volmer approximation, where the overpotentials are linked to the potential of the whole cell, and the reaction currents sum up to the total external current. Chemical reactions could be described by forward and backward reaction rates. The time evolution of all species in the system leads to production and consumption as a result of all types of reactions, together with the effect of shuttle on higher order polysulfides and the effect of precipitation/dissolution on lower order polysulfides.

The choice of cell components, its geometry and size, all contribute to its behaviour. The choice of electrolyte, for example, is well documented to affect many properties inside the cell, impacting on many model parameters. The magnitude of the polysulfide shuttle, and thus of k_s, is affected by electrolyte choice and additives.³³ The choice of electrolyte also impacts on the free solvation energies of ions, and thus implicitly on the electrode potentials, or the Nernst potentials of the various reactions, and ultimately on the cathode and cell voltage.³⁴ The saturation concentration of the various species, determining the initiation of precipitation/dissolution, is also expected to change as a result of electrolyte choice. Moreover, the actual mechanisms leading to charge/discharge can be affected, as different solvents promote different reaction paths.³⁴ Finally, if capacity fade due to reactions occurring at the anode surface is included in the model, the effect of electrolyte additives promoting a solid electrolyte interphase would also need to be considered.³⁵ In fitting the model parameters to a particular cell, many have to be inferred rather than obtained from independent experiments. As a result, many of these effects would be included automatically.

3 Results and discussion

The system parameters and their values are given in Table 1, where the subscript H denotes the high plateau reaction eqn (1a) and L the low plateau reaction eqn (1b).

Table 1 Model parameter values

Notation	Name	Units	Value
Physical constants
F	Faraday's constant	C mol^−1	9.649 × 10^4
M S8	Molar mass S^08	g mol^−1	32
N A	Avogadro number	mol^−1	6.0221 × 10^23
n e	Electron number per reaction	—	4
n S8, nS4, nS2, nS	Number of S atoms in polysulfide	—	8, 4, 2, 1
R	Gas constant	J K^−1 mol^−1	8.3145
ρ S	Density of precipitated sulfur	g L^−1	2 × 10^3
Cell design properties
a r	Active reaction area per cell	m^2	0.960
f H	Dimensionality factor H	g L mol^−1	0.7296
f L	Dimensionality factor L	g^2 L^2 mol^−1	0.0665
v	Electrolyte volume per cell	L	0.0114
m S	Mass of active sulfur per cell	g	2.7
Kinetic properties
E ^0H	Standard potential H	V	2.35
E ^0L	Standard potential L	V	2.195
i H,0	Exchange current density H	A m^−2	10
i L,0	Exchange current density L	A m^−2	5
Shuttle and precipitation parameters
	S^2− saturation mass	g	0.0001
k p	Precipitation rate	s^−1	100
k s	Shuttle constant	s^−1	0.0002
Operational parameters
I	External current	A	Variable
T	Temperature	K	298
Variables
E H, EL	Nernst potentials	V
i H, iL	Current contributions	A
η H, ηL	Overpotentials	V
V	Cell voltage	V
S^08, S^2−4, S^2−2, S^2−	Mass of dissolved sulfur	g
Sp	Mass of precipitated S^−2	g

---

3.1 Discharge and charge

Model predictions for the cell voltage during constant current discharge for two current rates are illustrated in Fig. 2(a). In the absence of precipitation and overpotential effects, the outputs of Mikhaylik and Akridge¹⁰ are retrieved (ppt no, BV no). In this case, the magnitude of the discharge current affects the cell behaviour only through its interplay with the shuttle: the higher the current, the higher the discharge efficiency in the high plateau, and correspondingly the higher the total discharge capacity. The value of the current has no effect on the height of the two voltage plateaus.

In reality, the magnitude of the discharge current is expected to impact on the discharge curve even in the absence of a shuttle effect; a cell cannot provide infinite instantaneous power due to losses from diffusion and charge transfer reactions. To account for the latter, as the discharge current is increased, greater losses should be caused by an increasing reaction overpotential. Predictions from the model including the voltage–overpotential relation in Fig. 2(a) indeed show that the accessible cell voltage is significantly lower than that predicted by the Mikhaylik model. The total capacity of the cell as defined by the state at which all available sulfur has reacted down the reaction chain to low order polysulfides does not change upon the addition of precipitation of end products and of current–overpotential relations. In practice, however, the capacity of a cell is defined by the charge available within a given voltage range. While the theoretical capacity of the cell remains the same as in the simpler model, in the model including kinetic limitations the voltage of the lower plateau decreases as the discharge current increases. As a result, the cell provides less usable capacity within a set voltage range, even before accounting for Ohmic losses.

The discharge curve predicted by the model with kinetic limitations but without precipitation in Fig. 2(a) (ppt no, BV yes) is still qualitatively different from those seen in typical experimental results, where the lower plateau is significantly flatter. The final set of discharge curves plotted in Fig. 2(a) (ppt yes, BV yes) exhibit the full model capabilities, including the effect of precipitation. The dip between the two voltage plateaus and the elevated flat voltage in the lower plateau are both revealed as effects of precipitation. The effect of precipitation is further analysed in Section 3.2.

Predictions of cell voltage during constant current charge are given in Fig. 2(b). For the simplest model ignoring both kinetic limitations and precipitation effects, the charging rate only affects the voltage at the end of charge, counteracting the effect of shuttle. During charge, the addition of a current-dependent reaction overpotential to the model by Mikhaylik and Akridge¹⁰ introduces an increase in cell voltage, effectively decreasing the charge capacity for a set voltage interval. The addition of precipitation further increases the cell voltage by introducing a large resistance associated with the dissolution of S_p. For a larger charging current, less S_p can be dissolved within the shorter time frame, increasing the cell voltage and leading to a charge capacity that is significantly reduced. As in the case of discharge, the charge voltage curve becomes similar to experimental data upon inclusion of precipitation effects. For charge, the presence of gradual dissolution removes the clear boundary between the two plateaus, as detailed in Section 3.2. In all cases presented, the charging current is strong enough to compensate for the shuttle, such that the cell can reach the nominal voltage for fully charged.

It should be noted that the choice of initial conditions for Fig. 2(b) is not unique. In comparing predictions for charge from a fully discharged state between a system where precipitation is allowed and one without precipitation, choosing the initial conditions is not straightforward. When fully discharged, an ideal cell following the present model would contain all sulfur in the form of S²⁻₂, S²⁻ and S_p in ratios influenced by the two precipitation parameters, k_p and . In the absence of precipitation, the fully discharged cell contains S²⁻ and S²⁻₂ only, in a 1 : 1 mass ratio. For a given total sulfur mass, these two conditions do not correspond to the same value of cell voltage. For comparison reasons, the same starting voltage was chosen for the predictions in Fig. 2(b). As a result, the predictions correspond to cells of different sulfur masses, with a variation of 2.898 × 10⁻⁴ g, or 0.364 × 10⁻³ Ah in their theoretical cell capacity, which was considered small enough to be ignored. This particular problem does not arise for discharge from a fully charged cell, where it is assumed that all sulfur is in the form of S⁰₈ for all model variants. A similarly small current-dependent variation in the total sulfur mass due to the way in which initial conditions are calculated arises when comparing predictions from the model including kinetic limitations to those from the model ignoring them. Either different values of i_H,0, i_L,0 but the same i_H, i_L, or different i_H, i_L, such as due to different current rates, but fixed i_H,0, i_L,0 cause a variation in η_H, η_L. For the discharge curves in Fig. 2(a), the maximum variation in m_S corresponds to 30.8 × 10⁻³ Ah.

3.2 Effect of precipitation

The presence of precipitation can significantly impact the shape of the discharge and charge voltage curves. In general, the occurrence of precipitation raises the voltage in the low plateau, compared to that given by a Nernst model, as seen for various values of the precipitation parameters in Fig. 3(a) and (b). The cell voltage in the low plateau traces the evolution of E_H in eqn (2a), if the assumption η_H ≈ 0 can be made. During discharge, as S²⁻ is gradually extracted from solution by precipitation, the low plateau reaction is driven forward and more S²⁻₄ is consumed. This in turn raises the Nernst voltage of the high plateau reaction, and thus also V, as observed also in Zhang et al.²⁶

The value of determines the state during discharge at which the precipitate starts forming, and thus affects the width of the voltage dip between the high and the low plateau. Higher values of allow for a later start of precipitation, an effectively faster initial rate of precipitation and a larger amount of S²⁻ in the electrolyte at any point during the discharge, as shown in Fig. 3(c). The latter has a significant effect on the value of the voltage in the low plateau: the smaller the amount of S_p, the lower the voltage, corresponding to a smaller departure from the Nernst voltage.

Fig. 3 Simulated constant current discharge and charge curves at 1.7 A (C/2) for various precipitation parameter values (a and b) and the corresponding evolution of the precipitated species (c and d). For the discharge, the x-axis in (c) corresponds to the grey interval in (a). As the precipitation of S²⁻ becomes more prevalent with increased k_p and decreased , the voltage in the lower plateau departs form the Nernst voltage, increasing and becoming flatter. Symbols distinguish between data sets rather than denote computed values.

Increasing the precipitation rate k_p has qualitatively similar effects on the voltage curve to decreasing . The effect of k_p is modulated by : the larger the latter, the less S_p is formed and the equilibrium is reached faster, lowering the importance of k_p. In the limit of infinitely fast precipitation, the cell voltage in the low plateau is determined by only. This cell voltage corresponds to the equilibrium voltage, similar to the concept of the open circuit voltage for intercalation Li-ion batteries.

In the absence of current overpotential and precipitation effects, I/k_s was shown to be a scaling factor for the voltage response of the cell in the model of Mikhaylik and Akridge.¹⁰ Similarly, in the absence of shuttle and overpotential effects, I/k_p can be considered the determining factor, due to the way in which shuttle and precipitation are described in competition with the applied current. As at least two of the three contributions, i.e. reaction overpotential, precipitation and shuttle, are usually present at any time during operation, the I/k_s and I/k_p factors are only indicative of the performance of a real cell, rather than serving as scaling factors. For relatively low overpotentials, I/k_s is indicative of the discharge voltage in the high plateau, while I/k_p serves this purpose in the low plateau.

As I/k_p is performance-determining for fixed in the low plateau, an experimental OCV curve for a Li–S cell can only be obtained for a constant current discharge at currents much lower than the rate of precipitation. A value of the latter is not readily available from literature, but could indicate a prohibitively low recomended value of I. Rest periods during discharge at rates as low as C/48 were found to lead to significant voltage recovery.³⁶ In the high plateau, moreover, shuttle would affect the discharge. In the model by Mikhaylik and Akridge,¹⁰ the voltage response to a C/30 discharge is considered to be the OCV. E⁰_H and E⁰_L are read directly from the data, at the points in the discharge capacity that correspond to those species concentration ratios that cancel the logarithmic terms in the Nernst equations, such that E_H = E⁰_H and E_L = E⁰_L respectively. The values thus obtained differ from simulation predictions,²¹ possibly due to automatically encompassing non-vanishing overpotentials from charge transfer and precipitation effects.

In order to analyse the effect of dissolution during charge, identical initial conditions are chosen for the four charging systems in Fig. 3(b), for an easier comparison. In reality, it is expected that, if charged from fully discharged, the different cells would start from different initial voltages and amounts of S_p. During charge, higher k_p and lower values result in a slower dissolution. Unlike for discharge, where the theoretical cell capacity is independent of precipitation rate, the charge capacity does depend on the system's ability to dissolve the precipitated species fast enough, as this ability impacts sulfur utilisation. During charge, a slow dissolution of S_p can act as a bottleneck, leading to the early formation of high plateau species, here S⁰₈, and to reduced capacity. The effect of this bottleneck on the cell voltage is already visible in Fig. 2(b) (ppt yes, BV yes) for charge at two different currents. The mechanism of the bottleneck becomes apparent when tracking the contribution of the two electrochemical reactions to the total current, shown in Fig. 4(a), and the evolution of species in the system, in Fig. 4(b). Thus, in the presence of slow dissolution, the boundary between the high and low plateaus becomes less sharp, leading to the disappearance of a distinct low plateau during charge. A slower dissolution rate during charge also corresponds to a larger initial jump in voltage, and a higher voltage during charge, as seen in Fig. 3(b).

Fig. 4 (a) Simulated reaction currents and (b) species evolution during a 1.7 A (C/2) constant current charge with high saturation concentration = 0.005 g, k_p = 100 1/s, (red △) and low saturation concentration = 0.0001 g, k_p = 100 1/s (blue ●). The corresponding voltage is illustrated in Fig. 3(b) with the same symbols. The slower dissolution acts as a bottleneck during charge, triggering an earlier onset of the high plateau reaction and eliminating the boundary between the plateaus. Symbols distinguish between data sets rather than denote computed values.

4 Conclusions

We describe a zero dimensional model for predicting the voltage response of a Li–S cell under the assumption of a two-step electrochemical reaction chain. By including kinetic limitations according to the Butler–Volmer approximation and a simple model of precipitation/dissolution, important features of the discharge and charge voltage curves are reproduced, such as the flat low voltage plateau and the dip between voltage plateaus during discharge, and the possible merging of the two plateaus when dissolution is the limiting phenomenon during charge. In the present model, the precipitating species does not participate in further reduction reactions, such that the theoretical discharge capacity of the cell does not change upon including precipitation effects. Dissolution, however, can become the limiting process for high charging currents, reducing the cell capacity.

Comparative to the existing one dimensional models,^14,19 the model presented here offers a tool for analysing the evolution of species concentration and the interplay between reaction and precipitation parameters with significantly lower computational requirements. Beyond the analysis of mechanisms inside a Li–S cell and their effects, this model is a suitable platform for predictive and diagnostic methods used by control engineers in evaluating the performance of Li–S cells in various applications.

The model also enables the interpretation of an open circuit voltage for Li–S, and indicates that such experimental data can be obtained from constant current discharge only at prohibitively low currents. The approach of using Nernst equations to model the current-free state can eliminate the need to define and obtain the open circuit voltage for a Li–S cell. In exchange, however, parameter values must be obtained from other modelling and experimental data, such as standard reaction potentials from simulations, but also precipitation and shuttle parameters.

The model can be easily expanded to include the voltage contribution from the anode and that from the concentration-dependent electrolyte resistivity, as considered in Zhang et al.,²⁶ for the purpose of comparison to experimental data. Other features that might play a role in the mechanism of Li–S cells, such as a longer chain of electrochemical reactions,¹² the existence of chemical reactions, that of parallel electrochemical reaction paths, or the precipitation/dissolution of other species can also be included. Each of these additions increases the number of parameters that must be obtained from experimental data, such that the need of any such improvement should be based on experimental evidence of the feature's importance in the studied cell.

Acknowledgements

The authors would like to thank the Engineering and Physical Sciences Research Council in the UK for funding this work under the Revolutionary Electric Vehicle Battery (REVB) project EP/L505298/1, and acknowledge useful discussions with Mark Wild and Laura O'Neill from OXIS Energy.

References

1. P. G. Bruce, S. A. Freunberger, L. J. Hardwick and J.-M. Tarascon, Nat. Mater., 2012, 11, 19–30.
2. J. Brückner, S. Thieme, F. Böttger-Hiller, I. Bauer, H. T. Grossmann, P. Strubel, H. Althues, S. Spange and S. Kaskel, Adv. Funct. Mater., 2014, 24, 1284–1289.
3. M. Hagen, D. Hanselmann, K. Ahlbrecht, R. Maça, D. Gerber and J. Tübke, Adv. Energy Mater., 2015, 5, 1401986.
4. Y.-S. Su, Y. Fu, T. Cochell and A. Manthiram, Nat. Commun., 2013, 4, 2985.
5. I. A. Hunt, Y. Patel, M. Szczygielski, L. Kabacik and G. J. Offer, J. Energy Storage, 2015, 2, 25–29.
6. A. Manthiram, Y. Fu, S.-H. Chung, C. Zu and Y.-S. Su, Chem. Rev., 2014, 114, 11751–11787.
7. D. Bresser, S. Passerini and B. Scrosati, Chem. Commun., 2013, 49, 10545–10562.
8. R. Xu, J. Lu and K. Amine, Adv. Energy Mater., 2015, 1500408.
9. R. M. Spotnitz, Electrochem. Soc. Interface, 2005, 14, 39–42.
10. Y. V. Mikhaylik and J. R. Akridge, J. Electrochem. Soc., 2004, 151, A1969.
11. Y. V. Mikhaylik and J. R. Akridge, J. Electrochem. Soc., 2003, 150, A306.
12. D. Moy, A. Manivannan and S. R. Narayanan, J. Electrochem. Soc., 2014, 162, A1–A7.
13. S. Risse, S. Angioletti-Uberti, J. Dzubiella and M. Ballauff, J. Power Sources, 2014, 267, 648–654.
14. K. Kumaresan, Y. Mikhaylik and R. E. White, J. Electrochem. Soc., 2008, 155, A576.
15. M. Ghaznavi and P. Chen, J. Power Sources, 2014, 257, 394–401.
16. M. Ghaznavi and P. Chen, J. Power Sources, 2014, 257, 402–411.
17. M. Ghaznavi and P. Chen, Electrochim. Acta, 2014, 575–585.
18. J. P. Neidhardt, D. N. Fronczek, T. Jahnke, T. Danner, B. Horstmann and W. G. Bessler, J. Electrochem. Soc., 2012, 159, A1528–A1542.
19. D. N. Fronczek and W. G. Bessler, J. Power Sources, 2013, 6–11.
20. A. F. Hofmann, D. N. Fronczek and W. G. Bessler, J. Power Sources, 2014, 259, 300–310.
21. L. Wang, T. Zhang, S. Yang, F. Cheng, J. Liang and J. Chen, J. Energy Chem., 2013, 22, 72–77.
22. S. Waluś, C. Barchasz, R. Bouchet, J.-C. Leprêtre, J.-F. Colin, J.-F. Martin, E. Elkam, C. Baehtz and F. Alloin, Adv. Energy Mater., 2015, 5, 1500165.
23. H. Park, H. S. Koh and D. J. Siegel, J. Phys. Chem. C, 2015, 119, 4675–4683.
24. G. Yang, S. Shi, J. Yang and Y. Ma, J. Mater. Chem. A, 2015, 3, 8865–8869.
25. P. Partovi-Azar, T. D. Kühne and P. Kaghazchi, Phys. Chem. Chem. Phys., 2015, 22009–22014.
26. T. Zhang, M. Marinescu, L. O'Neill, M. Wild and G. Offer, Phys. Chem. Chem. Phys., 2015, 17, 22581–22586.
27. Z. Deng, Z. Zhang, Y. Lai, J. Liu, J. Li and Y. Liu, J. Electrochem. Soc., 2013, 160, A553–A558.
28. V. Kolosnitsyn, E. Kuzmina and S. Mochalov, J. Power Sources, 2014, 252, 28–34.
29. Q. Wang, J. Zheng, E. Walter, H. Pan, D. Lv, P. Zuo, H. Chen, Z. D. Deng, B. Y. Liaw, X. Yu, X. Yang, J.-G. Zhang, J. Liu and J. Xiao, J. Electrochem. Soc., 2015, 162, A474–A478.
30. G. Scatchard, Chem. Rev., 1936, 19, 309–327.
31. K. S. Pitzer, J. Phys. Chem., 1973, 77, 268–277.
32. M. Cuisinier, P.-E. Cabelguen, S. Evers, G. He, M. Kolbeck, A. Garsuch, T. Bolin, M. Balasubramanian and L. F. Nazar, J. Phys. Chem. Lett., 2013, 4, 3227–3232.
33. N. Azimi, Z. Xue, L. Hu, C. Takoudis, S. Zhang and Z. Zhang, Electrochim. Acta, 2015, 154, 205–210.
34. H. Schneider, C. Gollub, T. Weiss, J. Kulisch, K. Leitner, R. Schmidt, M. M. Safont-Sempere, Y. Mikhaylik, T. Kelley, C. Scordilis-Kelley, M. Laramie and H. Du, J. Electrochem. Soc., 2014, 161, A1399–A1406.
35. M. Barghamadi, A. S. Best, A. I. Bhatt, A. F. Hollenkamp, P. J. Mahon, M. Musameh and T. Rüther, J. Power Sources, 2015, 295, 212–220.
36. C. E. Parfitt, PhD thesis, University of Warwick, 2012.

---

Footnote

† Electronic supplementary information (ESI) available: Matlab code solving the system of equations and producing the plots. See DOI: 10.1039/c5cp05755h

---