In this paper, the simulation and modeling of a polyethylene/clay nanocomposite were undertaken to predict the nanocomposite's dielectric behavior and to help design a nanocomposite material with optimum electrical properties for electrotechnical or electronic applications. A 3D simulation model using the finite elements method was employed in order to study the effective permittivity and electric field distribution of twophase nanocomposite materials for ordered and random distributions of inclusions in a lowloss host matrix such as polyethylene. The influence of the dispersion of reinforcing particles, and of the permittivity and radius of the inclusions, was analysed. The simulation results were compared with alternative, known theoretical solutions obtained from classical models, and were found to be in good agreement. The numerical results also indicate that for fixed volume fractions of nanoparticles the effective permittivity of the mixture, for ordered and random distributions, does not vary with the degree of dispersion. The variation of the effective permittivity with the particle radius is shown, using numerical data, to agree with the analytical modules.
1. INTRODUCTION
In recent years, nanocomposites have attracted a great deal of attention from many researchers. Several studies have shown that the incorporation of nanometric inclusions in a polymer matrix can often significantly improve its mechanical, dielectric and optical properties when compared to a pure polymer matrix
[1

9]
, provided that the particles are reasonably well dispersed. In the electrical insulation field, polyethylene is extensively used in medium/high voltage electrical cables due to its excellent dielectric properties, very low dielectric losses and high intrinsic breakdown strength. In order to improve these properties and be able to meet new needs, such as those relating to insulation systems in DC power cables, composite materials consisting of a can be used; this results in a possible increase in breakdown strength, dielectric endurance, thermal conductivity, and allows the electrical conductivity to be tailored to avoid space charge accumulation in DC applications.
The effective permittivity of composite materials generally depends on the material microstructure, which includes the volume fraction, as well as the shapes and types of component. For some specific conditions (i.e. when the mixture has a periodic welldefined structure), the effective permittivity can be calculated from the analytical solution of the field distribution. This has resulted in various analytical models, such as the commonly used laws of mixtures. However, if the structure is disordered, as may be the case for a real compounded composite, analytical models cannot be used to accurately estimate the material’s effective permittivity. However, it should be noted that in this case, it can still be measured experimentally.
Numerical methods have been developed and calculation capabilities expanded to evaluate the effective dielectric permittivity of composite materials. Numerical simulations have in fact grown to represent another method for dealing with a large proportion of the physical problems of composite materials.
In this paper, the effective permittivity of a polyethylenebased nanocomposite is calculated by numerical simulation using the Comsol Multiphysics software, which is based on the finite element method (FEM). The influence of dispersion as well as the variation of the permittivity and radius (or the volume fraction) of inclusions on effective permittivity is studied. The electric field and polarization distribution in the nanocomposite materials are also reported in this paper. A comparison between the numerical results and those of the analytical models is presented.
2. ANALYSIS
The effective dielectric permittivity of a multiphase material can be estimated either by analytical or numerical methods.
 2.1 Analytical models
