Modeling of electrical transport in Zinc Oxide varistors

The electrical transport in zinc oxide (ZnO) varistors is analyzed using microstructural material modeling. The fully three dimensional current distribution is computed by means of a nonlinear equivalent circuit model representing the assembly of current carrying grains and grain boundaries of the material. The investigation focuses on the phenomenon of current filamentation due to inhomogeneities of the varistor microstructure. Numerical results highlight the importance of 3-D percolation effects in the modeling of varistor currents as well as that of the grain bulk resistivity which so far has been neglected in previous studies.


Introduction
ZnO varistors are ceramic semiconductor materials that are characterized by a highly nonlinear electrical response.Their resistance drops remarkably with increasing voltage.Thus, a varistor based device behaving as an insulator under normal operating conditions, will draw large currents from the circuit as soon as an overvoltage event occurs (Meshkatoddini, 2011).ZnO varistors are, furthermore, capable of absorbing high amounts of energy (Gupta, 1990).For this reasons, they are most often applied as surge arresters for protection against temporary electrical overstresses, such as lightning and switching surges in power distribution lines, transformers and switches (Standler, 2002).Further applications of ZnO varistors include passive voltage stabilization for low power electronics (Eda, 1989).
The properties of ZnO varistors are a result of their crystalline microstructure.Varistor ceramics are composed of a large number of ZnO grains of different shape, size and orientation.Small amounts of doping additives such as Bi 2 O 3 or Mn 3 O 4 are usually present (Eda, 1984;Gupta, 1990).The nonlinear electrical response is attributed to the potential bar-riers forming at the boundaries between ZnO grains (Blatter and Greuter, 1986).These boundaries represent, in fact, Schottky-like microjunctions with an electrical response corresponding roughly to that of a back-to-back pair of Zener diodes.The grain boundary resistances determine varistor currents for voltages in the nonlinear breakdown voltage region of the electrical characteristics as well as in the linear pre-breakdown region.The resistance of the varistor at higher voltages becomes again Ohmic as it is essentially determined by grain bulk resistivity (Vojta et al., 1996).
Although this basic operation principle has been known for some time, the characterization of ZnO varistors for engineering applications remains a major challenge.Several numerical models have been proposed to describe electrical transport in varistors (Bartkowiak and Mahan, 1995;Vojta et al., 1996;Qingheng et al., 2002;Hofstaetter and Supancic, 2013).In these models, the varistor microstructure is first represented geometrically by a space filling mesh.The mesh cells correspond to the ZnO grains whereas cell interfaces (mesh edges in 2-D or mesh faces in 3-D) to the grain boundaries.Almost all of the models so far presented in the literature operate in 2-D planar geometry.They are based either on a Cartesian mesh representation (Vojta et al., 1996;Wen and Clarke, 1993) or on a Voronoi-type tessellation of the varistor sample (Bartkowiak et al., 1996a;Zhao et al., 2007).Given a microstructure representation, electrical transport can be modeled by means of an equivalent electrical circuit with voltage nodes and current branches associated with the grain centers and grain boundaries, respectively.The grain boundary resistances are taken into account by introducing diode-like elements with specified electrical characteristics (Vojta et al., 1996) in the current branches.The macroscopic electrical response is, then, obtained by the solution of the resulting circuit equations.It should be noted that varistor microstructure (grain sizes and shapes) as well as grain boundary resistances are subject to stochastic distributions which makes the modeling task in the general case very difficult.
In this work, we focus on two particular aspects of varistor modeling.The first concerns the influence of the 3-D geometry on modeling accuracy.Since percolation thresholds for 2-D and 3-D systems are in general very different, the commonly used 2-D models fail to reproduce properly the current distribution in the material and thus compute the varistor characteristic accurately.We perform a comparison between 2-D and 3-D varistor models using the approach of (Vojta et al., 1996) based on a Cartesian mesh representation.The second aspect is related to the effect of grain bulk resistivity.This effect involving both grain geometry and the topology of the grain network has been thoroughly neglected in previous studies.Finally, we present a simulation example featuring the current filamentation phenomenon using a 3-D modeling approach for a realistic varistor sample.

Varistor modeling
ZnO varistors are basically composed of polyhedral grains of various shapes and sizes, typically, varying in the range of 10-100 µm.The electrical current distribution in the material is strongly affected by its microstructure including grain shape, topology of the grain network as well as type and distribution of crystalline defects.The first modeling task consists, thus, in the construction of an appropriate geometrical model of the varistor microstructure.

