Multi-Scale Electrolyte Transport Simulations for Lithium Ion Batteries

1Dassault Systèmes, CB4 0WN Cambridge, United Kingdom 2Dassault Systèmes, IDEON Gateway, 22363 Lund, Sweden 3Mathematical Sciences, University of Southampton, SO17 1BJ, United Kingdom 4The Faraday Institution, Harwell Campus, Didcot, OX11 0RA, United Kingdom 5Department of Chemistry, University of Cambridge, CB2 1EW, United Kingdom 6Department of Engineering, University of Cambridge, CB2 1PZ,United Kingdom

Commercially available lithium ion batteries (LIB) for electric vehicles and consumer goods applications are typically based on Li ion chemistry with an organic liquid electrolyte. 1 An automotive roadmap for further development includes electrolyte chemistries and formulations that are non flammable, non-toxic, and environmentally friendly, without compromising the overall battery performance. 1 To achieve this goal, it is necessary to link the molecular composition and microstructure of the battery components to the charge transport characteristics of the electrochemical cell, battery pack, and even vehicle performance.
Most individual battery components are reasonably well understood. On molecular scales, the physical properties of individual electrolyte constituents and additives can be designed using experimental and computational techniques. [2][3][4][5][6][7][8][9] Commercially available electrolyte formulations generally contain the salt LiPF 6 , 10,11 which is combined with carbonate solvents. A mixture of additives is also included for different purposes, 12,13 such as redox shuttle, 7 HF scavenger, 8 or to help form the solid electrolyte interphase (SEI). 9 All additives impact the transport properties of the overall formulation, thus directly influencing the electrochemical and overall device performance.
This observation is the starting point for our contribution. We establish a direct computational link between the electrolyte chemistry and battery cell performance. To this extent, we present predictions of transport properties for liquid electrolytes starting with molecular dynamics simulations. These results are used in a previously established LIB model, 10,11 which quantitatively predicts measured discharge voltage curves across a wide range of typical operating conditions. We note in passing that we explicitly spell out battery cell throughout this paper to avoid confusion with the simulation cell in our molecular modeling work.
For a given electrolyte, transport coefficients can be measured directly as a function of temperature and salt concentration. 10,11,[14][15][16][17][18][19] These experiments are frequently used to parameterize electrochemical cell models. The qualitative description of lithium diffusion and electrolyte conductivity has also been achieved using different molecular dynamics techniques. [20][21][22][23][24][25][26] Simulations of battery cell behavior are commonly performed with Newman models and are usually combined with measuring some of a Present Address: Laboratory for Computational Science and Modeling, École Polytechnique Fédérale de Lausanne, 1015 Lausanne, Switzerland. z E-mail: felix.hanke@3ds.com the required material parameters. 11,17,18,[27][28][29] The models contain sufficient detail to allow investigations of different charge and discharge regimes, including fast, cold, and partial charging, all of which are required for the development of entire battery packs and battery management systems. However, it remains challenging to fully parameterize Newman models based exclusively on the physical properties of the cell components, whether calculated or measured directly.

