A thermo-chemo-mechanically coupled model for cathode particles in lithium–ion batteries

As the demand for lithium-ion batteries increases, a better understanding of the complex phenomena involved in their operation becomes crucial. In this work, we propose a coupled thermo-chemo-mechanical model for electrode particles of Li–ion batteries. To this end, we start with a general finite strain continuum framework for the coupled thermo-chemo-mechanical problem and then narrow it down to cathode active particles of Li–ion batteries, particularly to lithium manganese oxide particles. Electrochemical kinetics at the surface of the particle and also heat generation due to current exchange are taken into account. Next, the numerical treatment of the problem using the finite element method is presented. Specific line elements are needed to evaluate the flux of ions at the surface of the particle. Finally, the performance of the proposed model is evaluated using a few representative boundary value problems.


Introduction
Lithium-ion batteries are broadly used in portable consumer electronics, medical equipment, military devices, etc. Recently, the demand for lithium-ion batteries has increased dramatically due to the rise of electric and hybrid electric vehicles. Additionally, Li-ion batteries have the potential to be used for storing renewable energy derived from intermittent energy sources such as wind and sunlight, see Zubi et al. [61] and Diouf and Pode [20]. In comparison to other rechargeable batteries such as lead-acid or nickel-based batteries, Li-ion batteries provide higher specific energy and cycle life, which explains their popularity; however, due to safety reasons, they require a battery management system. With the increase in the demand for Li-ion batteries, researchers continue to enhance their performance and safety while decreasing the production costs.
A Li-ion battery cell consists of four main components: an anode, a cathode, a separator and an electrolyte. The anode is typically made of graphite, the cathode is normally a lithium metal oxide and the electrolyte is usually a lithium salt mixed in an organic solvent. The separator, which is usually made of polyethylene or polypropylene, has almost zero electrical conductivity and is placed between the electrodes to prevent short circuit while allowing the flow of lithium ions through its micropores. Despite the attempts to improve the battery capacity by replacing graphite anodes with alternative materials such as silicon-based anodes, graphite remains the preferred anode material for most industrial applications. Hence, the cathode active materials have been the focus of attention in enhancing battery performance. Some examples of cathode active materials are: lithium cobalt oxide (Li Co O 2 or LCO), lithium manganese oxide (Li Mn 2 O 4 or LMO), lithium nickel manganese cobalt oxide (NMC) and lithium nickel cobalt aluminum oxide (NCA), see Buchmann [8] for a lucid introduction. The reader is also referred to Whittingham [52]. The latter two cathode chemistries mentioned above have risen in popularity in the recent years, especially in electric vehicles. The charging process in Li-ion batteries is driven by an external voltage. During the charging of the battery, electrons move in the A. Nateghi · M.-A. Keip (B) Institute of Applied Mechanics, Chair of Material Theory, University of Stuttgart, Pfaffenwaldring 7, 70569 Stuttgart, Germany E-mail: marc-andre.keip@mechbau.uni-stuttgart.de Predicting the amount of heat generation during the operation of batteries is of key importance in designing thermal management systems. Numerous efforts have been made to predict and simulate heat generation in batteries. A general expression for energy balance in battery cells, which considered electrochemical reactions, phase transitions, mixing effects and Joule heating, was given in Bernardi et al. [4]. In their work, the temperature inside the cell was assumed to be uniform. Gu and Wang [27] presented a coupled thermal and electrochemical model using the volume-averaging approach, which allowed for a spatial variation of temperature within the battery cell. Heat generation in a one-dimensional cell with a porous electrode was simulated by Thomas and Newman [45] and compared with the experimental measurements. A three-dimensional thermal model for batteries with layered structure was given in Chen et al. [9]. Zhang et al. [57] used the expression given by Thomas and Newman [45] to calculate heat generation in LMO cathode particles under potentio-dynamic control. A different formulation of transport phenomena in Li-ion batteries, including heat transfer, was given by Latz and Zausch [33] based on non-equilibrium thermodynamics, see De Groot and Mazur [16]. See also the work of Hu and Shen [29] on variational formulation of coupled thermo-chemo-mechanical processes based on non-equilibrium thermodynamics. Hu et al. [28] presented a continuum model for hollow core-shell electrodes by coupling a one-dimensional electrochemical battery model with a two-dimensional thermal model. Werner and Weinberg [51] proposed, based on thermodynamics of irreversible processes, a coupled thermo-chemo-mechanical model for solid state batteries, which also takes phenomena such as phase transformations, thermal diffusion and thermoelectric effects into account. A coupled thermo-chemomechanical model for graphite, NMC, and LMO electrode particles was proposed by Masmoudi et al. [35] based on thermodynamics of irreversible processes. Their small strain formulation has been supplemented with numerical examples; however, an algorithmic implementation of the model has not been addressed.
Following the previous works by Dal and Miehe [15] and Miehe et al. [38], we propose a coupled thermochemo-mechanical model for electrode particles in Li-ion batteries. Our four-field formulation uses also elements of the previous works by Miehe et al. [36] on variational formulation of diffusion in elastic solids and Miehe et al. [37] on phase field modeling of brittle fracture in thermo-elastic solids. The model is used to simulate the response of LMO cathode particles during battery operation. This finite strain model can be used to simulate other electrode particles, provided that the electrochemical properties of the material and the necessary material parameters are known. In Sect. 2 of this work, a general framework for coupled thermo-chemo-mechanics is formulated. The governing equations of the coupled problem are given and the thermodynamical consistency of the model is discussed. An alternative approach to fulfill the second law of thermodynamics is discussed, which is related in the subsequent section of the paper to Onsager's reciprocal relations. In Sect. 3, the general framework of Sect. 2 is narrowed down to electrode particles of Li-ion batteries. To this end, the corresponding constitutive functions are specified and the electrochemical reaction kinetics and heat generation at the boundary of particles are discussed. In Sect. 4, the numerical implementation of the model with a specific application to electrode particles is presented. Finally, in Sect. 5, the model response is evaluated with the help of a few representative boundary value problems.

