Mathematical Modelling of EM Fields

The 3D Finite Difference Time Domain (FDTD) method involves simulating the propagation of an Electromagnetic (EM) wave over a series of time steps. Each of these time steps are extremely small, typically just a few pico seconds. When the wave has propagated then the frequency information can be extracted using a Fast Fourier Transform FFT).

Much of my work has involved the application of parallel processing to the 3D FDTD. The following is an introduction to the FDTD method.

1. Introduction

The processing power and memory capacity of modern computers increase by the year. This has made possible the simulation of electromagnetic field problems in the time-domain rather than in the frequency-domain. Another change in simulation techniques has been from continuous equations to discrete approximations. These discrete forms are usually easier to implement on a computer.

A good example of a time-domain simulation, using discrete equations, is the 3D FDTD method. It determines the frequency response over a wide spectrum of frequencies, whereas many other simulation methods require different models and/or techniques for different frequency spectra. Papers [3.1]-[3.12] outline the basic theory and application of the 3D FDTD method.

The 3D FDTD method derives directly from Maxwell's curl equations and is relatively simple to implement. Unfortunately, it requires large amounts of computer memory and processing time and, has, in the past, only been used with super-computers which have the processing power and memory capacity to apply it. With the arrival of high-speed desktop computers with large and cheap memory storage the method can now be fully exploited in the areas such as microstrip antenna modelling, analysis of microstrip circuits and in biological applications.

The 3D-FDTD method has two main advantages over empirical analysis. It provides a direct solution to Maxwell's equations without much complexity and takes into account both the electric and magnetic fields in a 3D model.

As the 3D FDTD method is time-based the results produced can also help to provide an understanding of EM wave propagation within the structures. Frequency-domain techniques often conceal how the EM waves propagate within the structure. For example, a microstrip antenna (or patch antenna) can be modelled as a transmission line, or as lumped parameters. This modelling can often hide the fact that the incident waves within the antennas head are reflected back and forward within the antenna. Electromagnetic radiation then leaks out from the ends of the patch. Results from 3D FDTD simulations allow the wave to be visualised, which helps in checking results.

This chapter discusses the theory of the 3D FDTD method and its application to electromagnetic wave propagation within microstrip antennas and PCBs. Chapters 8 and 9 apply the 3D FDTD method and show some 3D pulse visualisations.

2. Background

Yee [3.1] was the first researcher to propose a differencing scheme which is now known as the FDTD method. It has since been used to model microstrip circuits [3.2]-[3.7], to scattering problems [3.8] and in the simulation of electromagnetic radiation [3.9], [3.10]. Other researchers have applied it to the simulation of waveguides, to coaxial cable and similar structures [3.11]-[3.16], and to digital signal processing and ferrites [3.17]-[3.20]. It is also useful in areas, such as in Biomedical research to model electromagnetic radiation on human tissues and to radar wave simulations.

3. Simulation Steps

Figure 3.1 shows the main steps taken in a 3D FDTD simulation. Initially, a 3D model is made to represent the physical structure, including conductors, dielectics and boundaries. Next an applied pulse, normally either a sine-wave or a Gaussian pulse, acts as the input stimulus at all the sources. Then at increments of time the E and H fields are calculated. After each increment the input electric field amplitude is calculated and the E and H fields are again recalculated. This continues until the E and H fields within the system decay to zero.

After completing the simulation an FFT program extracts frequency information from the transient response. The location of the transient data depends on the required system response. For example to determine the reflection coefficient, the input and reflected waves at the sources are monitored. For a radiation pattern, points are taken in free-space around the structure.

 



Figure 3.1: FDTD method

4. Finite-Difference Time-Domain (FDTD) Method

The FDTD method uses Maxwell's equations which define the propagation of an electromagnetic wave and the relationship between the electric and magnetic fields, these are:

 

(3.1)

(3.2)

(3.3)

(3.4)

 

For a uniform, isotropic and homogeneous media with no conduction current Maxwell's curl equations then become:

(3.5)

(3.6)

 

By applying appropriate boundary conditions on sources, conductors and mesh walls an approximate solution of these equations can be found over a finite three-dimensional domain. Taking an example of the first equation in the i direction gives:

 

(3.7)

 

The central difference approximation is then used on both the time and space first-order partial differentiations to obtain discrete approximations. This gives:

 

(3.8)

 

Rearranging gives:

 

(3.9)

 

The half time-steps indicate that E and H are calculated alternately to obtain central differences for the time derivatives. In total there are six equations similar to Equation (3.9). These define the E and H fields in the x, y and z directions and are given in Equations (3.10) and (3.11). The permitivity (e ) and the permeability (m ) values in these equations are set to approximate values depending on the location of each of the field component.

 

(3.10)

(3.11)

5. Problem Conception

The structure simulated in Chapter 9 is a PCB with four electrical sources, as shown in Figure 3.2. It consists of a substrate layer, such as Duroid (relative permittivity, e_r, of 2.2) above a ground plane. A copper layer is formed by etching the top of the substrate to give the required pattern.

1. 3D gridding

A 3D grid is placed around the structure, as illustrated in Figure 3.3. The number of cells within the grid is normally selected with consideration to the simulation time limit and the amount of computer memory. As an example a linear grid placed around a microstrip antenna contained within a volume of 30´ 30´ 9.6 mm³ with a 100´ 100´ 12 grid gives a element volume of 0.3´ 0.3´ 0.8 mm³.



Figure 3.2: PCB with copper tracks

