A Multiscale Investigation on Electrolyte Systems of [(Solvent +Additive

19 hours ago - This website uses cookies to improve your user experience. By continuing to use the site, you are accepting our use of cookies. Read th...
0 downloads 0 Views 2MB Size
Subscriber access provided by Nottingham Trent University

C: Energy Conversion and Storage; Energy and Charge Transport

A Multiscale Investigation on Electrolyte Systems of [(Solvent +Additive) +LiPF] for Application in Lithium-Ion Batteries 6

Behnam Ghalami Choobar, Hamid Modarress, Rouein Halladj, and Sepideh Amjad-Iranagh J. Phys. Chem. C, Just Accepted Manuscript • DOI: 10.1021/acs.jpcc.9b04786 • Publication Date (Web): 15 Aug 2019 Downloaded from pubs.acs.org on August 15, 2019

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 55 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

A Multiscale Investigation on Electrolyte Systems of [(Solvent +Additive) +LiPF6] for Application in Lithium-Ion Batteries Behnam Ghalami Choobar1, Hamid Modarress1,*, Rouein Halladj1, Sepideh Amjad-Iranagh1 1

Department of Chemical Engineering, Amirkabir University of Technology (Tehran Polytechnic), Tehran, Iran . *Corresponding author. Email: [email protected] Tel: +982164543176 Fax: +982166405847

1 ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Abstract

To investigate how the electrolyte system of [(solvent+ additive) + lithium salt LiPF6] can promote the performance of lithium-ion batteries (LIBs), a multiscale method consisting of density functional theory (DFT) calculations at Pico-Nano scale level, molecular dynamics (MD) simulation at Nano scale level and pseudo 2-dimensional (P2D) electrochemical modeling at macro scale level were employed. The solvent used in the studied electrolyte systems was ethyl methyl carbonate (EMC) and the additives were: ethylene carbonate (EC), fluoroethylene carbonate (FEC) and succinonitrile (SN). The results indicated that at moderate temperatures the electrolyte systems [(EMC+SN 10)+LiPF6] and [(EMC+FEC 10)+LiPF6] represent higher performance compared with the conventional electrolyte [(EMC+EC 30)+LiPF6], since they create lower ohmic polarization in the LIB operation. Also the molecular configuration of the electrolytes were obtained by DFT calculations and their structural and transport properties were evaluated by MD simulations. To determine the operational properties of LIB such as voltage and polarizations at macro scale level, the evaluated properties by DFT and MD methods were inserted, as the basic parameters, for solving the continuity equations in P2D model. In this way we achieved an insight into relation between molecular structure and transport properties of electrolytes which can influence the performance of a LIB where it can be useful in designing novel LIBs.

2 ACS Paragon Plus Environment

Page 2 of 55

Page 3 of 55 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

1. Introduction By industrial development, the total power demand has increased constantly; such that total energy consumption reached from 64 billion watts in 1971 to 160 billion watts in 2016 and it is expected to triple by 20501. The provision of appropriate storage systems is essential to utilize the discontinuous energy sources, e.g. solar energy, wind and wave power. Up to now, lithium-ion batteries (LIBs) have assisted the technological revolution of portable electronic devices. Unfortunately application of current LIBs in medium- or large-scale systems like electrical vehicles is limited to short distances due to low energy densities, low charge rates and safety concerns2–4. The gaps between required and actual performance indicate the need for improving the present technology and also encourage the researchers to develop the new energy storage systems. The most important component of the LIBs is the electrolyte which plays a vital role in performance of the LIBs. These electrolytes, typically, composed of [(solvent + additive) + lithium salts)]. The best performance of LIBs depends on how effectively electrolyte interact with the electrodes3,5–7. The conventional electrolyte system used widely in LIBs is [(EMC + EC) + LiPF6] where ethyl methyl carbonate (EMC) is the solvent and ethylene carbonate (EC) is the additive. However this electrolyte system has some shortcomings such as high viscosity, low conductivity, reductive/oxidative instability and forming non-uniform solidelectrolyte interface (SEI) on the anode surface3,5–12. Therefore instead of EC other additives have been examined by researchers8–11. Our survey of literature revealed that the electrolyte system consisting of (solvent + additive) such as (ethyl methyl carbonate (EMC) + fluoroethylene carbonate (FEC)) 8,9 and (ethyl methyl carbonate (EMC) + succinonitrile (SN)) 10,11 have represented some favorable advantages; namely the former makes SEI with uniform morphology and the latter 3 ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

manifests higher oxidation potentials3,5,8,12. However, although some properties of these electrolyte systems such as oxidation potential and their tendency to form SEI have been measured experimentally, unfortunately no information about their transport properties have been provided. Therefore we studied these electrolyte systems by employing molecular dynamics (MD) simulation to evaluate the transport properties and then by inserting them, as the needed parameters, into the continuity equations, we calculated the voltage of LIBs for the specific current densities and in this way we examined the performance of LIBs to see which electrolyte system represent higher efficiency. Molecular modeling techniques have been widely employed to understand the principles of electrochemical and transport phenomena in LIBs. Quantum mechanical (QM) calculations, generally based on the density functional theory (DFT) 13, are capable to unravel the basic details of the interfacial structure and dynamics of electrolyte used in LIBs, such as lithium ion diffusion, solvent reorganization and decomposition reactions and also diffusion of ions and their insertion into the electrodes. Although QM calculations can effectively describe the behavior of complex electrochemical systems, their application is restricted because of high computational cost. On the other hand, MD simulations is a powerful tools for analysis of the microscopic properties of the systems containing thousands of atoms in the nanosecond time range. A number of MD and QM studies have focused on understanding the structural properties of electrolytes14–17, electrodes 18,19 and their interfacial behavior20–24. However, so far, few studies have been performed to obtain transport properties which are required as the basic parameters inherent in continuity equations, to evaluate lithium concentration and electrical potential distributions in the electrolyte, especially in close vicinity of electrodes20,25,26. 4 ACS Paragon Plus Environment

Page 4 of 55

Page 5 of 55 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Kumar et al. 14 studied the properties of EC and LiPF6 solution to be used in LIBs by applications of MD and DFT methods and by considering B3PW91 /6-311 ++ G (d, p) level of approximation for DFT calculations. Skarmoutsos et al. 15 investigated the solvation of Li cations in pure EC and in mixture with propylene carbonate (PC) and dimethyl carbonate (DMC) by MD and QM methods, where the results indicated cluster formation in tetrahedral arrangement of Li+ in pure EC and in mixture with the above mentioned solvents, which was due to strong tendency of Li+ for solvation by PC and EC molecules. Also analysis of radial distribution function (RDF) and self-diffusion coefficient evaluated by MD simulations revealed the preference of Li+ to interact with DMC molecules within its first solvation shell. You et al. 16 carried out MD simulation with OPLS-AA force field27 and provided complementary information regarding high solubility strength of organic solvents (PC and EC) for Li+ depending on dielectric constant of solvent where the obtained results indicated capability of OPLS-AA force field to predict the properties of the studied electrolyte systems. Boyer et al. 20 employed MD simulation to explore the structure of electrical double layer (EDL) of organic electrolyte (1.0 M LiPF6 in EC/DMC=50/50 (mole/mole) )in vicinity of graphite anode and it was found that when graphite was negatively charged, EC due to its cyclic molecular structure had a high population in vicinity of graphite electrode compared with that of DMC. This observation was used to explain the mechanism of decomposition and diffusion in early stages of SEI layer formation where the calculated differential capacitance of EDL (CEDL= 6.8 with the value of 4-5

2

2)

was consistent

as reported for similar electrolyte on graphite basal

electrode surface 21. By using reactive forcefield (ReaxFF28) in a MD simulation study, Yon et al. 23 predicted the formation of various compounds on SEI layer such as gases (C2H4, CO, CO2, CH4 and C2H6), inorganic compounds (Li2CO3, Li2O and LiF) and organic compounds (ROLi and ROCO2Li: R = -CH3 or -C2H5) 5 ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

on Si anode. ReaxFF has the advantage that its parameters can be adjusted for development of new enabler additives in mixture with various solvents for application to the LIBs and also for determination of electrolyte susceptibility to form a passivation layer on the anode electrode. The results of DFT calculations compared with experiments were used to explain the discrepancy of EC and PC behavior in the formation of the SEI layer, where the partial desolvation of the solvated Li+ on graphite active sites as a critical step, as well as occurrence of competitive solvation of Li+ by anion and solvent molecules determined the tendency of electrolyte to form a SEI layer. These results were consistent with the interpretations based on QM molecular orbital theory24. On the other hand, electrochemical model based on mass and charge conservation can provide information about the Li concentration and potential distribution across the cell dimension 29,30. In this case solving of corresponding continuity equations requires transport parameters depending on molecular structure of electrolyte and electrodes. These crucial parameters are usually measured either by carrying out experiments or calculated by MD simulation. Ravikumar et al. 25 determined the ionic conductivity of [EC+LiPF6] electrolyte for a range of LiPF6 concentrations (from 0.05 M to 4 M). Yoo et al.26 used MD simulation and DFT method to calculate the key parameters of the lithium-air battery continuity model. These parameters included electron transfer reaction rate constant, species diffusion coefficient and oxygen solubility and it was shown that, the battery performance was improved by increasing operational temperature, due to change in solubility of O2 in the electrolyte. Despite the number of reported works on calculating the structural and dynamic properties of electrolytes, as briefly reviewed above, there is still a severe need for an extensive and thorough QM and MD calculations to elucidate the transport properties of EC-free electrolytes 6 ACS Paragon Plus Environment

Page 6 of 55

Page 7 of 55 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

composing of (solvent + additive) mixture. The values of these properties can be used as the needed parameters to solve the continuity equations derived on the basis of Newman’s electrochemical model30 for LIBs operation. In this respect solving the pseudo-two dimensional (P2D) electrochemical model leads to obtain the profile of voltage versus current which can be used to understand the performance of LIBs. Considering an electrolyte system consisting of [(solvent + additive) + LiPF6] the needed parameters for solving continuity equations include the diffusion coefficient of Li+ and PF6- to pass through solvents, their transference number as well as ionic conductivity of electrolyte. The solvent mixtures investigated in this work were: (EMC+EC 30), (EMC+FEC 10) and (EMC+SN 10), where the numbers in parenthesis represent the weight percent of additives. The above mentioned parameters are calculated by MD simulation. The solvent and additive molecular structure and their corresponding atomic partial charges, as required in MD simulations, are determined by QM calculation31. The obtained MD results along with those obtained by solving continuity equations will be interpreted and discussed to elucidate the performance of the studied LIBs. In what follows the essential theoretical aspects of continuity equations in Newman’s electrochemical model as well as MD simulation are described.

2. Results and discussion 2.1 Molecular dynamics (MD) simulation protocol (System preparation) MD simulation of electrolytes consisting of [(solvents + additive) + LiPF6] were performed by employing GROMACS simulation package (version 5.0.7)32. In the course of simulation equilibration was achieved by minimization of energy in 20,000 steps followed by employing NVT ensemble simulation run for 1 ns and the NPT ensemble run for 20 ns to eliminate the close contacts between atoms and 7 ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

attain the equilibrated density for simulated electrolytes. Then, production run in NPT ensemble were performed for 10 ns. Nose-Hoover thermostat33,34 and Parrinello-Rahman barostat35 were used to maintain the simulation temperature and pressure with damping relaxation times of 1 and 2 ps, respectively. A Leap-frog integrator with a time step of 1 fs was used to integrate the Newton’s equations of motion. The cutoff radius of 1.2 nm was applied to vdW and particle mesh Evald (PME) summation36 for treating the electrostatics interactions. Dynamic loadbalancing (DLB)32 was used to overcome the load imbalance on the CPU and minimize the execution time. Snapshots of trajectories were stored every 20 fs to compute the time-averaged results for structural and transport properties. The detailed information of the simulation procedure is presented in Table S 8 of the Supporting Information. The current simulations were run on the HPE ProLiant DL380 Gen9 with Intel Xeon 2630 V3 and 32 GB memory on the 64-bit Ubuntu 16.04 LTS. Around 60 h were needed for a three-stage simulation including energy minimization, equilibrium attainment and property production. 2.1.1 Molecular configuration All quantum mechanical (QM) calculations to optimize the initial molecular configurations of solvent, additives and ions in addition to their respective atomic partial charges were performed by Gaussian 03 software 31. Geometry optimizations and potential energy calculations were performed for all of the molecules using density functional theory (DFT) method at B3LYP leves with 631++G(d,p) as the basis set. The optimized geometries were used as the inputs for further calculations based on the Gasussian-3 (G3) theory 37. The G3 total energy is given as37

8 ACS Paragon Plus Environment

Page 8 of 55

Page 9 of 55 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

E 0 G3

E QCISD T / 6 31G d E

G 3L

E SO

E plus

E HLC

(1)

E 2df , p

E ZPE

where E (QCISD (T)/6-31G (d)) is the correlated energy at the MP2(Full)/631G(d) optimized geometry , E (plus) is the effect of diffuse functions, E (2df, p) is the effect of polarization functions, E =3.6 is a correction for larger basis set effects, E (SO) is the spin-orbit correction, E (HLC) is a higher level correction, E (ZPE) is the zero point energy. The optimized configurations of solvent, additives and ions are shown in Figure 1. The results of QM calculations based on G3 level of theory are tabulated in Table S1 – S7 of the Supporting Information

(a)

(b)

(c)

9 ACS Paragon Plus Environment

(d)

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(e)

Page 10 of 55

(f)

Figure 1. Molecular snapshots of the ions, solvent and additive molecules: Ions: (a) Li+, (b) PF6-, and organic additive molecules: (c) EC, (d) FEC, (e) SN and solvent molecule: (f) EMC. The atoms are represented in spheres in colors as: carbon (gray), hydrogen (white), oxygen (red), nitrogen (blue), fluorine (light blue), lithium (pink) and phosphorus (yellow).

2.1.2 Forcefield parameterization and validation One of the essential requirements to obtain accurate results in MD simulation is utilizing an effective force fields. A number of force fields have been proposed in this respect which can be categorized as classical38 such as OPLS-AA force field27 and reactive force fields(ReaxFF)28. OPLS-AA force field27, which is commonly known as a suitable force field (FF) for MD simulation of electrolytes16,20,39–41 was employed in this work. This FF consists of six terms and is in the following form: U

1

" 2k

b ,i

ri

r0,i

1

" 2k

2

bonds

5

Ci 1 cos n i improper dihedral 2 N

""

i 1 j i 1

0,i

i

angles

"

N 1

(2)

2 ,i

qi q j e 2 4

r

0 ij

" "C

0,i

ij

cos 180

n i

proper dihedral n 0 12

4

n ,i

6

ij

ij

rij

rij

!

where the terms in this equation respectively represent: bond stretching, angle bending, improper dihedral torsions, proper dihedral torsions and the last two terms 10 ACS Paragon Plus Environment

Page 11 of 55 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

are for non-bonded interactions expressed as Columbic and Lenard-Jones. The Lenard-Jones parameters

ij

and

ij

are evaluated by Lorentz-Berthelot mixing38

rules for dissimilar atoms. The FF parameters for organic solvent (EMC) and additives( EC, SN, FEC) were taken from the LigParGen web server42and for PF6ion from the Reference14 and for Li+ parameters from OPLS database27. These parameters are tabulated in Table S1 – S7 of the Supporting Information. Before proceeding with the task of evaluating the accuracy of the employed OPLS-AA forcefield, the atomic partial charges of organic compounds used in this study should be rescaled to take the polarization into account. The rescaling was done by calculating the physical properties (density, viscosity, dielectric constant and dipole moment) of pure organic solvent: EMC, and additives: EC, FEC, SN were predicted by MD simulation and then compared with the reported available experimental data and theoretically calculated values according to the proposed model. The rescaling of atomic charges of electrolytes are necessary since the conventional forcefield (FF) OPLS-AA which is used in this work lacks the polar terms43–45. Therefore to compensate for this effect and to achieve closer agreement with experimental data for evaluating the physical properties by MD simulation, we found through application of a trial and error approach and by using a computational procedure to refine the results that; the atomic partial charges should be rescaled by an average scale factor of 0.6. The results are reported in Table 1. Table 1. Density, viscosity, dielectric constant and dipole moment of: EC, EEC, FMC and SN as evaluated by MD simulations and in comparison with the reported experimental data16,46–48.

EC

77

1.297

Simulation Viscosity Dielectric (Pa.s) constant ×103 1.43 80.00

EMC FEC SN

25 25 60

1.038 1.473 0.942

0.51 2.36 1.85

Solvent

Temperature (N')

Density (g/cm3)

2.54 62.8 43.22

Dipole moment

Density (g/cm3)

3.44

1.275

0.63 2.07 3.10

1.006 1.209 0.987

11 ACS Paragon Plus Environment

Experimental Viscosity Dielectric (Pa.s)× constant 103 1.9 84 (60 N' 0.65 2.96 4.1 65 2.4 56.09 (30 N'

Dipole moment

Ref

4.61

16,46

0.89 2.5 3.7

46 47 48

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

2.1.3 Molecular dynamics simulation results 2.1.3.1 Structural properties (i) Density The density of each electrolyte solution as calculated by MD simulations is given in Figure 2. For further verification of the forcefield, as rescaled by the partial charges, for application to the electrolyte mixture [(solvent + additive) + lithium salt], additional testing of physical properties were conducted by calculating density of the electrolytes as a function of salt concentration and the obtained results were compared with the reported theoretical and experimental studies. Ma et al. 8 and Logan et al.49 investigated the physical properties of (EMC+EC 30) and (EMC+FEC 10) as function of salt (LiPF6) concentration at different temperatures by proposing a model introduced as Advanced Electrolyte Model (AEM)50. AEM has been developed by using statistical-mechanics concepts based on molecular scale interactions which can accurately predicts the transport properties of various electrolyte systems. As shown in Figure 2 (a) and (b), the density results as obtained through MD simulations in this work are in good agreement with the calculated AEM values results with an average deviation of less than 1 %. The initial configuration of electrolytes constructed for density of 0.9 g/cm3 in a cubic simulation box by applying Periodic boundary conditions (PBC) in three directions and for varying salt molality (m) are represented in Table 2.

12 ACS Paragon Plus Environment

Page 12 of 55

Page 13 of 55

1.3 °

MD: 10 C MD: 20°C MD: 30°C °

AEM: 10 C °

AEM: 20 C

1.2

(g/cm 3)

AEM: 30 °C

1.1

1 0

0.5

1

1.5

2

m(mol/kg)

(a) 1.3 MD: 10°C MD: 20°C MD: 30°C AEM: 10 °C AEM: 20 °C

1.2

(g/cm 3)

AEM: 30 °C

1.1

1 0

0.5

1

1.5

2

m(mol/kg)

(b) 1.2 MD: 10°C °

MD: 20 C MD: 30°C

1.1

(g/cm 3)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

1

0.9 0

0.5

1

1.5

m(mol/kg)

13 ACS Paragon Plus Environment

2

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 14 of 55

(c) Figure 2. Density # g / cm 3 versus LiPF6 concentration (mol/kg) for the electrolyte systems of [(solvent + additive) +LiPF6]. The (solvent+additive) mixtures are: (a) (EMC+EC 30), (b) (EMC+FEC 10) (c) (EMC+SN 10) at Temperatures 10 N') 20 N' and 30 N'. Solid lines represent the evaluated densities calculated by AEM model49.

Table 2. MD simulation results of investigated electrolyte [(solvent + additive) + LiPF6] systems. Nm is the number of molecules in (solvent + additive) mixture.

(solvent +additive)

(EMC+EC 30) Nm = (250:130)

(EMC+FEC 10) Nm = (250:27)

(EMC+SN 10) Nm = (250:37)

Molality (mol /kg)

Initial Box Dimensions (nm3)

0.03 0.50 1.00 1.50 2.00 0.03 0.50 1.00 1.50 2.00 0.03 0.50 1.00 1.50 2.00

4.10×4.10×4.10 4.20×4.20×4.20 4.30×4.30×4.30 4.39×4.39×4.39 4.48×4.48×4.48 3.72×3.72×3.72 3.81×3.81×3.81 3.89×3.89×3.89 3.98×3.98×3.98 4.05×4.05×4.05 3.77×3.77×3.77 3.85×3.85×3.85 3.94×3.94×3.94 4.02×4.02×4.02 4.11×4.11×4.11

Equilibrium Volume (nm3) 10 N' 20 N' 30 N' 58.12 58.82 59.51 59.65 60.37 61.10 61.35 61.97 62.74 62.87 63.62 64.28 64.55 65.18 65.88 46.38 46.96 47.88 47.62 48.13 48.74 48.81 49.39 49.96 50.09 50.61 51.33 51.32 51.97 52.42 48.70 49.32 49.97 49.76 50.42 51.10 50.99 51.69 52.33 52.27 52.89 53.57 53.51 54.13 54.83

10 N' 0.029 0.529 1.029 1.479 1.929 0.036 0.663 1.293 1.857 2.427 0.034 0.501 0.944 1.398 1.800

Molarity (mol/L) 20 N' 0.028 0.523 1.018 1.462 1.911 0.037 0.656 1.278 1.838 2.397 0.034 0.494 0.932 1.381 1.779

30 N' 0.028 0.516 1.006 1.447 1.890 0.036 0.647 1.263 1.812 2.376 0.033 0.487 0.920 1.364 1.756

(ii) Radial distribution functions Additionally, the structural distributions of electrolyte systems were studied by radial distribution functions (RDFs). The RDF, g i

j

r , of atoms j in a molecule

in respect to atom i in another molecule chosen as the reference is defined as38:

14 ACS Paragon Plus Environment

Number of LiPF6 1 19 38 57 76 1 15 29 44 58 1 15 29 44 58

Page 15 of 55 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

gi

j

r

n r # 4 r 2dr 1

(3)

where, # is bulk number density and n r is number of j atoms in a molecule at a radial distance r from atom i in reference molecule. Figure 3 (a, b, c, d) shows the RDF for atom pairs i and j, where i is Li+ and j in Figure 3 (a, b, c, d) are: (a) carbonyl oxygen of EC, (b) carbonyl oxygen of EMC, (c) P of PF6- and (d) F of PF6- pair in the electrolyte systems [(EMC+EC 30) + LiPF6] as function of concentration at 20 N'3 The RDFs for Li+-carbonyl oxygen of EC and EMC ( Figure 3 (a) and (b)) exhibit a sharp peaks at the range of 0.20-0.23 nm which indicates the first solvation shell as reported in literature14,20,25 for similar electrolyte systems. Comparatively smaller peaks as observed at r

0.8 nm

r

0.4 nm

and

for carbonyl oxygen represent next solvation shells. Also as expected,

at large values of r the RDFs approach unity and due to increase in salt concentration which causes a decrease in the height of the sharp peaks as well as a shift in their position, as depicted in Figure 3 (a) and (b). Association of cations and anions can be understood by inspecting the RDF of Li+-PF6- pair as it is seen in Figure 3 (c) and Figure 3 (d) where by increase in the salt concentration the height of the first peak increases. This observation which indicates that at low concentrations most of the salt exists in the dissociated form and, as discussed in Section 3.2, the amount of salt in the associated form increases with increase in concentration. The RDFs for electrolyte systems [(EMC+FEC 10) + LiPF6] and [(EMC+SN 10) + LiPF6] are presented in Figure S 2 and Figure S 3 of the Supporting Information.

15 ACS Paragon Plus Environment

The Journal of Physical Chemistry

3

25

0.03 m 0.50 m 1.00 m 1.50 m 2.00 m

20

2.5

15

gLi-OEC(r)

10

2

0.03 m 0.50 m 1.00 m 1.50 m 2.00 m

15 10

5

1.5 0.2

30 25 20

gLi-OEMC(r)

3

0.3

2 5

0.2

0.3

1

1

0.5

0.0 0.0

0.5

0 0.0

1.0

0.5

1.0

r(nm)

r(nm) (a)

(b)

15

50

0.03 m 0.50 m 1.00 m 1.50 m 2.00 m

14

0.03 m 0.50 m 1.00 m 1.50 m 2.00 m

40 30

gLi-F(r)

13

gLi-P(r)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 16 of 55

5 4

10

3

5

2 1 0 0.0

0.5

0 0.0

1.0

0.5

r(nm)

1.0

r(nm) (d)

(c)

Figure 3. The radial distribution functions for i and j atoms (as defined in Equation (3)). i is Li+ and j atoms are: (a) carbonyl oxygen of EC, (b) carbonyl oxygen of EMC, (c) P of PF6- and (d) F of PF6- in electrolyte systems [(EMC+EC 30) + LiPF6] as a function of salt concentration at 20 N'3 The RDF curves in Figures (a) and (b) for higher values of

and

are respectively shown in the boxes.

2.1.3.2 Transport properties The transport properties which are evaluated by MD simulation and for solving P2D electrochemical model30 of LIB are as follows: (1) diffusion coefficient D ,which represents mass transport of species in an electrolyte mixture, (2) ionic 16 ACS Paragon Plus Environment

Page 17 of 55 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

conductivity $ of electrolyte, which represents the amount of charge carried by the ions and (3) transference number t Li

which shows the contribution of ions to

the total ionic conductivity. In what follows we present summary of the above mentioned properties. (i) Diffusion coefficients The diffusion coefficient of molecule i D i is related to the mean square displacement (MSD) of individual molecules as described by the Einstein equation38: Di

1 % &' 6N % i lim

Ni

"

R j ,i %

t

R j ,i t

(4)

2

j 1

where N i is the numbers of molecule of species i and R j ,i is position of

j

th

molecule of species i . The angle brackets denote an ensemble average computed over

t ( [0,% ] .

Due to the long-range interactions the diffusion coefficients as

calculated by Equation (4) is strongly dependent on size of simulation system 51,52. Yeh and Hummer52 developed the analytical correction to reproduce quantitatively the diffusion coefficient based on hydrodynamic theory for spherical particles as follow: D i'

D iMD

2.837 k BT 6 )' L

(5)

where D iMD is finite system size diffusion coefficients obtained by using Equation(4) and the second term is known as Yeh-Hummer correction where the other symbols in this equation are: k B Boltzmann’s constant, T absolute temperature, L length of the simulation box and )' shear viscosity of the system at defined temperature T. In addition, it should be noted the Yeh-Hummer correction does not explicitly 17 ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 18 of 55

depend on the size of molecules in mixture. However, all species of a multicomponent system experience identical finite-size effects52,53. In MD simulations, the required shear viscosity is computed from transverse current correlation function54, through the transverse momentum fields 32. In this method, 16 transverse-current autocorrelation functions corresponding to different k-vectros are considered, resulting in 16 values of ) k . The shear viscosity )' is calculated by fitting values of ) according to Equation (6). ) k

)' 1 a k 2