Four-field formulation of the coupled problem
This section outlines a general continuum formulation of coupled thermo-chemo-mechanics. We will start by introducing the primary field variables, their gradients and the material and spatial fluxes. We will then discuss the global and local balance laws of the coupled problem, the principle of irreversibility and the constitutive equations that fulfill the principle of irreversibility.

Primary field variables of the coupled problem
The coupled problem under consideration can be treated as a four-field problem. The primary field variables are the deformation of the solid ϕ : the species concentration and the chemical potential c : and the absolute temperature At any time t ∈ R + , the deformation ϕ is a one-to-one mapping of the points X ∈ B 0 of the reference configuration B 0 onto the points x ∈ B t of the current configuration B t . The concentration field c ∈ [0, 1] is a dimensionless measure of the amount of the diffusing species (lithium ions) in unit reference volume of the solid. It may be considered as the molar concentrationc of the diffusing species in unit reference volume of the solid divided by the maximum possible molar concentrationc max . The chemical potential μ is the energy-conjugate variable to concentration. The absolute temperature θ is a positive scalar field defined over B 0 , i.e., θ ≥ 0. See Fig. 1 for a visualization of the primary variables. The material gradients of deformation, (negative) chemical potential and temperature fields are denoted by F, M and F and defined by F := ∇ϕ(X, t), M := −∇μ(X, t), In order to circumvent penetration of matter, we require that the determinant of the deformation gradient be positive, i.e., J := det[F] > 0. The deformation gradient F, its cofactor cof[F] = J F −T and its determinant J serve as the mappings of infinitesimal line, area and volume elements, respectively: deformation field species concentration field chemical potential field temperature field Fig. 1 Primary field variables of the coupled thermo-chemo-mechanical problem. The boundary ∂B 0 of the reference configuration is decomposed into Dirichlet and Neumann boundaries for the deformation field ϕ, the chemical potential field μ and the temperature field θ. The prescribed values on the boundary are marked with an overbar

Stress tensors and flux vectors
Consider an arbitrary part P 0 ⊂ B 0 cut out of the reference configuration B 0 of the body under consideration. Let P t ⊂ B t denote the deformed configuration of P 0 at time t. Thus, we can write P t = ϕ t (P 0 ). Let also ∂P 0 and ∂P t express the boundaries of P 0 and P t , respectively. We will write the balance laws for this arbitrary part, assuming that solid material can not cross the boundary ∂P t of the spatial region P t as it deforms with the body.

