Conductivity simulations of field-grading composites

The electrical conductivity and the percolation threshold of field grading polymer composites intended for high voltage applications were examined with representative elementary volume simulation methods based on percolation threshold modeling (PTM) and electrical network modeling (ENM). Comparisons were made with experimental conductivity data for SiC-EPDM composites with spherical and angular particles, using different filler fractions and electrical field strengths. With a known conductivity of the filler particles (powder), the simulations could predict the percolation threshold and the composite conductivity as functions of the electrical field for a wide range of SiC-filler fractions. The effects of morphology, dispersion and filler shape were examined and the simulations were able to explain the experimental difficulty of reaching sufficient reproducibility when designing composites with filler fractions close to a percolation threshold. PTM of composites containing hard-core/soft-shell spheres revealed a y  =  (a  +  bx)(−1/c) relationship (R2  =  0.9997) between filler fraction and relative soft-shell thickness.


Introduction
One of the foremost challenges facing society today is the efficient generation, transmission and distribution of energy. High-voltage DC (HVDC) cable systems are important for power transmission over long distances, especially across seas and in densely populated regions where overhead cables cannot be used. By using higher voltage levels, the energy losses during transmission are reduced, but the insulation system is then exposed to higher electrical and thermal stresses, increasing the likelihood of electrical breakdown, especially at critical cable interfaces and triple points. Electric field-grading materials (FGM) are used in order to reduce the local enhancement of the electric field in these critical regions, by having a field-dependent resistivity. The material should be conductive at high fields and insulating at low fields. FGM are percolated composites based on inorganic semi-conductive particles, e.g. silicon carbide and zinc oxide, or conducting fillers, such as carbon black, typically in polymers such as EPDM rubber and silicone rubber [1].
Percolation theory is a critical concept in understanding the conductivity of a composite. It is generally accepted that, above a critical volume fraction of semiconducting fillers in an insulating matrix, the electrical conductivity increases dramatically due to the development of a continuous electron transport path [2]. The critical volume fraction φ is often referred to as the percolation threshold (φ p ) [2]. In order to manufacture field grading materials with a robust non-linear conductivity, it is important to understand the percolation phenom enon and the resulting dielectrical properties.
Theoretical percolation thresholds for idealized composites can be determined with numerical Monte-Carlo techniques, in either two or three spatial dimensions. A common approach is to simulate a small but representative brick-shaped section of (Some figures may appear in colour only in the online journal) Original content from this work may be used under the terms of the Creative Commons Attribution 3.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. the composite and examine at which filler fraction a continuous path is formed from one side of the cubic domain to the opposite side. Percolation simulations can be performed either on-lattice or off-lattice and the particles in the simulated composites can be solid, soft (possibly overlapping), or soft with a solid core. For spherical soft particles, off-lattice simulations predict a percolation threshold of φ p = 0.2895 [3,4]. Percolation for other geometries has also been studied, including soft cylinders [5], soft ellipsoids [6], various lattice structures [7][8][9], hard-core/soft-shell sphero-cylinders [10] and platelets [11].
If the conductivities of both the matrix and the filler are known in advance, the composite conductivity can be predicted quantitatively with various kinds of effective media modeling, for instance with finite element modeling (FEM) [12][13][14][15], explicit analytical mathematical expressions [13,[16][17][18][19][20] or on-lattice electrical network model (ENM) simulations [21][22][23]. Due to computer-resource limitations, previous on-lattice ENM simulations of field-grading composites [21][22][23] were performed on too small lattices (maximum 24 × 24 × 24 grid-points) for enabling systematic studies of particle clustering and other super structures, something which is addressed in this study. Further, the previous electrical network models were always restricted to a regular mesh, but in this paper an off-lattice ENM model is also presented.
The morphology of the composite has a documented impact on both the composite conductivity and φ p . Highaspect-ratio fillers typically give lower φ p values than their spherical equivalents [20], fillers oriented along the electrical field give a lower φ p than perpendicular fillers [13], and spherical particles typically give rise to a single threshold while angular particles may give rise to a double percolation threshold [24]. The composite conductivity often decreases with improved dispersion [25], and with increasing particle size [26,27].
The composite conductivity generally increases with increasing electrical field strength E but can either increase (at low E) or decrease (at high E) with temperature T. This can be explained by the conduction mechanisms in the composite. Under moderate electric fields the dominant electronic charge carriers are (i) hopping via localized states in the polymer bulk and (ii) double Schottky-barrier transport at the particle-particle interfaces. Both these processes are thermally activated and the conductivity increases with increasing temper ature [28][29][30]. Under high electric fields, the applied field influences the barrier heights and hot-electrons limit electron movements, leading to a reduced conductivity with increasing temperature [31,32].
The conductivity of a polymer composite is further influenced by the constituent properties, the polymer crystal fraction, saline content, additives, chemical reactions, pressure, packing, time, electrical field history and moisture content [33][34][35].
The aims of the present study were (i) to develop improved simulation tools for predicting the percolation threshold and the electrical conductivity of field-grading composites and (ii) to use the models to analyse how improved field grading composites might be constructed.

Models and methods
The electrical conductivity and the percolation threshold of polymer composites with known polymer and particle conductivities were studied with two meso-scale simulation techniques based on PTM and ENM. Both techniques assumed that the composite was globally homogeneous and could thus be represented by a sufficiently large, approximately repeating, sub-volume of the material, i.e. a representative elementary volume [36]. Composite morphologies were constructed either with on-lattice strategies, where the particle positions were restricted by a 3D grid, or with off-lattice strategies, where the particle positions were free.

Geometry assembly
On-lattice geometries with spherical and cubic fillers were generated with a Monte-Carlo algorithm on hexagonal and face-centered cubic (FCC) grids, respectively. A typical grid contained 50 × 50 × 50 grid points, where filled positions represented filler particles and empty positions represented the polymer matrix. Each filled grid position corresponded to one particle and adjacent particles were always in close contact. During the on-lattice Monte-Carlo process, new particles were successively added to the grid stochastically until the desired filler fraction was reached. Super-structures were assessed by using grid-filling probabilities depending on the number of adjacent particles. Dispersed and agglomerated morphologies were obtained by setting a low or high probability for adding new particles close to previous fillers, while string-shaped morphologies were obtained by increasing the probability of adding particles at positions with only one or two adjacent particles.
Off-lattice geometries with spherical fillers were constructed with a Monte-Carlo algorithm without any underlying grid. One by one, spherical particles were added at random positions within a brick-shaped periodic domain until the desired filler fraction was reached. Suggested particle positions were accepted only if they did not cause illegal collisions with previous particles [37,38]. The particles could either be solid or soft. In the latter case, the suggested positions were always accepted, since overlapping particles were then tolerated. Clustering was simulated by moving new spheres into close contact with their nearest neighbor, assuming that the nearest particle was positioned within a predefined radius.

Percolation threshold models (PTM)
The essence of PTM is to examine when an unbroken conductive path is formed through the whole representative volume of a simulated composite. The percolation threshold can be defined as the filler fraction when either two specific or two arbitrary opposite faces first become connected, yielding higher or lower values of φ p , respectively. Further, φ p can either occur when the surfaces themselves become connected (lower φ p ) or when an infinite path, considering periodic boundary conditions, is formed between them (higher φ p ) [4]. We define φ p as a connection between two specific faces. PTM determines only φ p and thus not directly the composite conductivity. However, φ p can be used in semi-empirical conductivity models [20] and the PTM methodology can be applied to much larger systems than is possible with FEM and to much more complex particle shapes than is possible with ENM.
In this study, the PTM methodology was applied to offlattice geometries containing either solid or soft spherical particles, which were added to the domain one by one until percolation was reached. The solid particles were surrounded by soft shells, the thicknesses of which were systematically varied in order to determine the effect of layer thickness on percolation. A connection between two particles was created if they were in close contact, if they overlapped or if their shells overlapped. The connection pattern between the particles was stored in a sparse connectivity matrix which was continuously updated. Particle clusters with a connection either to the voltage source or to the ground were tracked. Percolation was reached when the first particle cluster was connected to both ground and voltage. The average φ p -value of 100 simulated composites was always used.

Electrical network models (ENM)
Electrical network models were developed to predict the conductivity of composites containing many particles with a simple geometrical shape, such as a sphere or a cube. The conductivity of each simulated composite was calculated by transforming the geometry into an equivalent electrical circuit and solving with respect to the laws of Ohm and Kirchhoff. Each node in the electrical network was either a grid position (on-lattice) or a particle (off lattice). The total number of nodes (n) was thus either equal to the total number of meshpoints or to the number of particles. Each node was connected to its nearest neighbors and the total number of connections was denoted (m). A binary m × n matrix (A) was used to describe the connections between the nodes in the electrical network. A diagonal m × m matrix (C) was used to store the conductance of each connection. The boundary conditions for voltages and currents were stored in the n × 1 vector (b) and in the m × 1 vector ( f ), respectively. The current ( y ) over each connection was calculated with Kirchhoff's current law A T y = f, while the n node potentials (x) were computed with Ohm's law y = C(b − Ax). Linear networks with C independent of x were assumed and the equations could thus be solved simultaneously [39]: The shape of the domain was always a rectangular box with sides L x , L y and L z and with six faces. The upper face was connected to a voltage source (V Γ = V 0 ) and the bottom face to ground (V Γ = V 0 ), and the remaining faces had repeating boundary connections (V a = V b ). The non-normalized composite conductivity was calculated as the sum of all currents exiting the grounded face. Normalization was achieved by dividing by V 0 × L x × L y and multiplying by L z . The mean (or median) conductivity of 100 generated geometries was always reported.
For on-lattice geometries, the topology matrix A exhibited a band structure where the number a bands equalled the number of connections per node. A hexagonal grid has 12 bands, a FCC grid with only face contacts has 6, a FFC including edge contacts has 18 and a FCC including corner contacts has 26. Corner contacts can however be neglected [24]. Face contacts were modelled as being coupled in series σ = 2σ p σ m /(σ p + σ m ), and edge contacts as geometrical sums σ = sqrt(σ p σ m ) [24], where σ m and σ p are the matrix and particle conductivities. The composite conductivity was normalized with respect of the number of downward connections from each node, i.e. 3 for the hexagonal lattice and 5 for the FCC without corner contacts. Cubic particles were represented by the FCC lattice and spherical particles by the hexagonal lattice.
For off-lattice geometries with spherical fillers, the topology matrix was constructed by first finding all neighbors within a predefined radius of each sphere. If more than 10 neighbors remained, only the 10 closest were considered. The resistance of each link was calculated by summing the contributions from the matrix, the particle and the interlayer, taking into account the thickness of each phase. The conductivity of the network was calculated with equation (1) and normalized, while the conductivity of the remaining parts of the domain was approximated with Maxwell-Garnetts mixing rule, σ = [18]. The total composite conductivity was finally calculated as the sum of the two conductivity contributions multiplied by their respective volume fractions.

On-lattice results
Electrical network (ENM) simulations of composites containing mono-disperse spherical particles were carried out for composites generated with the on-lattice Monte-Carlo model, using hexagonal meshes with 50 × 50 × 50 grid-points. The filler fraction φ ranged from 0 to 50 vol.% and the conductivity of the particles was set to 10 9 times the matrix conductivity. Composites with well-dispersed, clustered and string-shaped morphologies were constructed and are compared in figure 1, where typical geometries with 10 vol.% particles are also inserted.
The composites with well-dispersed particle distributions showed higher φ p values than the materials forming clusters and strings. A similar dispersion effect was observed in experimental studies [25]. Agglomerates and especially stringshaped formations thus increase the probability of forming percolating paths through the composites. This is partly an effect of the higher aspect ratios of the super-structures than of the individual spherical particles. Well dispersed spheres have the highest φ p of all geometrical shapes [3,4]. The simulations stress the possibility of constructing field-grading composites with low filler fractions by using self-assembling particles that form string-shaped morphologies [40]. Suitable particles could probably be constructed by grafting techniques [41].
Well dispersed composites, with both spherical and cubic particles, were examined in detail with the on-lattice ENM model. Smaller meshes (30 × 30 × 30 grid points) were here used, since the mean and median values of σ still approximately coincided, indicating sufficiently large meshes. The particle conductivities were 10 3 , 10 6 and 10 9 times larger than the matrix conductivity. Similar simulations were also carried out for composites resembling EPDM rubber (σ m = 10 −15 S m −1 ) containing well dispersed powder of either spherical green (σ f = 3 · 10 −10 ) or angular black (σ f = 10 −7 ) silicon carbide (SiC) [42]. The results were compared with experimental data [24], and all the results are shown in figure 2. Four observations were made: (i) The standard deviations were always small far from the percolation threshold, but close to the threshold they became significant also for well dispersed high quality samples, and this should be taken into account when developing, testing and manufacturing real composites. (ii) The maximum standard deviation increased notably with increasing filler conductivity and thus indirectly with the electric field strength. (iii) Spherical fillers resulted in a single φ p value. Highly-conducting cubic fillers resulted in double thresholds, the lower φ p being a consequence of edge contacts while the higher φ p was a result of face contacts [24]. A meta-stable plateau, which can be investigated further for field grading purposes, was located between the two thresholds. (iv) The simulations gave φ p values approximately 10 vol.% lower than the experimental values, which was in agreement with previous findings [24], although the shapes of the curves were conceptually correct.
The main conclusion is that composites intended for highvoltage field-grading applications should have a value not below φ p , otherwise they will typically not be field grading, but nor should they have a value too close to φ p , otherwise the variations between different material sections will be too high and the overall composite quality will become unacceptably low.
The capability of the on-lattice ENM model to predict field-dependence was examined by varying the electrical field (E ) from 1 · 10 5 to 1 · 10 6 V m −1 for composites with black angular SiC in EPDM rubber. The conductivity of the EPDM was constant while the field-dependent response of the SiC powder was described by: The simulation curves were basically in good agreement with the experimental curves, although the latter were slightly more concave at low E-values. The intrinsic fieldgrading properties of the powder thus seem to dominate the field-grading behaviour of the composite, which leads to the conclusion that powder measurements combined with simulations can be used to predict the field dependence of composites. Powder measurements clearly contain vital information about both of the particle size effect and the Schottky-induced field-grading effect.

OFF-lattice results
PTM results for composites containing spherical inclusions, with or without inter-particle attractions within 10 sphere radii, were compared. The size of the domain was 50 sphere radii. The thickness of the conducting layer divided by the sphere radius was defined as d. In figure 4, the resulting φ p values are plotted against d for the two different attraction settings.
For thick conducting layers (d > 0.3) the φ p value was higher for composites with attractive forces, but for smaller d-values, the opposite trend was observed, since the particles must become closer in order to connect when d decreases. The experimentally measured φ p value was obtained at about d = 0.15 and d = 0.20 for composites with and without attractive forces, respectively. The standard deviations decreased with increasing attraction distances, indicating that the presence of agglomerates leads to a greater probability for undesired variations in manufactured composite materials.
The φ p curves were found to be accurately fitted by the function y = (a + bx) (−1/c) . For the morphologies without attractive forces, the values a = 1.695 ± 0.132, b = 2.032 ± 0.581, c = 0.463 ± 0.071 gave a fit of R 2 = 0.9997. By repeating the simulations with other interaction distances, we found that a, b and c could be approximated as convex third-order polynomial functions with respect of the inter-particle interaction distance, as seen in the inset in figure 4.
The effect of inter-particle attractions was also studied, for the same geometries, with off-lattice ENM. EPDM (σ m = 10 −15 S m −1 ), green spherical SiC (σ f = 3 · 10 −10 S m −1 ) and an interface with thickness d = 0.20 and conductivity σ f were used as composite constituents. The results are presented in figure 5. The particle contact conditions and the composite morphology have a significant impact on the composite conductivity. Clustered particle systems have a much lower φ p value than composites with dispersed particle distributions, and this has been confirmed experimentally [25].
The influence of the conductive layer thickness (d ) was also examined with off-lattice ENM, using green spherical SiC in EPDM rubber. The median of the normalized conductivity was calculated for filler fractions ranging between 0-30 vol.% for d = 0, 0.05, 0.1, 0.15, 0.20, 0.25 and 0.50. The results could also represent the effect of particle size, assuming a fixed d. The data are reported in figures 6. In figure 6(a) no inter-particle attractions are included, but in figure 6(b) clustering within 10 sphere radii has been applied.
The φ p value decreased with increasing d (or with decreasing particle size for a fixed d ), especially in the case of the samples without attractive forces, in agreement with experimental findings [26,27]. Since φ p is affected both by d and by the attraction distance, it would be difficult to distinguish between the effect of these two factors without doing simulations.
In figure 7, off-lattice ENM simulation results are presented for geometries containing randomly positioned soft-core particles, which can overlap each other. The constituents were spherical green SiC in EPDM. The simulations results are presented with corresponding experimental data [24]. The slope of the conductivity increase below the percolation threshold Figure 6. The impact of the soft-shell thickness (or particle size) of composites containing solid spherical fillers (a) with and (b) without inter-particle attractions, studied with off-lattice ENM simulations. The σ/σ m (median) is plotted against φ for seven shell thicknesses. in the simulations is higher than expected, but considering that the simulations are predictions without adjustable parameters, the agreement with experimental data is good. The processing of this kind of materials is most probably sensitive to the optimization process of the compounding, which also could explain deviation from the theoretical model.

Conclusions
Two representative elementary volume simulation methods based on PTM and ENM have been developed in order to predict the electrical conductivity and percolation threshold (φ p ) for field-grading composites. The effects of agglomeration, dispersion, filler fraction and particle shape have been systematically examined by creating various morphologies with Monte-Carlo techniques. PTM/ENM studies of dispersion, both with and without an underlying 3D-lattice, predicted that φ p would decrease with the presence of agglomerates, and this was in agreement with experimental findings. Off-lattice hard-core/soft-shell PTM revealed a (R 2 = 0.9997) dependence φ p = (a + bx) (−1/c) between φ p and the relative soft-shell thickness, where a, b and c were dependent on the particle interaction distance as convex third order polynomials. The dispersion simulations also predicted that low-filled field-grading composites with string-shaped morphologies could be constructed using partially attractive, partially repulsive, spherical fillers. On-lattice ENM simulations resulted in a single or double φ p for composites with spherical or cubic fillers, respectively, which was qualitatively in agreement with experimental measurements on SiC/EPDM composites. The field grading effect of SiC/EPDM was found to be dominated by the intrinsic properties of the SiC-powder. ENM-simulations also predicted comparatively large conductivity variations between materials with filler fractions close to φ p , even for well dispersed high-quality composites, an observation that may be of importance when designing highly filled composites for industrial applications. Off-lattice ENM soft-core simulations, without adjustable parameters, accurately predicted the composite conductivity of SiC/EPDM composites with spherical fillers, in agreement with experimental measurements above φ p .