(6)

where a is fitting parameter. Yeh and Hummer52 proved that the shear viscosity is independent of the system size. However, the evaluated viscosity by Equation (5) are almost constant as shown in Table S 9 of the Supporting Information. The evaluated diffusion coefficients ( D i' ) along with their average deviation (± 0.1) or their statistical confidence for infinite system size are reported in Table 3. These values have been obtained by applying Yeh-Hummer correction and by using the viscosity data as represented in Figure 4. Figure 5 illustrates, as an example, the plot of MSD for Li+ and PF6- versus time for various concentrations and for [(EMC+EC 30) +LiPF6] at 10 N'3 The finite system size diffusion coefficients ( D iMD ) has been determined by evaluating the slope of this plot (Figure 5), according to Einstein regime. Inspecting the D iMD values indicate an initial increase with salt concentration, and then they decrease at higher concentrations. However by performing Yeh-Hummer correction, the obtained D i' values (as it is seen in Table 3) show a decreasing trend with increasing salt concentration. That is, as expected, they have their maximum value for diffusion coefficient of ions at infinite salt dilution. However our search of literature indicated that the only 18 ACS Paragon Plus Environment

Page 19 of 55 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

reported experimental data for diffusion coefficient of Li+ is for 1 m (1 mol/kg) concentration of LiPF6 in electrolyte system of [(EMC + EC 30) +LiPF6] at 30 N' which is 2.66×10-10 m2/s55. As it is seen in Table 3 for this system at the same concentration and temperature, the obtained value by our MD simulation is 2.72×10-10 m2/s which is in close agreement with the experimental value. Also the viscosity values as obtained by our MD simulation, and illustrated in Figure 4, are in relatively good agreement with the experimental data for the studied electrolyte systems and also, as it is seen in Figure 4 (a) for the system of [(EMC + EC 30) + LiPF6] our viscosity results are well-fitted to the calculated values as obtained by AEM49 method. The ability of MD simulation to predict the diffusion and viscosity coefficients for electrolyte systems, as represented in this work, can be considered a significant progress compared to the previous MD simulations of electrolyte system, where polarization effects were ignored in application of the forcefield. Therefore it is expected that the high accuracy of these coefficients would lead to obtaining more accurate values for ionic conductivity compared with previously reported results14,20,25. This point will be discussed further in next section. Table 3. Computed diffusion coefficients D i' and their statistical confidence (±0.1) of ions (Li+ and PF6-) in the studied electrolyte systems [solvent + additive) +LiPF6] evaluated by MD simulations for molality concentration (m) of salt LiPF6 at different temperatures.