Stress tensors
The mechanical force per unit area exerted on P t by the rest of the body B t \P t is called the traction vector and is denoted by t. Based on Cauchy's stress theorem, at any point x ∈ ∂P t there exists a second-order tensor σ , called the Cauchy stress tensor, that relates the outward unit normal n to the traction vector t acting on P t : We define the referential traction vector T by scaling the spatial force tda to the referential area element d A such that the identity tda = T dA holds. The first Piola-Kirchhoff stress tensor P is then defined by forming a relationship, analogous to (6), between the referential outward unit normal N and traction vector T : The area mapping (5) 2 relates the infinitesimal area element d A = NdA to da = nda; thus, the relationship between the first Piola-Kirchhoff stress tensor P and the Cauchy stress tensor σ can be obtained as

Species flux vectors
The species flux h represents the amount of the species that flows outwards, per unit time, through a unit area element of the boundary ∂P t of the part P t . It is related to the outward normal n by a flux vector h in the form The material species flux, denoted by H and defined by the identity hda = H dA, represents the flow of the species, per unit time, through a unit area element of the boundary ∂P 0 . In analogy to (9), we assume that the material species flux H is related to the outward normal N to ∂P 0 by a material species flux vector H Using the area mapping (5) 2 , the relationship between the material species flux vector H and the spatial species flux vector h is obtained:

Heat flux vectors
The heat flux q is the outward flow of thermal energy, per unit time, through a unit area element of the boundary ∂P t . It is related to the outward normal n to the surface ∂P t by the heat flux vector q q(x, t; n) = q(x, t) · n.
The material heat flux Q is defined by the identity qda = QdA and represents the amount of thermal energy which flows outward, per unit time, through a unit area element of the boundary ∂P 0 . Similar to (12), the material heat flux Q is related to the outward normal N to ∂P 0 by a material heat flux vector Q With the area mapping (5) 2 , the relationship between the material heat flux vector Q and the spatial heat flux vector q can be written as 2.3 General equations of thermo-chemo-mechanics

Global balance equations
The general equations governing the coupled problem are written for the arbitrary part P 0 ⊂ B 0 cut out of the reference configuration B 0 . No mass production is assumed to take place. During the deformation the total mass of the solid constituent in the part P 0 should remain unchanged and equal to its initial value. Considering the part P 0 as a control volume, balance of species content demands that the temporal changes in the amount of the diffusing species be equal to the in-flux (negative out-flux) through the surface d dt Neglecting the inertial effects, balances of linear and angular momentum for the part P 0 can be written as respectively. In the above equations, γ (X, t) is the body force exerted per unit reference volume. Finally, the balance of energy reads d dt where V (X, t) :=φ(X, t) is the material velocity field, e(X, t) is the internal energy of the unit reference volume and r (X, t) is the heat source per unit reference volume.

Local balance equations
Using the definitions introduced in Sect.

Principle of irreversibility
The global form of the second axiom of thermodynamics for an arbitrary part P 0 ⊂ B 0 of the body is d dt where η(X, t) is the entropy per unit reference volume. The local form of the entropy imbalance readṡ Using the balance of energy (18) 4 and introducing the Helmholtz free energy ψ by the Legendre transformation e = ψ + ηθ, the above inequality takes the forṁ Multiplying the above inequality by the absolute temperature and defining G := − ∇θ θ , the local form of the second axiom is obtained in the form of the dissipation inequality This poses a constraint on the material modeling framework. It needs to be fulfilled locally such that the second axiom of thermodynamics is not violated. A stronger constraint is achieved by demanding the positiveness of some sets of the individual terms in the dissipation inequality. To this end, the dissipation inequality (22) is split up into an internal part and a part due to diffusion and heat conduction and each of these parts is demanded to be non-negative: An even stronger constraint may be achieved by splitting the dissipation inequality (22) into an internal part and two separate parts due to diffusion and heat conduction. Each of these parts is demanded to be non-negative separately:

