Simulating hindered grain boundary diffusion using the smoothed boundary method

Grain boundaries can greatly affect the transport properties of polycrystalline materials, particularly when the grain size approaches the nanoscale. While grain boundaries often enhance diffusion by providing a fast pathway for chemical transport, some material systems, such as those of solid oxide fuel cells and battery cathode particles, exhibit the opposite behavior, where grain boundaries act to hinder diffusion. To facilitate the study of systems with hindered grain boundary diffusion, we propose a model that utilizes the smoothed boundary method to simulate the dynamic concentration evolution in polycrystalline systems. The model employs domain parameters with diffuse interfaces to describe the grains, thereby enabling solutions with explicit consideration of their complex geometries. The intrinsic error arising from the diffuse interface approach employed in our proposed model is explored by comparing the results against a sharp interface model for a variety of parameter sets. Finally, two case studies are considered to demonstrate potential applications of the model. First, a nanocrystalline yttria-stabilized zirconia solid oxide fuel cell system is investigated, and the effective diffusivities are extracted from the simulation results and are compared to the values obtained through mean-field approximations. Second, the concentration evolution during lithiation of a polycrystalline battery cathode particle is simulated to demonstrate the method’s capability.

Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence.Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

Introduction
Defects, such as grain boundaries in polycrystalline solids, often have a significant effect on material properties.Grains and grain boundaries in materials arise from various processing techniques such as sintering, agglomeration of particles, or nucleation and growth.The unique configurations of atoms at grain boundaries, in which bonding structures differ from atoms in the intrinsic lattice, lead to differences in mechanical and/or transport behavior between grain boundaries and bulk regions.In many cases, materials can be designed to take advantage of these grain boundary effects [1], such as improving corrosion resistance [2,3], creep strength [3,4], ductility [5,6], and damage resistance [3,7,8] in metallic alloys.Grain boundaries have also been shown to increase transport rates of charge-carrying species in photovoltaics [9][10][11] and battery materials [12,13] as well as enhancing dopant diffusion rates in polycrystalline silicon [14].
However, some materials systems exhibit the opposite behavior, in which grain boundaries act to hinder transport.In the solid oxide fuel cell material yttria-stabilized zirconia (YSZ), space charge layers and oxygen vacancy depletion at grain boundaries decrease the diffusivity of oxygen ions at the grain boundaries compared to the bulk [15][16][17].The decrease in diffusivity at grain boundaries becomes particularly important for performance in nanocrystalline YSZ, in which the effects of the grain boundaries become more prominent as the grain boundary density increases [18].Battery cathode particles often consist of agglomerations of smaller particles that form dense, polycrystalline structures [19][20][21].During electrochemical charging and discharging, stress from volumetric changes may cause decohesion and cracks to form at the grain boundaries of the cathode particle which may obstruct ionic transport (if electrolyte does not penetrate) and negatively affect battery performance [19].Hydrogen diffusion in aluminum [22] and helium diffusion in tungsten [23] are hindered across grain boundaries.Hindered grain boundary diffusion has been simulated at the atomic scale in bicrystals via molecular dynamics [18,24] and the kinetic Monte Carlo method [25].Additionally, there are instances of hindrance relevant to the biomedical field; such is the case for medicine delivery to tumors where diffusion may be impeded at membranes in vasculature [26].
At the microstructure scale, recent studies have considered both hindered [27][28][29] and enhanced [27][28][29][30][31] grain boundary diffusion.In [30], a diffuse interface method was developed to simulate enhanced grain boundary diffusion, as discussed further below.In [27,28], both hindered and enhanced diffusion of hydrogen were considered in polycrystalline metals as part of parametric studies that also included hydrogen trapping.In [29], lithium in polycrystalline lithium lanthanum zirconium oxide was modeled to have hindered grain boundary diffusion at low temperatures and enhanced grain boundary diffusion at high temperatures, which was predicted by their molecular dynamics simulations.The effect of grain boundary diffusivity was modeled in [27][28][29] by solving a diffusion equation with diffusivity that is a function of phase-field order parameters that specify the grains.A grain boundary is represented implicitly by the diffuse interface where the order parameters of adjacent grains change in value.In the models presented in [27][28][29], while the diffuse interface enables accurate modeling with a simple (e.g.uniformly meshed) spatial discretization, its thickness in the simulation must be similar in magnitude to the physical width of the grain boundary, limiting the grain sizes of the microstructures that can be investigated [28].
In this paper, we follow previous work on enhanced grain boundary diffusion [30] by proposing a new diffuse interface approach based on the smoothed boundary method (SBM) to solve the diffusion equation with hindrance at the interfaces between grains.SBM is a method for applying boundary conditions to irregularly shaped domains described by a diffuse interface [32].By treating the hindrance as a boundary condition, the diffuse interface in our simulations can be decoupled from the physical width of the grain boundary, which becomes a model parameter.The proposed model is compared to a sharp interface model in one dimension.The effect of the numerical and physical parameters on the error of the model is explored.Finally, two case studies are presented.The first case study is oxygen diffusion in nanocrystalline YSZ, in which the effective diffusivities of the inhomogeneous structures are calculated and compared to the values obtained from mean field approximations.The second case study is the lithiation of a polycrystalline battery cathode particle that consists of primary particles (grains).
The proposed model will provide the dynamic, non-uniform concentration field that sets composition-dependent material properties and thus influences physical processes, such as reaction kinetics of electrochemically active materials or the stress state of a composite electrode.Additionally, the concentration information can be used to calculate effective properties of a microstructure.For example, the effective diffusivity of a polycrystalline microstructure can be determined from the flux in the system in which the grain boundaries are explicitly included without making simplifying assumptions required in mean field approximations such as the Hart [33] or Maxwell Garnett [34] approaches.