(solvent+ additive) Molality (m)

(EMC+EC 30)

(EMC+FEC 10)

0.03 0.50 1.00 1.50 2.00 0.03 0.50 1.00 1.50

D Li' * 0.1 10 oC 3.10 2.84 1.64 1.12 0.85 3.49 2.91 2.13 1.81

m2 + 1010 s

20 oC 3.23 3.05 2.21 1.40 1.12 3.70 3.10 2.45 2.03

30 oC 3.47 3.23 2.72 1.77 1.27 3.85 3.29 2.78 2.17

19 ACS Paragon Plus Environment

' * 0.1 D PF 6

10 oC 3.21 2.84 1.64 1.26 0.86 3.55 2.92 2.14 1.80

m2 + 1010 s

20 oC 3.24 3.03 2.21 1.40 1.24 3.75 3.13 2.47 2.03

30 oC 3.51 3.21 2.66 1.74 1.27 3.93 3.34 2.81 2.17

The Journal of Physical Chemistry

2.00 0.03 0.50 1.00 1.50 2.00

(EMC+SN 10)

1.44 3.64 3.01 2.20 1.76 1.41

1.70 3.72 3.29 2.54 2.10 1.65

2.10 3.97 3.70 3.21 2.43 1.97

1.45 3.47 3.01 2.00 1.75 1.42

