Simulation of Impedance Spectra for a Full Three-Dimensional Ceramic Microstructure Using a Finite Element Model

A method of characterizing electrically heterogeneous electroceramics for a full three-dimensional collection of randomly shaped grains is presented. Finite element modeling, solving Maxwell's equations in space and time is used to simulate impedance spectroscopy (IS) data. This technique overcomes several deficiencies associated with previous methods used to simulate IS data and allows comprehensive treatment of a full three-dimensional granular representation of ceramic microstructure without the requirement for equivalent circuits based on the Brickwork layer model (BLM) or the introduction of constant phase elements to describe any nonideality of the IS response. This is applied to a full three-dimensional ceramic microstructure with varying grain size and electrical properties to generate IS plots that highlight limitations of the BLM in data analysis.


I. Introduction and Background
I MPEDANCE spectroscopy (IS) is widely employed to deconvolute the intrinsic (bulk) and/or extrinsic (grain boundary, electrode effects, etc.) contributions to the electrical properties of electroceramics by measuring the impedance response over a frequency spectrum, 1,2 commonly from mHz to MHz. Since the late 1960s, extracting information from IS data, such as grain core (bulk) and grain-boundary capacitance and resistance values, has been done using some form of appropriate equivalent electrical circuit. 3 This usually consists of an arrangement of resistors and capacitors, connected in series and/or in parallel to model the IS response of the polycrystalline ceramic under investigation and to provide insight into the intrinsic and extrinsic properties. Identification of the correct form of the equivalent circuit is required for meaningful analysis of the system. 1,4 This is based on the likely physical processes that occur in the material and often requires some level of intuition. In many cases, to a first approximation, the grain-core (bulk) response is described in an equivalent circuit by a parallel combination of a resistor and a capacitor (RC). This combination results in an ideal arc in the complex impedance and electric modulus plane plots, Z* and M*, respectively, and an ideal Debye peak in spectroscopic plots of the imaginary components of impedance, Z″, and electric modulus, M″, with a full-width half maximum (FWHM) of 1.14 decades on a logarithmic frequency scale. 5 Due to heterogeneities associated with defects, impurities, and complex conduction processes, such an ideal response for the grain core is seldom obtained. This leads to a nonideal Debye-like response, i.e., a depressed arc in Z* and M* plots and a nonideal Debye response in Z″ and M″ spectra with a FWHM >1.14 decades. Such responses cannot be treated accurately using a simple RC circuit 4,6 and normally requires the addition of a constant phase element (CPE) to the equivalent circuit.
One of the first attempts to correlate the microstructure of an electroceramic with a mathematical combination of resistors and capacitors was proposed in the late 1960s by Bauerle. 3 This attributed the response of the ceramic to two parallel RC elements connected in series; one assigned to the grain core and the other to the grain boundary. This successfully represented an ion-conducting ceramic and modeled the dual-arc Z* plots obtained from experimental IS data. This simple model was further developed into a three-dimensional resistive boundary layer model 7 and later into the wellknown brick layer model (BLM) in the early 1980s. 8,9 The BLM is a general representation of a ceramic using the analogy of bricks surrounded by mortar to represent grain cores surrounded by grain boundaries. Nafe 10 in the mid-1980s developed this further to allow the possibility of current flow around the grain core, through the grain core or a combination of both, by summing pathways (where appropriate) in parallel. Using these approximations it is possible to convert bulk data such as resistances and capacitances into intrinsic material properties such as conductivity and permittivity for the grain core. However, due to the unknown geometry of the grain boundaries, this method is generally considered to be unreliable to extract grain-boundary conductivity and permittivity values.
The BLM method has also been incorporated into a finite difference pixel-based simulation to calculate the current distribution. [11][12][13] Here, a pixel consists of six orthogonal nested cubes, each being assigned an RC element with the properties of a grain core or grain boundary. These pixels form the points on which the conduction path can be calculated with the nested cubes allowing a 3D interconnectivity of the microstructure to be constructed that the previous BLM methods could not provide. This not only allows the treatment of current pathways, thus replicating IS data, but also permits the BLM to be used for grain-core volume fractions from zero to unity with no breakdown of the calculation. These BLM methods, however, all have two intrinsic limitations. First, they simulate grains as cubes or regular shapes. Studies by Kidner et al. 13 varying the imposed shape of simulated grains have shown that the cubic grain approximation is only applicable to micrometer-sized grains in ceramics and is no longer valid for nano-sized grains. Second, the pathway the models predict through the sample is dependent on the nested cube connections. As this is limited to six per pixel it cannot fully represent the complex conduction paths that are possible in a ceramics with irregular grain shapes.
An alternative approach to simulate the electrical response of an electroceramic is to use effective medium theory. 14,15 This is based on Maxwell's concept of an effective medium describing the ceramic as a collection of similarly shaped, coated spheres. Each one represents a grain core, which is coated with a shell to describe the grain boundary. The spheres are packed, either filling or partially filling an effective medium, which is then given the material properties of the grain core and grain boundary. The system can then be solved for the conductivity as a function of volume fraction of the grain cores. This method has been successfully used to determine the electrical response of heterogeneous ceramics and to extract values for grain-core and grain-boundary conductivity and permittivity. One major drawback is the model does not resemble the real microstructure of a ceramic. 13 Simulating the IS response of an electroceramic using finite element modeling can overcome the deficiencies associated with the methods discussed earlier. The finite element method (FEM) is a powerful tool widely used for numerical modeling and simulation in many areas of science and engineering. The idea of using the FEM to model IS data is not new, and has been successful in describing highly resistive grain boundaries [16][17][18] ; however, this approach has been limited to two-dimensional models and a comprehensive treatment of granular three-dimensional samples is still lacking.
Here, we present a finite element package, developed inhouse to simulate IS data for electroceramics using realistic microstructures. We then apply this FEM to various electrical microstructures to simulate IS spectra. Using an appropriate equivalent circuit to extract the resistance and capacitance of the electroactive grain and grain-boundary components, we apply a bricklayer method to estimate the corresponding conductivity and permittivity of the components. These values are then compared with the input data of the simulations to highlight the appropriate and limiting conditions of this method for data analysis.