Geometrical representation
The simplest geometrical model consists of a 2-D Cartesian mesh of square grains.This approach was used, e.g., in (Wen and Clarke, 1993;Vojta et al., 1996).Because of the simple geometry, it is not possible to take into account local grain size variations.Also, the resulting topology of the grain network is far from realistic.This model, however, allows for investigating current flow percolation effects in a system- For a more general geometrical description, an unstructured mesh representation is needed.Unlike the commonly used 2-D tessellations (Bartkowiak et al., 1996a;Zhao et al., 2007;Hu et al., 2010), however, we aim at a fully 3-D polyhedral grain model.The following 3-D tessellation procedure is proposed.We start with a simplicial tetrahedral mesh on the varistor sample.Such a mesh can be easily constructed using, e.g., a conventional FEM tool.The only constrain imposed on this 'primary' mesh is that the average edge length should coincide with the sought grain size of the material.
A polyhedral grain structure is constructed on the basis of the primary mesh by applying a midpoint subdivision procedure (Li et al., 1995).According to this procedure, new mesh nodes are introduced at the midpoints of every edge, face and cell of the primary mesh, respectively.Mesh cells of polyhedral shape (representing the ZnO grains) are, then, constructed by connecting these nodes as exemplarily shown in Fig. 1a.By construction, each node of the primary mesh corresponds to one cell of the polyhedral mesh.These nodes can be considered as the generating germs of the resulting grain structure.The grain boundaries are represented by the generally non-planar interfaces between the neighboring polyhedral cells.For illustration, the midpoint subdivision of a tetrahedron into 4 hexahedrons is depicted in Fig. 1b.
Note that, in general, grain cell convexity is not guaranteed.Also, grains with sharp angle corners and very small edges may occur.It is, however, not the purpose of this construction to generate a high mesh quality in the FEM sense.This procedure allows in the first place for a more realistic representation of the varistor microstructure than that of a simple Cartesian mesh model.Furthermore, grain size, grain number and (to some extent) grain shapes can be easily controlled.Tessellations for arbitrary material samples of, e.g., cubic (cf.Fig. 2) or cylindrical shape can be generated.Due to the robustness of modern mesh generators it is, furthermore, possible to construct large varistor models with several hundreds of thousands of grains.1).(b) Grain boundary equivalent circuit corresponding to the model of (Vojta et al., 1996).
Another strength of this approach is the ability to model microstructure inhomogeneities.This is important for the investigation of crystalline defects or inclusions and their impact on the varistor characteristic.Figure 3 demonstrates the model of a varistor with a non-uniform grain size distribution generated by this approach.The defective pattern is initially prescribed by means of a few small spheres properly located within the sample volume.The result is an nonuniform primary mesh and, therefore, a non-uniform polygonal grain tessellation with grain size distribution closely following the assumed defect distribution.

Grain boundary conductivity models
The J − V characteristic of a single grain boundary in the breakdown region can be described by the empirical relation where k is a parameter, α is the dynamic conductance, J is the current density and V the voltage drop across the boundary.In terms of an equivalent electrical circuit the contribution of a single grain boundary can be incorporated by introducing a nonlinear resistor between adjacent grains sharing this boundary (see Fig. 4a).This approach, however, does not count for the grain bulk resistivity.Furthermore, the reverse leakage diode currents which are responsible for the pre-breakdown varistor response are neglected.
A more general approach was proposed by (Vojta et al., 1996) and then later used by several authors (Zhao et al., 2007).The grain boundary characteristic is assumed to where V B is the barrier voltage, σ s and σ g are Ohmic conductivities, d is the distance between the centers of the two grains sharing the boundary and s is the nonlinearity varistor parameter.A closer inspection of Eq. ( 2) reveals that this is equivalent to two resistors connected in parallel across the boundary (see Fig. 4b).This configuration is able to reproduce qualitatively the varistor response in the full voltage range.Note that σ g is interpreted as the conductivity of the grain bulk material since it determines the characteristics at high voltages.
Referring to Fig. 4b, however, it is not possible to separate the contributions of grain bulk and grain boundary at lower voltages.Thus, neither the shape nor the size of the single ZnO grains can be explicitly accounted for in this approach.
The physical understanding that grain currents enter through its boundaries only, implies that bulk and boundary resistors should be connected in series.

