Permittivity effects of particle agglomeration in ferroelectric ceramic-epoxy composites using finite element modeling

The size dependence of the dielectric constants of barium titanate or other ferroelectric particles can be explored by embedding particles into an epoxy matrix whose dielectric constant can be measured directly. However, to extract the particle dielectric constant requires a model of the composite medium. We compare a finite element model for various volume fractions and particle arrangements to several effective medium approximations, which do not consider particle arrangement explicitly. For a fixed number of particles, the composite dielectric constant increases with the degree of agglomeration, and we relate this increase to the number of regions of enhanced electric field along the applied field between particles in an agglomerate. Additionally, even for dispersed particles, we find that the composite method of assessing the particle dielectric constant may not be effective if the particle dielectric constant is too high compared to the background medium dielectric constant.


I. INTRODUCTION
Dielectric materials for electrostatic capacitors with enhanced energy storage capability are crucial to the future of electrified transportation, power conditioning, and power conversion.A promising class of candidate materials are ceramic nanoparticles such as barium titanate (BaTiO 3 or BTO).BTO is a ferroelectric material commonly used in bulk as a dielectric, as it has a high dielectric constant of about 1500-2000. 1][4][5][6][7] For instance, some have observed a dramatic increase in dielectric constant with decreasing particle size down to a critical value, after which the dielectric constant decreases sharply. 2,6hough exciting, such claims remain unverified.This size-dependence also differs from a more commonly accepted behavior seen in sintered ferroelectric materials in which BTO permittivity begins to decrease and approach zero as the grain size drops below ∼ 1µm. 8,9For this reason, reliable measurements of the dielectric constant of nanosized BTO are essential to the development of improved BTO dielectrics.While we are primarily interested in BTO, the methods and results from this study can readily be generalized to other dielectric materials since we do not consider material properties beyond the dielectric constant.
Rather than attempting to isolate nanoparticles and measuring their dielectric constants, we can more easily arrive at the values indirectly.One approach is to embed particles in an epoxy composite and then measure the resulting dielectric constant.1][12][13][14] While this approach provides a measured value for the composite dielectric constant, additional steps are required to extract the particle dielectric constant.With the pure epoxy dielectric constant and volume fraction of BTO to epoxy known, the particle dielectric constant can be extracted from the composite's measured value using analytical or numerical models.
There are several well-established effective medium approximations (EMAs) that describe the macroscopic dielectric behavior of composites and mixtures of linear dielectric phases.Generally, the approximations relate the effective dielectric constant ε c of a binary composite to the particle dielectric constant ε p and the background medium dielectric constant ε m for a volume fraction of particles δ p .Several EMAs have been developed with various limitations, so we have focused on those that are most applicable to our system.
The Maxwell-Garnett EMA 15 is given by: The main distinction of this EMA is that it assumes all particles are spatially separated and noninteracting, so while perhaps not generally applicable, we used it to validate our finite element (FE) calculations for isolated particles. 16he Bruggeman EMA 17 is given by: This approximation was initially obtained from a model of a regular array of spheres, and is symmetric with respect to the choice of particle and background medium phases. 18While the Bruggeman EMA has the advantage that it does not assume spatially separated particles, it performs poorly when ε p ε m , as is likely true for BTO particles in epoxy.The Jayasundere-Smith approximation was originally derived from the Kerner equation for spheroidal particles, which includes an additional parameter describing the ratio of the average electric fields in the constituent phases along the applied field. 19By replacing this ratio with an analytical solution for electric field interactions between two spheroidal particles, the Jayasundere-Smith approximation is obtained: 20 Because of this modification to the Kerner equation, Jayasundere-Smith partially accounts for the effects of neighboring particles, which are likely important in real composites.Though EMAs are a useful starting point, we have pursued FE modeling both as additional verification and to investigate specific particle arrangements not well-described by any EMA.2][23][24][25] One key case of interest is particle agglomeration, which is common for nanoparticles due to the energy cost of their high surface-to-volume ratio.We have experimentally observed that the dielectric constant of a BTOepoxy composite is often higher than the maximum possible value predicted by EMAs, and we suspect particle agglomeration to be the cause.Sareni et al. modeled agglomerates in two dimensional composites and found that they produced a higher composite dielectric constant than both randomly and regularly arranged isolated particles. 22However, to the extent of the authors' knowledge, there has not been significant modeling done of three-dimensional particle agglomeration in dielectric composites.
In the present study, we use a FE model of BTO-epoxy composite capacitors to study various arrangements of particles, focusing on the effect of particle interaction and agglomeration.Results are compared to EMAs when applicable and interpreted using basic electrostatics.We also explain some features of our model using the method of images and provide insight into the range of dielectric constants for the embedded particles and background medium over which the composite technique will yield useful results.