The sharp interface model
Grain boundaries that hinder diffusion have their greatest effect at orientations perpendicular to the diffusion direction as that configuration positions the high-and low-diffusivity regions in series.Therefore, as a first step to develop a numerical model and assess its accuracy, we consider a one-dimensional system in which grain boundaries are placed perpendicular to the diffusion direction (i.e. the grain boundary normal is parallel to the diffusion direction).The flux in concentration C in the bulk of the grains is described in one dimension by Fick's first law, where x is the spatial coordinate and D bulk is the bulk diffusivity, which we assume to be constant.The evolution of the concentration in time is given by the diffusion equation, which, in one dimension with constant diffusivity, can be expressed as where t is time.
To account for the hindrance of diffusion at the grain boundary, a Neumann boundary condition is imposed at each grain boundary such that for an identical concentration gradient, the flux across the grain boundary is lower than the flux in the bulk, resulting in a sharper concentration gradient in grain boundaries than in the bulk at steady state.The other physical constants needed to parameterize the model for a given system are the grain boundary diffusivity, D gb , the grain boundary thickness , δ, and the grain size, d.
First, a sharp interface model for the flux boundary condition at grain boundaries will be described and used as a baseline.The flux at the grain boundary, J gb , is dependent on the concentration drop across the grain boundary and is given by: where ∆C gb is the concentration drop across the grain boundary, and κ is the degree of hindrance.
The degree of hindrance, κ, is the parameter that describes how strongly the grain boundary blocks chemical transport, and it is related to the diffusivity and thickness of the grain boundary.These properties may vary with grain boundary misorientation and have different values for high-symmetry grain boundaries such as twins or Σ-boundaries.They can also be altered at triple junctions, which can be identified by a triple product of the grain order parameters.In this work, however, for simplicity we will consider them to be constants for a given material.The assumption of constant grain boundary properties has been suggested as a reasonable approximation for sintered ceramic materials where high-angle grain boundaries predominate [29].To formulate κ, we consider diffusivity across grain bulk and grain boundaries as an analogy to electrical conductivity in a series circuit.Conductivity is a measure of how easily electricity can flow through a material whereas diffusivity across a grain boundary refers to how easily chemical species can flow across the grain boundary.Using this analogy, we formulate κ such that the degree of hindrance is the difference between the inverse of D gb and D bulk multiplied by δ such that the degree of hindrance becomes independent of δ when D gb = D bulk as expected.Therefore, our formulation of κ is The second term counteracts the constant bulk diffusivity assumed everywhere.The boundary conditions at the left and right boundaries of the domain are Dirichlet boundary conditions of C = 1 and C = 0, respectively.The concentration evolution is calculated using equation ( 2).This one-dimensional sharp interface model is sufficient for systems in which the grain boundaries are all perfectly perpendicular to the diffusion direction; however, a model that can allow for angled or curved grain boundaries is necessary for more complex systems.
We implement the sharp interface model using a one-dimensional finite difference discretization that will share the same time discretization (forward Euler), grid spacing (∆x), and time step (∆t = 0.05∆x 2 /D bulk ) as the SBM model described in the next subsection.The sharp interface model requires two points at each grain boundary location (one on each side of the grain boundary) for accurate evaluation of J gb in equation ( 3).The bulk flux in equation ( 1) is calculated on half points using centered differences, (J bulk , where i is the grid point index.At points in the bulk grains, the finite difference update is where n is the time step index, which is second-order accurate in ∆x [35].At grain boundary points, a second-order-accurate one-sided difference is used.For a point with index i on the left side on the grain boundary this is . Similarly, the formula for the point on the right side is

The SBM Model
SBM is a diffuse interface approach for setting internal boundary conditions using a field variable called a domain parameter [32].For a polycrystal, we use domain parameters, ϕ q , to describe the location of each grain in the system such that ϕ q = 0 outside Grain q and ϕ q = 1 inside Grain q, as in the phase-field modeling of grain growth [36].The interface, which represents a grain boundary, is indicated as a region of smooth transition between 0 and 1 described by a hyperbolic tangent function in one dimension, where x q is the location of the grain boundary and ζ is a numerical parameter that controls the thickness of the diffuse interface.The hyperbolic tangent profile for interfaces naturally arises as a solution to the Allen-Cahn equation [37].At a boundary between any two grains, we expect where Q is the number of grains in the system (behavior of phase-field models at junctions between three or more grains is more complex [38] and will be neglected here).Phase-field models enable the definition of grain inclination and misorientation using the gradients of the grain order parameters [36].Thus, while we here consider κ to be independent of orientation, our approach could be extended to include misorientation-dependent hindrance.
To include the effect of grain boundaries, the diffusion equation is modified to include a term that sets the Neumann boundary condition at the grain boundaries (i.e.regions where the gradient of the domain parameter is nonzero) [32]: where the flux in the bulk of the grains J SBM bulk is given by in which C q is the concentration field for grain q.See [32] for the derivation of equation ( 6).
The grain boundary flux J SBM gb takes the same form as J gb in equation ( 3), but now a separate concentration field must be solved for each grain to accurately capture ∆C gb .This concentration drop for the diffuse interface approach is calculated by such that ∆C SBM gb is the difference between the concentration of the grain considered and a weighted average of the concentrations of all the other grains.The use of this weighted average term also prevents any spurious effect from the concentration solution outside a grain located far from the grain under consideration.
For the one-dimensional simulations, we choose a domain which consists of three grains and two grain boundaries (figure 1(a)) to evaluate the error resulting from the interaction between grain boundaries.The first and last grain in the domain have lengths half that of the other grain(s).This choice is made based on the symmetry it provides, as it allows the bulk-tointerface ratio to remain constant regardless of the number of grains in the system.To illustrate the effect of grain boundary hindrance in this test system, the simulated concentration C = Q q=1 C q ϕ q is shown at steady state in figure 1(b) for three different values of the ratio between grain boundary and bulk diffusivity, D gb /D bulk .Lower ratios D gb /D bulk lead to larger gradients in concentration across the diffuse grain boundaries and smaller gradients in the bulk grains.
The SBM model was discretized using finite differences, namely a forward Euler tine stepping scheme and centered differences in space on a uniform grid.The bulk flux is calculated at half points as , where q denotes the grain index and i denotes the grid point index.The spatial derivative of the bulk flux is calculated at grid points using a centered difference of J SBM bulk values at half points (i + 1 2 , etc) and the term |∂ϕ q /∂x| is calculated at grid points using a centered difference of ϕ q values at full points (i, etc).To reduce the effect of rounding error in the comparison, all parameters were nondimensionalized for simulations using a length scale W equal to the grain size and a time scale τ = W 2 /D bulk , such that non-dimensional bulk diffusivity was equal to unity.Concentrations are scaled to fall within the range of zero to one.
For the three-dimensional structures, W was set to be equal to the smallest grain size (40 nm for YSZ and 100 nm for battery cathode primary particles).Two YSZ geometries are considered: a representative bulk polycrystalline microstructure consisting of a cube with side length 240 nm and a thin-film specimen with dimensions 240 nm × 240 nm × 40 nm.The numerical grid sizes for the YSZ simulations were 120 × 120 × 120 and 120 × 120 × 20 grid points for the cube structure and the thin-film structure, respectively.The cathode particle simulation consisted of a secondary (polycrystalline) particle with nominal diameter 1 µm inside a cubic simulation domain having a side length 1.2 µm that was discretized using 240 × 240 × 240 grid points.Before smoothing, grains in the cathode-particle domain with their center of mass located outside the sphere of diameter 1 µm were removed, leaving a rough surface similar to those of experimentally synthesized particles (e.g.see [19]).The particle is described by an additional domain parameter ψ = ϕ q with which the Dirichlet concentration boundary condition is set on the outside surface of the particle via SBM [32].The three-dimensional polycrystalline microstructures were originally generated using Dream.3D[39], and then the order parameters were smoothed and reduced following the description in [30].
The numerical parameter ratios are ζ/∆x (a measure of the resolution of the numerical interface) and N gb ζ/L x (a measure of the numerical interface volume fraction or the thickness of the numerical interface), where L x is the domain size and N gb is the number of grain boundaries.The physical parameter ratios are D gb /D bulk and V gb = N gb δ/L x , the latter of which is the physical grain boundary volume fraction.

Error analysis
The sharp and diffuse interface models were compared by examining one-dimensional simulations on uniform grids with the same values of ∆x and with finite difference approximations that have the same order of accuracy (second order with respect to ∆x and first order with respect to ∆t).By making the comparisons between the numerical solutions solved by consistent methods and parameters, we ensure that the error represents the effect of the diffuse interface, rather than a difference in accuracy of the numerical method and discretization.The metric chosen with which to calculate error is the average difference in bulk flux between the two methods, We consider the evolution of this error over the course of a simulation, and we compare sets of parameters based on the maximum value it reaches, ϵ max .

One-dimensional simulations and error analyses
The solutions from the sharp-interface model, equations ( 1)-( 3), and our model, equations ( 6)-( 8) are obtained and compared to assess the error introduced by the diffuse interfaces in our model.An initial parameter set of ζ/∆x = 1, N gb ζ/L x = 1/48, D gb /D bulk = 1/50 = 0.02, and V gb = 1/40 was chosen.The resolution of the numerical interface, ζ/∆x, corresponds to three points across the interface, defined approximately by the range 0.1 ⩽ ϕ ⩽ 0.9.
Figure 2(a) shows the transient behavior of the error over the course of this simulation.Initially zero, the error rapidly reaches a sharp initial peak, followed by a second, shallow peak, followed by a decline to near zero corresponding to steady state.(Note that ϵ = 0 here means that the simulations have identical fluxes; both fluxes still have truncation error associated with the finite difference discretization.)Figure 2(b) compares the sharp and SBM concentration profiles at the time corresponding to the initial error peak, t = 0.0543 τ .At this time, most of the gradient in concentration is still in the left-most grain and the first grain boundary, and the second grain is just beginning to increase to the steady-state profile, which is shown in figure 2(c).This time corresponds to the rapid evolution in concentration (concentration pileup at the hindering grain boundary) in the vicinity of the left-most grain boundary, which we hypothesize leads to the peak in error.The same explanation applies to the second peak, which corresponds to concentration pile up at the second grain boundary that is fast relative other times but still slower than at the first/left-most grain boundary (thus resulting in a smaller peak value).In additional simulations, we found that ϵ max was unaffected by changing the number of grains N gb when the numerical interface volume fraction N gb ζ/L x was held constant, supporting the idea that ϵ max corresponds to an interaction with a single grain boundary.
Figure 3 shows how ϵ max changes as a function of the numerical interface volume fraction N gb ζ/L x for three values of the ratio of diffusivities D gb /D bulk .Other parameters are held   (10).Data and fits for three values of D gb /D bulk are shown, namely, 0.004 (black circles, black line), 0.02 (blue squares, blue line), and 0.1 (red triangles, red line).Values for a, the fitted parameter in equation (10), and R 2 , the coefficient of determination, are displayed for each fit in the legend.constant from the simulation in figure 2. We find that ϵ max is second order in N gb ζ/L x , and the solid lines in figure 3 depict fits to the form where values of the fitting parameter a and the coefficients of determination R 2 are provided in the legend.The values of R 2 are close to unity, indicating good agreement with equation (10), which is consistent with an analytical prediction for SBM with Neumann-type boundary conditions [40].We find that the fitting parameter a increases as the ratio D gb /D bulk decreases.
Recalling figure 1(b), the ratio D gb /D bulk affects how much the gradient in concentration is localized to the grain boundaries rather than distributed through the bulk of the grains.This suggests that large gradients in concentration at the diffuse grain boundaries also affect ϵ.
The second-order relationship between error and N gb ζ/L x is predicted to hold for complex geometries in two or three dimensions [40], although the parameter a may change.Furthermore, the steady-state error, which was near zero (∼10 −12 ) in our one-dimensional simulations, is likely to be more significant for complex geometries due to interfacial curvature.Because we compare between simulations rather than to an analytical solution, our analysis does not consider error due to the finite difference discretization.Due to our use of centered differences and a forward Euler time discretization, we expect error from this source to vary proportionally to ∆x 2 and ∆t when other parameters are held constant [35].In simulations with complex geometries in two or three dimensions, care should be taken to ensure that the diffuse interface width ζ is small relative to the radii of curvature of microstructural features [32] and that ∆x/ζ is small enough to limit discretization error and avoid mesh anisotropy.

Case study 1-solid oxide fuel cell
To use our model to accurately investigate the transport properties of solid oxide fuel cell materials such as YSZ, the model must be extended to three dimensions.The generalization of equations ( 6) and (7) gives A cubic microstructure was generated with an average grain size of 50 nm with a range of 40-60 nm (comparable to those found in [15]) and an edge size of 240 nm in each direction.The microstructure was generated and smoothed in the manner presented in [30] to obtain domain parameter interfaces with hyperbolic tangent function profiles with thicknesses and grid spacing corresponding to N gb ζ/L x = 1/25 and ζ/∆x = 1, respectively.A thin-film structure having dimensions of 240 nm × 240 nm × 40 nm was similarly generated using the same set of parameters other than the thickness of the computational domain.Other physical parameters used for YSZ are D bulk = 2.17 × 10 −13 m 2 s −1 , D gb = 5.0 × 10 −15 m 2 s −1 , and δ = 1 nm [15].This parameter set corresponds to an expected error of less than 0.6% from figure 3.
Figure 4(a) shows the microstructure through plotting the sum of the square of the domain parameters.Figures 4(b)-(d) show the evolution of the concentration profile in the bulk of the YSZ.Sharp drops in concentration can be observed at the grain boundaries with a much more gradual concentration drop across the bulk of the grains.As time progresses and the system evolves towards steady state, the concentration profile becomes more diffuse in the direction of the concentration gradient for both the 40 nm and 240 nm thick YSZ structures.
One additional application of this hindered grain boundary diffusion model is to use the computed steady-state concentration profile to calculate the effective diffusivity of the microstructure using the expression where C in and C out are the concentration values for the Dirichlet boundary conditions at the left and right faces of the domain, respectively (C in = 1 and C out = 0), and J SBM is the simulated steady-state flux in the primary diffusion direction (pointing from left to right) averaged over the left or right face.For the steady-state concentration profile in figure 4(d), we find D eff = 1.17 × 10 −13 m 2 s −1 , whereas the experimentally reported value for effective diffusivity is 6.8 × 10 −14 m 2 s −1 [15], which is 42% less than the computed value.The discrepancy can be attributed to the uncertainty of the effective thickness of the grain boundary as the space charge layer thickness is greater than that of the grain boundary itself [15].Additionally, there is uncertainty of the actual particle size distribution because only the average grain size of 50 nm was reported in [15], and the 40-60 nm range was arbitrarily chosen when generating the microstructure.
There are multiple other models utilized to predict the effective diffusivity of polycrystalline solids in which grain boundary diffusivities differ from the bulk diffusivity.One such model is Hart's equation [33], which assumes a simplified geometry where all grain boundaries are parallel to the diffusion direction: In a three-dimensional structure, the grain boundary volume fraction, V gb , is calculated by multiplying δ by the total area of the grain boundaries as calculated by summing the area of triangular patches generated from a MATLAB isosurface of the domain parameters.Another mean field approach is Maxwell Garnett's equation [34], which assumes spherical grain boundaries: For the 240 nm thick YSZ microstructure, D Hart = 2.05 × 10 −13 m 2 s −1 and D MG = 1.16 × 10 −13 m 2 s −1 .The Maxwell Garnett prediction is very close to that of the numerical model because the spherical grain assumption is reasonable for an isotropic structure.The Hart model's assumption of parallel grain boundaries is not accurate for such a microstructure and overestimates the effective diffusivity.
Not only is the presented numerical model more comprehensively applicable than the existing mean field approaches (because it can dynamically predict the concentration profile), but it is also more robust when being used only to predict the effective diffusivity as it makes no assumptions about the geometry of grain boundaries.In the 240 nm thick microstructure, Maxwell Garnett's equation gives nearly the same result as the numerical model without the need-nor the cost-of a simulation.However, Maxwell Garnett's equation is only successful because the structure is isotropic.If we introduce an anisotropic microstructure in which grains are elongated in one direction, we expect the numerical model to become much more accurate than Maxwell Garnett's prediction.Furthermore, Maxwell Garnett's equation assumes that all grain boundaries possess the same properties.As mentioned earlier, our model is capable of relaxing such an assumption.
Another method to introduce anisotropy in the grain boundary structure is to consider a thin film with a thickness on the order of the grain size.In this case, there will be relatively fewer grain boundaries perpendicular to the diffusion direction.Figures 4(e)-(h) show the microstructure and concentration evolution of a 40 nm thick YSZ structure using the same average grain size of 50 nm.For this 40 nm thin film, we find the effective diffusivity value of D eff = 1.66 × 10 −13 m 2 s −1 from the simulation, while D Hart = 2.08 × 10 −13 m 2 s −1 and D MG = 1.32 × 10 −13 m 2 s −1 .The decrease in grain boundary volume fraction from the 240 nm to 40 nm structure (5.7% to 4.3%) results in an increase in D MG by 14%.However, the numerical model, which accounts for both the change in volume fraction and the introduced anisotropy, predicts an increase in D eff of 42%.Hart's equation predicts very little change between the 240 nm and 40 nm structures as the assumption of parallel grain boundaries is already at the extreme limit of anisotropy and the change in volume fraction matters little in the parallel orientation.It is evident that neither Hart's model nor Maxwell Garnett's model yield quantitatively accurate approximation in this case, indicating a need for numerical simulations.

Case study 2-cathode particle
Battery cathode particles often consist of an agglomeration of smaller, primary particles and thus are roughly spherical polycrystalline structures [19].Our numerical model for hindered grain boundary diffusion can be applied to this case, where cracks that form at the grain boundaries during cycling are the sources of diffusion hindrance.The model is parameterized for the nickel manganese cobalt (NMC) cathode material with D bulk = 6.1 × 10 −13 m 2 s −1 [21], an average grain size of 150 nm (with a range of 100 nm to 200 nm), and a cathode particle size of 1 micron [30].Because the width and effective diffusivity across the crack (or grain boundary with poor contact) are not well studied, an arbitrary choice is made for the degree of hindrance (κ = 1.0 × 10 11 s m −1 ).This parameter could be tuned and the simulation results compared against experimental data to extract an accurate κ.The simulation numerical parameters are the same as those used in the YSZ case study in section 3.2.
A boundary condition for lithium concentration (in terms of the site fraction) of C = 1 is set at the surface of the particle that is initially devoid of lithium.Figures 5(b) and (c) show the evolution of the concentration as the cathode particle is intercalated with lithium.Again, sharp concentration drops can be seen at the grain boundaries where the cracks (or grain boundaries with poor contact) are inhibiting the lithium transport.Potential uses of this numerical model for cathode particles are to investigate the cracks' effects on cycling performance and/or coupling the concentration field to a mechanical solver to study the dynamic concentration's effects on stress and crack growth.Such simulations have been demonstrated for the case where a liquid electrolyte infiltrates cracks [41], enhancing transport, but not for the case where cracks hinder transport, which is more likely for all-solid-state batteries where the solid electrolyte cannot flow into newly opened cracks [42].Future work could be made more accurate by implementing a variable diffusivity which depends on the local concentration as lithium diffusivity in NMC is known to be concentration dependent [21].

Conclusion
In summary, we developed a model that utilizes the SBM and accurately predicts the concentration evolution in polycrystalline structures with hindered grain boundary diffusion.The model developed in this work can explicitly consider geometries of grain boundaries, making it more robust than mean field approaches when calculating effective diffusivity.Furthermore, when diffusion plays a role in coupled-physics phenomena, the model can account for the details of the microstructure more accurately.For example, it can be applied to study the effects of grain morphology on chemical transport in systems such as solid oxide fuel cell and battery cathode materials, accounting for the interplay between reaction kinetics, mechanical response, and transport.Because the diffuse interface width in our model is decoupled from the width of the physical grain boundary, users can change the diffuse interface width to trade off numerical accuracy with computational cost.The error analyses conducted in this work will enable users to select proper parameter sets to accurately predict the diffusion behavior of polycrystalline systems with hindered grain boundary diffusion.Along with our previous work that focused on enhanced grain boundary diffusion, this work provides a framework for examining diffusional transport in polycrystalline systems.

Figure 1 .
Figure 1.One-dimensional SBM simulations.(a) Plots of the domain parameter, ϕ , for three grains with two grain boundaries as a function of position, x.(b) Plots of the simulated concentration field C for three different ratios of grain boundary diffusivity to bulk diffusivity, 0.004 (black line), 0.02 (blue line), and 0.1 (red line).

Figure 2 .
Figure 2. Comparison of a sharp interface simulation and an SBM simulation in one dimension.(a) Normalized error ϵ defined in equation (9) plotted vs. time.(b) Concentration profiles of the sharp interface and SBM simulations at time t = 0.0543 τ corresponding to the first peak in (a).(c) Concentration profiles corresponding to steady state.In (b) and (c), the SBM result is indicated by a solid black line and the sharp interface result is indicated by a dashed red line.

Figure 3 .
Figure 3. Calculated maximum error ϵmax in the one-dimensional simulations (indicby symbols) plotted as a function of the volume fraction of numerical interface N gb ζ/Lx, alongside their fits (indicated by solid lines) to the second-order polynomial in equation(10).Data and fits for three values of D gb /D bulk are shown, namely, 0.004 (black circles, black line), 0.02 (blue squares, blue line), and 0.1 (red triangles, red line).Values for a, the fitted parameter in equation(10), and R 2 , the coefficient of determination, are displayed for each fit in the legend.

Figure 4 .
Figure 4.The microstructure and dynamic concentration profiles for YSZ.The sum of the square of the order parameters, ∑ q ϕ 2 q , is plotted to show the grain boundaries of (a) the 240 nm thick YSZ structure and (e) the 40 nm thick YSZ structure.The evolving concentration field from the SBM simulation, C, is shown for (b) an early time, 250τ = 1.84 s, (c) a later time, 1000τ = 7.37 s, and (d) at steady state for the 240 nm YSZ structure and (f) an early time, 15τ = 0.111 s, (g) a later time, 40τ = 0.294 s, and (h) at steady state for the 40 nm YSZ structure.Each plot has one quarter of the domain removed to show the interior.

Figure 5 .
Figure 5.The microstructure and dynamic concentration profiles for a battery cathode NMC particle.(a) The sum of the square of the order parameters, ∑ q ϕ 2 q , showing the grain boundaries of NMC particle.The concentration, C, profiles are shown for (b) an early time, 150τ = 2.46 s, (c) a later time, 750τ = 12.3 s, and (d) at a time approaching complete lithiation, 1000τ = 49.2 s.Each plot has one quarter of the domain removed to see the interior.