Approach for the Electrolyte Transport Simulations
In electrochemical cell models, the charge transport properties of the bulk electrolyte are described using the ionic conductivity κ, the Li diffusion coefficient D + and the Li transference number t + which denotes how much of the total current in the solution is carried by Li ions. Our first task is to calculate these parameters as a function of concentration for electrolyte solutions.
About 200 solvent molecules and an appropriate number of ions for each salt concentration are packed randomly into a simulation box using the Amorphous Cell module in BIOVIA Materials Studio 30 and their structure is then optimized. 31 The COMPASSII force field 32 is highly suitable for condensed phases of organic compounds and will be used throughout. The Amorphous Cell setup provides an energetically favorable liquid-like structure at experimentally known densities of around 1.2 g/cm 3 as illustrated in Fig. 1.
Each structure is thermalized with a 100 ps NPT molecular dynamics (MD) run using the Nosé-Hoover thermostat and a Berendsen barostat at zero pressure. We then continue with a 5 ns MD run in the NVE ensemble to determine the equilibrium transport properties for each starting configuration. Each simulation is repeated independently five times and the results are averaged to get the final transport properties of the electrolyte at a given temperature and concentration. Fig. 1 also shows the Li-O radial distribution functions (RDF) obtained from our simulations, which have a well-known characteristic double-peak signature of the PC solvation shell around each Li atom at all concentrations studied.
For each MD simulation, the diffusion coefficients D +/− are calculated from the average mean square displacement (MSD) of the ions. The MSD is computed by averaging over all available time differences up to 2.5 ns. From this data, the bulk diffusion coefficient is established by fitting a line to the upper half of the MSD time interval, e.g. from 1.25 ns to 2.5 ns. Typical R 2 values for the linear fit are well above 0.99, which means that the long-term diffusion limit has been reached. From the diffusion coefficients D for both Li and PF 6 ions, the conductivity κ at concentration c is calculated using the Nernst-Einstein equation as where T is the temperature, e is the elementary charge, k B is the Boltzmann constant, and N A is Avogadro's number. The bulk electrolyte conductivity is obtained from summing over ionic species and the transference number is given by At very low concentrations, this procedure can yield significant error bars due to limited statistics. Since these concentrations have no technological relevance, our approach provides a pragmatic parameter-free method to predict charge and mass transport properties. To prepare our MD data for the macroscale Newman model, we fit the conductivity using the following Ansatz. Asymptotically, the diffusion coefficient for ions in solution decreases exponentially with the concentration, D ∼ exp(−Bc), where B is a constant, see e.g. Ref. 17. The Nernst-Einstein relation 1 links D to the conductivity with a standard Arrhenius factor providing the temperature dependence. We introduce the following analytic function to describe the overall conductivity as a function of concentration and temperature: The numbers A, B, and the activation energy E a are obtained by fitting to the calculated conductivity, accounting for the standard error from the simulations. R is the ideal gas constant. Eq. 3 along with the predicted transference number t + forms the basis of our macroscale description of the electrolyte transport properties. The Li diffusion coefficient is obtained from 1. The conductivity function introduced in Eq. 3 scales linearly with concentration in dilute electrolyte solutions. 33 It decays exponentially for large concentrations as observed experimentally. 17 Unlike the polynomial expansions often found in the literature, it provides the correct scaling in the high and low salt concentration limits, and contains the correct temperature dependence. Moreover, the three parameters required in the formula 3 to describe both the temperature and concentration dependence are substantially fewer than what is used for any polynomial fit function.
We now move up in scale and use our transport parameterization to model an entire battery cell. Calculations are carried out with pseudo two-dimensional electrochemical Newman model 11,[17][18][19][27][28][29] implemented in the Battery Library in Dymola, a Modelica-compliant solution for modeling and simulation of integrated systems. 34 These simulations consider the lithiation/delithiation through diffusion of ions in all battery components as well as transport phenomena such as charge transfer and conduction in the solid active material of the electrodes and heat exchange. The model for the active material, the kinetics of the Li intercalation processes, and the manufacturing-dependent electrode and separator tortuosities are adopted from Ref. 11. Both the diffusion coefficient in the solid active material and the exchange current density at the SEI are dependent on the local lithium ion concentration and temperature. The electrolyte is described using the atomistic approach developed above.
A thermal model considering heat generation and storage in the battery cell and dissipation into the surroundings describes the temperature of a battery cell in contact with a heat bath at ambient temperature T amb as mc batteryṪ = CCC (T amb − T ) +Q r +Q j , [4] where dots denote time derivatives, m is the mass of the battery, and c battery its specific heat capacity. The reaction heat currentQ r is calculated from summing over the product of overpotential and the local current between electrolyte and electrode particles, while the joule heat currentQ j describes ohmic losses in the solid active material and the electrolyte. The recently introduced Cell Cooling Coefficient (CCC) 35 denotes the overall ability of the cell to dissipate heat to its surroundings.