II. METHODOLOGY
Some terms used throughout are reiterated in Table I.FE calculations were performed using COMSOL Multiphysics version 5.2 and its AC/DC Module, 26,27 which supports electrostatics modeling.Our model, shown in Figure 1, consists of a cube of the epoxy background medium of side length .Thin copper electrodes of cross sectional area 2 and thickness t were placed on the top and bottom faces of the cube.Spheres of radius r, representing BTO nanoparticles, were placed inside the cube, replacing the epoxy material in those regions.The value of r and the number of particles were chosen to achieve the desired volume fraction, δ p , of BTO particles in epoxy.We assumed the particles are spheres rather than irregular spheroids, allowing us to neglect particle orientation and shape.
We chose dimensions = 100 mm and t = 1 mm, though the absolute dimensions of the model are irrelevant since the dielectric constant is dimensionless.However, we assumed a sufficiently large spatial scale that renders valid a classical continuum field theory that neglects discrete atomic or molecular effects.The electrode thickness and choice of copper were arbitrary since the electrodes were simply treated as perfect conductors.The only material property by which the epoxy and BTO materials were distinguished was their dielectric constant, which we assigned.The epoxy background medium dielectric constant ε m was chosen to be 4.5, our average measured value for the particular epoxy that we used (Fibre Glast System 2000 Epoxy Resin).For calculations in which the BTO particle dielectric constant ε p was fixed, we set ε p = 1500.This was within the range of observed bulk values and also large enough that increasing its value produced negligible changes in the composite dielectric constant ε c for all values of δ p and all particle configurations, as explained in Section III A.
The inner surface of the top electrode (where it makes contact with the epoxy) was held at 0 V, while the inner surface of the bottom electrode was held at 10 V. Periodic boundary conditions were imposed on the four remaining faces of the cube so as to model an ideal (infinite) parallel plate capacitor.The default interface condition for the electric displacement field D was applied to the boundaries within the cube between particles and epoxy and between copper and epoxy: 27 where σ f is the free surface charge density between adjacent regions 1 and 2 and n 2 is the normal pointing from region 2 to region 1.In our case, σ f = 0, so this boundary condition simply expresses the continuity of the normal component of the electric displacement field.
For the FE mesh, we first added a triangular mesh on two adjacent non-electrode faces of the cube.Then this mesh was copied to the opposite faces to comply with the periodic boundary conditions.A tetrahedral mesh was applied to the remaining entities in the model.The element sizes were chosen using the predefined options in COMSOL for convenience and consistency.Unless stated otherwise, the "Fine" mesh size was used in our calculations.This was found to be sufficient, as using finer meshes yielded values for ε c within 1% in all cases.
We used the default stationary study type from the electrostatics interface in COMSOL to solve Laplace's equation for dielectric materials in the entire model region, 26,27 where ε 0 is the permittivity of free space, V is the electric potential, and P is the electric polarization.
After obtaining a solution, we calculated the resulting ε c by first integrating the surface charge density over the inner electrode surfaces to determine the charge, Q, on the capacitor.The surface integral was performed over each electrode, and Q was taken to be the average of the absolute value of the two, though they were always within 0.1% of each other.We then found ε c from the capacitance, C, using the parallel plate capacitor equation where V is the potential difference, A is the plate area, and d is the plate separation.Alternatively, the stored energy U is the integral of the electric energy density over the whole volume and is directly proportional to the capacitance: where E is the magnitude of the electric field and ε is the dielectric constant, equal to ε p in the particles and ε m in the epoxy background medium.Then Equation ( 6) can be used to find ε c from C. We found that these two methods consistently yielded values for ε c within 1% (and usually within 0.001%) of each other, so we used the former.