15 °

MD: 10 C °

MD: 20 C MD: 30°C

Viscosity(Pa.s x103 )

Exp: 10°C °

Exp: 20 C

10

°

Exp: 30 C AEM: 10 °C AEM: 20 °C AEM: 30 °C

5

0 0

0.2

0.4

0.6

0.8

1

1.2

1.4

1.6

1.8

2

1.2

1.4

1.6

1.8

2

m(mol/kg)

(a) 5 °

MD: 10 C

4.5

°

MD: 20 C °

MD: 30 C

4

Viscosity(Pa.s x103 )

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 20 of 55

3.5 3 2.5 2 1.5 1 0.5 0 0

0.2

0.4

0.6

0.8

1

m(mol/kg)

(b)

20 ACS Paragon Plus Environment

1.73 3.96 3.33 2.57 2.11 1.70

2.11 4.03 3.67 3.22 2.47 1.98

Page 21 of 55

5 °

MD: 10 C

4.5

°

MD: 20 C °

MD: 30 C

Viscosity(Pa.s x103 )

4 3.5 3 2.5 2 1.5 1 0.5 0 0

0.2

0.4

0.6

0.8

1

1.2

1.4

1.6

1.8

2

m(mol/kg)

(c) Figure 4. Viscosity coefficient (Pa.s ×103) as a function of LiPF6 concentration for electrolyte systems: (a) [(EMC+EC 30) + LiPF6] (b) [(EMC+FEC 10) + LiPF6] (c) [(EMC+SN 10) + LiPF6] at temperatures 10 oC (green), 20 oC (red) and 30 oC (blue). Solid lines in Figure (a) are AEM results49, triangle symbols are experimental49and circles are MD simulation results.