The starting point of most analytical approaches to estimating the effective properties of heterogeneous media is the solution of the single inclusion problem for which a constant field Eo along the zdirection is applied at a distance from the inclusion. This approach has been detailed in many textbooks and review papers
[10

20]
. In spherical coordinates, the solution of the Laplace equation for a spherical inclusion of radius R is given by:
where ψ(r,θ) is the electrical potential, r is the radial coordinate, θ is the angle between the position vector and the zcoordinate, and ε
_{1}
and ε
_{2}
are the permittivities of the inclusion and the matrix, respectively. It should be noted that the same equations hold in steadystate AC conditions in which the potential would be a phasor and the permittivity can take complex values, including a possible conductivity term. The electrical field along the zdirection inside the inclusion according to (1) is given by:
A similar calculation can be made for the more general case of an ellipsoidal inclusion leading to:
where A
_{1}
is the depolarization factor along the ellipsoid principal axis parallel to the electrical field
[11]
. For spherical particles, A
_{1}
= A
_{2}
= A
_{3}
= 1/3 and (3) is identical to (2). By definition, the effective dielectric constant of a twocomponent heterogeneous linear material can be defined by
[21]
:
where ε
_{c}
is the effective permittivity, q
_{1}
and q
_{2}
are the volume fraction of the matrix and the inclusion, and the brackets denote averages over phase 1, phase 2 and over the material’s volume. An analytical calculation of the electrical field in a composite material can only be done if the minority phase is present in a small concentration and has regular shape inclusions. A number of results can be found in which an exact solution for several matrix systems with periodic arrangements of regular inclusions is obtained
[22
,
23]
. In the case of a dilute suspension of ellipsoidal shape inclusions with a permittivity ε
_{2}
in a continuum matrix of permittivity ε
_{1}
, it is possible to use the solution of the singleinclusion problem (equation (3)), assuming that the field E
_{o}
is equivalent to the average field in the matrix (phase 1). This leads to
[11]
:
A similar procedure leads, for randomly oriented inclusions, to:
where A
_{i}
is the depolarization factor for the ith axis of the ellipsoid. For spherical particles, A
_{1}
= A
_{2}
= A
_{3}
= 1/3. In the case of spheroids for which two axes are equal (a = b ≠ c), the analytical expressions for Ai for oblate spheroids (disklike spheroids with a = b > c) and prolate spheroids (needlelike spheroids with a = b < c) can be found in literature
[10]
.
Another approach, the effective medium approximation, also relies on the solution of the singleinclusion boundary value problem. The generalization of the Maxwell approximation for ellipsoidal inclusion leads to:
for a twophase composite consisting of a perfectly oriented ellipsoidal inclusion (phase 2) inside a matrix (phase 1). It can be shown that this equation is equivalent to equation (5), and the randomly oriented case would be equivalent to (6). This is also the exact solution of the coatedspheres model
[24]
.
A different approach in the effective medium approximation family is the selfconsistent approximation, which was originally developed by Bruggeman
[25]
. It leads to a slight modification of (7):
which in turn leads to a quadratic equation for the effective permittivity. For the case of the randomly oriented ellipsoidal, equation (8) can be written
[26]
as:
Finally, using a symmetric integration technique
[27]
, it can be shown that the Looyenga equation for randomly oriented ellipsoids, independent of their shape, is given by:
As mentioned previously, all the above equations also hold in steadystate AC conditions. In this case, the electric fields can be replaced by their respective phasors and the permittivity can be replaced by their complex representation. Accordingly, equations (5) to (10) can also be used to predict a composite dielectric response, i.e., the variation of the effective complex permittivity as a function of frequency, by replacing ε
_{1}
and ε
_{2}
by their frequencydependent complex representations.
 2.2 Numerical models
The effective permittivity of the dielectric mixture can be calculated once the electric field vector is known inside the material. In a purely electrostatic case, it can be calculated by solving the Poisson's equation given by:
where ε
_{r}
, ε
_{0}
and Ψ are the relative permittivity, the vacuum permittivity and the electrical potential, respectively. ρ is the charge density. For the neutral condition (ρ = 0), and if we take into account a possible conductivity σ and dielectric losses, then (11) can be written in the steady state more generally as
[28

30]
:
where the complex permittivity is given by:
Once the material microstructure and the properties of each phase are known, (11) or (12) can be numerically solved by the finite elements method (FEM) using a commercially available package like COMSOL Multiphysics. After the field distribution is numerically evaluated, the effective permittivity of the dielectric mixture can be calculated in several ways
[30

38]
. In this paper, the effective dielectric permittivity of the composites was calculated by using the averages of the electric field and dielectric displacement values. Therefore, the effective dielectric permittivity can be expressed as follows:
where <
D
> denotes the mean value of the electric displacement field and <
E
> is the mean value of the electric field over the material. Both averages were taken in the direction of the applied field.
3. RESULTS AND DISCUSSION
 3.1 Simulation setup