A. Dielectric constant sweep
To test the dependence of ε c on ε p , we swept over a range of values for ε p while holding ε m fixed at 4.5.Figures 2 and 3 show the results of the sweep for an evenly spaced cubic array of 64 particles at various δ p .
In each case here (and in other particle arrangements), ε c increases with increasing ε p until it reaches a plateau value.These values are summarized in Table II and compared to EMAs.The Maxwell-Garnett EMA predicts a lower ε c than FE, though they match to within 1% up to δ p = 0.3.Jayasundere-Smith and Bruggeman, on the other hand, predict a higher ε c than FE, as shown in Figure 3 for δ p = 0.1 and 0.2, with Jayasundere-Smith having more than double the discrepancy of Bruggeman.The discrepancy also tends to grow as δ p is increased.The asymptotic behavior can be explained physically.In the limit that ε p goes to infinity, the electric field inside the particle must go to zero in order for the polarization to remain finite (this follows from the definition of electric susceptibility).Therefore, as ε p increases, the particle behaves more and more like a perfect conductor.Once this state is reached and there is zero electric field inside the particle, increasing ε p further makes no difference.For all other calculations, ε p was fixed at 1500, well into the plateau region, to give an upper bound on ε c in each case.
The disagreement between our FE results and the various EMAs can be interpreted by considering the assumptions made in each.Maxwell-Garnett assumes spatially separated particles, so it is accurate for evenly spaced particles at low δ p but diverges from the FE results at higher δ p , where interactions between particles become significant.Bruggeman and Jayasundere-Smith do not assume spatially separated particles, so they predict a higher ε c when particle interactions are minimized.Since Jayasundere-Smith accounts for some of the particle interactions explicitly, it consistently produces a higher value than Bruggeman.It is also known that the Bruggeman EMA does not work well for ε p ε m , 18 consistent with Figure 3.An important observation to note is that due to the shape of curves in Figure 2, extracting ε p using our model may not actually work for modest values of δ p .Based on its bulk dielectric constant, nanosized BTO may reside on the plateau, in which case the uncertainty in the extracted ε p would be too large to be useful.Wada et al. encountered this issue when extracting a value for ε p using FE modeling (see Figure 13 of Ref. 2).This suggests that an alternative procedure may be necessary for high-permittivity inclusions.