Results
Our approach to the bulk electrolyte transport parameters is initially validated for LiPF 6 in propylene carbonate (PC). The results in Fig. 2 show excellent agreement between our parameter-free MD predictions for conductivity and diffusion coefficient with experimental results. 14 It is worth noting that the chiral nature of the PC molecule has no influence on the results. The calculations presented in Fig. 2 were performed with an even mixture of both enantiomers, but separate simulations with enantiopure PC yielded the same results. For high ionic concentrations, the Nernst-Einstein Equation 1 can overestimate the conductivity as it assumes independent charge carriers. Given the good correspondence between computation and measurement, this does not appear to be an issue here.
Encouraged by our single component analysis, we now turn to a multi-component commercial electrolyte used in Refs. 10 and 11.  we computed the anion/cation coordination numbers by integrating the radial density functions between P and Li. For technologically relevant concentrations below 1.5 mol/l the coordination remains below 2 up to a distance of 5 Å, suggesting that there is no aggregation and cluster formation at these concentrations and adding to the confidence in our transport coefficient calculations based on the formalism presented above. The results for the transport coefficients are shown in Fig. 3 along with the analytic form Eq. 3. Our predicted conductivity compares well with experimental results, slightly underestimating it at room temperature for the experimentally established concentration range between 0.5 and 1.5 mol/l. 10 The activation energy E a is found to be 14.2 kJ/mol, within chemical accuracy (4 kJ/mol) of the measured value of 17.1 kJ/mol. 10 Our prediction of the transference number t + is slightly overestimated when compared to the values of 0.35-0.4 found in Ref. 17. t + is more sensitive to correlated ionic motions than the other transport properties, 37 but given the overall close match of the conductivity and the thermal activation energy we are confident about the overall utility of our approach for liquid electrolytes.
Finally, we use the transport relations established in Fig. 3 to compute the discharge curves of battery cells using our extended Newman model. The parameters in Eq. 4 were obtained specifically for the commercially available battery used experimentally in Refs. 10 and 11. This pouch cell from the manufacturer Kokam has a capacity   of 7.5 Ah and a surface area of A s = 2 × (102 mm×107 mm). By matching to the measured discharge curve at current rate 1 C and temperature 263 K, we obtain a heat transfer coefficient of 2.3 W/m 2 /K. This corresponds to a cell cooling coefficient of 5 W/K which we use throughout the remainder of the simulations. The heat capacity was fixed at 1008 J/kg/K, as measured in Ref. 35 for this particular battery. Fig. 4 presents results for different temperatures and discharge rates, covering typical operating conditions in an automotive environment. The insets in each panel show the predicted battery temperatures, illustrating how operation at room temperature or higher and discharge rates of 1 C (e.g. the current required to fully deplete the cell in 1 h) or lower leads to very little heating of the battery. Our voltage traces correspond very well to experiments, suggesting that the underlying material parameters are of sufficient quality to explain the battery operation in those conditions. Fast discharge (up to discharge rates of 5 C) or operation down to a temperature of 263 K leads to substantial dissipation of heat inside the battery and a corresponding increase in its temperature. Test calculations at constant temperature underestimate the capacity by up to 25%, confirming that operation at low temperatures is highly influenced by heat generation and dissipation.
We note that separate simulations using the measured electrolyte transport properties from Ref. 10 are very similar to the results presented in Fig. 4. This corroborates our ansatz that an entire in-silico approach to the electrolyte formulation will be a successful strategy to identify suitable additive components going forward.

Discussion
In this contribution we have presented a seamless multi-scale link between the molecular composition of liquid electrolytes in commercially available batteries and electrochemical performance. Detailed molecular dynamics simulations establish the transport properties of the solution across a broad range of concentrations and temperatures. The Nernst-Einstein relation 1 along with Eq. 3 and four constants obtained directly from simulation data describe the bulk electrolyte without any phenomenological parameters. These functions are used in a Newman model and provide excellent agreement with experimentally observed discharge curves for across the typical operating conditions.
Our workflow is fully automated for a given electrolyte composition, salt concentration, and temperature. This opens a route toward the computational exploration of chemical space for the electrolyte, for example to establish the impact of individual additives on the macroscopic cell performance. The detailed optimization of the electrolyte is generally performed at a late stage in the battery cell development process, once the basic electrode chemistry has been established. Changes in the electrolyte formulation then can be implemented on engineering time scales, which are typically much shorter than the lead times involved in novel materials discovery.
In the near to intermediate future, the roadmap for automotive battery development calls for substantial focus on liquid electrolytes, with respect to flammability, degradation, and environmental impact. 1 The workflow presented here can speed up these developments considerably, for example by developing the electrolyte chemistry based on a fully parameterized reference battery model.
Finally, we highlight that a key ingredient to the electrochemical Newman model calculations is the operating temperature, which can be substantially higher than the ambient temperatures for lowtemperature operation or fast discharge rates. Here we use a simple setup to obtain the temperature dependence qualitatively. The thermal properties of the immediate environment of the battery therefore need to be an integral part of describing the batteries operating conditions if fully predictive models are to be successful.