A conservative and efficient model for grain boundaries of solid electrolytes in a continuum model for solid-state batteries

A formulation is presented to efficiently model ionic conduction inside, i.e. across and along, grain boundaries. Efficiency and accuracy is achieved by reducing it to a two-dimensional manifold while guaranteeing the conservation of mass and charge at the intersection of multiple grain boundaries. The formulation treats the electric field and the electric current as independent solution variables. We elaborate on the numerical challenges this formulation implies and compare the computed solution with results from an analytical solution by quantifying the convergence towards the exact solution. Towards the end of this work, the model is firstly applied to setups with extreme values of crucial parameters of grain boundaries to study the influence of the ionic conduction in the grain boundary on the overall battery cell voltage and, secondly, to a realistic microstructure to show the capabilities of the formulation.


Introduction
The technology of solid-state batteries (SSBs) has made remarkable progress in recent years towards their usage in real-world applications.However, some key challenges still remain and are part of current research [1].One of the challenges is to tune the grain boundaries inside polycrystalline solid electrolytes to obtain desired properties.They occur in polycrystalline solid electrolytes between the grains of oxides (e.g., garnet LLZO) or sulfides (e.g., LPS) [2].Grains can be defined as the geometric domains where atoms are periodically arranged [2], and therefore, grain boundaries are the locations where this periodicity ends.The role of grain boundaries in terms of their influence on the overall cell performance is part of the discussion in the scientific community.Among others, the following physical effects are observed at grain boundaries: The deposition and growth of lithium filaments [3], a different ionic conductivity inside of the grain boundary compared to the grains [4,2], the grain boundary as an electron conductor in the solid electrolyte [5] and in the SEI [6], the grain boundary as an ionic resistor [4], or cracking along the grain boundary [2].In the past years, results have been reported to modify and thereby improve the properties of grain boundaries [7,8,9,10].A model that resolves the ionic conduction along the grain boundaries is needed to best profit from the possibility of tuning the properties of grain boundaries.Such a model is also demanded in the literature [2].Such a model allows quantifying the competing transport mechanisms through the grains and through the grain boundaries [4].Especially for sulfides with high grain boundary conductivity, which could reach the magnitude of the bulk conductivity [11], such a model becomes inevitable.Different types of models have been reported to capture the effect of grain boundaries, ranging from atomistic models [4,12,13] to continuum models which discretize the grain boundaries by a phase field [14].However, a continuum model that geometrically resolves the grain boundaries sharply, i.e., not by a smeared phase field, and is applicable to realistic microstructures, is still missing.One challenge is that the thickness of the grain boundary, which is reported to be less than 10 nm [15,16], is significantly smaller than the typical length scale of a realistic microstructure of a battery cell, and therefore, a standard numerical discretization would exceed the limits of available computational resources (c.f.[17]).Another challenge comes with the need to incorporate the intersection of multiple grain boundaries at one point while maintaining conservation properties.In this paper, we propose a continuum model that resolves the ionic conduction property of grain boundaries across and along the grain boundary together with a mass and charge-conserving formulation at the intersection of grain boundaries (Section 2).Afterwards, the model is analyzed in terms of convergence, and conservation properties before it is applied to realistic microstructures (Section 3).The model is parametrized using experimental studies that differentiate the contribution of measured ionic conductivities into the bulk and the grain boundary [18,19].Many sources indicate that the ionic conductivity of the grain boundary is usually smaller than the conductivity of the solid electrolyte [4,2,19].However, it was shown that the electrochemical properties of grain boundaries can be tuned using surface modification by coating the grains with another solid electrolyte [20] or modifications in the lattice structure [21].Thus, we varied the ionic conductivity of the grain boundary in a wide range to quantify its influence on the effective conductivity of the solid electrolyte.

An efficient and conservative model for grain boundaries within SSBs
In this paper, we present a novel model to include the effect of ionic conduction across and along the grain boundaries into a microstructure-resolved continuum model for solid-state batteries.
The effect of grain boundaries on the entire cell performance in terms of ion conduction can only be studied if the grain boundaries are geometrically resolved within the solid electrolyte to also account for potentially preferred conduction paths along the grain boundaries.Therefore, a set of equations describing the conservation of mass and charge in the electrodes, the electrolyte, and the current collectors, as well as inside the grain boundaries, is required.