B. Two-particle arrangement
To investigate the effect of particle agglomerations, we examined what happens when two separated particles are brought closer together.We considered two arrangements: the "parallel" arrangement where the particles are separated along the direction of the applied electric field, and the "perpendicular" arrangement where the particles are separated perpendicular to the applied electric field.We doubled the size of the model along one direction to accommodate the two particles, but all other dimensional and electrostatic conditions remained the same.In each case the particles are centered on the long axis of the box, as shown in Figure 4.
The results for the two-particle arrangement are shown in Figure 5 for δ p = 0.1, which displays the same qualitative behavior as other volume fractions.As shown in Figure 4, the particle separation here is taken between the nearest points on the boundaries of the particles rather than their centers, and is normalized to the particle radius.Decreasing the separation to zero in the parallel case causes an increase in ε c , while doing so in the perpendicular case causes a decrease.However, the magnitude of the change in the parallel case is about ten times greater than in the perpendicular case.
These results are consistent with Jayasundere and Smith's analysis of a similar FE model for two spheres (separate from their EMA), in which they found that the particles interact most strongly when aligned along the applied field. 20While the Jayasundere-Smith EMA includes these interactions, it does not account for the slight reduction in ε c in the perpendicular case, which may partly explain why it tends to predict a higher ε c .As the separation increases, both cases approach the Maxwell-Garnett EMA for spatially separated particles.Note that maximum spatial separation occurs at a separation of 1.47 particle radii in Figure 5, and larger values lead to proximity to the electrodes as discussed later in this section.The composite dielectric constant intersects with the Bruggeman and Jayasundere-Smith EMAs at lower separations, where the interaction between particles becomes much more pronounced.
The reasons for these changes in ε c can be observed in the electric field results, shown in Figures 4 and 6 for the parallel and perpendicular arrangements, respectively.In the parallel case, the field between the particles is strongly enhanced as the separation is reduced, which increases ε c .For perpendicular fields, we observe that there are lobes of electric field pointing in opposite directions surrounding the two particles.As the particle separation is decreased, these lobes begin to cancel out in the region between the particles, decreasing ε c slightly.
As shown in Figure 5, the two-particle behavior is symmetric about a separation of 1.47 particle radii, which is where the particle separation is exactly double the distance between each particle and the closest electrode.Essentially, the particle interacts with its image on the opposite side of the electrode, as we discuss in Section III D.

C. Simple cubic arrangement
Building on the model for two particles, we studied a simple cubic (sc) array of particles to investigate the effects of agglomerates.As shown in Figure 7, we considered an evenly spaced array of particles, referred to as the "spread" cubic arrangement, and an array of particles with zero separation between adjacent particles, referred to as the "tight" cubic arrangement, as well as intermediate spacings.Figure 8 shows the results of varying particle separation in a sc arrangement of 64 particles.Beginning with the spread cubic arrangement, we decreased the separation between adjacent particles to zero, reaching the tight cubic arrangement.This results in an increase in ε c .As shown in the results in Figure 9, there is an increased magnitude of electric field in the regions between the particles as they are brought together, as was the case in the parallel two-particle arrangement.
These results agree with our results for two particles.As the particle separation in the cubic arrangement decreases, particles are brought closer together both parallel and perpendicular to the applied field, and these two effects compete with one another.However, since the effect of bringing particles together along the applied field is an order of magnitude greater than that of bringing them together in the perpendicular direction, the net effect is an increase in ε c , as we observe.This behavior suggests that the agglomeration of particles within the composites tends to increase ε c .
While Figure 9 may suggest that the electrostatic energy is stored entirely in the background medium, a non-negligible portion is stored in the particles themselves.The relative amounts of energy stored inside each region can be calculated using Equation (7).There is almost no electric field inside  the particles compared to the background medium, but because ε p ε m , the dielectric constant factor in the energy density balances the storage of the energy more than one might initially expect.For example, for 64 particles at δ p = 0.1, the relative amount of energy stored in particles is 4.2% in the tight cubic arrangement and 0.3% in the spread cubic arrangement.It is unlikely that real agglomerates of BTO nanoparticles are in sc packing; however, this approximation is still quite useful and we found that the tight cubic arrangement sets an upper bound on ε c for a given number of particles and δ p .While particles may be arranged more tightly in a packing such as hexagonal close-packed (hcp), we have seen that particle interactions along the direction of the applied field are most significant, and these are maximized in sc.We compared an hcp arrangement of 64 particles (four close-packed planes) to a tight sc arrangement of 64 particles for δ p = 0.1.The composite dielectric constant was consistently about 10% lower in hcp than in sc, even as both structures were rotated about a direction perpendicular to the applied field to account for some variation in orientation.Previous studies found that sc gave a consistently higher ε c than both body-centered cubic and face-centered cubic (fcc). 21,28Sareni et al. did find hcp to yield a higher ε c than that of sc in some cases, though not dramatically so. 21We note that agglomerates may form in fcc packing, which has the same packing efficiency as hcp.Both of these allow for a higher maximum δ p than sc (0.74 versus 0.52), however, we did not focus on agglomeration in the high δ p regime.
Finite element computations with realistic numbers of particles would exceed the limits of most computers, so we are interested in how the computed ε c depends on the number of particles used in the simulation.For the tight cubic arrangement, we fixed δ p and varied the number of particles used.This can be thought of as taking a single BTO particle and subdividing it into more and more particles arranged in a tight cubic array.The results of this process for various δ p are shown in Figure 10.In each case, ε c increases with increasing number of particles.
The increase in ε c in going from a single particle to eight (2 × 2 × 2 array) is expected because of the large magnitude of the electric field between particles, as in Figure 9.If we assume that ε c is scaled up by some constant multiplicative factor c each time the particles are subdivided, then we would expect the scaling with the number of particles to follow a power law with exponent log 8 c.The results in Figure 10 do appear to follow power laws for large number of particles, with R 2 > 0.95 for each value of δ p if we ignore the first three points in our fit.For the first three agglomerate sizes (1, 8, and 27), all or all but one of the particles are located at the edge of the agglomerate, so we would expect the dielectric constant to be less than that predicted by this simple power law.In contrast to the tight cubic arrangement, we found no dependence on number of particles in the spread cubic arrangement, for reasons discussed below in Section III D.