II. Three-Dimensional Finite Element Approach
When an alternating voltage is applied across a material it generates a time-varying electric field, causing the propagation of charge carriers such as electrons, holes, or ions that generate a current through the material. The temporal evolution of an electrical field can be established by solving Maxwell's equations in time and space. Here, we highlight how this can be achieved within a finite element framework. We assume the material properties to be isotropic, linear, and time independent allowing the problem to be simplified. We also assume that, in the mHz-to-MHz frequency range, inductive effects are negligible when compared with the capacitive behavior. This allows the relationship of the timedependent electric displacement to the electric field to be simplified to where E ! ðr; tÞ is the local electric field and D ! ðr; tÞ is the electric displacement at time t and position r. As no time dispersion is considered, the electric permittivity e(r) is a function only of position.
Using the differential form of Maxwell's continuity equation: where j ! is the current density and q is the charge density. The current density can be written as: where r is the conductivity, j ! c is the differential form of Ohms law, and j ! d the displacement current density. Assuming the isotropic case, such that D ! ðr; tÞ ¼ eðrÞ E ! ðr; tÞ, the current density can be written such that impedance, Z, can be written as: with a real conductivity, r(t), and imaginary susceptibility, ixe(r), where i = √À1 and x is the angular frequency. Using the differential form of Gauss's law that: where u is the electric potential and e o is the permittivity of free space. Combining Eqs. (4) and (5), we transform Maxwell's Eq (2) into r j ! ðr; tÞ ¼ Àr ðrðrÞruðr; tÞ þ eðrÞ o ot r/ðr; tÞÞ ¼ 0 Using a time-domain finite element method (TDFEM) allows us to approximate the electric potential, u(r,t), in Eq. (6) as a function of space and time. This permits the current density to be calculated by integrating over the whole sample, and thus in turn allows simulation of the IS response of an electroceramic. Implementing Dirichlet boundary conditions at the electrode-air interface fixes the electric potential. We assume that the displacement currents crossing the free surface of the material are zero by using Neumann boundary conditions. A powerful aspect of this approach is that it allows the complete microstructure of the electrochemical system including contacts, grain boundaries, and grain cores to be created, meshed, and analyzed for their influence on IS data. Each grain and grain-boundary phase can be assigned its own unique time constant (i.e., intrinsic conductivity and permittivity) and we can therefore create heterogeneity within the ceramic microstructure. This model can then be calculated using the FEM package, without the requirement of an equivalent circuit consisting of a combination of R, C, and CPE elements. It should be noted that while the electrical response of the system can be solved using this package, it cannot account for any change due to a chemical or charge transfer process.

