1 Introduction
The Finite-Difference Time-Domain (FDTD) method is a computational electromagnetic technique for solving for the electric and magnetic fields in arbitrary spatial domains in the time
domain. In contrast to techniques such as the Finite Element Method (FEM) and the Method of
Moments (MoM), this technique is straightforward to understand and is simple to program. A
rudimentary 2D TMz code is included in Section §6 and is used to illustrate the main features
of the method.
The aim of this assignment is to use the provided FDTD code in a series of numerical investigations, and compare quantitatively its predictions against theory, which the student is expected
to research independently after the completion of the taught part of the EMAP module. A
formal report is not required, but your assignment report needs to answer all the assignment
questions, in a self-contained manner.
2 The basics behind the FDTD algorithm
2.1 Defining the Lattice
The basic FDTD method (in Cartesian coordinates) makes use of a regular lattice of interleaved
electric and magnetic field components as originally proposed by Yee [1]. In the case of a 2D
TMz lattice1
, it is possible to derive the following from Maxwell’s equations:
∗ and σ are the magnetic loss (Ω/m) and electric conductivity (S/m) respectively. The
Yee algorithm uses second-order central difference approximations to discretise the spatial and
temporal partial differentiation operators. If ∆x, ∆y (m) are the dimensions of a lattice cell in
the x and y directions respectively, and ∆t (s) is the time step, it is useful to adopt the notation
for a field component U (U may be either E or H) given by
U(x, y, t) = U(i∆x, j∆y, n∆t) = U|
n
i,j .
Consider now the (i, j)th TMz lattice cell, as shown in Fig. 1. Using the above notation, it is
possible to form the update equations [2, p73ff] for the various field components, given by
1A TMz lattice only contains the field components Hx, Hy and Ez. The converse is a TEz lattice which only
contains the field components Ex, Ey and Hz.
(3)
The current term Jsource can be used to excite the lattice. The coefficient matrices Ca(·), Cb(·),
Da(·) and Db(·) are used to incorporate different materials within the lattice and are are given
2.2 Excitation
It is always necessary to excite the lattice in some fashion, and the specific method by which this
is done is problem dependent. One way is to specify the term Jsource according to a predefined
time sequence; alternatively a given field component can be directly specified in a similar fashion
(in the TMz case it is usual to specify a single Ez component). As the simulation progresses,
the field will be observed to propagate outwards from the source.
In many cases it is desirable to obtain the response of a system at a fixed frequency3
. Although it
would seem logical to use sinusoidal excitation to determine this, in practice is is usually better
to estimate the impulse response of the system using a wideband pulse such as a Gaussian
derivative pulse [3, p88] given by
(4)
which provides a unit peak amplitude at t− m = ±s. Plotting the impulse response can provide
a very good qualitative understanding of how the propagating electromagnetic wave interacts
with objects in the problem domain.
If the time-harmonic response of the system is required, this can be straightforwardly determined
from the impulse response. For example, if the time-harmonic response at frequency f is required
for the Ez component, the time-harmonic electric field Ez(f) is given by [4, p169]
Ez(f) = X
T
n=0
[Ez|
n
cos(2πfn∆t) − jEz|
n
sin(2πfn∆t)] (5)
In practice, Ez(f) can be calculated by maintaining a separate complex field buffer which is
incrementally determined by adding the new contribution at each time step.
2.3 Spatial Step Size
In the FDTD method it is necessary to select a spatial step size of approximately λ/20 in the
most electromagnetically ‘dense’ material in the solution domain (i.e. in the region with the
greatest value of ǫr).
2.4 Time Step Size
To ensure stability, it is necessary to select a time step size that is less than or equal to the
Courant limit. In the case of a uniform mesh in 2D with cell size ∆, the Courant limit is given
by [3, p70]
∆tCourant =
∆
u
√
2
2The number of iterations required is problem dependent. Usually it is necessary to perform sufficient iterations
such that any transients have decayed to an acceptable level.
3This is often referred to as the time-harmonic response.
4
where u is the speed of light in the most electromagnetically ‘dense’ material in the lattice. In
practice a time step of ∆t = 0.95∆tCourant is used to ensure any finite precision rounding errors
do not cause numerical instability [5, p31].
2.5 Absorbing Boundary Conditions
The FDTD lattice is, by default, terminated on its periphery by a perfect conductor which
acts to reflect any outwardly propagating fields (can you figure out why?). However, this is not
appropriate for problems which have open boundaries, in which any outwardly propagating fields
should be absorbed. In these cases it is necessary to modify the material coefficient matrices
and/or the update equations in the vicinity of the boundaries to minimize any reflections. The
development of high performance absorbing boundaries has been an area of active research for
some time, and boundaries such as the Uniaxial Perfectly Matched Layer (UPML) [4, p212]
and Convolutional Perfectly Matched Layer (CPML) [4, p225] can achieve very high levels of
performance. These boundaries can, however, be complex to implement.
A much simpler boundary condition is the Absorbing Boundary Condition (ABC) discussed in
[3, pp82-83]. In the case of the boundary at +x, the new value of a tangential field component
i.e., the new field on the boundary is a function only of the old field on the boundary and the
field one lattice cell in from the boundary. In the case of the TMz case being considered here,
this condition need only be specified for Hx at the ±y boundaries and Hy at the ±x boundaries
— the remaining field components are established by the update equations.
3 Assignment
Your report should be in three sections each providing an answer to the following questions.
The assessment criteria for each section are, (a) a clear description of what you have done: 10%,
(b) presentation of your simulation results in an appropriate form for interpretation and discussion: 20%, brief summary of relevant theory researched (including citations of key references,
but avoiding giving an unnecessary tutorial), implementation and calculation of corresponding
theoretical predictions: 30%, (c) discussion of numerical and theoretical results: 30%, and (d)
drawing meaningful conclusions: 10%.
1. Investigate the behaviour of a diffracting PEC knife-edge in the vicinity of its incidence
shadow boundary and in its shadow region. Ensure that the simulation has converged and
examine the field strength in logarithmic units along a vertical and a horizontal receiver
path on the side of the knife-edge opposite to the transmitter. Compare your observations quantitatively against a theoretical model of knife-edge diffraction which you have
researched (Hint: A knife-edge can be specified by setting pec blocks = [150 1 151
225]. Also, plot field values in dBs.) [50 marks]
2. Investigate the behaviour of the field in the shadow region of three cascaded PEC knifeedges defined by the PEC blocks command by,
pec blocks = [150 1 151 225
225 1 226 275
5
300 1 301 250];
Comment on the magnitude of the field compared to the single knife-edge case. [20 marks]
3. Compare the field distribution in the shadow region of three cascaded PEC knife-edges to
the predictions made by the ‘fast’ empirical multiple knife-edge models of Bullington [6],
Epstein-Peterson [7], and Deygout [8]. Discuss the relative accuracy of these three models.
[30 marks]
4 Submission details
You need to submit a short report, no more than 8 sides of A4 excluding figures, in 11-point
Sans Serif font (e.g. Arial), single line spacing and 1.5 cm margins all round. The report should
have a cover and feedback sheet which can be downloaded from the module’s Canvas page and
completed with your student ID number clearly visible on all pages.
Please ensure that all material included from the literature is adequately referenced to avoid
any potential plagiarism penalties.
The report should be submitted in Acrobat PDF format, on Canvas, by the deadline set on the
EMAP module’s Canvas page.
5 Statement of expectations
An excellent report will provide sufficient information to enable the assessor to reproduce its
results independently, will have thoroughly researched the state-of-the-art literature for the appropriate theory to compare with numerical simulation results, will produce insightful discussions
and will draw scientifically sound conclusions.
A failing report will generate numerical simulation results which cannot be verified independently
for corectness, will simply cite relevant literature without justification of its appropriateness, and
will limit its discussions to straight-forward observations.
6 Code Listing — fdtd 1
% fdtd_1
%
% TMz FDTD Code for EE4D Assignment #2
% Code written by Dr. M. J. Neve
%
clear all; % Clear all variables from memory
% Problem parameters
id = ’example1’; % Experiment identification string
f = 1e9; % Frequency (Hz)
samples_per_wavelength = 20; % Samples per wavelength - 20 is good
lattice_size_in_wavelengths = 20; % Lattice size in wavelengths
tmax = 10e-9; % Maximum simulation time (s)
want_abc = 1; % Should absorbing boundary condition be used? 1-yes, 0-no
want_plot_excitation = 1; % Plot excitation waveform?
want_plot_pulse = 1; % Plot pulse waveform at final iteration?
want_plot_time_harmonic = 1; % Plot time harmonic response at frequency f?
want_plot_lattice = 0; % Plot lattice?
want_plot_geometry = 1; % Overplot geometry?
6
want_plot_movie = 0; % Plot/create a movie/animation?
ncontours = 100; % Number of contours
% Constants
c = 2.997925e8; % Speed of light in vacuum (m/s)
eps_0 = 8.854e-12; % Permittivity of free space (F/m)
mu_0 = 4*pi*1e-7; % Permeability of free space (H/m)
eta_0 = sqrt(mu_0/eps_0); % Intrinsic impedance of free space (ohms)
% Excitation
s = 3e-10;
m = 4*s;
xs_idx = 20; % Location of source (Ez field indices)
ys_idx = 200;
% Define PEC blocks
% This section can be used to define PEC blocks
% Each row defines a separate PEC block
% Format is [xmin_idx ymin_idx xmax_idx ymax_idx]
% Indices are specified relative to locations of the Ez component -
% if Ez included components on all sides of cell included
%
%pec_blocks = [200 1 201 200];
pec_blocks = [];
% Preliminary calculations
lambda = c / f; % Wavelength in free space
lx = lattice_size_in_wavelengths * lambda; % x lattice size in m
ly = lx; % y lattice size in m
nx = samples_per_wavelength * lattice_size_in_wavelengths; % Number of samples in x direction
ny = nx;
dx = lx / (nx - 1); % delta x
dy = ly / (ny - 1); % delta y
dt = 0.95 * dx / (c * sqrt(2)); % Time step according to Courant limit
nt = ceil(tmax / dt); % Number of time steps required
% Preliminary calculations are now complete
% Preallocate field storage
ez = zeros(nx - 1, ny - 1); % TMz field components
hx = zeros(nx - 1, ny);
hy = zeros(nx, ny - 1);
ei_z = zeros(size(ez)); % Time harmonic buffer at frequency f
% Preallocation material constants to free space
ca = ones(size(ez));
cb = ones(size(ez)) * dt / (eps_0 * dx);
dax = ones(size(hx));
day = ones(size(hy));
dbx = ones(size(hx)) * dt / (mu_0 * dx);
dby = ones(size(hy)) * dt / (mu_0 * dx);
% Process any defined PEC blocks
[n_pec_blocks, temp] = size(pec_blocks); % n_pec_blocks is the number of PEC blocks
for ii = 1:n_pec_blocks
xmin_idx = pec_blocks(ii,1);
ymin_idx = pec_blocks(ii,2);
xmax_idx = pec_blocks(ii,3);
ymax_idx = pec_blocks(ii,4);
ca(xmin_idx:xmax_idx,ymin_idx:ymax_idx) = -1.0;
cb(xmin_idx:xmax_idx,ymin_idx:ymax_idx) = 0.0;
% dax,dbx,day,dby do not change because magnetic loss of PEC is still 0
end
% Define excitation waveform
t = [0:nt-1] * dt;
v = -exp(0.5)*(t - m) / s .* exp(-(t - m).^2/(2*s^2));
% Main time step loop
for ii = 1:nt
7
ez = ca .* ez + cb .* (hy(2:nx,:) - hy(1:(nx-1),:) + hx(:,2:ny) - hx(:,1:(ny-1)));
ez(xs_idx,ys_idx) = v(ii);
if (want_abc)
hx(:,1) = hx(:,1) * (1 - c * dt / dx) + hx(:,2) * c * dt / dx;
hx(:,ny) = hx(:,ny) * (1 - c * dt / dx) + hx(:,ny-1) * c * dt / dx;
hy(1,:) = hy(1,:) * (1 - c * dt / dx) + hy(2,:) * c * dt / dx;
hy(nx,:) = hy(nx,:) * (1 - c * dt / dx) + hy(nx-1,:) * c * dt / dx;
end
hx(:,2:ny-1) = dax(:,2:ny-1) .* hx(:,2:ny-1) + dbx(:,2:ny-1) .* (ez(:,2:ny-1) - ez(:,1:(ny-2)));
hy(2:nx-1,:) = day(2:nx-1,:) .* hy(2:nx-1,:) + dby(2:nx-1,:) .* (ez(2:nx-1,:) - ez(1:(nx-2),:));
% Update time harmonic buffer
ei_z = ei_z + ez * (cos(2*pi*f*(ii-1)*dt) - j * sin(2*pi*f*(ii-1)*dt));
end
% The remaining code is for plotting/visualisation purposes only
if (want_plot_excitation)
f0 = figure;
plot(t*1e9,v);
xlabel(’t (ns)’);
ylabel(’v (V)’);
title(sprintf(’FDTD: %s: Excitation waveform’, id));
print(f0, ’-depsc2’, sprintf(’%s_excitation.eps’, id))
end
if (want_plot_pulse) % Transpose of data needed to get matrix in correct orientation
f1 = figure;
[ch,ch]=contourf(ez’,ncontours);
set(ch,’edgecolor’,’none’);
axis equal;
%caxis([-0.2 0.2]);
colorbar;
xlabel(’x sample index’);
ylabel(’y sample index’);
title(sprintf(’FDTD: %s: Pulse response’, id));
if (want_plot_geometry)
hold on;
for ii = 1:n_pec_blocks
x1 = (pec_blocks(ii,1) - 1);
y1 = (pec_blocks(ii,2) - 1);
x2 = pec_blocks(ii,3);
y2 = pec_blocks(ii,4);
plot([x1 x2 x2 x1 x1],[y1 y1 y2 y2 y1],’w’);
end
plot(xs_idx, ys_idx, ’wo’);
hold off
end
print(f1, ’-depsc2’, sprintf(’%s_pulse.eps’, id))
end
if (want_plot_time_harmonic)
f2 = figure;
[ch,ch]=contourf(real(ei_z)’,ncontours);
set(ch,’edgecolor’,’none’);
axis equal;
%caxis([-2 2]);
colorbar;
xlabel(’x sample index’);
ylabel(’y sample index’);
title(sprintf(’FDTD: %s: re(ei_z), f = %5.2f GHz’, id, f/1e9));
if (want_plot_geometry)
hold on;
for ii = 1:n_pec_blocks
x1 = (pec_blocks(ii,1) - 1);
y1 = (pec_blocks(ii,2) - 1);
x2 = pec_blocks(ii,3);
8
y2 = pec_blocks(ii,4);
plot([x1 x2 x2 x1 x1],[y1 y1 y2 y2 y1],’w’);
end
plot(xs_idx, ys_idx, ’wo’);
hold off
end
print(f2, ’-depsc2’, sprintf(’%s_timeharmonic.eps’, id))
end
if (want_plot_lattice)
f3 = figure;
hold on;
for ii = 1:nx-1
for jj = 1:ny-1
% hx
x1 = (ii - 1) * dx;
y1 = (jj - 1) * dx;
x2 = ii * dx;
y2 = (jj - 1) * dx;
plot([x1 x2],[y1 y2],’b’);
end
end
hold off
print(f3, ’-depsc2’, sprintf(’%s_lattice.eps’, id))
end
if (want_plot_movie)
f4 = figure;
frame = 0;
mesh(real(ei_z));
set(gca,’nextplot’,’replacechildren’);
for ii = 0:20:340
frame = frame + 1;
theta = ii * pi/180;
[ch,ch]=contourf(real(ei_z * exp(j * theta))’,ncontours);
set(ch,’edgecolor’,’none’);
axis equal;
caxis([-3 3]);
colorbar;
xlabel(’x sample index’);
ylabel(’y sample index’);
title(sprintf(’FDTD: %s: re(ei_z), f = %5.2f GHz’, id, f/1e9));
if (want_plot_geometry)
hold on;
for ii = 1:n_pec_blocks
x1 = (pec_blocks(ii,1) - 1);
y1 = (pec_blocks(ii,2) - 1);
x2 = pec_blocks(ii,3);
y2 = pec_blocks(ii,4);
plot([x1 x2 x2 x1 x1],[y1 y1 y2 y2 y1],’w’);
end
plot(xs_idx, ys_idx, ’wo’);
hold off
end
drawnow;
mov(frame) = getframe;
end
end
9
References
[1] K. Yee, “Numerical solution of initial boundary value problems involving Maxwells equations
in isotropic media,” IEEE Trans. Antennas Propagat., vol. 14, no. 3, pp.302-307, 1966.
[2] A. Taflove and S. C. Hagness, Computational electrodynamics: the finite-difference timedomain method. Boston: Artech House, 2005.
[3] D. B. Davidson, Computational Electromagnetics for RF and Microwave Engineering, 2nd
ed. Cambridge, 2011.
[4] U. S. Inan and R. A. Marshall, Numerical Electromagnetics: The FDTD Method. Cambridge,
2011.
[5] A. C. M. Austin, “Interference Modelling for IndoorWireless Systems using the FiniteDifference Time-Domain Method,” Ph.D. dissertation, Department of Electrical and Computer Engineering, The University of Auckland, New Zealand, 2011.
[6] K. Bullington, “Radio propagation for vehicular communications,” IEEE Trans. Vehicular
Tech., vol. 26, no. 4, pp.295-308, 1977.
[7] J. Epstein and D. W. Peterson, “An Experimental Study of Wave Propagation at 850 MC,”
Proc. IRE, vol. 41, no. 5, pp.595-611, 1953.
[8] J. Deygout, “Multiple knife-edge diffraction of microwaves,” IEEE Trans. Antennas and
Propagat., vol. 14, no. 4, pp.480-489, 1966.