D. Method of images
As mentioned in previous sections, the interaction between particles and the electrodes in our model may be understood via the method of images, which is commonly used in electrostatics to find the electric field near the surface of a conductor.If we consider a single BTO particle at the center  of the cube, where the electrodes are held at 0 V and 10 V, we can remove the top electrode at 0 V and electrostatically image the system about that plane.This results in a box that is twice as long in the direction of the applied field with electrodes at -10 V (top) and 10 V (bottom), now containing the original particle and a new image particle.The center plane remains at 0 V by symmetry, and the field solution in the original region is also the same.The composite dielectric constant is unchanged because both the applied potential and the electrode separation have doubled.
As shown in Figure 11, we indeed obtain the same solution in these two cases.This means that the field between a particle and an electrode can be framed as the field between that particle and an image particle in a box twice as large.In this sense, our model is effectively periodic in the direction of the applied field as well as the two directions that are given explicit periodic boundary conditions.However, the period in this direction is twice the height of the cubical region.
The method of images accounts for the symmetric behavior observed in the two-particle arrangement, as discussed in Section III B. It also explains why we find that ε c in the spread cubic arrangement does not depend on the number of particles.Since every particle in the spread cubic arrangement is the same distance from its nearest real particle as it is from its nearest image particle in each direction, one particle is mathematically equivalent to eight, 27, and so on.