III. Results and Discussion
(1) Model Setup For our model design we base our technique on previous granular structure generation for magnetic simulations. 19 We first distribute an array of seed points, representing the centers of the grain cores within a cube. A Voronoi tessellation is then performed to generate a three-dimensional structure, filling the volume of the box. The arrangement of the seed points defines the structure of the system so that, if the points are distributed upon a regular grid and tessellated, a regular arrangement of cubic bricks is generated. Postprocessing on the Voronoi tessellation is then performed, eliminating any extremely small surfaces, to allow the structure to be discretized with tetrahedron elements with no complications. If required, the volumes can be shrunk toward their center point from thin gaps between individual volumes. These are then filled with prism elements that can represent very thin volume regions and can be assigned their own distinct material properties. This method allows for much higher volume ratios of differing thicknesses to be calculated. This then allows a volume ratio of the grain-core domain to the grain-boundary domain to be assigned, which we denote here as V gc:gb .
Applying material parameters of conductivity and permittivity to these regions, we solve using our FEM package to simulate the frequency response of a defined sample. Typically, we consider a frequency range 0.01 Hz-1 GHz using a potential of 100 V applied on a contact material with conductivity of 10 k/Sm.

(2) Response of Idealized Microstructures
To verify the capabilities of the FEM package we designed, simulated, and compared two simple structures; a simple layered structure shown in Fig. 1(a), and an encased structure depicted in Fig. 1(b). Each model is based on a cube with lateral dimensions of 1 lm and meshed using a combination of 250 000 tetrahedron and prism elements.
We first consider a simple layered system. A cube is divided into two distinct layers where the conductivity, r, is selected so that each region represents a grain core (gc) or grain boundary (gb) with r gc = 100 l/Sm and r gb = 0.1 l/Sm, respectively. The permittivity is held constant with a relative permittivity of e r = 100 for both the core and boundary regions. The thicknesses of the gc and gb are chosen to be identical, forming equal volumes of each material and V gc:gb = 1 as shown in Fig. 1(a). Here, the boundary makes up 50% of the total height of the cube.
A general analytical formulation can be written to describe this system using a BLM model with two RC elements in series, such that the total resistance R and capacitance C of the system are given by the combined total thickness, l, and cross-sectional area, A, of the different phases: Using the intrinsic dimensions and material properties, the analytical values of resistance and capacitance for the graincore and -boundary phases are calculated to be R gc = 5 GΩ and C gc = 1.77 fF, and R gb = 5 TΩ and C gb = 1.77 fF, respectively.
Using the FEM package we can solve this model to generate IS data that can be plotted and analyzed in all four immittance formalisms allowing values of resistance and capacitance to be extracted. The simulated IS data are shown here in the form of a Z* plot, Fig. 2(a), with the associated M″ spectroscopic plot in Fig. 2(b). To extract the values of resistance and capacitance for both the grain core and boundary, we apply an equivalent circuit fit of two parallel RC elements connected in series. The extracted values using this method give R gc = 5 GΩ and C gc = 1.75 fF and R gb = 4.99 GΩ and C gb = 1.76 fF for the grain-core and boundary phases, respectively. This shows excellent agreement with the values predicted analytically. Using the FEM IS data, the Z* arc for the layered structure also has no measureable depression angle and therefore exhibits a near-ideal Debye-like response with a measured FWHM of 1.15 decades in the corresponding M″ spectroscopic plot. It should also be noted that the two peaks in the M″ spectroscopic plot are of equal height indicating equal volume fractions of the two phases as expected from the BLM where the permittivity of the grain-core and -boundary phases is assumed to be the same.
We now compare this to an encased structure, where the grain-boundary material surrounds the grain core as shown in Fig. 1(b). We maintain the same volume ratio of V gc: gb = 1. To achieve this requires a grain-core cube of length of 0.794 lm to be encased by a boundary layer with a thickness of 0.103 lm. This gives a total lateral thickness of the boundary as 0.206 lm or 20.6% of the total thickness of the system. From the IS data in Figs. 2(a) and (b) a number of differences are observed. Although the volume ratio of the two materials is maintained, the grain-boundary response in the Z* plot is reduced from 5 to 2.25 TΩ. The M″ spectrum shows a drop in the peak height of the grain-boundary response and an increase in the peak height of the grain-core response indicating a volume ratio more in the order of V gc:gb = 5 than the true fraction. There is also an associated increase in the FWHM of the M″ Debye peak for the  grain-core response to 1.23 decades signifying an apparent electrical homogeneity in the sample. We extend this study to increased volume ratios where for simplicity we focus on two other ratios, that of V gc:gb = 10 and V gc:gb = 100 with dimensions shown in Table I.
As the volume fraction is increased to V gc:gb = 10 the Z* arc associated with the grain-boundary response remains less than half that of the layered structure, Fig. 2(c). The M″ Debye peak associated with the grain-core response, Fig. 2(d), however, decreases in FWHM to 1.18 decades indicating a more homogenous grain-core electrical response than for V gc:gb = 1, while the M″ peak height becomes comparable with that obtained from the layered structure. The trend continues for V gc:gb = 100 where although the Z* arc associated with the grain-boundary response remains lower for the encased model, Fig. 2(e), the grain-core M″ spectroscopic response shown in Fig. 2(f) is comparable for both models and is that of a near-ideal Debye-like response.
To explain this apparent electrical heterogeneity at low volume ratios we use the power of the FEM model to plot current density through the various models, as shown in Fig. 3. Figures 3(a)-(c) show the current density through the layered structures for V gc:gb = 1, 10, and 100, respectively. The current density is homogenous through each layer, indicating a linear flow of current through each phase, and as such resulting in the Debye-like responses in the simulated IS data. This does not occur for the encased structure. Figure 3(d) highlights that the current does not have to pass through the highly resistive grain boundaries whose normals are orthogonal to the current flow to reach the lower contact. These areas can be avoided, resulting in a near-zero current density within these regions. This leads not only to an increase in the current density through the grain-core regions (by over a factor of 3) but also in electrical heterogeneous behavior of the current flow. The IS data of the grain core thus results in an apparent electrically heterogeneous response with a FWHM in the M″ spectrum that exceeds 1.14 decades, even though the grain-core properties are homogeneous. This is also true at larger volume ratios, Figs. 3(e) and (f). These effects, however, are reduced as the grain boundary contributes less to the system as whole and thus smaller departures of the grain-core response from an ideal Debye-like response are observed.
To observe the effect of this change in microstructure on the extraction of the material properties we use the simulated IS data and the known lateral dimensions of the grain-core and -boundary regions to allow for the effects of sample geometry to obtain the intrinsic material properties. The parameters extracted from the layered structures for conductivity and relative permittivity (shown in Fig. 4 for a range of volume ratios) are all within 1% of the input parameters for both grain core and boundary. When the encased structure is analyzed in a similar way there is a strong dependency on the volume ratio. The resistance and capacitance of the core and boundary regions at large volume ratios agree with the true value; however, as the volume of the boundary region is increased, not only is the grain-boundary resistance underestimated by up to~50% but the grain-core resistance is also underestimated by up to~40%, Figs. 4(a) and (b). The permittivity for both the grain core and boundary is underestimated by 10%, Figs. 4(c) and (d). These results highlight how a small change in microstructure can significantly influence IS data and give rise to potential issues with extracting material properties using the BLM.