Objective free energy and thermodynamical forces
In order to close the problem, we need to specify the constitutive functions such that the principle of irreversibility is fulfilled. Assuming a local theory of grade one, the Helmholtz free energy ψ is assumed to depend on the primary variables ϕ, c and θ and their first gradients, The stored energy functionψ must remain invariant under a rigid body motion of the form ϕ + = Rϕ + c superimposed on the current configuration. Here, R(t) ∈ SO(3) denotes an arbitrarily chosen rotation tensor of the special orthogonal group and c(t) is a translation vector. Based on this requirement, the stored energy functionψ can not depend on ϕ and depends on the deformation gradient F through the right Cauchy-Green tensor C := F T F, which remains invariant under this rigid body motion. This conclusion together with a standard Coleman-Noll procedure applied on the internal dissipation inequality (23) 1 or (24) 1 gives the reduced set of arguments for the stored energy and the constitutive equations for the thermodynamical forces Note that by making the stored energy function dependent on the right Cauchy-Green tensor C, the balance of angular momentum is automatically fulfilled.

Temperature evolution equation
With e = ψ + ηθ and the constitutive equations given in (27), balance of energy (18) 4 takes the form of a temperature evolution equation where the following definitions are used for the sake of compactness:

Thermo-chemically coupled dissipation potential
A thermodynamically consistent setting can be achieved by demanding that the internal part of the dissipation inequality (23) 1 and the part due to diffusion and heat conduction (23) 2 be non-negative. Let us repeat these inequalities: Following a standard Coleman-Noll procedure and using the constitutive equation (27), the internal dissipation inequality (30) 1 can a-priori be fulfilled. What remains is to restrict the material modeling framework such that the dissipation due to diffusion and heat conduction (30) We assume that the flux vector J is specified by an objective dissipation potential which depends on the driving force I at a given state {F, c, θ} of deformation, concentration and temperature. For a convex, positive and normalized φ, the dissipation inequality (31) is a-priori fulfilled if the flux vector J is specified by the following constitutive equation: Convexity of the dissipation potential in I implies that the Hessian H := ∂ 2 II φ is positive semi-definite. Moreover, H is a symmetric tensor due to symmetry of second derivatives. Let us divide H into four blocks as below: where we a and c are themselves symmetric. Positive semi-definiteness of H implies a constraint on the blocks a, b and c. For a positive semi-definite block a 0, the tensor given in (34) is positive semi-definite if and only if the Schur complement H/a of a is also positive semi-definite, that is With (27) and (33) at hand, we can summarize the constitutive equations in the following box:

Decoupled chemical and thermal dissipation potentials
A simpler formulation with decoupled thermal and chemical dissipation potentials is obtained by demanding that the dissipation inequalities (24) be satisfied. Let us repeat these inequalities: Following a standard Coleman-Noll procedure and making use of the constitutive Eq. (27), the internal dissipation inequality (37) 1 can a-priori be fulfilled. We need to restrict the material modeling framework such that the dissipation inequalities due to diffusion (37) 2 and heat conduction (37) 3 are also satisfied. To this end, the material species flux vector H is assumed to be specified by the objective chemical dissipation potential which depends on the (negative) gradient M of chemical potential at a given state {F, c, θ} of deformation, concentration and temperature. Assuming a convex, positive and normalized chemical dissipation potential φ c in terms of M , the dissipation inequality (37) 2 is a-priori satisfied if the material species flux vector H is specified by the constitutive equation Similarly, we assume the material heat flux vector Q to be specified by an objective thermal dissipation potential which depends on the vector G at a given state {F, c, θ} of deformation, concentration and temperature. For a convex, positive and normalized φ t , the dissipation inequality (37) 3 is a-priori fulfilled if the material heat flux vector Q is specified by the following constitutive equation: With (27), (39) and (41) at hand, we can summarize the constitutive equations in the following box:

Boundary conditions of the coupled problem
The boundary ∂B 0 of the body in the reference configuration is split into Dirichlet and Neumann boundaries for deformation ϕ, chemical potential μ and temperature θ fields such that Given are the Dirichlet and Neumann boundary conditions for the deformation for the chemical potential and for the temperature withφ,μ andθ being the prescribed deformation, chemical potential and temperature and withT ,H and Q being the prescribed traction, species outflux and heat outflux. The governing equations of the coupled problem, which are used to derive the weak form, are summarized in the following box: Note that in the numerical treatment of our four-field setting, the constitutive equation for chemical potential (48) 2 has been treated as a global equation. The reason is that a three-field formulation in terms only of the primary field variables {ϕ, μ, θ} will lead to complications in evaluating the species flux at the boundary of the particles using the Butler-Volmer equation, see Sect. 3.3; thus, we consider the concentration c as a primary (global) variable and regard the constitutive equation for chemical potential as an additional global equation.