Geometric definitions
In Figure 1, an SEM image of grains (bright domains) and grain boundaries (dark lines) of LLZO is shown.One geometric characteristic of grain boundaries is that more than two can intersect at one line (in the two-dimensional picture at one point).An abstract geometric schematic containing all relevant domains, interfaces, and boundaries of a full battery cell (c.f. Figure 2) includes: The current collector on the anode side Ω cc,a , the anode Ω a , the solid electrolyte Ω el , the cathode Ω c , and the current collector on the cathode side Ω cc,c .At the intersection of two domains (Ω i , Ω j ) interfaces are defined as Γ i−j .The outer boundary is split into real physical boundaries at the current collectors Γ cc-o and modeling boundaries Γ cut that are required to limit the size of the geometry in the lateral direction.Additionally, grain boundaries are defined within the solid electrolyte Ω gb with the interface between them and the solid electrolyte Γ gb-el .They intersect at one location marked with '*'.Note that this location is a point in two-dimensional setups and lines and points in three-dimensional setups.The electrodes Ω ed = Ω a ∪ Ω c and the current collectors Ω cc = Ω cc,a ∪ Ω cc,c are summarized to simplify the notation.

Governing equations in a continuum formulation
We introduced the set of conservation equations of charge and mass for SSBs in a continuum form in our previous work [22].The governing equations are with the overvoltage η = Φ ed − Φ el − Φ(c).The used symbols are listed in Table 1.Note that we consider the ionic conductivity κ el as constant and isotropic inside the grains of the solid electrolyte within this work.However, including an anisotropic and inhomogeneous ionic conductivity would be conceptually easy.This set of governing equations is extended by equations to model the transport of ions in the grain boundaries and the exchange of ions between the grains and the grain boundaries.We model the grain boundaries as single-ion conductors, such that the charge density ρ gb and mass density, expressed in terms of the concentration c gb inside the grain boundaries, are linked by the charge number z, the transference number t + = 1, and the Faraday constant F : ρ gb = zF c gb .Therefore, it is sufficient to only model the conservation of charge inside of the grain boundaries.For the derivation of the conservation of charge inside of the grain boundaries, we begin with the general equation for the conservation of charge ρ gb with the electric current i gb .The constitutive equation for the electric current is given by to model ion conduction inside of the grain boundary [2].No accumulation of charge occurs due to the assumption of electro-neutrality and a transference number of one, i.e.
∂ρ gb ∂t = 0.In our previous work [17], we showed that thin layers in SSBs can be modeled as two-dimensional manifolds if the transport phenomena that are orthogonal to the layer can be approximated by a priori knowledge.Often, this is given if both the local curvature of the layer and the thickness of the layer are significantly smaller than a typical length scale of the surrounding geometry.A similar concept is used in this work to model the conservation of charge inside of the grain boundaries.Here, we assume that the normal component of the electric field E n = E • n gb = ∇Φ • n gb is constant throughout the thickness of the grain boundary while it may vary along the grain boundary (see Figure 3 for the definition of n gb ).This leads to the following set of equations with the Nabla operator ∇ Γ being evaluated on curved surfaces with the Neumann boundary of a grain boundary ∂Γ gb h .The surface Γ gb is defined in the center of Ω gb .
From geometric considerations, it becomes obvious that no external fluxes across the outer edges of the grain boundary occur, i.e., ī = 0 on Γ gb ∩ Γ cut and Γ gb ∩ Γ ed .A model for the electric current between the grain boundaries and the solid electrolyte is formulated by a simple linear kinetics law with the contact resistance r c between the solid electrolyte and the grain boundary domain.However, any other kinetics law is applicable as well.A geometric characteristic of grain boundaries is that more than two can intersect in one line (see the location of the * symbol in Figure 2).The conservation of charge needs to be guaranteed at these locations as well.When modeling the grain boundaries as a two-dimensional continuum, the formulation of conservation of charge at these locations is equivalent to Kirchhoff's circuit laws in electrical circuits: (1) the electric potential at the tip of all branches of the grain boundaries is equal, and (2) all currents into and out of the branches have to sum up to zero.This results in the following two constraints wherever n grain boundaries intersect with the normal vector n i on the tip of the branches of the grain boundaries Γ gb i as defined in Figure 3. exemplary slope of the electric potential (blue line) and of the electric current (red line) is added to Figure 3 to visualize the constraints at the intersection point.While the electric potential is equal (denoted by the circular dotted line), the currents need to sum up to zero (denoted by the same length of the red arrows).Note that only the contribution of the flux out of the grain boundary is constrained in a three-dimensional setup.