12 0.03 m 0.50 m 1.00 m 1.50 m 2.00 m

10

8

MSD (nm2 )

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

6

4

2

0 0

1

2

3

4

5

6

7

(ns)

(a)

21 ACS Paragon Plus Environment

8

9

10

The Journal of Physical Chemistry

12 0.03 m 0.50 m 1.00 m 1.50 m 2.00 m

10

8

MSD (nm 2 )

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 22 of 55

6

4

2

0 0

1

2

3

4

5

6

7

8

9

10

(ns)

(b)

Figure 5. MSD vs time for ions in for electrolyte system [(EMC+EC 30) + LiPF6] with various molality concentration of LiPF6 at

"

5N' ions are: (a) Li+ (b) PF6-.

The decreasing trend of diffusion coefficients with concentration of lithium salt concentration, as it is seen in Table 3, has been attributed to the competitive Li+ solvent and Li+

PF6- interactions, as depicted in Figure 3, which lead to

formation of three different types of salt solvated complexes, namely, solventseparated ion pairs (SSIPs), contact ion pairs (CIPs) and aggregates (AGGs)56. Formation of these complexes has been observed by experimental studies in the system of [propylene carbonate (PC) + LiBF4]57 where in dilute electrolyte, there is a significant number of free solvent molecules, but by adding more Lithium salt, the solvated complexes (CIPs and AGGS) tend to form. These complexes have less mobility and carry less charge in comparison to SSIPs. Also the diffusion coefficients of Li+ and PF6- in the electrolyte systems are almost the same (as it is seen in Table 3). This can be ascribed to the fact that organic polar solvents have higher tendency to solvate the Li+ rather than PF6-. However due to higher charge density of Li+ compared to the dispersed charge of PF6-, the transport of Li+ 22 ACS Paragon Plus Environment

Page 23 of 55 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

becomes dominant and the change in the diffusion coefficient, in this case as resulted from solvent interference, is attributed to what is known as vehicular diffusion mechanism5,25,56. The dependence of diffusion coefficients on solvent composition is clearly seen in Table 3 and also it is observed that Li+ and PF6- diffusion coefficients increase when additive EC is replaced by additives SN and FEC. This behavior is due to presence of carbonyl group in EC molecular structure, which protrudes from the ring and can solvate Li+ ions more tightly than EMC. This causes Li+ readily, form a stable solvation SSIP complex with EC and as results retards the diffusion of Li+ would occur. This effect can be used to interpret the difference in diffusion coefficients of Li+ in the electrolyte systems of [(EMC + SN 10) + LiPF6] and [(EMC+FEC 10) + LiPF6]. Although in the former electrolyte system Li+ has a greater diffusion coefficients to Li+ at low salt LiPF6 concentration, but due to smaller number of FEC molecules (27 as it is seen in Table 2 compared with the number of SN molecule 37) in the latter electrolyte system, the diffusion coefficients of Li+ in the latter electrolyte system is greater than in the former electrolyte system at the higher concentration of LiPF6. The reason for this behavior can be explained by considering the fact that at low number of FEC molecules, as mentioned above, less SSIP complexes are formed between Li+ and FEC and then more free Li+ are present to diffuse easily through this electrolyte system. Inspecting the results in Table 3 indicates that the diffusion coefficients increase with temperature rise and this behavior is, of course, due to increase in thermal motion of ions. But this increase is almost the same for cations and anions. This behavior can be attributed to the capability of mixed solvents to separate cations and anions from each other and depends on variation of dielectric constant of the 23 ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 24 of 55

mixed solvent with temprature58,59. Although it is expected that for a binary solvent system, the dielectric constant should represent a decreasing trend with temperature, but it was found that the dielectric constant of binary solvents containing EMC, at higher weight fractions, remains almost constant with temperature variation60. For instance, the dielectric constant of binary solvent (EMC+FEC 10) indicates a small change from 11.10 to 10.64 respectively, for temperature variation of 20 N' and 30 N'3 The above mentioned values for dielectric constant for this binary mixture were calculated by using a mixing rule as suggested by Oster61 and from the reported dielectric values of pure EMC and FEC6. Therefore, according to these results, temperature variation has a negligible effect on dielectric of binary solvents and thereby on their capability to dissociate the salt into ions and also on diffusion coefficient of ions. (ii) Ionic conductivity Ionic conductivity is an important electrolyte property which quantifies the ability of the electrolyte to conduct the ions. It can be calculated from the collective mean square displacement of ions in the solution based on Einstein relation38 as represented in the following form: $ %

e2 % &' 6%V k T B lim

"" i

z i z j Ri % t

Ri t

Rj % t

Rj t

(7)

j

where $ % denotes ionic conductivity, z i is charge of species i , e is charge of electron, V is simulation box volume. kB, T and Ri are respectively Boltzmann constant, Kelvin temperature and position of species i and j at time t and % t .The ionic conductivity can be estimated by calculating the slope of the collective MSDs versus time % for % & ' . However, this straightforward method is not suitable due to computational burden and noisy nature of collective MSDs. 24 ACS Paragon Plus Environment

Page 25 of 55 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Therefore, an alternative procedure as proposed by Borodin et al. 62 and used by other researchers 63,64 as well can be applied which remarkably can reduces the complexity of the problem and computational burden. In this procedure the ionic conductivity can be expressed as summation of two terms62:

$ $ NE $cc

(8)

where, $ NE and $cc are defined as: $ NE %

$cc %