E. Randomized arrangement with agglomerates
To capture the effects of varying size and number of agglomerates as well as disordered particle positions, we modeled composites with particles randomly interspersed between agglomerates.Let N p denote the total number of particles in a given simulation.The agglomerates were treated as n × n × n tight cubic arrays of particles.We first placed the desired number of agglomerates N a in an evenly spaced cubic array and randomly displaced their positions by up to half of the agglomerate separation in each direction using a uniform random distribution.Then we placed randomly dispersed particles between the agglomerates with equal probability over the remaining space of the cube, excluding the volume of the agglomerates.Examples of the resulting configurations are shown in Figure 12 for N p = 125 particles.For each FE calculation, we consider 10 random configurations to obtain a reasonable average and uncertainty.
We first used this model to examine N p = 512 dispersed particles with no cubic agglomerates (N a = 0).These results are summarized in Table III and compared to EMAs.The average ε c was 2% and 5% higher than the spread cubic arrangement (and the Maxwell-Garnett EMA) for δ p = 0.1 and 0.2, respectively.This difference is because any random deviation of all particles from the spread cubic arrangement will result in a decrease in the separation between some particles parallel to the applied field, increasing ε c .For both the Bruggeman and Jayasundere-Smith EMAs, however, there is better agreement than in the spread cubic case from Section III A. The perfectly uniform spacing of particles in the spread cubic arrangement is extremely unlikely in practice, so the degree of particle interactions is more realistic in the randomized case.We also separately verified that changing the number of randomly dispersed particles at a fixed δ p does not change ε c significantly, at least up to N p = 1000, though the uncertainty tends to decrease with increasing number of particles.Therefore our model for randomly dispersed particles could reasonably be compared to measured ε c , assuming spherical particles and no agglomeration in the samples.
We define the degree of agglomeration as the fraction of the total particles in agglomerates (as opposed to dispersed).For N a n × n × n agglomerates out of N p total particles, this is given by For example, eight 2 × 2 × 2 agglomerates out of 512 total particles corresponds to Fixing N p = 512, we considered each possible degree of agglomeration f a , as well as each possible number of agglomerates N a to achieve that degree (with the restriction that an agglomerate can be no smaller than n = 2).The results of these calculations are shown in Figure 13 for δ p = 0.1 and 0.2.The composite dielectric constant increases roughly linearly with f a for both δ p = 0.1 and 0.2.This can be understood in terms of the number of regions of enhanced electric field between nearest-neighbor particles along the direction of the applied field in an agglomerate (as in Figure 9a).For an n × n × n agglomerate, there are n 2 (n − 1) such regions, which scales roughly as n 3 , the size of the agglomerate.Assuming that each region causes a constant increase in ε c , then we would expect ε c to increase linearly with the agglomerate size.
Furthermore, we can plot the same ε c data from Figure 13 versus the total number of enhanced field regions N e in each arrangement, which is given by for N a n × n × n agglomerates.This relationship is shown in Figure 14.The composite dielectric constant does depend approximately linearly on N e , as suspected.The deviations from linearity seem to be related to the number of agglomerates in each case.For example, for both δ p values there is a dip between N e = 180 (N a = 1, n = 6) and N e = 294 (N a = 1, n = 7) at N e = 256 (N a = 64, n = 2).This implies that there is an additional increase in ε c when more of the enhanced field regions are directly adjacent to one another, as they are for a single agglomerate.Though our results for two particles suggest that interactions between lateral nearest-neighbor particles would cause a slight decrease in ε c , this may be overshadowed by interactions between next-nearest-neighbor particles displaced partially along the applied field.We also observe from Figure 13 that increasing N a for the same f a causes a decrease in ε c , as the agglomerate is divided and the number of regions of enhanced field is reduced.Therefore a single agglomerate surrounded by dispersed particles provides an upper bound on ε c for a given δ p and f a .This model cannot be compared to measured values as easily as the zero agglomeration case due to the assumptions about agglomerate shape and particle number, but it could be useful in checking the validity of measured ε c .For example, a measured dielectric value that greatly exceeds the model value for a single agglomerate might indicate an anomaly with the sample beyond the presence of nanoparticle agglomerates.