Numerical realization of the model
First, we summarize the numerical discretization of the novel continuum model.Afterwards, the solution strategy applied to the discretized system of equations is discussed.

Discretization of the model and solution of the nonlinear system
For the discretization of the novel model in space, we apply the finite element method and, in time, the one-step theta method.The resulting nonlinear, algebraic system is solved by the iterative Newton-Raphson scheme.We propose to include the constraints formulated in Equations (6a) and (6b) by treating both the electric potential Φ and the electric current i as independent solution variables in the domain of the grain boundary Γ gb .In the remaining domains, the electric potential Φ and the concentration c are considered as solution variables.This results in more unknowns inside the domain of the grain boundary (i.e., the electric potential and the vector-valued electric current per node).However, we believe that this is affordable as the number of nodes in the grain boundary is significantly smaller than the total number of nodes if the grain boundary is resolved as a two-dimensional manifold.The weighted residual in the domain of the grain boundary consists of the sum of the standard residual R gb,std , and the residua originating from the constraints of the electric potential R gb,constr Φ and of the electric current R gb,constr i : R gb = R gb,std + R gb,constr Φ + R gb,constr i = 0.The standard residual is given by multiplication of the balance equation (Equation (4a)), the constitutive equation (Equation (4b)), and the homogeneous Neumann boundary (Equation (4c)) with an arbitrary test function w and subsequent integration over the entire domain or boundary, respectively Applying Gauss divergence theorem and the product rule of the divergence to the first term leads to where we make use of w Φ = 0 on ∂Γ gb \∂Γ gb h .The constraints (Equations (6a) and (6b)) are enforced by introducing Lagrange multipliers The constraint matrix C i is chosen, such that the flux constraint (Equation (6b)) is satisfied.Therefore, the constraint is split into a slave and a master side and reorganized with n s 1 being the largest component of the normal vector n s (n 1 s = max(n x s , n y s , n z s )), and n s 2 , n s 3 the other components In the remaining domains, only the electric potential and the concentration are considered unknown, leading to the following discretized form of the weighted residual The residua R gb and R bulk are discretized in space using the finite element method, and the resulting nonlinear, algebraic system is solved using the Newton-Raphson scheme (see Appendix A for the derivation of the linear system of equations in the domain of the grain boundaries and Appendix B for the remaining domains).Thereafter, two linearized systems of equations (K gb ∆Ψ gb = −R gb , K bulk ∆Ψ bulk = −R bulk ) arise The systems of equations for the grain boundaries and the bulk domains are coupled by the coupling flux defined in Equation ( 5).Therefore, the systems of equations contain the additional contributions

Solution of the linear system
Due to robustness and efficiency, this linear system of equations (Equation ( 13)) is solved monolithically [23].
Iterative solvers are required as soon as the size of the linear system of equations exceeds a certain threshold, which is often the case for realistic microstructures.This monolithic system could include a zero-block on the main-diagonal sub-block K gb ΦΦ if the source term s ρ gb is not a function of the electric potential, i.e.
∂sρ gb ∂ Φgb = 0.Many iterative solvers fail for such saddle-point problems, i.e. systems that contain a zero sub-block on the main diagonal (c.f.[24]).However, the determinant of the entire matrix K is non-zero if proper boundary conditions are applied, such that a solution exists.
In this work, we adapted an approach that employs the Block-Gauss-Seidel (BGS) algorithm (see Algorithm 1, adapted from [25]).Therefore, the linear system of equations (Equation ( 13)) is split into physically meaningful sub-blocks, which, in this case, are the bulk domains (i.e., the electrodes, solid electrolyte, and current collectors) and the domain of the grain boundary Application of the BGS algorithm results in multiple smaller systems of linear equations that only require inverting the main-diagonal blocks.We know that the main diagonal block of the grain boundary domain could have a saddle-point structure.Therefore, the system of equations of this sub-block is solved with a direct solver (in this case UMFPACK [26]) while the other sub-blocks are solved using an algebraic multigrid solver (c.f.[27]) together with a row and column-based equilibration.This is affordable, as for realistic microstructures, the number of unknowns in the grain boundary block is significantly smaller (for our example in Section 3.3: 6, 200 nodes•4 dofs node = 24, 800 dofs) compared to the number of unknowns in the other domains (157, 000 nodes • 2 dofs node = 314, 000 dofs), and thus, the computationally greater amount of workload of direct solvers does not weigh heavy.Note, that several other methods are available to circumvent the solution of a saddle-point system, e.g., by making use of the Schur-complement (c.f.[24] or [28]).However, the focus of this work is not on the assessment and comparison of different solvers.