e2 lim % &' 6%V k T B

"

2e 2 lim % &' 6%V k T B

""

z i2 R i % t

Ri t

(9)

2

i

i

z i z j Ri % t

Ri t

Rj % t

Rj t

(10)

i,j

The first term represents Nernst-Einstein conductivity, which is based on ionic mobility and represent the ideal ionic conductivity of an electrolyte in the absence of ions dynamic correlations. $ NE can be directly calculated from diffusion coefficients and corresponds to ionic conductivity of dilute aqueous KCl solution, as can be obtained from Walden plot12. By combining Equations (4) and (9), the Nernst-Einstein conductivity for an electrolyte consisting of cations (c) and anions (a) can be expressed as $ NE

e2 N a z a2 D a V k BT

(11)

N c z c2 Dc

The second term, $cc , expresses the cross-correlated conductivity contribution between the different ions for

i - j

. Because of correlated motion and association

between oppositely charge ions, this term is negative64. Therefore, an . % ratio defined as the degree of uncorrelated motion or degree of dissociation, can be proposed as62:

25 ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

. %

"" i

z i z j Ri % t

Ri t

Rj % t

Rj t

Page 26 of 55

(12)

j

"

z i2 R i % t

Ri t

2

i

where for . % . %

1.0 the ion motion is completely uncorrelated, whereas for

0.0 cations and anions move together as ion pairs and then give rise to zero

conductivity. In this method, the key assumption is that; the value of . % converges in short times whereas the dynamics may still have some subdiffusive contributions64. To validate this assumption, we calculated . % by using 5 % of our trajectories at different time periods (See Figure S 4 in the Supporting Information). The ionic conductivity values and lithium salt (LiPF6) concentrations along with the reported experimental data 49 are shown in Figure 6 and also are tabulated in Table S 10 in the Supporting Information . The mean square deviations (MSDs) were calculated for MD and AEM results of electrolyte [(EMC+EC 30) + LiPF6] according to Equation (13). MSD

1 " $ $ exp n

(13)

2

Comparison of results in Table S 10 in the Supporting Information shows an excellent agreement between the experimental49 and the ionic conductivity obtained through MD simulation with MSD of 0.27 which is close to MSD of wellaccepted advance electrolyte method (AEM). These results as calculated through the degree of uncorrelated motion ( . % ) are illustrated in Table 4. In this table, it is seen that by increasing Li salt concentration, the ionic conductivity initially increases, which is due to tendency of Li+ to form SSIP complexes, and then decreases due to formation of CIPs and AGGs species. This behavior which is

26 ACS Paragon Plus Environment

Page 27 of 55

confirmed by decrease in . % , as it is seen in Table 4, and is in excellent match with the results obtained using pulsed field gradient nuclear magnetic resonance (PFG-NMR) measurement 65. At infinite dilution of lithium salt, as it is seen in Table 4, the . % values are close to unity which indicates complete dissociation of lithium salt and uncorrelated motion of Li+ and PF6- through binary solvent. However the increase in viscosity (see Figure 4) at higher concentrations of lithium salt is another reason for decrease in ionic conductivity and diffusion coefficient. 15 MD: 10°C

Ionic conductivity(mS/cm)

MD: 20°C MD: 30°C Exp: 10°C Exp: 20°C

10

Exp: 30°C AEM: 10 °C AEM: 20 °C AEM: 30 °C

5

0 0

0.2

0.4

0.6

0.8

1

1.2

1.4

1.6

1.8

2

1.4

1.6

1.8

2

m(mol/kg)

(a) 15 MD: 10°C °

MD: 20 C

Ionic conductivity(mS/cm)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

MD: 30°C

10

5

0 0

0.2

0.4

0.6

0.8

1

1.2

m(mol/kg)

27 ACS Paragon Plus Environment

The Journal of Physical Chemistry

(b) 15 MD: 10°C MD: 20°C

Ionic conductivity(mS/cm)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 28 of 55

MD: 30°C

10

5

0 0

0.2

0.4

0.6

0.8

1

1.2

1.4

1.6

1.8

2

m(mol/kg)

(c) Figure 6. Ionic conductivity versus LiPF6 concentration for the studied electrolytes systems [(solvent + additive) + LiPF6]. (Solvent + additive) are: (a) (EMC+EC 30) (b) (EMC+FEC 10) (c) (EMC+SN 10) at temperatures 10 oC (green), 20 oC (red) and 30 oC (blue). Solid lines in Figure (a) are AEM results49, triangles represents experimental data49and circles are the MD simulation results. Table 4. The degree of uncorrelated motion .

%

for the studied

electrolyte systems, [(solvent + additive) + LiPF6] (solvent + additive)

(EMC+EC 30)

(EMC+FEC 10)

(EMC+SN 10)

. %

Molality (m) 10 N' 0.91 0.43 0.55 0.42 0.37 0.84 0.52 0.59 0.43 0.35 0.81 0.55 0.63 0.47 0.39

0.03 0.50 1.00 1.50 2.00 0.03 0.50 1.00 1.50 2.00 0.03 0.50 1.00 1.50 2.00

28 ACS Paragon Plus Environment

20 N' 0.94 0.54 0.54 0.46 0.44 0.82 0.54 0.56 0.44 0.36 0.82 0.54 0.67 0.44 0.41

30 N' 0.96 0.56 0.50 0.44 0.46 0.81 0.56 0.54 0.44 0.37 0.82 0.55 0.48 0.45 0.37

Page 29 of 55 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Although the studied electrolytes containing SN and FEC have lower . % than the one containing EC, their conductivity is higher (See Figure 6), which is due to high diffusivity of ions in these systems, as reported in section 3.2.1. Lower . % values of these two electrolytes arise from low dielectric constant of solvent mixture and when these additives are used in binary mixture with the solvent (EMC), the mixture would possess lower polarity and as a result weaker capability to separate or dissociate lithium salt into ions (see Table 1), as has been pointed out earlier in Section 3.2.1. As it is seen in Figure 6, raising temperature causes an increase in ionic conductivity of electrolytes. However, . % remains approximately constant (see Table 4) and does not show any significant trend of variation with temperature. These results indicate that the temperature variation does not affect ion pairing in the studied electrolyte systems which is due to almost negligible variations of dielectric constants with temprature60. In other words, since the ionic conductivity depends on ionic diffusivity, it is expected that, it should indicate an increasing trend similar to diffusivity as a result of temperature rise. But as stated earlier ionic conductivity consists of two contributions arising from correlated and uncorrelated motions of ions, therefore, temperature rise may have a reverse effect on these contributions, that is, increasing one causes a decrease in the other. As a result . % would remain constant with temperature variation.

Walden plot is a primary method based on ion-pairing concept12 and is used to assess the dissociative capability of a solvent to separate or dissociate a dissolved species into ions. This plot represents the variation of molar conductance / S cm 2 / mol

$ mS cm 1 +10

3

c mol cm

3

of aqueous solution of KCl, as an

ideally strong electrolyte which completely dissociates into ions, versus inverse 29 ACS Paragon Plus Environment

The Journal of Physical Chemistry

viscosity coefficient of solution. The other electrolyte solutions have a deviation from KCl and as it is shown in Figure 7, all the / values lie in an area between ideal KCl line (blue solid line) and the line of 10 % ionization (red discrete line). This plot can be used for verification of degree of dissociation and the deviation from Nernst-Einstein conductivity as represented in Table 4. It is interesting, that the studied electrolyte systems containing SN and FEC indicate larger deviation than the electrolyte system containing EC (See Figure 7). Therefore we may be misled by this plot to believe that the electrolyte system containing EC, as an additive, would be a better choice to be used in LIBs. But in section 3.3, by solving Newman’s electrochemical model we will show that the electrolyte containing FEC and SN would represent more effective performance if used in LIBs. Therefore it may be stated that the Walden plot loses its credibility for mixed solvent systems and perhaps a criterion based on solvent-co-solvent rule correlation, as reported in literatures12,66 would be more effective in this respect. 103 Ideal KCl Line 10 % Ionisation Line