The morphology of polyethylene/clay nanocomposite in the nanometric scale observed by scanning electronic microscopy (SEM) is shown in
Fig. 1
. As can be seen, two components are presented in the nanocomposite. This twophase nanocomposite consists of 5 wt% of nanoclay particles dispersed in a polyethylene matrix. A comparison between the morphology of polyethylene with 5 wt% nanoclay, PE/OMMT (
Fig. 1(a)
), and polyethylene with 5 wt% nanoclay and 10 wt% of compatibilizer, PE/OMMT/PEMA (
Fig. 1(b)
), shows that the density and size of aggregates were decreased in the compatibilized nanocomposite. This improvement in the nanoclay dispersion is due to the presence of the polar compatibilizer, PEMA
[8]
. In this paper, the dielectric permittivity of the polyethylene matrix and that of the nanoclay reinforcing fillers are taken to be as follows: ε
_{1}
=2.3 and ε
_{2}
=4.4, for the matrix and the filler, respectively. A singleinclusion 3D model cell was used as a first approach, as shown in
Fig. 2
. The inclusion was assumed to have a spherical form to simulate the shape of a nanoclay particle, and the radius of this was later calculated in order to meet the requirement of a 5% volume fraction of the particles. The bottom face of the cube was set to a constant potential (Ψ = 1 V), and the opposite face was set to ground (Ψ = 0). The other faces of the cube were set to periodic conditions.
SEM micrograph of (a) PE/OMMT and (b) PE/OMMT/PEMA nanocomposites.
Unit cell model of a singleinclusion twocomponent periodic composite material.
 3.2 Ordered distribution
In order to study the effect of the dispersion of nanoclay particles on dielectric properties, such as electric field distributions and effective permittivity, four geometries (with the same volume fraction q=0.05) were drawn with 1, 8, 27 and 64 spheres corresponding to particle normalized radii of 0.229, 0.114, 0.076 and 0.057 (
Figures 3(a)

3(d)
).
Electric field distribution in the nanocomposite with (a) 1 particle, (b) 8 particles, (c) 27 particles, and (d) 64 particles (ordered distribution).
The surface plots of the electric field distribution in the 3D model were obtained by FEM simulations, and are shown in
Fig. 3
.
As can be predicted by the singleinclusion solution, a field enhancement is present at the interface between the particle and the matrix on the bottom and top sides in the zdirection, the field is almost constant and smaller within the inclusion. The maximum value of the electric field E
_{max}
increases as the degree of dispersion of the nanoclay particles increases, and the highest value is obtained when the number of spheres is 27.
The polarization vector (in C/m
^{2}
) distributions for all 4 geometries are shown in
Figure 4
, along with an enhancement at the matrixparticle interfaces in the zdirection. The highest value of polarization was found at the surface of the particles in the cell with 64 spheres. It is evident from the images that the polarization increases as the degree of dispersion of the nanoclay particles increases.
Polarization field distribution in the nanocomposites with (a) 1 particle, (b) 8 particles, (c) 27 particles, and (d) 64 particles (ordered distribution).
 3.3 Random distribution
In the above section, the uniform reinforcing particles are considered to be orderly dispersed in a polyethylene matrix. In real materials, nanoclay particles or any reinforcing filler are distributed more or less randomly in the polymer matrix with agglomerate size distributions.
Figures 5(a)