Model problem: electrode particles of Li-ion batteries
In what follows, we reduce the continuum formulation of Sect. 2 to the specific case of electrode particles in lithium-ion batteries. To this end, we present the specific choice of constitutive functions for the stored energy and the dissipation potentials. Then, we comment on the kinetics of electrochemical reaction at particle interface and finally, we discuss heat generation due to current exchange at the boundary of particles.

Constitutive functions for the stored energy
A multiplicative decomposition of the deformation gradient into an elastic part F e , a chemical part F c and a thermal part F t is considered: where the chemical and thermal parts are volumetric deformations in terms of the chemical and thermal expansions Here, c 0 and θ 0 are the initial concentration and temperature, respectively. is a swelling parameter and α is the coefficient of linear thermal expansion. The stored energy function is assumed to additively decompose into elastic, chemical and thermal contributions The elastic part ψ e of the stored energy function depends on the elastic part F e of the deformation gradient and is of compressible Neo-Hookean type Here, γ is the shear modulus and the parameter ζ > 0 depends on the Poisson ratio ν by the relationship ζ = 2ν 1−2ν . The purely chemical contribution to the stored energy function is given by where A is a material constant. Finally, the purely thermal contribution to the stored energy function has the form where the heat capacity C is assumed to be constant.

Constitutive functions for the dissipation potentials
A simple choice for the dissipation potential function discussed in Sect. 2.3.6 would be a quadratic function in terms of M and G, where a = a T and c = c T , while the constraint (35) holds and the block a is positive semi-definite a 0. Note that a, b and c generally depend on deformation, concentration and temperature. The choice of a quadratic dissipation potential function leads to a linear dependence of J on I in the form J = H · I, where H is a symmetric positive semi-definite tensor independent of I. Using the representation (34) of H, we can write The tensor b couples species flow to temperature gradient and heat flow to the gradient of chemical potential. These couplings are called the Soret effect and the Dufour effect, respectively, see De Groot and Mazur [16] and Rahman and Saghir [41]. Note that the symmetry of H is in line with Onsager's reciprocal relations, see De Groot and Mazur [16]. We choose the following specific form for the dissipation potential function (55) where M, L and K are material parameters and C is the right Cauchy-Green tensor. Positive semi-definiteness of H imposes a constraint on the material parameters M, L and K . Since c ∈ [0, 1] and the inverse of the right Cauchy-Green tensor C −1 is symmetric and positive definite, based on (35), the following constraints need to be fulfilled such that the second law of thermodynamics is not violated: In Sect. 5.1, we will comment on the response of the coupled thermo-chemo-mechanical model using the above formulation. Setting L = 0 in (57) and thus, b = 0 in (55) and (56), brings us the decoupled setting discussed in Sect. 2.3.7, where the following chemical dissipation potential function characterizes the species flow at a given state of deformation represented by the right Cauchy-Green tensor C and concentration c. Here, M > 0 is identified as the species mobility, see Miehe et al. [36]. The thermal dissipation potential function governs the heat flow and has a similar structure to the chemical dissipation potential in (59). It is a convex at a given state of deformation and concentration, see Miehe et al. [37]. Here, the parameter K > 0 represents the thermal conductivity.

Electrochemical reactions at the boundary of particles
During discharging of the battery, lithium ions migrate from anode through the electrolyte to cathode and the electrons flow in the external circuit in the same direction, whereas during charging, lithium ions migrate from cathode through the electrolyte to anode. An external voltage is needed for the charging process. The Following the approach of Fuller et al. [24] and Newman and Thomas-Alyea [39], the rate of the electrochemical reaction is related to the surface overpotential η s , which serves as the driving force of the reaction. Surface overpotential is the difference between the potential V of the working electrode and the open-circuit potential U ocp ; thus, we have η s = V − U ocp . The open-circuit potential U ocp is the potential of an electrode of the same material which is in equilibrium with the electrolyte. It is a function of the state of charge, which is expressed by the stoichiometric ratio y = 1 − x of lithium in the electrode Li y Mn 2 O 4 . This can be considered as equivalent to the unitless concentration c of lithium ions in the electrode. Following the work of Doyle et al. [22], the open-circuit potential of LiMn 2 O 4 is given in terms of the state of charge by the curve of Fig. 2.
Additionally, similar to Zhang et al. [55] we can assume that, due to the small size of the particle, the potential within the particle is homogeneous and equal to the applied potential U applied under potentiodynamic stimulus. We can then summarize The current density i at the surface of the particle is related to the surface overpotential η s by the Butler-Volmer equation Here, i 0 is the exchange current density, β is the symmetry factor, R is the gas constant and F is Faraday's constant. The exchange current density i 0 is given by wherec amb is the molar concentration of lithium ions in the electrolyte and k is the reaction rate constant. Finally, the flux of lithium ions through the particle surface is given bȳ We will assume that concentration of lithium ions in the electrolytec amb remains almost constant during the operation of the battery. Hence, we consider the multiplier k(c amb ) (1−β) as a model parameter and denote it by H in the forthcoming sections.  [46] has been performed using the sum of sines with 4 terms in MATLAB