The first grid point in the z-direction lies on the top of the ground plane. Normally, there are fewer cells in the z-direction because there are very few discontinuities in this direction. A discontinuity causes reflections in the electromagnetic wave and they have a great effect on the frequency characteristics of the simulated model. Thus, to provide higher accuracy around discontinuities, a non-linear grid is sometimes placed around them. A fine grid is placed around discontinuities and a course grid where there are no discontinuities.

2. Permittivity and permeability

The calculation of the magnetic fields involves permeability. As conductors are assumed to have zero thickness, the value of m _r is always taken as 1 (thus the permeability m will be 4p ´ 10^-7 H m^-1).

The calculation of electric fields uses permittivity which varies depending on whether the field is within the substrate or in the surrounding air. The permittivity in the medium above the substrate is e _{r1e 0}, and within the substrate it is e _{r2e 0}, (where e ₀ is 8.854´ 10^-12 F.m^-1). At the interface between the air and the substrate, the approximate relative permittivity is taken to be the average of the two, that is:

 

(3.12)



Figure 3.3: 3D gridding

3. Input signal

The input signal can be of any shape, but, it is normally a Gaussian pulse. This type of pulse has a frequency spectrum that is also Gaussian and thus has the advantage of providing frequency information from DC up to a desired cut-off frequency. The form of the input signal in a continuous form is:

(3.13)

 

where t₀ is the pulse delay and T relates to the pulse width. Written in a discrete form gives:

(3.14)

 

where n is the current time-step, m the pulse delay time-step and x the width of the pulse in time-steps. Figure 3.4 shows Gaussian pulses with pulse widths of 5, 10 and 20 time-steps. Each pulse has been delayed by 30 time steps.

The width of the Gaussian pulse sets the required cut-off frequency. Figure 3.5 shows the relative power of a Gaussian pulse width of 5, 10 and 15 time steps. In can be seen that the thinner the pulse the larger its signal bandwidth.

The pulse width is normally chosen to have at least 20 points per wavelength at the highest frequency significantly represented in the pulse. In most cases in this thesis the pulse width is 11 time-steps, which gives a bandwidth of 20 GHz.

Initially in the simulation, all the electric and magnetic fields are set to zero. The Gaussian pulse applied at the source has only a field component which is perpendicular to the ground plane (that is, E_z). Thus, E_y and E_x, at the source, are always zero. A change in the electric field at the source with respect to time causes a change in the magnetic field in the x-direction. Thus, the wave propagates in the y-direction, as shown in Figure 3.6.



Figure 3.4: Gaussian pulse



Figure 3.5: Gaussian pulse



Figure 3.6: Propagation of the wave

4. Conductors

The 3D FDTD method assumes perfect electrical conductors. Thus, the tangential electric field components that lie on the conductors are assumed to be zero. Figure 3.7 shows that the E field components on the conductor will be zero in the x- and y- direction.



Figure 3.7: Conductor treatment

5. Boundary walls

There is a limit to the size of physical grid applied around the model. To reduce the requirements for a large grid an absorbing wall is placed on the six mesh boundary walls. The electric field component tangential to the ground plane is zero and the tangential electric fields on the other five mesh walls are calculated so that a wave propagating against them does not reflect back. For the structure simulated in this thesis the pulses are normally incident on the mesh walls. This leads to simple approximations for continuous absorbing boundary conditions. The tangential fields on the absorbing boundaries then obey the one-dimensional wave equation in the direction normal to the mesh wall. For the normal y-direction wall the one-dimensional wave equation may be written as:

 

(3.15)

 

This equation is Mur's [3.21] first approximate absorbing boundary condition and in a discrete form it is:

 

(3.16)

 

where E₀ represents the tangential electric field on the mesh wall and E₁ the electric field one grid point within the mesh wall. Similar equations can also be derived for the other four absorbing mesh walls.

The method, unfortunately, does not take into account fringing fields that are propagating tangential to the walls. Thus the absorbing boundary must be placed well away from any fringing fields.

6. Maximum time step

The maximum time step that may be used is limited by the stability restriction of the finite difference equations [3.2]. This is given by:

 

(3.17)

 

where c is the speed of light (300 000 000 m.s^-1) and D x, D y and D z are the dimensions of the unit element. Table 3.1 gives two example time steps for different element sizes.

Table 3.1: Example time intervals

Model size (mm3)	nx, ny, nz elements	D x (mm)	D y (mm)	D z (mm)	D t (picoseconds)
30´ 10´ 10	100, 100, 12	0.3	0.3	0.3	0.68
80´ 80´ 50	100, 100, 10	0.8	0.8	0.5	1.25

6. Extracting Frequency Data

A fourier transform extracts frequency information from the transient response. Figure 3.8 shows an example of the electric field at a source. It can be seen that initially the Gaussian pulse is applied at the input. Then, after the pulse reaches the head of the antenna, a pulse reflects back to the source from the interface between the antenna head. Reflected pulses also return back from the resonance with the antenna head. The reflection coefficient is then the ratio of the reflected wave divided by the applied wave. In general, the scattering parameters S_jk may be obtained using a fourier transform on the transient waveforms, thus:

 

(3.18)

 



Figure 3.8: Applied and reflected wave

7. Improvements to the FDTD Method

Improvements can be made to the FDTD method which can improve accuracy, such as the sub-gridding method around discontinuities [3.22] and a modified frequency domain Finite-Difference Method that condenses nodes and uses an gifs/Image principle [3.23] can be devised. Other researchers have incorporation of static field solutions and Z-transforms into the FDTD method [3.24]-[3.26].

Improvements can also be made to the boundary conditions and the modelling of sources [3.27]-[3.29]. Simulation time can be reduced by using parallel processing methods [3.30],[3.31]. These will be discussed in more detail in chapters 6 and 7.

8. References