Remarks on alternative formulations
Alternative formulations exist to enforce the constraint on the electric current using the finite-element method.One possibility is to change the function space of the shape functions to Hermite shape functions (c.f.[29] for details).These shape functions include the derivative at the nodal values, which scales with the electric current, as an additional unknown.Thereby, constraints on the electric current could easily be applied.However, this change of the function space requires a fundamental reformulation of finite element codes and is often not applicable.Another possibility is given by enforcing the constraint in standard formulations using Lagrangian shape functions.However, the derivative of a quantity includes, in general, all nodes of an element, and therefore, the constraint affects not just the nodes at the interface Γ gb-el but all elemental nodes, which adds additional hurdles to the implementation.

Numerical examples
We first analyze the novel model and show the correctness of the model and its implementation for a geometrically simplified setup.Afterwards, the influence of ionically conducting grain boundaries on the effective conductivity of solid electrolytes is discussed for a simple geometry.Finally, results for grain boundaries in realistic microstructures are shown.

Analysis of the novel model
The correctness of the model is shown by a convergence study with an analytical solution as a reference and the evaluation of the conservation of charge at the intersection between three grain boundaries.

Convergence study on simple grain geometries
An analytic solution can be found for certain geometric setups and boundary conditions.This analytic solution serves as a reference for a spatial convergence analysis.Here, a three-dimensional setup without units is chosen where three planar grain boundaries are modeled between solid electrolyte grains (see Figure 4).They intersect at one line (depicted in the quasi-two-dimensional setup in Figure 4 as a point).The coordinates x 1 , x 2 , and x 3 are introduced along the grain boundaries.The electric potential in the grains is fixed to zero Φ se = 0. Thus, the equation for the conservation of charge inside the grain boundaries reduces to a one-dimensional differential equation with the second term being the exchange current between the grain boundary and the grains.The analytic solution of the electric potential Φ gb,ana is given by with the constants c 1 i , c 2 i .Consequently, the current i gb,ana is By enforcing the constraints defined in Equations (6a) and (6b) and applying Dirichlet boundary conditions for the electric potential at all ends of the T-shape (Φ 1 = 0, Φ 2 = 0.1, Φ 3 = 4, see circles in Figure 4), the constants c 1 i , c 2 i can be computed.The analytic solution of the electric potential and the electric current within the three grain boundaries Γ 1 -Γ 3 are plotted in Figure 5 for κ gb = 1, and √ κ gb t gb r n = √ 10.The positive direction of the current is defined in the positive direction of the respective coordinate.As expected, the currents sum up to zero, i.e., i 1 n 1 +i 2 n 2 +i 3 n 3 = 0, with i i denoting the current of each grain boundary Γ i at the intersection.Note that the direction of the normal vectors of the three grain boundary domains w.r.t. the coordinate has to be considered when computing the sum.A spatial convergence study is performed based on this analytic solution.Therefore, the relative L2-norm ϵ of the deviation of the numerical solution of the electric potential Φ gb from the analytic solution is computed as The development of the relative L2-norm for a decreasing edge length of the hexahedral elements used to mesh the three-dimensional geometry is shown in Figure 6 together with two lines representing linear and quadratic convergence, respectively.As expected, quadratic convergence is observed, which confirms a correct formulation and implementation of the constraints.