Distributed grain conductance model
A more accurate model is obtained by employing a conductance matrix to describe current flow within the ZnO grains.
In order to illustrate the idea we consider a single hexagonal grain as shown in Fig. 5a.Neglecting leakage resistances, the model of (Vojta et al., 1996) reduces to a star network of resistors.Applying an arbitrary voltage V 1 at one of the boundaries while keeping the rest at a fixed potential V 0 , results in an equal current distribution among the remaining boundaries.This is true independently from the number of boundaries or the actual grain shape.As can be easily verified by simulations, however, most of the current should flow through the boundaries closest to the excited one.For the varistor problem, where the macroscopic characteristic depends strongly on the current flow pattern within the material, this approach may lead to inaccurate results.We propose the elementary equivalent circuit shown in Fig. 5b.The grain bulk resistance is described by the conductance matrix G of a N-port network, where N is the number of grain boundaries.The grain boundary resistances are connected in series at each port of the network.They consist of a leakage resistor connected in parallel to a nonlinear one obeying the J − V characteristic (Eq.1).The conductance matrix can be determined numerically by 3-D FEM simulations for arbitrarily shaped grains.The elementary equivalent circuits constructed by this principle are, then, assembled into a global equivalent circuit which can be solved numerically for the current distribution within the varistor.The drawback of this approach is that a large number of FEM simulations may be needed to obtain conductance matrices for each individual grain within the material.

Influence of 3-D geometry
We consider the simplest varistor model using a Cartesian grain structure with microjunction characteristics given by Eq. ( 2).The 3-D model consists of a 15 × 15 × 15 mesh of cubic grains with a side length of 10 µm.The 2-D geometry case is emulated by employing a mesh with 15×15×1 grains of the same size.The situation is illustrated in Fig. 6a.For consistency, the global varistor characteristic computed with the two models for a uniform distribution of microjunction properties is shown in Fig. 6b.In this case, all grain boundaries are described by the same parameters, σ s = 10 −3 S / m, σ g = 1 S / m, s = 8 and V B = 3.2 V.The current flow between top and bottom electrodes is uniform, thus, 2-D and 3-D models yield the same varistor characteristic.
The influence of 3-D geometry becomes obvious when inhomogeneity is introduced into the model.We adopt the approach of (Vojta et al., 1996) and (Bartkowiak et al., 1996a) where a number of grain boundaries uniformly distributed within the sample volume are assumed to be electrically inactive.This is implemented by setting the diode parameters for these boundaries to s = 0 and σ s = 0 S / m. Figure 7 shows the results obtained when this procedure is applied for a fixed portion of inactive grain boundaries (20 %).A sequence of models is considered which is obtained by adding grain layers in the z direction to the 2-D model (single layer in the z direction) while keeping the grain size constant.The 3-D varistor characteristic is shifted toward lower voltages as the number of layers is increased.The lower switching voltage in the 3-D case compared to the 2-D one can be explained by the higher probability for 3-D currents to find a conductive path to the bottom electrode.By further increasing the number of layers (see Fig. 7c) this effect saturates and the varistor characteristics does not change any more.At this point, the percolation threshold for the current paths in 3-D geometry has been reached.In general, this phenomenon depends on grain size distribution and on the topology of the grain network.This simple example, however, already demonstrates the necessity of 3-D modelling in varistor simulations.

Influence of grain bulk conductivity
The same experiment is repeated for the distributed grain conductance model described in Sect.2.3.The nominal grain boundary parameters are assumed to k = 4.055 × 10 −18 S / m 2 , α = 45 and σ s = 10 −3 S / m.Again, the same shifting of the switching voltage towards lower values is observed, when 20 % of the boundaries are assumed to be inactive.When a sufficient number of 3-D layers has been added the varistor characteristic converges to its final shape (see Fig. 8).Compared to the results obtained by the Vojta approach (Fig. 7), this characteristic appears much smoother in the transition region between breakdown and upturn voltages.This behavior is due to the influence of grain bulk  conductivity which is properly taken into account in the distributed conductance model.