Heat generation due to current exchange
Neglecting side reactions, we assume, similar to Masmoudi et al. [35], that heat generation within the particle can be evaluated based on the current exchange at the interface of the particle. Two contributions to the generated heat are considered, namely the irreversible resistive heating and the reversible entropic heating. The resistive heating is denoted by R r and is given in the form of the following integral over the surface of the particle where i is the current density determined by the Butler-Volmer equation (63) and η s is the surface overpotential (62), see Sect. 3.3. The entropic heating R e is given by  Fig. 3, to the experimental measurements has been performed using the sum of sines with 4 terms in MATLAB. This specific choice has been made to avoid the complexities of implementing a piecewise function in our numerical treatment. The fit is given by The overall heat generation due to current exchange would be the sum of (66) and (67) which will be applied within the particle as a local heat source. To this end, we will use the following equation to evaluate the local heat source r per unit reference volume We will assume a homogeneous distribution of the local heat source r in B 0 . Thus, one can write where |B 0 | denotes the total volume of B 0 .
Note that in (74), the heat source r is frozen at its previous value r n at time t n . This approximation simplifies the numerical treatment of the problem dramatically. We only need to evaluate the current heat source r and save it as a history variable to be used in the next time step.

Space discretization
We construct the weak form of the problem using a standard Galerkin method. For this purpose, the timediscrete forms of the governing equations are multiplied by arbitrarily chosen test functions out of the admissible function spaces listed below: The weak forms of the governing equations are then obtained as for all admissible test functions δϕ, δμ, δc and δθ. Using Gauss's theorem, the Cauchy-like relationships (7), (10) and (13), and the Neumann boundary conditions (45) 2 , (46) 2 and (47) 2 , the above equations can be rewritten as Note that in Eq. (77) 3 , the species flux at the boundary of the particleH is evaluated using the Butler-Volmer equation (65). For this reason, in addition to discretization of the referential volume B 0 , a discretization of the Neumann boundary ∂B h 0 using surface elements is required. Let us consider a finite element discretization of the bulk of B 0 and the Neumann boundaries ∂B t 0 , ∂B h 0 and ∂B q 0 as follows: Here, N is the number of the bulk finite elements, whereas N t , N h and N q are the numbers of the surface finite elements used for the discretization of traction, species flux and heat flux boundaries. The following set of approximations have been used for the primary variables and the corresponding gradients within the bulk Note that these shape functions should be compatible with the shape functions used for the bulk elements.
Introducing the above approximations into the weak form (77) leads to the residual expressions The above coupled system has to be solved for the nodal displacements d ϕ , concentrations d c , chemical potentials d μ and temperatures d θ . Let us introduce the following shorthand notation for the global vectors of residuals and generalized nodal displacements The problem under consideration is then to find d such that R(d) = 0. A standard Newton-Raphson scheme can be used to solve this system of nonlinear equations. The vector of generalized nodal displacements is updated by until a convergence criterion is met, e.g., |R| < tol. Here, K is the tangent of the Newton-Raphson scheme.

Evaluation of the model response
In what follows, we will show the results of a few numerical simulations obtained using the model proposed in the previous sections. In Sect. 5.1, we will comment on the results obtained by the coupled formulation of the dissipation potential, which was discussed in Sect. 2.3.6. This academic study does not directly apply to modeling of electrode particles in Li-ion batteries. The latter topic will be treated later in Sects. 5.2 and 5.3, where we will discuss the response of two cathode particles with different geometries under a potentio-dynamic loading.
K Thermal conductivity 10.0 mN/s