Conservation at non-orthogonally intersecting grain boundaries
A geometry including non-orthogonal intersections of the grain boundaries is investigated in another threedimensional example, again without units.The solid electrolyte grains are modeled as three intersecting cylindrical grains (see circles in Figure 7).The electric potential in the solid electrolyte grains is fixed to zero  Φ se = 0.At the outer ends of the grain boundaries, the electric potential is fixed (Φ 1 = 0, Φ 2 = 0.1, Φ 3 = 2, see Figure 7) to different values and the ionic conductivity of the grain boundaries is set to κ = 0.1.For this setup, the sum of the currents at the intersection of the grain boundaries can be evaluated to quantify the fulfillment of the constraints.Therefore, the relative sum of the electric currents into the grain boundaries ϵ = isum itot is computed at the intersection of the grain boundaries with i sum = i n i • i i , and the total current i tot = i |i i |.For the outlined case, the electric current is shown in Figure 8 and the relative sum is ϵ = isum itot = 6.8 • 10 −9 at the intersection, which shows the fulfillment of the charge conservation at the intersection within the expected numerical tolerances.

Assessment of the influence of the grain boundaries
An artificial geometry with regularly arranged grains, as shown in Figure 9, allows finding analytical results for extreme cases of the ionic conductivity inside of the grain boundary.These results serve as a basis to gain deeper insights and enable comparison with the solution obtained from the simulation.We analyze the difference in electric potential ∆Φ SE through the solid electrolyte, i.e., from the leftmost point of the solid electrolyte to the rightmost point, for a given current i through the solid electrolyte.The difference in electric potential is a measure of the effective conductivity of the solid electrolyte.We distinguish between two extreme cases: (1) the ionic conductivity inside the grain boundaries approaches zero, and (2) the ionic conductivity inside the grain boundaries approaches infinity.In the first extreme case, conduction inside the grain boundaries is unfavored, and therefore, the shortest conduction path is exclusively in the z-direction, i.e., it is strictly orthogonal across the grain boundaries, which are normal to the z-direction.For this case, the difference in electric potential can be estimated by a series of resistors of n SE = 12 solid electrolyte grains (r se = n se tse κse ), n gb = 11 grain boundaries (r gb = n gb t gb κ gb ), and 2n gb contact resistances (r gb-se = 2 n gb r c ) scaled by the current i In the second extreme case, conduction in the grains is unfavored.Therefore, the shortest conduction path is through the grain boundaries because their resistance approaches zero in the limit of infinite ionic conductivity.Thus, the remaining resistance originates from two solid electrolyte grains (r se = 2 tse κse ) and two contact resistances (r gb-se = 2r c ) The solution of both equations is plotted in Figure 10a in dashed lines for different values of the ionic conductivity in the grain boundary κ gb .The other parameters are constant as defined in Table 2.It can be seen that both curves are separated from each other, and especially if the focus of the investigation lies between the extreme cases, a model that resolves conduction along the grain boundary becomes vital.Vertical lines exemplarily indicate the values for the conductivity of LLTO inside the grains (dark grey) and inside the grain boundary (light grey).They highlight that the incorporation of conduction inside grain boundaries is important for realistic material parameters.The difference in electric potential of evaluations  of the novel model with different values for the ionic conductivity of the grain boundary κ gb are added to Figure 10a.It can be seen that the proposed model is able to cover both extreme cases with a smooth transition between them.This transition is within a relevant range of the ionic conductivity of typical solid electrolyte materials and grain boundaries and, therefore, requires resolving charge transport also along the grain boundary.Based on these two extreme cases, a deeper understanding of the ionic conduction in the grain boundary can be derived.Therefore, we modify single parameters from Equations (19a) and (19b) and plot the results in Figure 10b together with the results of evaluations with the novel model.It can be seen that the characteristic shape remains the same for all combinations while the magnitude and slope change: • Reduction of the contact resistance (blue curve, reduction by one order of magnitude) reduces the total resistance.This is relevant within a large range of the ionic conductivity of the grain boundaries.For small ionic conductivities of the grain boundaries, the total resistance is dominated by this small ionic conductivity, and therefore, the influence of the contact resistance becomes negligible in these regions.
• Reduction of the ionic conductivity in the solid electrolyte (red curve, reduction by two orders of magnitude) increases the total resistance.However, this is only relevant for high ionic conductivities of the grain boundary as for small ionic conductivities, the total resistance is dominated by this small ionic conductivity.
• Reduction of the current (yellow curve, reduction by factor 2.5) parallelly shifts the difference in electric potential.For this setup, the current is simply a linear amplification.
• Reduction of the thickness of the grain boundary (purple curve, reduction by one order of magnitude) reduces the total resistance for small ionic conductivities of the grain boundary while it remains the same for high ionic conductivities.We can conclude that the novel model is capable of resolving all variations, including the extreme cases.