IV. CONCLUSIONS
In this study we used a finite element model to investigate the effects of BTO particle agglomeration in BTO-epoxy composites.Our results for dispersed particles tend to agree with the Jayasundere-Smith, Bruggeman, and Maxwell-Garnett EMAs, and the discrepancies are reasonable given the assumptions made by each EMA.However, our results suggest the composite method for extracting useful values of the BTO particle dielectric constant may not be effective for high-permittivity particles at low volume fractions due to the asymptotic behavior of the composite dielectric constant.For interacting particles, we found that bringing particles together generally increases the composite dielectric constant due to an increase in the electric field between particles in the direction of the applied field.
At a fixed volume fraction, the composite dielectric constant due to a simple cubic agglomerate increases as a power law of the number of particles.This result suggests that, in addition to the simplified geometry, the small number of particles in our model makes direct comparison to experiment difficult, unless the power law relationship extends to realistic numbers of particles.We also observed that particles in close proximity to electrodes have the same effect as particles in close proximity to other particles, which was explained using the method of images.Modeling varying size and number of agglomerates surrounded by dispersed particles showed that the composite dielectric constant increases roughly linearly with the number of regions of enhanced electric field between particles and that the composite dielectric constant at a given degree of agglomeration was maximized for a single agglomerate.
This study highlights some challenges to extracting the dielectric constant of discrete ferroelectric nanoparticles within composites.Most notably, the plateau effect for high particle dielectric constant indicates that even if the composite dielectric constant is measured with reasonable precision, the corresponding uncertainty in the particle dielectric constant from modeling may be too large to be informative.Estimates made using this method should therefore be met with some caution.While this technique may be more challenging to apply to high permittivity ferroelectrics like BTO, it could be more easily applied to lower permittivity materials.Also, the plateau is pushed to larger particle dielectric constants for larger volume fractions, so it would be desirable to examine composites with higher volume fractions.However these can be difficult to synthesize while maintaining dispersed particles.A natural extension of our work would be modeling more realistic agglomerates based on experimentally observed packings of nanoparticles.This may allow for an alternative procedure for extracting single particle dielectric constants where some agglomeration is accounted for, thus enabling measurements at higher volume fractions of particles.We also note that agglomeration may be preferable for capacitor applications, as long as electrical breakdown can be avoided, due to the enhanced energy storage within the agglomerates.We expect these results will be useful for understanding and modeling the performance of dielectric composite materials.

FIG. 2 .
FIG. 2. Composite dielectric constant ε c versus particle dielectric constant ε p for various volume fractions, δ p .An evenly spaced cubic array of 64 particles was used.Dotted lines indicate Maxwell-Garnett EMA and solid lines indicate FE results.

FIG. 4 . 7 Kaufman
FIG. 4. Natural logarithm of the z-component of the electric field in V/m (along the applied field) in the parallel arrangement at separation of (a) 0.1 radii and (b) 1.47 radii.δ p = 0.1.The perpendicular components of the electric field (not shown) are very small.

FIG. 6 .
FIG. 6. y-component of the electric field in V/m (perpendicular to the applied field) in the perpendicular arrangement at separation of (a) 0.1 radii and (b) 1.47 radii.δ p = 0.1.The parallel component of the electric field (not shown) is approximately the applied electric field.

FIG. 9 .
FIG. 9. Electric field magnitude in V/m for a cubic array of 64 particles at δ p = 0.1 and a separation of (a) 0.1 radii and (b) 1.47 radii, where 1.47 radii corresponds to the spread arrangement.Notice the difference in color scale between (a) and (b).

FIG. 10 .
FIG. 10.Log-log plot of ε c versus number of particles in a tight cubic arrangement for various δ p fit to power laws (dotted lines), ignoring the first three points.Fitted values of c range from 1.03 to 1.19.

FIG. 11
FIG. 11. z-component of the electric field in V/m (along the applied field) for a (a) single particle and (b) the equivalent arrangement from the method of images.The slight decrease in field amplitude around the midplane in (b) may be an artifact of the calculation.

FIG. 13 .FIG. 14 .
FIG.13.Composite dielectric constant ε c versus degree of agglomeration f a for N p = 512 particles and varying number of agglomerates N a at δ p = 0.1 and 0.2.Each point is the average of 10 random configurations.Error bars indicate ± one standard deviation.

TABLE I .
Important term symbols and definitions.

TABLE II .
Plateau values of ε c (at ε p = 1500) for various δ p , calculated by FE and EMAs.Parentheses indicate percent difference from FE result.

TABLE III .
Average ε c and standard deviation from randomized FE model for N p = 512 with N a = 0, compared to EMAs.Parentheses indicate percent difference from FE result.