A finite strain electro-chemo-mechanical theory for ion transport with application to binary solid electrolytes

A finite strain electro-chemo-mechanical theory for ion transport with application to binary solid electrolytes

Journal of the Mechanics and Physics of Solids 125 (2019) 681–713 Contents lists available at ScienceDirect Journal of the Mechanics and Physics of ...

4MB Sizes 0 Downloads 5 Views

Journal of the Mechanics and Physics of Solids 125 (2019) 681–713

Contents lists available at ScienceDirect

Journal of the Mechanics and Physics of Solids journal homepage: www.elsevier.com/locate/jmps

A finite strain electro-chemo-mechanical theory for ion transport with application to binary solid electrolytesR Markus Ganser a,b,∗, Felix E. Hildebrand a, Marc Kamlah b, Robert M. McMeeking c,d,e a

Robert Bosch GmbH, Corporate Sector Research and Advanced Engineering, Stuttgart 70465, Germany Institute for Applied Materials Karlsruhe Institute of Technology, 76344 Eggenstein-Leopoldshafen, Germany c Materials Department and Department of Mechanical Engineering, University of California Santa Barbara, CA 93106, USA d School of Engineering, University of Aberdeen, King’s College, Aberdeen AB24 3UE Scotland, United Kingdom e INM-Leibniz Institute for New Materials, Saarbruecken 66123, Germany b

a r t i c l e

i n f o

Article history: Received 2 August 2018 Revised 9 December 2018 Accepted 4 January 2019 Available online 11 January 2019

a b s t r a c t An electro-chemo-mechanical formulation of ion transport in solid electrolytes, in particular for binary systems, is presented. Starting with conservation laws and the second law of thermodynamic, we state a consistent Helmholtz-energy-based framework taking electrostatics, component transport and nonlinear elastic mechanical interaction into account. With the help of finite strain continuum mechanics, we include the effect of geometry changes on ion transport. Changes of local concentration cause swelling and shrinkage and hence stress assisted diffusion. Further coupling originates via an osmotic pressure. Since binary systems are of special interest in battery applications, we formulate both, a fully resolved and an electroneutral model for ion transport. The latter turns out to be an extended version of Newman’s concentrated solution theory taking mechanical effects into account. We demonstrate the importance of these mechanical effects by means of double layers adjacent to blocking electrodes and concentration profiles during galvanostatic charging. Further, we investigate the effect of external deformation as, e.g. found in dendrite growth. © 2019 Elsevier Ltd. All rights reserved.

