Antennas and Electromagnetic Wave Propagation Assignment No.22019
Prof Costas Constantinou
Department of Electronic, ElectricalSystems Engineering
Contents
University of Birmingham
c.constantinoubham.ac.uk
1 Introduction 2
2 The basics behind the FDTD algorithm 2
2.1 DefiningtheLattice ……………………………. 2 2.2 Excitation…………………………………. 4 2.3 SpatialStepSize ……………………………… 4 2.4 TimeStepSize ………………………………. 4 2.5 AbsorbingBoundaryConditions ……………………… 5
3 Examples 5
3.1 example1PropagationinFreeSpacetmax 10ns . . . . . . . . . . . . . . . 5 3.2 example2PropagationinFreeSpacetmax 50ns . . . . . . . . . . . . . . . 5 3.3 example3PropagationinthePresenceofaPECObstacle . . . . . . . . . . . 7
4 Assignment 10
5 Submission details 10
6 Statement of expectations 11
7 Code Listingfdtd 1 11
This assignment was originally devised by Dr Michael Neve of the University of Auckland and has been
slightly modified for the purposes of the present module.
1
1 Introduction
The FiniteDifference TimeDomain FDTD method is a computational electromagnetic tech nique 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 7 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 investiga tions, 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 selfcontained 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 Maxwells equations:
Hx1Ez Hx
t Hy t Ez t
0r y
1 Ez Hy 0r x
1 Hy
o0or x y
whereandare the magnetic loss m and electric conductivity Sm respectively. The Yee algorithm uses secondorder 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
Ux,y,tUix,jy,ntU ni,j.
Consider now the i,jth 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
Ez n0.5Ca i0.5,j0.5Ez n0.5 i0.5,j 0.5 i0.5,j 0.5
Cb i0.5,j0.5 Hy ni,j0.5 Hx ni0.5,j Hx ni0.5,j1
Hy ni1,j0.5 1
Jsourceni0.5,j0.5
Hx Ez
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.
2
Ez i12,j32
Hy i1,j12 Hy i,j12
Hx i12,j1
Ez i12,j12
Ez i12,j12
Hx i12,j i,j
Figure 1: i, jth TMz lattice cell.
Hx n1 i0.5,j 1
Da i0.5,j1Hx ni0.5,j1Db i0.5,j1 Ez n0.5
and Hy n1 i,j 0.5
Da i,j0.5Hy ni,j0.5Db i,j0.5 Ez n0.5
The current term Jsource can be used to excite the lattice. The coefficient matrices Ca , Cb , Daand Dbare used to incorporate different materials within the lattice and are are given by
Cai,j1i,jt 2oi,j
i0.5,j 0.5
Ez n0.52 i0.5,j 32
Ez n0.5. 3 i0.5,j 0.5
i0.5,j 0.5
1 i,j t2oi,j
Cbi,j1i,jt 2oi,j
toi,j
t
1 i,j 2i,j
1 Dai,j t
i,j 2i,j
t
i,jandDbi,j t.
1 i,j 2i,j
It should be noted that the indices in the coefficient matrices correspond to the locations of the field components that are being updated. Although appearing cumbersome, these update
3
E
E E
equations can be programmed in a straightforward fashion. Initially, all field components are initialized to zero, and the field components updated in the order Ez 1 followed by Hx 2 and Hy 3. These calculations are then repeated in sequence until sufficient number of iterations have been performed2.
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
e12 tm2 2s2
v0t stme 4
which provides a unit peak amplitude at t ms. 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 timeharmonic response of the system is required, this can be straightforwardly determined from the impulse response. For example, if the timeharmonic response at frequency f is required for the Ez component, the timeharmonic electric field zf is given by 4, p169
T
zf Ez n cos2fnt jEz n sin2fnt 5
n0
In practice, zf 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 or.
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 u2
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 timeharmonic response. 4
where u is the speed of light in the most electromagnetically dense material in the lattice. In practice a time step of t0.95tCourant 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 andor 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, pp8283. In the case of the boundary at x, the new value of a tangential field component
n1 is given by Nx
n1n 1 ctctn
Nx Nx xxNx1
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 boundariesthe remaining field components are established by the update equations.
3 Examples
A rudimentary FDTD code fdtd 1 has been written in MATLAB and is included in Section 7. Various examples using this code will be investigated in this section.
3.1 example1Propagation in Free Space tmax10 ns
This example is for the code included in Section 7. The source is located at 20,200, and the total simulation time is 10 ns. The excitation waveform is plotted in Fig. 2, and the pulse response at 10 ns is shown in Fig. 3.
The propagating wavefront is clearly seen in Fig. 3, together with evidence of a small reflection which has arisen due to the imperfect nature of the ABC on the lefthand side of the computa tional domain. It is not meaningful to extract the timeharmonic response from this result, as the solution had not reached steady state. This is now considered as example2.
3.2 example2Propagation in Free Space tmax50 ns
The magnitudes of the fields in Fig. 5 are noticeably smaller than those in Fig. 3, as all propa gating fields have encountered the ABC on the periphery of the computational domain at least once. However, the residual field is still of appreciable magnitude, and the only way to reduce
5
1 0.8 0.6 0.4 0.2 0 0.2 0.4 0.6 0.8 1
FDTD: example1: Excitation waveform
0 1 2 3 4 5 6 7 8 9 10 t ns
Figure 2: example1Excitation waveform.
Figure 3: example1Pulse response at t10 ns.
6
v V
FDTD: example2: Excitation waveform
1
0.8
0.6
0.4
0.2
0
0.2
0.4
0.6
0.8
1
0 5 10 15 20 25 30 35 40 45 50
t ns
Figure 4: example2Excitation waveform.
these is to use a higher performance absorbing boundary such as the UPML or CPML. The timeharmonic response in Fig. 6 shows a dominant cylindricallypropagating wave, although some ripple is present and is due to the presence of reflections from the boundaries of the computational domain.
3.3 example3Propagation in the Presence of a PEC Obstacle
A PEC obstacle has been defined with vertices at 150,100 and 300,250. This is done by
including the following definition for pec blocks:
pecblocks150 100 300 250;
As in example2, tmax50 ns, and the pulse response at 50 ns is plotted in Fig. 7 and the time harmonic response in Fig. 8 the excitation waveform is the same as in example2. The pulse response in Fig. 7 is somewhat complex, as a result from waves reflecting from and diffracting around the PEC block. The effects of reflection can also be seen in Fig. 8 with the presence of a standing wave between the source and the box, and a significantly reduced field amplitude behind the box as a result of diffraction4.
4An animation of the timeharmonic fields can also be useful in visualizing these effects. This can be done by enabling the option want plot movie1; in the header, and then executing the command moviemov,n where n is an integer specifying the number of times the movie should be played.
7
v V
350
300
250
FDTD: example2: Pulse response
4 x 10
10
8
6
4
200 2
150
100
50
50 0 50 100 150 200 250 300 350 400 450
0
2
4
6
350
300
250
200
150
100
50
10
8
6
4
2
0
2
x sample index
Figure 5: example2Pulse response at t50 ns.
FDTD: example2: reei , f1.00 GHz z
y sample index y sample index
50 0 50 100 150 200 250 300 350 400 450
x sample index
Figure 6: example2Timeharmonic response.
8
350
300
250
200
150
100
50
350
300
250
200
150
100
50
10
8
6
4
2
0
2
FDTD: example3: Pulse response
4 x 10
1
0
1
50 0 50 100 150 200 250 300 350 400 450
x sample index
Figure 7: example3Pulse response at t50 ns.
FDTD: example3: reei , f1.00 GHz z
y sample index y sample index
50 0 50 100 150 200 250 300 350 400 450
x sample index
Figure 8: example3Timeharmonic response.
9
4 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 discus sion: 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 conclusions: 10.
5
1.
2.
3.
4.
Implement the WalfischBertoni propagation model 6 for a typical UK scenario of cellular telephony coverage at 900 MHz, for an area such as Selly Oak. In this scenario, the model parameters shown in Fig. 1 of 6 are R500 m, d25 m, h10 m, H5 m, and hm10 m. Use the theoretical model in 6 to determine the attenuation function relative to free space in dB at a mobile receiver centred at street level. Note that the model assumes that the buildings are knifeedge obstacles of zero thickness. 25 marks
Extend the freespace simulation space of example1 so that it corresponds to the prop agation scenario of 1 by selecting the correct lattice size, frequency, source location, simulation convergence time, etc. Avoid having the source too close to the simulation boundaries, and discuss why this is a sensible precaution. Calculate the field value at the location which corresponds to the mobile receiver in 1, assuming the complete absence of the ground. 25 marks
Simulate the WalfischBertoni propagation scenario using the FDTD code, assuming that the buildings can be modelled as perfectly conducting knifeedge block of thickness equal to 1 lattice interval and plot the field strength in the entire domain in logarithmic units. to calculate the field at the receiver point, and thus the attenuation relative to free space in dB. Compare this value to the one predicted using the WalfischBertoni model in 1 and discuss the source of any discrepancy. 25 marks
Repeat 3 but now replace the knifeedges with rectangular blocks of 8.5 m thickness, centred on the location of each knifeedge. Discuss the observed difference if any between the attenuation function obtained in this case and the corresponding estimates obtained in 1 and 3. 25 marks
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 modules 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 12:00 noon on 10 February 2020.
10
6 Statement of expectations
An excellent report will provide sufficient information to enable the assessor to reproduce its results independently, will have thoroughly researched the stateoftheart literature for the ap propriate 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 straightforward observations.
7 Code Listingfdtd 1
fdtd1
TMz FDTD Code for EE4D Assignment 2
Code written by Dr. M. J. Neve
clear all;Clear all variables from memory
Problem parameters
idexample1;
f1e9;
samplesperwavelength20; latticesizeinwavelengths20; tmax10e9;
wantabc1; wantplotexcitation1; wantplotpulse1; wantplottimeharmonic1; wantplotlattice0; wantplotgeometry1; wantplotmovie0; ncontours100;
Constants
c2.997925e8;
eps08.854e12;
mu04pi1e7;
eta0sqrtmu0eps0;
Excitation
s3e10;
m4s;
xsidx20;Location of source Ez field indices ysidx200;
Experiment identification string
Frequency Hz
Samples per wavelength20 is good
Lattice size in wavelengths
Maximum simulation time s
Should absorbing boundary condition be used? 1yes, 0noPlot excitation waveform?
Plot pulse waveform at final iteration?
Plot time harmonic response at frequency f?
Plot lattice?
Overplot geometry?
Plotcreate a movieanimation?
Number of contours
Speed
Permittivity of free space Fm
Permeability of free space Hm
Intrinsic impedance of free space ohms
of light in vacuum ms
pecblocks200 1 201 200; pecblocks;
Preliminary calculations
lambdacf;Wavelength in free space
lxlatticesizeinwavelengthslambda;
lylx; ylatticesizeinm
nxsamplesperwavelengthlatticesizeinwavelengths;Number of samples in x direction nynx;
Define
This
Each
PEC blocks
section can be used to define PEC blocks row defines a separate PEC block
Format is xminidx yminidx xmaxidx ymaxidx
Indices are specified relative to locations of the Ez component
if Ez included components on all sides of cell included
x lattice size in m
11
dxlxnx1; deltax
dylyny1; deltay
dt0.95dxcsqrt2;Time step according to Courant limit ntceiltmaxdt;Number of time steps required
Preliminary calculations are now complete
Preallocate field storage
ezzerosnx1, ny1;TMz field components
hxzerosnx1, ny;
hyzerosnx, ny1;
eizzerossizeez;Time harmonic buffer at frequency f
Preallocation material constants to free space caonessizeez;
cbonessizeezdteps0dx;
daxonessizehx;
dayonessizehy;
dbxonessizehxdtmu0dx; dbyonessizehydtmu0dx;
Process any defined PEC blocks npecblocks, tempsizepecblocks; for ii1:npecblocks
npecblocks is the number of PEC blocks
xminidxpecblocksii,1;
yminidxpecblocksii,2;
xmaxidxpecblocksii,3;
ymaxidxpecblocksii,4;
caxminidx:xmaxidx,yminidx:ymaxidx1.0; cbxminidx:xmaxidx,yminidx:ymaxidx0.0;
dax,dbx,day,dby do not change because magnetic loss of PEC is still 0
end
Define excitation waveform
t0:nt1dt;
vexp0.5tms . exptm.22s2;
Main time step loop for ii1:nt
ezca . ezcb . hy2:nx,:hy1:nx1,:hx:,2:nyhx:,1:ny1;
ezxsidx,ysidxvii;
if wantabc
hx:,1hx:,11cdtdxhx:,2cdtdx; hx:,nyhx:,ny1cdtdxhx:,ny1cdtdx; hy1,:hy1,:1cdtdxhy2,:cdtdx; hynx,:hynx,:1cdtdxhynx1,:cdtdx;
end
hx:,2:ny1dax:,2:ny1 . hx:,2:ny1dbx:,2:ny1 . ez:,2:ny1ez:,1:ny2; hy2:nx1,:day2:nx1,: . hy2:nx1,:dby2:nx1,: . ez2:nx1,:ez1:nx2,:;
Update time harmonic buffer
eizeizezcos2pifii1dtjsin2pifii1dt;
end
The remaining code is for plottingvisualisation purposes only
if wantplotexcitation f0figure;
plott1e9,v;
xlabelt ns;
ylabelv V;
titlesprintfFDTD: s: Excitation waveform, id; printf0, depsc2, sprintfsexcitation.eps, id
end
if wantplotpulseTranspose of data needed to get matrix in correct orientation
12
f1figure; ch,chcontourfez,ncontours; setch,edgecolor,none; axis equal;
caxis0.2 0.2;
colorbar;
xlabelx sample index; ylabely sample index; titlesprintfFDTD: s: Pulse if wantplotgeometry
hold on;
for ii
x1
y1
x2
y2
1:npecblocks
pecblocksii,1
pecblocksii,2
pecblocksii,3;
pecblocksii,4;
response, id;
1;1;
plotx1 x2 x2 x1 x1,y1 y1 y2 y2 y1,w; end
plotxsidx, ysidx, wo;
hold off end
printf1, depsc2, sprintfspulse.eps, id end
if wantplottimeharmonic f2figure;
ch,chcontourfrealeiz,ncontours; setch,edgecolor,none;
axis equal;
caxis2 2;
colorbar;
xlabelx sample index;
ylabely sample index;
titlesprintfFDTD: s: reeiz, f5.2f GHz, id, f1e9; if wantplotgeometry
hold on;
for ii
x1
y1
x2
y2
1:npecblocks
pecblocksii,1
pecblocksii,2
pecblocksii,3;
pecblocksii,4;
1;1;
plotx1 x2 x2 x1 x1,y1 y1 y2 y2 y1,w; end
plotxsidx, ysidx, wo;
hold off end
printf2, depsc2, sprintfstimeharmonic.eps, id end
if wantplotlattice f3figure;
hold on;
for ii1:nx1
for jj
hx
x1
y1
1:ny1
ii1dx;
jj1dx;
x2iidx;
y2jj1dx; plotx1 x2,y1 y2,b;
end end
hold off
printf3, depsc2, sprintfslattice.eps, id end
if wantplotmovie f4figure;
frame0;
meshrealeiz;
13
setgca,nextplot,replacechildren; for ii0:20:340
frameframe1;
thetaiipi180; ch,chcontourfrealeizexpjsetch,edgecolor,none;
axis equal;
caxis3 3;
colorbar;
xlabelx sample index;
ylabely sample index; titlesprintfFDTD: s: reeiz, f if wantplotgeometry
theta,ncontours;
5.2f GHz, id, f1e9;
y2 y2 y1,w;
hold on;
for ii
x1
y1
x2
y2
1:npecblocks pecblocksii,11; pecblocksii,21; pecblocksii,3; pecblocksii,4;
plotx1 x2 x2 x1 x1,y1 y1
end
plotxsidx, ysidx, wo; hold off
end
drawnow;
movframegetframe;
end end
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.302307, 1966.
2 A. Taflove and S. C. Hagness, Computational electrodynamics: the finitedifference time domain 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 Finite Difference TimeDomain Method, Ph.D. dissertation, Department of Electrical and Com puter Engineering, The University of Auckland, New Zealand, 2011.
6 J. Walfisch and H. L. Bertoni, A theoretical model of UHF propagation in urban environ ments, IEEE Trans. Antennas Propagat., vol. 36, no. 12, pp.17881796, 1988.
14
Reviews
There are no reviews yet.