Figures 5(d)
show four geometries and electric field distribution plots for different numbers of randomized nanoclay particles. The volume fraction of nanoclay particles was set at 5% and dispersed in a unit cell.
The corresponding effective permittivity and normalized maximum electric field E
_{max}
/E
_{0}
(where E
_{0}
=1 kV/mm) of the nanocomposite materials are plotted in
Figs. 6
and
Figs. 7
, respectively. As can be seen from
Fig. 6
, the effective permittivity for ordered distributions does not vary significantly with the number of particles. This is due to the fact that the volume fraction of nanoclay particles was constant and set at 5%. On the other hand, the effective permittivity of the random distribution was found to be higher than that of ordered distributions. This change is due to a decrease of the distance between inclusions in the random distribution.
Electric field distribution in the nanocomposites with (a) 1 particle, (b) 8 particles, (c) 27 particles, and (d) 64 particles (random distribution).
Figure 7
shows the normalized maximum value of the electric field in the zdirection as a function of the particle dispersion. At a low particle number, it can be seen that the normalized maximum electric field is almost identical for the ordered distribution as for the random cases, and we can also see it increases as the quality of dispersion is improved. For particle numbers higher than 27, the normalized maximum value of the electric field in the random distribution decreases slightly as the quality of dispersion increases, while in the ordered distribution, the electric field is found to increase considerably with an increase in the number of particles.
Comparison of effective permittivity for random and ordered nanoparticle distributions.
Comparison of normalized maximum electric field for random and ordered nanoparticle distributions.
 3.4 Effect of the permittivity of the inclusion on effective permittivity
In this section, the two geometries of 64 spheres of the ordered and random distributions presenting the distribution of the electrical potential simulated by FEM are shown in
Figs. 8(a)
and
8(b)
to show the effect of a variation of the dielectric permittivity of inclusions on the effective permittivity of the ordered and random nanocomposites. The resulting effective permittivities obtained from FEM are compared with those calculated from equation (7), the MaxwellGarnett model, equation (9), the Brug In this section, the two geometries of 64 spheres of the ordered and random distributions presenting the distribution of the electrical potential simulated by FEM are shown in
Figs. 8(a)
and
Figs. 8(b)
to show the effect of a variation of the dielectric permittivity of inclusions on the effective permittivity of the ordered and random nanocomposites. The resulting effective permittivities obtained from FEM are compared with those calculated from equation (7), the MaxwellGarnett model, equation (9), the Bruggeman model, and equation (10), the Looyenga model.
Electric potential distribution in the nanocomposite with 64 particles: (a) ordered distribution and (b) random distribution. The red lines present the electrical field stream lines.
Numerical and theoretical results of the effective permittivity as a function of the dielectric permittivity of inclusions.
Figure 10
shows the effect of varying the permittivity of the inclusion ε
_{2}
on the normalized maximum value of the electrical field in the zdirection. As expected, for ordered and random distributions, the maximum value of the electrical field increases as the value ε
_{2}
increases.
Normalized maximum electric field for ordered and random nanoparticle distributions as a function of the dielectric permittivity of inclusions.
 3.5 Effect of the radius (and volume fraction) of the inclusion on effective permittivity