(3) Response of Complex Microstructures
Although complex phenomena are observed in these simple structures, the idealized brick-shaped grains do not represent the true complexity of ceramic microstructures. In the BLM, the connectivity between neighboring brick-shaped grains is exactly 6 facing and 12 edge paths between adjacent grains. Realistic grains do not obey this rule precisely due to their complex structural geometry and so the number of grain boundaries that are transversed changes, which therefore affects the impedance spectra.
We now expand our analysis from a simple two-component system, to systems that incorporate 512 individually addressable grains and associated grain boundaries arranged in an 8 9 8 9 8 configuration. To highlight the significance of microstructure we considered three designs. The first sys-  tem uses a simple layered structure of two materials to form a specific volume ratio. The second is based on the encased structure, featuring a regular brick-mortar system as shown in Fig. 5(a). The second builds upon the complexity of this using irregular shaped grains created by a Voronoi tessellation of random seed points. We maintain the average volume of each region at 1 lm 3 , however introduce a distribution with a standard deviation of 0.5 lm 3 . The typical volume spread in the Voronoi model is shown in Fig. 5(c). These individual volumes are then shrunk down around their individual centers to form the desired volume ratio, where the remaining volume is discretized to account for the grain boundaries. Due to the relatively low number of grains in each model compared with experimental samples, the results are averaged over 10 simulations. Simulated IS data of the layered, brick-mortar, and Voronoi models for volume ratios of V gc:gb = 1, 10, and 100, respectively, are shown in Fig. 6. As with the simple study, the layered structure shows a much larger grain-boundary response in the Z* plot due to the current path being forced through the whole layer. A reduction in the magnitude of the grain-boundary Z* arc is again observed when the grainboundary material surrounds the grain core. In all cases shown here, the reduction is of the order of a third, corresponding to the fraction of grain-boundary regions that are avoided. The Z* response of the grain boundary for the complex microstructure of the Voronoi model also decreases by a further 5%, whereas the core remains relatively unaffected. Figure 6(b) highlights the effect on the M″ spectrum, where the peak height of the grain-core response increases with the complexity of the microstructure and therefore the current pathways. As the volume ratio increases from V gc:gb = 1 to 100, the M″ peak associated with the grain core for the brick-mortar and Voronoi models begins to converge to that of the layered structure, indicating, as before, a more ideal Debye-like response. As shown in Fig. 7, at high volume ratios the grain-core response in the brick-mortar system shows a near-ideal Debye response where the FWHM is 1.15 decades. This is in contrast to the Voronoi model structure which still indicates significant nonideality with a FWHM of 1.20 decades. As the grain-boundary region becomes larger, the M″ Debye peak FWHM increases to 1.25 and 1.32 for the brick-mortar and Voronoi models, respectively. Current density plots at low frequency are shown in Figs. 8(a), (b), and (c) for V gc:gb = 1, 10, and 100, respectively. It is clear from these that a combination of irregular shaped grains with higher resistance grain boundaries influences the associated current flow around these areas, resulting in a large nonlinear response of the current density in the grain cores. At each volume ratio, an inserted image of a typical grain exhibiting significant nonlinearity of the current is shown. At equal fractions of grain core and boundary, the current density through a single electrically homogenous grain can be as large as 90% of the total. A similar but reduced effect (60% and 45% of the total current) is also observed at larger volume fractions (V gc:gb = 10 and 100, respectively).
We can use this FEM to predict the effect of microstructure on IS data and therefore comment on the confidence of using the BLM for such systems. We follow the standard experimental procedure to analyze the results of the model. First, we create a cross section of the model, similar to the images in Fig. 8, and use a line-scan method to estimate the percentage thickness of grain boundary and grain-core thickness for the system. Various slices through the model were used, finding that on average these values agreed well with  Table I for the various volume fractions. Using these values we correct for the geometry and extract the conductivity and permittivity of the grain-core and -boundary components. As shown in Fig. 9, the microstructure affects both the grain-core and -boundary response. At high grain-boundary volumes (boundary accounts for 20% of the thickness) a deviation of over 60% from the expected (input) values is obtained for the conductivity and over 10% for the relative permittivity. As before, this converges to the expected grain-core material properties as the volume fraction of grain core is increased. At values of V gc: gb > 20, the extracted values are all within 10% of the expected values. Thus, the use of the BLM method is shown to predict the correct values to within 10% for homogenous granular structures, when the grain boundary contributes 1% of the total thickness of the sample. Care should be taken with the grain-boundary material properties; however, as even at this ratio, the extracted properties are overesti-  show the grain-core and -boundary permittivity, respectively. A line for the eye is overlaid for each data set. mated by~12% for the conductivity and~20% for the permittivity.

IV. Conclusions
A fast and efficient FEM framework has been developed and used to allow a comprehensive study of IS data for threedimensional heterogeneous ceramics. In the model presented here we incorporate contacts, grain boundaries, and grain cores to replicate the microstructure of realistic ceramics; however, the flexibility of this code allows us to simulate virtually any microstructure such as a porous, nano-sized, and multiple-phase electroceramics. We show that an electrically homogenous grain core can give rise to an apparently heterogeneous IS response due highly resistive grain-boundary regions. Using the BLM (based on input materials parameters where the grain-boundary resistivity is four orders of magnitude larger than the grain-core resistivity but the permittivity of the two phases is the same) with its associated equivalent circuit to extract material properties (conductivity and permittivity) from IS data can lead to potential discrepancies of up to 60% of their true values based only on changes in microstructure.