102

(S cm 2 /mol)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 30 of 55

101

100

10-1

10-2 10-2

10-1

100

101 -1

(Poise-1 )

30 ACS Paragon Plus Environment

102

103

Page 31 of 55 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Figure 7. Walden plot12 for the studied electrolyte systems [(solvent + additive )+LiPF6] where (solvent + additive ) are represented by symbols: (EMC+EC 30) circles, (EMC+FEC 10) stars and (EMC+SN 10) triangles and temperatures are represented by colors: 10 N' blue , 20 N' green and 30 N' red.

(iii) Transference number The transference number, t i , represents the fraction of the current, I i , that is carried by an ion i in a solution2 and is defined as: ti

Ii I0

(14)

where I 0 stands for the total current. In a dilute solution, the transference number can be evaluated from the diffusion data as obtained by PFG-NMR measurements through following Equation29: ti ,p

zi

"

zj

2

(15)

ci Di 2

cjDj

j

where, c j , D i and z i are molar concentration, diffusion coefficient and ion charge for the ionic species i , respectively. This dimensionless property is also called the transport number which measures the relative contribution of species to the overall mass transport in electrolytes. Transport number given by Equation (15) can be calculated by using D i ,NMR as measured by PEG-NMR29. The transport number is expected to be about 0.5 due to equality of diffusion coefficients of Li+ and PF6ions58,59. The transference number of lithium ions t Li is usually less than its transport number, because cations (Li+) and anions (PF6-) have correlated motion. Therefore to take the correlative properties of ions pairs into account and calculate more accurate transference number, we will use the following Equation67:

31 ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

$ PF

t f ,Li

6

$ PF

PF6

6

where $ i

j

PF6

+ $ Li

+ $ Li

$ PF

Li

$ Li2

Li 6

PF6

Page 32 of 55

(16)

PF6

2 + $ Li

PF6

represent the conductivity of pair correlative ions in the electrolyte

system. The pair correlative conductivities can be evaluated through following equations: $ Li

Li

$ PF

2$ Li

lim

PF6

N

e2 % &' 6%V k T B lim

z i z j Ri % t

Ri t

Rj % t

Rj t

(17)

i 1 j 1

lim

%

N

""

e2 % &' 6%V k T B

%

PF6

6

e2 % &' 6%V k T B

%

N

N

""

z i z j Ri % t

Ri t

Rj % t

Rj t

z i z j Ri % t

Ri t

Rj % t

Rj t

(18)

i 1 j 1 N

N

""

(19)

i 1 j 1

where N and N in the above equation stand for the numbers of Li+ and PF6-, respectively. The other symbols in Equations (17) to (19) have been defined in previous section. Equation (16) can be rearranged in the following form: t f , Li

02 4 12 1

(20)

41 1 1 0

where, 0 and 1 , describe the non-ideality of solution, and are given as follows:

0 1

$ PF6

2$ PF6

Li

PF6

$ Li

$ Li $ PF

6

PF6

(21) Li

(22)

Li

$ Li

Li

Table 5 reports the transport number t p ,Li

along with the ion pair correlated

transference number t f ,Li , as respectively can be calculated by using equation (15) and (20). As it is seen in this table, since the same value of diffusion coefficients for both ions for the studied electrolyte systems has been inserted in Equation (15), 32 ACS Paragon Plus Environment

Page 33 of 55 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

the transport number of Li+ is obtained as 0.5 by MD simulation (See Table 3). However it has been reported by other researchers58,59 that t p ,Li as evaluated by MD simulation does not change with temperature and with salt concentration variations, whereas experimental measurement of diffusion coefficients by PEGNMR indicated that the evaluated transport number of lithium ion t p ,Li increases with salt concentration68,69. On the other hand, it has been shown that the transference number t f ,Li behaves differently with salt concentration change58. Inspecting t f ,Li as obtained in this work and reported in Table 5 for Li+ by MD simulation, indicates that, its variation is in the range of 0.1-0.4 which is the same range as reported by other researches2,5,69. Also the decrease of t f ,Li with salt concentration is due to reduction of degree of dissociation69. It should be noted that small tolerance in t f ,Li by concentration variation, as it is seen in Table 5, is related to the high sensitivity of t f ,Li values due to changes of 0 parameter67. Therefore the proper value for 0 according to Equation (20) as calculated for 0 2 t f ,Li 2 0.4 , should be 0.95 2 0 2 1 (See Figure S 5 in the Supporting Information).

33 ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Table 5. The transference number t f

(solvent + additive)

(EMC+EC 30)

(EMC+FEC 10)

(EMC+SN 10)

, Li

Page 34 of 55

and transport number t p , Li for the studied electrolyte systems.

10 N' 0.37 0.27 0.30 0.30 0.12 0.37 0.32 0.37 0.10 0.16 0.44 0.34 0.35 0.37 0.19

0.03 0.50 1.00 1.50 2.00 0.00 0.50 1.00 1.50 2.00 0.03 0.50 1.00 1.50 2.00

t p ,Li

t f ,Li

Molality (m)

20 N' 0.36 0.26 0.32 0.22 0.13 0.50 0.28 0.39 0.26 0.13 0.37 0.35 0.25 0.13 0.15

30 N' 0.35 0.31 0.31 0.33 0.13 0.31 0.25 0.22 0.21 0.18 0.41 0.31 0.26 0.28 0.18

10 N' 0.48 0.50 0.50 0.47 0.50 0.50 0.49 0.50 0.50 0.50 0.51 0.50 0.52 0.50 0.50

20 N' 0.50 0.50 0.50 0.50 0.47 0.50 0.50 0.50 0.50 0.50 0.50 0.50 0.50 0.50 0.50

30 N' 0.52 0.50 0.51 0.50 0.50 0.49 0.50 0.50 0.50 0.50 0.49 0.50 0.50 0.50 0.50

3.2 P2D electrochemical model In this section we will demonstrate that, how application of pseudo twodimensional (P2D) electrochemical model30 provides an effective approach to deal with operational properties of the LIBs at macro scale level. These properties which can be evaluated by considering variation of voltage and current density in charge-discharge operation, as occurs in a LIB, are generally influence by some undesirable effect such as voltage drop, polarization, ohmic resistance and hysteresis energy loss which significantly reduce the performance of LIB. To examine the efficiency of a LIB, which is intended to promote its performance, by newly designed electrodes or electrolytes, a great deal of experimental measurements are required which their provision is time consuming and expensive. Therefore it is essential to calculate the desired properties by employing a 34 ACS Paragon Plus Environment

Page 35 of 55 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

proficient model. Fortunately, in pervious sections of this manuscript we succeeded to evaluated properties of electrolyte by using DFT and MD methods at Pico-Nano scale level. These properties are needed as the basic parameters for P2D macro scale model. Therefore were are, now, in a position to employ P2D model, to make a comparison between our newly designed electrolytes, and decide application of which electrolyte would provide the LIB with more efficient performance. The electrochemical model was initially introduced by Doyle et al.30 based on porous electrode theory, concentrated solution theory and Butler-Volmer kinetics70,71. This model can be employed to study the performance of LIBs by using pseudo-two dimensional (P2D) electrochemical model as proposed by Newman30. Figure 8 shows a schematic of P2D electrochemical model for LIB operating as a discharging battery and consists of three sections: (i) anode: (negative electrode 0 2 x 2 L neg ), (ii) separator: (electrolyte, L neg 2 x 2 Lsep L neg ), (iii) cathode: (positive electrode L neg Lsep 2 x 2 Lcell ), where Lcell