1. Introduction Among coupled transport theories, the transport of ionic species in a deformable, solid host is particularly interesting both with respect to possible applications and from a modeling perspective. In terms of application, such solid electrolytes enable state of the art energy storage (e.g. solid state batteries (Yang et al., 2017), energy conversion (e.g. fuel cells (Kreuer et al., 2004)) and actuation or sensing (e.g. ion-exchange polymer-metal composites (IPMC) (Pugal et al., 2010)). In terms of modeling, electrical, chemical as well as mechanical aspects have to be considered with particular emphasis on the coupling between these physical domains. The coupling between electrical and chemical domains has thereby become its own discipline in the form of electrochemistry, whereas the coupling with mechanics is less well established. A direct experimental validation of the existence and impact of this coupling are IPMC, sometimes referred to artificial muscles, that are used R ∗

Dedicated to the memory of Christian Miehe Corresponding author at: Robert Bosch GmbH, Robert Bosch Campus 1, 71272 Renningen, Germany. E-mail address: [email protected] (M. Ganser).

https://doi.org/10.1016/j.jmps.2019.01.004 0022-5096/© 2019 Elsevier Ltd. All rights reserved.


M. Ganser, F.E. Hildebrand and M. Kamlah et al. / Journal of the Mechanics and Physics of Solids 125 (2019) 681–713

Fig. 1. Graphical representation of the coupling between of multicomponent ion transport and mechanics in a Solid State Electrolyte. Mechanical response in the form of geometry change and stresses is triggered by swelling or external forces. The ions are subject to this change of geometry and mechanical driving forces give rise to stress driven diffusion.

as actuators and sensors (Pugal et al., 2010) and show very strong mechanical responses, usually bending of the material (Enikov and Seo, 2005; Stalbaum et al., 2015). With the rise of solid electrolytes in lithium ion batteries and the related challenges of solid separators in contact with lithium metal anodes, the “electro-chemo-mechanical” coupling has gained increased interest over the last decade. Oftentimes, this coupling is incorporated by extending an existing electrochemical model. Here, we take the approach of a fundamental derivation of the fully coupled transport theory using the framework of nonlinear continuum mechanics to include a full picture of the electrical, chemical and mechanical response and their coupling, which goes beyond multicomponent mass transport in fluids (Krishna and Wesselingh, 1997). Despite the broad range of applications of solid electrolytes, we will present our proposed framework in the context of solid state batteries to fix ideas and to motivate modeling choices. In this context, solid state electrolytes (SSE) present a promising opportunity for increased capacity and safety. Besides the substitution of flammable liquids, solid electrolyte separators demonstrate the potential for suppressing the growth of dendrites on metal anodes (Yang et al., 2017), Lin et al. 2017, thus enabling the use of lithium metal as the anode material with its class leading energy density and low electrochemical potential (Lin et al., 2017). A broad class of electrolytes, ranging from inorganic ceramics to polymers, are currently under investigation (Goodenough and Singh, 2015). Polymer systems like poly(ethylene oxide) (PEO) serve as a solvent for binary ion transport, with mobile anions and cations. The dynamics of multicomponent transport within this host material contribute significantly to the cell performance. Several polymer modifications have been suggested in order to increase electrochemical properties as well as mechanical strength (Golodnitsky et al., 2015; Tang et al., 2016), e.g. via cross linking (Khurana et al., 2014) or block/triblock polymers Singh et al., 2007(Niitani et al., 2005; Stone et al., 2012; Young et al., 2014), respectively. The class of single ion conducting electrolytes, such as oxidic or sulfidic ceramics, represent a subproblem of a binary electrolyte with stationary anions and is not the focus of this work. The purely electrochemical behavior of binary electrolytes is usually treated using the well-established concentrated solution theory (Newman and Thomas-Alyea, 2004) which was initially derived for liquid electrolytes. This theory has also been used to describe binary solid state electrolytes (SSE) such as PEO or PS-PEO and the respective parameterization have been published for these materials (Doyle and Newman, 1996; Ma, 1995; Pesko et al., 2017; Timachova et al., 2018). Since the concentrated solution theory is the standard approach for transport description in the battery modeling community, we will rewrite our theory in this form to highlight the additional terms due to the coupling with mechanics. Also, we will present a rule for the conversion of parameters between this theory and our proposed framework. Within this work, we establish a consistent model of ion transport in a deformable host in the spirit of Bucci et al. (2016). We exploit concepts of finite strain continuum mechanics, consider rigorously electrostatic influences and use standard notation of electrochemistry thereby making a connection to well established results in that area. With the help of a thermodynamic framework of Helmholtz energy and a series of consecutive steps deriving the physics and electrochemistry, a clear and consistent picture of the interdependencies of electrochemistry and mechanics is drawn, leading to a combination of ideas from multicomponent mass flow (in contrast to porous media theory (Ehlers et al., 2009; Leichsenring et al., 2017)) under the influence of an electric field with finite deformation theory for the mechanical response of the host. Both a treatment based on the Gauss law for a system containing unbalanced charges (see also Salvadori et al., 2014; Salvadori et al., 2015a; Salvadori et al., 2015b) and an electroneutral version for a binary electrolyte are presented, thereby establishing a link to the well known concentrated solution theory (Newman and Thomas-Alyea, 2004) similar to Monroe and Delacourt (2013) by introducing local electroneutrality and by linking the commonly used transport parameters conductivity, transference number, diffusivity and thermodynamic factor to mobilities and chemical potentials. Non-idealities usually associated with a concentrated solution are thereby accounted for. Throughout the derivation, specific care is taken in considering coupling mechanisms between ion transport and mechanics. These coupling mechanisms include swelling, stress driven diffusion and the change of geometry on which the transport takes place, see Fig. 1. We are convinced that swelling is present in binary solid electrolytes since the salt, which is redistributed due to concentration polarization during charging, makes up a significant fraction of the electrolyte. These coupling mechanisms between transport and mechanics are also present in IPMC (Enikov and Seo, 2002) and hydrogels (Hong et al., 2008), where small molecules diffuse in a network of long chain molecules and where a change of local component concentration leads to heterogeneous swelling or shrinkage of the host material and consequently to mechanical stresses. These phenomena can, in turn, provide an additional driving force for

M. Ganser, F.E. Hildebrand and M. Kamlah et al. / Journal of the Mechanics and Physics of Solids 125 (2019) 681–713


Fig. 2. An undeformed geometry provides the material description where all balance equations are solved. The deformed geometry as well as spatial field variables are linked to material counterparts via the deformation gradient F.

diffusion. Similar relations between swelling and transport are also found in the context of battery active materials (Bohn et al., 2013; Bower et al., 2011; Christensen and Newman, 2006; Cui et al., 2012; Huttin and Kamlah, 2012; Klinsmann et al., 2016; Salvadori et al., 2018; Walk et al., 2014; Zhang and Kamlah, 2018). The electro-chemo-mechanical coupling gains specific importance if inhomogeneous external forces are applied on the solid electrolyte, as is the case when lithium inhomogeneously deposited on an adjacent lithium metal anode leads to non-planarities (Ferrese and Newman, 2014) or dendrite growth (Brissot et al., 1999; Harry et al., 2016). Such stresses and their coupling with ion transport (and interfacial reaction kinetics) play an important role when considering the morphological stability of the lithium/SSE interface (Tikekar et al., 2016). Swelling of active particles, infiltration of stiff active particle into the separator due to prestresses at cell assembly or even by an undesired damage of the battery also lead to deformation and stresses within the electrolyte. All these examples underline the need for a consistent, fully coupled electro-chemo-mechanical framework. Thermodynamic consistent transport models have to follow, especially in a coupled setting, two important modeling paradigms: The first important ingredient is the description of the constitutive functions (stress, chemical potential, etc.) which drive e.g. mechanical displacement or mass transport. To derive individual thermodynamic restrictions for these quantities, we utilize the Coleman-Noll argument (Coleman and Noll, 1963). As the resulting framework consistently links e.g. free energy, stresses, fluxes etc., it is straight forward to derive consistent interdependencies of different physical effects of e.g. swelling and stress driven diffusion or chemical diffusion and osmotic pressure. The second ingredient describes the cross-relations between different driving force/mass flux couples and is commonly satisfied using the Onsager theory (Landstorfer and Jacob, 2013; Latz and Zausch, 2011; Monroe and Newman, 2009), the Maxwell-Stefan formalism (Krishna and Wesselingh, 1997; Monroe and Delacourt, 2013) or the Poisson-Nernst-Planck equation. The structure of the work is as follows. In Section 2, we introduce field variables in the context of finite deformation. Section 3 deals with fundamental conservation laws considering chemical, electrical and mechanical contributions. Based on energetic considerations, the relationship of Helmholtz energy and constitutive material properties, and hence their interdependencies is derived. Section 4 describes a specific Helmholtz energy and discusses the consecutive material laws. The set of equations is closed by a thermodynamically consistent formulation of mass fluxes. In Section 5, we consider a binary electrolyte and introduce for this system electroneutrality to obtain an extended set of transport equations similar to those in concentrated solution theory (Newman and Thomas-Alyea, 2004). In Section 6, we present three numerical examples to highlight and investigate the influence of mechanics in (1) the double layer regime; (2) the polarization of concentration and (3) the complex transport behavior arising in the case of external deformation. An appendix provides the reader with more context for large strain continuum mechanics, Maxwell-Stefan theory, the derivation of transport parameters and steady state current. 2. Field variables For the description of the fully coupled electro-chemo-mechanical system we distinguish between a solid host α = 0 and N species with index α = 1 . . . N that can move relative to the host. The electroneutral host provides the mechanical strength whereas the (ionic) mobile components act as carriers of charge. 2.1. Kinematics We employ the framework of large strain continuum mechanics, see e.g. (Holzapfel, 20 0 0; Truesdell and Noll, 2004), and describe the deformation of the host material at time t by a mapping ϕt between a reference configuration B0 with material coordinates X ∈ B0 and deformed configuration Bt with spatial coordinates x ∈ Bt , such that

ϕt :=

B0 × T → Bt (X , t ) → x = ϕ (X , t ).



M. Ganser, F.E. Hildebrand and M. Kamlah et al. / Journal of the Mechanics and Physics of Solids 125 (2019) 681–713

Accordingly, the material and spatial velocity of the host material are given by

V0 = ϕ˙ t (X ) :=

∂ ϕ (X , t ) ∂t


v0 = ϕ˙ t (ϕt−1 (x )).


The displacement of the host is given by u = x − X . We emphasize that the host displacement will serve as reference to which the species transport will be defined. To describe local stretches and rotations, we make use of the deformation gradient F defined by


∂ ϕt ∂x = = Grad(x ), ∂X ∂X

Fi j = Grad j (xi ) =

∂ xi . ∂Xj


The determinant J = det F thereby determines the mapping between material volume elements dV ⊂ B0 and their deformed spatial counterpart dv ⊂ Bt , i.e.

dv = detF =: J > 0, dV


where the positivity constraint follows to avoid interpenetration of matter. Note that all quantities and equations can be formulated both in the referential and the spatial setting. Since the physical interpretation is more intuitive in the spatial setting, we will usually formulate governing equations in the spatial form and transform them to the material formulation, which is the natural setting for a numerical treatment. From here on, whenever possible, we use lower case symbols for the spatial quantities (e.g. x) and upper case symbols or diamond superscript (e.g. X, ρα ) for the material quantities. 2.2. Chemical species We assume that N species are distributed each as a continuum in the solid matrix consistent with a mixture theory. If the number of moles of species α in an undeformed volume element dV is given by dNα and its molar mass is Mα , we can define the material molar concentration and the partial material density respectively as

cα =

dNα , dV

ρα =

dNα Mα = cα Mα dV

∀α = 1..N.


The spatial quantities follow similar definitions and are linked to the material concentration and density by Eq. (5), yielding

cα =

dNα c = α, dv J

ρα =

dNα Mα ρ = cα Mα = α dv J

∀α = 1..N.


The spatial concentration cα describes the amount of ions of type α per deformed area dv. We assume that at each point of the host material species α moves with velocity Vα . Recall that the velocity of the host is defined via Eq. (2) and generally vα = ϕ˙ t . The overall densities of the mixture, including the host material in the undeformed and deformed configurations, respectively, are given by  ρtot =


α =0

ρα ,

ρtot =


ρα .


α =0

2.3. Charged species Since species α = 1.N can carry charge, whereas the host material is assumed to be charge neutral, the material and spatial charge densities for each species are

ρqα = F zα cα ,

ρq α = F z α c α ,


where zα is the charge number associated with ion α . The total charge density is computed via

ρq =


α =1

ρqα ,

ρq =


α =1

ρq α .


3. Balance and field equations In this section, we derive all required field equations to describe ion transport in a solid host material. Using the framework of finite strains, we identify the balance equations in the spatial configuration and formulate them in the reference frame with the help of the transformations given in Appendix 8.1.

M. Ganser, F.E. Hildebrand and M. Kamlah et al. / Journal of the Mechanics and Physics of Solids 125 (2019) 681–713


3.1. Balance of mass For the balance of mass, we distinguish between the host (α = 0) and the chemical species (α = 1 . . . N). The host material is assumed to have constant mass which yields (Holzapfel, 20 0 0)

ρ˙ 0 = 0.


To derive the balance of a single species α , we consider at time t an arbitrary control volume Ptα ∈ Bt with species velocity vα . The only change in mass of α in that volume occurs due to a mass source rα associated with chemical reaction and we can write

d dt


cα d vα =


rα d vα

∀α = 0.


To rewrite this in terms of the control volume Pt tied to the motion of the host, we use the Reynolds theorem (97), i.e.

d dt


cα d v +

∂ Pt

cα (vα − v0 ) · nda =


rα d v

∀α = 0.


The second term on the left hand side denotes the mass flux j α = cα (vα − v0 ) through the boundary ∂ Pt . Hence, we describe naturally the mass flux of a species relative to the velocity of the host v0 . That distinguishes our description from multicomponent transport in liquids which usually utilize the molar mixture velocity (Krishna and Wesselingh, 1997), barycentric velocity (Dreyer et al., 2015) or the solvent velocity (Monroe and Delacourt, 2013). Now, we apply a transformation to the material setting, the divergence and localization theorem, to obtain the local form of the balance of species mass and total mass, respectively, as

c˙ α + Div(J α ) = Rα

∀α = 0

 ρ˙ tot + Div(J tot ) = 0,


∂ (J ) with the definition of the divergence operator Div(J α ) = ∂ Xα i . The material species and total molar flux are J α = JF−1 j α i

 and J tot = N α =1 J α . The material concentration source is Rα = Jrα and since a chemical reaction does not produce mass, we N have α =1 Rα = 0, which gives the second equation. 3.2. Balance of charge and Gauss law As we assume that charge is exclusively transported by ions without involvement of electrons, the balance of charge can be directly derived from the balance of mass

  ρ˙ qα + Div J qα = Rqα

∀α = 0

  ρ˙ q + Div J q = 0,

(14) N

with charge flux J qα = F zα J α and source Rqα = F zα Rα . Since charge cannot be generated in reactions, we have α =1 Rqα = 0,  which —together with the definition of current J q = N α =1 J qα — yields the total balance of charge on the right hand side. Note that Eq. (14) is not independent of the balance of mass as seen in Eq. (13). As a result, only N − 1 mass fluxes Jα are independent for a given total current flux J q . The Gauss law relates the spatial electric displacement d to the total charge density ρ q . The spatial and material representations are

div(d ) = ρq

Div(D ) = ρq

(15) JF−1 d.

with the material electric displacement D = Further we follow (McMeeking and Landis, 2005) and introduce a surface charge ωSC = D · N as the jump of electric displacement. Since we assume an electroneutral host, this surface charge d will be assumed to be zero. The electric field e is related to the magnetic field by curl(e ) = dt b = 0 due to Maxwell’s equations, where we have assumed that the change of magnetic field b in time can be neglected. To satisfy Maxwell’s equation1 , we define the electric potential  such that

e = −grad() = −F−T Grad() = F−T E


with E = −Grad() the material electric field. The electric field is related to the electric displacement by d = 0 e + p = 0 r e, where  0 and  r are the vacuum and relative permittivity, respectively. By inserting Eq. (16) in Eq. (15), we obtain the well known Poisson equation

−div(0 r grad() ) = ρq


which forms the acoustic limit of the full Maxwell equations. In the case of electroneutrality, the electric potential is decoupled from individual mobile components and Ohm’s law

j q = κe = −κgrad()


in terms of the total current j q becomes valid. In such circumstances the current is driven by an electric field through a conductivity κ as shown in Eq. (18). 1

We make here use of the identity curl(grad(ξ ) ) = 0

∀ξ .


M. Ganser, F.E. Hildebrand and M. Kamlah et al. / Journal of the Mechanics and Physics of Solids 125 (2019) 681–713

3.3. Conservation of linear momentum To state the balance of linear momentum, we sum the rates of the partial linear momenta of host and species and equate E C them with mechanical (e.g. gravity), electrical and chemical body forces of each species bM α , bα , bα and effective tractions tM , tE , tC acting on the host material and obtain

 N  N   d  ρα v α d v α = bM + bEα + bCα dv + t M + t E + t C da. α dt Ptα Pt ∂ Pt α =0


α =0

We have to carefully apply the time derivative on left hand side. Consider a snapshot at time t and choose the component domain Ptα such that it aligns with host domain Pt := Pt0 . Since both domains move with a different velocity, the increments α . To resolve that issue and exchange integral and time derivative, we apply Reynolds theorem and describe differ Pt+δt = Pt+ δt the state in the reference geometry P0 and with use of Eq. (101) obtain for the left hand side of Eq. (19)

d dt

ρα vα dvα =





d (ρ  Vα ) + Div(Mα Vα J α )dV dt α


ρα V˙ α + Grad(Mα Vα ) · J α + Mα Vα Rα dV, 



where we utilized for the second line with Eq. (13) the balance of mass. The body force BVα corresponds to an inertia force of component α within the host material.  E  E E The total electric body force BEtot = N α =0 Jbα is expressed as a Piola-type Maxwell stress tensor Btot := Div P . The same procedure is used for the chemical body forces. Further, mechanical, electric and chemical traction are balanced via a jump in mechanical stress, see (McMeeking and Landis, 2005)

∂ Pt

t M + t E + t C da =

∂ P0

T M + T E + T C dA = −

 ∂ P0

PM  · N dA =

Div PM dV.



Consequently, an expression for the balance of momentum of the host material is

ρ ϕ¨t =  0

BM tot



+ Div P + P + P




α =1



N N M M where we have introduced the effective total body force BM α =0 Bα = α =0 Jbα and the total Piola stress P. The inertia, tot = advection and reaction contribution of the mobile components can be regarded as an additional body force BVα acting on the host material. 3.4. Conservation of angular momentum We elaborate the tedious derivation for the balance of angular momentum in Appendix 8.2 and mention only the result:



α =1

Mα J α  Vα = FPT −


α =1

Mα Vα  J α .


Without component motion this aligns with the standard result PFT = FPT . 3.5. Conservation of energy Consider the balance of energy in an arbitrary volume Pt ∈ Bt . The change of total internal energy Eint , as well as kinetic energy Ekin , is equal to the power supply in Pt due to heat LH , mechanics LM and electricity LE :

 d  int E + E kin = LH + LM + LE dt


The left hand side of Eq. (24) describes the rate of internal energy eα and kinetic energy

E int :=

N   α α =0 Pt

ρα e α d v α


E kin :=

N   α α =0 Pt

1 ρα vα · vα dvα . 2

1 2 ρα vα

· vα of each species


M. Ganser, F.E. Hildebrand and M. Kamlah et al. / Journal of the Mechanics and Physics of Solids 125 (2019) 681–713


Note since each component α moves with its own velocity vα , we have to again apply Reynolds transport theorem through Eq. (97) and obtain N  d int E = dt

α =0



d dt


ρα e α d v +

 ∂ Pt

ρα eα (vα − v0 ) · nda

d (ρ  eα ) + Mα Grad(eα ) · J α + Mα Rα eα − Mα c˙ α eα dV, dt α

α =0 P0


where we have used Eq. (103) to obtain the second version, involving the transformation to the material configuration, the divergence theorem and balance of mass. The same procedure applies for the kinetic energy N  1 d kin  E = BVα · Vα − Vα · Vα Mα Rα dV. dt 2 P0


α =0

where we recall the definition of BVα in Eq. (20). On the right hand side of Eq. (24), we observe energy production due to heat sources and heat fluxes. In terms of heat we treat host and components as one system and get the contribution

LH :=

rH dv −



∂ Pt

j · nda,


where we have introduced a homogenized heat flux jH and an external heat source rH , e.g. an electromagnetic heating. Note that heat generated due to chemical reactions is an internal process of conversion within the internal energy. The reference description of Eq. (28) is

LH =



RH − Div J H dV.


The rate at which mechanical work is done on the host material and the species is N  

LM :=

α =0 Pt

bM α · vα dv +

∂ Pt

t M · v0 da,


where only the host material can supports traction forces. The individual body forces act on the host material as well as on each mobile component. The result converts to

LM =


bM α · (vα − v0 )dv +

α =0 Pt



α =0 Pt

 bM α · v0 dv +

∂ Pt

t M · v0 da,

N    1 T M ˙ F Bα · J α + BM tot + Div (P ) · V0 + P:FdV, cα


P0 α =0

where we have made use of Eq. (21), the divergence theorem and the notation of double contraction, meaning A:B = Ai j Bi j to obtain the version with respect to the reference configuration. The electrical energy supply due to transport of charged species at an electric potential  by means of a current j q reads

 LE := −

∂ Pt

 jq · nda =



e · j q − div j q dv.


Transformation to the reference formulation yields

LE =



E · J q + Div D˙ dV =


E · J q + E · D˙ dV +

∂ P0

D˙  · N dA


where we have used Div J q + D˙ = 0, a consequence of Eqs. (14) and (15). The last term can be identified as a surface

charge ωSC = D · N , see e.g. (McMeeking and Landis, 2005). Since all free charges are carried by ions, we assume that the change of surface charge is zero and hence ωSC = 0 on ∂ P0 . Combining Eqs. (27), (29), (31), (33) and using balance of momentum through Eq. (22), we obtain the local material form for the balance of energy N N    d   ρα eα = RH − Div J H + P:F˙ + E · D˙ + E · J q + Mα eα c˙ α dt

α =0

α =1


α =1

 N    1 ˜ α + Grad(Mα eα ) · J α − B Mα eα − Mα Vα · Vα Rα 

α =1



V ˜ α = − 1 FT BM with the extended specific body force B α + Bα . Additional to the thermal and mechanical contributions, we cα see a power supply due to electric interaction, e.g. Joule heating, species transport and due to chemical reaction.


M. Ganser, F.E. Hildebrand and M. Kamlah et al. / Journal of the Mechanics and Physics of Solids 125 (2019) 681–713

3.6. Balance of entropy We proceed similarly with the balance of entropy and sum over the component entropies ηα . We locally assume equal absolute temperature ≡ α for each constituent and obtain

   N  H d  rH j ρα η α d v α = dv − · nda + dt Ptα Pt ∂ Pt Pt α =0

ρ0 δV dv


where δ V is the dissipation density. Applying Reynold’s theorem through Eq. (97) on the left hand side, divergence theorem on the right hand side and switching to the material frame gives N d   ρα η α = dt

α =0

ρ0 δV + −



   1 (RH − Div J H + J H · Grad( )

[Grad(Mα ηα ) · J α − Mα ηα c˙ α + Mα ηα Rα ]

α =1


as an expression for the rate of total entropy. Note that Eq. (36) is a representation of the Clausius-Duhem inequality (Holzapfel, 20 0 0). 3.7. 2nd law of thermodynamics We now introduce the total Helmholtz energy ψ˜ and the species Helmholtz energy ψ α consisting of the Legendre transform components of the internal energies eα

ψ˜ (F, cα , , D ) =


α =0

ρα ψα =


α =0

ρα (eα − ηα ).


Combining Eqs. (34) and (36), making use of Eq. (37) and solving for ρ0 δV , we obtain the so-called full dissipation inequality, also known as the Clausius-Planck inequality as

      N ˜  ∂ ψ ∂ ψ˜ ∂ ψ˜ ∂ ψ˜ ρ0 δV = − ηα − ˙ + P − : F˙ + Mα ψα −  c˙ α + E − · D˙ ∂ ∂F ∂ cα ∂D α =0 α =1 N 


α =1

˜ α + Grad(Mα ψα ) + F zα E − (Mα ηα )Grad( ) · J α − B

− Mα ψα −


Grad( ) · J H

1 Mα Vα · Vα Rα . 2


The left hand side has to be non-negative with 0 ≥ ρ0 δV due to the second law of thermodynamics. Following the Coleman, D ˙ , F˙ , c˙ α ˙ and therefore require the parenNoll argument (Coleman and Noll, 1963), we assume independence of the rates thetical terms in Eq. (38) to vanish. This yields the consistency conditions N 

α =0

ηα = −

∂ ψ˜ , ∂


∂ ψ˜ , ∂F

μα := Mα ψα =

∂ ψ˜ , ∂ cα


∂ ψ˜ . ∂D


where we have introduced the chemical potential (per mole) as μα := Mα ψ α of species α , which is also sometimes referred to as the convective potential (Biot, 1977), a measure for the energy required to insert a mole of species α into a fixed volume dV. Consequently the total Helmholtz energy of Eq. (37) is a combination of the host energy and the sum of the chemical potentials of the components

ψ˜ = ρ0 ψ0 +


α =1

cα μα .


The second term on the right hand side is normally identified as Gibbs energy G. However this form of Gibbs energy is formulated for a system in which internal energy is a function only of extensive variables and Euler’s theorem for homogeneous functions utilized. This is not the case in the context of solid mechanics where the internal energy is expressed as a   2 function of deformation which is not extensive. Thus N α =1 cα μα cannot, in this case, be identified as a Gibbs energy. 2

Approaches to extend the Gibbs energy by deviatoric contributions are found in e.g. (Baker et al., 2016; Goyal and Monroe, 2017).

M. Ganser, F.E. Hildebrand and M. Kamlah et al. / Journal of the Mechanics and Physics of Solids 125 (2019) 681–713


Further, since ψ˜ is composed from of the individual ψ α , the definition in Eq. (39) has an implicit structure and can be resolved via N ∂  ρ ψ  ∂ cα β =0 β β

Mα ψα =


β =0


∂ψβ =0 ∂ cα


β =1


∂μβ ∂ψ0 =−  , ∂ cα ∂ cα


where the last representation can be seen as a Gibbs-Duhem like equation at constant temperature and deformation. We continue with the evaluation of the dissipation inequality (38). If all equations of Eq. (39) are fulfilled, Eq. (38) reduces to



α =1


˜ α + μα + F zα  + Mα ηα Grad( ) · J α Grad B

Grad( ) · J H −


 1 μα − Mα Vα · Vα Rα 2

α =1


 where we have reformulated the Joule heating E · J q = − N α =1 Grad() · F zα J α . Eq. (42) constitutes a thermodynamic constraint for the molar fluxes Jα (and therefore species velocity vα ), the heat flux JH and the reactions Rα . Under the assumption of each process being independent, the first term on the right hand side of Eq. (42) relate molar flux to chemo-electrical driving forces, inertial and body forces. Similarly, when independent of other terms, the second term leads to the well known Fourier’s law of thermal conduction and the third term guides the direction of reactions. Assuming isothermality (Grad( ) = 0), quasistatic conditions (ϕ˙ t ≈ 0), the absence of body forces (BM α ≈ 0) and further assuming that inertia related forces are small (BVα ≈ 0). We obtain the thermodynamic constraints



α =1

Grad(ωα ) · J α




μα R α ,


α =1

where we have introduced the electrochemical potential

ωα = μα + F zα .


The first of Eq. (43) provides a constraint for all mass fluxes Jα as discussed in the upcoming sections. The second of Eq. (43) denotes a requirement for chemical reactions within the bulk, e.g. degeneration of electrolyte. This will not play a role in the following sections, where Rα = 0 will be assumed. 4. Material description To complete the model, we provide a definition of an appropriate energy functional, derive the constitutive equations and prescribe thermodynamically consistent mass fluxes. 4.1. Constitutive Helmholtz energy We assume that the Helmholtz free energy for the electro-chemo-mechanical system can be written with additive mechanical, electrical and chemical contributions as

ψ˜ (F, cα , D ) = ψ˜ M (F, cα ) + ψ˜ E (F, D ) +


α =1

ψ˜ αC (F, cα ),


where we have chosen that the Helmholtz free energy of the host at reference state is zero. 4.1.1. Electrostatics For the sake of simplicity, we use an isotropic formulation of electrostatics in the spatial domain consistent with Eq. (17), which leads to

ψ˜ E J


1 d · d. 20 r


The reference description is then given as

ψ˜ E =

1 ( )−1 :(D  D ) 20 r


r = JF−1 r F−T


We introduce the material description of the relative electric permittivity r to keep the notation concise. With Eq. (47) we do not include dissipative materials (Rosato and Miehe, 2014) or anisotropic polarization (McMeeking and Landis, 2005) and assume a concentration independent permittivity.


M. Ganser, F.E. Hildebrand and M. Kamlah et al. / Journal of the Mechanics and Physics of Solids 125 (2019) 681–713

Fig. 3. The deformation gradient is composed of purely isotropic swelling Fs due to the pressure of ions and elastic deformation Fe .

4.1.2. Mechanics We describe the mechanical response of the material by a hyperelastic relation and assume that local changes of concentration lead to volume changes in the sense of swelling. We thus assume a multiplicative split of the deformation gradient into swelling and elastic contributions (Fig. 3)

F = Fe Fs ,


where Fe describes the elastic response of the material and Fs corresponds to a volumetric stress free expansion or compression. We use a model of linear, isotropic swelling given by 1

Fs = Js 3 I


Js = 1 + β s


α =1

  α cα − cαref ,


ref denotes the unswollen state (i.e. the concentration at which the material was fabriwhere the reference concentration cα cated). The factor β s thereby accounts for non-ideality in swelling and can be determined experimentally or by atomistic  ref simulations (Cui et al., 2012). Due to the requirement of Js > 0, the reference concentration has to obey N α =1 α cα < 1. The impact of swelling depends strongly on the size of the ionic species which is described by the molar volume α .3 We assume a Neo-Hookean material law for the description of the elastic response of the host. The mechanical energy,  which depends solely on Fe = FF−1 s (cα ), then becomes

1 2

1 2

M ψ˜ M (F, cα ) := ψ˜ NH (Fe ) = γNH (I1 − 3 ) + λNH (log I3 )2 − γNH log I3


with the Lamé constants γNH , λNH and the invariants

I1 = tr FTe Fe ,

I3 = Je = det Fe


see e.g. Ogden (2011). Note, that the Lamé constants can be converted from Young’s modulus E or shear modulus G and the Poisson ratio ν by

Eν 2 Gν E = and γNH = = G. (52) (1 + ν )(1 − 2ν ) 1 − 2G 2 (1 + ν ) Hereby, E, G and ν have their conventional meaning relating infinitesimal strain relative to the current configuration and

λNH =

change in Cauchy stress. 4.1.3. Chemistry The chemical contribution to the Helmholtz energy is assumed to be

ψ˜ αC

= cα μ0α + R




f α cα fαref cαref

dcα .


with the gas constant R and the standard chemical potential μ0α which is independent of concentration and deformation. In concentrated solutions, a correction factor, the so called fugacity or activity coefficient, fα is introduced to incorporate non-idealities in the formulation (Atkins et al., 2006; Bazant, 2013). In dilute solutions, an ideal gas response is obtained with fα = 1. Definition (53) is formulated very similarly to classical mixture theories in fluids and solids, but there are two features we want to highlight. In contrast to conventional usage in solids, where vacancy motivated chemical potential is a function  Bohn et al. (2013), we follow a concept used for fluids and postulate that the chemical energy of material concentration cα depends on the amount of ions per deformed volume, i.e. on the spatial concentration cα . Therefore, material compression leads to a decrease of spatial volume and hence an increase of effective spatial concentration (Fig. 4). The importance of this assumption and its implications will be illustrated in Example 6.3. The second feature is the integration domain of the integral in Eq. (53) which represents the entropic contribution to the chemical energy. We choose it to be zero at reference ref . This will be beneficial for the mechanical response in the upcoming section. concentration cα 3

Note that with Eq. (49) there is an underlying assumption that the components and the host undergo the same volumetric deformation.

M. Ganser, F.E. Hildebrand and M. Kamlah et al. / Journal of the Mechanics and Physics of Solids 125 (2019) 681–713


Fig. 4. Consider an isolated system: The effective concentration cα changes from the undeformed (t1 ) to the deformed state (t2 ). Therefore, the chemical energy increases and gives rise to chemical stresses tC also known as osmotic stress or pressure.

4.1.4. Constitutive equations Following Eq. (39), we can deduce the stress tensor, the chemical potential and the electric field from Eqs. (47), (50) and (53). The Piola stress tensor then follows as

P = PM + PE + PC

  ⎧ M F−T ⎨P = F−1 s : (λNH log I3 − γNH ) e + γNH Fe  1 with PE = F−T E  D − 2 (E ·D )I ⎩PC = N −F−T R  cα 1 + ∂ ln fα dc . α α =1 ∂ ln cα cref




Note that the electric contribution is the Maxwell stress tensor as derived in Rosato and Miehe (2014) or McMeeking and Landis (2005)4 and reduces to an electric body force in the material and spatial representations


BEtot = Div PE = ρq F−T E

bEtot = ρq e.


The chemical contribution PC originates from the dependency of concentration on deformation in the chemical energy  ). Its equivalent Cauchy stress (Jcα = cα

σ = C


α =1

−R cα − cα + ref

cα cαref

 ∂ ln fα dc I . ∂ ln cα α


is purely hydrostatic. It corresponds to the osmotic pressure pC = − 31 tr σ C , which is usually known from a setting of a membrane permeable for a solution but not for a solvent. The jump in concentration across the membrane yields diffusion of solvent through it to minimize chemical energy. Diffusion through the membrane stops when a critical pressure across the membrane is reached that balances the chemical driving force (Hong et al., 2008). An ideal solution with fα = 1 and N ref = 0 gives the familiar result for an ideal gas pC = cα α =0 R cα (Atkins et al., 2006). In contrast to the case of fluids, we C ref ) = 0 is present at the reference concentration, thereby avoiding formulate Eq. (53) such that a stress free state with σ (cα further complexity with osmotic pressure having to be balanced in the reference configuration by a mechanical response. Now we evaluate the chemical potential through

μα = μCα + μM α with

 μCα =

∂ ψ˜ αC f α cα 0 ∂ cα = μα + R ln fαref cαref M ∂ ψ˜ s M μM α = ∂ cα = Je β α p .


The chemical part μCα is in alignment with standard chemical potentials with an energetic and an entropic contribution (Bazant, 2013; and Thomas-Alyea, 2004). The mechanical contribution μM α due to the elastic part of the pressure  Newman  pM = − 13 tr J1 PM FT leads to stress assisted diffusion (Larche and Cahn, 1984). The factor Je is a consequence of large strain considerations. Concentration dependent properties of elasticity and electrostatics would yield additional contributions to the chemical potential (Cui et al., 2012). Finally we evaluate the electric field from Eq. (39) and get


1 ( )−1 D, J0 r


which gives together with Eq. (15) the Poisson equation

−Div(0 r Grad() ) = ρq .


The description in the material configuration of Eq. (59) aligns with the well known spatial Poisson equation, see Eq. (17).


The authors in McMeeking and Landis (2005) formulate the mechanical response with Cauchy stresses σ =

1 PFT . J


M. Ganser, F.E. Hildebrand and M. Kamlah et al. / Journal of the Mechanics and Physics of Solids 125 (2019) 681–713

4.2. Mass flux Motivated by the requirement of the reduced dissipation inequality (43), we postulate a linear ansatz for the mass flux (Bucci et al., 2016; Ferguson and Bazant, 2012; Hong and Wang, 2013; Hong et al., 2008) as

Jα = −


β =1

Mαβ Grad

 ωβ ,


where we introduce a mobility matrix Mαβ , which we require to be positive semi-definite to guarantee non-negativity of Eq. (43). The mobility matrix

Mαβ :=

M11 M21 .. .

M12 M22 .. .

⎞   ⎠  

··· ··· .. .


aα · Mαβ · aβ ≥ 0

∀ aα , aβ ∈ R d



consists of N × N entries where each entry can be again a Rd×d tensor to represent anisotropic properties. A deeper discussion of the implementation of anisotropy within this theory can be found in Appendix 8.3. The mobility matrix describes the influence of the driving force Grad(ωβ ) on the mass flux Jα . A low mobility is a result of high internal friction with respect to the host (on-diagonal terms) and with respect to component-component interaction (off-diagonal terms). They are related to drag coefficients and binary diffusion coefficients in the Maxwell-Stefan-diffusion model as discussed in Appendix 8.4. We would like Eq. (60) to be consistent with the spatial relation

jα = −


β =1

Mαβ grad

 ωβ ,


correlating the mass flux on the deformed geometry with a spatial gradient of the electrochemical potential and obtain5

Mαβ = JF−1 Mαβ F−T .


The property Mαβ is the spatial mobility matrix, a parameter closely related to measurable quantities and Eq. (63) effectively informs the material configuration of changes of transport paths due to deformation. Eventually, we solve the material mass flux in Eq. (60) with the material mobility of Eq. (63) which is equivalent to the desired spatial mass flux of Eq. (62). If the domain is deformed, the material mobility Mαβ becomes anisotropic due to the conversion of Eq. (63), even if isotropic material properties are assumed. This anisotropy is solely related to deformation of transport paths. The requirement that the material mobility matrix has to obey a positive semi-definite structure leads along with Eq. (63) to the spatial mobility matrix having to be positive semi-definite. This property can be checked e.g. by the presence of non-negative eigenvalues. Therefore, the determinant of Mαβ yields a non-negative value. 4.3. Generic mass flux and current The driving force for mass transport in Eq. (60) is partially due to the chemical potential and partially driven by gradients of the electric potential

J α = J Cα + J Eα

J Cα = −



β =1 Mαβ Grad

J Eα = −κα Grad()



where we have introduced a mass conductivity

κα = F


β =1

Mαβ zβ ,


which defines the impact of the electric field on mass transport.6 A similar split for the electric current J q is

Jq =

J Cq


J Eq


J Cq = −


α =1 κ α Grad(μα ) 

J Eq = −κ Grad()


with the ionic conductivity

κ =


α =1

F zα κα


5 We deduce from j α · nda = J α · N dA the relation J α = JF−1 j α , see Holzapfel (20 0 0). Further, the material and spatial gradient operator are linked via Grad(· ) = ∂∂Xx grad(· ) = Fgrad(· ). 6 In the case of uncharged species—as found e.g. in active particles (Bohn et al., 2013)—the conductivities from Eq. (64) are zero. Consequently, the electric portion J Eα drops out and only mass transport due to chemo-mechanics remains.

M. Ganser, F.E. Hildebrand and M. Kamlah et al. / Journal of the Mechanics and Physics of Solids 125 (2019) 681–713


From Eq. (65) we can see that κ is directly related to the mobilities. Further, we introduce a so-called transference number tα as

tα J Eq = F zα J Eα

tα κ = F zα κα

tα = F zα κα (κ )−1 ,


as seen in Lai and Ciucci (2011) and Gouverneur et al. (2018) for binary systems. Physically, the transference number tα represents the fraction of the overall current transported by species α if there is only migration and no diffusion mechanism is present (Grad(μα ) = 0) (Evans et al., 1987). The transference number in Eq. (68) is therefore a representation of mobilities as discussed in Appendix 8.5 and thus a well defined material property. Note that transference numbers are in general   anisotropic, their sum is the identity tensor ( N α =1 tα = I ) and, based on our derivation, we find that there is no thermodynamic restriction demanding positivity of the transference number. In fact, negative entries in tα are possible when one component is more strongly affected by the chemical potential of another species, i.e. the component under consideration is dragged along by one or more other components.7 Recent experimentally work (Gouverneur et al., 2018; Pesko et al., 2017) confirms the theoretical result. We also refer to Appendix 8.6 which gives further insights on how to check thermodynamic consistency of transport parameters. In the following we assume that at least one species is charged (zα = 0). In such systems, the conductivity and transference number are measurable characteristic parameters and we can rewrite the current as

J q = −κ

N  1  t Grad(μα ) − κ Grad(). F zα α


α =1

Next, we rewrite the electric part of the mass flux of Eq. (69) as

J Eα =

 1  tα J q − J Cq , F zα


which gives us the possibility of decomposing the species mass flux into three components

J α = (I − tα )J Cα − tα

N  zβ C 1  J + t J zα β F zα α q


β =α =1

namely, a self diffusion, a drag diffusion and a migration portion, respectively. Note that Eqs. (69) and (71) are generic versions of Newman type equations for transport in concentrated solutions (Newman and Thomas-Alyea, 2004). (6 )

(48 )

 /J = c  / ( J J ) Care has to be taken when evaluating the gradient of the chemical potential of Eq. (57). With cα = cα e s α due to the effect of compression this gives rise to

Grad(μα ) = ξαc with the prefactors

1 ξα = 1− J c

  ∂μCα Grad(cα ) + ξ p α Grad pM + ξαJe Grad(Je ) ∂ cα

 ∂ ln Js , ∂ ln cα

ξ p = Je β s ,

ξαJe = β s α pM −

cα ∂μCα . Je J ∂ cα



The first term on the right hand side of Eq. (72) drives classical diffusivity, the middle term is the mechanical driving force due to pressure gradients, whereas the last term specifies a second order effect due to volume change. In other words, if the volume changes due to elastic deformation, the ions are packed closer and the chemical energy increases. If the elastic deformation is inhomogeneous, the third driving force in Eq. (72) originates. 5. Binary system - link to classical concentrated solution theory Ionic charge in binary electrolyte systems is carried by a positively charged cation and a negatively charged anion which we denote by α = {−, +}. We assume here that z+ = 1 and z− = −1. In solid electrolytes, these ions are embedded in a solid host material, whose kinematics are described by the deformation map ϕt ; this is fundamentally different from liquids where a third mobile component, the solvent, is present. Therefore, we have only two mobile species and obtain the binary mobility matrix

Mαβ |α , β = {+, −}



M+ = M±


M± . M−


with three independent mobilities (that can each be functions of concentration, stress, etc.). The diagonal contains the cation mobility M+ and the anion mobility M− . The off-diagonal term M ± in addition to the fugacity stems from the non-ideality "


0.1 Consider an isotropic binary electrolyte with mobilities Mαβ = 0.2

eigenvalues, but the transference number is negative with

tiso +


− 17 .


0.2 1

(α = +, β = − ). This obeys thermodynamically consistency with positive


M. Ganser, F.E. Hildebrand and M. Kamlah et al. / Journal of the Mechanics and Physics of Solids 125 (2019) 681–713

of concentrated solutions and can be associated with the motion of clusters of both cations and anion with unlike charge that cause the constituting cations to be influenced by the anion potential and vice versa. In the settings of IPMC, this effect is known as electro-osmosis drag where water is shuttled along with hydrogen (Zhu et al., 2013; 2011). Another motivation for the off-diagonal term is a counting argument. There are four constitutive quantities, namely the chemical potential and the three mobilities. This is the same number of parameters as for the standard parametrization of binary ion transport (Newman and Thomas-Alyea, 2004), see Appendix 8.5. In the following, we will present three approaches computing the motion of the mobile species {c+ , c− , ρq }. The first is a straight forward application of the previously derived theory for binary electrolytes. The second version is reformulated in order to introduce the balance of charge and serves as the transition to the third, the electroneutral version that constitutes simplifications of the first two. All of these formulations describe transport and, hence, the underlying structure of the mechanics with the Piola stress tensor P from Eq. (54), the material mobilities Mαβ from Eq. (63) and swelling from Eq. (49), is unaffected and remain the same in all versions.  , c  , } 5.1. Full formulation in terms of {ϕt , c+ −

We state the balance of linear momentum, add to it mass conservation for positive and negative ions derived from Eq. (62) and include the Gauss law as


J + = −M+ Grad(ω+ ) − M± Grad(ω− ) J − = −M± Grad(ω+ ) − M− Grad(ω− )

 c˙ + + Div J + = 0    c˙ − + Div J − = 0 Div(P ) = 0 Div(D ) = ρq

D = −

 0 r Grad



where the electrochemical potentials are given by

ω + = μ+ + F 

ω− = μ− − F .



 , c , }. The The equations in (75) constitute the system of field equations necessary for the solution of the fields {ϕt , c+ −  , c and boundary conditions (BC) for ϕ , c , c and . initial boundary value problem requires initial conditions for c+ − t + − Typical electrochemical boundary values (BV) are those for blocking electrodes given by

−N · J + |∂ B0 = −N · J − |∂ B0 = 0,

˜, |∂ B0 = 


˜ is a prescribed, position dependent electrical potential on the boundary. Also typical are galvanostatic BC specified where  by

−N · J + |∂ B0 =

J app , F

−N · J − |∂ B0 = 0,


where a boundary current per undeformed area Japp is applied and only cations leave the system. The boundary current is usually given by reaction kinetics such as the Butler-Volmer equation (Ganser et al., 2019). Note that for global electroneu trality, it is necessary that ∂ P J app dA = 0. To within an additive constant the electric potential is fully described due to its 0 coupling to ion transport. Under galvanostatic charging given by Eq. (78), the system from Eq. (75) will eventually obtain a steady state where the anions are in equilibrium and their mass flux vanishes, i.e. J − = 0. This gives an expression of a steady state Ohm’s law as

 SS Jq


= −(κ )SS Grad




  −1  (κ )SS = F 2 M+ − M± M− M± .


Expression (79) aligns with the findings of Balsara and Newman (2015), see Appendix 8.7 and shows the importance of off-diagonal mobilities. The steady state conductivity would only be governed by the cation mobility if M± is assumed to be zero. It is worth noting that effects from mechanics and the chemical potential do not enter the calculation explicitly. In the case of an isotropic system (Mαβ = Mαβ I ), we can deduce from a positive steady state conductivity requirements on the entries of the binary mobility matrix from Eq. (74). Therefore, we write the steady state conductivity in the spatial form as

κ iso



 1 1  F(κ )SS FT = M+ M− − M2± I . J M−


The expression in the brackets is the determinant of the spatial mobility matrix. Since both, the determinant and the steady state conductivity, are non-negative, we conclude M− ≥ 0. Due to the structure of the determinant, one can further conclude that M+ has to be non-negative, too. 5.2. Towards balance of charge To reformulate Eq. (75), we recall the cation mass flux from Eq. (71) as


= − t− M+ + t+ M± Grad(μ+ ) + t+ M− + t− M± Grad(μ− )W +

t+ J . F q


M. Ganser, F.E. Hildebrand and M. Kamlah et al. / Journal of the Mechanics and Physics of Solids 125 (2019) 681–713


where we identify an ambipolar mobility matrix as mentioned e.g. in Ferguson and Bazant (2012) as

Mamb := t− M+ + t+ M± = t+ M− + t− M± = F 2 M+ M− − M± 2

  −1 κ .


Note that in contrast to the treatment in Newman and Thomas-Alyea (2004), all concentration dependencies are hidden in the mobilities.   with the field variable total charge density ρ via c = c − ρq and Further, we exchange the anion concentration c− q − + F obtain from Eq. (69) a statement for the charge flux

Jq = −

κ  F

t+ Grad(μ+ ) − t− Grad(μ− ) − κ Grad().


After some rearrangement and using  = ωF+ − μF+ we obtain


 c˙ + + Div J + = 0   ρ˙ q + Div J q = 0 Div(P ) = 0 Div(D ) = ρq

J + = −Mamb Grad(μ+ + μ− ) + tF+ J q      J q = − κF t+ − I Grad(μ+ + μ− ) − κF Grad(ω+ ) D = −

 0 r Grad



where the current is driven by the electrochemical potential of the cation, see Eq. (76). With respect to the BC, we exchange the boundary flux for the anion flux in Eq. (78) to give

−N · J q = J app .


5.3. Electroneutral version in terms of {ϕt , c , ω+ } Ion polarization occurs only in the double layer close to interfaces and results in strong electric potential gradients. The width of this domain can be approximated by the Debye length (Dreyer et al., 2015) and is usually of the order of nanometers. Away from this region, electrostatic forces are dominant, drawing unlike charges together, and therefore electroneutrality prevails. The relationship   c+ ≈ c−

(86)  c+


is an adequate approximation for simplifying the transport equations and the individual concentrations and can be  = c . Consequently, the charge density is zero ρ  = 0 and Poisson’s law decouples both replaced by a salt concentration c± q from mass transport (Landstorfer and Jacob, 2013).  and c by c , it is beneficial to replace  by the electrochemical potential of the cations Besides the replacement of c+ −

ω+ = μ+ + F z+ ,


because it is directly accessible in probe measurements (Newman and Thomas-Alyea, 2004). Consequently, we reduce the  , c , } to {ϕ , c , ω }. Conservation of charge, as given in the form of an extended Ohm’s set of field variables from {ϕt , c+ + − t law in Eq. (84), then serves as the governing equation for the electrochemical potential. Besides eliminating a field variable and thus saving computational effort, another advantage of the electroneutrality assumption is the reduction in the number of required parameters. This is significant especially since individual ion properties are very difficult to measure. With the molar volume of the salt ± = + + − , the chemical potential of the salt μC± := μC+ + μC− (Newman and Thomas-Alyea, 2004) and considering Eq. (72), we obtain


 c˙ + + Div J + = 0

J + = −D Grad(c ) +


t+ J F q

− ξ p ± G Grad pM − ξ±Je G Grad(Je )

Div J q = 0


J q = Dq Grad(c )− κF Grad(ω+ ) + ξ p ± Gq Grad p 

(88) + ξ±Je Gq Grad(Je )

Div(P ) = 0 with contributions of diffusion, migration, stress-assisted diffusion and volumetric effects. The prefactors are given by

D = ξ±c Mamb

Dq =

ξ±c  F

∂μC± (c ) , ∂c

I − t+

 ∂μC± (c )  κ, ∂c

G = Mamb =

Gq =



 1 I − t+ κ , F

∂μC± (c ) ∂c


D ,




M. Ganser, F.E. Hildebrand and M. Kamlah et al. / Journal of the Mechanics and Physics of Solids 125 (2019) 681–713

with the Fickian diffusion coefficient D formulated in the material setting. Choosing the chemical potential μC± to be of the form of Eq. (57) with mean salt activity coefficient f± =


( f− )( f+ ) (Atkins et al., 2006), we obtain

∂μC± (c ) 2R ∂ ln f± (c ) = 1+ , ∂c c ∂ ln c


where expression in brackets is usually called the thermodynamic factor. The set of Eq. (88) constitutes a system allowing for the solution of {ϕt , c , ω+ } if suitable initial and boundary conditions are given. The set (88) thereby constitutes an extended version of the classical Newman theory (Doyle et al., 1993; Newman and Thomas-Alyea, 2004) with which it coincides for the case of ϕt = X (no deformation) and ± = 0 (no swelling). Under these conditions the last two terms on the right-hand side of the mass and current flux in Eq. (88) drop out. Other than the parameters associated with swelling in Eq. (49) due to the salt molar volume, i.e.  ± and the reference salt concentration cref , no further parameters are required to couple mechanics and ion transport in this formulation. An overview of the relationships of mobilities, transport parameters and binary diffusion coefficients D can be found in f ± (c ) Appendix 8.5. Using the parametrization {D , κ , t+ , 1 + ∂ ln∂ ln }, one has either to check the eigenvalues of the mobility c matrix Mαβ or, in a binary isotropic system, ensure that conductivity, steady state conductivity, diffusivity and the thermodynamic factor are non-negative and the transference number is bounded between zero and one, see Appendix 8.6. Finally, we observe that the osmotic stress in Eq. (56) depends solely on the thermodynamic factor with

pC = 2R



 ∂ ln f± dc. ∂ ln c


6. Examples To highlight features of our different model formulations, we will present a set of numerical studies. We choose a polymer electrolyte with mobile anions and cations and assume that the electrochemical properties remain unchanged while we examine the influence of mechanical material properties on the response of the electrolyte. We utilize the finite element solver COMSOL Multiphysics® . Conservation of mass, the Poisson equation and balance of charge were implemented as general form partial differential equations and discretized with quadratic Lagrange elements. The ”structural mechanics” and ”nonlinear structural materials” module were used for solving the deformation field with cubic Lagrange elements. All equation were solved fully coupled and the field variables involved were normalized to optimize accuracy and convergence. 6.1. Double layer in solid electrolytes We start the section with an example using the full model given by Eq. (75) that treats anions and cations separately and can thus be applied to situations where electroneutrality is violated and where the ions behave differently such as in a double layer. This effect is characterized by local anion/cation polarization and is most pronounced in the case of blocking electrodes. Related studies are found for liquid electrolyte systems (Dreyer et al., 2015; Landstorfer et al., 2016), solid electrolytes with fixed anions (Braun et al., 2015) and references cited therein. Double layer properties can be obtained experimentally e.g. via impedance spectroscopy measurements (Marzantowicz et al., 2008). Material: We will restrict our investigation to the impact of mechanics and thus neglect surface dynamics and electrode response. Material parameters, which are consistent with measurable properties for electrolytes (Edman et al., 20 0 0; ref 3 Singh et al., 2007), are chosen to be κ iso = 0.01S/m, Diso = 5 × 10−12 m2 /s, tiso + = 0.3 and cα = 10 0 0mol/m . We assume an ideal mixture such that f+ = f− = 1. These parameters transform with Eq. (118) to M+ = 1.1 × 10−12 mol 2 , M− = 1.5 ×−12 mol 2 /(mJs ) and M± = 7.6 × 10−13 mol 2 /(mJs ), showing that the anion is more mobile than the cation. The electric permittivity is set to r = 10 (Fontanella and Wintersgill, 1989) and the molar volume of the salt is ± = 140cm3 /mol (Lundgren et al., 2015). Although ion sizes for some non-aqueous solvents are reported (Marcus and Hefter, 2004), there is a lack of such information for solid state electrolytes and therefore we assume − = 10+ . As the influence of mechanics is the focus of this work, we vary Young’s modulus in the range 10MPa to 10 0 0MPa. The lower bound corresponds to a PEO material (Kelly et al., 2016; Maitra et al., 2007) whereas values in the middle of this range are e.g. found in Block-Copolymers (Singh et al., 2007). The host material is assumed to be compressible with a Poisson ratio of ν = 0.4, similar to polystyrene (McKeen, 2009). Boundary value problem (BVP): We consider a blocking electrode (Fig. 5a) of width L = 0.5μm with no-flux conditions for anions and cations on the boundaries. We assume that the right boundary is electrically grounded (r = 0) and the left boundary subjected to l = 10mV . The system is mechanically fixed on both sides. Results: A stationary solution having a double layer is reached within 1ms. As can be seen, electroneutrality holds for most of the domain (Fig. 5b). The potential drops occur very locally at the two electrodes within the double layer (Fig. 5c) which extends to a width of 0.5nm, comparable to the results in Braun et al. (2015) and the Debye length λD = 0.1nm8 Fig. 6a 8

The Debye length represents the length scale of the diffuse layer and is given by λD =


r 0 F2c

(Landstorfer et al., 2016).

M. Ganser, F.E. Hildebrand and M. Kamlah et al. / Journal of the Mechanics and Physics of Solids 125 (2019) 681–713


Fig. 5. (a): Setup of the boundary value problem for blocking electrodes solved with Eq. (75). (b) illustrates the drop of electric potential over the whole domain and (c) zooms into the double layer regions where the electric potential decreases and the chemical potential increases. In systems with a high elastic modulus, mechanics has an influence on the result.

Fig. 6. The figures show the concentration (a) and the pressure profile (b) in the double layer for various material stiffnesses. The higher the elastic modulus, the more pronounced is the pressure gradient which acts as a mechanical driving force and yields a lower concentration of the anions. The impact on the smaller cations is less.

shows the concentration profiles of anions and cations within this double layer. In the region close to the negative electrode the concentration of anions decreases whereas the cations accumulate. The concentration of cations is such that their chemical potential gradient balances the electric field (Fig. 5c). Comparison of the results with different Young’s moduli reveals that the anion concentration depends strongly on the mechanical stiffness: The higher the material stiffness, the more pronounced is the non-symmetric response of the ions (Fig. 6a). The reasons for this is that the large size of the anions leads to swelling that cannot be countered by depletion of the smaller cations. Consequently, a pressure gradient builds up and increases the mechanical contribution to the chemical potential (Fig. 5c). Although both ions are affected by this driving force, the influence on bigger ions outweighs that on smaller ones and results in a considerably lower increase of anion concentration at the electrode for high Young’s modulus. The cation concentration, however, increases slightly as the modulus is increased. This is a result of the lower drop of the anion concentration. The size of the double layer region stays nearly the same when the modulus is varied. Discussion: In summary we see that the suggested full model of Eq. (75) is well suited to describe double layers and can also predict the influence of the mechanical properties of the electrolyte on these electrochemical phenomena. Note that mechanical effects on the double layer are more pronounced in solid electrolytes than in liquid electrolytes, where the solvent is mobile and can relax such constraints. Finally we point out that the electroneutral transport model from Eq. (88) cannot capture the potential drop associated with blocking electrodes, whereas only a full model with dedicated Poisson equation such as Eq. (75) can accurately describe the situation, specifically at the boundaries. However, in contrast to the case of blocking electrodes, charge separation is negligible when flux boundary conditions are used such as for galvanostatic charge and discharge (see Example 6.2). In this situation, the discontinuities in electric potential adjacent to the electrode are features of the redox reactions occurring at the interfaces between the electrolyte and the electrodes. As such they are accounted for by models like Butler-Volmer reaction kinetics or Tafel equations. As a result, the electroneutral model serves as an appropriate approximation when flux BC are invoked. The blocking electrode BC is therefore the condition that most strongly requires an approach that accounts for local charge separation. 6.2. Current through a solid electrolyte When a battery is charged, cations pass from the positive electrode (cathode) to the negative electrode (anode) through the electrolyte. When it occurs at a given current density this process is known as galvanostatic charging. In a binary electrolyte, both anions and cations act as carriers of charge, but only one kind of ion (usually the cation) enters and leaves the system. The other (usually the anion) accumulates against what are effectively blocking electrodes. The cation is locally


M. Ganser, F.E. Hildebrand and M. Kamlah et al. / Journal of the Mechanics and Physics of Solids 125 (2019) 681–713

Fig. 7. The graphs show the converted mobilities (a), drag coefficients (b) and chemical stresses (c) converted from PEO polymer electrolyte (Ma, 1995). The data set is thermodynamically consistent although negative drag coefficients are present. This is a consequence of the negative transference number (shown inset (b)). The chemical stress from Eq. (56) is of the order of 1Mpa and is computed from the thermodynamic factor shown inset.

attracted to the polarizing anions due to strong electrostatic forces (electroneutrality) and eventually a salt concentration gradient builds up. These profiles dictate the performance of a battery and are thus one of the classical tasks for electrochemical simulations to predict. Locally low concentrations are not favorable due to their influence in slowing reaction kinetics on the interfaces, while locally high concentrations can lead to phase separation. The concentration gradient is usually governed by the diffusion coefficient, the transference number and the thermodynamic factor. In this study, we emphasize that mechanical properties can also be key to controlling the salt concentration gradient and thus the performance of the electrolyte. Material: We utilize isotropic parameters of the PEO polymer with NaCF3 SO3 salt (Ma, 1995) at a 85°C operating temperature, see Fig. 7a. This polymer creeps rapidly at elevated temperatures and its mechanical response may be closer to a liquid electrolyte than to a solid. Nevertheless, PEO is also used in composite systems such as PS-PEO with increased mechanical strength as shown in Example 6.3. Beside improving mechanical properties, such compositing strategies generally reduce the volume fraction of the conducting phase, and transport properties are thus degraded. This means that in reality, one cannot improve, e.g., mechanical stiffness without negatively impacting electrochemical transport. In what follows, we disregard this coupling as we are interested in highlighting the effect of mechanics on electrochemistry which is easier to see if electrochemical properties remain unchanged. The mobilities of PEO with NaCF3 SO3 in Fig. 7a were computed with Eq. (118) and show that anions are slightly more mobile than cations through the concentration range. Strong cross interactions between the ions are present due to the high off-diagonal mobility M ± . This is confirmed by the drag coefficient K ± in Fig. 7b obtained with Eq. (114). The negative drag coefficient K−0 is a consequence of the negative transference number. Although the negative sign is non-intuitive, the system is thermodynamically consistent. Fig. 7c shows the osmotic pressure of Eq. (56) as a function of salt concentration. Being up to 2MPa it is small but not negligible. In addition to the electrochemical properties, we use the molar volume ± = 80m3 /mol (Xiao and Trentaine, 1997), Poisson ratio ν = 0.4 and Young’s moduli in the range 10MPa to 500MPa. For the full description as given in Eq. (75), we utilize r = 10, f+ = f− and − = 10+ . These assumptions are irrelevant to the electroneutral model given in Eq. (88), where only the salt properties are required. BVP: We utilize both the full model without electroneutrality from Eq. (75) and the extended Newman model from Eq. (88) to study the complex interplay of electrostatics, chemistry and mechanics. The charging process can be treated as a 1D problem where we assume a separator length of L = 30μm. We apply a current of J app = 0.5mA/cm2 and set  = 0 for the full and ω+ = 0 for the electroneutral model on the right, see Fig. 8. Again, we fix the left and the right boundary mechanically. Results - general: Focusing first on the electrochemical response, we carry out a simulation with fixed Young’s modulus of E = 10MPa. From the results we see that a concentration gradient builds up due to migration of anions in the first 100s, see Fig. 8a. A stationary state is reached when the chemical driving forces on the anions equal the electric driving forces and hence J − = 0. Comparison of the extended Newman model and the full formulation yield the same concentration and electrochemical potential profile within numerical accuracy. In the full formulation, only a small double layer develops with differences in anion and cation concentration of ≤ 0.01mol/m3 . Note that spatial and material quantities are almost identical in this example because mechanical deformation is small. Results - full model: To gain some further understanding, it is instructive to focus on the driving forces within the electrolyte and analyze the drop of the components of the electrochemical potential ωα = μCα + μM α + zα F , i.e.  = right − left throughout the separator. Such a potential drop serves as a first approximation to the internal driving forces α 9 due to grad(ωα ) ≈ ω L . For E = 10MPa, we measure ω− = −8.2mV and hence a relatively high potential drop for the


We use the spatial description to motivate the findings and avoid therefore the formal material anisotropy as discussed in Eq. (63).

M. Ganser, F.E. Hildebrand and M. Kamlah et al. / Journal of the Mechanics and Physics of Solids 125 (2019) 681–713


Fig. 8. Setup of the boundary value problem of galvanostatic charging for (a) the full model, Eq. (75) and (b) the reduced model, Eq. (88). (c): A concentration gradient builds up during charging for E = 10MPa. Both, model (75) and model (88) lead to the same result.

Fig. 9. (a): Steady state concentration gradient which occurs during charging for three different Youngs’moduli. (b): Mechanical pressure pM , osmotic pressure pC and total pressure p within the electrolyte.

anions and ω+ = 13.2mV for the cations. The values can be confirmed via an analytic expression using the steady state assumption (J − = 0) in Eq. (75) which leads to

grad(ω+ ) = −

M− M+ M− − M2±


j+ ,

grad(ω− ) =

M± M+ M− − M2±


j+ .


app /F is constant for a 1D-problem. As soon as off-diagonal mobilities are present, one can The cation mass flux || j SS + || = J observe a drop in the anion electrochemical potential Grad(ω− ) = 0, although J − = 0. To understand the influence of mechanics, we repeat the simulations with E = 100MPa and E = 300MPa. We observe that concentration gradients are less pronounced with increasing mechanical stiffness E, see Fig. 9a. In general, the concentration polarization has two effects on mechanics: (1) The osmotic pressure pC increases sightly for an increase of salt concentration and (2) the swelling mechanism leads to a pressure gradient (Fig. 9b). The mechanical pressure pM thereby shows different characteristics for different moduli. It has a reversed slope with respect to the concentration polarization   for E = 10MPa. This is a consequence of the osmotic pressure pC which couples with the mechanical stress via Div PM + PC = 0. For high moduli E = 300MPa, the picture changes such that the sign of the pressure gradient follows the sign of the concentration gradient because the swelling mechanism becomes dominant. Local volume change is observed to be small with J = 1 ± 0.025. Fig. 10a shows the potential drops with respect to the mechanical moduli E. The electrochemical potential for the cation ω+ is composed of , μC+ and μM + . The chemical contributions for both ions decrease for increasing stiffness and at the same time as there is a less pronounced salt concentration polarization, see Fig. 9a. The mechanical contribution μC− and μC+ is in tandem with the mechanical pressure in Fig. 9b. However, due to the small size of the cation, μM + and therefore its contribution to the electrochemical potential is small. The effect on the anion is substantially more pronounced, see Fig. 10a. There is a small decrease of ω+ , which comes from the concentration dependencies of the mobilities. Integrating the left equation in Eq. (93) we find that

ω + =


grad(ω+ )dv = j +


M− M+ M− − M2±



where the mobilities are concentration dependent. The integral can be interpreted as a macroscopic steady state resistance. Using the concentration profile from Fig. 9a we obtain a smaller value from the integral and therefore a smaller ω+ for an increase of E. Further, we see a small increase of  for increasing E in Fig. 10a. To explain its behavior, we consider the


M. Ganser, F.E. Hildebrand and M. Kamlah et al. / Journal of the Mechanics and Physics of Solids 125 (2019) 681–713

Fig. 10. (a): Potential drops across the electrolyte of the individual contributions for various Young’s moduli. (b): The cation mass flux J + for different stiffnesses at X = 12 L is composed of diffusion, migration, stress driven diffusion and volume change, see Eq. (88).

case of constant mobilities and therefore get from Eq. (62) the estimation


+  const.

∼ grad(ω+ ) = Grad


   μC+ + Grad μM + +F Grad () 



The gradient in electric potential has to compensate the change of chemical potential. Results - Extended Newman formalism: Now, we move our focus to the solution using the extended Newman formalism as given in Eq. (88) and probe the contribution on the mass fluxes in the middle of the separator (X = 12 L). The cation mass flux of −5.2mol/(ms2 ) is composed of diffusion, migration, stress-driven diffusion and a contribution from the effect of volume change, see Fig. 10b. The fractions due to migration and the effect of volume change stay almost constant for the moduli investigated. The stress driven contribution gains importance for higher material stiffness and, therefore, the diffusion contribution is less pronounced. The latter can be seen from the smaller gradient of concentration in Fig. 9a. Note that probing at a fixed location with varying concentration yields slightly different contributions due to the concentration dependent transport properties. Positivity of the migration originates from the negative transference number (Fig. 7b). Discussion: This example illustrates several aspects of the complex mechanisms occurring during galvanostatic charging. Both models, Eqs. (75) and (88), align and no pronounced double layer effects are observed. The parametrization (Ma, 1995) shows strong interactions between anions and cations by means of cross mobilities M ± . They have a direct influence on the steady state current of Eq. (79) and potential drop of Eq. (93), and yield Grad(ω− ) = 0 although J − = 0. Young’s modulus has a great impact on the concentration profile. The decrease of chemical potential drop with increase of Young’s modulus goes together with an increase of electric potential drop. The mechanical contribution μM + is found to be small. Overall, the electrochemical potential drop μ+ is barely influenced by a change of Young’s modulus. The mechanical driving forces in Eq. (88) due to swelling show a pronounced influence on ion transport. Especially for stiffer systems, the influence of stress-driven diffusion is on the order of the diffusion mechanism and must be accounted for. Although we see a strong impact on the concentration distribution in the simulation, the absolute value of mechanical stress is small and therefore not in the regime of inelastic deformation. For soft materials osmotic pressure leads to small strains and stresses which induce (small) transport contributions in the direction opposing the concentration profile. However, for stiffer materials, the importance of osmotic pressure is negligible within the parameter set investigated. Transport effects due to volume change are small in all cases. 6.3. Mechanically deformed electrolyte For our last example, we study a two dimensional problem where the solid electrolyte undergoes an inhomogeneous external deformation. This can be caused for example by inhomogeneous deposition on the metal interface, i.e. the growth of a dendrite intrusion penetrating the separator. The change of effective transport paths is an obvious outcome, but also a stress field in the host will build up which couples with ion transport. With this example, we want to demonstrate the importance of the large deformation formalism for correctly predicting ion transport in a deformed solid electrolyte and we want to show that external deformation has a strong influence on the concentration distribution within a SSE. Setup: We use the full model as given in Eq. (75) to solve this problem. The formulation includes various aspects of coupling between mechanical deformation and electrochemistry. An analysis of the different contributions in the fully coupled system is difficult since all effects interfere with each other. Therefore, we consider in addition to the full model two simplifications to analyze the effects and the importance of the individual coupling mechanisms: Setting A Setting B Setting C

βs = 0 βs = 0 βs = 1

(No swelling) (No swelling) (With swelling)

μC ( c  ) μC ( c ) μC ( c )

(No spatial diffusion) (Spatial diffusion) (Spatial diffusion)

M. Ganser, F.E. Hildebrand and M. Kamlah et al. / Journal of the Mechanics and Physics of Solids 125 (2019) 681–713


Fig. 11. The graphs show the converted mobilities (a), drag coefficients (b) and chemical stresses (c) converted from PS-PEO block-copolymer electrolyte (Timachova et al., 2018). The data set is thermodynamically consistent although negative mobilities and drag coefficients between ions are present. The chemical stress from Eq. (56) is computed from the thermodynamic factor shown inset.

Fig. 12. Boundary conditions for Example 6.3. Mechanical deformation with displacement BC is applied at Phase 3. P1 and P2 mark the probe points for Fig. 13.

Both Setting A and Setting B exclude the swelling mechanism and therefore no contribution of stress driven diffusion is considered. Setting A excludes in addition the contribution of the change of volume to the chemical energy. Hence, the effect in Eq. (72) drops out and no osmotic pressure is present. In the extended Newman model in Eq. (88), Setting A translates p to ξαc = ξα = ξαJe = 0 and only the change of transport path due to deformation is taken into account. Setting C includes all features. Material: For this calculation we use the electrochemical data measured for a PS-PEO block-copolymer (Timachova et al., 2018). We assume a homogeneous material and define the concentration with respect to the overall PS-PEO volume. The resulting mobilities are depicted in Fig. 11a and show a low cation mobility in the concentration range of interest. The mobile anion has a small drag coefficient (Fig. 11b). The negativity of the off-diagonal mobility goes along with a negative drag coefficient between anion and cation. Note that this negative mobility is not related to a negative transference number. Again, these data are admissible from a thermodynamic point of view. The chemical stress is an order of magnitude higher compared to the parametrization of Example 6.2. We assume cref = 1500mol/m3 and that the mechanical properties are governed by the PS-phase with E = 10 0MPa (Singh et al., 20 07; Stone et al., 2012) and ν = 0.4. BVP: We consider a rectangular domain of size 50μm × 30μm with periodic BC left and right (Fig. 12). The top edge is fixed in space and grounded ( = 0). The bottom edge will be deformed into a Gaussian-like shape10 with an amplitude of 5μm and a wavelength of 40μm. Further, we impose a current density Japp on the top and bottom edges. To highlight the results of our theory, we will successively employ six temporal phases: Phase Phase Phase Phase Phase Phase

1 2 3 4 5 6

Current flow Relaxation Deformation Relaxation Current flow Relaxation

t t t t t t

= [0, 250]s = [250, 500]s = [500, 505]s = [505, 750]s = [750, 10 0 0]s = [10 0 0, 1250]s

J app J app J app J app J app J app

= 0.175m A/cm2 =0 =0 =0 = 0.175m A/cm2 =0

undeformed undeformed deformation is imposed deformed deformed deformed

Phase 1 serves as a comparison datum for a current flow on the deformed geometry in Phase 5. The relaxation phases are used to obtain a steady state subsequent to transients. In Phase 3, we apply a deformation much faster than the diffusion 10 We utilize the Wendland function uBC = [X, A(1 − X/R ) )4 (1 + 4X/R )] ∀X ≤ R with amplitude A and radius R to define the protrusion with compact support of size 2R.


M. Ganser, F.E. Hildebrand and M. Kamlah et al. / Journal of the Mechanics and Physics of Solids 125 (2019) 681–713

Fig. 13. Spatial concentration c and electrochemical potential ω+ measured at the tip (P1) and the flank (P2).

Fig. 14. Color contour plots showing the salt concentration and the electrochemical potential in the electrolyte immediately after deformation (end of Phase 3). The spatial concentration shows a peak close to the protrusion tip. The electrochemical potential shows no gradient for Setting A, small effect for Setting B and a strong one for Setting C.

time to highlight the effect of mechanics. Since we obtain via Phase 4 a steady state once more, the features in Phase 5 and Phase 6 are independent of the rate of deformation as long as a purely elastic response is considered. The outcome would coincide with a slowly penetrating protrusion. Results: Fig. 13 shows the evolution of electrochemical potential ω+ and spatial concentration c at two probes attached to the boundary of the material. One probe is associated with the tip of the protrusion (X P1 = [0, 0]) and the other is on the flank (X P2 = [−15μm, 0]), see Fig. 12. During application of current in Phase 1, a concentration gradient builds up. Therefore, the concentration at probes P1 and P2 increases in all settings, see Fig. 13. Since there is no deformation, the results for tip probe P1 and the flank probe P2 coincide. The concentration gradient is not as pronounced for the relatively stiff material and, therefore, Setting C has the least increase of salt concentration. The electrochemical potential is linked to the salt concentration and hence has the same characteristics. Afterwards, in Phase 2, the salt distribution relaxes to the reference concentration cref for all settings. In Phase 3, the volume above the tip (P1) is compressed with Je ≈ 0.8, where swelling effects are small (Je ≈ 1) due to the fast deformation process. Therefore, the same amount of ions is located in a smaller (compressed) volume. Fig. 14 shows contour plots of the salt concentration and electrochemical potential immediately after deformation. If no swelling is ref = 1.25cref . Indeed, this happens in Setting A (Fig. 13a), considered, the salt concentration increases to a value of cα ≈ 01.8 cα α but the electrochemical potential stays constant because of its material description (Fig. 14a). For the other two settings, this value is not reached. The diffusion mechanism due to increased chemical potential counteracts the concentration build-up ref . The from the beginning of the deformation. The peak of salt concentration at P1 in Fig. 13b increases only up to cα ≈ 1.08cα mechanical pressure locally reaches 50MPa and therefore, in Setting C, stress-driven diffusion plays a major role, increases ω+ locally and therefore grad(ω+ ) dramatically, and squeezes ions away from the tip. The snapshot, Fig. 14c, taken directly after deformation, shows depletion in the region of high pressure. An opposite effect can be seen in P2. This region experiences

M. Ganser, F.E. Hildebrand and M. Kamlah et al. / Journal of the Mechanics and Physics of Solids 125 (2019) 681–713


Fig. 15. Color contour plots showing the salt concentration and mechanical pressure in the deformed electrolyte for the relaxed state in equilibrium (end of Phase 4). The concentration in Setting A does not relax, Setting B shows complete relaxation and in Setting C the ions are pushed away from high pressure regions.

Fig. 16. Color contour plots showing the salt concentration and the electrochemical potential in the electrolyte during galvanostatic charging (Phase 5). Setting A shows still a peak at the concentration tip, an almost linear concentration polarization occurs in Setting B and Setting C shows, even during concentration polarization, a depletion of salt at the tip.

a small tension and thus an increase of volume (Je ≈ 1.02). Therefore, the ions occupy more space and the concentration decreases. During relaxation after deformation (Phase 4, Fig. 15), no change of salt concentration occurs in Setting A since the chemical potential defined with respect to the material concentration is ”not aware” of the deformation. Its value stays constant as seen in Fig. 13a and consequently the gradient, which drives transport, is zero. Setting B, in contrast, equilibrates the concentration profile to a value slightly above cref . This is a consequence of the no-flux-boundary conditions and conservation of mass in the deformed computational domain. Since  is prescribed on the boundary, the change of 50mol/m3 affects also the electrochemical potential which increases 5mV. The time constant to reach steady state is comparable to Phase 2. The response in Setting C shows a salt depletion at the tip with c ≈ 0.87cref (Fig. 15c) and therefore goes even below the reference concentration. The high pressure in this region makes the tip an energetically less favorable place for the ions. On the other hand, the region of negative pressure experiences a small increase of concentration. In Phase 5, we solve for galvanostatic charging the deformed geometry (Fig. 16). All of the three settings account for the deformed transport path associated with Eq. (63), i.e. the electrochemical potential at the tip (P1) needed for transport is 6 to 9mV less compared to that at P2 (Fig. 13). A concentration polarization with salt depletion on the tip (P1) and ion accumulation on the flank (P2) originates and, roughly speaking, is superposed on the steady state of Fig. 15. Therefore, the very high concentration in Setting A is observed. Setting B, on the other hand, shows a moderate increase of salt concentration at P1, an effect of the shorter transport path (Fig. 16b). The gradient of concentration in the electrolyte is almost homogeneous and therefore, one finds the highest concentration in the region furthest away from the counter electrode, i.e. in P2. Setting C shows a similar increase of concentration at both probes, P1 and P2, but the squeezing effect is still dominant and the concentration difference between P1 and P2 stays at 300mol/m3 . Nevertheless, due to the additional contribution of pressure in the electrochemical potential, the difference of ω+ at P1 and P2 is smaller compared to the other two settings. Discussion: These examples lay out the complex coupling mechanisms for ion transport in an inhomogeneously deformed solid electrolyte and illustrates the various contributions to the coupling in our model. Each setting (the simplified Setting A and Setting B as well as the full Setting C) with their different assumptions shows not only qualitative, but also quantitative differences in terms of the mechanical and electrochemical response. The missing relaxation mechanism after deformation in Setting A is not compatible with the understanding of mobile ions within an electrolyte. Setting B neglects the property of ion volume which introduces a penalty term for high pressure regions. Only Setting C provides a consistent coupling with considerations of the change of transport path, appropriate relaxation and stress driven diffusion. Therefore, we can predict a concentration and electrochemical potential profile along a deformed separator. Further analysis should incorporate reaction kinetics, i.e. of the Butler-Volmer variety (possibly also including mechanics), to investigate the effect


M. Ganser, F.E. Hildebrand and M. Kamlah et al. / Journal of the Mechanics and Physics of Solids 125 (2019) 681–713

of the electrochemical potential ω+ on the interface on the mass flux through the interface to obtain a full picture of the electro-chemo-mechanical interactions. 7. Conclusion We presented a rigorous formulation of ion transport in a solid electrolyte. Using the concept of finite strain, the Coleman-Noll procedure and the Helmholtz energy, we derived in a consistent manner the driving forces for species transport, the mechanical response of the host material and their coupling mechanism. We combined therein the transport description in the deformed geometry but solved it —as common in continuum mechanics— in a reference description. The influence of mechanics is broad. Most notable is the consideration of ion size by means of swelling and eventually stress driven diffusion. Further, we considered the effect of the change of geometry on mass flux and in the chemical potential. The latter gives rise to an osmotic pressure and is closely related to the thermodynamic factor. Using electroneutrality assumptions for a binary electrolyte, we derived an extended version of concentrated solution theory (Newman and Thomas-Alyea, 2004) taking diverse mechanical coupling mechanisms into account. We identified the connections among mobilities, drag coefficients and transport parameters by means of conductivity, diffusivity, transference number and thermodynamic factor from concentrated solution theory. Further, the relation to the steady state current was shown and thermodynamic consistency was discussed. The off diagonal terms of the mobility matrix and hence the anioncation interaction were found to be significant. Several examples illustrated the impact of mechanics on ion transport by using both a transport model with full considerations of electrostatics and one with the assumption of electroneutrality. Only the first is capable of resolving the double layer at a blocking electrode. In the case of galvanostatic charging, the electroneutrality condition holds throughout the domain and both models align. Maxwell stresses therefore are negligible and also the osmotic pressure was found to be small for relatively stiff materials. Using electrochemical transport parameters from PEO based systems and extending them by mechanical stiffness gave us the possibility of investigating the impact of Young’s moduli on double layers (in blocking electrodes) as well as the salt concentration and electrochemical potential (galvanostatic charging). The salt concentration profile is strongly affected by the stiffness of the host material although only small pressure gradients occur. The electrochemical potential, however, alters only little. Finally we studied the response of a deformed electrolyte. Reasonable simulation results were only obtained with comprehensive mechanical considerations. External deformation thereby influences the electrochemistry significantly. The concentration profile along the interface varies strongly and suggests that a depletion of salt occurs in the region of highest intrusion. This understanding should be considered in studies of dendrite formation. Acknowledgments The results presented were partly achieved during a visit of four months in the group of Prof. Robert McMeeking at the University of California, Santa Barbara. This stay was financially supported by a scholarship of the Karlsruhe House of Young Scientists (KHYS). Markus Ganser and Felix Hildebrand are very grateful to Ulrich Sauter and Matthias Hanauer for their valuable comments. Appendix The aim of this section is to provide the reader with some background of finite deformation continuum mechanics, multicomponent transport by means of the Maxwell-Stefan theory, its relation to our derivation and conversion of transport parameters in binary electrolyte systems. We further discuss requirements of the transport parameters for thermodynamic consistency and the definition of the steady state current in a binary electrolyte. 8.1. Continuum mechanics formalities We utilize Reynold’s theorem and the framework of finite deformation continuum mechanics (Holzapfel, 20 0 0; Truesdell and Noll, 2004) for a mathematically consistent description of multicomponent flow within a deforming host material. 8.1.1. Reynold’s theorem Let Ptα ⊆ Bt be a control volume occupied by species α moving at velocity vα . Postulating that Ptα coincides at time t with a control volume Pt associated to a host material moving at velocity v0 = ϕ˙ t , volume and surface integrals coincide at time t:


y d vα =


yd v


∂ Ptα

y˜ · ndaα =

∂ Pt

y˜ · nda,


where y is a scalar or vector field defined on the domain Pt and y˜ a vector field defined on the boundary ∂ Pt , respectively. However, when considering the rates, we have to take the different time-dependencies of the integration domains Ptα and Pt into account. To describe the relative motion as seen in ❷ of Fig. 17, we use Reynold’s theorem

d dt


y d vα =

d dt


yd v +

∂ Pt

y ( (vα − v0 ) · n )da

M. Ganser, F.E. Hildebrand and M. Kamlah et al. / Journal of the Mechanics and Physics of Solids 125 (2019) 681–713


Fig. 17. Transformation steps: ❶ Host domain Pt and species domain Ptα underly a different motion during a timestep δ t. ❷ With the Reynold’s theorem, one considers the change of concentration within the host material with a mass flux jα . ❸ All quantities can be expressed with respect to a material, undeformed configuration P0 .

Fig. 18. Anisotropic material is defined in the intermediate configuration M∗αβ is defined on an intermediate configuration, where the orientation is linked to the material frame but the stretches and hence the length of the effective transport paths to the spatial frame.

d dt



 yd v +

1 ∂ Pt cα

y ( j α · n )da,


where we have introduced the mass flux j α = cα (vα − v0 ) with concentration cα . The second term on the right hand side denotes the influx of volume quantity y into the host domain due to a relative velocity. 8.1.2. Finite strain transformation Each spatial point in Pt can be linked to a reference point in P0 ⊆ B0 . These material points in P0 are usually associated with an undeformed material. The deformation gradient F describes the mapping of tangential and normal vectors via t = FT and n = F−T N . The integrals in the material and spatial configuration coincide for properly transformed quantities (see ❸ in Fig. 17) and we obtain


yd v =

Y dV


∂ Pt


y˜ · nda =

∂ P0

Y˜ · N dA,


where the corresponding material properties are Y = Jy and Y˜ = Jy˜ F−T . The latter is known as Nanson’s relationship. Mass fluxes transform therefore with J α = JF−1 j α . Since P0 is a reference description, it is independent of time and thus the integral operator and time derivatives are interchangeable

d dt


Y dV =


d Y dV = dt

Y˙ dV.



8.1.3. Reynolds theorem in finite strains Combining Eq. (98) with Eq. (97) and using

y cα


reference description of matter via

d dt

y d vα =



∂ P0

1 Y (J α · N )dA. cα


Y˙ dV +

Jy Jc α


Y  cα

gives us the transformation to link moving species within a


    The divergence theorem ∂ P Y · N dA = P Div Y˜ dV yields a fully volumetric description with 0 0 d dt


y d vα =


Y˙ + Div

Y  J α dV. 


With the chain rule, we can further resolve the divergence to

d dt


 Pα t

y d vα =


Y˙ + Grad


Y Jα +

1 Y Div(J α )dV. cα



M. Ganser, F.E. Hildebrand and M. Kamlah et al. / Journal of the Mechanics and Physics of Solids 125 (2019) 681–713

Using balance of mass with Eq. (13) yields eventually

d dt

y d vα =



Y˙ + Grad


Rα c˙  Y J α +  Y − α Y dV. cα cα cα


Eq. (103) gives a straightforward expression to write a time dependent species integral Ptα in material coordinates. To do so, convective flux, reaction and change of mass terms have to be taken into account. 8.2. Balance of angular momentum The balance of angular momentum states that the time derivative of the angular momentum L is equal to the total torque M:

d L=M dt





α =0 Ptα x × ρα vα dvα 


ECM + oEC dvα + ∂ Pt x × t ECM da α =0 Ptα x × b


 ECM M E C with the body force bECM = N = t M + t E + t C . Further, an angular momentum due α =0 bα + bα + bα and traction force t to dipoles has to be considered. Following (McMeeking and Landis, 2005), we use oE = oEi = i jk σkEj . Although a deviatoric contribution of the chemical stress is unlikely, we set oC = i jk σkCj , and yields with the electrostatic contribution oEC = oE + oC . The rate of angular momentum computes to

1  N   d d L= (X × ρα Vα ) + Div  (X × ρα Vα )  J α dV dt cα P0 dt α =0




α =0 P0

i jk X j


  d   α ρα Vk + Divm Mα Vkα Jm + Mα i jk δ jmVkα Jmα dV dt


and the applied torque is

M= =


α =0 Pt N  

α =0 P0

i jk x j bECM dv + k

 ∂ Pt

x × (σ M n )da +


i jk σkECj dv

  M  i jk X j BECM + Divm Pkm + i jk δm j PklM Fml + i jk PklEC Fjl dV. k

Combining Eqs. (105) and (106) gives



α =0 P0


i jk X j


α =0 P0



   M d   α − Divm Pkm dV ρ V + Divm Mα Vkα Jm − BECM k dt α k

  i jk PklECM Fjl − Mα Vkα Jαj dV.


The left hand side is equal to  zero due to balance of linear momentum of Eq. (20) and BECM = BM + Div PC + PE . The right hand side results in the symmetry condition



α =0

Mα J α  Vα = FPT −


α =0

Mα Vα  J α .


Under the assumption stated over Eq. (43), this reduces to PFT = FPT . 8.3. Anisotropic mobilities In order to include anisotropic material properties as found e.g. in laminar structures (Golodnitsky et al., 2015; Timachova et al., 2018) in our transport framework, we link transport paths with the material configuration and use properties like the effective length of a transport path from the deformed configuration. For illustration purposes, imagine a set of parallel lines on a material which are associated with the material points. Their shape and length will change with deformation. Now consider an ion traveling through the material without being able to cross the lines. The ion has then to follow the material structure (the lines are borders) but needs to move the spatial length of these lines (deformed geometry). To realize this, we decompose the deformation gradient into F = RU with a right stretch tensor U and a rotation tensor R.11 This defines an intermediate configuration, where the material was stretched with U but not rotated. It serves for the definition of the 11

The orthogonal rotation tensor R has the properties det(R ) = 1 and RT = R−1 and the symmetric right stretch tensor is U =

√ FT F (Holzapfel, 20 0 0).

M. Ganser, F.E. Hildebrand and M. Kamlah et al. / Journal of the Mechanics and Physics of Solids 125 (2019) 681–713


transport properties, e.g. anisotropic mobility matrix M∗αβ . If the structure rotates, the effective spatial mobility Mαβ does so via

Mαβ = RM∗αβ RT .


We get with Eq. (63) the material mobility matrix as

Mαβ = U −1 M∗αβ U −T .


The material mobility matrix is ultimately influenced only by stretches. The mass flux stated in the material frame takes length changes of the transport path into account, the orientation, however, is governed by a path associated with the material. If an electrolyte is isotropic, we can reduce the formulation by the introduction of an isotropic mobility M∗αβ = Mαβ I , which results in

and Mαβ = J FT F

Mαβ = Mαβ RRT = Mαβ I




Therefore, the spatial mobility matrix coincides with the intermediate definition. Nevertheless, the material mobility behaves anisotropic due to deformation. This anisotropy is required because of the nature of finite deformation where the partial differential equations are solved in the undeformed configurations. 8.4. Link to Maxwell-Stefan theory From Eq. (62), we can directly deduce the Maxwell-Stefan diffusion model (Krishna and Wesselingh, 1997), which is also the foundation of Newman’s derivation (Kim and Srinivasan, 2016; Monroe and Delacourt, 2013; Newman and ThomasAlyea, 2004). For simplicity, we restrict ourself to isotropic mobilities Mαβ = Mαβ I . Taking the inverse of the mobility matrix and multiplying with cβ on both sides of Eq. (62) gives

cβ grad

N N      ˜ αβ (vα − v0 ) = K Kαβ vα − vβ , ωβ =

α =1

with the definition

˜ αβ = −cβ M−1 cα K αβ


α =0


M−1 αβ =

M11 M21 .. .

M12 M22 .. .

··· ··· .. .

⎞−1   ⎠  .  



The transformation between the two drag coefficient matrices is

˜ αβ = K


− γ =α =0 Kαγ Kαβ N +1×N +1

Note that Kαβ ∈ R N×N

and Mαβ ∈ R

α=β α = β

and Kαβ =

⎧ 0 ⎪ ⎨− N ⎪ ⎩−

˜ N ˜ γ =0=0 Kγ β

Kαγ γ =0=0

˜ αβ K

α=β α = β = 0 β = α = 0



˜ αβ ∈ RN×N considers friction with the host but no self or background friction. The matrices K

have non-zero diagonal elements. These diagonal elements correspond to interaction with the underlying

˜ −1 media. All reformulations maintain symmetry. The transformation of drag coefficients to mobilities is Mαβ = cα cβ K αβ . The Maxwell-Stefan diffusion coefficients Dαβ are related to the second drag coefficients with

Kαβ = R

cα cβ cT Dαβ



where cT is the total concentration. Further, a host concentration c0 is necessary within this description. For dilute solutions cT c with no component-component interaction (Mαβ = 0 ∀α = β ), we obtain the correlation Mββ = R cβ D0β ∀β = 1..N. The 0

Maxwell-Stefan diffusion coefficients are similar to self-diffusion coefficients obtained by experiments (Kim and Srinivasan, 2016; Timachova et al., 2018) or MD simulations (Borodin et al., 2006). If cross dependencies are present, the conversion within the concepts of mobilities and Maxwell Stefan coefficients is not as trivial and we refer to Appendix 8.5 for a binary electrolyte. 8.5. Conversion of transport parameters, mobilities and drag coefficients Binary electrochemical systems are usually characterized by macroscopic measurements yielding isotropic transport pa∂ ln f± (c ) iso rameter {κ iso , tiso + ,D ,1 + ∂ ln c }. They are linked to mobilities with Eqs. (67), (68), and (89) by

κ iso = F 2 (M+ + M− − 2M± ),


M. Ganser, F.E. Hildebrand and M. Kamlah et al. / Journal of the Mechanics and Physics of Solids 125 (2019) 681–713

tiso + = Diso =

M+ − M± M+ + M− − 2M± M+ M− − M2± 2R ξ±c M+ + M− − 2M± c

 ∂ ln f± (c ) . ∂ ln c



The spatial properties transform into the material description with

 −1 κ = J FT F κ iso ,

D = J FT F

t+ = tiso + I,


Diso .


Taking the inverse of Eq. (116), we deduce mobilities from measured parameters as

Mγ = M± =

Diso c


2 R ξ± 1 +

 + κ iso ∂ ln f± (c ) ∂ ln c

Diso c

2 R ξ

c ±

1 + ∂ ln∂ lnf±c(c )

 − κ iso

tiso γ



γ = {+, −}

iso tiso + t− . F2


The first part in Eq. (118) describes the influence of diffusion. The second equation reveals the diffusion and migration contribution onto the mobility. The ratio of both for the anion is given by the Newman number, see Eq. (126). Further, if we assume f± = f+ = f− , we can compute from the fugacity the thermodynamic factor by integration as

f± (c∗ ) = exp

c∗ 0

 ∂ ln f± (c ) 1 dc . ∂ ln c c


Eqs. (118) and (119) therefore yield the transformation rules to convert transport properties used in formulation (88) to the properties used in Eq. (75). Also drag coefficients are used to describe transport. Using Eq. (114) with the inverse of Eq. (113) gives the correlation between binary drag and diffusion coefficients for a binary system as

M+ M±

M± M−


1 K0− + K± K± c2 (K0+ K0− + K± (K0+ + K0− ) )

K± . K0+ + K±


The on-diagonal mobilities Mαα stem from a complex interplay between friction of host as well as the adjacent component. If the drag coefficients are positive then the mobilities are also positive. Combining Eqs. (116), (120), and (115) provides us with

F 2 (K0+ + K0− ) (K0+ K0− + K± (K0+ + K0− ) ) K0 − D0+ = = K0+ + K0− D0+ + D0−

κ iso = tiso +

Miso amb =


c2 cT c D0+ D0− cT c = =: D K0+ + K0− R c0 D0+ + D0− R c 0


where we have linked the isotropic ambipolar mobility Miso amb from Eq. (82) and the thermodynamic diffusivity D (Kim and Srinivasan, 2016; Newman and Thomas-Alyea, 2004). It is worth noting that only the conductivity depends on the component-component friction K ± . Further, the transference number can only be negative if one of the drag coefficients K0− or K0+ are negative. It is important to note, that measurements of electrochemical transport parameters do not usually take the mechanical response into account. As seen by the derivation, conductivity and transference number are not influenced by mechanics. The steady state conductivity from Eq. (80) is another measurable parameter without direct impact of mechanics. The diffusion coefficient is, however, usually identified by experiments where quite large concentration gradients occur. In such measurements, we expect already a strong influence of mechanical driving forces which assist the diffusion. The measured relaxation time is usually only associated to an effective diffusion coefficient. Atomistic simulations could resolve this problem. Further, we want to point out that a third component would introduce three new transport parameters due to the nondiagonal entries. In general, having N species requires N (N + 1 )/2 + 1 independent measurable parameters to account for all possible interactions between ions and the chemical potential. 8.6. Thermodynamically consistent transport parameters We have seen that the mobility matrix Mαβ has to be a positive semi-definite form for thermodynamic consistency in Eq. (61). The eigenvalues, therefore, can be checked for validation. However, assuming a binary isotropic electrolyte, one can see consistency more easily. From Eq. (82) we deduce that the determinant is positive for positive conductivity, diffusivity

M. Ganser, F.E. Hildebrand and M. Kamlah et al. / Journal of the Mechanics and Physics of Solids 125 (2019) 681–713


and thermodynamic factor:

  c det Mαβ = Diso κ iso 2 R ξ c


∂ ln f± (c ) ∂ ln c


−1 .


Note that a positive determinant does not necessarily yield a positive definite matrix. However, we know that the diagonal elements are positive due to positive steady state conductivity in Eq. (80). If the symmetric mobility matrix has a strongly diagonally dominant structure, it has a positive definite form. This feature is given with 0 < tiso + < 1 as seen in Eq. (116). Therefore, we can conclude that if

0 < κ iso ,


κ iso


0 < Diso ,


0 < tiso + < 1,


∂ ln f (c ) , ∂ ln c


the second law of thermodynamics is automatically fulfilled. If not, the eigenvalues of Mαβ have to be computed and checked for non-negativity. Therefore, the transference number tiso + cannot necessarily take any real value as proposed in Monroe and Delacourt (2013).

8.7. Relation to steady state current in a binary electrolyte Finally, we want to discuss the steady state current in more detail. Other than in Eq. (79), we derive the steady state with the model of Eq. (84) and follow the spirit of Balsara and Newman (2015). In the case of steady state, the diffusive and electric mass flux of the anions cancel out and we get

J − = −Mamb Grad(μ+ + μ− ) −

t− SS ! J = 0. F q


This gives a correlation between the steady state current J SS q and the chemical potential. Rearrangement yields

J SS q =

1 κ −t− Grad(μ+ + μ− ) , Ne F


where we have introduced the so called Newman number

Ne =

iso 2 (tiso ) − κ

F 4 det Mαβ


(M− − M± )2 det Mαβ


= ξc

κ iso tiso ∂μ± − . F 2 Diso ∂ c


Here we have used Eq. (82) for Mamb and assumed isotropic transport properties. Eq. (126) is in alignment with Balsara and Newman (2015). The Newman number is non-negative due to the thermodynamic constraint det Mαβ ≥ 0. The right hand side of Eq. (125) corresponds to the diffusion part of the electric current in Eq. (84). Consequently, we obtain  SS J SS q = − (κ ) Grad

ω  +



(κ )SS =

κ 1 + Ne


with the steady state conductivity (κ )SS . This is consistent with the result of in Eq. (79). The total current can then be computed via

 SS Jq


1  E SS J . 1 + Ne q


In the case of a single ion conductor with t− = 0 this yields Ne = 0 and the electric current drives the whole mass transport. Another limiting case is when no cross relations between the ions are present. Then, the right hand side of Eq. (128) reduces to the positive transference number t+ and is in alignment with Balsara and Newman (2015) and Bruce and Vincent (1987). Separation of ionic current in chemical and electric contribution in Eq. (66) gives

 C SS Jq


Ne  E SS J 1 + Ne q


which shows with Ne ≥ 0 that the chemical part of the ionic current always flows in the opposite direction of the electrical induced current and therefore reduce the effective transport. Note that no assumptions to electroneutrality or dilute limit (e.g. fα = 1) are required for this result. The Newman number was recently used to characterize transport properties, especially the transference number (Pesko et al., 2017).


M. Ganser, F.E. Hildebrand and M. Kamlah et al. / Journal of the Mechanics and Physics of Solids 125 (2019) 681–713

8.8. Symbols

Table 1 Symbols [material frame; spatial frame; unit; description]. Geometric properties X x Grad( · ) grad( · ) N n dV dV dA da Pt P0 Bt B0 Fields cα cα

ρα  ρtot ρqα ρq

ρα ρ tot ρqα ρq

Jα J Eα , J Cα Jtot Jq J Eq , J Cq Rα Rqα D E P

jα j tot jq rα rqα d e p


RH JH Eint , Ekin LH , LM , LE eα

rH jH

ψα ηα ψ˜ ψ˜ M , ψ˜ E , ψ˜ αC μα  ωα δV

m 1/m 1 m3 m2

Spatial coordinates Spatial derivative Normal vector Volume element Area element Arbitrary control volume Body of the host material

mol/m3 kg/m3 kg/m3 As/m3 As/m3 mol/m2 s mol/(m2 s) mol/(m2 s) A/(mol2 ) A/mol2 mol/(m3 s) As/(m3 s) As/m2 V/m As/m2 As/m2 J/m3 s J/m3 J J/s J/kg J/kg JK/kg J/m3 J/m3 J/mol V J/mol J/m3 s

Concentration Partial mass density Total mass density Partial charge density Charge density Total molar flux of species α Molar flux due to electric and chemical potentials Overall molar flux Total ionic current Ionic current due to electric and chemical potentials Concentration sink or source Charge source Electric displacement Electric field Electric polarization Surface charge Heat source Heat flux Internal and kinetic energy in Ptα Power supply due to heat, mechanics and electrostatic Specific internal energy of species α Specific Helmholtz energy of species α Specific entropy of species α Total Helmholtz energy Mechanical and electrical and chemical Helmholtz energy Total chemical potential Electric potential Electrochemical potential Dissipation density

Table 2 Symbols [material frame; spatial frame; unit; description]. Mechanics:


u F J Fe , Fs Je , Js Vα V0 E C BM α , Bα , Bα BVα ˜α B

TM , TE , TC PM PE PC P . . Parameter:

α ±

E C bM α , bα , bα

tM , tE , tC

σC σ p pC , pM

m m 1 1 1 1 m/s m/s N/m3 N/m3 N/m3 N/m2 N/m2 N/m2 N/m2 N/m2 N/m2 N/m2

Mapping from material to spatial frame Displacement vector Deformation gradient Jacobian of deformation gradient Elastic and swelling part of J Jacobian of elastic deformation gradient Velocity of species α (same for spatial and material frame) Host velocity Mechanical, electrical, and chemical body forces Body force related to component inertia Extended specific body force Mechanical, electrical and chemical traction Mechanical Piola stress tensor, see Eq. 34 Maxwell stress tensor Stress tensor due to chemical forces Total Piola (material) and Cauchy (spatial) stress tensor Hydrostatic pressure Osmotic pressure and pressure due to mechanical stresses

m3 /mol m3 /mol

Molar volume of species α Molar volume of salt (continued on next page)

M. Ganser, F.E. Hildebrand and M. Kamlah et al. / Journal of the Mechanics and Physics of Solids 125 (2019) 681–713


Table 2 (continued)

βs zα


cα Mαβ


 r ref

M+ , M− , M± Mamb ˜ αβ Kαβ , K

M+ , M− , M ±

fα f ± (c ) 1+ ∂ ln∂ ln c

κ (κ )SS κα



tα D

γNH , λNH

1 1 1 m3 /mol mol2 /(mJs)

tα D

E, G

Jm 1 1 1/(m) 1/(m) mol/(V s m) 1 m2 /s MPa 1


Factor to account for non-ideality in swelling Charge number Relative electric permittivity Reference concentration Mobility matrix between species α and β Mobility of cation, anion and cross mobility Ambipolar mobility matrix Drag coefficient Fugacity or activity coefficient Thermodynamic factor Ionic conductivity Steady state conductivity Mass conductivity Transference number of species α Diffusion coefficient lamé constants Young’s modulus, shear modulus Poisson ratio

Table 3 Symbols [material frame; spatial frame; unit; description]. Other: t T

dNα I1 , I2 , I3 ξαc , ξ p , ξαJe I Constants F

Value 964,885 0 8.854 · 10−12 R 8.314 Boundary conditions Japp uBC

s s K 1 1 1 1

Time Time domain Temperature Amount of species α in volume dV Invariants Prefactors defined in Eq. (73) Identity matrix

As/mol F/m J/K mol

Faraday constant Electric permittivity of vacuum Universal gas constant

A/mol2 m

Current density applied on interface Prescribed boundary displacement

References Atkins, P.W., Paula, J.D., Walters, V., 2006. Physical Chemistry, 8th W H Freeman, New York. OCLC: 153577665. Baker, D.R., Verbrugge, M.W., Bower, A.F., 2016. Thermodynamics, stress, and stefan-maxwell diffusion in solids: application to small-strain materials used in commercial lithium-ion batteries. J. Solid State Electrochem. 20 (1), 163–181. doi:10.1007/s10008-015-3012-7. Balsara, N.P., Newman, J., 2015. Relationship between steady-state current in symmetric cells and transference number of electrolytes comprising univalent and multivalent ions. J. Electrochem. Soc. 162 (14), A2720–A2722. doi:10.1149/2.0651514jes. Bazant, M.Z., 2013. Theory of chemical kinetics and charge transfer based on nonequilibrium thermodynamics. Acc. Chem. Res. 46 (5), 1144–1160. doi:10. 1021/ar300145c. Biot, M., 1977. Variational lagrangian-thermodynamics of nonisothermal finite strain mechanics of porous solids and thermomolecular diffusion. Int. J. Solids Struct. 13 (6), 579–597. doi:10.1016/0 020-7683(77)90 031-2. Bohn, E., Eckl, T., Kamlah, M., McMeeking, R., 2013. A model for lithium diffusion and stress generation in an intercalation storage particle with phase change. J. Electrochem. Soc. 160 (10), A1638–A1652. doi:10.1149/2.011310jes. Borodin, O., Smith, G.D., Henderson, W., 2006. Li + cation environment, transport, and mechanical properties of the liTFSI doped N -methyl- N alkylpyrrolidinium + TFSI - ionic liquids. J. Phys. Chem. B 110 (34), 16879–16886. doi:10.1021/jp061930t. Bower, A., Guduru, P., Sethuraman, V., 2011. A finite strain model of stress, diffusion, plastic flow, and electrochemical reactions in a lithium-ion half-cell. J. Mech. Phys. Solids 59 (4), 804–828. doi:10.1016/j.jmps.2011.01.003. http://linkinghub.elsevier.com/retrieve/pii/S0 0225096110 0 0 044. Braun, S., Yada, C., Latz, A., 2015. Thermodynamically consistent model for space-charge-layer formation in a solid electrolyte. J. Phys. Chem. C 119 (39), 22281–22288. doi:10.1021/acs.jpcc.5b02679. Brissot, C., Rosso, M., Chazalviel, J.N., Lascaud, S., 1999. Dendritic growth mechanisms in lithium/polymer cells. J. Power Sources 81–82, 925–929. doi:10. 1016/S0378- 7753(98)00242- 0. Bruce, P.G., Vincent, C.A., 1987. Steady state current flow in solid binary electrolyte cells. J. Electroanal. Chem. Interfacial Electrochem. 225 (1–2), 1–17. doi:10.1016/0 022-0728(87)80 0 01-3. Bucci, G., Chiang, Y.M., Carter, W.C., 2016. Formulation of the coupled electrochemical mechanical boundary value problem, with applications to transport of multiple charged species. Acta Mater. 104, 33–51. doi:10.1016/j.actamat.2015.11.030. http://linkinghub.elsevier.com/retrieve/pii/S1359645415300811. Christensen, J., Newman, J., 2006. Stress generation and fracture in lithium insertion materials. J. Solid State Electrochem. 10 (5), 293–319. doi:10.1007/ s10 0 08-0 06-0 095-1. Coleman, B.D., Noll, W., 1963. The thermodynamics of elastic materials with heat conduction and viscosity. Arch. Ration. Mech. Anal. 13 (1), 167–178. doi:10.1007/BF01262690. Cui, Z., Gao, F., Qu, J., 2012. A finite deformation stress-dependent chemical potential and its applications to lithium ion batteries. J. Mech. Phys. Solids 60 (7), 1280–1295. doi:10.1016/j.jmps.2012.03.008. http://linkinghub.elsevier.com/retrieve/pii/S0 0225096120 0 0622. Doyle, M., Fuller, T.F., Newman, J., 1993. Modeling of galvanostatic charge and discharge of the lithium/polymer/insertion cell. J. Electrochem. Soc. 140 (6), 1526. doi:10.1149/1.2221597.


M. Ganser, F.E. Hildebrand and M. Kamlah et al. / Journal of the Mechanics and Physics of Solids 125 (2019) 681–713

Doyle, M., Newman, J., 1996. Comparison of modeling predictions with experimental data from plastic lithium ion cells. J. Electrochem. Soc. 143 (6), 1890. doi:10.1149/1.1836921. Dreyer, W., Guhlke, C., Mueller, R., 2015. Modeling of electrochemical double layers in thermodynamic non-equilibrium. Phys. Chem. Chem. Phys 17 (40), 27176–27194. doi:10.1039/C5CP03836G. Edman, L., Doeff, M.M., Ferry, A., Kerr, J., De Jonghe, L.C., 20 0 0. Transport properties of the solid polymer electrolyte system p(EO n liTFSI. J. Phys. Chem. B 104 (15), 3476–3480. doi:10.1021/jp993897z. Ehlers, W., Karajan, N., Markert, B., 2009. An extended biphasic model for charged hydrated tissues with application to the intervertebral disc. Biomech. Model. Mechanobiol. 8 (3), 233–251. doi:10.1007/s10237- 008- 0129- y. Enikov, E. T., Seo, G. S., 2002. Large deformation model of ion-exchange actuators using electrochemical potentials. Pp. 199–209 doi:10.1117/12.475165. http://proceedings.spiedigitallibrary.org/proceeding.aspx?articleid=882864. Enikov, E.T., Seo, G.S., 2005. Analysis of water and proton fluxes in ion-exchange polymer-metal composite (IPMC actuators subjected to large external potentials. Sens. Actuators, A 122 (2), 264–272. doi:10.1016/j.sna.2005.02.042. http://linkinghub.elsevier.com/retrieve/pii/S09244247050 0150 0. Evans, J., Vincent, C.A., Bruce, P.G., 1987. Electrochemical measurement of transference numbers in polymer electrolytes. Polymer 28 (13), 2324–2328. doi:10.1016/0032- 3861(87)90394- 6. Ferguson, T.R., Bazant, M.Z., 2012. Non-equilibrium thermodynamics of porous electrodes. J. Electrochem. Soc. 159 (12), A1967–A1985. Ferrese, A., Newman, J., 2014. Mechanical deformation of a Lithium-Metal anode due to a very stiff separator. J. Electrochem. Soc. 161 (9), A1350–A1359. doi:10.1149/2.0911409jes. Fontanella, J.J., Wintersgill, M.C., 1989. Low frequency dielectric properties of polyether Electrolytes. In: Polymer Electrolyte Reviews 2. ELSEVIER. http: //www.dtic.mil/dtic/tr/fulltext/u2/a196494.pdf. Ganser, M., Hildebrand, F.E., Klinsmann, M., Hanauer, M., Kamlah, M., McMeeking, R.M., 2019. An Extended Formulation of Butler-Volmer Electro-Chemical. React. Kinet. Includ. Influence. Mech. Submitted to JES. Golodnitsky, D., Strauss, E., Peled, E., Greenbaum, S., 2015. Review-on order and disorder in polymer electrolytes. J. Electrochem. Soc. 162 (14), A2551–A2566. doi:10.1149/2.0161514jes. Goodenough, J.B., Singh, P., 2015. Review—solid electrolytes in rechargeable electrochemical cells. J. Electrochem. Soc. 162 (14), A2387–A2392. doi:10.1149/ 2.0021514jes. Gouverneur, M., Schmidt, F., Schönhoff, M., 2018. Negative effective li transference numbers in li salt/ionic liquid mixtures: does li drift in the ”Wrong” direction? PCCP 20 (11), 7470–7478. doi:10.1039/C7CP08580J. Goyal, P., Monroe, C.W., 2017. New foundations of newman\’s theory for solid electrolytes: thermodynamics and transient balances. J. Electrochem. Soc. 164 (11), E3647–E3660. doi:10.1149/2.0611711jes. Harry, K.J., Higa, K., Srinivasan, V., Balsara, N.P., 2016. Influence of electrolyte modulus on the local current density at a dendrite tip on a lithium metal electrode. J. Electrochem. Soc. 163 (10), A2216–A2224. doi:10.1149/2.0191610jes. Holzapfel, G.A., 20 0 0. Nonlinear solid mechanics: A continuum approach for engineering. Wiley, Chichester, New York. Hong, W., Wang, X., 2013. A phase-field model for systems with coupled large deformation and mass transport. J. Mech. Phys. Solids 61 (6), 1281–1294. doi:10.1016/j.jmps.2013.03.001. http://linkinghub.elsevier.com/retrieve/pii/S0 0225096130 0 0550. Hong, W., Zhao, X., Zhou, J., Suo, Z., 2008. A theory of coupled diffusion and large deformation in polymeric gels. J. Mech. Phys. Solids 56 (5), 1779–1793. doi:10.1016/j.jmps.2007.11.010. http://linkinghub.elsevier.com/retrieve/pii/S0 0225096070 02244. Huttin, M., Kamlah, M., 2012. Phase-field modeling of stress generation in electrode particles of lithium ion batteries. Appl. Phys. Lett. 101 (13), 133902. doi:10.1063/1.4754705. Kelly, T., Ghadi, B.M., Berg, S., Ardebili, H., 2016. In situ study of strain-dependent ion conductivity of stretchable polyethylene oxide electrolyte. Sci. Rep. 6, 20128. doi:10.1038/srep20128. Khurana, R., Schaefer, J.L., Archer, L.A., Coates, G.W., 2014. Suppression of lithium dendrite growth using cross-linked polyethylene/poly(ethylene oxide) electrolytes: a new approach for practical lithium-metal polymer batteries. J. Am. Chem. Soc. 136 (20), 7395–7402. doi:10.1021/ja502133j. Kim, S.U., Srinivasan, V., 2016. A method for estimating transport properties of concentrated electrolytes from self-diffusion data. J. Electrochem. Soc 163 (14), A2977–A2980. doi:10.1149/2.0541614jes. Klinsmann, M., Rosato, D., Kamlah, M., McMeeking, R.M., 2016. Modeling crack growth during li extraction in storage particles using a fracture phase field approach. J. Electrochem. Soc. 163 (2), A102–A118. doi:10.1149/2.0281602jes. http://jes.ecsdl.org/content/163/2/A102. Kreuer, K.D., Paddison, S.J., Spohr, E., Schuster, M., 2004. Transport in proton conductors for fuel-cell applications: simulations, elementary reactions, and phenomenology. Chem. Rev. 104 (10), 4637–4678. doi:10.1021/cr020715f. Krishna, R., Wesselingh, J., 1997. The maxwell-stefan approach to mass transfer. Chem. Eng. Sci. 52 (6), 861–911. doi:10.1016/S0 0 09-2509(96)0 0458-7. Lai, W., Ciucci, F., 2011. Mathematical modeling of porous battery electrodes: revisit of newman’s model. Electrochim. Acta 56 (11), 4369–4377. doi:10.1016/ j.electacta.2011.01.012. http://linkinghub.elsevier.com/retrieve/pii/S0 0134686110 0 0739. Landstorfer, M., Guhlke, C., Dreyer, W., 2016. Theory and structure of the metal-electrolyte interface incorporating adsorption and solvation effects. Electrochim. Acta 201, 187–219. doi:10.1016/j.electacta.2016.03.013. http://linkinghub.elsevier.com/retrieve/pii/S0013468616305278. Landstorfer, M., Jacob, T., 2013. Mathematical modeling of intercalation batteries at the cell level and beyond. Chem. Soc Rev. 42 (8), 3234. doi:10.1039/ c2cs35050e. Larche, F., Cahn, J., 1984. The interactions of composition and stress in crystalline solids. J. Res. Natl. Bur. Stand. 89 (6), 467. doi:10.6028/jres.089.026. Latz, A., Zausch, J., 2011. Thermodynamic consistent transport theory of li-ion batteries. J. Power Sources 196 (6), 3296–3302. doi:10.1016/j.jpowsour.2010. 11.088. http://linkinghub.elsevier.com/retrieve/pii/S0378775310020793. Leichsenring, P., Serdas, S., Wallmersperger, T., Bluhm, J., Schröder, J., 2017. Electro-chemical aspects of IPMCs within the framework of the theory of porous media. Smart Mater. Struct. 26 (4), 045004. doi:10.1088/1361- 665X/aa590e. http://stacks.iop.org/0964- 1726/26/i=4/a=045004?key=crossref. 50831928960c1723ffb03383e645a160. Lin, D., Liu, Y., Cui, Y., 2017. Reviving the lithium metal anode for high-energy batteries. Na. Nanotechnol. 12 (3), 194–206. doi:10.1038/nnano.2017.16. Lin, D., Liu, Y., Pie, A., Cui, Y., 2017. Nanoscale perspective: materials designs and understandings in lithium metal anodes. Nano. Res. doi:10.1007/ s12274- 017- 1596- 1. Lundgren, H., Scheers, J., Behm, M., Lindbergh, G., 2015. Characterization of the mass-transport phenomena in a superconcentrated liTFSI:acetonitrile electrolyte. J. Electrochem. Soc. 162 (7), A1334–A1340. doi:10.1149/2.0961507jes. Ma, Y., 1995. The measurement of a complete set of transport properties for a concentrated solid polymer electrolyte solution. J. Electrochem. Soc. 142 (6), 1859. doi:10.1149/1.2044206. Maitra, M., Sinha, M., Mukhopadhyay, A., Middya, T., De, U., Tarafdar, S., 2007. Ion-conductivity and young’s modulus of the polymer electrolyte PEO ammonium perchlorate. Solid State Ion. 178 (3–4), 167–171. doi:10.1016/j.ssi.2007.01.006. http://linkinghub.elsevier.com/retrieve/pii/S01672738070 0 0306. Marcus, Y., Hefter, G., 2004. Standard partial molar volumes of electrolytes and ions in nonaqueous solvents. Chem. Rev. 104 (7), 3405–3452. doi:10.1021/ cr030047d. Marzantowicz, M., Dygas, J., Krok, F., 2008. Impedance of interface between PEO:liTFSI polymer electrolyte and blocking electrodes. Electrochim. Acta 53 (25), 7417–7425. doi:10.1016/j.electacta.2007.12.047. http://linkinghub.elsevier.com/retrieve/pii/S0013468607014958. McKeen, L.W., 2009. Table of Poisson’s Ratios. In: The Effect of Creep and Other Time Related Factors on Plastics and Elastomers. Elsevier, pp. 373–381. doi:10.1016/B978- 0- 8155- 1585- 2.50013- 3. McMeeking, R.M., Landis, C.M., 2005. Electrostatic forces and stored energy for deformable dielectric materials. J. Appl. Mech. 72 (4), 581. doi:10.1115/1. 1940661. http://AppliedMechanics.asmedigitalcollection.asme.org/article.aspx?articleid=1415426.

M. Ganser, F.E. Hildebrand and M. Kamlah et al. / Journal of the Mechanics and Physics of Solids 125 (2019) 681–713


Monroe, C.W., Delacourt, C., 2013. Continuum transport laws for locally non-neutral concentrated electrolytes. Electrochim. Acta 114, 649–657. doi:10.1016/ j.electacta.2013.10.006. http://linkinghub.elsevier.com/retrieve/pii/S0013468613019208. Monroe, C.W., Newman, J., 2009. Onsager’S shortcut to proper forces and fluxes. Chem. Eng. Sci. 64 (22), 4804–4809. doi:10.1016/j.ces.20 09.05.0 09. http: //linkinghub.elsevier.com/retrieve/pii/S0 0 092509090 03285. Newman, J.S., Thomas-Alyea, K.E., 2004. Electrochemical Systems, 3rd edition. J. Wiley, Hoboken, N.J. Niitani, T., Shimada, M., Kawamura, K., Dokko, K., Rho, Y.H., Kanamura, K., 2005. Synthesis of li[sup +]ion conductive PEO-PSt block copolymer electrolyte with microphase separation structure. Electrochem. Solid-State Lett. 8 (8), A385. doi:10.1149/1.1940491. Ogden, R.W., 2011. Non-linear elastic deformations. Dover Publications. OCLC: 868280590. URL http://public.eblib.com/choice/publicfullrecord.aspx?p= 1897319. Pesko, D.M., Timachova, K., Bhattacharya, R., Smith, M.C., Villaluenga, I., Newman, J., Balsara, N.P., 2017. Negative transference numbers in poly(ethylene oxide)-based electrolytes. J. Electrochem. Soc. 164 (11), E3569–E3575. doi:10.1149/2.0581711jes. Pugal, D., Jung, K., Aabloo, A., Kim, K.J., 2010. Ionic polymer-metal composite mechanoelectrical transduction: review and perspectives. Polym. Int. 59 (3), 279–289. doi:10.1002/pi.2759. Rosato, D., Miehe, C., 2014. Dissipative ferroelectricity at finite strains. variational principles, constitutive assumptions and algorithms. Int. J. Eng. Sci. 74, 162–189. doi:10.1016/j.ijengsci.2013.08.007. http://linkinghub.elsevier.com/retrieve/pii/S0 0207225130 01262. Salvadori, A., Bosco, E., Grazioli, D., 2014. A computational homogenization approach for li-ion battery cells: part 1 - formulation. J. Mech. Phys. Solids 65, 114–137. doi:10.1016/j.jmps.2013.08.010. https://linkinghub.elsevier.com/retrieve/pii/S0 0225096130 01634. Salvadori, A., Grazioli, D., Geers, M., 2015. Governing equations for a two-scale analysis of li-ion battery cells. Int. J. Solids Struct. 59, 90–109. doi:10.1016/j. ijsolstr.2015.01.014. https://linkinghub.elsevier.com/retrieve/pii/S0 0207683150 0 0165. Salvadori, A., Grazioli, D., Geers, M., Danilov, D., Notten, P., 2015. A multiscale-compatible approach in modeling ionic transport in the electrolyte of (lithium ion) batteries. J. Power Sources 293, 892–911. doi:10.1016/j.jpowsour.2015.05.114. https://linkinghub.elsevier.com/retrieve/pii/S0378775315010241. Salvadori, A., McMeeking, R., Grazioli, D., Magri, M., 2018. A coupled model of transport-reaction-mechanics with trapping. part i - small strain analysis. J. Mech. Phys. Solids 114, 1–30. doi:10.1016/j.jmps.2018.02.006. http://linkinghub.elsevier.com/retrieve/pii/S0022509617311067. Singh, M., Odusanya, O., Wilmes, G.M., Eitouni, H.B., Gomez, E.D., Patel, A.J., Chen, V.L., Park, M.J., Fragouli, P., Iatrou, H., Hadjichristidis, N., Cookson, D., Balsara, N.P., 2007. Effect of molecular weight on the mechanical and electrical properties of block copolymer electrolytes. Macromolecules 40, 4578– 4585. doi:10.1021/ma0629541. Stalbaum, T., Pugal, D., Nelson, S.E., Palmre, V., Kim, K.J., 2015. Physics-based modeling of mechano-electric transduction of tube-shaped ionic polymer-metal composite. J. Appl. Phys. 117 (11), 114903. doi:10.1063/1.4914034. Stone, G.M., Mullin, S.A., Teran, A.A., Hallinan, D.T., Minor, A.M., Hexemer, A., Balsara, N.P., 2012. Resolution of the modulus versus adhesion dilemma in solid polymer electrolytes for rechargeable lithium metal batteries. J. Electrochem. Soc. 159 (3), A222–A227. doi:10.1149/2.030203jes. Tang, M., Wang, H., Lee, Y.G., Takeda, Y., Yamamoto, O., Xu, J., Yuan, A., Imanishi, N., 2016. Electrochemical and mechanical properties of polyolefin hard segment with polyethylene oxide conductive phase block copolymers. Solid State Ion. 289, 188–193. doi:10.1016/j.ssi.2016.03.014. http://linkinghub. elsevier.com/retrieve/pii/S0167273816300169. Tikekar, M.D., Archer, L.A., Koch, D.L., 2016. Stabilizing electrodeposition in elastic solid electrolytes containing immobilized anions. Sci. Adv. 2 (7). doi:10. 1126/sciadv.1600320. e1600320–e1600320. Timachova, K., Villaluenga, I., Cirrincione, L., Gobet, M., Bhattacharya, R., Jiang, X., Newman, J., Madsen, L.A., Greenbaum, S.G., Balsara, N.P., 2018. Anisotropic ion diffusion and electrochemically driven transport in nanostructured block copolymer electrolytes. J. Phys. Chem. B 122 (4), 1537–1544. doi:10.1021/ acs.jpcb.7b11371. Truesdell, C., Noll, W., 2004. The Non-Linear Field Theories of Mechanics. Springer Berlin Heidelberg, Berlin, Heidelberg doi:10.1007/978- 3- 662- 10388- 3. Walk, A.C., Huttin, M., Kamlah, M., 2014. Comparison of a phase-field model for intercalation induced stresses in electrode particles of lithium ion batteries for small and finite deformation theory. Eur. J. Mech. A. Solids 48, 74–82. doi:10.1016/j.euromechsol.2014.02.020. http://linkinghub.elsevier.com/retrieve/ pii/S09977538140 0 0539. Xiao, C., Trentaine, P.R., 1997. Apparent molar volumes of aqueous sodium trifluoromethanesulfonate and trifluoromethanesulfonic acid from 283 k to 600 k and pressures up to 20 MPa. J. Solution Chem. 26 (3), 277–294. doi:10.1007/BF02767999. Yang, C., Fu, K., Zhang, Y., Hitz, E., Hu, L., 2017. Protected lithium-metal anodes in batteries: from liquid to solid. Adv. Mater. 29 (36), 1701169. doi:10.1002/ adma.201701169. Young, N.P., Devaux, D., Khurana, R., Coates, G.W., Balsara, N.P., 2014. Investigating polypropylene-poly(ethylene oxide)-polypropylene triblock copolymers as solid polymer electrolytes for lithium batteries. Solid State Ion. 263, 87–94. doi:10.1016/j.ssi.2014.05.012. http://linkinghub.elsevier.com/retrieve/pii/ S0167273814002070. Zhang, T., Kamlah, M., 2018. A nonlocal species concentration theory for diffusion and phase changes in electrode particles of lithium ion batteries. Continuum. Mech. Thermodyn. doi:10.10 07/s0 0161-018-0624-z. Zhu, Z., Asaka, K., Chang, L., Takagi, K., Chen, H., 2013. Multiphysics of ionic polymer-metal composite actuator. J. Appl. Phys. 114 (8), 084902. doi:10.1063/ 1.4818412. Zhu, Z., Chen, H., Chang, L., Li, B., 2011. Dynamic model of ion and water transport in ionic polymer-metal composites. AIP Adv. 1 (4), 040702. doi:10.1063/ 1.3668286.