Current filamentation and varistor characteristics
In the following, a more realistic varistor simulation is discussed.A cylindrical sample of radius 50 µm and height 200 µm with an average grain size of 10 µm is considered.The microstructure geometry is constructed by applying the tessellation procedure described in Sect.2.1.The resulting 3-D grain model consists of a total of 1737 ZnO grains (see the inset of Fig. 9).Furthermore, the model contains 10 081 grain boundaries corresponding to the number of nonlinear resistors employed in the equivalent electrical circuit model.
To model the grain boundary microjunctions we apply the Vojta approach (Eq.2) with nominal parameters, σ g = 1 S / m, σ s = 10 − 6 S / m, V B = 3.2 V and s = 8, respectively.Inhomogeneity is introduced into the model by assuming a certain number of inactive microjuctions randomly distributed among the available grain boundaries.Figure 10 shows the current distribution within a varistor containing 10 % inactive boundaries for three different voltages applied at the top and bottom electrodes.Depicted is the current flowing through each grain normalized to the total varistor current.For visualization purposes, only grains carrying a current above a certain threshold are shown.An inhomogeneous current distribution is observed already at low voltages (see Fig. 10a).In this case, a single grain can carry up to 50 % of the total varistor current.In the breakdown region, this phenomenon becomes substantially stronger.As shown in Fig. 10b, nearly the total varistor current flows along a single narrow path of width no larger than twice the grain size.As the voltage is increased above the upturn voltage, the current distribution becomes again more homogeneous (see Fig. 10c).The varistor response in this region is essentially Ohmic, thus, as expected the effect of inactive boundaries in the current flow distribution is less pronounced.
The influence of inhomogeneous current flow on the overall varistor performance is illustrated in Fig. 9. Three cases with 5, 10 and 20 % inactive boundaries, respectively, are considered.The J − V characteristic for 5 % inactive boundaries corresponds essentially to the single grain boundary response (Eq.2).Increasing the number of inactive boundaries causes the switching woltage to shift towards lower values while the effective dynamic conductance of the material is decreased.The transition from fully nonlinear to Ohmic varistor response is, however, not trivial.As already observed in (Bartkowiak et al., 1996b) using 2-D simulations, the J −V characteristic of a varistor containing inactive grain boundaries manifests a double-knee pattern featuring a second inflection in the current-voltage curve in the region between switching and upturn voltages, respectively.

Conclusions
The simulation of electrical transport in ZnO varistors requires 3-D modeling.This was demonstrated by a direct comparison of the electrical varistor characteristics computed by 2-D and 3-D simulations.The simulations show a systematic shift of the switching voltage toward lower values in the 3-D case.This can be explained by the intrinsically inhomogeneous current flow pattern through the grain www.adv-radio-sci.net/12/29/2014/Adv.Radio Sci., 12, 29-34, 2014 microstructure which implies that the percolation behavior of the current paths in 3-D space is different from the 2-D one.Furthermore, we propose a novel equivalent circuit approach employing a distributed grain conductance model.This approach allows to properly take into account the influence of grain bulk resistivity on the macroscopic varistor response for arbitrarily shaped ZnO grains.A procedure for the fully 3-D geometrical representation of varistor microstructures based on a particular 3-D tessellation method is introduced.Using this approach, simulations for a realistic varistor sample were performed.They show the experimentally observed current filamentation phenomenon resulting from microstructural inhomogeneities which strongly influence the electrical performance of ZnO varistors.

Figure 1 .
Figure 1.(a) Basic procedure for the construction of a polygonal grain tessellation on the basis of a simplicial mesh.(b) Example barycentric decomposition of a tetrahedron.

Figure 3 .
Figure 3. Modeling of a nonuniform grain size distribution: (a) Geometrical description of the defective pattern by small spheres located within the volume of a cubic varistor sample.3-D cut (b) and cross-sectional view (c) of the resulting grain structure.

Figure 4 .
Figure 4. (a) Equivalent electrical circuit representation of a grain boundary according to Eq. (1).(b) Grain boundary equivalent circuit corresponding to the model of(Vojta et al., 1996).

Figure 5 .
Figure 5. (a) Grain boundary resistors for a hexagonal grain and current distribution resulting from the model of(Vojta et al., 1996).(b) Equivalent circuit representation in the distributed grain conductance approach.

Figure 6 .
Figure 6.(a) Cartesian mesh grain structure and corresponding equivalent circuit.(b) Varistor characteristics resulting form 2-D and 3-D simulations, respectively, when a uniform distribution of the microjunction properties is assumed.

Figure 7 .
Figure 7. Varistor characteristics obtained when 20 % of the microjuctions are assumed to be inactive.A sequence of 3-D models is considered which is obtained by adding grain layers in the z direction: (a) 1 layer.(b) 3 layers.(c) 10 layers.(d) 15 layers.

Figure 8 .
Figure 8. 3-D Varistor characteristics resulting from the distributed grain conductance model when 20 % of the microjuctions are assumed to be inactive.The 3-D models consist of (a) 3 grain layers and (b) 10 grain layers, respectively.

Figure 9 .
Figure 9. Grain structure representation of a cylindrical varistor sample and the computed characteristics for 5, 10 and 20 % inactive boundaries, respectively.