Model response of the fully coupled formulation
To evaluate the model response of the fully coupled formulation discussed in Sect. 2.3.6, we consider two thermo-chemo-mechanical loading scenarios on a square block with the side length of 5µm. The block is mechanically restricted such that it can only expand in horizontal direction, see Fig. 4a. Plane strain conditions are assumed. The material parameters used for this study are listed in Table 1 and do not represent any real material. For the discretization of the sample, elements with bi-linear shape functions for all of the unknown fields are used. Since no species flux boundary conditions are considered, there is no need to use surface elements at the boundaries. In the first numerical experiment, a linearly increasing chemical potential μ =μ is applied on the right side of the block. See Fig. 4b for the loading function. No species flux is assumed to take place on all of the other sides. Furthermore, no heat flux through the boundaries is permitted. The problem is indeed a one-dimensional problem and the primary fields vary only in X direction. For three different values of the coupling parameter L = {0.0, 0.5, −0.5} µm 2 /s, introduced in Sect. 3.2, the values of species concentration c, chemical potential μ and temperature θ are plotted on the mid-line of the block in Fig. 5 at the final stage of loading t f = 100 s. Let us consider the case when L = 0.0, shown in the graphs by solid black lines. Since chemical potential is increased on the right side of the block, its gradient points in the X direction. Thus, M points in the negative X direction, which by (56) creates a flux vector H also in the negative X direction. This leads to an increase of the concentration field within the block, which can be observed in Fig. 5a. Furthermore, by the temperature evolution Eq. (48) 4 , we can expect an increase of the temperature within the block, observed in Fig. 5c. For a positive L, an additional term in the negative X direction will be added to the heat flux vector Q in (56). This term will lead to a temperature decrease on the right side of the block shown in the temperature graph. The opposite effect takes place for a negative L.
In the second numerical experiment, we apply a linear temperature rise θ =θ on the right side of the square, see Fig. 4. No species flux is permitted at the boundaries. Heat transfer is also prevented on the top, bottom and  . 7 a Boundary conditions applied on the circular particle. Due to symmetry, only one quarter of the particle is discretized. b The electric potential profile used by Zhang et al. [57] is applied on the particle. This induces a species flux through the boundary and simulates a charging-discharging cycle left sides of the block. For the above-mentioned values of the coupling parameter L, the species concentration c, chemical potential μ and temperature θ plots on the mid-line of the block are given in Fig. 6 at the final stage of loading t f = 100 s. Even though the temperature graphs of Fig. 6c seem identical for L = 0.5 µm 2 /s and L = −0.5 µm 2 /s, the corresponding concentration graphs show a totally different pattern. For a positive L, an additional term in negative X direction will be added to the species flux vector H in (56). This term will lead to a concentration decrease on the right side of the block and consequently, to a concentration rise on the left side observed in Fig. 6a. The opposite effect is observed for a negative L.

