Examensarbete 30 hp
Oktober 2019
Electromagnetic Homogenization-
simulations of Materials
Teknisk- naturvetenskaplig fakultet UTH-enheten Besöksadress: Ångströmlaboratoriet Lägerhyddsvägen 1 Hus 4, Plan 0 Postadress: Box 536 751 21 Uppsala Telefon: 018 – 471 30 03 Telefax: 018 – 471 30 00 Hemsida: http://www.teknat.uu.se/student
Abstract
Electromagnetic Homogenization-simulations of
Materials
Julia Törnqvist
This thesis aims to determine the distribution of the relative permittivity for random mixtures of material using electromagnetic simulations. The algorithm used in the simulations is the FDTD method which solves Maxwell's equations numerically in the time-domain. The material is modeled as randomly shaped particles with radius 12 ± 10 micrometre in x- and y-direction and radius 3 ± 1 micrometre in z-direction. The scattering parameters from the transmitted and reflected electric field when a plane wave interacts with the material are measured. The relative permittivity is determined from the scattering parameters using the iterative Baker-Jarvis method. The simulations shows that both the distribution and the value of the relative permittivity is low when the particles have non conducting layers to force interruptions to prevent percolation, a conducting path between the particles. The most important result is of the kind where the simulations do not have any boundaries to prevent percolation. These simulations reflects how the relative permittivity distributes in real measurements. It is established that the value of the relative permittivity has a large distribution and also that percolation occurs because of the periodic structures.
Ämnesgranskare: Jan Isberg
Sammanfattning
När elektromagnetiska fält infaller på olika blandningar av material kan fälten bli modifierade. Olika bland-ningar av material kan absorbera, dämpa, blockera eller omdirigera elektromagnetisk strålning. Elektromag-netiska simuleringar görs ofta för att bestämma radarmålarea för delar av, eller hela flygfarkoster. Därför är det viktigt att ha kunskap om materialens egenskaper så som permittivitet och permeabilitet.
Acknowledgements
I would like to thank my supervisors Torleif Martin and Tomas Lundin at Saab for being supportive and for guiding me throughout this work. Thanks for all the valuable inputs and answers to my questions and for guiding me in the right directions whenever needed.
I would also like to thank Jan Isberg at the department of Electricity at Uppsala University for being subject reviewer for this thesis.
Contents
1 Introduction 6
1.2 Background . . . 7
1.3 Thesis problem description . . . 7
1.4 Purpose and goals . . . 7
1.5 Outline . . . 7 2 Waveguide Theory 9 2.1 Maxwell’s Equations . . . 9 2.2 Constitutive Relations . . . 10 2.2.1 Electromagnetic Materials . . . 10 2.3 Waveguides . . . 11 2.4 TEM waves . . . 13 2.5 Coaxial line . . . 13 3 Material Characterization 15 3.1 Estimation of material parameters . . . 15
3.1.1 Scattering parameters . . . 16
3.1.2 Nicolson-Ross-Weir Algorithm . . . 16
3.1.3 Baker Jarvis Method. . . 17
4 Simulation tools 19 4.1 Finite-Difference Time-Domain method . . . 19
4.2 Simulation set-up . . . 20
5 Simulations 22 5.1 Simulation specifications . . . 23
5.2 Cluster-Cluster aggregation model . . . 23
5.3 Particle model based on ellipsoids. . . 25
5.3.1 Modeling a single particle . . . 25
5.3.2 Modeling of several particles . . . 26
Appendix A 47
Appendix B 48
Chapter 1
Introduction
Determine the electromagnetic properties of a material is challenging with theories based on Maxwell’s equations, waveguides and material characterization. The material parameters are a macroscopic model of microscopic material distributions. Beginning with Maxwell’s equations that can be applied to both micro-and macroscopic calculations the electric micro-and magnetic fields are defined [10]. The fields and also the electric and magnetic flux forms constitutive relations which appears in different forms depending on the material and its permittivity and permeability [2].
The permittivity and permeability of a material can be found when a sample is brought under illumination of an incident electromagnetic field and the reflected and transmitted fields are measured. This can be done using waveguides or using a free space measurement setup with antennas [7] [19]. Figure1.1shows an example of a measurement setup to measure the transmission and reflection scattering parameters using a coaxial cable.
Figure 1.1: Diagram of Experimental Fixture [19].
The permittivity can be determined by measuring the scattering parameters using the transmission/reflection method. Estimating the scattering parameters can be done with several different methods such as the Nicolson-Ross-Weir algorithm [8] and with Baker-Jarvis method [4].
1.2
Background
Different mixtures of material modifies the electromagnetic fields interacting with them. The material can absorb, block, attenuate or simply redirect the electromagnetic radiation. Electromagnetic simulations are made to determine radar cross-section (RCS) for parts of, or entire aircrafts. It is therefore important to have knowledge of how the materials behave when exposed to external electromagnetic fields and that can be possible by exploring the material parameters permittivity and permeability µ.
1.3
Thesis problem description
From previous measurements carried out at Saab it has been observed that for material mixtures containing particles with high conductivity, the spread of permittivity is large. The idea is to recreate similar measure-ments through electromagnetic simulations to investigate whether the simulations is consistent with previous results.
1.4
Purpose and goals
This thesis is carried out at Saab Aeronautics in Linköping and aims to determine the relative permittivity for random mixtures of material by using electromagnetic simulations. The purpose is to study the spread of the relative permittivity depending on the shapes of the particle and the random distribution of particles in material mixtures with controlled content of particles.
1.5
Outline
In chapter 2 and chapter 3 the waveguide theory and material characterization are explained in detail. Chapter 4 is about the simulation tools. Scattering parameters are explained more in detail together with the Finite-Difference Time-Domain (FDTD) method which is a solver for electromagnetic wave equations in time domain.
Figure 1.2: Particles in a material shown on a microscopic scale.
The beginning of chapter 5 includes the modeling step by step to the final model where each particle is based on an ellipsoid with random dimensions in each direction. The second part in chapter 5 describes the test simulations where the shapes of the particles is determined in the shape analysis together with the distri-bution analysis where the distridistri-bution of the sizes of the particles is investigated. The last part of chapter 5 contains the main simulations where the different cases depends on the particle content and the size of the FDTD volume. To prevent percolation, a network of conductive pathways created by the particles, a non conducting layer around each particle is introduced to eliminate connection. To prevent percolation to repeat periodically a gap through the FDTD volume is also introduced.
Chapter 2
Waveguide Theory
2.1
Maxwell’s Equations
Faraday’s law of induction, Ampére’s law and Gauss’ laws for the electric and magnetic fields are the fun-damental blocks of Maxwell’s equations2.1, here formulated in SI unit convention. [11]
∇ × E = −∂B ∂t ∇ × H = J +∂D ∂t ∇ · D = ρ ∇ · B = 0 (2.1)
Where E is the electric field [V/m] and B is the magnetic flux density [Wb/m2], or [T]. H is the magnetic field [A/m] and D is the electric flux density [C/m2]. The quantities of ρ and J describes the volume charge density [C/m2] and the electric current density [A/m2] respectively.
Maxwell’s equations can be applied to both microscopic and macroscopic calculations. Applying the equa-tions on macroscopic calculaequa-tions results in two new auxiliary fields, the electric flux density, also called the displacement field D and the magnetic field H. With the newfound fields the behaviour of materials can be described in large-scale [10].
Using the Fourier transform on Maxwell’s equations combined with the constitutive relations the equations can be expressed in frequency domain instead of time domain.
∇ × E = −jωµH ∇ × H = jωE
∇ · E = 0 ∇ · H = 0
(2.2)
2.2
Constitutive Relations
Constitutive relations describes the relation between the flux densities D and B and the fields E and H [5]. The constitutive relations can appear in different forms depending on the material. The constitutive relations for linear material are not a function of the applied field whereas the parameters for non-linear material are. Homogeneous materials are uniformed throughout and the constitutive relations are not a function of position as they are for inhomogeneous materials. Isotropic materials has the same properties in all directions and non-isotropic materials has different properties in different directions. The constitutive relations for non-isotropic materials vary with the direction of the applied field. When the constitutive re-lations depends on the frequency, the materials are dispersive [2]. For vacuum the constitutive relations are expressed as
D = 0E, B = µ0H (2.3)
where the constants 0 = 8.854 × 10−12 [F/m] is the permittivity in vacuum and µ0 = 4π × 10−7 [H/m] is the permeability in vacuum. The speed of light and the characteristic impedance can be derived from the quantities of permittivity and the permeability in vacuum.
c0= 1 √ µ00 = 3 × 10 8m/sec Z0= r µ0 0 = 377Ω (2.4)
For homogeneous isotropic dielectric and magnetic materials the constitutive relations is described as
D = 0rE, B = µ0µrH (2.5)
where r and µr are the relative permittivity and permeability respectively.
2.2.1
Electromagnetic Materials
Dielectric materials are electric insulators which means that when applying a voltage there will be no cur-rent flowing through the material since there is no free charges that can travel. Dielectrics can however be polarized when when applying a external electric field which gives the material the ability to store energy. Polarization means that the negative charges, the electrons, shift from their initial position towards the pos-itive voltage and creates potential energy. The ability to store energy from an external field is represented by the dielectric constant also called the permittivity, [16]. The permittivity of a material is defined as
where 0is the permittivity of free space and ris the relative permittivity. The relative permittivity depends on several different factors such as temperature and frequency.
Magnetic materials has the ability to repel or attract other substances when being exposed to external mag-netic fields. The permeability, µ, describes the degree of magnetization in materials from external magmag-netic fields. Magnetization is the ability to maintain the magnetic fields within materials. The permeability of a material is defined as
µ = µ0µr (2.7)
where µ0 is the permeability of free space and µr is the relative permeability.
The electric and magnetic properties of a material can be determined by the permittivity and permeability. The permittivity and permeability describes the interaction between an applied electromagnetic field and material. Equation2.8describes the relative permittivity using frequency dependent complex parameters.
r(ω) = (ω) 0 = 0(ω) − j00(ω) 0 = 0r(ω) − j00r(ω) (2.8)
For material with conductivity σ the relative permittivity is = r+ σ
iω0 = r− j σ
ω0 (2.9)
Where the real part of the above equation 0 is the relative permittivity which express the energy storage, the exchange of energy between the material and field. The imaginary part 00is the loss factor and describes the dissipation of energy which occurs when the material absorbs electromagnetic energy from a continually changing field.
In a similar way, the complex permeability is defined by
µr(ω) = µ(ω) µ0 = µ0(ω) − jµ00(ω) µ0 = µ 0 r(ω) − jµ00r(ω) (2.10)
2.3
Waveguides
An overview of the different modes can be seen in the list below. In this thesis only the TEM mode will be used. • Ez = 0, Hz = 0, TEM modes • Ez = 0, Hz6= 0, TE or H modes • Ez6= 0, Hz = 0, TM or E modes • Ez6= 0, Hz6= 0, hybrid, HE or EH modes
Figure 2.1: Example of waveguide structures [11].
Electromagnetic waves can be described in a simple way by assuming a uniform plane wave propagating in the z-direction [11]. In a lossless material and with a fixed frequency the electric and magnetic fields are
E(z, t) = E(z)ejωt
H(z, t) = H(z)ejωt (2.11)
Deriving Maxwell’s equations with respect of the direction of propagation the electric and magnetic fields become
E(x, y, x, t) = E(x, y)ejωt−jβz
H(x, y, x, t) = H(x, y)ejωt−jβz
(2.12)
where β is the propagation wave number and are calculated by
β = 1 c0 p ω2− ω2 c = ω √ µ (2.13)
2.4
TEM waves
Transverse electromagnetic (TEM) mode is when both the electric and magnetic field lines are transverse to the direction of propagation. The electric component Ez= 0 and the magnetic component Hz= 0 are zero in z-direction. The TEM mode is derived from Maxwell’s equations and thorough derivations can be found in the book of Collin [3].
The electric and magnetic fields in the TEM mode are determined by equation 2.14 and 2.15 where the transverse electric and magnetic fields are perpendicular to each other and propagating in the direction of z.
∇T× ET= 0
∇T· ET = 0
(2.14)
HT = 1
Zˆz × ET (2.15)
Where Z is the characteristic impedance and is determined by the permittivity and permeability such as
Z =r µ
(2.16)
2.5
Coaxial line
Coaxial lines, can be seen in figure 2.2, are most widely used TEM transmission line and often used in a measurement set-up as shown in figure1.1. The line has two concentric conductors with an inner radius a and an outer radius b. The space between the conductors are filled with a dielectric material [11].
The characteristic impedance for the coaxial line can be calculated with equation2.17. Z = Z0 2πln b a (2.17) Where Z0=r µ0 r (2.18)
Chapter 3
Material Characterization
3.1
Estimation of material parameters
There are several different methods for measure and estimate the material parameters permittivity and permeability. Some of these methods are resonant methods, free space methods, reflection methods and transmission/reflection methods. Each method has limitations such as specific frequencies, applications and materials [18].
In this thesis the transmission/reflection method is used to determine the material parameters. Since using the TEM mode a suitable waveguide for the simulations is to use the coaxial cable to carry TEM waves in a given direction. As described in section 2.5 the transmitted and reflected fields are simulated using an infinite layer with thickness d. The principle is to measure the transmitted and reflected electromagnetic fields when a material is under illumination from an incident electromagnetic field [7]. Figure 3.1 shows a sample of a material under illumination of an incident electric field E1+ and the reflected E1− and the transmitted E2+ electric fields.
Scattering parameters are the amplitude from the measured signals. From the scattering parameters the per-mittivity and permeability can be obtained using different conversion methods. Two methods are Nicolson-Ross-Weir (NRW) method and Baker Jarvis (BJ) method. These will be applied in this thesis.
3.1.1
Scattering parameters
Scattering parameters (S-parameters) describes the behaviour of linear electrical networks. The parameters consists of either the reflection or the transmission coefficient of a component, normally a voltage, and are represented as complex numbers [12]. The scattering matrix describes the network and relates the ampli-tudes of the incident wave voltages to the reflected or transmitted ones. For a two port network, seen in figure3.2, the scattering matrix is described as3.1[12].
Figure 3.2 V− 1 V2− =S11 S12 S21 S22 V+ 1 V2+ (3.1)
Where V1,2+ are the amplitude of the incident voltage waves on the two ports and V1,2− are the amplitude of the reflected voltage waves. The reflected signals are S11, S22and the transmitted signals are S12, S21.
3.1.2
Nicolson-Ross-Weir Algorithm
The Nicolson-Ross-Weir (NRW) algorithm [8] is a fast, non iterative method based on reflection and trans-mission. NRW uses the scattering parameters S11 and S21 to determine the permittivity and permeability. S11 is the incoming wave divided with the reflected wave at port 1. S21 is the ratio of the transmitted wave at port 2 and the incoming wave at port 1 [19]. The equations for the NRW algorithm are as follows: The reflection coefficient is
Γ = χ ±pχ2− 1 (3.2) where χ = S 2 11− S212 + 1 2S11 (3.3)
T = S11+ S21− Γ
1 − (S11+ S21)Γ (3.4)
The complex permeability and permittivity are calculated according to
µr= 1 Λqλ12 0 − 1 λ2 c 1 + Γ 1 − Γ (3.5) r= λ2 0 µr 1 λ2 c − 1 Λ2 (3.6)
where λ0is the free space wavelength in vacuum and λcis the cutoff wavelength. The normalized wavelength Λ is defined by 1 Λ2 = − 1 2πlln( 1 T) 2 (3.7)
where l is the length of the sample.
The NRW method has an instability since equation3.7has an infinite number of roots when the frequency corresponds to a multiple of one-half wavelength in the sample. The phase ambiguity can be corrected by
ln(T ) = ln(T ) + j(2πn), n = 0, ±1, ±2, ... (3.8)
The advantages are that the NRW method is fast and non-iterative. The disadvantages are that the sample of the material should be short and it is not suitable for low loss materials [18].
3.1.3
Baker Jarvis Method
Baker Jarvis (BJ) method is an iterative algorithm where all the scattering parameters are taken into account compared to the Nicolson-Ross-Weir method where only S11 and S21 are used. The iterative BJ method is described by equation3.9.
1
2[(S12+ S21) + β(S11+ S22)] =
T (1 − Γ2) + βΓ(1 − T2)
1 − T2Γ2 (3.9)
and l as the length of the sample. Γ is the reflection coefficient defined by Γ = γ0 µ0 − γ µ γ0 µ0 + γ µ (3.10)
where γ0, γ are the propagation constants for air and for the material and µ0, µ are the permeability of vacuum and the permeability of the material.
Chapter 4
Simulation tools
4.1
Finite-Difference Time-Domain method
Finite-Difference Time-Domain method (FDTD) solves Maxwell’s equations numerically in the time-domain, where the partial derivatives are replaced by finite-differences. The version of FDTD solver used in this thesis has been developed at FOI (Swedish Defence Research Agency) as part of different research projects [9] and since 2014 shared with Saab as part of an agreement between Saab and FOI. The version used has periodic boundary conditions in x- and y-directions for simulations of infinite periodic structures.
The details about the FDTD-algorithm can be found in [9]. The electromagnetic fields are updated in the time domain at discrete time steps. The scattered fields propagate to an outer boundary condition – UPML (Uniaxial Perfect Boundary Condition). This type of boundary condition absorbs the scattered waves more or less independently of the incidence directions.
A typical length of a FDTD-simulation is in the order of approximately 1000-30 000 time steps, depending on the size and properties of the scattering object. Resonant structures require longer runs and this is because the energy in the system must decay sufficiently before Fourier transformation of the results into the frequency domain can be performed. The simulations performed in this thesis are either 12 000 or 100 000 time steps.
Figure 4.1: Visualization of the Yee-cell [9].
4.2
Simulation set-up
Chapter 5
Simulations
There are several ways to model and generate particles depending on material, shape, distribution and how many particles that will form an aggregate. Since the material to be modeled is only known from photos
1.2, 5.1, the modeling of particles are made in many steps in order to rule out different methods and make new approaches in order to get closer to the reality.
When the particles in figure 5.1 are embedded in a matrix material, experiments show that there is little to no percolation. Percolation occurs when the particles are connected to each other and creates a path for the current to travel. The idea is to model particle shapes similar to the ones in the figure and study the spread of the relative permittivity with controlled content of particles. Since the relative permittivity are expressed as a complex value the occurrence of percolation can be found while studying the imaginary part of the relative permittivity.
5.1
Simulation specifications
The parameters used in the simulations can be found in the list below. The relative permittivity and con-ductivity for the particles and the background material are fixed parameters given from the company Saab. The shape and size of the particles will vary within the range in the different simulations.
• Relative permittivity of background material: 3 • Relative permittivity of particles: 5
• Conductivity of background material: 0 S/m • Conductivity of particles: 104 S/m
• Permeability: 4π × 10−7 H/m • Frequency: 0-18 GHz
• Cell dimension in FDTD volume: 10−6m
• Dimension of small FDTD volume: 61 × 61 × 201µm • Dimension of large FDTD volume: 121 × 121 × 201µm • Range of particle dimension: 4-50 µm
• Content of particles: 10%, 15% and 20%
5.2
Cluster-Cluster aggregation model
The first approach was to simulate clusters of aggregates with a predefined fractal dimension. The fractal dimension can be calculated using the fractal scaling equation which is described in appendixA. This model is based on the description by Thouy and Jullien [14] and is also well described in a previous work [13]. Therefore only the principle will be described in this section. A Matlab-script was provided with this model and the aggregate consisted of spherical particles. The script were modified to represent the aggregate with cubes since the FDTD volume consists of a three-dimensional grid. This method require the input of a specific number of particles (np) and a fractal dimension (Df ∈ [1,3]) [13].
The principle of the cluster-cluster (CC) model is to merge two clusters of particles together where the starting point is to create a cluster of two particles. The next step is to merge the previous cluster with a new one consisting two particles. The merging continues until the number of particles chosen by the user are merged together and results in one final cluster.
Figure 5.2: Visualization of a particle cluster consisting of 32 particles and with a fractal dimension of 1.5.
5.3
Particle model based on ellipsoids
Since the shape of the particle aggregate in the CC model was to different from the ones in figure5.1, the new approach became to model particles based on the shape of an ellipsoid.
5.3.1
Modeling a single particle
The particle is modeled using the equation of an ellipsoid which is described in appendix B. A three-dimensional grid is applied trough the ellipsoid to be able to extract the coordinates of the grid located within or on the surface of the ellipsoid. By extracting the coordinates one can represent the ellipsoid with cubes. The criteria to extract the coordinates needs to be less or equal to one considering the ellipsoid equation.
The advantage of using this method is that the dimension of the particle can be chosen and changed in an easy way to imitate the shape of an already known particle. An example of a particle are shown in figure
5.4and5.5where the radius in x-, y- and z- direction is 14, 8 and 4 respectively.
Figure 5.5: Visualization of a particle represented by cubes.
5.3.2
Modeling of several particles
Each particle is modeled with its center at (0,0,0) in the coordinate system. When generating several particles the previous one needs to be moved to a new point inside the FDTD volume to prevent the next generated particles to gather at point (0,0,0). The translation of a particle is made using a vector with random x-, y-and z- coordinates. This vector is added to the coordinates of the particle to move it from the origin. The development to model several particles consisted of two steps to gradually design a functional model. The steps are:
• Identical shaped particles.
• Identical shaped particles rotated in space.
The first approach with identical shaped particles was mainly to prove that the particles become distributed in the FDTD volume. The approach where the particles is rotated were applied to not have identical particles to model a realistic material.
When the first particle is modeled it is rotated in the xyz-plane using rotation matrices described in appendix
C. The rotation around each axis occurs randomly with an angle between 0◦ and 360◦. The rotated particle is moved to a random point inside the FDTD volume to prevent the next modeled particle to be at the same coordinates as the previous ones.
There are two criteria for the generated particle that must be fulfilled:
2: If the particle moves across the boundary of the FDTD volume in ± z-direction, the excess part of the particle will be removed.
Figure5.6and5.7shows an example of when the above criteria are fulfilled. One can observe the periodically movement at the boundary to the right and left in the figures. At the top of the figures one can observe that the excess part of the particles have been removed. The matrix material prevents observations on how the particles are distributed inside the FDTD volume and will be removed in the future, but only for future visualizations.
Figure 5.7: Visualization of the boundary criteria for the particles.
5.3
Initial Simulations
5.3.1
Shape analysis
The shape analysis consists of four cases with particles shaped as spheres, needles and discs with two different sizes. Each case are simulated 25 times where the volume content of the particles is approximately 10%. With regard of the boundary conditions and the randomly distributed particles in the small FDTD volume, the volume content can differ for the same amount of particles in each simulation.
•
Case 1: Spheres
The particles have the shape of spheres with 6µm radii. The number of particles are 84 and the range of the volume content is 10 ± 0.51%
•
Case 2: Discs
The radii of the disks in x- and y-direction is set to 12µm and in z-direction 3µm. The number of particles are 43 and the range of the volume content is 10 ± 0.66%.
•
Case 3: Needles
The needle shaped particles have 3µm radii in x- and y-direction and the radius in z-direction is 12µm. The volume content has a range of 10 ± 0.25% for a number of 169 particles.
•
Case 4: Larger discs
The radii in x- and y-direction is 24µm and in z-direction the radius is 3µm. This type of disk will be larger and thinner compared to case 2. The number of particles are 11 with a volume content range of 10 ± 0.64%.
Figure 5.8: Visualization of particles shaped as spheres (left) and as discs (right).
Each case are simulated 25 times and the relative permittivity can be represented over the frequency range 0-18 GHz. The simulations are made in time-domain an the results are transformed with the Fourier trans-form to the frequency-domain. The transmitted electric field is represented by the scattering parameter S21
3.1which is used in Baker-Jarvis method 3.9to calculate the relative permittivity. The relative permittiv-ity is shown in figure 5.10 and the figure contains the complex permittivity estimated using Baker-Jarvis (BJ) method and also Maxwell-Garnet (MG). The homogenization is an approximation of a heterogeneous structure using a mixture formula [6]. The relative permittivity for the matrix material is 3, and 5 for the particles. The value of the permittivity given by the homogenized version should be a value between the ones for the materials. The value for the relative permittivity estimated by Maxwell-garnet remain the same for all 25 simulations as long as the particle content is the same. Further explanations of the mixture formula can be seen in [6].
Figure 5.10: The complex relative permittivity for simulation number 2 from the case with needles. Further visualizations of the relative permittivity will be in form of histograms to also get an overview of the distribution. Each bin in the histogram represent the mean value of the real part of epsilon for each simulation such as the one shown in figure5.10.
Figure5.11shows the different distributions of epsilon for the different shaped particles. The mean value of both the real part and the imaginary part of epsilon can be seen in table5.1. The mean values are calculated based on the values represented in the histograms. It can be observed that the real part of epsilon gets a higher value when the particles are small and shaped as needles. The reason epsilon gets a higher value is due to higher charge accumulation at particle tips and edges. If several particles are in connection with each other, then longer dipoles are formed which increases the real part of the permittivity.
Table 5.1: Mean values of the real part of epsilon for the cases in the shape analysis.
Case 1 2 3 4
Mean value of Re() 4.4876 5.8025 9.6179 4.9190
The red bins in the histogram indicates that there is percolation in the simulations. Equation 2.9 shows that the imaginary part of the relative permittivity is dependent on the frequency. Here the 1/f dependence indicates high homogenized conductivity for low frequencies and the mean value entails a high value of the imaginary part. The limit for the mean value of the imaginary part is set to 10 to indicate percolation of the particles when the mean value is higher. The percolation can become very high and occurs when the particles are connected to each other and creates a path for the current to travel. It indicates how high the electrical losses are. Figure 5.12shows an example of the percolation from the case with needle shaped particles. The left graph shows simulation number 4 with percolation and the right graph shows simulation number 2 without percolation. Observe that the graphs have a logarithmic scale on the y-axis to be able to also show the real part of epsilon with a value of approximately 10 and 9 respectively.
Figure 5.12: Example of how the imaginary part of epsilon looks for simulations with percolation (left) and without percolation (right).
5.3.2
Distribution analysis
The distribution analysis starts with two cases where the size and shape of the particles varies. In the first case the shapes from the shape analysis are distributed in which each shape contributes to one quarter of the 10% volume content of particles. The second case distributes randomly shaped particles. As in the previous analysis each case are simulated 25 times with the particles volume content approximately at 10%.
•
Case 1: One quarter of each shape
The spheres, needles and discs are identical to case 1-4 in the shape analysis. The distribution is 20 spheres, 11 discs, 40 needles and 3 larger discs. The range of the volume content is 10 ±0.55%
•
Case 2: Randomly shaped particles
The shape of each particle is random. The radii in x- and y-direction are in the range of 12 ± 6µm and the radius in z-direction is in the range of 6 ± 3µm. With this conditions the particles can vary in shape from thin discs to long needles. With 22 particles the range of the volume content is 10 ±0.79%.
Figure 5.13: Visualization of randomly shaped particles.
The real part of the relative permittivity for the two cases can be seen in figure5.14. By represent the values of epsilon in histograms one can easily observe the distribution. For the first case the distribution is quite large. The red bins in the histogram corresponds to the simulations where percolation occur. Compared to the second case where there is no percolation, the distribution is more narrow and epsilon is also lower. The values of epsilon can be seen in table5.2. The reason for the higher value of epsilon in case 1 is due to the amount of particles shaped as needles that gives a higher real part of epsilon as observed in the shape analysis.
Figure 5.14: Distribution of the real part of epsilon for the cases in the distribution analysis.
Table 5.2: Mean values of the real part of epsilon for the cases in the distribution analysis.
Case 1 2
It was decided to continue to simulate randomly shaped particles and have an other range for the radii of the particles to get more variation of the sizes and shapes. The new dimensions of the radii is set to 12 ± 10µm in x- and y-direction and 3 ± 1µm in z-direction. With the new dimensions, two new cases are developed and simulated 25 times each, where the first case has a content of particles of 10% and the second has a content of 15%.
•
Case 3: 10% randomly shaped particles
The shape of each particle is random and 40 particles is generated to fill the FDTD volume to approx-imately 10%. The range of the particle content varies with 10 ± 0.96 %.
•
Case 4: 15% randomly shaped particles
The shape of each particle is random and 65 particles is generated to fill the FDTD volume to approx-imately 15%. The content of particles varies in a range of 10 ± 0.92 %.
Figure 5.15: Distribution of the real part of epsilon in the distribution analysis with different particle content.
Table 5.3: Mean values of the real part of epsilon for the new cases in the distribution analysis.
Case 3 4
Mean value of Re() 6.6322 10.2110 Mean value of Im() 13.4154 119.6729
5.4
Main Simulations
The main simulations contains several different approaches. As observed in the initial simulations, the perco-lation occurs in several cases. Since the amount of percoperco-lation is not expected due to previous measurements where the occurrence is quite low, the following simulations are also made to study how the percolation can be prevented. One idea to decrease the percolation is to create a non conducting layer around the particles to prevent any connection between to avoid creating a conducting path for the current. The layer is set to have the same permittivity as the background material and a thickness of 2 µm. The simulations consists of particles with and without a non conducting layer for comparison. An other approach to avoid the percolation is to expand the FDTD volume in x- and y-direction to get the particles more distributed to avoid larger connections between the particles. Figure5.16visualizes the small and the large FDTD volumes.
Figure 5.16: Visualization of the small FDTD volume (left) and the large FDTD volume (right). The drawback with the non conducting layer is the forced break between each one of the particles. The reality shows that the particles have some amount of connection and the simulations are not fully trustwor-thy. An other approach will be to have conductivity through the whole particle and instead introduce a gap through the FDTD volume. The gap extends through both yz- and xz-direction in a few simulations and only in yz-direction on other. The thickness of the gap is either set to 2 µm or 5 µm and it has the same properties as the background material. This is reasonable due to the periodicity of the boundary conditions and this allows some particles to create a path for the current but not create it periodically. Figure 5.17
visualizes how the non conducting layer looks like together with a visualization of the gap.
Figure 5.17: Visualization of particles with a non conducting layer (left) and the FDTD volume with a gap extending in yz-direction (right).
5.4.1
Particle content of 10%
At a particle content of 10% all of the aforementioned ideas are tested. Each case are simulated 128 times and the range for the content of particles is within 10±0.5%. This built the following cases:
•
Case 1: The small FDTD volume, particles with non conducting layer
The particles have a non conducting layer with a thickness of 2µm simulated with the small FDTD volume. The number of particles 50.
•
Case 2: The small FDTD volume, particles without non conducting layer
The particles does not have the non conducting layer and the number of particles are 45, simulated with the small FDTD volume.
•
Case 3: The large FDTD volume, particles with non conducting layer
The particles have a non conducting layer with thickness of 2µm simulated in the large FDTD volume. The number of particles are 200.
•
Case 4: The large FDTD volume, particles without non conducting layer
•
Case 5: The large FDTD volume, one gap
The number of particles are 180 and are without the non conducting layer. The simulation is with the large FDTD volume. In this case a gap with a thickness of 2 µm is introduced in yz-direction to try to prevent the percolation since the boundary condition is periodic.
Table 5.4: Mean value of the real part of epsilon for different cases with a particle content of 10%.
Case 1 2 3 4 5
Mean value of Re() 5.2037 6.6291 5.2042 7.7139 6.9019
Mean value of Im() 0.0255 28.8102 0.0232 2.3563 0.5847
The first two graphs shows case 1 and case 2, the second two shows cases 3 and 4. Case 5 is shown in the last graph5.18. There can be observed that the distribution is small when having the non conducting layer on the particles with the large FDTD volume corresponding to the smallest distribution of epsilon with value 1. It can be noted that the value of epsilon also is lower when having the layer to prevent the connection between the particles. Simulating without the non conducting layer gives a very large distribution that extends to approximately 7.
5.4.2
Particle content of 15%
For the 15% particle content there is a change in the cases from the previous section. The small FDTD volume is only used when simulating particles with the non conducting layer since there was observed al-ready in the previous section that percolation occurred when using the small FDTD volume and particles without the layer. Since there were so few simulations with percolation present for the large FDTD vol-ume in the previous section it was decided to simulate a case with a particle content of 15% with particles without the layer. This case will also be simulated with more time steps to increase the time of the simu-lation to see if there is any difference in the result. Here the range of the particle content are within 15±0.5%.
•
Case 1: The small FDTD volume, particles with non conducting layer
The small FDTD volume is used and the particles have the non conducting layer with a thickness of 2 µm. The number of particles are 60.
•
Case 2: The large FDTD volume, particles with non conducting layer
The particles have the non conducting layer of 2 µm and are simulated with the large FDTD volume. The number of particles are 240.
•
Case 3: The large FDTD volume, one gap
The number of particles are 270 and they are without the non conducting layer. The simulations are with the large FDTD volume. In this case there are one gap through the volume in yz-direction to prevent a periodic connection between the particles.
•
Case 4: The large FDTD volume, 12 000 time steps
The particles does not have the non conducting layer and the large FDTD volume is used. The number of particles are 270 and the time steps in the simulations are set to 12 000.
•
Case 5: The large FDTD volume, 100 000 time steps
This case is identical to case 4 except that the time steps for the simulations are increased to 100 000.
Table 5.5: Mean value of the real part of epsilon for different cases with a particle content of 15%.
Case 1 2 3 4 5
Mean value of Re() 7.0351 6.9679 12.9142 14.8331 18.5268 Mean value of Im() 0.4801 0.3631 10.4513 57.9200 47.8740
Figure 5.19: Distribution of the real part of epsilon for particles with a non conducting layer for the small and large FDTD volume together with a case of distribution when introducing a gap.
Introducing a gap in the third case gives interesting results. The mean value of epsilon increases sharply to almost 13 and have doubled in value for a particle content of 15% compared to the same case in previous section. By not having any other restrictions the distribution for this case become very high and the red bins that indicating percolation is near 50%.
Measurements performed at Saab does not show the amount of percolation the simulations does. A reason-able explanation is that the transmitted electric field in the simulations has not decayed sufficiently. Figure
5.21shows the transmitted electric field for a case that indicates no percolation and a case with percolation. In the first graph one can see that the field is back to zero and have completely decayed. The second graph shows that the time steps are to few for the electric field to decay sufficiently and be back to zero. This means that the red bins in the histograms in this section can not fully be explained to be percolation.
Figure 5.21: Graphs showing the transmitted electric field sufficiently decayed (left) and not sufficiently decayed (right).
5.4.3
Particle content of 20%
For the particle content of 20% the cases are the same as in the previous section with a particle content of 15%. Only one new case is introduced where a second gap is added through xz-direction to prevent percolation. The range of particle content are within 20±0.5%.
•
Case 1: The small FDTD volume, particles with non conducting layer
The number of particles are 90 and they have the non conducting layer with a thickness of 2 µm. For the simulation the small FDTD volume is used.
•
Case 2: The large FDTD volume, particles with non conducting layer
The number of particles are 630 and they have the 2 µm thick non conducting layer. The large FDTD volume is used in this simulation.
•
Case 3: The large FDTD volume, one gap
The number of particles are 360 and the simulations are made with the large FDTD volume. In this case there are one gap through the volume in yz-direction to prevent a periodic connection between the particles since they do not have the non conducting layer.
•
Case 4: The large FDTD volume, two gaps
The number of particles are 360 and the simulations are made with the large FDTD volume. In this case a second gap is introduced through the volume in xz-direction. The two gaps form a cross to prevent a periodic connection between the particles since they do not have the non conducting layer.
•
Case 5: The large FDTD volume, 12 000 timesteps
The particles does not have the non conducting layer and the large FDTD volume is used. The number of particles are 360 and the time steps in the simulations are set to 12 000.
•
Case 6: The large FDTD volume, 100 000 timesteps
This case is identical to case 5 except that the time steps for the simulations are increased to 100 000.
The mean values of the relative permittivity can be seen in table 5.6. As seen in the previous sections the real part of epsilon gets a higher value when the contents of particles increases. For a particle content of 20% epsilon has indeed increased for a higher particle content of 20 %. The first two graphs in figure5.22
corresponds to case 1 and 2 where the particles have the non conducting layer. The second graph simulated with the large FDTD volume shows that the distribution is fairly small and extends to approximately 1.5 compared to the first graph simulated with the small FDTD volume where the distribution is quite large.
Table 5.6: Mean value of the real part of epsilon for different cases with a particle content of 20%.
Case 1 2 3 4 5 6
Mean value of Re() 9.3372 8.4110 22.6558 20.7778 16.1896 17.2073
Mean value of Im() 1.9555 0.2099 48.2153 37.7851 414.3614 421.4701
indicates percolation throughout all simulations. Similar to the corresponding cases in previous section it can not be established that there are in fact percolation since the electric field has not decayed sufficiently for all simulations.
Figure 5.22: Distribution of the real part of epsilon for particles with a non conducting layer for the small and large FDTD volume together with the cases containing one and two gaps respectively.
Chapter 6
Conclusions
Increasing the particle content was challenging since the probability for percolation was higher than expected. The probability of percolation depends, among other things, on periodic structures. Since there are periodic boundary conditions, it can form long connections between the particles which increases the dipole moment resulting in percolation. It is possible to increase the particle content using a larger FDTD volume in the simulations to prevent percolation. It was also discovered that truncated simulation time can be perceived as percolation and the time steps needs to be set carefully.
The course of action to obtain realistic particles was challenging since no previous models could be used. The modeling of particles using ellipsoids does not reflect the existing material completely but is similar enough to get useful results. It can be confirmed that the shapes of the particles affects the relative permittivity which increases for particles with sharper edges.
The simulations shows that for a given content of particles, the spread of the relative permittivity is relatively large. When the particle content gets higher, the spread of the relative permittivity increases.
6.1
Future work
Appendix A
The pattern of fractals have self-similarity and scaling the fractals will repeat the exact same pattern over and over. The connection between the self-similarity and the fractal dimension is described by the fractal scaling equation [13]:
NL∼ L−Df (A.1)
Where NL is the number of sticks with length L needed to measure the length of an object with fractal dimension Df. The fractal dimension can be calculated using the combination of equations below.
NL∼ L−Df
log NL∼ −Dflog L
Df = −∆ log NL ∆ log L
Appendix B
The shape of an ellipsoid is described by equationB.1when using the Cartesian coordinate system with the origin in the centre of the ellipsoid [17].
(X − Xc)2 Xr2 + (Y − Y c)2 Y r2 + (Z − Zc)2 Zr2 = 1 (B.1)
Appendix C
Elemental rotation of a vector is a rotation around one of the axes and can be applied in several dimensions. Rotating a vector in three dimensions can be described by the equations C.1,C.2 andC.3. The rotations occur around the x-, y- and z-axes with angles α, β and γ respectively [1].
Rx(α) = 1 0 0 0 cos α − sin α 0 sin α cos α (C.1)
Counterclockwise rotation around the x-axis with an angle α.
Ry(β) = cos β 0 sin β 0 1 0 − sin β 0 cos β (C.2)
Counterclockwise rotation around the y-axis with an angle β.
Rz(γ) = cos γ − sin γ 0 sin γ cos γ 0 0 0 1 (C.3)
References
[1] H. Anton, C. Rorres Elementary Linear Algebra, 11th ed. Wiley, 2014.
[2] C. A. Balanis Advanced Engineering Electromagnetics, John Wiley & Sons, 1989. [3] R. E. Collin Foundations of Microwave Engineering, 2nd ed. John Wiley & Sons, 2001
[4] J. Baker-Jarvis, M. D. Janezic, B. F Riddle, R. T. Johnk, P. Kabos, C. L. Holloway, R. G. Geyer, C. A. Grosvenor Measuring the Permittivity and Permeability of Lossy Materials: Solids, Liquids, Metals, Building Materials and Negative-Index Materials, 2005, NIST Technical Note 1536.
[5] J. D. Jackson Classical Electrodynamics, 3rd ed. John Wiley & Sons, 1962.
[6] G. Kristensson, C. Larsson, T. Martin, D. Sjöberg, N. Wellander Mikromodellering av radarabsorbenter, II, 297769-LB784962, 2010-03-31.
[7] E. Lindenholst, F. Olsson TEM Cell design for Material Characterization, Lund University, 2016. [8] A. M. Nicolson, G. F. Ross Measurement of the Intrinsic Properties of Materials by Time-Domain
Tech-niques, IEEE Transactions on Instrument and Measurement, 19(4):377-382, Nov 1970.
[9] T. Martin Broadband Electromagnetic Scattering and Shielding Analysis using the Finite Difference Time Domain Method, PhD thesis No. 669, Linköping University, Linköping, 2001.
[10] K. A. Milton, J. Schwinger Electromagnetic Radiation: Variational Methods, Waveguides and Acceler-ators., Springer, August 2005.
[11] S. J. Orfanidis Electromagnetic Waves and Antennas, www.ece.rutgers.edu/ orfanidi/ewa, 2002-2016. [12] D. M. Pozar Microwave Engineering, 4th ed. John Wiley & Sons, 2012.
[13] A. Roth Image analysis to estimate the fractal dimension of soot aggregates, Lund University, 2018. [14] R. Thouy, R. Jullien A cluster-cluster aggregation model with tunable fractal dimension, Journal of
Physics A (Mathematical and General), 27:2953-2963, 1994.
[15] A. Vicente, G. Dip, C. Junqueira The step by step development of NRW method, 2011, 10.1109 / IMOC.2011.6169318.
[16] A. A. Ward Dielectric materials for advanced applications, 2016, 10.13140/RG.2.1.3481.5600.
[17] The MathWorks, Inc., Ellipsoid, 1994-2019, https://se.mathworks.com/help/matlab/ref/ellipsoid.html [18] Application Note Measurement of Dielectric Material Properties, Rohde & Schwarz,
RAC-0607-0019_1_5E