The distributions of the electric potential simulated by FEM for the two geometries for the singleinclusion case, of the ordered and random distributions, are shown in
Fig. 11
. In the simulation, the permittivity of the polymer matrix ε
_{1}
is still assumed to be 2.3, and the permittivity of the inclusion ε
_{2}
is set at 4.4. To study the effect of nanoclay loading on the effective permittivity ε
_{c}
, the normalized radius of the inclusion r/r
_{0}
(where r
_{0}
=1 μm) was varied from 0.010 to 0.400, corresponding to a volume fraction of nanoparticles ranging from 0.0004% to 26.81%.
Electric potential distribution in the nanocomposite with 1 particle: (a) ordered distribution and (b) random distribution. The red lines represent the electrical field stream lines.
The effective permittivity of the nanocomposites as a function of the radius of nanoclay particles is plotted in
Fig. 12
. It can be observed that the effective permittivity increases with the nanoclay particle radius, and it is evident from this figure that the effective permittivity of the ordered and random distributions obtained from FEM is quite similar to that calculated with the MaxwellGarnett model, and is also very close to those of the Bruggeman and Looyenga models. These results confirm the validity of our simulation method.
Variation of the effective permittivity with normalized radius of the spherical nanoclay particles.
4. CONCLUSIONS
A 3D simulation model using the finite elements method was developed in order to study the effective permittivity and electric field distribution of polyethylene/clay nanocomposite materials for electrical applications. An enhancement of the electric field and polarization was observed as the degree of dispersion of the nanoclay particles increased; however, the effective permittivity of the nanocomposites was not affected by improving the quality of dispersion of nanoclay particles in the host matrix. The numerical results indicate that the MaxwellGarnett model is appropriate for evaluating the effective permittivity of ordered distributions, while the Bruggeman Symmetry model remains the most suitable for calculating the effective permittivity for random distributions. It was observed that in both distributions the maximum value of the electrical field and the effective permittivity ε
_{c}
increase as the value ε
_{2}
increases. Finally, this numerical model can be extended to design nanocomposite materials with optimum dielectric properties for electrotechnical or electronic applications.
Acknowledgements
The authors are grateful for financial support from the National Sciences and Engineering Research Council of Canada (NSERC).
Hoyos M
,
Garcia N
,
Navarro R
,
Dardano A
,
Ratto A
,
Guastavino F
,
Tiemblo P
2008
Journal of Polymer Science Part B: Polymer Physics
[DOI: http://dx.doi.org/10.1002/polb.21464]
46
1301 
Osman M. A
,
Rupp J.E.P
,
Suter U. W
2005
Polymer
[DOI: http://dx.doi.org/10.1016/j.polymer.2004.11.112]
46
1653 
Utracki L
,
Kamal M
2002
Arabian Journal Science & Engineering Special Issue
27
43 
Kornmann X
,
Lindberg H
,
Berglund L. A
2001
Polymer
[DOI: http://dx.doi.org/10.1016/S00323861(00)003463]
42
1303 
Han Y
,
Wang Z
,
Li X
,
Fu J
,
Cheng Z
2001
Current Trends in Polymer Science
[DOI: http://dx.doi.org/10.1016/S13601385(00)018306]
6
1 
Kawasumi M
,
Hasegawa N
,
Kato M
,
Usuki A
,
Okada A
1997
Macromolecules
[DOI: http://dx.doi.org/10.1021/ma961786h
30
6333 
Vaia R. A
,
Teukolsky R. K
,
Giannelis E. P
1994
Chemistry of Materials
[DOI: http://dx.doi.org/10.1021/cm00043a025]
6
1017 
Zazoum B
,
David E
,
Ngô A
2013
Journal of Nanotechnology
2013
10 
David E
,
Frechette M
,
Zazoum B
,
DaranDaneau C
,
Ngô A. D
,
Couderc H
2013
Journal of Nanomaterials
Torquato S
2002
Random Heterogeneous Materials: Microstructure and Macroscopic Properties
Interdisciplinary Applied Mathematics
Springer
[DOI: http://dx.doi.org/10.1007/9781475763553]
Banhegyi G
1986
Colloid and Polymer Science
[DOI: http://dx.doi.org/10.1007/BF01410321]
264
1030 
Cret R
,
Cret L
2004
Journal of Optoelectronics and Advanced Materials
6
1045 
Tuncer E
,
Serdyuk Y. V
,
Gubanski S. M
2002
IEEE Transactions
9
809 
Peng J. H
,
Yang J. J
,
Huang M
,
Sun J
,
Wu Z. Y
2009
Frontiers of Materials Science in China
[DOI: http://dx.doi.org/10.1007/s1170600900152]
3
38 
Koledintseva M. Y
,
Patil S. K
,
Schwartz R. W
,
Huebner W
,
Rozanov K. N
,
Shen J
,
Chen J
2009
Dielectrics and Electrical Insulation
IEEE Transactions
16
793 
Barber P
,
Balasubramanian S
,
Anguchamy Y
,
Gong S
,
Wibowo A
,
Gao H
,
Ploehn H. J
,
Zur Loye H. C
2009
Materials
[DOI: http://dx.doi.org/10.3390/ma2041697]
2
1697 
Blanchard C
,
Portí J. A
,
Morente J. A
,
Salinas A
,
Navarro E. A
2007
J. Appl. Phys.
[DOI: http://dx.doi.org/10.1063/1.2779216]
102
064101 
Karkkainen K. K
,
Sihvola A. H
,
Nikoskinen K. I
2000
Geoscience and Remote Sensing, IEEE Transactions
38
1303 
Sihvola A. H
,
Kong J. A
1988
Geoscience and Remote Sensing, IEEE Transactions
26
420 
Beran M
1968
M. Beran, N. York
Emets Y. P
1998
Journal of Experimental and Theoretical Physics
[DOI: http://dx.doi.org/10.1134/1.558701]
87
612 
Tuncer E
,
Gubanski S. M
,
Nettelblad B
2001
J. Appl. Phys.
[DOI: http://dx.doi.org/10.1063/1.1372363]
89
8092 
Hashin Z
,
Shtrikman S
1962
J. Appl. Phys.
[DOI:http://dx.doi.org/10.1063/1.1728579]
33
3125 
Bergman D. J
1982
Annals of Physics
[DOI: http://dx.doi.org/10.1016/00034916(82)901762]
138
78 
Lal K
,
Parshad R
1973
J. Phys. D: Appl. Phys.
[DOI:http://dx.doi.org/10.1088/00223727/6/11/311]
6
1363 
Looyenga H
1965
Physica
[DOI: http://dx.doi.org/10.1016/00318914(65)900455]
31
401 
Tuncer E
,
Nettelblad B
,
Gubaski S. M
2002
J. Appl. Phys.
[DOI: http://dx.doi.org/10.1063/1.1505975]
92
4612 
Tuncer E
,
Serdyuk Y
,
Gubanski S
2001
Comparing Dielectric Properties of Binary Composite Structures Obtained with Different Calculation Tools and Methods, Electrical Insulation and Dielectric Phenomena
2001 Annual Report. Conference on. IEEE
665 
Brosseau C
,
Beroual A
2001
J. Phys. D: Appl. Phys.
[DOI: http://dx.doi.org/10.1088/00223727/34/5/307]
34
704 
Tuncer E
,
Gubański S. M
,
Nettelblad B
2001
J. Appl. Phys.
[DOI: http://dx.doi.org/10.1063/1.1372363]
89
8092 
Jebbor N
,
Bri S
2012
Journal of Electrostatics
[DOI: http://dx.doi.org/10.1016/j.elstat.2012.05.007]
70
393 
Venkatesulu B
,
Jonsson B.L.G
,
Edin H
,
Norgren M
2013
Dielectrics and Electrical Insulation, IEEE Transactions
20
177 
Fa X
,
Rcas X
,
Cret C
,
Petreus X. R. D
,
Palaghit X. X
2010
Modeling and Simulation of Dielectric Mixtures Using Finite Elements Method, Design and Technology in Electronic Packaging (SIITME)
2010 IEEE 16th International Symposium
305 
Cret R
,
Darabant L
,
Farcas C
,
Turcu A
2011
Considerations about the influence of some factors related to the geometric characteristics of inclusions on effective permittivity of dielectric mixtures, Advanced Topics in Electrical Engineering (ATEE)
2011 7th International Symposium
1 
Karkkainen K
,
Sihvola A
,
Nikoskinen K
2001
Geoscience and Remote Sensing
IEEE Transactions
39
1013 
Sihvola A. H
,
Kong J. A
1988
Geoscience and Remote Sensing, IEEE Transactions
26
420 
Nilsson F
,
Gedde U. W
,
Hedenqvist M. S
2011
Composites Science and Technology
[DOI: http://dx.doi.org/10.1016/j.compscitech.2010.11.016]
71
216 