A circular disk under potentio-dynamic loading
We evaluate the model response by considering a circular LiMn 2 O 4 cathode particle of radius r = 5µm. The exterior of the particle is assumed to be free of traction. Additionally, it is assumed that the temperature at the boundary of the particle remains fixed since it is in constant contact with the electrolyte. We assume that the temperature of the electrolyte does not change. Due to symmetry, only a quarter of the particle is discretized by bi-linear bulk elements and linear surface elements. See Fig. 7a for the geometry of the particle and the applied boundary conditions. Plane strain conditions are considered in the simulation. A species flux through the boundary is induced by an electric potential applied on the particle. The applied electric potential follows the sweeping profile given in Zhang et al. [57], which simulates a charging-discharging cycle. See Fig. 7b for the applied electric potential. The initial values of the concentration and temperature are chosen to be c 0 = 0.996 and θ 0 = 298.15 K = 25 • C (room temperature), respectively. The material parameters given in Table 2 are chosen as realistically as possible under the assumptions of our modeling framework. The coupling parameter L is set to zero. The governing equations of the problem (48) are solved using the finite element method to evaluate the four unknown fields shown in Fig. 1. Flux of species H induced by the applied potential is shown in Fig. 8a. In order to compare this result with the flux profile given in Zhang et al. [57], the computed values of H are multiplied by the maximum molar concentrationc max = 2.37 × 10 4 mol/m 3 . Figure 8b shows a comparison of the fluxes obtained by the simulation with the flux profile of Zhang et al. [57]. A relatively good agreement can be observed between the two results. Li-ion concentration versus time is shown in Fig. 9a at the center and on the surface of the particle. Due to the small size of the particle, variations of concentration within the particle are relatively small. This can be observed better in Fig. 9b, which shows the variations of Li-ion concentration within the particle in the radial direction. At some time instances, there are very small concentration gradients within the particle, while at some other instances, concentration of Li-ions can not be considered as homogeneous. The generated resistive and entropic heat at the particle interface are plotted in Fig. 10 in the first half of the loading cycle. As it can be observed in the figure, the entropic heat is more pronounced compared to resistive heat. The temperature evolution equation (48) 4 is solved by treating the overall generated heat R e + R r as a homogeneously distributed source of heat inside the particle, determined by Eq. (72), and using the time  The flux shown on the left is multiplied byc max to obtain the molar species fluxH . It is then compared to the results obtained of Zhang et al. [57] discretization approach discussed in Sect. 4.1. It turns out that the heat generated due to current exchange at the boundary of the particle as well as the heat generated due to diffusion and deformation do not have a major effect on the temperature inside the particle. Temperature within the particle remains almost constant and equal to the temperature of the surrounding environment. Hence, the dominant factor that determines the temperature inside the particle is the surrounding temperature. This observation can be attributed to the small size of the particle and hence, to the fast dissipation of the generated heat to the surroundings of the particle. If instead of the Dirichlet temperature boundary condition at the boundary of the particle, shown in Fig. 7, a Neumann boundary conditionQ = 0 for heat flux was considered, the temperature within the particle would increase dramatically. See Fig. 11 for the results of the simulation with zero flux boundary condition. This temperature rise is due to deformation, diffusion and heat generation caused by electric current exchange. These results underline once more the significance of heat dissipation to the surrounding environment without which, temperature rise is extreme.

A two-dimensional particle geometry under potentio-dynamic loading
In a second study, we consider a realistic two-dimensional geometry of a LiMn 2 O 4 cathode particle shown in Fig. 12. This geometry is a two-dimensional reduction of the scanning electron microscope image provided in Bohn [5]. Similar boundary conditions to the previous example are assumed on the exterior of the particle,  66) and (67), are plotted versus time in the first half of the loading cycle. Integration over the entire surface of the particle is carried out with the help of the line finite elements. The overall heating R e + R r is applied as a homogeneously distributed heat source inside the particle determined by Eq. (72) i.e., the boundary of the particle is traction-free, the surrounding temperature is constant and a species flux on the boundary is induced by applying the electric potential profile shown in Fig. 7b. The initial conditions for concentration and temperature inside the particle are c 0 = 0.996 and θ 0 = 298.15 K (room temperature), respectively. The coupled problem is solved using the finite element method. The geometry has been discretized using 598 bulk elements and 100 surface elements. The material parameters of Table 2 are used for this problem. The coupling parameter L is assumed to be equal to zero. Figure 13 shows the contour plots of the species concentration at different time instances in the first half of the charging-discharging cycle and also at the final time t = 4000 s. As it can be observed in the figure, concentration of Li-ions decreases during the first half of the loading, leading to a decrease in the volume of the particle. The initial volume of particle is recovered at the end of the loading cycle. The temperature rise within the particle is very small and can be ignored as far as the ambient temperature does not change. Once more, the ambient temperature turns out to govern the temperature evolution within the particle. This can be attributed to the small size of particle and fast dissipation of the produced heat to surrounding environment. In this work, we presented a four-field formulation of coupled thermo-chemo-mechanical problems with application to the simulation of cathode particles in lithium-ion batteries. First, the governing equations of the coupled problem were discussed in a general framework. An alternative way of satisfying the second law of thermodynamics using Onsager's symmetric positive-definite matrix of coefficients was explained. This formulation enables us to take the Soret and Dufour effects into account. Subsequently, the coupled formulation was applied to the simulation of electrode particles in lithium-ion batteries. The generated heat due to current exchange at the boundary of the particle was considered as a volume heat source within the particle. Thereafter, time and space discretization of the problem were discussed. The response of the model was evaluated under a potentio-dynamic loading applied on different particle geometries. Temperature evolution within the particle was found to be very minor. The ambient temperature was found to be the most important factor governing the temperature distribution within the particle.