35 ACS Paragon Plus Environment

L neg Lsep L pos .

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 36 of 55

Li+

Discharge

Figure 8. Schematic diagram of P2D electrochemical model for a LIB operating as discharging battery. The electrolyte is [(EMC+EC)+LiPF6], where c s x , r , t in the circles indicates the concentration of Li. The darker colors are used to show higher concentration of Li.

The following chemical reactions occur in LIB: Oxidation in anode: L x C 6 & Li x z C 6 zLi

(23)

ze

where lithiated carbon LixC6 loses z electrons and produce z Li+ (z lithium ions). Reduction in cathode: L y z MO 2 zLi

ze & Li y MO 2

(24)

where M is a transition metal. In Equation (24) lithium ion (Li+) in the electrolyte takes z electrons from cathode surface and produces z Li (z lithium atom) which would intercalate into the structure of Lithium metal oxide. Thus the overall cell reaction is: L y z MO 2 Li x C 6 & Li y MO 2 Li x z C 6

(25)

The separator consists of an electrolyte [(solvent + additive) + LiPF6] in mixture with a polymer gel such as PP/PE/PP72,73. 36 ACS Paragon Plus Environment

Page 37 of 55 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

As a result of chemical reaction when LIB is operating as a discharge battery, the concentration of Li+ in the electrolyte in the vicinity of anode increases while the amount of lithium atom (Li) in the anode decreases. In the cathode during discharge the concentration of Li+ in the electrolyte in the vicinity of the cathode decrease while the amount of Li atom in the cathode decreases. Therefore a concentration gradient is developed in the electrolyte across the LIB cell which causes diffusion of Li+ through electrolyte from the anode toward the cathode. Since Li+ diffusion flow is sufficiently fast, the concentration gradient intensifies with passage of time and then leads to concentration segregation and polarization of electrolyte across the cell where this would influence the performance of cell by controlling the rate of cell chemical reaction. When LIB cell is operating as a charging battery, the reverse of reactions (23) to (25) would occur. That is: LiyMO2, positive electrode, will act as the anode and LixC6, negative electrode, as the cathode. As the performance of LIB is controlled by transport properties of Li+ along with other transport properties such as ionic conductivity of electrolyte, it is necessary to use a cell model to study the influence of these properties. In the Supporting Information we presented the governing equations in the pseudo-two dimensional electrochemical model as proposed by Newman30. This model along with ButlerVolmer equation70,71 (Equation (S 10), in the Supporting Information) will be employed to evaluate jLi the rate of cell chemical reaction at different concentration of salt as well as the cell voltage (Equation (26)) which is considered as the most important factor and determines the efficiency of LIB in used to justify the performance of LIB when operating as discharging/charging battery29,30,73: V t

s

Lcell , t

s

0, t

Rf I app t A

(26)

37 ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

where

s

Lcell , t and

s

0, t , as expressed by Equations (S 3 and S4) are the

potentials of solid phase electrode at x

Lcell and x 0 at time t, respectively. In

Equation (26) R f is an empirical constant for contact resistance of the cell and I app t

is the applied current.

Therefore, to investigate the behavior of the studied electrolytes we need to solve simultaneously the continuity equations, Equations (S 1) to (S 12) for one dimensional LIB according to P2D electrochemical model due to Newman29 and as represented in Figure 8. The purpose of this task is to evaluate the performance of the electrolyte systems [(solvent-additive) + LiPF6] under different operating conditions for charge and discharge processes of LIB. The electrodes of this LIB model are LixC6 (anode) and LiyMn2O4 (cathode) where the oxidation/reduction reactions follow the same mechanism as explained in Section 2.1. For solving the equations of P2D electrochemical model and calculating the required properties the COMSOL Multiphysics 5.2a software74 was used. The data inserted into the equation were temperature, salt concentration and the transport properties, as the parameters, which were obtained by MD simulation. These parameters were: the diffusion coefficients, ionic conductivities and transference number. The other parameters used, were chosen as the realistic values and state-of-the-art for LIB compositions, as summarized in the Table S 11 of the supporting Information. Time-dependent solver and a direct solver of MUMPS74,75 with absolute tolerance of 0.001 were utilized to solve relevant equations. The geometry was discretized automatically by employing a physicscontrolled mesh and refined until sufficient accuracy was achieved with finer size. Then, the discharge-charge cycle was applied to determine the effect of electrolyte properties on performance of LIB. As shown in Figure 9 (a), the cycle applies to: 38 ACS Paragon Plus Environment

Page 38 of 55

Page 39 of 55

100 s of discharge at 5 C (as defined in the Table S 11 in the supporting Information) for current density (as shown by curve A to B), 20 s at open circuit (as shown by curve C to D), then 100 s of charge at 5 C current density (as shown by curve E to F) and finally open circuit conditions (as shown by curve G to H). Cell voltage and total polarization profiles are depicted in Figure 9 (b) and (c), respectively. The cell with conventional electrolyte [(EMC+EC 30) + LiPF6] experiences higher polarization than the electrolyte systems containing SN and FEC as the additives. Figure 9 (b) and (c) shows more evidently the performance of these electrolyte system, whereas Walden plot, as shown in Figure 7, was not capable to manifest the difference in performance of these additives and even could lead to erroneous interpretation in regard to behavior of additives in the studied electrolytes. 100

4.6 A

F

B

80

4.5

60

4.4

40

4.3

Cell Voltage (V)

Current Density (A/m2 )

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

20 D

0

G

H

C

-20

[(EMC+EC 30) + LiPF6 ] [(EMC+FEC 10) + LiPF 6 ] [(EMC+SN 10) + LiPF6 ]

4.2 A

H

4.1 E 4

-40

3.9

-60

3.8

C

-80

3.7 E

-100 0

50

100

G

D

B

F 150

200

250

3.6

300

0

Time(s)

50

100

150

Time(s)

(a)

(b)

39 ACS Paragon Plus Environment

200

250

300

The Journal of Physical Chemistry

0.4 [(EMC+EC 30) + LiPF6 ]

B 0.3

[(EMC+FEC 10) + LiPF 6 ] [(EMC+SN 10) + LiPF6 ]

0.2

Total Polarization (V)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 40 of 55

C

0.1 A

D

0

H E

-0.1 -0.2

G

-0.3 F

-0.4 0

50

100

150

200

250

300

Time(s)

(c) Figure 9. Discharge-charge cycle performance of one dimensional Li-ion battery (based on P2D electrochemical model) operated with different electrolytes at temperature 20 N'3 (a) The applied current density (b) Cell voltage (c) Total polarization.

Figure 10 (a, b, c) have been selected as examples to show the effect of cell variations including Li+ concentration in the electrolyte, Li concentration in active material of electrodes (LiyMn2O4 and LixC6 ) and local current density on the performance of LIB by color-bars. The results illustrated in Figure 10 corresponds with Figure 9 for electrolyte systems containing (EMC +EC 30). Similar results for electrolyte systems containing (EMC+FEC 10) and (EMC+SN 10) are depicted respectively in Figure S 6 and Figure S 7 of the Supporting Information. In Figure 10 (a and b) the intercalation process5,6 during discharge-charge cycle is demonstrated. In discharge, the metallic lithium in the anode (LixC6) loses an electron and changes into lithium ion (Li+) and as result Li concentration in the anode decreases, whereas Li+ concentration in the electrolyte increases. As it is seen in Figure 10 (a), from A to B, in the region 0