代写data程序、代做Java，Python编程

日期：2023-07-27 09:04

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.

10