Application of the model to a realisitc microstructure
The focus of this work is to establish a formulation to include ion conduction in intersecting grain boundaries into a continuum model for SSBs.Hence, for completeness, we show its applicability to a realistic and, therefore, complex microstructure.

Microstructure
The grains are approximated as spherical particles (see Figure 11a  cathode consists of active material particles (mean of diameter distribution µ c = 2.3, and standard deviation σ c = 0.05 [30], evaluated in µm) and solid electrolyte grains (µ el = 1.831, σ el = 0.548 [31], evaluated in µm) which follow a log-normal distribution for the diameter d with the probability density function Both phases a volumetric ratio of vc vc+v el = 0.48.The thickness of the composite cathode is 30µm.The anode is modeled as a planar lithium metal foil with a thickness of 10µm.A thin current collector foil with a thickness of 2µm is attached to both electrodes.The separator is located between the electrodes with a thickness of 20µm.It consists of the same solid electrolyte grains as in the composite cathode.We construct a planar grain boundary where two solid electrolyte grains intersect.Thereby, a connected network of grain boundaries is created (see Figure 11b).The average porosity of the porous separator and the porous composite cathode is ϵ = 0.15.The lateral edge length is 36 µm to limit the size of the computational domain.

Materials
The active material of the cathode is NMC622, the solid electrolyte is LLTO, the anode is lithium metal, and the current collectors are copper and aluminum, respectively.The material parameters are summarized in Table 3.

Boundary and initial conditions
A discharge scenario is simulated.Therefore, the initial concentration in the electrodes is chosen to represent a charged state, i.e c c,0 = 21, 000 mol m 3 and c a,0 = ρa Ma = 76, 900 mol m 3 .A constant voltage of Φ = 0 V is enforced at the current collector on the anode side, and a constant current is applied to the current collector at the cathode side, such that the battery cell is discharged with a c rate of 0.1C until the cut-off voltage between both current collectors of ∆Φ = 2.7 V is reached.

Results
The electric potential in the grain boundaries at the end of discharge is shown in Figure 12.A decrease in  preferred compared to conduction in the grains.Furthermore, it can be seen that the in-plane current does not remain constant but changes over time.This can be attributed to the inhomogeneous lithiation of the cathode, which results in an inhomogeneous electronic conductivity in the cathode (see Equation ( 27)) and an inhomogeneous distribution of the equilibrium potential at the interface between cathode and solid electrolyte and, therefore, to shifted optimal conduction paths during discharging.If, thereby, the impedance in the cathode is increased, conduction in the solid electrolyte and, hence, along the grain boundaries become more favored and vice versa.Thus, the current along the grain boundary changes over time.Additionally, the cell voltage is shown in Figure 13b together with a zoom into the in-plane current for the standard value of the ionic conductivity in the grain boundary to highlight the dependence of the state of charge of the in-plane current.
In Figure 14, the geometrically resolved magnitude of the in-plane current inside of the grain boundaries is shown for different values of the ionic conductivities at the end of the discharge in a slice close to the anode.Note that the dominant direction of the current is out of the shown plane, i.e., in the z-direction of Figure 11a.An increase in current is observed for higher conductivities (indicated by the different color bars in Figure 14).Moreover, an increased current is also observed at the intersection between more than two grain boundaries, marked with red circles in Figure 14.This is caused by the current that is merged from two grain boundaries into one grain boundary; therefore, this increase is not observed in non-intersecting grain boundaries (purple circles in Figure 14).This local increase in current could be unfavored as it could initiate the development of dendrites or local plating [32,33,34].Beyond this physical insight, this also highlights the necessity of consistent constraint enforcement at the intersections to properly resolve the increase in current there.

Summary
A modeling approach is presented to incorporate the ionic conduction along grain boundaries into a continuum model for solid-state batteries that geometrically resolves the microstructure.Based on a formulation to represent transport in thin layers of solid-state batteries [17], the grain boundaries are reduced to a twodimensional manifold.This reduction raises the question of how to guarantee the conservation of mass and charge at locations where more than two grain boundaries, modeled as two-dimensional manifolds, intersect.In terms of ionic conduction inside of the grain boundaries, this means a unique electric potential and a net current of zero.These constraints are enforced by treating the electric potential and the electric current as independent unknowns within the system of equations.This formulation comes with some numerical challenges like the solution of a saddle point system, which we discuss in this work together with solution strategies for them.We show the fulfillment of the formulated constraints as well as the convergence of the numerical error.Afterwards, we discuss extreme cases for the limit of infinite and zero ionic conductivity of the grain boundaries in terms of their influence on the voltage drop within the solid electrolyte.We observe that the novel formulation is able to cover both extreme cases as well as a smooth transition between them.Finally, we show the applicability of the model to realistic microstructures, extract the current along the grain boundaries, and find an increase in the magnitude of the flux at the intersection of grain boundaries, which could be attributed to degradation mechanisms in further studies.Furthermore, the results of impedance measurements could be classified into contributions from grain boundaries and bulk domains.The outlined formulation is exemplarily shown for ionic conduction in the grain boundaries.However, recent publications indicate that also electronic conduction [5] may occur in grain boundaries of solid-state batteries.With the model established in this work, this and also various other transport phenomena in grain boundaries can be incorporated into the model together with advanced interface kinetics to also model, e.g., lithium deposition in the grain boundaries.

Funding
We gratefully acknowledge support by the Bavarian Ministry of Economic Affairs, Regional Development and Energy [project "Industrialisierbarkeit von Festkörperelektrolytzellen"] and the German Federal Ministry of Education and Research [FestBatt 2 (03XP0435B)].

A Discretization of the equations in the grain boundary
For discretizing Equations (8a) to (8c) the unknowns Φ gb , i gb and the test functions w Φ , w i , polynomial shape functions that are organized in the matrix N , i.e., [Φ gb , i gb , w Φ , w i ] T = N Φgb , îgb , ŵΦ , ŵi T are employed.In this work, linear Lagrangian polynomials are used to form the space of the shape functions.The spatially discretized forms of the weighted residua are Rgb,std = The residual Rgb = Rgb,std + Rgb,constr Φ + Rgb,constr i is reorganized w.r. with This nonlinear, algebraic system of equations is iteratively solved for the unknowns organized in a vector Ψgb = [ Φi , Φm , Φs , îi , îm , îs , λΦ , λi ] T by the Newton-Raphson scheme The subscripts denote the interior, master, or slave side matrices.The Lagrange multipliers, as well as the slave side values, are removed from the global system of equations by static condensation, and the final system of linear equations is with the condensed vector of unknowns

B Discretization of the equations in the bulk domains
For discretizing Equation (10), the unknowns Φ bulk , c and the test functions w Φ and w c are as well discretized with polynomial shape functions N , i.e., [Φ bulk , c, w Φ , w c ] T = N Φbulk , ĉ, ŵΦ , ŵc T .Again, linear Lagrangian polynomials are used to form the space of the shape functions.
The linearized system of equations for the Newton-Raphson scheme in the bulk domains is given as

C Block Gauss-Seidel algorithm for saddle-point systems
The Block Gauss-Seidel algorithm (Algorithm 1) is used to iteratively solve the linear system of equations K∆x = R until norm (E r KE c ∆x − E r R) < ϵ, with a tolerance ϵ Thereby, only systems of equations that contain the n blocks on the main diagonal need to be solved.If this sub-system contains a saddlepoint structure, a direct solver is employed.The full system of equations is equilibrated by row (E r ) and column (E c ) multiplication to improve its condition.We chose both matrices as diagonal matrices containing the reciprocal of the largest value of the respective row or column within each sub-block.adapted for NMC622 -β-LPS from [22] and [30] interface anodesolid electrolyte Γ an-el exchange current density i 0 8.87 A m 2 ( The open circuit potential [35] of NMC622 is shown in Figure 15.Lithiation state χ Open circuit potential Φ 0 in V Figure 15: Open circuit potential of NMC622 as a function of the lithiation state based on [35].

Figure 1 :
Figure 1: SEM image of the cross-section of LLZO including grains and grain boundaries (M.Balaish, Department of Chemistry, Technical University of Munich).

Figure 2 :
Figure 2: Schematic sketch of the geometry of a battery cell including all domains Ω i , boundaries Γ i , and interfaces between domains Γ i−j .

2 Figure 3 :
Figure 3: Schematic of the intersection of three grain boundaries (black lines) with the definition of the normal vectors and exemplary slopes of the electric potential (blue line) and the electric current (red line), including the constraints at the intersection.The arc illustrates the equal electric potential at this location, and the red arrows sum up to zero, indicating the constraint on the electric current.

Figure 4 : 3 Figure 5 :
Figure 4: Geometry for the convergence study.The grain boundaries occur at the intersection of the solid electrolyte domains (colored rectangles).The edge size of each grain is a = 4 or 2a.The coordinates x i are defined along the grain boundaries.The electric potential is set as a boundary condition at the ends of the grain boundaries marked with a circle.

Figure 6 :
Figure 6: Convergence of the relative L2-norm (blue line).Model evaluations are marked with a cross.The dashed lines represent linear and quadratic convergence rates, respectively.

Figure 7 :
Figure 7: Geometry with non-orthogonally intersecting grain boundaries.The grain boundaries occur at the intersection of the solid electrolyte domains (colored domains with a radius of r = 1).The electric potential is set as a boundary condition at the ends of the grain boundaries marked with circles.

Figure 8 :
Figure 8: Current in the grain boundary.The jump at the intersection, s.t. the net current sums up to zero, is visible.Note that the geometry is distorted for visualization compared to Figure 7.

10 − 8 10 −6 10 − 4 10 −2 10 0 10 2 10 −1 10 0
ionic conductivity in S m diff. in el.pot. in solid electrolyte in V (a) The extreme cases (infinite conductivity and zero conductivity in the grain boundaries) are shown by dashed lines, and evaluations of the proposed model are represented by the green line.The ionic conductivity of the grains (dark grey) and the grain boundary (light grey) are indicated by vertical lines for LLTO[18].10−8 10 −6 10 −4 10 −2 10 0 10 2 el.pot. in solid electrolyte in V (b) Variation of parameters that influence the extreme cases.The dashed black line represents the default values and the colored lines the variations of the parameters in the following manner: blue: rc < rc 0 , red: κse < κse 0 , yellow: i < i0, purple: t gb < t gb 0 , where the index 0 represents the default value.The evaluations with the novel model are symbolized by the dotted lines.

Figure 10 :
Figure 10: Evaluation of the influence of the ionic conductivity of the grain boundary on the voltage drop in the solid electrolyte.
, blue: solid electrolyte, grey: cathode active material, red: grain boundaries, silver: anode, silver/brown: current collectors).The composite Battery cell with porous separator and composite electrode.(b)Network of grain boundaries.

Figure 11 :
Figure 11: Geometry of battery cell with grain boundaries.The grain boundaries occur at the intersections of the solid electrolyte grains and are modeled as planar surfaces.

Figure 12 : 2 ( 2 (
Figure 12: Electric potential in the grain boundaries at the end of discharge.

Figure 13 :
Figure 13: In-plane current over time.

Figure 14 :
Figure 14: Electric current inside of slices through grain boundaries for different values of the ionic conductivity of the grain boundary (from left to right: small to large, κ gb = 1.88 • 10 −4 , 10 −2 , 1 S m ) at the end discharge.Note the different limits of the color bar that indicate the different orders of magnitude of the current for the different ionic conductivities.The red circles indicate the intersection between three grain boundaries.The purple circle indicates non-intersecting grain boundaries.

Table 1 :
List of symbols.

Table 2 :
Parameters for the setup with the regular grains.
d∂Γ, and G = Γ gb κ gb N T ∇ Γ N dΓ.The test functions ŵ have arbitrary values, s.t. each term has to be individually zero, i.e.
Algorithm 1 Block Gauss-Seidel with saddle-point structureA ← E r KE c ▷ Equilibration c ← E r R y i ← 0 ∀ i ∈ n h ← 0 while r > ϵ do h ← h + 1 for i ∈ {1, . . ., n} do ▷ Loop over main-diagonal blocks g ← c i for j ∈ {1, . . ., n} do ▷ Loop over off-diagonal blocks if j < i then g ← g − A ij y h−1 j else if j > i then g ← g − A ij y h