DEVELOPMENT OF COMPLEX ELECTROMAGNETIC PROBLEMS USING FDTD SUBGRIDDING IN HYBRID COMPUTATIONAL TECHNIQUES Khairan N Ramli Raed A Abd Alhameed Pet PDF
DEVELOPMENT OF COMPLEX ELECTROMAGNETIC PROBLEMS USING FDTD SUBGRIDDING IN HYBRID COMPUTATIONAL TECHNIQUES Khairan N Ramli Raed A Abd Alhameed Pet PDF
DEVELOPMENT OF COMPLEX ELECTROMAGNETIC PROBLEMS USING FDTD SUBGRIDDING IN HYBRID COMPUTATIONAL TECHNIQUES Khairan N Ramli Raed A Abd Alhameed Pet PDF
DEVELOPMENT OF COMPLEX
ELECTROMAGNETIC PROBLEMS
USING FDTD SUBGRIDDING
IN HYBRID COMPUTATIONAL
TECHNIQUES
No part of this digital document may be reproduced, stored in a retrieval system or transmitted in any form or
by any means. The publisher has taken reasonable care in the preparation of this digital document, but makes no
expressed or implied warranty of any kind and assumes no responsibility for any errors or omissions. No
liability is assumed for incidental or consequential damages in connection with or arising out of information
contained herein. This digital document is sold with the clear understanding that the publisher is not engaged in
rendering legal, medical or any other professional services.
ENGINEERING TOOLS, TECHNIQUES
AND TABLES
DEVELOPMENT OF COMPLEX
ELECTROMAGNETIC PROBLEMS
USING FDTD SUBGRIDDING
IN HYBRID COMPUTATIONAL
TECHNIQUES
KHAIRAN N. RAMLI,
RAED A. ABD-ALHAMEED
AND
PETER S. EXCELL
EDITORS
New York
Copyright © 2014 by Nova Science Publishers, Inc.
All rights reserved. No part of this book may be reproduced, stored in a retrieval system or transmitted
in any form or by any means: electronic, electrostatic, magnetic, tape, mechanical photocopying,
recording or otherwise without the written permission of the Publisher.
For permission to use material from this book please contact us:
Telephone 631-231-7269; Fax 631-231-8175
Web Site: http://www.novapublishers.com
Independent verification should be sought for any data, advice or recommendations contained in this
book. In addition, no responsibility is assumed by the publisher for any injury and/or damage to
persons or property arising from any methods, products, instructions, ideas or otherwise contained in
this publication.
This publication is designed to provide accurate and authoritative information with regard to the subject
matter covered herein. It is sold with the clear understanding that the Publisher is not engaged in
rendering legal or any other professional services. If legal or any other expert assistance is required, the
services of a competent person should be sought. FROM A DECLARATION OF PARTICIPANTS
JOINTLY ADOPTED BY A COMMITTEE OF THE AMERICAN BAR ASSOCIATION AND A
COMMITTEE OF PUBLISHERS.
Additional color graphics may be available in the e-book version of this book.
Preface vii
Chapter 1 Numerical Solution of Maxwell Equations
Using FDTD and MoM 1
K. N. Ramli, R. A. Abd-Alhameed
and P. S. Excell
Chapter 2 FDTD Technique for Field Truncation 17
K. N. Ramli, R. A. Abd-Alhameed
and P. S. Excell
Chapter 3 Surface Kernel Solution of the Method of Moments 45
K. N. Ramli, M. Lashab, K. H. Sayidmarie
and R. A. Abd-Alhameed
Chapter 4 Quasi-Static Finite-Difference Time-Domain
Subgridding Technique 113
K. N. Ramli, R. A. Abd-Alhameed
and P. S. Excell
Chapter 5 Effect of Anisotropic Magneto-Chirality
on Microstrip Resonator Characteristics 143
C. Zebiri, F. Benabdelaziz and M. Lashab
Chapter 6 Interaction of EM Fields to the Human Body
Using Hybrid Computational Method 207
K. N. Ramli, R. A. Abd-Alhameed
and P. S. Excell
vi Contents
numbers of time steps, in the order of tens of millions, has been eased to tens
of thousands by employing quasi-static methods.
The book also illustrates the principle of the equivalent surface boundary
employed close to the antenna for MoM-FDTD-SGFDTD hybridisation. It
shows the benefit of using hybrid methods due to their ability to analyse a
system of multiple discrete regions by employing the principle of equivalent
sources to excite the coupling surfaces. A human body model is developed to
interact with a short range RFID antenna. The near and far field radiation
pattern are explored for which the cumulative distribution function of antenna
radiation efficiency is presented. The field distributions show reasonable and
stable results at 900 MHz. This work helps deeper investigation of the
phenomena in the interaction between electromagnetic fields and human
tissues.
FDTD method is further used at powerline and mobile communication
frequency to model the interaction between the EM wave and electrically
small objects. The method is also employed for lumped-element biological cell
model in which Floquet periodic boundary conditions are imposed to introduce
a periodic replication effect of the biological cell. Furthermore, the hybrid
MoM-FDTD technique permits different parts of an antenna to be computed.
The Equivalence Principle describes the fields between the different parts
accurately. The method facilitates precise computation of antenna in proximity
with lossless and lossy dielectric media.
In: Development of Complex Electromagnetic … ISBN: 978-1-61122-013-1
Editors: K. N. Ramli et al. © 2014 Nova Science Publishers, Inc.
Chapter 1
ABSTRACT
This chapter briefly describes the historical background of the
numerical solution for Maxwell equations using FDTD and MoM
techniques in solving complex electromagnetic scattering problems. It
also presents the state of the art of the original contributions of the
analysis methods towards further developments of the hybrid
computational technquiesin Electromagnetics.
2 K. N. Ramli, R. A. Abd-Alhameed and P. S. Excell
problem space of the computational domain. In the past, MoM has been
applied for aircraft field analysis [15], dual-reflector antenna and feeding horn
[16] and microstrip antennas [17]. Moreover, FDTD has been employed for
patch antennas [18], sloped interfaces [19-20] and SAR calculation [21].
Generally, the algorithm used by Yee was described by the electric field
component which was spatially and temporally offset from the magnetic field
component to acquire the update equations. These equations were used in a
leap-frog manner to propagate the electric and magnetic fields ahead in time.
The equations provide the present fields in terms of the past fields all over the
computational domain. After Yee’s publication, the approach was widely used
with different endeavour [23-27].
The boundaries of the computational domain in FDTD need to be
carefully treated when simulating problems in open regions. Spurious
4 K. N. Ramli, R. A. Abd-Alhameed and P. S. Excell
reflections will generally occur from the termination of the grid. The problem
can be solved by means of the well known method called the absorbing
boundary condition (ABC). It is generally meant to absorb any outgoing
propagating waves without ideally producing spurious reflections. The ABCs
was first proposed in 1971 by Merewether [28] to solve the open region
difficulties. The development chronicle to magnify the practicability study of
the technique was continued in the literature by [29-33] which were based of
nonmaterial type.
In contrast, Berenger presented a new idea in 1994 called the perfectly
matched layer (PML) ABC which was based on material category [34]. The
state of the art of Berenger’s PML contributes to notably better precision when
compared to the other ABCs in the written works [35, 36] for broad assortment
of applications.
The main handicap of FDTD lies in the truth that only consistent grids can
be used. Accordingly, the geometry resemblance in FDTD is restricted to
staircase-shaped boundaries which lead to a large number of computer
memory requirements and the CPU time particularly when dealing with
curvature geometries with fine features [37]. The total number of cells in the
computational domain grows significantly due to a global fine mesh. Another
FDTD weakness is the presence of error due to numerical dispersion [38, 39].
In this case, many scientists were prompted to examine the subgridding
scheme as an approach to parry the problem. A variety of methods have been
proposed to boost the efficiency of FDTD technique such as non-uniform
meshing [40], sub-cellular technique [41], non-orthogonal meshing [42],
alternative direction implicit (ADI) method [43-45], higher-order technique
[46, 47], hybrid method [48-52] and subgridding method [53-59]. In FDTD
subgridding technique, the smaller size components in a structure is filled with
fine grids and the remaining of the space is represented by coarse grids. The
fields on the boundary between coarse and fine grids are basically unknown in
nature. They are predicted by using spatial and temporal interpolations.
The regions of the coarse and fine grids are computed by the FDTD
method and are kept in time step which satisfy the Courant stability condition.
Consequently, the stable subgridding algorithm can refine the mesh locally.
Hence, the accuracy of the solution can be improved without increasing the
computational efforts significantly.
Numerical Solution of Maxwell Equations Using FDTD and MoM 5
REFERENCES
[1] A. Taflove and S. C. Hagness, Computational electrodynamics: The
finite-difference time-domain method, 3rd ed., Boston, MA: Artech
House, 2005.
[2] K. Tap, P. Pathak and R. Burkholder, “Complex source beam - moment
method procedure for accelerating numerical integral equation solutions
of radiation and scattering problems,” IEEE Transactions on Antennas
and Propagation, vol. 62, no. 4, pp. 2052 - 2062, 2014.
[3] D. J. Ludick, R. Maaskant, D. B. Davidson, U. Jakobus, R. Mittra, and
D. De Villiers, “Efficient analysis of large aperiodic antenna arrays
using the domain Green's function method,” IEEE Transactions on
Antennas and Propagation, vol. 62, no. 4, part 1, pp. 1579 - 1588, 2014.
Numerical Solution of Maxwell Equations Using FDTD and MoM 9
[4] X.-M. Pan and X.-Q. Sheng, “Improved algebraic preconditioning for
MoM solutions of large-scale electromagnetic problems,” IEEE
Antennas and Wireless Propagation Letters, vol. 13, pp. 106-109, 2014.
[5] M. A. Echeverri Bautista, M. A. Francavilla, F. Vipiana, and G. Vecchi,
“A hierarchical fast solver for EFIE-MoM analysis of multiscale
structures at very low frequencies,” IEEE Transactions on Antennas and
Propagation, vol. 62, no. 3, pp. 1523-1528, 2014.
[6] M. Shafieipour, I. Jeffrey, J. Aronsson, and V. I. Okhmatovski, “On the
equivalence of RWG method of moments and the locally corrected
nyström method for solving the electric field integral equation,” IEEE
Transactions on Antennas and Propagation, vol. 62, no. 2, pp. 772-782,
2014.
[7] N. Hendijani, J. Cheng and R. J. Adams, “Combined field integral
equation using a constraint-based Helmholtz decomposition,” IEEE
Transactions on Antennas and Propagation, vol. 62, no. 3, pp. 1500-
1503, 2014.
[8] J. Alvarez, L. D. Angulo, M. R. Cabello, A. Rubio Bretones, and S. G.
Garcia, “An analysis of the leap-frog discontinuous Galerkin method for
Maxwell's equations,” IEEE Transactions on Microwave Theory and
Techniques, vol. 62, no. 2, pp. 197-207, 2014.
[9] R. Xiong, B. Chen, Y. Yi, H.-L. Chen, and H.-F. Zhao, “Two-step
method for the FDTD analysis of long apertures having depth,” IEEE
Antennas and Wireless Propagation Letters, vol. 12, pp. 39-42, 2013.
[10] A. C. Lesina, A. Vaccari and A. Bozzoli, “A modified RC-FDTD
algorithm for plasmonics in Drude dispersive media nanostructures,” 6th
European Conference on Antennas and Propagation (EUCAP), pp. 687-
690, 2012.
[11] G. M. Noetscher, Y. Xu and S. N. Makarov, “Accuracy of point source
models with coincident phase centers in a cubic FDTD grid for arbitrary
source orientation,” IEEE Antennas and Propagation Society
International Symposium (APSURSI), pp. 1-2, 2012.
[12] V. Nayyeri, M. Soleimani and M. Dehmollaian, “Implementation of a
PEMC boundary condition in the 2-D FDTD technique,” IEEE Antennas
and Propagation Society International Symposium (APSURSI), pp. 1-2,
2012.
[13] M. Weiwei, S. Dong and W. Xianliang, “UPML-FDTD parallel
computing on GPU,” International Conference on Microwave and
Millimeter Wave Technology (ICMMT), vol. 4, pp. 1-4, 2012.
10 K. N. Ramli, R. A. Abd-Alhameed and P. S. Excell
Chapter 2
2.1. INTRODUCTION
Over the past few years, finite-difference time-domain (FDTD) method
[1-11] have become increasingly prevalent in the computational
electromagnetic problems due to its simplicity, efficiency, robustness and
versatility scheme for highly complex configuration in the computational
18 K. N. Ramli, R. A. Abd-Alhameed and P. S. Excell
B
E J m
t (2.1)
D
H Je
t (2.2)
FDTD Technique for Field Truncation 19
Jm H
(2.3)
Je E (2.4)
B=H (2.5)
D=E (2.6)
H
t
1
E H
(2.7)
*H is the magnetic losses which may exist inside the medium. Inserting
(2.4) and (2.6) to (2.2) and dividing by gives:
E 1
H E
t (2.8)
E is the electric losses which may exist inside the medium. In Cartesian
coordinates, equations (2.7) and (2.8) yield the following six scalar equations:
20 K. N. Ramli, R. A. Abd-Alhameed and P. S. Excell
H x 1 E y E z
H x
t z y (2.9)
H y 1 E z E x
H y
t x z (2.10)
H z 1 E x E y
H z
t y x (2.11)
Ex 1 H z H y
Ex
t y z (2.12)
E y 1 H x H z
Ey
t z x (2.13)
Ez 1 H y H x
Ez
t x y (2.14)
(i , j , k ) (i x , j y , k z ) (2.15)
Let F denote any function of discrete space and time given by:
F Fi 1 2 , j , k Fi 1 2, j , k
n n
O ( x) 2
x x (2.17)
n 1 2 n 1 2
F Fi , j , k Fi , j , k
O ( t ) 2
t t (2.18)
In equation (2.17), O(x) 2 is the error term that represents all the
remaining terms in a Taylor series expansion. It is known as a central finite
difference scheme in space with second-order accuracy. Similarly, (2.18) is
second-order accurate in time. Applying Yee's finite-difference scheme to
(2.9) gives:
E n n
y i , j , k 1 2 E y i , j , k 1 2
n 1 2 n 1 2
H x i, j ,k H x i, j ,k 1
z
t i, j ,k E n n
z i , j 1 2, k E z i , j 1 2,k n
i, j , k H x i , j ,k
y (2.19)
n
The H x i , j ,k
field component in (2.19) is evaluated at time step n.
n
However, the value of H x i , j ,k
at time step n is not available and hence the
following interpolated approximation is used:
n 1 2 n 1 2
n
H x i , j ,k H x i , j ,k
H x i , j ,k
2 (2.20)
n 1 2
By substituting equation (2.20) in (2.19), leaving H x i , j ,k
on the left
hand side and passing the all remaining terms to the right, assuming cubical
FDTD cells are used, the finite difference updating equation for the magnetic
and electric field components can be derived as:
22 K. N. Ramli, R. A. Abd-Alhameed and P. S. Excell
E n n
n1 2 n1 2 y i, j ,k 1 2 Ey i, j ,k 1 2
H x i, j ,k Da, H x H x i, j ,k Db, H x
i , j ,k i , j ,k n n
Ez i, j 1 2,k Ez i, j 1 2,k
(2.21)
Ez n n
n 1 2 n 1 2 i 1 2, j , k E z i 1 2, j , k
H y i , j ,k Da , H y H y i , j ,k Db, H y
i , j ,k n n
i, j ,k
E x i , j ,k 1 2 E x i , j , k 1 2
(2.22)
Ex n n
n 1 2 n 1 2 i , j 1 2,k Ex i , j 1 2,k
H z i , j ,k Da , H z H z i , j ,k Db , H z
i , j ,k i , j ,k n n
E y i 1 2, j ,k E y i 1 2, j ,k
(2.23)
H z n1 2 H z n1 2
n 1 n i , j 1 2,k i , j 1 2, k
E x i , j , k Ca , E x Ex i , j ,k Cb, Ex
i , j ,k i , j ,k n 1 2 n 1 2
H y i , j ,k 1 2 H y i , j ,k 1 2
(2.24)
H x n 1 2 H x n 1 2
n 1 n i , j , k 1 2 i , j , k 1 2
Ey Ca , E y Ey Cb , E y
i , j ,k n 1 2 n 1 2
i , j ,k i , j ,k i , j ,k
H z i 1 2, j ,k H z i 1 2, j ,k
(2.25)
H n1 2 H n 1 2
n 1 n y i 1 2, j ,k y i 1 2, j , k
E z i , j ,k Ca , E z E z i , j ,k Cb, Ez
H x i ,j112 2,k H x in,j112 2,k
i , j ,k i , j ,k n
(2.26)
It can be seen that the coefficients on the left hand side are referred to as
Yee's updating coefficients. The electric field coefficients are given by:
t i , j ,k t
C a i , j , k 1 i , j , k 1
2 i , j ,k 2
i, j ,k
(2.27)
FDTD Technique for Field Truncation 23
t i , j , k t
Cb p 1
i , j ,k 2 i , j ,k
i , j ,k p (2.28)
i, j , k t i, j , k t
Da 1 1
i, j ,k 2 i, j ,k 2 i, j,k
(2.29)
t i, j ,k t
Db p 1
i , j ,k i , j ,k p 2i , j ,k
(2.30)
Figure 2.1. The electric and magnetic field components distribution on the FDTD
lattice [16].
Figure 2.2. Relationship between field components: (a) within a quarter of a unit cell,
(b) on a plane [16].
FDTD Technique for Field Truncation 25
Figure 2.3. Space-time chart of the Yee algorithm in a leapfrog arrangement [11].
The field quantities are solved with a leapfrog arrangement where a half
time step separates the solutions of the electric and magnetic fields as shown
in Figure 2.3. The time-stepping of the FDTD algorithm is continued until the
desired late time pulse response or steady state behaviour is reached. Figure
2.4 depicts the time-stepping FDTD algorithms flow chart. Apparently, a
suitable size of the time step should be chosen properly to avoid the late time
instability of the algorithms, after determining the spatial resolution based on
the geometrical features and the operating frequency. It can be seen that a
leapfrog arrangement between the electric and magnetic field vector
components is used to implement the time step of the FDTD algorithm. The
grids of the electric and magnetic fields are displaced half a cell between them
indicating that the computer must work through them in turn.
26 K. N. Ramli, R. A. Abd-Alhameed and P. S. Excell
Figure 2.4. The flow chart of the time stepping FDTD algorithm [13].
FDTD Technique for Field Truncation 27
Figure 2.5. Variation of the numerical phase velocity with wave propagation angle in
two-dimensional FDTD grid [11].
dispersion error, when the grid size and the angle of propagating wave are
increased as illustrated in Figure 2.5. Numerical dispersion is observed to be
approximately 1% when the cell size is exactly /10.
The time increment t must satisfy the CFL criterion [11, 13] to ensure
the stability of the computation process, which sets the relation between the
time and space resolution for three-dimensional FDTD given by:
1
vmax t
1 1 1
c
x y z 2
2 2
(2.31)
h
t
vmax d (2.32)
The second hard source provides a wideband Gaussian pulse with finite
DC content. The pulse is centred at the time step no and has a 1/e characteristic
decay of ndecay time steps. It is simply given by:
Ez
n
Eoe
n n o / n decay 2
is
(2.34)
Ez
n
E0 e
n n0 / ndecay
2
sin 2f 0 n n0 t
is
(2.35)
The plane wave source excitation model was used by Yee in 1966 [15].
An arbitrary incident plane wave is used for the excitation purposes. The
30 K. N. Ramli, R. A. Abd-Alhameed and P. S. Excell
incident plane wave is modelled and approximated within the FDTD space
lattice. Two methods of excitation of the plane wave on the space lattice is
discussed here, namely the initial condition technique and the total/scattered
field formulation technique. The initial condition method was used by Yee in
the implementation of the FDTD method to represent an incident plane wave
with special applications in the RCS simulations. The values of electric and
magnetic field vector components at the zero time step of the incident plane
wave throughout the problem space are preset in sign and magnitude. As a
result, this will give the desired polarisation of the incident plane wave.
However, this method has a drawback due to the fact that it is a non-compact
wave source. It requires a large number of additional free space cells to
contain initial conditions of long duration pulses or continuous sinusoids for
oblique incident angles of the propagating wave. This will unfortunately
increase the size of the problem space which one must avoid. Hence, the
technique is limited to special usage only. The most popular method used in
many FDTD software is the total/scattered field formulation method [17, 18]
due to the fact that it allows the modelling of the FDTD with long-duration
pulsed or sinusoidal illuminations for arbitrary plane wave propagation
directions.
Figure 2.6. Total and scattered field zoning for a generic scattering case [11].
FDTD Technique for Field Truncation 31
Figure 2.6 presents the total and scattered field zoning for a generic
scattering case. Basically, it shows two distinct regions of interest namely the
total field and scattered field region. Both of the regions are separated by a
non-physical virtual surface applied numerically with a special treatment to
include the incident wave excitation and to split the problem space into the
total field and scattered field regions. As a result, the scattered field vector
values can be computed in the scattered field region without the presence of
the incident field. Furthermore, the arbitrary incident plane wave with different
oblique incidence angles using incident-field array (IFA) excitation scheme
[11] can be modelled efficiently. The IFA is generally an FDTD-based lookup
table from which the incident-field values are covered on the FDTD lattice in
the propagating direction.
Figure 2.7. Total and scattered field components for 1-D FDTD grid.
n1
Ez tot,i Ez tot,i
L
n
L
t
Hy
n1 2
tot,iL 1 2
Hy
n1 2
scat,iL 1 2
(2.36)
Subscripts tot and scat stand for total and scattered fields respectively. It is
clear that the above equation is inconsistent, so Hy,tot at (iL-1/2) must be
subtracted from Hy,tot at (iL+1/2) to advance the value of Ez,tot at iL. To correct
this updating equation, the vector function Hy,inc at (iL-1/2) is added as an
excitation wave source of the FDTD algorithm. The boundary Ez updating
equation is given by:
n1
Ez tot,i Ez tot,i
L
n
L
t
Hy
n1 2
tot ,iL 1 2
Hy
n1 2
scat,i L 1 2
t
Hy
n1 2
inc,iL 1 2
(2.37)
Subscript inc stands for incident field. The added term is assumed to be a
known function such as sinusoidal wave or Gaussian pulse for plane wave
representation, while the rest of the terms of the right hand side are assumed
stored in computer memory from the previous updating time step. By
following proper modifications of equation (2.37), the updating equations for
the other three special magnetic and electric boundary field components are
expressed as follows:
Hy
n1 2
scat,iL 1 2
Hy
n1 2
scat,iL 1 2
t
n n
Ez tot,i Ez scat,i 1
L L
t
n
Ez inc,i
L
(2.38)
n 1
Ez tot,i Ez tot,i
R
n
R
t
Hy
n 1 2
scat,i R 1 2
Hy
n 1 2
tot ,i R 1 2
t
Hy
n 1 2
inc,iR 1 2
(2.39)
Hy
n1 2
scat,iR 1 2
Hy
n1 2
scat,iR 1 2
t
n n
Ez scat,i 1 Ez tot,i
R R
t
n
Ez inc,i
R
(2.40)
Figure 2.9. Location of Ey() and Ez() components in planes i = io and i = i1.
34 K. N. Ramli, R. A. Abd-Alhameed and P. S. Excell
n 1 n 1 n 1 2
Ey {E y }(2.25) Cb, E y H z ,inc i 1 2, j ,k
i0 , j , k i0 , j , k i0 , j , k 0
(2.41)
n 1 n 1 n 1 2
E z i , j ,k {E z i , j , k }( 2.26) Cb , E z H y ,inc
i0 , j , k i0 1 2, j , k
0 0
(2.42)
Figure 2.10. Location of Ex() and Ez() components in planes j = jo and j = j1.
n 1 n 1 n 1 2
Ez i1 , j , k
{E z }
i , j , k ( 2.26 )
Cb , E z H y ,inc
i1 , j , k i1 1 2, j , k
1
(2.44)
n 1 n 1 n 1 2
Ex i, j {E x i , j ,k }( 2.24) Cb , Ex H z ,inc i , j
0 ,k i , j0 , k 0 1 2, k
0
(2.45)
FDTD Technique for Field Truncation 35
n 1 n 1 n 1 2
Ez i, j {E z i , j ,k }( 2.26) Cb , E z H x ,inc i , j
0 ,k i , j0 , k 0 1 2, k
0
(2.46)
n 1 n 1 n 1 2
E x i , j ,k {E x i , j ,k }( 2.24) Cb , E x H z ,inc i , j 1 2,k
i , j1 , k
1 1 1
(2.47)
n 1 n 1 n 1 2
Ez i , j1 , k
{E z }
i , j , k ( 2.26 )
Cb , E z H x ,inc i , j 1 2,k
i , j1 , k
1 1
(2.48)
Figure 2.11. Location of Ex() and Ey() components in planes k = ko and k = k1.
n 1 n 1 n 1 2
E x i , j ,k {E x i , j ,k }( 2.24) Cb , E x H y ,inc
i , j ,k0 i , j , k 0 1 2
0 0
(2.49)
n 1 n 1 n 1 2
Ey {E y }( 2.25) Cb, E y H x ,inc i , j ,k
i , j , k0 i , j ,k0 i , j , k0 0 1 2
(2.50)
Ex (i = io + ½, …, i1 - ½; j = jo, …, j1; k = k1) is given by:
36 K. N. Ramli, R. A. Abd-Alhameed and P. S. Excell
n 1 n 1 n 1 2
E x i , j ,k {E x i , j ,k }( 2.24 ) Cb , E x H y ,inc i , j ,k 1 2
i , j , k1
1 1 1
(2.51)
n 1 n 1 n 1 2
Ey {E y }( 2.25) Cb, E y H x,inc i , j ,k 1 2
i , j , k1 i , j , k1 i , j , k1 1
(2.52)
The consistency condition for the magnetic field components located 0.5
outside of each face of the total/scattered-field interface are given by analogy.
The Hy and Hz components can be determined from Figure 2.9. Hy (i = io - ½; j
= jo, …, j1; k = ko + ½, …, k1 - ½) is given by:
n 1 2 n 1 2 n
Hy {H y }( 2.22) Db, H y E z ,inc i , j ,k
i0 1 2, j , k i0 1 2, j ,k i0 1 2, j , k 0
(2.53)
n 1 2 n 1 2 n
Hz i0 1 2, j , k
{H z }
i0 1 2, j , k ( 2.23)
Db , H z E y ,inc i , j ,k
i0 1 2, j , k 0
(2.54)
n 1 2 n 1 2 n
Hy {H y } Db, H y Ez ,inc i , j , k
i1 1 2, j , k i1 1 2, j , k ( 2.22 ) i1 1 2, j , k 1
(2.55)
n 1 2 n 1 2 n
Hz i1 1 2, j , k
{H z i 1 2, j ,k }( 2.23) Db , H z E y ,inc
i1 1 2, j , k i1 , j , k
1
(2.56)
n 1 2 n 1 2 n
H x i, j {H x i , j }
1 2, k ( 2.21)
Db , H x E z ,inc i , j
0 1 2, k i , j0 1 2, k 0 ,k
0
(2.57)
FDTD Technique for Field Truncation 37
n 1 2 n 1 2 n
Hz i , j0 1 2, k
{H z }
i , j0 1 2, k ( 2.23)
Db , H z E x ,inc i , j
i , j0 1 2, k 0 ,k
(2.58)
n 1 2 n 1 2 n
H x i , j 1 2,k {H x i , j 1 2,k }( 2.21) Db , H x E z ,inc i , j ,k
i , j1 1 2, k
1 1 1
(2.59)
n 1 2 n 1 2 n
Hz i , j1 1 2, k
{H z i , j 1 2,k }( 2.23) Db , H z E x ,inc i , j ,k
i , j1 1 2, k
1 1
(2.60)
n 1 2 n 1 2 n
H x i , j ,k {H x i , j , k } ( 2.21) Db , H x E y ,inc i , j ,k
0 1 2 0 1 2 i , j , k 0 1 2 0
(2.61)
n 1 2 n 1 2 n
Hy {H y } Db, H y Ex ,inc i , j ,k
i , j , k 0 1 2 i , j , k0 1 2 ( 2.22 ) i , j , k0 1 2 0
(2.62)
n 1 2 n 1 2 n
H x i , j , k 1 2 {H x i , j , k 1 2 }( 2.21) Db, H x E y ,inc
i , j , k1 1 2 i , j , k1
1 1
(2.63)
In 3-D, all six Cartesian field vector components are separated to realize
the following twelve modified version of Maxwell’s equations [22]:
FDTD Technique for Field Truncation 39
H xy ( E zx E zy )
o *y H xy
t y (2.65)
H xz ( E yx E yz )
o z* H xz
t z (2.66)
H yz ( E xy E xz )
o z* H yz
t z (2.67)
H yx ( E zx E zy )
o x* H yx
t x
(2.68)
H zx ( E yx E yz )
o x* H zx
t x (2.69)
H zy ( E xy E xz )
o *y H zy
t y
(2.70)
E xy ( H zx H zy )
o y E xy
t y (2.71)
E xz ( H yx H yz )
o z E xz
t z (2.72)
E yz ( H xy H xz )
o z E yz
t z (2.73)
E yx ( H zx H zy )
o x E yx
t x (2.74)
40 K. N. Ramli, R. A. Abd-Alhameed and P. S. Excell
E zx ( H yx H yz )
o x E zx
t x (2.75)
E zy ( H xy H xz )
o y E zy
t y (2.76)
PML matching conditions and lattice structures are basically similar to the
TM case.
The standard time-stepping of the Yee algorithm cannot be used since the
attenuation of the outgoing source waves in the PML layer is very fast.
Consequently, an explicit exponential time step algorithm is suggested in the
PML layer [27]. From equations (2.70) and (2.71), the updating Hzy and Exy in
the PML medium can both be expanded by the expressions:
*y
*y i , j ,k t
i , j ,k t i , j ,k E n n
xy i , j 1 2, k E xy i, j 1 2,k
n 1 2 i , j ,k
n 1 2 1 e
H zy e
H zy
i, j ,k i, j ,k *y i, j , k p E xz n E
n
i , j 1 2, k xz i , j 1 2, k
(2.77)
y i , j ,k t
y i , j ,k t
i , j ,k H zx n 1 2 H zx n 1 2
n 1 i , j ,k n 1 e i , j 1 2,k i , j 1 2, k
E xy e
E xy
i , j ,k i , j ,k y i , j ,k p H zy
n 1 2
H zy
n 1 2
i , j 1 2, k i , j 1 2, k
(2.78)
oc ln g
o g
ln R0g
p p
2 p g N 1 (2.79)
FDTD Technique for Field Truncation 41
CONCLUSION
The fundamental concepts of the FDTD computational method were
briefly explained in this chapter. Some of them deal with the solution of
Maxwell's equations by means of finite differences, the derivation of FDTD
updating equations for electric and magnetic field vector components, and the
explanation on the accuracy and stability which necessitate the FDTD
computational method. In addition, the excitation source inside the FDTD
lattice was also explained. The general idea of the ABC to the specific
Berenger’s PML in 2-D and 3-D were well treated. Last but not least, Fortran
source code was written successfully in 2-D FDTD with subgridding
(SGFDTD) and 3-D MoM-FDTD-SGFDTD hybrid technique for different
applications in electromagnetic scattering problems.
REFERENCES
[1] W. G. Whittow, C. C. Njoku and Y.C. Vardaxoglou, “Fine scale
simulations of patch antennas with heterogeneous substrates,” USNC-
URSI Radio Science Meeting (Joint with AP-S Symposium), pp. 223,
2013.
[2] R. Xiong, B. Chen, Y. Yi, H. -L. Chen and H. -F. Zhao, “Two-step
method for the FDTD analysis of long apertures having depth,” IEEE
Antennas and Wireless Propagation Letters, vol. 12, pp. 39-42, 2013.
[3] S. K. Kamepally, B. P. Kumar and C. S. Paidimarry, “FDTD estimation
for accurate specific absorption rate in a tumor,” Annual International
Conference on Emerging Research Areas and International Conference
on Microelectronics, Communications and Renewable Energy
(AICERA/ICMiCR), pp. 1-5, 2013.
[4] C. Guiffaut, A. Reineix and B. Pecqueux, “New oblique thin wire
formalism in the FDTD method with multiwire junctions,” IEEE
42 K. N. Ramli, R. A. Abd-Alhameed and P. S. Excell
[28] J. P. Berenger, “Perfectly matched layer for the FDTD solution of wave-
structure interaction problems,” IEEE Transactions on Antennas and
Propagation, vol. 44, pp. 110-117, 1996.
[29] J. P. Berenger, “Improved PML for the FDTD solution of wave-structure
interaction problems,” IEEE Transactions on Antennas and
Propagation, vol. 45, pp. 466-473, 1997.
In: Development of Complex Electromagnetic … ISBN: 978-1-61122-013-1
Editors: K. N. Ramli et al. © 2014 Nova Science Publishers, Inc.
Chapter 3
K. N. Ramli1, M. Lashab2,
K. H. Sayidmarie3 and R. A. Abd-Alhameed4
1
Faculty of Electrical and Electronics Engineering,
Universiti Tun Hussein Onn Malaysia,
Parit Raja, Batu Pahat, Johor, Malaysia
2
Electronics Department, University of Skikda, Skikda, Algeria
3
Department of Communication Engineering,
College of Electronic Engineering, University of Mosul, Mosul, Iraq
4
Mobile and Satellite Communications Research Centre,
Bradford University, Bradford, UK
ABSTRACT
The classical methods known as asymptotic or high frequency, such
as physical optics (PO) or geometrical theory of diffraction (GTD), can
present an approximate or almost exact solution for electromagnetic
scattering to some given antenna shape. The numerical methods known as
the full wave lead generally to an exact solution for any given shape of
antenna, these can be listed as finite element method (FEM), finite-
difference time-domain method (FDTD) and moment method (MoM).
The moment method is based on the digitalization of the electric field
integral equation (EFIE) or the magnetic field integral equation (MFIE),
obtained from the vector solution of Maxwell equation and for given
46 K. N. Ramli, M. Lashab, K. H. Sayidmarie et al.
3.1. INTRODUCTION
Great research work has been carried out to improve the numerical method
applied to electromagnetic scattering, especially for very complicated
geometries or for stratified and inhomogeneous medium. The moment method
is still one the most powerful method, which could give very exact solution for
many antenna structures [1-3], but for some cases where the geometry of the
antenna is complicated or its medium is inhomogeneous, in this case a hybrid
solution is needed [4, 5]. The choice of the numerical solution depends on the
antenna shape and the type of medium where the antenna is placed. The
Surface Kernel Solution of the Method of Moments 47
moment method (MoM) is used for very simple structures, where surface
integral equations are used [6, 7].
The finite element method (FEM) and the finite-difference time-domain
(FDTD) are used for very complicated geometries where volumetric integral
equations are used.
The computing time is usually O(N2) for the moment method and O(N3)
for the other numerical methods, where N is the number of unknowns. For
very large structures (over 100 wavelength), the computing time and the
memory storage became very big, a hybrid with approximate methods such as
physical optics (PO) is desirable [8, 9].
The moment method was first introduced by Harrington in 1968 [10],
since then many effort were made to improve both the precision and the
computational time spent. A good choice of basis and testing function can
improve the precision and the computing time, for instance sinusoidal
functions known as entire domain, when used good results are obtained but the
computing time is a bit longer than ordinary functions [11-15].
Several research work have been carried out to reduce the computing time
of the classical moment method, such as the fast multipole (FMPM) [16] and
improved moment method based on wavelets [17-21].
These two methods have the ability to reduce the computing time by
simply reducing the impedance matrix coefficients.
.E
(3.3)
48 K. N. Ramli, M. Lashab, K. H. Sayidmarie et al.
.H 0 (3.4)
We define E and H respectively the electric and magnetic field, J is
the current density and the electric charge density.
The medium is considered homogenous and isotropic as shown in Figure
3.1, the permittivity and the permeability are considered constant with
respect to curl and divergence action.
Figure 3.1. Homogenous and isotropic structure fed by electric or magnetic field.
where, Einc is the incident electric field and Escat is the scattered electric
field. Using Maxwell’s equations previously given, the scattered electric field
can be expressed by the following equation.
Escat j. A
(3.6)
where A is the magnetic potential vector, and is the electric potential
scalar. They are obtained after substitution of equation (3.6) into Maxwell’s
equations also by using the unit vector into the expression gauges de Lorentz
[24], given by:
A . J (r ' ).G(r , r ' ).dS '
S (3.7)
1
(r ' ).G(r, r ' ).dS '
S
(3.8)
where G(r,r') is Green’s function in the space domain, and is the charge
density. They are expressed by the following relations:
jk r r '
e
G (r , r ' )
4 . r r '
(3.9)
1
..J
j. (3.10)
n̂ is the unit normal vector to the surface of the structure. Finally, the integral
equation of the electric field can be given as:
j
[k 2
J (r ' ).G(r , r ' ) ' J (r ' ).G(r, r ' ).]dS '
Einc (r ) S
0 (3.12)
where and are respectively the gradient with respect to source point
(x,y,z) and the divergence with respect to the observation point (x,y,z).
Equation (3.12) is defined as the electric field integral equation (EFIE). The
temporal dependence ejt is omitted here, but it is considered present in all this
manuscript. The solution of the integral equation leads to the density of current
in which the directivity and the gain of the antennas are obtained.
J s (r )
nˆ H inc (r ) nˆ J s (r ' ) ' G(r , r ' ).dS '
2 S
(3.14)
Here the integral is called principal value, where the case (r = r’) is
excluded. This type of integral is called Cauchy integral [25].
Surface Kernel Solution of the Method of Moments 51
The moment method (MoM) is very well known technique for resolving
linear equations [4]. For antenna analysis, the moment method converts the
integral equation of the field into a matrix equation, by digitalizing this later,
the matrix equation is given as [3, 22]:
g L( f ) (3.15)
where g is the excitation source, L is a linear integro-differential operator,
and f is the unknown function, which is written as the sum of N independent
functions f n , called the basis functions with unknowns amplitudes an .
N
f an f n
n 1 (3.16)
By injecting equations (3.16) into (3.15) and owing to the fact that L is a
linear operator, equation (3.15) can be written as:
N
g an L( f n )
n 1 (3.17)
N
wm , g an wm , L ( f n )
n 1 (3.18)
52 K. N. Ramli, M. Lashab, K. H. Sayidmarie et al.
g m Z mn . an (3.19)
Equation (3.19) often written in a similar manner to the Ohm’s law, using
the matrix so called impedance, defined by the following scalar product:
Z mn wm , L ( f n )
(3.20)
w1, L( f1 ) w1, L( f 2 )
Z mn w2 , L( f1 w2 , L( f 2
(3.21)
Also the unknown amplitudes and the scalar product of the test function
and excitation are given as:
a1 w1 , g
an a2 gm w2 , g
; (3.22)
The exact solution depends on the choice of the basis function and the test
function. In this case, these have to be linear and independent.
Surface Kernel Solution of the Method of Moments 53
Since the path is 1-D, the integral equations of the electric field (TM case)
and magnetic (TE case) respectively equation (3.12) and (3.14) are written in
the following form.
a/2
EZinc ( x)
4
a / 2
J ( x).H 0( 2) (k. x x' ).dx
(3.24)
54 K. N. Ramli, M. Lashab, K. H. Sayidmarie et al.
J ( x) k a/2
H Zinc ( x) j J ( x' ).H1( 2) (k x x' ).dx'
2 4 a / 2 (3.25)
( 2)
where H (2)
0 and H 1 are Hankel’s functions zero and first order of the
second type respectively, defined as given bellow:
2 j ( kr 3 / 4 )
H 1( 2 ) e
kr (3.27)
The integral equations (3.24) and (3.25) are solved by applying the
moment method as illustrated by equation (3.18). Equations given below are
respectively for the TM and TE case.
N
a
x n 1
n 1
n ( x xm ),
4
xn
( x xn ).H 0( 2 ) ( k x x ' ).dx '
( x xm ), EZinc ( x )
(3.28)
N
m, n k xn1
a
n 1
n ( x xm ), (
2
j
4 xn
( x xn ).H1( 2) (k x x' ).dx'
( x xm ), H Zinc ( x)
(3.29)
1 m n
m, n
0 m n
x xm 1
, x xx
x
m 1 m
x xm 1
( x xm ) , x xx
x
m m 1
0, elsewhere
(3.30)
N
a
xn
n [ ( x xn 1 ).H 0( 2 ) (k xm x ).dx
n 1 4. x x n 1
x n 1
( xn 1 x ).H 0( 2 ) ( k xm x ).dx ] E Zinc ( xm )
xn
(3.31)
N
k
a (
xn
n m,n j [ ( x xn 1 ).H 1( 2 ) (k xm x ).dx
n 1 4. x x n1
56 K. N. Ramli, M. Lashab, K. H. Sayidmarie et al.
x n 1
( xn 1 x ).H 1( 2 ) ( k xm x ).dx ].) H Zinc ( xm )
xn
(3.32)
The matrix elements are written in explicit form and given below by
polarization case TM and TE respectively.
x
(xxn1).H0(2) (k xm x)dx (xn1 x).H0(2) (k xm x).dx
xn1
4.
,n
n
ZmTM
xn1 xn
x (3.33)
ZmTE,n m,n
k xn
(x xx1).H1(2) ( xm x ).dx (xn1 x).H1(2) (k xm x ).dx
xn1
j
4.x xn1 xn
(3.34)
The integral equations (3.33) and (3.34) are written in Matlab code taking
into account the presence of singularity when xn xm . The singularity is due
to discontinuity of Green function [28, 29], in all this chapter we have been
trying to avoid this singularity using the technique of contournement related to
the choice the integration path. The excitation as shown in Figure 3.2,
projected in the x-axis leading to following equations:
depends on geometry of the object and angle of excitation. The RCS is defined
in the literature by the following equation:
2
Ereflechis
( ) lim r 2r 2
Eincident
(3.37)
2
By taking the magnitude of the incident electric field Eincident 1 , and
2 j jk .r
H0( 2) (kr) e , for far field r , is the observation angle.
.kr
By substituting in equation (3.31) and (3.32) respectively for the TM and TE
case, the following relations are obtained [22]:
2
k 2 N
TM
( ) x ane jk . xm cos( )
4 i 1
(3.38)
16
Current Density Jz(mA)
14
12
10
4
-0.5 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.5
Length of the Patch (L/Lmax)
Figure 3.4. Current density of the path, with length 2λ and 4λ, incident angle φi = 90°,
TM polarization, F = 0.3 GHz.
58 K. N. Ramli, M. Lashab, K. H. Sayidmarie et al.
4.5
3.5
2.5
-0.5 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.5
Length of the Patch (L/Lmax)
Figure 3.5. Current density of the path, with length 2λ and 4λ,incident angle φi= 90°,
TE polarization, F = 0.3 GHz.
2
k 2 N
( )
TE
x an .cos( ).e jk . xm cos( )
4 i 1
(3.39)
where is the angle between the vector observation r’ and the vector source
r. The bi-static RCS is represented in decibel using the following formula:
dB 10 . log 10 ( ( )) (3.40)
Figure 3.6 and Figure 3.7 show the bi-static RCS respectively for TM and
TE case. It is remarkable that the RCS is more important for a perpendicular
direction to the patch, and is minimal for parallel direction to the patch. It is
also minimal for longer patch with respect to the wavelength. The results are
in good agreement with the work presented in [25].
Surface Kernel Solution of the Method of Moments 59
EiZinc (r )
4 C
J z (r ' ).H 0( 2) (k r r ' ).dC '
(3.41)
20
Bistatic RCS (dB)
15
10
-5
0 20 40 60 80 100 120 140 160 180
Observation Angle (degres)
Figure 3.6. Bi-static RCS of the patch, length 2λ and 4λ, incident angle φi = 90°, TM
polarization, F = 0.3 GHz.
60 K. N. Ramli, M. Lashab, K. H. Sayidmarie et al.
15
10
Bistatic RCS (dB)
-5
-10
-15
-20
-25
0 20 40 60 80 100 120 140 160 180
Observation angle (degres)
Figure 3.7. Bi-static RCS of the patch, length 2λ and 4λ, incident angle φi = 90°, TE
polarization, F = 0.3 GHz.
J c (r ) k
H Zinc (r ) j J c (r ' ).H1( 2) (k r r ' ).cos( r , r ' ).dC '
2 4 C (3.42)
Surface Kernel Solution of the Method of Moments 61
The index C insinuates that the integral should be done on the contour of
the cylinder, also r , r ' is the angle between the vector (r r ' ) and the normal
vector to the surface n̂ . It is noteworthy that the current Jc, is a tangential
current to the contour and not in the z-axis. Equations (3.41) and (3.42) can be
written by just following the contour C, with changing variables from
Cartesian coordinate to polar coordinate. The integral equations related to
circular cylinder are given as:
2
EZinc (r ( ))
4 0
J Z ( ' ).H0( 2) (k r ( ) r ' ( ' ) ). r ' ( ' ) .d '
(3.43)
J ( ) k 2
H Zinc ( ) j J ( ' ).H1( 2 ) (k r r ' ) cos( r ,r ' ). r ' ( ' ) .d '
2 4 0 (3.44)
The structure of Figure 3.8 presents an axial symmetry, by the fact that the
2-D problem is brought to a 1-D problem, where the only variable is the angle
theta (θ). The integral equations (3.43) and (3.44) are digitalized by the
moment method using equation (3.18). The following equations are
respectively for the cases TM and TE.
N
n1
a
n1
n ( m ),
4 n
( n ).H 0( 2) (k R ). r ' ( ' ) .d '
( m ), EZinc (r ) (3.45)
N
m,n k n 1
a
n 1
n ( m ), (
2
j
4 n
( n ). cos( r , r ' ).H 1( 2 ) ( k R ). r ' ( ' ) .d '
( m ), H Zinc ( r )
(3.46)
with,
2
R r ( ) r ' ( ' ) , n (n 1). , N
and
62 K. N. Ramli, M. Lashab, K. H. Sayidmarie et al.
1
m [(m 1) ]
2
The test functions are Dirac delta and the basis functions are rectangular
type ( m ) , defined as follow:
(3.47)
N n 1
Ei ( r ( m ))
4 i 1
an H 0( 2) (k r ( m ) r ' ( ' ) ). r ' ( ' ) .d '
n
(3.48)
N
m .n k n1
H Zinc ( m ) a n
2
j cos( r , r ' ).H1( 2 ) (kR ). r ' ( ' ) .d '
4 n
n 1 (3.49)
n1
Z mTM,n
4 n
H 0( 2) (k r ( m ) r ' ( ' ) ). r ' ( ' ) .d '
(3.50)
k n1
Z mTE, n m.n j cos( r , r ' ).H1( 2 ) ( kR). r ' ( ' ) .d '
2 4 n (3.51)
The radiation pattern given in Figure 3.14 shows the total pattern in both
cases. However, it is clear that the maximum gain for TM case is facing to the
feeding point, whereas the maximum gain for the TE case is in the opposite
side of the feeding.
0
0 50 100 150 200 250 300 350
Observation Angle(degres)
Figure 3.10. Current density of circular cylinder, cylinder diameter = 1.5λ, incident
angle φi = 0°, TM polarization, F = 0.3 GHz.
0
0 50 100 150 200 250 300 350
Observation Angle (degres)
Figure 3.11. Current density of the circular cylinder, cylinder diameter = 2λ, incident
angle φi = 0°, TE polarization, F = 0.3 GHz.
Surface Kernel Solution of the Method of Moments 65
4
Bi-static RCS (dB)
-2
-4
-6
0 50 100 150 200 250 300 350
Observation Angle (Degrees)
Figure 3.12. Bi-static RCS of the circular cylinder, cylinder diameter=1.5λ, incident
angle φi = 0°, TM polarization, F = 0.3 GHz.
15
Bi-static RCS (dB)
10
-5
-10
0 50 100 150 200 250 300 350
Observation Angle (Degrees)
Figure 3.13. Bi-static RCS of the circular cylinder, cylinder diameter=2λ, incident
angle φi = 0°, TE polarization, F = 0.3 GHz.
66 K. N. Ramli, M. Lashab, K. H. Sayidmarie et al.
a b
Figure 3.14. Radiation pattern for the circular cylinder, (a) TM case (b) TE case.
For simplicity, here in this section the moment method is applied for 3-D
structures having symmetry of revolution known as bodies of revolution
(BOR), as shown in Figure 3.15. Conical horn, parabolic reflector, circular
waveguide and many others are examples of BOR structures. The solution of
the integral equation of such structures can be brought to 2-D cases even
though the structure is 3-D [27, 31]. The integral equation for 3-D structure
and body of revolution, solved by the moment method is given below as the
inner product [29]:
j
Einc ,W
0 S S '
[k 2 (W .J ).G (r , r ' ) ( ' J )(.W )G (r , r ' )]dS '.dS
(3.54)
where W and J are respectively testing and basis functions as vector form,
written in the tangent system of coordinate. The total current on the surface of
the structure is written as [29]:
J (t, ) J (t , ).e
jn
J (t , ) t
n (3.55)
where, J t and J are the two current density components, associated to the
tangential unity vectors. The term e jn corresponds to the Fourier series
expansion with n is the mode number.
= ×( − ) (3.56)
=− ×( − ) (3.57)
68 K. N. Ramli, M. Lashab, K. H. Sayidmarie et al.
The fields inside the volume can be set to arbitrary values, since the
equivalent currents are realized for any discontinuity. The internal fields are
set to zero for simplicity of the problem, = 0 and = 0.
a b
Figure 3.16. Surface equivalence principle, (a) Original problem (b) Replaced by
surface of electric current and magnetic currents.
1
K ( J (r )) J (r ) nˆ J (r ' ) ' G(r , r ' ).ds' nˆ H i (r )
2 S
(3.58)
2 0 1 j 2 m1
J t (t , ) a .
n0
t
n
t
j ,n (t , ) cmt , n mt , n (t , )
m0 n0 (3.59)
2 0 1 j 2 m1
J (t , ) an . j , n (t , ) cm , n m , n (t , )
n 0 m0 n 0 (3.60)
where ( m,n , m,n ) and ( j ,n , j ,n ) are respectively, the mother wavelet and
t t
j 2 m 1
J t (t , ) cmt , n mt , n (t , )
m0 n0 (3.61)
j 2 m 1
J (t , ) cm , n m , n (t , )
m 0 n 0 (3.62)
W , K ( J ( r )) W , nˆ H i (r )
(3.63)
1 t t
2 sW . J ( r ).ds
s s
[( ˆ
n W ) J (r ' )].G (r , r ' )ds '.ds
W t .(nˆ H i )ds
s (3.64)
1
2 sW .J (r ).ds s s[(nˆ W ) J (r ' )].G (r , r ' )ds'.ds
W .(nˆ H i )ds
s (3.65)
Knowing that,
t
nˆ W W tˆ , and nˆ W W ̂
t
Also,
(W t nˆ) J (W tˆ ) ( Jt tˆ' Jˆ ' )
(3.66)
(W nˆ) J (W tˆ) ( Jt tˆ' Jˆ ' )
(3.67)
n t
n H i
1 0 0 Hi tˆ H tiˆ
0 H ti Hi
(3.68)
Surface Kernel Solution of the Method of Moments 71
The inner product in equation (3.64) and (3.65) can be expressed as:
W t .J W t tˆ.( J t tˆ' J ˆ ' ) W t J t (3.69)
W .J W tˆ.( J t tˆ' J ˆ ' ) W J (3.70)
1
G(r, r ' ) Rˆ .( jk ).G(r , r ' )
R (3.71)
where
Rˆ R / R
,
'
R r r ' ( ' ) 2 ( z z ' ) 2 4 ' sin 2 ( )
2
By replacing equations (3.69) and (3.70) into equations (3.64) and (3.65),
and adding the summation in (3.61) and (3.62):
1 t t
T T T
Wq J p . .dt (Wqt J tpˆ tˆ'Wqt J p ˆ ˆ ' ).I G1. ' .dt '.dt
p
0 2 0 0
T
Wqt H i tˆ.tˆ.I G 2 .dt
0 (3.72)
T 1 T T
p
0 2
Wq J p . .dt (Wq J tp tˆ tˆ'Wq J p tˆ ˆ ' ).I G1. ' .dt '.dt
0 0
T
Wq H tiˆ .ˆ .I G 2 .dt
0 (3.73)
72 K. N. Ramli, M. Lashab, K. H. Sayidmarie et al.
2
I G1 G ( r , r ' ).e j . 'd '
0 (3.74)
1 2
IG 2
2 0
e j d
(3.75)
where the index is the mode number and T is the maximum arc length of the
antenna (for a sphere T .r , r is the radius). The indices (p,q) correspond
respectively to the combinations (m p , n p ) and ( mq , nq ) . Equations (3.72)
and (3.73) can be written as a matrix equation as [37]:
Z tt Z t cmt , n H 1
.t .
Z Z cm , n H 2
(3.76)
Each element of the matrix (3.76) is given based on equations (3.72) and
(3.73):
1 t t
Z ttpq .Wq J p .dt W J q
t t
p .ˆ tˆ'.I G1. ' dt ' dt
t 2 t t'
(3.77)
1
Z ttpq p , q , T (t , t ).(t , )
2 (3.78)
where T(t,t) is the second term under the double integral, (t,, ) is an
operator of changing variables from(t, t ' ) to ( , ' ) , is a variable related
the wavelet belonging to domain [0,1], finally p and q are wavelets of
different rank. Equation (3.78) can be written in more explicit form as
introduced by Wang [38] and Tretiakov [39].
Surface Kernel Solution of the Method of Moments 73
1 1 1
Z ttpq q [ p pT (t , t ) (t , ) d .D( )]D( ' )d '
0 2 0
(3.79)
Z tpq W q
t
J p .ˆ ˆ '.I G1. ' dt ' dt
t t' (3.80)
Z tpq q , p , T (t , ).(t ,
(3.81)
and,
Z pq.t W
q J tp .tˆ tˆ'.I G1. ' dt ' dt
t t' (3.82)
Z pq.t q , p , T ( , t ).(t ,
(3.83)
Also,
1
pq
Z .Wq J p .dt W
q J p .tˆ ˆ '.IG1. ' dt ' dt
t 2 t t'
(3.84)
1
pq p , q ,
Z T ( , ).(t , )
2 (3.85)
H q1 q , H t I G 2 .(t , )
(3.86)
H q2 q , H I G 2 .(t , )
(3.87)
74 K. N. Ramli, M. Lashab, K. H. Sayidmarie et al.
Once the matrix equation (3.76) is solved the unknowns [cmt ,n , cm , n ] which
are the components of the current density are obtained. To find the co-
polarization and cross-polarization, a transformation from tangential
coordinate to polar coordinate have to be done as in [40]. The wavelet
employed is constructed from Haar orthogonal wavelet with vanishing
moment N=7. The lowest resolution level is chosen as 2j = 27 = 128. Since 128
wavelets are involved, a system with 128 128 matrix elements is
generated.The radiation pattern given in Figure 3.18 and Figure 3.19 are
obtained for given ratio focal distance over length of diameter, with F/D = 0.5,
the diameter = 10λ, the incident angle φ=45° and the frequency F = 30 GHz.
The figures show a gain of almost 10 dB is obtained, and the co-polarization is
maximum at the front side with its very narrow main lobe (20° width).
In reverse, the cross polarization is null at the front of the antenna and
maximum at its sides. The radiation pattern given in Figure 3.20 and Figure
3.21 show the rectangular and polar radiation pattern respectively. with F/D =
0.5, the incident angle of the excitation is 45°, and the diameter = 5λ and the
frequency F = 30 GHz. The figures show a gain of almost 15 dB is obtained
with the main lobe of the co-polarization is larger than last one (40° width).
The obtained results were compared to previously obtained results in the
literature [16, 41] and they were in good agreement.
-10
Radiation Pattern (dB)
-20
-30
-40
-50
-60
-70
-150 -100 -50 0 50 100 150
Observation Angle - (degres)
Figure 3.18. Rectangular radiation pattern of the reflector antenna, D=10λ at frequency
F=30 GHz.
Surface Kernel Solution of the Method of Moments 75
Figure 3.19. Polar radiation pattern of the reflector antenna, D=10λ at frequency F=30
GHz, co-polarization.
10 Copolarization
CrossPolarization
0
Radiation Pattern (dB)
-10
-20
-30
-40
-50
-60
Figure 3.20. Rectangular radiation pattern of the reflector antenna, D=5λ at frequency
F=30 GHz.
76 K. N. Ramli, M. Lashab, K. H. Sayidmarie et al.
Figure 3.21. Polar radiation pattern of the reflector antenna, D=5λ at frequency F=30
GHz, co-polarization.
assume a uniform surface current distribution across the cross section, (e.g.,
Djordjevic et al. [50], Burke and Poggio [51] and Richmond [52]).
Hence, surface resistive losses and reactive effects that may be augmented
by the non-uniform surface current will not be correctly predicted.
This problem is particularly significant for resonant coiled electrically-
small antennas, such as the normal-mode helical antenna (NMHA), in which
the surface current distribution has a critical effect on the efficiency and Q-
factor. A moment-method (MoM) formulation uses two orthogonal basis
functions on the surface of the wire including its ends was thus developed to
investigate this problem in details. The work detailed a more generalised
theory and results of the work done by the present authors [47]. It should be
noted that there are more advanced commercial codes now available which, in
particular, implement patch modelling more effectively (e.g., FEKO [53], CST
[54], HFSS [55] and IE3 [56]). In general MoM [57-63] has the ability to solve
complex antenna pattern in frequency and time domain.
Basically, the original motivation for this work was to assess the degree of
benefit that would be obtained if an antenna of this type were to be realised in
high-temperature superconductor. Electrically-small antennas have a low
radiation resistance that is easily swamped by ohmic loss resistance, resulting
in a low efficiency. Superconductors have the potential to remove much of the
loss and hence raise the efficiency significantly. There is then the possible
disadvantage that the inherent Q-factor of the antenna may become very high:
whether this is a real disadvantage depends on the nature of the system into
which the antenna is proposed for deployment. To quantify the reduction in
loss, and hence improvement in efficiency, which might accrue from the use of
superconductor, it is necessary to quantify the surface loss Ps measured in
W/m2:
J s2
Ps
s (3.88)
where Js is the surface current density (A/m) and s is the surface conductivity
(). The self-resonant helix had already been identified as a convenient design
of electrically-small antenna in which quite interesting results were reported
for example, broad band V-helical antenna [64], circular NMHA [65], double
pitch NMHA [66] and multiple pitches NMHA [67], however, for realisation
in high-temperature superconductor, the superconducting element may be left
electrically isolated. The detailed quantification of Js in this particularly
78 K. N. Ramli, M. Lashab, K. H. Sayidmarie et al.
complex case was thus the main original objective of the work. The very
detailed modelling procedure that has been developed has much wider uses,
particularly in the accurate modelling of normally-conducting NMHAs, which
see extensive use in mobile telecommunications.
Complete validation of the predictions of the procedure poses
considerable difficulties, since it would require measurement of the surface
current distribution on wire. This matter is an important topic for future work,
but an adequate degree of validation can be claimed for the results that have
been presented in this work from this type of modelling process.
I j
j .L ( J j ) (E i )
(3.89)
where Ij is a basis function for the surface current Jj, Ei is the incident electric
field strength and L is the integro-differential operator given by:
L ( J ) ( j A ) tan (3.90)
where A and are the vector and scalar potentials. If a set of testing functions
Wm is defined, equation (3.89) may be rewritten as:
I j
j Wm , L ( J j ) Wm , E i for j 1, 2,... N
(3.91)
where
Surface Kernel Solution of the Method of Moments 79
Wm , L ( J j ) ( W
s s'
m .L ( J j )) ds ' ds Z mj
(3.92)
Wm , E i ( W
s'
m .E i ) ds Vm
(3.93)
where ds and ds are the differential areas on the wire surface for the source
and the observation points respectively, m = 1, 2, ... N is the index of the
testing function and Z and V are the conventional abbreviations for the
interaction matrix and excitation vector terms in the Method of Moments.
1
Z mj j J .W
s s'
j m
k2
(.J j )(.Wm ) g ( R)ds' ds
(3.94)
where g(R) is the free-space Green function [69-71] and is given by the
expression:
e jkR
g (R)
4R (3.95)
R is the distance between the observation and source points on the wire
surface. Singular integral occurs when R = 0 (i.e., g(R) ). The coordinates
for a point on the surface of the helix wire can be given by:
P
z , a sin cos
2 (3.98)
where:
P
tan
2b (3.101)
where a is the radius of the helix wire, b is the radius of the helix, P is the
pitch distance between the turns, is the azimuth angle of the circumferential
cross-section wire and is the pitch angle. Equations (3.99), (3.100) and
(3.101) are the exact coordinates of helix geometry. Defining two orthogonal
directions on the surface of the helix wire as shown in Figure 3.22, the unit
vectors of the curvilinear surface patches in both directions are:
Figure 3.22. Basic geometry of the helical antenna driven by a voltage source at its
centre. The directions of the orthogonal basis or test functions are shown on the right
and represent the source or observations points on the wire surface and its ends.
Surface Kernel Solution of the Method of Moments 81
aˆ sin cos aˆx cos cos aˆy sin aˆz (3.102)
aˆ cs sin cos cos sin sin aˆ x
sin sin cos cos sin aˆ y cos cos aˆ z
(3.103)
where â and âcs are the unit vectors in the axial and circumferential
surfaces of the wire respectively as shown in Figure 3.22. The differential
lengths in both directions are:
dcs ad (3.105)
2
P '2
b'' b
2 (3.106)
2
b ' (b a cos( )) 2 (3.107)
x r , b r cos
(3.108)
where 0 ra. Hence, the unit direction vectors of the basis function on the
end surface can be expressed as:
82 K. N. Ramli, M. Lashab, K. H. Sayidmarie et al.
where â r and â ce are the unit vectors in the radial and the circumferential
directions on the end surface of the wire respectively as shown in Figure 3.22.
The differential area on the end surface can be given by:
f
f for 0 and 1 2
f 1
(3.114)
f ' f ' for 0 and 1 2
(3.115)
matrix elements can be found. As an example, the impedance element for basis
and test function in the axial direction can be stated as follows:
1
Z ' j aˆ .aˆ f ( ) f ( ' ) k
s s'
' 2
f ( ) f ( ' ) g ( R)ds' ds
(3.116)
where:
The other self and mutual impedance elements for all other basis
directions can be obtained in a similar way. The magnitude of R is written as:
2
P 2
(b a cos( ))
2
(3.120)
Let say there is another axial length 1 of the curvilinear patch, then:
2
P 2
1 (b a cos( )) 1
2
(3.121)
It can be clearly seen that remains constant along the axial length. In
other words:
84 K. N. Ramli, M. Lashab, K. H. Sayidmarie et al.
f ( )
1 1 (3.122)
Moreover, when the surfaces of the patches are very close in proximity to
each other still there will be pseudo-singularity. The condition occurs when:
r p puˆ p tuˆ t
(3.126)
Surface Kernel Solution of the Method of Moments 85
r1s r1 r p
(3.127)
H( , , ' , ' )e jkR( , , ', ')
Z ' ' j
4R( , , ' , ' )
d ' d ' dd
' '
J ( , , p' , t ' )e jkR( , , p ',t ')
4R( , , p' , t ' )
dp' dt ' dd
p' t '
J ( , , p' , t ' )e jkR( , , p ',t ')
j
p' t '
4R( , , p' , t ' )
dp' dt ' dd
(3.129)
86 K. N. Ramli, M. Lashab, K. H. Sayidmarie et al.
k 1l 1
4R(i , j , p'k , t 'l ) t 'l p ' k
jkR (i , j , p ',t ')
n n J (i , j , p ' , t ' )e
j dt ' dp' w j wi
i 1 j 1 p 't ' 4R(i , j , p ' , t ' )
(3.130)
The last term of equation (3.130) can be solved using analytical solution:
For simplicity, the above equation can be reduced to the following format:
a1 a 2 t
I
p t v 2 ( p p1 ) 2 (t t1 ) 2
dtdp
(3.132)
where a1, a2, v, p1 and t1 are constants. The analytical solution of the equation
is given by the expression:
distribution over the wire surface and wire ends as for the dipoles are used. A
computer program was written to implement the analysis given in the previous
section. The surface patch subdivision was automatically generated by the
program, subject to the number of the basis functions in both orthogonal
directions. The impressed field Ei is modelled by a delta-function voltage
source at the centre of the dipole and the helix whereas in loop it was placed at
= 0.
a b
Figure 3.24. The geometry models of two parallel dipoles and loop antennas including
the directions of the basis or test functions used; (a) dipoles, (b) loops. A simple axial
excitation (in the -direction) was considered. Thus the impressed field can be given
by:
1
Ei ( c )aˆ
2a (3.134)
where a is the radius of the wire antenna and c is the half axial length of the
dipole or the helix. The antenna wire for all geometries was assumed to be
perfectly conducting and surrounded by free space. Several examples were
used to investigate the surface current distribution of dipoles, loops and
NMHA as predicted by the formulation, as follows.
The response of the input impedance of two parallel dipoles of 50 cm
length and 5 mm wire radius separated by 15 mm is presented in Figure 3.25.
In this example, both dipoles were centrally fed as presented in equation
(3.47).
The axial and circumferential lengths were subdivided by 16 curvilinear
patches of equal lengths in both directions. The attachment basis modes
between the wire ends and the wire surface were placed subject to the
88 K. N. Ramli, M. Lashab, K. H. Sayidmarie et al.
Figure 3.25. The input impedance versus frequency for two parallel dipoles separated
by 15 mm; solid line: NEC [51], (‘ooo’ Resistance, ‘xxx’ Reactance present work).
Surface Kernel Solution of the Method of Moments 89
Figure 3.26. Input impedance at 300 MHz operating frequency of two parallel dipoles
of 50 cm length and 5 mm wire radius versus the separated distance between them;
(‘***’, ‘+++’: Present work), (‘ooo’ and ‘xxx’: NEC).
a b
Figure 3.27. The normalised magnitudes of the axial and circumferential surface
current components of the antenna geometry given in Figure 3.24 separated by 15 mm
versus at = 6.25 cm for axial and = 4.6845 cm for circumferential from the
bottom of the dipoles for different operating frequencies: ‘xxx’ 100 MHz, ‘***’ 300
MHz, ‘ooo’ 500 MHz; (a) axial component, (b) circumferential component.
90 K. N. Ramli, M. Lashab, K. H. Sayidmarie et al.
a b
Figure 3.28. The normalised magnitudes of the axial and circumferential surface
current components versus at similar location as Figure 3.27, for different separated
distances between the dipoles at 300 MHz (equivalent half wavelength dipoles with
0.005 wire radius): (a) axial component, (b) circumferential component. (‘ooo’: 15
mm, ‘xxx’: 30 mm, ‘***’: 50 mm, ‘+++’: 100 mm, ‘□□□’:500 mm).
a b
Figure 3.29. The normalised magnitudes of the axial and circumferential surface
current components versus at similar locations as Figure 3.27, for different separated
distances between the dipoles at 300 MHz (equivalent to half wavelength dipoles with
0.01wire radius): (a) axial component, (b) circumferential component. (‘ooo’: 30 mm,
‘xxx’: 50 mm, ‘***’: 100 mm).
92 K. N. Ramli, M. Lashab, K. H. Sayidmarie et al.
Figure 3.30. Power loss ratio of two parallel half wavelength dipoles for various
separation distances; ‘ooo’: wire radius = 0.005λ, ‘xxx’: wire radius = 0.01.
a b
Figure 3.31. The normalised magnitudes of the axial and circumferential surface
current components versus at = 0: axial component and = 33.75:
circumferential component, for a single loop antenna for different operating
frequencies. The loop radius and wire radius are 3 cm and 5 mm respectively.(a) axial
component, (b) circumferential component. (‘xxx’: 400 MHz, ‘ooo’: 600 MHz, ‘+++’:
800 MHz, ‘***’: 900 MHz).
Surface Kernel Solution of the Method of Moments 93
a b
Figure 3.32. The normalised magnitudes of the axial and circumferential surface
current components versus at = 0: axial component and = 33.75:
circumferential component, for two parallel loops separated by 15 mm and each has a
radius of 3 cm and 5 mm wire radius versus different operating frequencies. (a) axial
component, (b) circumferential component. (‘+++’: 400 MHz, ‘xxx’: 600 MHz, ‘***’:
800 MHz, ‘ooo’: 900 MHz).
a b
Figure 3.33. The normalised magnitudes of the axial and circumferential surface
current components versus at = 0: axial component and = 33.75:
circumferential component, for two half wavelength parallel loops each of radius
0.0796 wavelength and 0.013 wavelength wire radius, versus the separation distance
‘d’. (a) axial component, (b) circumferential component. (‘+++’: d = 3a, ‘xxx’: d = 6a,
‘***’: d = 20a, ‘ooo’: d = 100a).
a b
Figure 3.34. The normalised magnitudes of the axial and circumferential surface
current components of a half wavelength single turn helix antenna versus at different
positions from the first end of the helix. The helix radius, wire radius and pitch
distance are 0.0796, 0.013 and 3a respectively. (a) axial component: ‘***’:0.031,
‘ooo’:0.125, ‘+++’:0.25 (b) circumferential component: ‘***’:0.0156,
‘ooo’:0.109, ‘+++’:0.234.
Surface Kernel Solution of the Method of Moments 95
a b
Figure 3.35. The normalised magnitudes of the axial and circumferential surface
current components of the same antenna geometry given in Figure 3.33, except that the
operating frequency is half than that used in Figure 3.33 (i.e. operating wavelength
2). (a) axial component: ‘***’:0.031, ‘ooo’:0.125, ‘+++’:0.25 (b) circumferential
component: ‘***’:0.0156, ‘ooo’:0.109, ‘+++’:0.234.
a b c
Figure 3.36. The normalised magnitudes of the axial and circumferential surface
current components of a half wavelength two turns helix antenna versus for different
pitch distances. The helix radius and wire radius are 0.04 and 0.006 respectively.
(‘ooo’: P = 3a, ‘+++’: P = 5a, ‘xxx’: P = 7a, ‘***’: P = 9a)(a) axial component: taken
at 0.031 from the bottom end of the helix. (b) axial component: taken at the centre of
the helix. (c) circumferential component: taken at 0.023 from the bottom end of the
helix.
pointed inside the helix for all locations presented. These results might help to
approximate the equivalent of these variations into one curve that might be
taken along all the local length of the helix to assess the total loss power in the
axial direction. Similar variations were observed on the circumferential
component except that the strong effect has shifted from = 90 to = 135
for the locations presented. However, considering the strong effect locations of
these currents regardless of the locations of their maxima and taking into
account the similarities in these variations it may be concluded that this
indicates the approximate contribution of this current to the total loss power.
However, for the same antenna geometry given in Figure 3.34, the
currents were computed at half the operating frequency (i.e. operating
wavelength 2) as shown in Figure 3.35. It is clearly shown that the strong
effects of these current variations are mostly similar to that presented in Figure
3.34. Moreover, a strong correlation can be found in the variations of the
circumferential currents for different locations for this example. However, a
one turn helix is not sufficient to permit comment on the currents variations of
a multi-turn helix antenna, thus the following examples are considered. The
normalised magnitudes of the axial and circumferential surface current
components of a half wavelength of two, three, four and five turns helix
antenna versus for different pitch distances are shown in Figures 3.36, 3.37
and 3.38.
a b c
Figure 3.37. The normalised magnitudes of the axial and circumferential surface
current components of a half wavelength three turns helix antenna versus for
different pitch distances. The helix radius and wire radius are 0.0265 and 0.004
respectively. (‘ooo’: P = 3a, ‘+++’: P = 5a, ‘xxx’: P = 7a, ‘***’: P = 9a); (a) axial
component: taken at 0.0294 from the bottom end of the helix. (b) axial component:
taken at the centre of the helix. (c) circumferential component: taken at 0.022 from
the bottom end of the helix.
Surface Kernel Solution of the Method of Moments 97
a b
Figure 3.38. The normalised magnitudes of the axial surface current components of a
half wavelength four turns helix antenna versus for different pitch distances. (a) four
turns: the helix radius and wire radius are 0.02 and 0.003 respectively; (b) five turns:
the helix radius and wire radius are 0.015 and 0.002 respectively; (‘ooo’: P = 3a,
‘+++’: P = 5a, ‘xxx’: P = 7a, ‘***’: P = 9a).
It is very clear that the axial and circumferential components were non-
uniform even if the pitch distance between the helix turns was nine times the
radius of the wire. Another interesting point is that the peak values of the axial
component were pointed inside the helix (i.e. around = 180) for all helices
presented. This is clearly shown in the variations of these currents at the feed
points on Figures 3.36(b), 3.37(b), 3.38(a) and 3.38(b) and the helix turns in
Figures 3.36(a) (first turn) and 3.37(a) (first turn). The similarities of these
variations permit the approximate calculation of the effective power loss in
that particular direction. It was also observed that the maximum variations of
the circumferential component on the first turn confined between = 0 and
180 as shown in Figures 3.34(b), 3.35(b), 3.36(c) and 3.37(c). However, the
maximum ratio of the axial component to the circumferential component for
all pitch distances was found between 15:1 and 40:1 for all helices more than
one turn.
General comments on the trend of these results are to predict the accurate
or approximated equivalent power losses that are associated with non-uniform
variations along the wire surfaces. Hence this will follow to affect the
radiation efficiency of this kind of antennas.
98 K. N. Ramli, M. Lashab, K. H. Sayidmarie et al.
CONCLUSION
The application of moment method in space domain has been presented in
this chapter, for 1-D, 2-D and 3-D antennas. From the 1-D antenna, we can
understand the very basic principles of the moment method since only one
variable involved. The application of this case can be just for very thin wire
antenna. It has been noticed that the current at the edges is vanishing due to the
current has greater value at the edges. In the second case of 2-D antenna, an
example of infinite cylinder has been studied. This can have an application for
thin wire antenna but infinite, the most important concept presented here are
the electric integral equation (EFIE) and the magnetic integral equation
(MFIE) manipulation. For each type of polarization we have to choose the
corresponding integral equation, also the concept of digitalizing the integral
equation is another important task to understand. In the last section, the 3-D
antennas have been treated, but only for body of revolution (BOR), where a 3-
D problem is brought to 2-D problem.
An application has been presented concerning the reflector antenna fed by
a dipole antenna, working at 30 GHz. The basis and testing functions are
chosen as wavelets mainly to reduce the computation time and the size of the
impedance matrix.
Finally, the radiation pattern of the reflector antenna was given for
different size of the diameter, the co-polarization is maximum in front of the
antenna and the cross-polarization is null in front of the antenna.
The surface current distributions on structures with closely spaced parallel
wires, such as dipoles, loops and helical antennas, can be computed by using
the method of moments with a general surface patch formulation. The current
distribution varies substantially from the common assumption that it is
uniform around the wire cross-section.
Transverse (circumferential) currents are shown to be present: they are
relatively weak on thin wires (around 0.01 wire radius) excited by axial
component parallel to the local axis of the wire. The effect is still significant
when the wire separation distance is relatively large.
In spite the strong variations of the axial and circumferential currents, it
was found that the input impedance and the average value of the axial surface
current are in reasonable good agreement with the results of thin wire codes
such as NEC using an extended kernel solution. The power loss ratio resulting
from use of non-uniform surface current, compared with the conventional
uniform assumption of two parallel dipoles, shows a significant increase of
power loss when they are closely separated. However, these current variations
Surface Kernel Solution of the Method of Moments 99
will dominate the radiation efficiency when predicting the accurate total power
loss on these types of antennas and this can be important in some applications,
e.g., highly resonant antennas and antennas realised in superconducting
materials. As a matter of interest it was computed that the maximum ratio of
the variations of the axial component to the circumferential component on a
half wavelength helix of a few turns for different pitch distances was between
15:1 to 40:1. This behaviour is expected as the NMHA has a hybrid of dipole
and loop behaviour. The modelling method employed a two-dimensional
electric surface patch integral equation formulation solved by independent
piecewise-linear basis function methods in the circumferential and axial
directions of the wire. A similar orthogonal basis function was used on the end
surface and appropriate attachments with the wire surface were employed to
satisfy the requirements of current continuity. The results were stable and
showed good agreement with less comprehensive earlier work by others.
REFERENCES
[1] M. Amir, S. Bedra, S. Benkouda, and T. Fortaki, “Bacterial foraging
optimisation and method of moments for modelling and optimisation of
microstrip antennas,” IET Microwaves, Antennas and Propagation, vol.
8, no. 4, pp. 295-300, 2014.
[2] D. J. Ludick, R. Maaskant, D. B. Davidson, U. Jakobus, R. Mittra, and
D. De Villiers, “Efficient analysis of large aperiodic antenna arrays
using the domain Green's function method,” IEEE Transactions on
Antennas and Propagation, vol. 62, no. 4, part 1, pp. 1579 - 1588, 2014.
[3] S. R. Zang and J. R. Bergmann, “Analysis of omnidirectional dual-
reflector antenna and feeding horn using method of moments,” IEEE
Transactions on Antennas and Propagation, vol. 62, no. 3, pp. 1534-
1538, 2014.
[4] M. N. O. Sadiku, Numerical techniques in electromagnetics, New York:
CRC Press, 2001.
[5] Q. Ye, Electromagnetic scattering by numerical methods for large
structures, PhD thesis, Canada: Manitoba University, 2000.
[6] E. V. Jull, Antenna handbook – Theory, applications and design, Ch1,
London, 1999.
[7] C. P. Davis, Understanding and improving moment method scattering
solutions, M. Sc. Thesis, Birmingham: Birmingham University, 2004.
100 K. N. Ramli, M. Lashab, K. H. Sayidmarie et al.
APPENDIX A
To find the tangent surface equation at a point on the curved surface can
be stated in the following. For simplicity, assume each point on the curved
surface S for all 3-D axes can be represented by two variables as u and v given
by:
x f u , v , x f u , v (A.1)
y g u , v , y g u , v (A.2)
z h u , v , z h u , v (A.3)
a fu gu hu (A.4)
b fv gv hv (A.5)
L
M ab
(A.6)
N
The above notation can be extended using the cross product as in the
following matrix:
L g u hv g v hu
M f h f h
v u u v (A.7)
N f u g v f v hu
D L2 M 2
N2 (A.8)
The tangent plane should satisfy the following formula:
L M N
x y z Lx My Nz (A.9)
D D D
Thus, the tangent surface for a point at (xo, yo, zo) can simply be expressed
by the following:
L x x M y y N z z 0 (A.10)
F x , y , z 0 (A.11)
The above equations are also applicable to direct surface equation such as
that given by the unit sphere at the origin:
x2 y2 z2 1 0 (A.12)
Surface Kernel Solution of the Method of Moments 107
F
x x F y y F z z 0 (A.13)
x x , y , z y x , y , z
z x , y , z
Again the above derivation can be applied to simple surfaces such as the
one given by the following example. Surface equation is z = 4x3y2 + 2y and the
tangent at point (1, -2, 12) can be derived as:
f x , y 4 x 3 y 2 2 y (A.15a)
f
12 x 2 y 2 (A.15b)
x
f
8x3 y 2 (A.15c)
y
f
48 (A.15d)
x x , y
108 K. N. Ramli, M. Lashab, K. H. Sayidmarie et al.
f
14 (A.15e)
y x , y
48 x 1 14 y 2 z 12 0 (A.15f)
48 x 14 y z 64 (A.15g)
Finally, the tangent surface of helix wire can be derived as follows. Given
the coordinates of the helix by the following:
P
z a sin (A.18)
2
x x
, ,
L x
y
M
(A.19)
N , ,
z
x
,
,
Surface Kernel Solution of the Method of Moments 109
APPENDIX B
The expression for C in equation (3.133) in Chapter 3 is given by:
1
1
C log(t1 t (v2 p2 2 pp1 p12 t 2 2tt1 t12 ) 2 ) p p p1
2
1
1 1
log(v2 p2 2 pp1 p12 ) v atan( (2 p 2 p1) / v) t / ((t t1)2 ) 2
2 2
1 1 1 1
log((2t 2 4tt1 2t12 2(v2 ) 2 ( p p1 (v2 ) 2 ) 2((t t1)2 ) 2 (( p p1 (v2 ) 2 )2
1 1 1 1 1 1
1
2(v2 ) 2 ( p p1 (v2 ) 2 ) t 2 2tt1 t12 ) 2 ) / ( p p1 (v2 ) 2 )) p1 t(v2 ) 2 / ((t t1)2 ) 2
2
1 1 1 1
log((2t 2 4tt1 2t12 2(v2 ) 2 ( p p1 (v2 ) 2 ) 2((t t1 )2 ) 2 (( p p1 (v2 ) 2 )2
1 1 1 1 1 1
1
2(v2 ) 2 ( p p1 (v2 ) 2 ) t 2 2tt1 t12 ) 2 ) / ( p p1 (v2 ) 2 )) t(v2 ) 2 / ((t t1)2 ) 2
2
1 1 1 1
log((2t 2 4tt1 2t12 2(v2 ) 2 ( p p1 (v2 ) 2 ) 2((t t1 )2 ) 2 (( p p1 (v2 ) 2 )2
1 1 1 1 1
1
2(v2 ) 2 ( p p1 (v2 ) 2 ) t 2 2tt1 t12 ) 2 ) / ( p p1 (v2 ) 2 )) t / ((t t1 )2 ) 2
2
1 1 1 1
log((2t 2 4tt1 2t12 2(v2 ) 2 ( p p1 (v2 ) 2 ) 2((t t1 )2 ) 2 (( p p1 (v2 ) 2 )2
1 1 1 1 1
1
2(v2 ) 2 ( p p1 (v2 ) 2 ) t 2 2tt1 t12 ) 2 ) / ( p p1 (v2 ) 2 )) p1 t1 / ((t t1 )2 ) 2
2
1 1 1 1
log((2t 2 4tt1 2t12 2(v2 ) 2 ( p p1 (v2 ) 2 ) 2((t t1 )2 ) 2 (( p p1 (v2 ) 2 )2
1 1 1 1 1 1
1
2(v2 ) 2 ( p p1 (v2 ) 2 ) t 2 2tt1 t12 ) 2 ) / ( p p1 (v2 )2 )) p1 t1 (v2 ) 2 / ((t t1 )2 )2
2
1 1 1 1
log((2t 2 4tt1 2t12 2(v2 )2 ( p p1 (v2 ) 2 ) 2((t t1 )2 ) 2 (( p p1 (v2 )2 )2
1 1 1 1
2(v2 ) 2 ( p p1 (v2 ) 2 ) t 2 2tt1 t12 ) 2 ) / ( p p1 (v2 )2 ))t
1
log( p p1 (v2 p2 2 pp1 p12 t 2 2tt1 t12 ) 2 ) t1
1 1 1
1
log( p p1 (v2 p2 2 pp1 p12 t 2 2tt1 t12 ) 2 ) t1(v2 ) 2 / ((t t1)2 ) 2 ...
2
110 K. N. Ramli, M. Lashab, K. H. Sayidmarie et al.
1 1 1 1
... log((2t 2 4tt1 2t12 2(v2 ) 2 ( p p1 (v2 ) 2 ) 2((t t1)2 ) 2 (( p p1 (v2 ) 2 )2
1 1 1 1 1
1
2(v2 ) 2 ( p p1 (v2 ) 2 ) t 2 2tt1 t12 ) 2 ) /( p p1 (v2 ) 2 )) t1 /((t t1)2 ) 2
2
1 1 1 1
log((2t 2 4tt1 2t12 2(v2 ) 2 ( p p1 (v2 ) 2 ) 2((t t1)2 ) 2 (( p p1 (v2 ) 2 )2
1 1 1 1
2(v2 ) 2 ( p p1 (v2 ) 2 ) t 2 2tt1 t12 ) 2 ) /( p p1 (v2 ) 2 )) p1
1
1 1
D (2 p 2 p1 )(v 2 p 2 2 pp1 p12 t 2 2tt1 t12 ) 2
4 2
1
1
log(p p1 (v 2 p 2 2 pp1 p12 t 2 2tt1 t12 ) 2 )v 2
2
1
1
log(p p1 (v 2 p 2 2 pp1 p12 t 2 2tt1 t12 ) 2 )t12
2
1
log(p p1 (v 2 p 2 2 pp1 p12 t 2 2tt1 t12 ) 2 )t 2 t1 p
1 1
1
log(t1 t (v 2 p2 2 pp1 p12 t 2 2tt1 t12 ) 2 ) t1 p t1tp1 /((t t1) 2 ) 2
2
1 1 1 1 1
log((2t 2 4tt1 2t12 2(v 2 ) 2 ( p p1 (v 2 ) 2 ) 2((t t1) 2 ) 2 (( p p1 (v 2 ) 2 ) 2 2(v 2 ) 2
1 1 1 1 1
1
( p p1 (v 2 ) 2 ) t 2 2tt1 t12 ) 2 ) /( p p1 (v 2 ) 2 ) t 2 2tt1 t12 ) 2 )) t1tp1 /((t t1 )2 ) 2
2
1 1 1 1 1
log((2t 2 4tt1 2t12 2(v 2 ) 2 ( p p1 (v 2 ) 2 ) 2((t t1)2 ) 2 (( p p1 (v 2 ) 2 )2 2(v 2 ) 2
1 1 1 1
1
( p p1 (v 2 ) 2 ) t 2 2tt1 t12 ) 2 ) /( p p1 (v 2 ) 2 )) t12 p1 /((t t1) 2 ) 2
2
1 1 1 1 1
log((2t 2 4tt1 2t12 2(v 2 ) 2 ( p p1 (v2 ) 2 ) 2((t t1) 2 ) 2 (( p p1 (v2 ) 2 ) 2 2(v2 ) 2
1 1 1 1
1
( p p1 (v 2 ) 2 ) t 2 2tt1 t12 ) 2 ) /( p p1 (v 2 ) 2 )) t12 p1 /((t t1) 2 ) 2
2
1 1 1 1 1
log((2t 2 4tt1 2t12 2(v2 ) 2 ( p p1 (v2 ) 2 ) 2((t t1)2 ) 2 (( p p1 (v2 ) 2 )2 2(v2 ) 2 ) 2
1 1 1 1
1
2(v 2 ) 2 ( p p1 (v2 ) 2 ) t 2 2tt1 t12 ) 2 ) /( p p1 (v 2 ) 2 )) t1 p1
2
1 1
1
log((v2 p2 2 pp1 p12 t1t(v2 ) 2 /((t t1)2 ) 2 ...
2
Surface Kernel Solution of the Method of Moments 111
1 1 1 1 1
... log((2t 2 4tt1 2t12 2(v2 ) 2 ( p p1 (v2 ) 2 ) 2((t t1 )2 )2 (( p p1 (v2 ) 2 )2 2(v2 )2
1 1 1 1 1
1
( p p1 (v2 ) 2 ) t 2 2tt1 t12 ) 2 ) / ( p p1 (v2 ) 2 )) t1t(v2 ) 2 / ((t t1)2 ) 2
2
1 1 1 1 1
log((2t 2 4tt1 2t12 2(v2 ) 2 ( p p1 (v2 ) 2 ) 2((t t1 )2 ) 2 (( p p1 (v2 ) 2 )2 2(v2 ) 2
1 1 1 1 1
1
( p p1 (v2 ) 2 ) t 2 2tt1 t12 ) 2 ) / ( p p1 (v2 ) 2 )) t12 (v2 ) 2 / ((t t1 )2 ) 2
2
1 1 1 1 1
log((2t 2 4tt1 2t12 2(v2 ) 2 ( p p1 (v2 ) 2 ) 2((t t1 )2 )2 (( p p1 (v2 )2 )2 2(v2 )2
1 1 1 1 1
1
( p p1 (v2 ) 2 ) t 2 2tt1 t12 )2 ) / ( p p1 (v2 )2 )) t12 (v2 )2 / ((t t1 )2 )2
2
1 1 1 1 1
log((2t 2 4tt1 2t12 2(v2 )2 ( p p1 (v2 )2 ) 2((t t1 )2 )2 (( p p1 (v2 )2 )2 2(v2 )2
1 1 1
1
( p p1 (v2 )2 ) t 2 2tt1 t12 )2 ) / ( p p1 (v2 )2 )) t1v atan( (2 p 2 p1) / v)
2
In: Development of Complex Electromagnetic … ISBN: 978-1-61122-013-1
Editors: K. N. Ramli et al. © 2014 Nova Science Publishers, Inc.
Chapter 4
QUASI-STATIC FINITE-DIFFERENCE
TIME-DOMAIN SUBGRIDDING TECHNIQUE
ABSTRACT
A new approach to the modeling of electromagnetic wave
propagation and penetration in and around electrically small objects is
presented. The traveling electromagnetic wave from a source is simulated
by the finite-difference time-domain solution of Maxwell’s equations,
and a sub-gridding technique is imposed at points of interest in order to
observe the electromagnetic field at high resolution.
The computational burden caused by the requirement for a large
number of time steps has been ameliorated by implementing the state-of-
the-art quasi-static approach. The method is demonstrated by finding the
induced electromagnetic fields near a buried pipeline that runs parallel to
400-kV power transmission lines; results are presented and discussed.
114 K. N. Ramli, R. A. Abd-Alhameed and P. S. Excell
4.1. INTRODUCTION
The sharing of route by overhead high voltage power transmission lines
(OHTL) and buried utility pipelines has become quite common nowadays [1].
Most of the countries in the world typically follow the same procedure to
distribute the energy to rural and urban territories. The primary reason of this
tendency is the restriction imposed by government and private organisations
concerning the environmental effect due to the development of the facilities
needed on a given location of interest.
The numerous constraints become more alarming as new construction
reaches high density urban areas. This sharing of way prompts the question to
the public concern as to the effect of OHTL has on underground pipelines. In
some situations, the underground utility pipelines may be very close in
proximity to OHTL. It is thus necessary to take into account the
electromagnetic fields induced above the pipelines [1]. Electromagnetic
interference can be generated in the pipelines due to the electromagnetic
induction between the underground pipelines and OHTL when they are close
to each other in the vicinity. In general, the interference can be composed into
two; namely conductive interference and inductive interference. Conductive
interference is the potential increase of the ground in the vicinity of the
pipelines. It happens from a large current injected into the ground from the
transmission line particularly due to the lightning strikes between the tower
and the overhead transmission line. In contrast, inductive interference is the
voltage generated in the pipelines due to the induction of the electromagnetic
field of the OHTL.
This work presents the development of a new approach of modelling the
source excitation and the penetration of structures by continuous propagating
electromagnetic plane waves. The technique incorporates the solution of time-
dependent Maxwell’s equations and the initial value problem as the structures
are illuminated by the plane waves. The propagation of waves from source
excitation is simulated by solving a finite-difference Maxwell's equation in the
time domain. The ultimate objective of research in this area is to access the
appropriateness of the method in determining the amount of electromagnetic
penetrating fields between OHTL and underground utility pipeline. The three-
phase OHTL are modelled as the AC sources and the pipeline as the dielectric
material. In this case, the pipeline is defined as the fine grid and the residual
spaces as the coarse grid in the computational spaces. The fields between these
two grids are unknown in nature and have to be calculated. Interpolation
algorithm is thus required between the grids. The aim of the present work is to
Quasi-Static Finite-Difference Time-Domain … 115
develop the general code for solving the electric and magnetic fields within
arbitrary metal or dielectric structures, while maintaining a boundary of
uncertainty low reflection level in two-dimensional approach.
In this work, finite-difference time-domain (FDTD) method has been used
due to its easiness and potentiality to treat complex geometry structures in the
huge calculation region [2-12]. This method of solving Maxwell’s differential
equations was first proposed in two-dimensional problems [13] and then
utilized in three-dimensional applications [14]. However, the standard FDTD
method is incompetent if the details of the geometry need to be modelled due
to a global fine mesh. As a result, the total number of cells increases
dramatically. The time step must be reduced to fulfill the Courant stability
condition causing the computational time to increase significantly. The
discretisation of time step is crucial for accurate determination of the scheme
and has to be small enough to resolve different dielectric or metal structures.
A frequency domain integral equation such as method of moments is well
suited for modelling complex of antennas in free space [15]. In other words, its
strength is solving PEC structures effectively. Generally, it employs a method
of weighted residuals.
All such techniques begin by establishing trial functions with one or more
adjustable parameters, and the residuals are obtained from the differences
between trial and true solutions. The parameters are found using minimisation
to give the best fit of the trial function. In contrast, the time domain FDTD
technique is best suited for modelling electromagnetic fields inside and outside
inhomogeneous media, in particular the ground.
The presence of arbitrary inhomogeneous objects inside the computational
domain does not seriously impact the number of unknowns to be determined.
Many researchers in the past have been prompted to investigate
subgridding technique as an analysis approach in the interaction between
source and scatterer [16-21]. In general, this technique is used to condense the
lattice at the point of interest locally and does not require any analytical
formula to be taken into account, and hence it is appropriate for objects of any
shape. The residual of the space is filled with coarse grids. The fields on the
boundary between coarse and fine grids are coupled using spatial and temporal
interpolations. The regions of the coarse and fine grids are computed by the
FDTD method and are kept in time step. A stable subgridding algorithm can
refine the mesh locally and improve the accuracy of the result without
increasing the computational efforts significantly. It is hence very useful for
FDTD code.
116 K. N. Ramli, R. A. Abd-Alhameed and P. S. Excell
The method mentioned above has two main advantages relative to the
other modelling approaches. First, it is simple to implement for complicated
dielectric or metal structures due to arbitrary electrical parameters can be
assigned to each cell on the grid.
Second, the entire computational spaces need not to be discretised with a
fine grid as it put unreasonable burden on the computer processing time. The
ultimate objective of research in this area is to access the appropriateness of
the method in determining the amount of electromagnetic penetrating fields
between transmission lines and underground utility pipeline. The three-phase
transmission lines are modelled as the AC sources and the pipeline as the
dielectric material. In this case, the pipeline is defined as the fine grid and the
residual spaces as the coarse grid in the computational spaces. The fields
between these two grids are unknown in nature and have to be calculated.
Interpolation algorithm is thus required between the grids. The aim of the
present work is to develop the general code for solving the electric and
magnetic fields within arbitrary metal or dielectric structures, while
maintaining a boundary of uncertainty low reflection level in two-dimensional
approach.
48] applied the same knowledge. In this case, the dimension of the object of
study is at a fraction of the wavelength. An attempt to gain the scheme of
using the same technique but at much higher frequencies of 900 MHz and
1800 MHz has been made by See et al. in 2003 [49-50].
The theoretical method discussed by the authors in [42-44] has been
realised in the present work to approximate the quasi-static FDTD
subgridding. In general, two conditions must be satisfied before applying the
quasi-static formulation:
where:
The ‘’ symbol in equation (4.1) basically refers to vector dot product.
From this equation with the two stated conditions satisfied, the scaling
relationship can be deduced [43-45]:
f j f
E tissue f
E tissue f
f j f (4.2)
where:
Quasi-Static Finite-Difference Time-Domain … 119
E tissue f is the resultant internal electric field (V/m)
f is the scaling internal electric field (V/m)
E tissue
f is the frequency of interest (Hz)
f is the scaling frequency (Hz)
ω is the angular frequency of interest (s-1)
ω is the scaling angular frequency (s-1)
is the conductivity of the object (S/m)
is the scaling conductivity of the object (S/m)
f f
E tissue f E f
f f tissue
(4.3)
(a) Without subgrid, and (b) With subgrid, as depicted in Figure 4.2 (refers to
the x symbol outside the grid region).
The problem space was excited by sinusoidal wave and Gaussian pulse at
1800 MHz. The electric fields at the same point for case 1 and case 2 were
observed and compared as illustrated in Figures 4.3(a).
Figure 4.1. Subgridding model with field components at the main and fine grid.
a b
Figure 4.2. The observed field was located inside/outside subgrid area: (a) Without
subgrid, (b) With subgrid.
Quasi-Static Finite-Difference Time-Domain … 121
Figure 4.3. The observed fields: (a) The electric field inside subgrid and at normal grid,
(b) The magnetic field inside subgrid and at normal grid.
Sinusoidal Gaussian
Inside Subgrid Outside Subgrid Inside Subgrid Outside Subgrid
Number of
Ez at Ez at Ez at Ez at
Time Steps Ez at Ez at Ez at Ez at
Normal Normal Normal Normal
Subgrid Subgrid Subgrid Subgrid
Grid Grid Grid Grid
1 -0.1400 -0.1400 -0.1900 -0.1900 -0.1900 0.0087 -0.1900 0.0087
1000 -0.1523 -0.1527 -0.2088 -0.2095 -0.2088 0.0191 -0.2088 0.0135
2000 0.1067 0.1154 0.1881 0.2014 0.1881 0.0061 0.1881 0.0043
3000 0.1994 0.2113 0.3301 0.3483 0.3301 0.0035 0.3301 0.0025
4000 -0.0135 -0.0091 0.0038 0.0106 0.0038 0.0025 0.0038 0.0018
5000 -0.2440 -0.2477 -0.3494 -0.3550 -0.3494 0.0019 -0.3494 0.0013
6000 -0.1781 -0.1794 -0.2483 -0.2504 -0.2483 0.0014 -0.2483 0.0010
7000 0.0899 0.0980 0.1624 0.1747 0.1624 0.0012 0.1624 0.0008
8000 0.1869 0.1983 0.3110 0.3285 0.3110 0.0010 0.3110 0.0007
9000 -0.0233 -0.0192 -0.0111 -0.0049 -0.0111 0.0008 -0.0111 0.0006
10000 -0.2520 -0.2559 -0.3616 -0.3676 -0.3616 0.0007 -0.3616 0.0005
a
Sinusoidal Gaussian
Inside Subgrid Outside Subgrid Inside Subgrid Outside Subgrid
Number of
Hy at Hy at Hy at Hy at
Time Steps Hy at Hy at Hy at Hy at
Normal Normal Normal Normal
Subgrid Subgrid Subgrid Subgrid
Grid Grid Grid Grid
1 0.1000 0.1000 0.0067 0.0067 0.0039 0.0039 0.0055 0.0055
1000 0.1428 0.1441 0.0571 0.0588 0.0571 0.0588 0.0571 0.0588
2000 -0.7502 -0.7757 -1.1335 -1.1675 -1.1335 -1.1675 -1.1335 -1.1675
3000 -1.0701 -1.1052 -1.5602 -1.6070 -1.5602 -1.6070 -1.5602 -1.6070
4000 -0.3365 -0.3496 -0.5820 -0.5995 -0.5820 -0.5995 -0.5820 -0.5995
5000 0.4586 0.4693 0.4781 0.4925 0.4781 0.4925 0.4781 0.4925
6000 0.2316 0.2356 0.1755 0.1808 0.1755 0.1808 0.1755 0.1808
7000 -0.6922 -0.7160 -1.0563 -1.0879 -1.0563 -1.0879 -1.0563 -1.0879
8000 -1.0271 -1.0609 -1.5028 -1.5479 -1.5028 -1.5479 -1.5028 -1.5479
9000 -0.3028 -0.3149 -0.5371 -0.5532 -0.5371 -0.5532 -0.5371 -0.5532
10000 0.4859 0.4975 0.5145 0.5300 0.5145 0.5300 0.5145 0.5300
b
124 K. N. Ramli, R. A. Abd-Alhameed and P. S. Excell
The electric fields in subgrid region (Ezg) and at normal grid (Ez) for case 1
were found to be identical to each other to confirm the proof of concept. The
electric fields Ez with and without subgrid for case 2 were also found to be
identical to each other. A similar explanation also applies for the magnetic
fields for both case 1 and 2.
The results in Figure 4.4 illustrate the stability of the simulation inside the
problem space. The electric field remained at 0.23 V/m when using different
values of subgrid cells to justify that the results converge with the mesh size in
the computational domain.
A source code was written to implement the design and analysis of the
interaction between overhead high voltage power transmission lines and
buried utility pipeline. Fortran 90 was used as a programming language
platform. The work was devoted to 2-D TM case.
Figure 4.4. Electric field distribution for different numbers of subgrid cells in one main
FDTD cell.
Quasi-Static Finite-Difference Time-Domain … 125
Figure 4.5. Outline of standard circuit 400 kV steel lattice transmission high voltage
suspension towers with normal span of 300 m (low height construction design, not to
scale) [51, 52].
Figure 4.5 illustrates the cross section and the dimension of a common
corridor in which a buried utility pipeline runs parallel to a 400 kV overhead
power transmission line. It is designed with low height construction. The
height from the ground to the bottom conductors is 17 m. The distance from
the overhead ground wire to the earth surface is 32 m. Phase A conductors at
the top were collocated horizontally with a separation of 3.0 m. The bottom
conductors of phase B, phase B to C and phase C were collocated horizontally
with a separation of 1.75 m, 3.5 m and 1.75 m between two adjacent
conductors respectively. The three phase steel lattice transmission high voltage
suspension tower was designed with six cables. These cables were used as the
source signal which propagates inside the problem space. Each of the two
cables carries the same phase of the AC current.
The general equations of phase A, phase B and phase C cables were given
respectively by the expressions:
Phase A sin(2ft )
(4.4)
126 K. N. Ramli, R. A. Abd-Alhameed and P. S. Excell
2
Phase B sin 2ft
3
(4.5)
4
Phase C sin 2ft
3
(4.6)
where f is the frequency (Hz) and t is the time (s). The pipeline was separated
at a distance of 100 m from the steel lattice suspension towers and buried 2 m
beneath the surface of the earth. It was made from metal with a very high
conductivity of 4.75 106 S/m. The radius of the pipeline was 25 cm. The soil
in the common corridor was designed to be inhomogeneous. It was modelled
with different relative permittivity by means of random number generator. It
was known that the relative permittivity of soil varies from 1 to 5 at 460 kHz,
where as the conductivity is kept at 2.0 10-3 S/m [53].
Figure 4.6 represents a histogram that indicates the number of occurrance
that soil with their respective value of random relative permittivity from 1 to 5.
Figure 4.7 depicts the cumulative distribution function of soil relative
permittivity.
The plot is basically based on equation (7). The representation of Figures
4.6 and 4.7 clearly indicates that the soil was designed as arbitrarily
inhomogeneous.
x1
1
CDF x dx
1 4
(4.7)
1
R g 0 1 x1 1
4 (4.8)
x1 4 Rg 0 1 1
(4.9)
was modelled with r = 3.0 while its conductivity remained the same value as
before. The computational region at the coarse grids was discretised at a
spatial resolution of 2,609 cells per wavelength (y = z = 25 cm).
Subgridding involves local mesh refinement in the pipeline and some part of
the ground in order to determine the propagation of the waves inside that area
while observing the change in the electric and magnetic fields. The
computational space for main region was 521 cells × 185 cells (130.25 m ×
46.25 m) as shown in Figure 4.8. The subgrid computational space was 40
subgrid cells × 40 subgrid cells as illustrated in Figure 4.9.
In other words, there were many colours that represent random
distribution of the soil. In contrast, the non-colour (white area) depicts the
pipeline. The values for the other parameters were summarised in Table 4.2.
The distribution of ground surrounding the pipeline was generated by
using random number to simulate the inhomogeneity of the media. The fine
grids was discretised at a spatial resolution of 10,435 cells per wavelength (y
= z = 6.25 cm). In other words, the ratio of the coarse to the fine grids was
4:1. The length of the coarse grids remained at 3.83 × 10-4 of the wavelength.
The length of the fine grids remained at 9.58 × 10-5 of the wavelength. The
induced EM fields section above the pipeline were observed for 30 cells × 20
cells (7.5 m × 5 m). The Courant stability condition for 2-D case is given by:
h
t
c 2 (4.10)
Figure 4.8. The main region in the computational domain for 400 kV steel lattice
transmission high voltage suspension towers model.
Quasi-Static Finite-Difference Time-Domain … 129
Figure 4.9. The subgrid region in the computational domain for 400 kV steel lattice
transmission high voltage suspension towers model.
Parameter Measurement
Source frequency 460 kHz
Coarse grids 25 cm
Fine grids 6.25 cm
Refinement factor 4
Time step 0.4 ns
Number of time steps 21,160
Number of cycles 4
Subgrid spatial resolution 40 subgrid cells 40 subgrid cells
Induced EM fields spatial resolution 30 cells 20 cells
130 K. N. Ramli, R. A. Abd-Alhameed and P. S. Excell
According to this equation, the time step was set at 0.4 ns to satisfy the
Courant stability condition. The simulation was run for 21,160 time steps to
allow for the wave to fully traverse the spatial domain for 4 cycles.
Figure 4.10. Three-phase sinusoidal sources driven from 400 kV steel lattice
transmission high voltage suspension towers.
Figure 4.11. The amplitude of EM fields plotted against time inside subgrid region: (a)
Electric field Ezg, (b) Magnetic field Hyg and Hxg.
The distribution of the electric field Ez, magnetic fields Hy and Hx through
the simulated FDTD computational space were given in Figures 4.12, 4.13 and
4.14 respectively.
In subgrid region, the distribution of Ezg, Hyg and Hxg were demonstrated
in Figures 4.15, 4.16 and 4.17 respectively. Here, the fields inside the metallic
pipeline were also found to be zero. The reason for this phenomenon was due
to the excess electrons at the surface of the metal preventing any incoming
132 K. N. Ramli, R. A. Abd-Alhameed and P. S. Excell
propagating waves from penetrating the pipeline. It was shown that the electric
field distribution surrounding the pipeline alters from 7.9 × 10-7 V/m (-122
dBV/m) to 7.9 × 10-4 V/m (-62 dBV/m) which shows good conformity with
[54].
The difference in relative distance of each phase from the nearby pipeline
can create phase imbalance in the transmission line. Under fault condition, the
currents on the faulty phases of transmission lines were high causing induced
AC voltage on the pipeline. The induced field will not contribute to shock
hazard in normal condition.
Figure 4.12. The electric field Ez distribution in the main FDTD grid.
Figure 4.13. The magnetic field Hy distribution in the main FDTD grid.
Quasi-Static Finite-Difference Time-Domain … 133
Figure 4.14. The magnetic field Hx distribution in the main FDTD grid.
Figures 4.18, 4.19 and 4.20 illustrate the induced EM fields for Ez, Hy and
Hx respectively. As can be seen, the current amplitude induced on the pipe
varies (Hy magnetic fields in figure 15(b)) from 1.8 × 10-8 A/m (-155 dBA/m)
to 1.3 × 10-7 A/m (-138 dBA/m), and (Hx magnetic fields in figure 15(c)) from
5.6 × 10-8 A/m (-145 dBA/m) to 4.0 × 10-7 A/m (-128 dBA/m). It is
noteworthy that these currents do not produce any sudden risk to a person
nearby.
Figure 4.15. The electric field Ezg distribution in the subgrid section.
134 K. N. Ramli, R. A. Abd-Alhameed and P. S. Excell
Figure 4.16. The magnetic field Hyg distribution in the subgrid section.
Figure 4.17. The magnetic field Hxg distribution in the subgrid section.
Figure 4.18. The induced electric field Ez at 1.75 m above metallic pipeline.
Quasi-Static Finite-Difference Time-Domain … 135
Figure 4.19. The induced magnetic field Hy at 1.75 m above metallic pipeline.
Figure 4.20. The induced magnetic field Hx at 1.75 m above metallic pipeline.
CONCLUSION
An approach to model the interaction between overhead transmission lines
and underground utility pipeline at power-line frequency has been presented.
This uses the FDTD technique for the whole structure of the problem
combined with subgridding method at the object of interest particularly at the
underground pipeline. By implementing a modified version of Berenger’s
136 K. N. Ramli, R. A. Abd-Alhameed and P. S. Excell
PML, the reflection on the boundary layers inside the spatial FDTD
computational region has been successfully decreased, although it is
surrounded by lossy penetrable media. The computational burden due to huge
number of time steps in the order of tens of millions has been eased to tens of
thousands by employing the method called quasi-static approximation scheme.
The use of inhomogeneous soil in the common corridor permits a non-trivial
proximity region of authentic ground properties to be simulated. Profound
investigation of the interaction between electromagnetic fields and natural or
utility arrangement with different electrical characteristics at different level of
spatial resolution can be assisted by such tools. The combination of frequency
scaling SGFDTD approach with arbitrary inhomogeneous dielectric volume,
and the modified Berenger’s PML paves a way as a good candidate model of
EM fields interaction modelling for complex geometries.
REFERENCES
[1] A. W. Peabody and A. L. Verhiel, “The effects of high-voltage AC
transmission lines on buried pipelines,” IEEE Transactions on Industry
and General Applications, vol. IGA-7, 1971.
[2] R. Xiong, B. Chen, Y. Yi, H.-L. Chen, and H.-F. Zhao, “Two-step
method for the FDTD analysis of long apertures having depth,” IEEE
Antennas and Wireless Propagation Letters, vol. 12, pp. 39-42, 2013.
[3] S. K. Kamepally, B. P. Kumar and C. S. Paidimarry, “FDTD estimation
for accurate specific absorption rate in a tumor,” Annual International
Conference on Emerging Research Areas and International Conference
on Microelectronics, Communications and Renewable Energy
(AICERA/ICMiCR), pp. 1-5, 2013.
[4] W. G. Whittow, C. C. Njoku and Y. C. Vardaxoglou, “Fine scale
simulations of patch antennas with heterogeneous substrates,” USNC-
URSI Radio Science Meeting (Joint with AP-S Symposium), pp. 223,
2013.
[5] A. C. Lesina, A. Vaccari and A. Bozzoli, “A modified RC-FDTD
algorithm for plasmonics in Drude dispersive media nanostructures,” 6th
European Conference on Antennas and Propagation (EUCAP), pp. 687-
690, 2012.
[6] G. M. Noetscher, Y. Xu and S. N. Makarov, “Accuracy of point source
models with coincident phase centers in a cubic FDTD grid for arbitrary
Quasi-Static Finite-Difference Time-Domain … 137
[51] A. J. Pansini, Power transmission and distribution, 2nd ed. Georgia: The
Fairmont Press Inc., 2005.
[52] G. A. Goulty, Visual amenity aspects of high voltage transmission.
Somerset, England: Research Studies Press Ltd., 1990.
[53] W. M. Middleton and M. E. V. Valkenburg, Reference data for
engineers: Radio, electronics, computer, and communications, 9th ed.
Boston, 2002.
[54] G. M. Amer, “Novel technique to calculate the effect of electromagnetic
field of H.V.T.L. on the metallic pipelines by using EMTP program,”
18th International Conference and Exhibition on Electricity Distribution
(CIRED), Turin, Italy, pp. 1-5, 2005.
In: Development of Complex Electromagnetic … ISBN: 978-1-61122-013-1
Editors: K. N. Ramli et al. © 2014 Nova Science Publishers, Inc.
Chapter 5
EFFECT OF ANISOTROPIC
MAGNETO-CHIRALITY ON MICROSTRIP
RESONATOR CHARACTERISTICS
ABSTRACT
The effect of a chiral bianisotropic substrate on the complex resonant
frequency of a rectangular microstrip resonator is studied on the basis of
the integral equation formulation. The analysis is based on numerical
resolution of the integral equation using the Galerkin procedure for the
moment method in the spectral domain. This work aims to study the
effect of the chirality of a bianisotropic substrate upon the resonant
frequency and the half-power bandwidth. The effect of a chiral-substrate
bianisotropy on the surface waves of the microstrip resonator is also
studied. The effective technique used to formulate the characteristic
equations of the surface waves in a medium equipped with a complex
anisotropy is presented. The equations concerning an evaluation of the
cut-off frequencies are given in more detailed forms. A simple
approximate formula for estimating the wave number of the surface mode
TM0 and TE1 are obtained. An estimated maximum value of chiral slab
144 C. Zebiri, F. Benabdelaziz and M. Lashab
5.1. INTRODUCTION
In recent years, the method of moments (MoM) [1-6] has been adopted in
the analysis of radiating elements such as resonators [7] and antennas [8, 9].
Microstrip resonators appeared during the fifties, and they were especially
developed during the seventies [10, 11]. They are small in size and have
simplicity, easy manufacturing and easy handling. They are robust when they
are fixed on rigid surfaces, and are also very effective in terms of resonance,
polarization, input impedance and radiation patterns [12-14]. The major
disadvantages of microstrip resonators are their weak purity of polarization,
and narrow half-power bandwidth which is typically a few percent. However,
increasing substrate thickness and reduced relative permittivity make it
possible to improve resonator efficiency by up to 90%, and to increase
bandwidth by up to 35%, provided that the surface waves are minimized [15].
The resonant frequency of microstrip resonator is strongly related to the
dielectric of the substrate, and remains a very important factor to determine in
resonator design. Some works have mentioned that many uniaxial anisotropic
substrates have been used for printed circuit resonators [16, 17]. However,
studies concerning the resonant frequency calculation of the microstrip
resonator, implemented on a substrate of uniaxial anisotropy [18], have shown
that this kind of substrate does not affect the quality factor, the half-power
bandwidth and obviously the resonant frequency. For these reasons, other
studies have focussed on cases of positive and negative anisotropy [18-20], or
bianisotropic media [21]. In the present work, the Spectral Domain Approach
(SDA) is applied, and intensively used in the analysis and design of planar
Effect of Anisotropic Magneto-Chirality ... 145
1
B H E (5.1)
c0
1
D E H (5.2)
c0
where the permittivity and permeability tensors are of uniaxial anisotropy and
the magneto-electric elements are respectively expressed as
146 C. Zebiri, F. Benabdelaziz and M. Lashab
t 0 0 t 0 0
0 t 0 0 t 0
0 0 z 0 0 z
, ,
0 0
j 0 0
0 0 0
(5.3)
~
E~e
S s , z
ES s , z ~h
1 E~
F (s , rs ) ~
x
ES s , z s Ey
z 1
j 0 0 ~
Ez
t s
z
1 ~
0 0 z Hz
s
(5.4)
~
H~e
S s , z
HS s , z ~ h
1 H~
F (s , rs ) ~
y
HS s , z s Hx
1
0 z 0 ~
Ez
s ~
0 z 1 Hz
j 0
t s z
(5.5)
Effect of Anisotropic Magneto-Chirality ... 147
y
F s , rs x , s x xˆ y yˆ , rs xxˆ yyˆ (5.6)
y x
~ ~
h in (5.4) and (5.5) denote the TM and TE waves respectively, E S and H S are
expressed according to the TE and TM modes.
~
e e
E z s , z A e e j z z B e e j z z (5.7)
~
h h
H z s , z A h e j z z B h e j z z (5.8)
where
ze 2 02 t t 2 s2 t
(5.9)
z
zh 2 02 t t 2 s2 t
(5.10)
z
s , ze and zh . These are respectively the free space propagation for TE and
TM modes. According to the previous equations, the following tangential
components of the fields are obtained
~
~
E e s , z
ES s , z ~ h
j z z
j z z
e A s e B s (5.11)
E s , z
148 C. Zebiri, F. Benabdelaziz and M. Lashab
~
~
H e s , z
HS s , z ~ h
j z z j z z
e g s A s e h s B s
H s , z
(5.12)
where
z diag ze zh (5.13)
A s j 12
s
z
t
j B
0
e
z
e 1
s2
0 z Ah
T
(5.14)
B s j 12
s
z
t
j B
0
e
z
e 1
s2
0 z B h
T
(5.15)
g s diag
0 t
j 0 j zh
j 0 j z
e
0 t
(5.16)
h s diag
0 t
j 0 j zh
j 0 j z
e
0 t
(5.17)
For an electric and non magnetic medium having biaxial anisotropy with
regard to the permittivity, the previous expressions are well detailed in [18-
20].
Ne
G s
1
j 0
1
diag e z ze , h 02 t
sin z d (5.18)
D D
where
Ne
1
e
2
0 t t t s2
z
(5.19)
z
D e ze z t cos zh d j 02 t t t s2 j z 0 t sin ze d
z
(5.20)
D h zh cos zh d j z t j 0 sin zh d (5.21)
N
J r M 0
J rs an xn s bm (5.22)
n 1 0 m1 J ym rs
d F ~
s s , rs G s J s 0 (5.23)
B1 N N B
2 N M a N 1
0
B3 M N B (5.24)
4 b M 1
M M
det(B()) 0 (5.25)
The resonator is designed to operate near its resonant frequency, and all its
characteristics are estimated around this frequency. Equation (5.25) is called
the characteristic equation for the complex resonant frequency f f r if i ,
where f r is the resonant frequency and f i is the dampling factor that
characterises the losses of the radiating antenna. The quality factor and the
half-power bandwidth are defined in [31, 32] as
Q f r (2 f i )
(5.26)
BW 1 Q (5.27)
The effect of the two chiral constitutive parameters has been studied.
These are first the chirality and the magnetic uniaxial anisotropy of the
substrate upon the complex resonant frequency and secondly the half-power
bandwidth of the microstrip resonator. In order to validate the results (in both
isotropic and chiral cases), we have made a comparison with the isotropic case
Effect of Anisotropic Magneto-Chirality ... 151
presented in the literature. The figures presented below show the normalized
real and imaginary parts of the resonant frequency, and the half-power
bandwidth, compared with the measurements in the literature with x=z=2.35
and x=z=7. The medium considered is defined as follows: ==-1, 0 and 1
respectively (magneto-electric elements), t=0.8, z=1 (positive anisotropy of
the permeability), t=1.2, z=1 (negative anisotropy). Uniaxial anisotropy is
obtained by changing t and keeping z constant.
G k s
d 0
d
j
diag 02 t 1 s2 ,t 02
z
(5.28)
Some authors such as Chew [19], have shown that only one basis function
(N=0, M=1) is needed to obtain excellent convergence of the results. In this
case, the current distribution on the patch conductor is given by
152 C. Zebiri, F. Benabdelaziz and M. Lashab
0 0
J rs b1 1 sin y
b b (5.29)
J r 2
b
y1 s
and the characteristic equation for the resonant frequency driven from the
fourth element of the matrix equation (5.24) is given as
d ~
d
B 11
2 ~
t 1 y2 J y1 s J y1 s 0 (5.30)
j
4 s 0
z
The function Jy1 has a Fourier transform given by the analytical expression
~ sin x a
cos y b
J y1 s b 2 2
2 2 y b2 2
(5.31)
x
4
02 t I1 I2 0 (5.32)
z b2
where
cos2 y
y2 cos2 y
I1 d y
, I 2 dk y
(5.33)
2 2
2 2 2 2 2 2
0 y 0 y
Using the integration along the contour of the spectral domain, the integral
equations of (5.33) is resolved analytically as
1
I1 , I2
4
(5.34)
Effect of Anisotropic Magneto-Chirality ... 153
c
fr (5.35)
2b z t
where c is the speed of light in the vacuum. The resulting formula gives an
idea of the effect of the parameters z and t upon the antenna resonant
frequency. It is clear that this depends on the permittivity along the optic axis
and the perpendicular component of the permeability in the limit of a thin
substrate, hence the resonant frequency depends only on the components z
and t.
1,00 0,025
[20]
Normalized frequency
[18]
Normalized frequeny
0,98
0,020
0,96
0,015
0,94
[20] 0,010
0,92 [18]
0,005
0,90
0,88 0,000
0,00 0,05 0,10 0,15 0,20 0,00 0,05 0,10 0,15 0,20
(a) d/b (b) d/b
5,4
[20]
[18]
4,5
3,6
BW %
2,7
1,8
0,9
0,0
0,00 0,05 0,10 0,15 0,20
(c) d/b
Figures 5.2. Normalized resonant frequency, a=1.5 cm, b=1 cm, r=7, (a) Real part, (b)
Imaginary part, (c) Half-power bandwidth.
154 C. Zebiri, F. Benabdelaziz and M. Lashab
1,00 0,063
[20]
0,98
Normalized frequeny
Normalized frequeny
0,054
0,96
0,045
0,94
0,92 0,036
0,90
[20] 0,027
0,88
0,018
0,86
0,009
0,84
0,82 0,000
0,02 0,06 0,10 0,14 0,18 0,02 0,06 0,10 0,14 0,18
(a) d/b (b) d/b
0,15 [20]
0,12
BW (%)
0,09
0,06
0,03
0,00
0,02 0,06 0,10 0,14 0,18
(c) d/b
Figures 5.3. Normalized resonant frequency, a=1.5 cm, b=1 cm, r=2.35, (a) Real part,
(b) Imaginary part, (c) Half-power bandwidth.
The effect of chirality is remarkable only for thick layers, whereas for
those infinitely small, no effect is perceived, which is confirmed by
the asymptotic form expression of the resonant frequency given by
Equation (5.35).
Effect of Anisotropic Magneto-Chirality ... 155
Normalized frequeny
Normalized frequeny
0,04
t=1.2, z=1
1,02
0,03
0,96
0,02
0,90
[20]
t=z=1
0,84
0,01 t=0.8, z=1
t=1.2, z=1
0,78 0,00
0,00 0,05 0,10 0,15 0,20 0,00 0,05 0,10 0,15 0,20
0,12
0,09
BW (%)
0,06
[20]
0,03 t=z=1
t=0.8, z=1
t=1.2, z=1
0,00
0,00 0,05 0,10 0,15 0,20
(c) d/b
Figures 5.4. Normalized resonant frequency, a=1.5 cm, b=1 cm, r=2.35, (a) Real part,
(b) Imaginary part, (c) Half-power bandwidth.
It can be seen that from Figures 5.4, in the case of permeability values less
than 1, for example μt=0.8, that the effect of uniaxial anisotropy of the
permeability clearly results in a large increase in the real part of the resonance
frequency, but an unimportant increase in the imaginary part of the resonance
frequency. Therefore, the bandwidth undergoes a slight reduction.
Permeability values greater than one, such as μt=1.2, cause variations in the
opposite sense with respect to the preceding case. The chirality is apparent
156 C. Zebiri, F. Benabdelaziz and M. Lashab
only in thick layers, whereas the effect of the permeability remains whatever
the thickness of the support. According to these results one can predict the
desired values of three parameters: fr, fi and bandwidth of the resonator, as
illustrated by Figure 5.5.
Normalized frequeny
1,050
0,045 t=0.8, z=1
0,975
0,036 t=0.8, z=1
0,900 0,027
0,018
0,825
t=0.8, z=1 0,009
t=0.8, z=1
0,750 0,000
0,00 0,05 0,10 0,15 0,20 0,00 0,05 0,10 0,15 0,20
(a) d/b
(b) d/b
0,16 [20]
t=1.2, z=1
0,14
t=1.2, z=1
0,12
t=0.8, z=1
0,10 t=0.8, z=1
BW (%)
0,08
0,06
0,04
0,02
0,00
0,00 0,05 0,10 0,15 0,20
(c) d/b
Figures 5.5. Normalized resonant frequency, a=1.5 cm, b=1 cm, r=2.35, (a) Real part,
(b) Imaginary part, (c) Half-power bandwidth.
Magnetized ferrites belong to the complex material class which have proved to
have potential applications as substrates in the field of MIC and APC.
Measurements have confirmed that the resonant frequencies of the microstrip
structures printed on ferrite substrates can be fixed to suit a range of
applications, simply by adjusting the polarization of an external magnetic field
[33]. However, ferrites can be employed to reduce the microstrip resonator
[34, 35] and to design circularly polarized antennas by the application of a
simple feed [36].
The possibility of employing chiral materials as substrates for the design
of MIC and APC, was reported by Lindell [27]. On the other hand, Pozar [37]
underlined the serious disadvantages of the use of these materials as substrates,
because of the losses due to the surface wave excitation and the significant
appearance of poles during numerical calculations. However, Toscano et al.
[38] and Zebiri et al. [39, 40] have recently proved that chiral substrates can
advantageously be employed to increase the bandwidth of microstrip antennas.
This result indicates that more attention must be paid to detailed studies of the
surface wave excitation in the microstrip antennas over chiral substrates. This
study evaluates the cut-off frequencies of TE0 and TM1 modes of a resonator
printed on an anisotropic medium in order to optimize these structures
according to their applications.
These resonators appeared during the 1950s and especially were
developed during the 1970s [10, 11]. They combine small size, simplicity,
facility of manufacture and practical implementation. Moreover, they are
easily matched to plane and non-plane surfaces, and they exhibit a great
robustness when assembled on rigid surfaces. They are also very effective in
terms of resonance, polarization, input impedance and radiation diagram [12,
41].
The major disadvantages of microstrip resonators lie in their low purity of
polarization. As regards antennas we can consider the reduction of their
bandwidth which is typically a few percent [41]. However, increasing
substrate thickness and reduced relative permittivity make it possible to
improve resonator efficiency by up to 90%, and to increase bandwidth by up to
35%, provided that the surface waves are minimized [15]. The idea of using
chiral materials as substrates and superstrates in the design of printed antenna
was first presented by Engheta [42] and the term “chirostrip“ was then
invented. In the literature, it is shown that the power of the surface wave can
generally be reduced when a chiral substrate is employed for antennas
intended for printed circuitry [43].
158 C. Zebiri, F. Benabdelaziz and M. Lashab
The surface wave modes are guided by the interface separating two
different dielectric media [44]. In the design of microstrip circuits as well as in
their applications, the study of surface waves is justified by one of the two
following reasons; to minimize their effects or to use them in certain
applications [45]. It is initially essential to take into account the propagation of
modes in the chiral layer. We study the propagation of these modes along z-
axis. For such propagation, the interface of the structure can be considered as a
particular case of a chiro-waveguide [46].
In [47] it is proved that in all chiro-waveguides the proper modes are
hybrid [46], but for the case of the selected chiral, decomposition into TE and
TM modes is possible [39, 40, 48, 49]. Their study would be also made during
the evaluation of the improper integrals appearing in the formulation of the
problem of a microstrip structure by the integral method. It would be also
necessary, in certain situations, to determine as a preliminary the propagating
modes of the surface waves which appear in a microstrip structure. The
spectral method and the tensorial Green’s function are adopted as an analysing
tool. We can mention as an example the case where the integration of the
characteristic matrix elements deduced from the Galerkin’s method is
performed on the wave numbers real axis [50, 51].
Moreover, if we are interested in the fields radiated by a microstrip
structure without neglecting the radiation due to the surface waves, the
determination of their wave numbers is necessary [52]. However, the majority
of the studies relating to this subject deal with the problem of the surface
waves as being a procedure of solving two transcendent equations for which
the expressions appear in the denominators of the Green’s function [32, 52,
53].
The study in the case of simple structures, made up of one or two
dielectric layers, does not present difficulties [54, 55]. However, for the
structures that consist of more than two layers, the study becomes increasingly
complicated, requiring tedious calculations. As a result, the resolution of the
resulting characteristic equations and the asymptotic behavior of the solutions
for low operating frequencies and low dielectric-layer thicknesses is more
complicated. Traditional methods do not lead, in general, to simple
expressions without significant algebraic calculation.
Effect of Anisotropic Magneto-Chirality ... 159
nc0
fc , n=0, 1, 2, … (5.36)
4h r 1
where c0 is the speed of light in vaccuo, and h and εr are respectively the
thickness and the permittivity of the substrate. Odd values of n give the cut-off
frequencies for the TEn modes, and even values of n give cut-off frequencies
for the TMn modes. For the TE1 mode, the computed values of the ratio h c1
are 0.217 for Duroid (εr=2.32) and 0.0833 for alumina (εr=10), using the
following equation
h n
1 (5.37)
c 4 r 1
Thus, the low order mode TE1 is excited at 41 GHz for a Duroid substrate
of thickness 1.6 mm, and at approximately 39 GHz for an alumina substrate of
thickness 0.635 mm. The thickness of the substrate is chosen so that the h 0
1
ratio is lower than the h c (λ0 is the wavelength in open space at the
operating frequency), which gives [57]
c
h (5.38)
4 fu r 1
where fu is the highest frequency in the operation band. Note that h should be
chosen as high as possible, under the constraint indicated in (5.38), so that the
160 C. Zebiri, F. Benabdelaziz and M. Lashab
0 .3 c
h (5.39)
2f u r
The TM0 mode has no cut-off frequency and is always present to a certain
extent. The excitation of the surface wave TM0 mode becomes appreciable
when h/λ>0.09 (εr≅2.3) and when h/λ>0.03 (εr≅10). Generally, suppressing
the TM0 mode requires lower permittivity and smaller substrate size.
Unfortunately, decreasing ε tends to increase the antenna size, whereas
decreasing h thus leads to a limitation of the antenna effectiveness, hence a
reduction of its bandwidth. The reduction of surface waves, for microstrip
antennas printed on ferrites, has been discussed by many researchers [59-63].
In [64-66], electromagnetic band gap (EBG) structures, also known under the
name of photonic crystals, were found to reduce and in certain cases to
eliminate the surface waves, leading to increases in directivity and bandwidth.
In this context, Yang [67] was the first who proposed antennas with large gain,
which could be carried out with a unique radiant element printed on a 2-D
photonic bandgap (PBG) material.
D e ze z t cos ze d j 02 t t t 2 j z 0 t sin ze d 0
z
(5.40)
For TE modes
0 0
e2 2
2
lim z 0 t t tz (5.42)
z j 02
z 1t
d
(5.43)
z 1 0d
z j 02 d
z 1 (5.44)
z
162 C. Zebiri, F. Benabdelaziz and M. Lashab
0 0
2
lim z 0 t t
h2 2 t
z
(5.45)
1
z 1 0 d (5.46)
t d
The wavenumbers in Equations (5.43) and (5.46) are general formulas for
the anisotropic case with chiral, they can be shown to equal those obtained by
Peixeiro et al. [55]. The development of an approximate formula for the
wavenumbers of TM0 and TE1 modes gives for the TM0 mode
1 1
z j 02 t 1 d (5.47)
1 0d z t
1
z 1 0 d (5.48)
t d
The chirality effect on the two wavenumbers is similar, for the case of
very thin layers, while bringing the ratio 1/1+0d closer to 1-0d. According
to the equation z2 02 s2 . By replacing z in the Equation (5.42), an
approximate expression can be found as
2d2 1 2
1 2
z t
2 2 0
(5.49)
sp 0 z 1 0 d 2
Effect of Anisotropic Magneto-Chirality ... 163
1
02 d 2 2
TM 0 1 2
1
12
(5.50)
0
z 1 0 d 2 z t
1 d 2
2
1
2 0
(5.51)
sp
0 t d
0
1
1 d 2 2
TE 0 1 0
(5.52)
1 0 t d
Equations (5.50) and (5.52) are useful as a good initial estimate for the
pole location, in a standard routine, to seek the true solution, as shown in
Figure 5.7. The approximation of the pole location in the integration path is
harmful for effective numerical evaluations of the integrals. It is advantageous
to subtract the pole singularity and to reinstate it thereafter analytically. In the
case of a thick substrate, several poles may exist, and the analytical evaluation
of the integrals around the half-circles becomes quite complicated if two or
several poles are very close to each other [68]. The total number of the poles is
determined by the operating frequency and the dimensional parameters of the
substrate [52, 68]. In the isotropic case, it is possible to predict the number of
poles, such as in [69], if we have 0 d r 1 2 then the denominator of
Gh(s) does not have any zero, while that of Ge(s) has only one.
In the case of a structure having only one isotropic dielectric layer with
z=t=r, z=t=1 and =0 the equality (5.51) reduces to the simple expression
given by [63, 70]
r 1
2
1 0 d
2 2 2 2
(5.53)
sp r2
0
164 C. Zebiri, F. Benabdelaziz and M. Lashab
0 d
2 02 (5.54)
1 0d
It is sufficient to put z=0 in Equations (5.40) and (5.41), to deduce the
cut-off frequencies of TE and TM modes. For TM modes
nc 1
f cTM , n is even. (5.55)
4 d
t t z 2
t
For TE modes
c 1
n 1
t t 2 t
f TE
arctan ,
z
n
c
d 2 t
2 2
t t z
is odd. (5.56)
nc 1
f cTM (5.57)
4 d
t t z
t
2n 1 c 1
f cTE
(5.58)
4 d t
t t z
nc 1
f cTM (5.59)
4 d
t z
t
2n 1 c 1
f cTE (5.60)
4 d t 1
nc 1
f cTM (5.61)
4 d r 1
2n 1 c 1
f cTE (5.62)
4 d r 1
The previous equations have the same forms of those developed in [55,
57]. The mode of the low order TE1 is excited at 40.77 GHz in a duroid
(r=2.32) substrate of thickness 1.6 mm without the effect of the chirality [57].
On the other hand, using the same substrate with the consideration of a
positive (=+1) or negative (=-1) coefficient of chirality, the TE1 mode for
the same structure is excited at 138.40 GHz or 192.75 GHz respectively.
Similarly, the low order TE1 is excited at approximately 39.33 GHz for an
isotropic 0.635 mm thick alumina (r=10) substrate [57], but at 50.75 GHz or
166 C. Zebiri, F. Benabdelaziz and M. Lashab
c 1
n 1
t t
2 t
d arctan (5.63)
z
f 2 t 2 2
t t
z
300 x=posititive
250
x=négative
200
150
100
50
0
0,85 0,90 0,95 1,00 1,05 1,10 1,15
Figure 5.8. Effect of the chiral constitutive parameters on the TE0 mode cut-off
frequency versus different ratios. TE0 is excited at 40.77 GHz in a Duroid (r=2.32)
substrate of thickness 1.6 mm [57].
Effect of Anisotropic Magneto-Chirality ... 167
600
isotropic case [57]
550
x=z/t, t=2.35
450
x=t/z, z=2.35
400 x=z/t, t=1.
350 x=t/z, z=1.
300
x=posititive
250
x=négative
200
150
100
50
0
0,85 0,90 0,95 1,00 1,05 1,10 1,15
x
Figure 5.9. Effect of the chiral constitutive parameters on the TM1 mode cut-off
frequency versus different ratios, in a Duroid (r=2.32) substrate of thickness 1.6 mm
[57]. isotropic again.
M N
Z in I xmVxm I ynV yn (5.64)
m1 n1
where Ixm and Iyn are the coefficients to be determined, and m and n define the
number of basis functions along the directions x and y, respectively.
d
1
i x x f y y f
Vxm
4 2 J
xm e G
0
zx dzd x d y (5.65)
d
1
i x x f y y f
V yn
4 2
J yne G
0
zy dzd x d y (5.66)
z sin ze d y J y
~
t 1
V 2 d x d y
p
(5.67)
4 0 z zeT
m
Effect of Anisotropic Magneto-Chirality ... 169
z sin ze d x J x
~
t 1
V 2 d x d y
p
(5.68)
4 0 z zeT
n
T t z cos ze d j z ze sin ze d 1 0 e tan ze d (5.69)
z
where the surface currents J on the patch are expanded into a finite series of
basis function Jxn and Jym. For digital convergence of Equations (5.67)-(5.69),
the steps discussed in [57, 80, 81] are followed.
Figures 5.9 and 5.10 show that the input impedance is directly related to
the chirality, which leads us to consider this effect in seeking a better
adaptation. The amplitude parameters raised vary depending on the chirality,
while in [72], where the case of chiral bi-isotropic is considered, these
parameters are constant in amplitude, which is advantageous for the circuit
design.
1,2 =-1
isotropic case
normalized input impedance (real part)
=1
1,0
0,8
0,6
0,4
0,2
0,0
1,16 1,17 1,18 1,19 1,20 1,21
f (GHz)
Figure 5.9. Chirality effect on the real part of the input impedance (a=7.62 cm,
b=11.43 cm, x0=1.52 cm, y0=0.385 cm, εx=εz=2.64).
170 C. Zebiri, F. Benabdelaziz and M. Lashab
1,25 =-1
0,50
0,25
0,00
-0,25
-0,50
-0,75
-1,00
1,16 1,17 1,18 1,19 1,20 1,21
f (GHz)
Figure 5.10. Chirality effect on the imaginary part of the input impedance (a=7.62 cm,
b=11.43 cm, x0=1.52 cm, y0=0.385 cm, εx=εz=2.64).
shape of the magneto-electric elements of the chiral. This case of study of the
TM or TE surface waves cut-off frequencies, with respect to the constitutive
parameters of the medium, ended in results which show clearly that the
possibility of excitation of the surface waves is more favoured when the
dielectric become increasingly thick, compared to components with traditional
dielectrics. Thus chirality presents interesting effects on the resonant
frequency, bandwidth and input impedance.
on the resonant frequency and the bandwidth of the rectangular microstrip are
investigated. The idea of using chiral materials as a substrate and superstrate in
printed antenna circuits was first introduced by Engheta [42]. It has been
shown that the power of the surface waves can be generally reduced when a
chiral substrate is employed in a printed microstrip circuit [90]. The theoretical
concept of the complex resonant frequency considered here is formulated in
terms of an integral equation, by using the vector Fourier transform [91, 92].
This part of work is presented as follows. Firstly, the bianisotropic material is
presented with its magneto-electric tensors, with a chiral material as substrate
and superstrate, Green’s tensors are deduced and the integral equation is
obtained in the spectral domain. Secondly, the numerical results are presented
for different values of the magneto-electric elements. In this contribution, the
intention is to show and observe the effects of chiral elements on the resonant
frequency and half-bandwidth of a microstrip antenna.
1
B μH ξE (5.70)
c0
1
D εE ηH (5.71)
c0
where E and H are the electric and magnetic field vectors, D and B are
the electric and magnetic induction vectors, c0 is the speed of light in vacuum
and , , and are the constitutive tensors, expressed as
Effect of Anisotropic Magneto-Chirality ... 173
μt 0 0 εt 0 0 0 ξ 0
μ μ0 0 μt 0 , ε ε0 0 εt
0, ξ η j ξ 0 0
0 0 μ z 0 0 ε z 0 0 0
(5.72)
N ie N ce e e
N κ z ,1
0
G κs
1 Die Dce
sin κ z ,1d1 (5.73)
jωε
0
N i
h
N ch 2
κ0
D i
h
Dch
εt ,1
κ 02 εt ,1 μt ,1 κ s2
ε z ,1
Ne (5.74)
κ ze,12
κ ze,2
Nie κ z cos κ z,e 2 d2 j εt,2
sin κ z,e 2 d2 (5.75)
Nce
κ 0ξ 2
κ ze,2
κ z εt,2 jκ0ξ2 sin κ z,e 2 d2 (5.76)
κz
Nih cos κ z,h2 d2 j h
κ z,2
μt,2 sin κ z,h2 d2 (5.77)
Nch j
κ0 ξ 2
h
κ z,2
sin κ z,h2 d2 (5.78)
Die κ z εt,1 cos κ z,e 1d1 jκ ze,1 sin κ z,e 1d1 cos κ z,e 2 d2
(5.79)
ε
sin κ
e
j t,1 κ ze,2 cos κ z,e 1d1 j
ε
κ κ
κ
ε t,2 sin κ z,e 1d1 z,1 z
e
e
z,2 d2
t,2 z ,2
κ0ξ1
Dce
κez,1
κz εt,1 jκ0ξ1 sin κz,e1d1 cos κz,e 2 d2
εt,1 κ0ξ2
εt,2 κz,e2
κz εt,2 jκ0ξ2 sin κz,e 2 d2 cos κz,e1d1
ε ξ κ ξ ξ κe ε κe
κ0 t,1 1 0e 1 e 2 κz εt,2 jκ0ξ2 jξ2 ez,1 jξ1 t,1 z,e2 sin κz,e 2d2 sin κz,e1d1
εt,2 ξ2 κz,1κz,2 κz,2 εt,2 κz,1
(5.80)
Effect of Anisotropic Magneto-Chirality ... 175
1
κz,h1 cos κz,h1d1 jκz μt,1 sin κz,h1d1 cos κ h
d
z,2 2
(5.81)
Dih κ κh μ
μt,1
j μt,2 z h z,1 cos κz,h1d1
κz,2
j t,1 κz,h2 sin κz,h1d1 sin κz,h2 d2
μt,2
ξ 2 μt,1 ξ1 μt,2 jκ0 ξ 2 κ z μt,2
κ0
j
κ z,h2 μt,2
sin κ z,h1d1 sin κ z,h2 d2
(5.82)
Dc
h
μt,1 κh
hz,1 ξ2 sin κ z,h2 d2 cos κ z,h1d1 ξ1 cos κ z,h2 d2 sin κ z,h1d1
κ z,2
κ z κ 02 κ s2 (5.83)
It is clear from (5.73) that the Green’s dyadic is diagonal in the (TM, TE)
representation in the presence of uniaxial anisotropic chirality can be
simplified when deriving Green’s dyadic for a multilayer anisotropic media.
Furthermore, if the substrate and superstrate layers are isotropic the
coefficients with ‘c’ indices are cancelled. The dyadic elements given in (5.73)
are the same as given in [20] and [97]. The surface current J on the patch can
be expanded into a finite series of basis functions
N
J r M 0
J rs an xn s bm (5.84)
n 1 0 m1 J ym rs
where an and bm are the mode expansion coefficients to be sought. The integral
equation describing the electric field on the patch is [98]
dκ F κ ,r G κ J κ 0
~
s s s s s (5.85)
176 C. Zebiri, F. Benabdelaziz and M. Lashab
F κ s ,rs rs 0
is the kernel of the vector transforms, where
κ κy
F κ s ,rs
rs 0
x
~
, and J κ s is the Fourier transformation of the
κ x
κ y
surface current J. Using the well known Galerkin procedure of the method of
moments, the integral equation in (5.85) is discretized and brought into the
following matrix equation
( Β1 ) N N ( B2 ) N M (a) N 1
( Β ) 0
( Β4 ) M M (b) M 1
(5.86)
3 M N
where (B1)NN, (B2)NM, (B3)MN and (B4)MM are the elements of the digitalized
matrix equation. (a)N1 and (a)M1 are respectively vectors of the unknown
~ ~
expansion coefficients of patch surface currents J xn s and J ym s . A
nontrivial solution of (5.86) is derived by seeking a complex frequency f=fr+ifi,
where fr and fi are respectively the real and imaginary parts of the complex
resonant frequency and 2fi/f is the half-power bandwidth of the antenna [18,
99, 100]. Owing to the fact that the narrow bandwidth of the resonator should
be exact, the obtained results are validated with those in [20] and [89]. Table
5.1 presents a comparison between the calculated frequencies and those
measured by Bahl [89] and calculated in [20]. The superstrate and the
substrate are considered in the isotropic case. The values are closer to those
given in [89] than to those given in [20]. Table 5.1 shows a better convergence
of our results.
In Figure 5.12, the effect of the electric uniaxial anisotropy is presented by
comparison with the isotropic case [20] (a rectangular patch trapped between a
superstrate and a substrate; a=6 cm, b=5 cm, d1=0.1 cm, t,1=z,1=2.35), for
various values of the permittivity (t,2=z,2=1.5, 2.35, 4, 10) and superstrate
thickness.
Effect of Anisotropic Magneto-Chirality ... 177
[20] [20]
Normalized real frequeny 1,002 t,2/z,2=0.8, z,2 =1.5 1,002 t,2/z,2=0.8, z,2 =2.35
1,000
t,2/z,2=0.8, z,2 =4. 1,00
0,975 0,94
0,970 0,93
0,965 0,92
0,0 0,5 1,0 1,5 2,0 2,5 3,0 0,0 0,5 1,0 1,5 2,0 2,5 3,0
[20]
0,990
0,996
0,992
0,975
0,0 0,5 1,0 1,5 2,0 2,5 3,0 0,0 0,5 1,0 1,5 2,0 2,5 3,0
[20] [20]
1,014 t,2/z,2=0.8, z,2 =1.5
1,020
t,2/z,2=0.8, z,2=2.35
t,2/z,2=1.2, z,2 =1.5 t,2/z,2=1.2, z,2=2.35
Normalized BW
Normalized BW
1,008
z,2/t,2=0.8, t,2 =1.5 z,2/t,2=0.8, t,2=2.35
1,005
z,2/t,2=1.2, t,2 =1.5 z,2/t,2=1.2, t,2=2.35
1,002
0,996 0,990
0,0 0,5 1,0 1,5 2,0 2,5 3,0 0,0 0,5 1,0 1,5 2,0 2,5 3,0
(a) d2/d1 (b) d2/d1
[20] [20]
1,02
t,2/z,2=0.8, z,2 =4. 1,04
t,2/z,2=0.8, z,2 =10
t,2/z,2=1.2, z,2 =4.
Normalized BW
1,02
t,2/z,2=1.2, z,2 =10
Normalized BW
1,01
0,96
0,98
0,94
0,97
0,0 0,5 1,0 1,5 2,0 2,5 3,0 0,0 0,5 1,0 1,5 2,0 2,5 3,0
In this region, the components t,2 hold the bandwidth close to BW0 or
higher. For dielectric permittivity (z,2, t,2>t,1=z,1) and higher thicknesses
(3d1>d2>2.6d1), the components z,2 and t,2 have a tiny effect on the
bandwidth, but for much higher values (z,2 or t,2=10), the component t,2 has
an important effect. In Figure 5.15, the effect of the uniaxial permeability is
presented in comparison with the isotropic case, for various values of the
permeability and thickness of the superstrate. According to Figure 5.15, while
increasing the value of the permittivity (t,2,z,2>r,1), the effect of the two axial
components of the permeability (µt,2, µz,2) becomes more negligible, whereas
for weak permittivity (1<t,2,z,2r,1), the component µz,2 has a stronger effect
on the real part of the resonant frequency than the component µt,2.
Figure 5.16 shows that for low permittivity (1<t,2,z,2r,1), the
component of the permeability µz,2 for (d2<1.5d1) is the most influential factor
on the imaginary resonant frequency, whereas for thick layers (d2>2.6d1), the
component µt,2 becomes more important than µz,2. For higher permittivity
(t,2,z,2>r,1), one can distinguish two regions. For the first one (d2<1.5d1), the
anisotropy in this case has a negligible effect (0.5% maximum), while in the
Effect of Anisotropic Magneto-Chirality ... 181
second region (d2>1.5d1), the component µt,2 has a visible effect (up to 2.5%).
Figure 5.17 shows that for (d2<d1), the influencing component on the
bandwidth is µz,2. Unfortunately, this effect remains weak in this region,
whereas for thick layers (d2>d1), the component µt,2 becomes more important
than µz,2, and the bandwidth in this case reaches a maximum of 3%. For higher
thicknesses of the superstrate (d2>d1), the effect of permeability on the
imaginary frequency and bandwidth is more important than the permittivity,
whereas for real frequency the latter constitutive elements have a reverse
effect.
0,994
0,985
0,993
0,0 0,5 1,0 1,5 2,0 2,5 3,0 0,0 0,5 1,0 1,5 2,0 2,5 3,0
1,000 1,00
Normalized real frequeny
0,0 0,5 1,0 1,5 2,0 2,5 3,0 0,0 0,5 1,0 1,5 2,0 2,5 3,0
(c) d2/d1 (d) d2/d1
0,992 0,980
0,0 0,5 1,0 1,5 2,0 2,5 3,0 0,0 0,5 1,0 1,5 2,0 2,5 3,0
(c) d 2 /d 1 (d) d 2 /d 1
Normalized BW
1,015
1,010
z,2/ t,2=0.8, t,2 =1. 1,010 z,2/t,2=0.8, t,2 =1.
1,005 z,2/ t,2=1.2, t,2 =1. 1,005
z,2/t,2=1.2, t,2 =1.
1,000
1,000 0,995
0,990
0,995
0,985
0,0 0,5 1,0 1,5 2,0 2,5 3,0 0,0 0,5 1,0 1,5 2,0 2,5 3,0
1,02
1,03 t,2/z,2=1.2, z,2 =1.
1,01 z,2/t,2=0.8, t,2 =1. 1,02
z,2/t,2=0.8, t,2 =1.
1,01
1,00 z,2/t,2=1.2, t,2 =1. 1,00
z,2/t,2=1.2, t,2 =1.
0,99
0,99 0,98
0,97
0,98 0,96
0,95
0,97 0,94
0,0 0,5 1,0 1,5 2,0 2,5 3,0 0,0 0,5 1,0 1,5 2,0 2,5 3,0
1,04
1,08
1,02
1,02
1,00
0,96
0,98
0,90
0,96
0,84
0,94
[20] r,2 =4. 0,78 [20] r,2=10.
0,92
0,90 2=-1. 0,72
2=-1.
0,66
2= 1.
0,88
0,86 0,60
2= 1.
0,0 0,5 1,0 1,5 2,0 2,5 3,0 0,0 0,5 1,0 1,5 2,0 2,5 3,0
(c) d2/d1 (d) d2/d1
Figure 5.18. The chirality effect of a superstrate on the normalized real resonant
frequency; a=6 cm, b=5 cm, d1=0.1 cm, r,1=2.35, (a) r,2=1.5, (b) r,2=2.35, (c) r,2=4,
(d) r,2=10.
that the negative chirality element is advantageous over the positive one, due
to the possibility of bandwidth cancellation for the case when d2≈0.5d1.
4
3
imaginary frequeny
imaginary frequeny
3
2
2
Normalized
Normalized
1
1
0
0
-1
[20] r,2 =1.5 -1
[20] r,2 =2.35
-2
-2 2=-1. -3 2=-1.
-3
2= 1. -4 2= 1.
0,0 0,5 1,0 1,5 2,0 2,5 3,0 0,0 0,5 1,0 1,5 2,0 2,5 3,0
3 4
imaginary frequeny 3
2
imaginary frequeny
2
Normalized
1 1
Normalized
0 0
-1
-1
-2
[20] r,2 =4. -2 [20] r,2=10.
-3
-3 2=-1. -4 2=-1.
-5
2= 1.
-4
2= 1. -6
0,0 0,5 1,0 1,5 2,0 2,5 3,0 0,0 0,5 1,0 1,5 2,0 2,5 3,0
(c) d2/d1 (d) d2/d1
Figure 5.19. The chirality effect of a superstrate on the normalized imaginary resonant
frequency; a=6 cm, b=5 cm, d1=0.1 cm, r,1=2.35, (a) r,2=1.5, (b) r,2=2.35, (c) r,2=4,
(d) r,2=10.
4
[20] r,2 =1.5
5
[20] r,2=2.35
Normalized BW
Normalized BW
2= 1. 3 2= 1.
2
2
1
1
0 0
0,0 0,5 1,0 1,5 2,0 2,5 3,0 0,0 0,5 1,0 1,5 2,0 2,5 3,0
(a) d2/d1 (b) d2/d1
6,0 11
Normalized BW 5,4
[20] r,2 =4. 10 [20] r,2 =10.
Normalized BW
4,8
2= -1. 9
2= -1.
4,2 8
3,6 2= 1. 7
2= 1.
6
3,0
5
2,4
4
1,8
3
1,2 2
0,6 1
0,0 0
0,0 0,5 1,0 1,5 2,0 2,5 3,0 0,0 0,5 1,0 1,5 2,0 2,5 3,0
Figure 5.20. The chirality effect of a superstrate on the normalized bandwidth; a=6 cm,
b=5 cm, d1=0.1 cm, r,1=2.35, (a) r,2=1.5, (b) r,2=2.35, (c) r,2=4, (d) r,2=10.
1,000
1,00
0,99
0,975
0,98
1 =-1., 2 =-1. 1 =-1., 2 =-1.
0,97 1 =+1., 2 =-1. 0,950 1 =+1., 2 =-1.
0,96 1 =-1., 2 =+1. 1 =-1., 2 =+1.
1 =+1., 2 =+1.
r,2 =1.5 r,2 =2.35
1 =+1., 2 =+1.
0,95 0,925
0,0 0,5 1,0 1,5 2,0 2,5 3,0 0,0 0,5 1,0 1,5 2,0 2,5 3,0
(a) d 2 /d 1 (b) d 2 /d 1
1,00 0,95
0,90
isotropic case [40]
0,95
0,85
1=0.,2=-1. [40]
0,80
1=0.,2=+1. [40]
1=-1., 2=-1. 1=-1., 2=-1.
0,75
0,90 1=+1., 2=-1. 1=+1., 2=-1.
0,70
1=-1., 2=+1. 1=-1., 2=+1.
1=+1., 2=+1.
r,2 =4. 0,65
1=+1., 2=+1.
r,2 =10.
0,85 0,60
0,0 0,5 1,0 1,5 2,0 2,5 3,0 0,0 0,5 1,0 1,5 2,0 2,5 3,0
r,2 =1.5
r,2 =2.
1=-1., 2=-1.
Normalized bandwidth
5 1=-1., 2 =-1.
Normalized bandwidth
4,0
1=+1., 2=-1.
3,5 1=+1., 2 =-1.
1=-1., 2=+1. 4
3,0 1=-1., 2 =+1.
1=+1., 2=+1.
1=+1., 2 =+1.
2,5 3
2,0 isotropic case [40] isotropic case [40]
1=0., 2=-1. [40] 2 1=0., 2 =-1. [40]
1,5
1=0., 2=+1. [40] 1=0., 2 =+1. [40]
1,0
1
0,5
0,0 0
0,0 0,5 1,0 1,5 2,0 2,5 3,0 0,0 0,5 1,0 1,5 2,0 2,5 3,0
(a) d2/d 1 (b) d 2/d 1
Normalized bandwidth
5,25
1 =+1., 2 =-1. 1 =+1., 2 =-1.
10
4,50 1 =-1., 2 =+1. 1 =-1., 2 =+1.
3,75 1 =+1., 2 =+1. 8 1 =+1., 2 =+1.
0,00 0
0,0 0,5 1,0 1,5 2,0 2,5 3,0 0,0 0,5 1,0 1,5 2,0 2,5 3,0
(c) d 2 /d 1 (d) d 2 /d 1
1,05
r,2 =2.35
r,2 =10
1,00 0,9
0,95
0,8
0,90
0,85 0,7
Figure 5.24. Combined effect of substrate and superstrate on the normalized resonant
frequency versus superstrate thickness; a=6 cm, b=5 cm, d1=0.5 cm, t,1=z,1=2.35,
z,2=1.5, 2.35, 4, 10.
1,06
r,2 =1.5
Normalized real frequency
1,05
1,02
1,00
1,00
0,98 0,95
0,96
1=2=0. 0,90
0,94
1=2=0.
0,92 1=-1., 2=-1. 1=0., 2=+1.
0,85 1=-1., 2=-1. 1=0., 2=+1.
0,90 1=+1., 2=-1. 1=0., 2=-1.
1=+1., 2=-1. 1=0., 2=-1.
0,88 0,80
1=-1., 2=+1. 1=-1., 2=+1.
0,86
1=+1., 2=+1. 0,75
1=+1., 2=+1.
0,84
0,0 0,5 1,0 1,5 2,0 2,5 3,0 0,0 0,5 1,0 1,5 2,0 2,5 3,0
r,2 =10.
Normalized real frequency
1,0
1,00
0,9
0,95
0,90 0,8
0,85 0,7
0,80
1=-1., 2=-1.
1=+1., 2=-1.
0,6 1=-1., 2=-1.
0,75
1=-1., 2=+1. 0,5 1=+1., 2=-1.
0,70 1=2=0.
1=2=0. 1=-1., 2=+1.
0,65 1=+1., 2=+1. 0,4
1=0., 2=+1. 1=0., 2=+1.
0,60 1=+1., 2=+1.
1=0., 2=-1. 0,3 1=0., 2=-1.
0,55
0,0 0,5 1,0 1,5 2,0 2,5 3,0 0,0 0,5 1,0 1,5 2,0 2,5 3,0
Figure 5.25. Combined effect of substrate and superstrate on the normalized bandwidth
versus superstrate thickness; a=6 cm, b=5 cm, d1=0.5 cm, t,1=z,1=2.35, z,2=1.5, 2.35,
4, 10.
Effect of Anisotropic Magneto-Chirality ... 191
Figures 5.25(a-d) show that the chirality effect of the substrate on the
bandwidth is symmetric in this case, but the combined effect (ξ1 and ξ2) is
asymmetric on the bandwidth. This is more apparent for low permittivity εr,
important ratio d2/d1 and for a positive magneto-electric element of the
superstrate, ξ2=+1, Figures 5.25(a-b). The asymmetry effect of the superstrate
chirality on the bandwidth becomes important when increasing both the
permittivity εr and the thickness of this layer, with a negative magneto-electric
element of the superstrate ξ2=-1, Figure 5.25(d), whereas in the case of a
positive magneto-electric element of the superstrate, the chirality of the
substrate is negligible.
The symmetry and asymmetry of the bi-isotropic chiral media, having
scalar constitutive parameters, have been reported by some authors [134, 135].
In [134], the authors found that the chiral substrate produces additional
longitudinal asymmetric and transverse symmetric field components, this is
the case of isotropic study. This effect could significantly change the
properties of microwave devices built on a chiral substrate, and with respect to
the asymmetry, it is necessary to evaluate the asymptotic form of the Green
tensor. This form is obtained analytically just below, when d1→ 0, and d2 are
maintained constant to show the asymmetry effect of this chiral layer,
according to the Green tensor given in (5.73) as the following
κe κ ξ
κ02εt,1 κ s2 κ z cos κ z,e 2 d2 j z,2 0e 2 κ z ε t,2 jκ0ξ 2 sin κ z,e 2 d2
ε κ
t,2 z,2
d
G κs 1 e
0
jωε κ z εt,1 cos κ z,2 d2 j κ z,2 sin κ z,2 d2 j
e
e
εt,1 κ 0ξ
εt,2 κ z,e2
2
κ z εt,2 jκ0ξ 2 sin κ z,2 d2
e
0 κ02
(5.87)
It can be noticed from the previous expression that the chirality anomaly is
observed as
κ
G e κ s f e0 κ z ε t,2 ξ 2 j κ 0 ξ 22 sin κ z,e 2 d2 (5.88)
κ z ,2
The term εt,2 ξ2 sin κ z,2 d2 e
is significant for important values of t,2, and
this last result can explain some previous results. The anomaly can be
explained also by the fact that the fields are divided into the symmetric and
asymmetric field components. The symmetric field components were only
192 C. Zebiri, F. Benabdelaziz and M. Lashab
slightly affected by the chirality of the media, but the asymmetric field
components were entirely due to the chirality of the media. This can be clearly
noticed by rewriting the equations, TE and TM immitance in symmetric and
asymmetric forms
g s diag j 2 02 t 0 e 2 j ze 2 2 0 t e 2
j 0 j
h
z
(5.89)
0
z 0 z
0 t
asymmetric part symmetric part asymmetric part
h s diag j 2 02 t 0 e 2 j ze 2 2 0 t e 2
j
0 j z
h
(5.90)
0
0z
z
0 t
asymmetric part symmetric part asymmetric part
CONCLUSION
Numerical results for resonant frequency and bandwidth were presented
for different values of constitutive parameters. The asymmetric effect of a thin
chiral substrate is negligible in the first case examined, a monolayer microstrip
antenna, but the effect of the superstrate is always important. Another original
result was obtained, not previously discussed in the literature, namely the
relationship between the permittivity (exactly t,2) and the chiral constitutive
parameter. It is shown that by increasing the permittivity, the asymmetrical
Effect of Anisotropic Magneto-Chirality ... 193
effect of the chiral becomes high, and for complex values of the magneto-
electric parameters , (in this case a non-diagonal matrix) has shown that the
transverse component occurs in the chiral asymmetry, which is not the case for
bi-isotropic case study where , have a scalar form. In addition to what is
found in the literature, it is shown that the layer plays an important role in the
chiral asymmetry. Bianisotropy and specially the asymmetry of the substrates
should always be taken into account in the design of the microstrip resonators.
Furthermore, to improve the resonator parameters the chiral is probably most
convenient for this task, since this artificial medium can be modelled as
required. The effects of a uniaxial permittivity, uniaxial permeability and
chirality of the superstrate on the complex resonant frequency and bandwidth
have been studied using the integral equation formulation. Fast numerical
convergence is obtained using sinusoidal basis functions to expand the current
on the patch, with the calculation being carried out on a PC with 1.86 GHz
Intel Pentium M and 1,024 GB RAM. The determinant cancellation of the
Green’s tensor is obtained in about 42 s. Numerical results for the resonant
frequency and bandwidth have been presented for various constitutive
parameter configurations. The most important parameters for a substrate-
superstrate are t,2, µz,2, in contrast to the monolayer case where they are z,1,
µt,1, whereas the chirality elements have an important effect for any case. The
advantage of positive uniaxial permittivity anisotropy is summarized by the
fact that the resonant frequency and the bandwidth are subject to an increase.
For all presented behaviour of the chiral medium in a rectangular microstrip
antenna substrate-superstrate configuration, it can be concluded that the
material is suitable for the conception of radiating resonators or selective filter
resonators.
REFERENCES
[1] X. -M. Pan and X. -Q. Sheng, “Improved algebraic preconditioning for
MoM solutions of large-scale electromagnetic problems,” IEEE
Antennas and Wireless Propagation Letters, vol. 13, pp. 106-109, 2014.
[2] M. A. Echeverri Bautista, M. A. Francavilla, F. Vipiana and G. Vecchi,
“A hierarchical fast solver for EFIE-MoM analysis of multiscale
structures at very low frequencies,” IEEE Transactions on Antennas and
Propagation, vol. 62, no. 3, pp. 1523-1528, 2014.
194 C. Zebiri, F. Benabdelaziz and M. Lashab
Chapter 6
INTERACTION OF EM FIELDS
TO THE HUMAN BODY USING HYBRID
COMPUTATIONAL METHOD
ABSTRACT
A new hybrid method of moments (MoM)/finite-difference time-
domain (FDTD), with a subgridded finite-difference time-domain
(SGFDTD) approach is presented. The method overcomes the drawbacks
of homogeneous MoM and FDTD simulations, and so permits accurate
analysis of realistic applications. As a demonstration, it is applied to the
short-range interaction between an inhomogeneous human body and a
small UHF RFID antenna tag, operating at 900 MHz. Near-field and far-
field performance for the antenna are assessed for different placements
over the body. The cumulative distribution function of the radiation
efficiency and the absorbed power are presented and analyzed. The
algorithm has a five-fold speed advantage over fine-gridded FDTD.
208 K. N. Ramli, R. A. Abd-Alhameed and P. S. Excell
6.1. INTRODUCTION
Rigorous determination of electromagnetic fields within arbitrary,
anisotropic, and inhomogeneous dielectric bodies is important for researchers
exploring the effects of microwaves upon living tissue. Electromagnetic fields
are induced inside any biological system, such as the human body, when it is
illuminated by an electromagnetic wave, which is also scattered externally.
The human body is an irregularly shaped heterogeneous conducting medium
whose permittivity and conductivity vary with the frequency of the incident
wave, so the distributions of the internal and scattered electromagnetic fields
depend on the body’s physiological parameters and geometry, as well as on
the polarization and frequency of operation. Due to its mathematical simplicity
and ease of implementation, FDTD [1-12] numerical solutions have been
widely adopted to solve complex stratified dielectric objects, such as human
bodies [13-15] and biological tissues [16-17].
The conventional uniform mesh FDTD becomes unfavourable due to
excessive computation time and memory when it is applied to electrically
small objects which require high spatial resolution. To circumvent this
problem, the non-uniform mesh and subgridding schemes in FDTD were
proposed. However, these methods may produce spurious solutions or suffer
from instability [18], so the FDTD method has been used in conjunction with
other numerical techniques in order to tackle this limitation. In 1993, Aoyagi
et al. [19] used the Yee algorithm with the scalar wave equation to reduce the
computations needed to model a vivaldi antenna, while Cangellaris et al. [20]
used a hybrid spectral-FDTD method to analyze propagation in anisotropic,
inhomogeneous periodic structures. Wang [21] introduced a hybrid ray-FDTD
method and used it to investigate scattering from a cavity with a complex
termination and wave penetration through inhomogeneous walls. In 1994,
Mrozowski [22] introduced a hybrid FDTD-PEE (partial eigenfunction
expansion) method to speed up the FDTD method when solving shielded
structure problems. In addition, finite element and finite volume methods have
recently been combined with FDTD, [23, 24], for accuracy in handling curved
geometries and systems with fine features. Hybrid methods operating entirely
in the time domain have been reported in the literature [25, 26], but the time-
domain MoM is not at the state of maturity and flexibility of the frequency-
domain version. Time-domain MoM does have the advantage of generating
information over a wide frequency band, and does not need an iterative
procedure to couple with FDTD, but it requires very long run-times when
treating a junction with more than two wires [27], unlike the frequency domain
Interaction of EM Fields to the Human Body … 209
(a)
Interaction of EM Fields to the Human Body … 211
(b)
J if nˆ H if
(6.1)
M if E if nˆ
(6.2)
Here n̂ is the outwardly directed unit vector normal to the surface from
the source region. Hif and Eif are equivalent to the forward-scattered magnetic
and electric fields respectively from the source region on the equivalent
surface Sc, and Jif and Mif are the corresponding electric and magnetic source
currents respectively on this surface. These currents are then treated as sources
in the FDTD computational province, propagating fields to the scatterer by
using the E and H curl equations given by the expression:
212 K. N. Ramli, R. A. Abd-Alhameed and P. S. Excell
B
E M if
t (6.3)
D
H J if
t (6.4)
J ib H ib nˆ (6.5)
M ib nˆ E ib (6.6)
where Hib and Eib are the back-scattered fields computed at Sc . Note that n̂ is
as above, directed outwards from the source region. Jib and Mib are the electric
and magnetic equivalent surface currents at S c .
Now, the voltage back-scattered in the source region (the excitation for the
MoM) can be evaluated using either of the following equations, defined by the
reciprocity theorem in the same way as in [33, 37]:
Vb dS
Sa
a ( J ms E ib )
(6.7)
V b J ib E ms M ib H ms , dS c
(6.8)
dS
Vb c ( J ib E ms M ib H ms )
S c
(6.9)
where
Interaction of EM Fields to the Human Body … 213
1
E ib j A( r ) V ( r ) F (r )
(6.10)
A( r ) dS c J ib G ( r , r ' )
S c
(6.11)
V (r )
j
dS c 's J ib G ( r , r ' )
S c
(6.12)
F (r ) dS
c M ib G ( r , r ' )
S c
(6.13)
G(r, r) is the free space Green function, given by the expression:
jk r r '
e
G (r , r ' )
r r'
(6.14)
n
S c
Vb (J
k 1
ib k E ms ( rk , r ' ) M ibk H ms ( rk , r ' )) a k
(6.15)
where rk is the position vector of the centre of the cell surface and ak is the
surface area of the cell surface. Therefore J ib and M ibk are the equivalent
k
214 K. N. Ramli, R. A. Abd-Alhameed and P. S. Excell
surface currents at the centre of the surface cell n. Since the excitation voltages
are known, the MoM can be executed to compute the new currents and the
procedure can be repeated until the steady state solution is achieved.
Figure 6.3. Magnitude of Ey and Ez electric field components along z axis at y = 7.2
cm: Near-field Ey(‘ooo’), Ez(‘xxx’), Far-field Ey and Ez (‘ ’).
216 K. N. Ramli, R. A. Abd-Alhameed and P. S. Excell
The field distribution over an x-z plane 7.2 cm distant from the sensor for
the near-field and far-field technique is shown in Figure 6.4 (a) and (b)
respectively. The plane size considered here was 20 cm 16 cm for the x and z
axes respectively. When the far-field and near-field techniques were checked
for comparison of one antenna geometry i.e. with Dt was set at 16.8 cm for
both techniques, the fields were found identical to each other. Both methods
show good stability and the results were convergent within four iterations.
However, the total field components were found to vary between ±2% when
compared to results computed using NEC software.
(a)
(b)
Figure 6.4. A basic. Distribution of the Ez and Etotal field components in dB at 7.2 cm
away from the sensor, (a) using near field method and (b) using far field method. (left:
Ez, right: Etotal).
To validate the accuracy of the method [40], where Figure 6.5 depicts the
calculation model used, and to compare the result with existing published
results [15, 41], a COST244 model [15, 41] was adopted as a human phantom.
This model represents a homogeneous cubical human head of size 200 200
Interaction of EM Fields to the Human Body … 217
200 mm3, with relative permittivity r = 41.5, conductivity = 0.95 S/m and
volume density = 1000 kg/m3.
A 160 mm long half-wavelength dipole, operating at 900 MHz, was used
as a radiation source, with a continuous wave with an input power of 250 mW
exciting the dipole. The distance between the antenna and the human head was
15 mm and the origin of coordinates was set at the centre of the surface of the
human head model, on the same horizontal level as the antenna feeding point.
The FDTD cell size and time step in the analysis area were 2.5 mm and 3.3 ps
respectively. The problem space, cell sizes and number of PML layers were
(127127127 cells), 2.5 mm and 6, respectively. The size of the equivalent
Huygens surface is 6668 cells, equivalent to 15 mm 15 mm 170 mm.
The subgridded volume is 1616 16 cells, equivalent to 40 mm 40 mm
40 mm. It should be noted that this subgridded volume is placed between two
dielectric media, i.e. free space (40540 mm3) and the human head
(403540 mm3) in order to evaluate the effectiveness of the proposed method
when applied to inhomogeneous tissues. Two and four subgridding factors
were used for each cell of the main FDTD. The dielectric properties for the
subgridded cells located on the boundary were averaged and assigned on each
field component, as shown in Figure 6.6. The variations of the unaveraged
specific absorption rate (SAR) along the y-axis with and without subgridding
are shown in the Figure 6.7.
These unaveraged SAR values without subgridding, with subgridding
factor of 2 and with subgridding factor of 4 were observed at every 1.25 mm,
0.625 mm and 0.3125 mm points respectively in the ascending order along the
y-axis. The outcomes were promising for the case of without subgridding and
both with subgridding factors applied in this example due to the convergence
of the results into one line.
The averaged SAR of 6.8 W/kg was computed by numerical simulation
(in this case, it was averaged) over 10 gm of human tissue. The value of 6.8
W/kg (averaged over 10 gm of tissue) when compared with [15, 41] were very
close. In order to prove the proposed technique is computationally viable, the
situation shown in Figure 6.5 was also modeled using FDTD with coarse and
fine grids, and using FDTD/MoM methods. The required computational times
for each method are compared in Table 6.1.
In the coarse grid FDTD model, a 2.5 mm cell size was used while in the
fine grid FDTD model, 1.25 mm and 0.625 mm cell sizes were adopted. It
should be noted that all simulations were performed on an Intel Core-i7-2000
desktop with 3.4 GHz CPU and 16 GB of RAM.
218 K. N. Ramli, R. A. Abd-Alhameed and P. S. Excell
Figure 6.6. Distribution of field components on the interface boundary model of two
different mediums of the proposed subgridding model.
Interaction of EM Fields to the Human Body … 219
Figure 6.7. The variations of the unaveraged SAR along y-axis of the model shown in
Figure 6.5 [15, 41].
220 K. N. Ramli, R. A. Abd-Alhameed and P. S. Excell
(a) (b)
Figure 6.8. RFID antenna with equivalent Huygens box: (a) Horizontal polarization,
(b) Vertical polarization.
A six cell perfectly matched layer (PML) is used to terminate the FDTD
space, and the distance between the antenna Huygens surface box and the
Interaction of EM Fields to the Human Body … 221
human body is 12.0 mm (or 2 cells). The human body model employed in this
work was developed by Mason et al. in 2000 [43]. The near-field and far-field
radiation of the antenna at 900 MHz for different locations has been analyzed
in order to investigate the performance of the antenna in proximity to the
human body. The number in Figure 6.9 indicates the location of the RFID
antenna.
Figure 6.9. Location of the antennas, represented by black dots, in proximity to the
human body.
Figure 6.10. Subgrid cells of 10×10×10 FDTD cells are taken inside the human body:
(a) Front location, (b) Back location (subgridded).
number of samples. Figure 6.10 shows the location of the antenna in front and
at the back of human body respectively, in which subgrid cells of 10×10×10
FDTD cells were taken inside the human body for near-field analysis. The
near-field and far-field of the antenna are calculated and analyzed at various
locations in order to build up a realistic picture of the antenna performance
while in close proximity with the human body. Figure 6.11 shows the electric
field distributions using an equivalent dB scale in the immediate
neighbourhood and the interior of the body model, for the vertical polarization
state, in x-y, x-z and y-z planes.
Figure 6.12 shows the electric field distribution inside the subgridded
region for x-y, x-z and y-z planes (represented in the equivalent dB scale).
Additional simulations were performed using the same antenna placements for
horizontally polarized antenna, and these were found to be comparable with
the vertical state. The electric field distributions obtained in the neighbourhood
of the human body model were almost the same regardless of whether the
antenna tag was horizontally or vertically polarized. It is interesting to note
that a similar observation was also produced within the subgridded region. The
electric fields were very strong when the antenna was located close to the
body, as represented by the red colouring in the plots.
Figure 6.13 illustrates the far-field radiation pattern for the antenna (in
horizontal polarization state) at the front and back of the human body model. It
should be noted that the underlying computation was normalized to 1 W input
power. The variation in the far-field patterns implies that the field distributions
are more concentrated in the direction facing the normal antenna axis, and
away from the body.
The field magnitude is reduced by between 10 dB to 20 dB, which appears
to be due to the tailing effect of the body. Once more, the field distributions
for vertical polarization were quite similar to the horizontal case, and hence
are not reproduced here.
It can be concluded from Figure 6.14 that the cumulative distribution
function (CDF) has similar curves for horizontal and vertical polarization in
these locations. Moreover, it is apparent that the standard deviation of the
radiation efficiency when the mobile is located in the front is less compared to
the back.
The antenna achieves better radiation efficiency of 43% mean percentage
value for both horizontal and vertical polarization comparing front and back
location.
Interaction of EM Fields to the Human Body … 223
(a)
(b)
Figure 6.11. Electric field distribution (dB scale) for vertical polarized antenna placed
at: (a) Front of the body, (b) Back of the body.
224 K. N. Ramli, R. A. Abd-Alhameed and P. S. Excell
(a)
(b)
Figure 6.12. Electric field distribution (dB scale) inside subgrid region for vertically
polarized antenna placed at: (a) Front of the body, (b) Back of the body.
Interaction of EM Fields to the Human Body … 225
Figure 6.13. Far-field pattern for horizontally polarized antenna placed at: (a) Front of
the body, (b) Back of the body; ‘o-o-o’: E, ‘x-x-x’: E.
a)
b)
Figure 6.14. Cumulative distribution function of radiation efficiency and the ratio of
Pabsorbed / Pradiated for: (a) Horizontal polarized antenna, (b) Vertical polarized antenna.
226 K. N. Ramli, R. A. Abd-Alhameed and P. S. Excell
This follows from a loss in the tissue at the front, and a greater degree of
absorbed power at the back of the human body, as can be seen from the plots.
Taking the previous simulation results into account, we may conclude that
altering the position of the antenna on the front will not result in a large
dispersion in the radiation efficiency of the antenna.
CONCLUSION
A hybrid MoM-FDTD-SGFDTD approach has been presented and
adopted for modelling the interaction of the human body with a short-range
RFID antenna. In order to ensure that optimal quality of results was achieved
using a limited sample size of 32 points, each point was investigated using two
both horizontal and vertical polarizations. The MoM technique was used to
produce the electric and magnetic fields created by the antenna in free space,
and the equivalent electric and magnetic sources on the Huygens surface were
produced at each time step. The FDTD technique was then applied for the
whole structure of the problem combined with subgridding at the object of
interest, namely inside the human body near the excitation source. The human
body model was designed to be inhomogeneous at close proximity to the
antenna. The near-field and far-field distributions were incorporated into the
study to heighten the understanding of the impact on human tissue both facing
the antenna, and not directly facing the antenna. The cumulative distribution
function of the radiation efficiency and the ratio of absorbed to radiated power
of the antenna at these locations were also computed. The results support the
conclusion that there was a clear improvement in the front of the human body
model compared to the back position. The combination of hybrid MoM-
FDTD-SGFDTD method approach with an arbitrary inhomogeneous human
body model encourages the development of this new EM field interaction
modeling approach. This study robustly supports the computing of the total
power dissipated and the SAR distribution inside human tissue.
REFERENCES
[1] R. Xiong, B. Chen, Y. Yi, H. -L. Chen and H. -F. Zhao, “Two-step
method for the FDTD analysis of long apertures having depth,” IEEE
Antennas and Wireless Propagation Letters, vol. 12, pp. 39-42, 2013.
Interaction of EM Fields to the Human Body … 227
Chapter 7
7.1. INTRODUCTION
The problem of modelling electromagnetic field interaction with structures
whose dimensions are a few wavelengths has been the subject of many studies
and numerical simulations for bioelectromagnetic applications. The methods
used for such numerical simulations generally fall into one of three categories,
method of moment (MoM) [1-12], finite-element method (FEM) [13-23] and
finite-difference time-domain (FDTD) [24-36].
Method of moment does an excellent job of analyzing unbounded
radiation problems and they excel at analyzing perfect electric conductor and
homogeneous dielectrics. However, it is not very effective when applied to
arbitrary configurations with complex geometries or inhomogeneous
dielectrics. As for finite-element method, it is a very versatile technique
because it allows the analysis of complex inhomogeneous structures.
However, it requires unaffordable computation resources to solve extremely
large problem, due to discretising the entire problem space and lack of
adequate free-space boundary conditions that can be placed close to the
objects modelled.
In contrast, the finite-difference time-domain, since its first introduction
by Yee [37], has been used extensively over the last decade for
bioelectromagnetic dosimetry – the numerical assessment of electromagnetic
fields coupled to complex heterogeneous biological entities at powerline
frequency and radio frequency. The use of FDTD is very attractive for this
application, due mainly to its algorithmic simplicity, capable to solve
complicated heterogeneous geometries and ease of modularization and
parallelisation of the algorithm on massively parallel computers. Moreover, it
is also capable to model electromagnetic wave interaction problems requiring
the solution of considerably more than 108 field vector unknowns. At this level
of complexity, it is possible to develop detailed, 3-D models of complete
engineering systems.
The classical form of FDTD method requires extremely small time-step
sizes when modelling electrically-small regions (much smaller than a
wavelength). Thus, it can become impractical due to the unaffordable
computation times required. This problem can be solved by implementing a
quasi-static approximate version of FDTD. This approach is based on
transferring the working frequency to a higher frequency, to reduce the
number of time steps required. Then, the generated internal field at the higher
frequency can be scaled back to the frequency of interest [37, 38-41]. It should
be noted that this approach is only valid if the size of the interacting structure
Quasi-Static FDTD Scheme for Electrically-Small Regions … 235
in the problem space is 10 times or more smaller than the wavelength and |σ +
jωε| >> ωεo [35, 38], where σ and ε are the conductivity and permittivity of the
tissues respectively, is the radian frequency, and εo is the permittivity of free
space. In order to prove the validity of the quasi-static approach, this chapter
will develop a 3-D FDTD program to directly model a single homogeneous or
multi-layered sphere inside a lossless or lossy problem space. The sphere was
excited by a plane wave, which is replaced by an equivalent surface. The edge
of the problem space was truncated by a properly modified Berenger’s
perfectly matched layer (PML) absorbing boundary condition (ABC), with
matching impedance condition (σ/ε = σ*/ for lossless or lossy media, whereas
σ/εo = σ*/o for free space) [42-45]. An optimum value of the g factor is
properly chosen to provide efficient low reflection termination of the
computational space for different penetrable media problem at the interface
layer.
The spherical Bessel function of the first kind, Bn(1) = Jn(r), indicates a
standing wave.
The spherical Bessel function of the second kind, Bn(2) = Yn(r),
indicates a standing wave as well and is also called the spherical
Neumann function or Weber function.
The spherical Bessel function of the third kind, Bn(3) = Jn(r) + iYn(r),
indicates an inward traveling wave and is also called the Hankel
function of the first kind, Hn(1)(r).
The spherical Bessel function of the fourth kind, Bn(4) = Jn(r) - iYn(r),
indicates an outward traveling wave and is also called the spherical
Hankel function of second kind, Hn(2)(r).
For the innermost layer of the sphere the spherical Bessel function of first
kind, Jn, is the only admissible solution, while for the surrounding medium the
fields will consist of incident waves and scattered waves in order to satisfy the
radiation condition at infinity. Hence, the general solution of the electric and
magnetic field vectors in pth region of concentric sphere can be formulated as
follows:
Quasi-Static FDTD Scheme for Electrically-Small Regions … 237
2 n 1 p 1
1 3 3
E p E o e i t i n an mo1n ibnp ne1n np mo1n i np ne1n (7.1)
n 1 nn 1
i nn 1 b m
k E eit
2n 1 1 1 3 3
Hp p o n p
ianp no1n np me1n i np no1n (7.2)
p n 1
n e1n
where a, b, , are the expansion coefficients and m , n are the spherical
vector wave functions as given by Stratton [47].
B nz k p r Pn1 cos
z 1 cos ˆ
mo
e 1n sin
sin
(7.3)
B nz k p r
sin
Pn1 cos cos ˆ
z nn 1 z
Bn k p r Pn1 cos cos rˆ
sin
no
e 1n
k pr
1
k p r k p r
k p rBnz k p r
1
sin
Pn cos cos ˆ
(7.4)
1
k p r sin k p r
k p rBnz k p r Pn1 cos
cos ˆ
sin
and electrical conductivity of the pth layer, and the positive value of the square
root is used.
238 C. H. See, A. S. Abdullah, J. M. Noras et al.
The n-th term of Pn1 cos is the associated first order Legendre
polynomials, properties for which can be found in literature [66] or might be
computed using the follow recursive relations:
Pn1 cos 1
Let: Qn and Rn Pn cos then:
sin
Qn
2n 1cosQn 1 nQn 2 and
n 1
Rn n cosQn n 1Qn 1
Figure 7.2. M-layered concentric spherical system. Shaded region represents pth layer
of the system with the radius of rp and electrical parameters of p, p, and p.
E
p E p 1
E
p E p 1 (7.6)
H
p H p 1
H
p H p 1
Substituting m and n functions into equations (7.1) and (7.2) and
applying the four relationships in equation (7.6), the expansion coefficients
may be written as follows [51, 67]:
anp
1
n
p
np jnp p hnp np 1 anp 1 np hnp 1 p hnp np 1 np 1 (7.7a)
240 C. H. See, A. S. Abdullah, J. M. Noras et al.
np
1
pn
p jnp np 1 np jnp 1 anp 1 p jnp np 1 np hnp 1 np 1 (7.7b)
bnp
1
n
p
np np jnp 1 hnp np 1 bnp 1 np np hnp 1 hnp np 1 np 1 (7.7c)
np
1
pn
jnp np 1 p np jnp 1 bnp 1 jnp np 1 p np hnp 1 np 1 (7.7d)
where:
jnp jnp k p rp
hnp hnp k p rp
pn jnp np hnp np
k p 1 p
p
k p p 1
np
1
k p rp k p rp
k p rp jnp
1
2n 1
n 1 jnp1 njnp1
np
1
k prp k prp
k p rp hnp
1
2n 1
n 1hnp1 nhnp1
anp Q11p Q12p anp 1
p (7.8)
n Q p Q p p 1
21 22 n
Quasi-Static FDTD Scheme for Electrically-Small Regions … 241
bnp U11p U12p bnp 1
p (7.9)
n U p U p p 1
21 22 n
From [43], for a M-layered sphere, the above equations can be stated as
follows:
M
a1n T11 T12 an
1 (7.10)
n T T M
21 22 n
M
bn1 V11 V12 bn
1 (7.11)
n V V M
21 22 n
M 1 M 1
where Tij Qijp and Vij U ijp (7.12)
p 1 p 1
Using the fact that n1 n1 0 (i.e., the Hankel function, in which its
divergent at the origin at the innermost sphere, must be absent for this region),
and a1 b1 1 . Equations (7.8) and (7.9) yield:
M M
a1n T11 T12 1
M (7.13)
0
T T
21 22 n
bn1 V11 V12 1
M (7.14)
0
V V
21 22 n
242 C. H. See, A. S. Abdullah, J. M. Noras et al.
anp 1 Q11p 1 Q12p 1 anp
p 1 (7.15)
n Q p 1 Q p 1 p
21 n
22
bnp 1 U11p 1 U12p 1 bnp
p 1 (7.16)
n U p 1 U p 1 p
21 n
22
reflection (matched across a planar boundary) for all normally incident waves
if the matching condition in equation (7.17) is satisfied for that medium [42].
(7.17)
o o
The idea of matching condition in free space for the PML will be extended
to lossless and lossy media with the following modified matching condition.
(7.18)
r o r o
where and * are the electric conductivity and magnetic conductivity of the
medium respectively. o and o are the free space permittivity and permeability
and the relative permittivity and permeability are r and r. In order to reduce
the reflection on the interface layer, an optimum value of the geometric
grading factor g has been selected, by using an empirical expression as in [44]:
c ln g
ln R ( 0 ) (7.19)
2x g N 1
where x is the spatial increment of FDTD mesh, R(0) is the normal reflection
coefficient, N is the number of the cells in the PML thickness, c is the velocity
of EM waves in free space or lossy medium.
demonstrated when eight grid cells are considered between the PML and the
structure.
Figure 7.3. Basic structure of the problem for FDTD computational domain.
Figure 7.4. 3-D view of the problem for FDTD computational domain.
246 C. H. See, A. S. Abdullah, J. M. Noras et al.
Figure 7.5. Huygens surface in free space FDTD computational domain (logarithmic
scale).
Quasi-Static FDTD Scheme for Electrically-Small Regions … 247
Figure 7.6. Ex component at yz-plane for single-layer sphere excited by plane wave of 1
V/m at 20 MHz (logarithmic scale).
Figure 7.7. Electric fields distribution along various central axes for single-layer sphere
excited by plane wave of 1 V/m at 60 Hz in free space.
Quasi-Static FDTD Scheme for Electrically-Small Regions … 249
Figure 7.8. Hy component at xz-plane for single-layer sphere excited by plane wave of
1 V/m at 20 MHz and 60 Hz (logarithmic scale).
Figure 7.9. Magnetic field Hy component distribution along various central axes for
single-layer sphere excited by plane wave of 1 V/m at 60 Hz and 20 MHz in free
space.
250 C. H. See, A. S. Abdullah, J. M. Noras et al.
Figure 7.10. Ex component at yz-plane for three-layers sphere excited by plane wave of
1 V/m at 20 MHz in free space (logarithmic scale).
Figure 7.11. Electric fields distribution along various central axes for three-layer
sphere excited by plane wave of 1 V/m at 60 Hz in free space.
Quasi-Static FDTD Scheme for Electrically-Small Regions … 251
Figure 7.13. Ex component at yz-plane for single-layer sphere excited by plane wave of
1 V/m at 20 MHz in lossless media (logarithmic scale).
Figure 7.14. Electric fields distribution along various central axes for single-layer
sphere excited by plane wave of 1 V/m at 60 Hz in lossless medium (logarithmic
scale).
254 C. H. See, A. S. Abdullah, J. M. Noras et al.
Figure 7.15. Hy component at xz-plane for single-layer sphere excited by plane wave of
1 V/m at 60 Hz and 20 MHz in lossless media (logarithmic scale).
Figure 7.16. Magnetic field distribution along various central axes for single-layer
sphere excited by plane wave of 1 V/m at 60 Hz and 20 MHz in lossless medium.
Quasi-Static FDTD Scheme for Electrically-Small Regions … 255
Figure 7.18. Ex component at yz-plane for double-layers sphere excited by plane wave
of 1 V/m at 30 GHz in lossy medium (logarithmic scale).
Figure 7.19. Electric field (Ex) distribution along various central axes for double-layers
sphere excited by plane wave of 1 V/m at 2450 MHz in lossy medium.
Quasi-Static FDTD Scheme for Electrically-Small Regions … 257
Figure 7.20. Electric field (Ez) distribution along x-axes for double-layer sphere excited
by plane wave of 1 V/m at 2450 MHz in lossy medium.
Figure 7.21. Hy component at xz-plane for double-layers sphere excited by plane wave
of 1 V/m at 2450 MHz and 30 GHz in lossy medium (logarithmic scale).
258 C. H. See, A. S. Abdullah, J. M. Noras et al.
Figure 7.22. Magnetic field distribution along various central axes for double-layer
sphere excited by plane wave of 1 V/m at 2450 MHz and 30 GHz in lossy medium.
CONCLUSION
By implementing the frequency scaling approach, the number of FDTD
time steps can be reduced. The reflection on the interface layers inside the
FDTD computation domain, has also been successfully reduced in lossless and
lossy penetrate media. The accuracy of the FDTD scaling approach with the
models of homogenous and layered spheres in free space, lossless and lossy
media, was verified. The numerical results were in good agreement with the
analytical ones. This lending its support to the validity of using the scaled-
frequency FDTD method to obtain induced electric fields at power-line and
mobile communication frequencies, as long as the quasi-static conditions are
valid.
REFERENCES
[1] Zang SR; Bergmann, JR. “Analysis of omnidirectional dual-reflector
antenna and feeding horn using method of moments,” IEEE
Quasi-Static FDTD Scheme for Electrically-Small Regions … 259
[24] Xiong, R; Chen, B; Yi, Y; Chen, H-L; Zhao, H-F. “Two-step method for
the FDTD analysis of long apertures having depth,” IEEE Antennas and
Wireless Propagation Letters, vol. 12, 39-42, 2013.
[25] Kamepally, SK; Kumar, BP; Paidimarry, CS. “FDTD estimation for
accurate specific absorption rate in a tumor,” Annual International
Conference on Emerging Research Areas and International Conference
on Microelectronics, Communications and Renewable Energy
(AICERA/ICMiCR), 1-5, 2013.
[26] Whittow, WG; Njoku, CC; Vardaxoglou, YC. “Fine scale simulations of
patch antennas with heterogeneous substrates,” USNC-URSI Radio
Science Meeting (Joint with AP-S Symposium), 223, 2013.
[27] Lesina, AC; Vaccari, A; Bozzoli, A. “A modified RC-FDTD algorithm
for plasmonics in Drude dispersive media nanostructures,” 6th European
Conference on Antennas and Propagation (EUCAP), 687-690, 2012.
[28] Noetscher, GM; Xu, Y; Makarov, SN. “Accuracy of point source models
with coincident phase centers in a cubic FDTD grid for arbitrary source
orientation,” IEEE Antennas and Propagation Society International
Symposium (APSURSI), 1-2, 2012.
[29] Nayyeri, V; Soleimani, M; Dehmollaian, M. “Implementation of a
PEMC boundary condition in the 2-D FDTD technique,” IEEE Antennas
and Propagation Society International Symposium (APSURSI), 1-2,
2012.
[30] Liu, D-C; Chang, H-C. “The second-order condition of FDTD method at
sloped dielectric interfaces by averaging the permittivities along surface
normal,” IEEE Transactions on Antennas and Propagation, vol. 60, no.
11, 5259-5267, 2012.
[31] Weiwei, M; Dong, S; Xianliang, W. “UPML-FDTD parallel computing
on GPU,” International Conference on Microwave and Millimeter Wave
Technology (ICMMT), vol. 4, 1-4, 2012.
[32] Guiffaut, C; Reineix, A; Pecqueux, B. “New oblique thin wire formalism
in the FDTD method with multiwire junctions,” IEEE Transactions on
Antennas and Propagation, vol. 60, no. 3, 1458-1466, 2012.
[33] Shyroki, DM. “Modeling of sloped interfaces on a Yee grid,” IEEE
Transactions on Antennas and Propagation, vol. 59, no. 9, 3290-3295,
2011.
[34] Kunz, KS; Lee, KM. “A three-dimensional finite-difference solution of
the external response of an aircraft to a complex transient EM
environment I: The method, and its implementation,” IEEE Transaction
on Electromagnetic Compatibility, vol. 20, 328-333, 1978.
262 C. H. See, A. S. Abdullah, J. M. Noras et al.
Chapter 8
COMPUTATION OF ELECTROMAGNETIC
FIELD INSIDE A TISSUE USING
QUASI-STATIC AND LUMPED-ELEMENT
FDTD SCHEME
8.1. INTRODUCTION
Research into possible mechanisms of interaction of electromagnetic (EM)
fields with biological tissues and cells in culture has motivated a growing need
for accurate models describing the EM behaviour of cells exposed to these
fields. Therefore, several numerical models have been set up in order to study
the interaction between EM fields and biological entities, at tissue level, cell
level and ionic level. In this area, the most frequently used technique for
computing the EM field is the finite-difference time-domain (FDTD) method
[1-12], due to its independence from the material parameters.
The standard FDTD method requires extremely small time-step when
modelling electrically-small regions (much smaller than a wavelength): this is
especially the case when modelling biological cells, since they have maximum
dimensions of a few tens of micrometers. Thus, it can become impractical due
to the unaffordable computation times required. This problem can be solved by
implementing a quasi-static approximate version of FDTD. This approach is
based on transferring the working frequency to a higher frequency, to reduce
the number of time steps required. Then, the generated internal field at the
higher frequency can be scaled back to the frequency of interest [13-16].
Cells are surrounded by thin membranes, typically about a few
nanometres thick [17]. They are the major barrier in the cell, separating the
inside of the cell from the exterior medium. This structure will allow cells to
selectively interact with their environment. Therefore, the cell membrane has
been identified as the primary target of the action of the EM field on biological
structures. Since the thickness of the membrane is about 1000 times smaller
than the biological cell, and if the standard FDTD procedure were to be blindly
applied to model detail in the membrane within a complete cell model, then
this will take some millions of iterations to complete one cycle of simulation.
This again will cause excessive computation time. To overcome this problem
in standard FDTD, the lumped element finite-different time-domain (LE-
FDTD) method [18-22] was implemented in such a way to model the
behaviour of the membrane, based on the Hodgkin-Huxley (HH) model [23-
27] on the surface of the biological cell.
This chapter presents a new approach to model and analyse of the HH
membrane model [23-27]. The HH model is represented as an electrical circuit
on the surface of the biological equivalent spherical cells. For the sake of
simplicity, the analysed structure has been represented with spherical cells and
then Floquet periodic boundary conditions [28-31] have been applied to the
border of the analysed structure in order to mimic the presence of the
Computation of Electromagnetic Field Inside a Tissue … 267
surrounding cells. Despite the cellular tissues are not perfectly periodic and the
living cells are not precisely spheres, this approximation allows the modelling
of biological tissue with only a small part of the structure and alleviates the
problem of huge requirement of computer resources for the simulation of a
complete tissue. Since the external medium of the biological tissue is lossy
fluid, the modified Berenger’s perfectly matched layer (PML) absorbing
boundary condition (ABC) is used to truncate the computation grid [32-35], in
order to reduce the reflections on the interface layers. It should be noted that
Berenger’s PML ABC is more accurate than the Mur ABC [36], used in other
recent work [16].
A further difficulty is the limited extent of studies on the dielectric
properties of cell tissues [37]; thus, the complex permittivity of each cell tissue
is not clearly established for radio frequencies. However, in this study, an
analytical method for estimating the electrical properties of cell tissues in the
RF band [38] will be adopted throughout the analysis. Earlier work only
considered two media (water and membrane) [16], but the procedure adopted
here enables the tissue model to consist of three media (lossy medium,
membrane and cytoplasm). In addition, a mass of connected biological tissue
is simulated by creating an equivalent stack of compact cell structures (such as
spherical and cylindrical with interstices, and cubical). The total electric fields
along the central axes of rows of these spherical, cubical and cylindrical cell
structures will be investigated.
Ey) which are not on the edge of the surface will be firstly to be discussed, and
subsequently, the updating equations of the edged tangential electric
component (Ez) are demonstrated. Consider the four surfaces planes at i = i0, iN
and j = j0, jM, in which the Floquet periodic boundary condition should be
applied as shown in Figure 8.1.
Figure 8.1. Geometry used in the analysis of 3-D infinite periodic structure illuminated
by a normal incident plane wave.
(a)
(b)
Figure 8.2. (a) Location of Ez() and Ey() components in plane i = i0 and iN, (b)
Location of Ez() and Ex() components in plane j = j0 and jM.
270 C. H. See, R. A. Abd-Alhameed, M. J. Ngala et al.
From Figure 8.2(a), the 3-D updating equations for the tangential electric
components (Ey and Ez) which are not located on the edge of surface plane i =
i0 and iN, can be written as:
n 12
i, j, k p 1 / 2
1
n
n 1 t H x i , j , k 1 / 2 H x
2
E y ( i , j , k ) E y (i , j , k )
n
n
1
n
1
(7.1a)
H z i 1 / 2, j , k H z i 1 / 2, j , k
2 2
n 12
t H y i 1 / 2, j, k H yn 1 / 2 iN 1 / 2, j, k
Ez i, j, k Ez i, j, k
n 1 n
n
1
n
1
(8.2a)
H x i, j 1 / 2, k H x i, j 1 / 2, k
2 2
n 12 n
1
t H z i, j 1 / 2, k H z i, jM 1/ 2, k
2
E n 1
i, j, k E (i, j, k )
n
x x
n
1
n
1
(8.3a)
H y i, j, k 1/ 2 H y i, j, k 1 / 2
2 2
n 12
t H y i 1/ 2, j, k H y i 1/ 2, j, k (8.4a)
n 1 / 2
En 1
i, j, k E i, j, k
n
z z
n
1
n
1
H x i, jM 1/ 2, k H x i, j 1/ 2, k
2 2
As for the case of the edged tangential electric components (Ex, Ey, Ez) on
plane i = i0, iN and j = j0, jM and due to the Ex and Ey components are on the
edge of the absorbing boundary conditions as seen in Figure 8.2(a) and (b),
therefore, they are assumed to be updated by the ABC updating equations. It
should be noted that the only edged Ez tangential components will be
considered here for the periodic boundary condition. According to Figure 8.3,
due to the normal incidence plane wave, hence, four equations of Ez tangential
components can be updated simultaneously:
n 12
t H y i 1/ 2, j , k H yn1 / 2 iN 1 / 2, j , k
E z i, j , k Ez i, j , k
n1 n
(8.5a)
n
1
n
1
H x i, jM 1/ 2, k H x i, j 1/ 2, k
2 2
Cells are surrounded by thin membranes. It is the major barrier in the cell,
separating cell contents from outside. Since the cell is separated from its
environment, the membrane must be able to accommodate the required cell
functions. Therefore, the membrane has to act as a selective barrier, allowing
nutrients to pass in but keeping out many harmful substances to the cell, and as
a dynamic barrier medium, constantly adapting to changing environmental
conditions (different concentrations of ions).
The dimensions of a biological cell are around a few tens of micrometres
while the thickness of the membranes is in the scale of a few nanometres,
strongly depending on the type of the tissue. Depending on the type of the cell,
voltages in the range of 20-200 mV can be obtained across the membrane.
When the cell is in a resting state, the current across the membrane averages
zero, but more generally it depends on the variation on the membrane voltage
[23, 24].
Hodgkin and Huxley (HH) gave a general description of the time course
of the current which flows through the membrane of a squid giant axon when
the potential difference across the membrane is suddenly changed from its
steady state [23, 24]. The results in [23], suggest that the behaviour of a
membrane may be represented by the electrical circuit shown in Figure 8.4.
Current can be carried through the membrane either by charging the
membrane capacitance or by movement of ions through the nonlinear
conductance in parallel with the membrane capacitance. The equations that
describe the HH model are [23]:
dVm
I Cm g k n 4[V , t ](Vm Ek )
dt
g Na m3[V , t ]h[V , t ] (Vm ENa ) (8.6)
gl (Vm EL )
dn [V , t ]
a n (V )(1 n[V , t ]) bn (V ) n[V , t ] (8.7)
dt
dm [V , t ]
a m (V )(1 m[V , t ]) bm (V ) m[V , t ] (8.8)
dt
Computation of Electromagnetic Field Inside a Tissue … 273
dh [V , t ]
a h (V )(1 h[V , t ]) bh (V ) h[V , t ] (8.9)
dt
where:
V Vm Vr (8.10)
The element definitions of the HH model of the above equations are given
as follows:
Table 8.1 quotes the constant value of all the components in HH model, as
stated in [23]. The percentage of the open potassium channels is n4 and the
percentage of the open sodium channels is m3h. It also needs to be highlighted
that the quantities appearing in equations (8.7) to (8.9) are defined as [23]:
Cm Ek gk ENa gNa EL gL
0.01 F/m2 72 mV 360 S/m2 -55 mV 1200 S/m2 50 mV 3 S/m2
274 C. H. See, R. A. Abd-Alhameed, M. J. Ngala et al.
0 . 01 (V 10 )
a n (V ) (8.11)
e (( V 10 ) / 10 ) 1
0 .1(V 25 )
a m (V ) (8.12)
e (( V 25 ) / 10 ) 1
1
bh (V ) ((V 30 ) / 10 )
(8.16)
e 1
d
( a b ) (8.17)
dt
where:
n , m , h
After solving the first order partial differential equation of (8.17), the
expression for ψ can be obtained:
a e
( a b ) t
(t ) 1 (8.18)
a b K
where K is the constant that defines the initial condition. As can be observed,
if the time t is much greater than the time constant τ = 1/(aψ+bψ), the parameter
ψ can be considered and replaced to the following constant:
Computation of Electromagnetic Field Inside a Tissue … 275
a
(t ) t (8.19)
a b
n m h
0.3177 0.0529 0.5961
cell’s membrane is normal to the x-axis. Therefore, the following equations are
exploited for this analysis in FDTD computational domain.
t
(8.22)
Cbm,n, p i, j ,k m
i , j ,k
i, j,k t Cm g n4 t g m3hmt g t
1 k m Na l m
2 pi, j ,k 2n pi, j ,k 2n pi, j ,k 2n pi, j ,k
i , j ,k n
1
(8.23)
Ccm,n, p n p
i , j ,k
i, j,k t Cm gk n4mt gNam3hmt gl mt
1
2 n pi, j,k 2n pi, j,k 2n pi, j,k 2n pi, j ,k
i , j ,k
Figure 8.6. Electric field along the centre of the lossy medium.
GK n 4 Ek GNa m3hE Na GL EL
Vr 60.27 mV
GK n 4 GNa m3h GL
Figure 8.7 shows the comparison between the FDTD model and the
analytical computed membrane voltage versus time. As can be seen, the
analytical and FDTD results are indistinguishable in the steady state. Since the
HH model contains a large capacitor, it causes the transient of the HH circuit
to require hundreds of thousands of time steps to reach the steady state when
simulation considers the dimensions of a realistic biological cell. To speed up
the computation time, the actual FDTD simulation is performed in two steps.
In the first step (known as DC simulation), simulation is performed with the
capacitor in the circuit and then removed and its effects on the system is
energized by the step voltage (DC) sources. Once the steady state is reached,
the electric and magnetic fields on each FDTD discretisation cell are used as
the initial values for the second step. During the transient simulation, the
capacitor and the sinusoidal source are activated. In short, to speed up the
computation time, the value obtained in the steady state without excitation can
be used as the initial value of the membrane voltage. This procedure is similar
to the approach used by the SPICE computer program [24] for the transient
and AC analysis of electronic circuits.
Computation of Electromagnetic Field Inside a Tissue … 279
Figure 8.7. Steady state membrane’s voltage verse number of time steps (the normal is
directed along X).
In the second test case, it was shown that the method above gives the same
results as when complete transient simulation is performed. The HH model
was again inserted in one of the FDTD grid cells, with dimension of 10 cm.
Then, it was excited by a plane wave at a frequency of 10 MHz. Two initial
conditions were assumed in the simulation; firstly the membrane voltage was
set at zero, and then at 60.27 mV. Figure 8.8 depicts the comparison of the two
initial conditions, leading to the conclusion that the computation time can be
speeded up by assigning the resting voltage on the cell’s membrane at the
initial state of the simulation. In a third test case, all the parameters in equation
(8.11) were set to zero, except for n = 1 and gk = 2,000,000 S/m2, so the HH
model behaves like a resistor whose value is:
1
R 50
gk s 2
In this case, the grid cell’s size is 0.1 mm, and a plane wave excitation was
used propagating in the z-direction and polarised in the x-direction at 40 GHz.
Figure 8.9 shows the distribution of voltage and current on the membrane. The
ratio between the peak values of voltage and current on the membrane was
50 Ω.
280 C. H. See, R. A. Abd-Alhameed, M. J. Ngala et al.
Figure 8.8. Comparison between membrane voltage with different initial conditions.
GHz and in the other there was no excitation. It should be noted that the plane
wave in the first case is propagating along the z-axis and electric field is
parallel to the x-axis. Figure 8.13 shows the numerical noise in the centre of
the structure without excitation. Figure 8.14 shows that the dominant
component (Ex) is corrupted by the numerical noise along the z-axis of the
structure, whereas Figure 8.15 shows that the noise is deterministic and the
harmonic excitation can be retrieved. This can be done by subtracting the
numerical noise in Figure 8.13 from the excitation field in Figure 8.14.
Figure 8.12. Total electric field along central axis of the structure.
Figure 8.13. The field versus time in the central FDTD cell of the structure in lossy
medium without excitation.
284 C. H. See, R. A. Abd-Alhameed, M. J. Ngala et al.
Figure 8.14. The corrupted excitation field in the centre of the FDTD cell of the
structure.
In general, the radius of the each biological cell in a tissue was 10 µm.
The biological cell consists of three layers, which are cytoplasm, membrane
and extracellular medium and the dielectric properties of these were obtained
from [38]. The equations of the effective conductivity and effective
dielectric permittivity of the material of biological cells are given in the
following two forms:
n
k k
0 (8.24)
k 1 1 k
2 2
n
k 2 k2
0 (8.25)
k 1 1 k
2 2
t0
0 Tr
0.51 cos t sin t 0t
2
f (t ) r r
(8.26)
sinr t t
Tr
2
where Tr is the period of the ramped cosine which is about 3 source cycles.
The PML, shown in Figure 8.19, was six FDTD elements wide, the
grading factor g was 10.1383 and the grid structure was effectively extended
to infinity in the x- and y-directions, by imposing the Floquet boundary
condition along the x- and y-axes. The Floquet periodic boundary condition
plays an important role to mimic the presence of an extended 3-D structure of
biological cells, simulating connected tissue. This can be easily imagined in
two dimensions, as shown in Figure 8.18. The FDTD problem space was
220×20×20 FDTD elements of size 1 µm while a discretization time step δt of
1.3 fs was chosen to drive the FDTD computation, to meet the requirements of
the Courant stability criterion.
Before implementing the Hodgkin-Huxley model into the simulated
structure, the effect of moving the Floquet boundaries gradually away from the
simulated structure was studied. Figures 8.20 and 8.21 depict the field
distribution through the centre of the simulated structure at 10 GHz with
varying locations of the Floquet boundaries, where Ncell is the number of
FDTD elements between the Floquet boundaries and the boundaries of the
biological cells, in the x and y directions. Figure 8.22 shows the field
distribution on the xz-plane of the simulated structure for the case of Ncell =
10. When the Floquet boundaries are exactly adjacent to the simulated
structure (Ncell = 0), the strongest coupling effect between cells can be
Computation of Electromagnetic Field Inside a Tissue … 289
obtained. The highest induced field on the membrane and lowest induced field
in the cytoplasm of the cell can be observed. Conversely, when the Floquet
boundaries are far away from the simulated structure (Ncell = 10), the lowest
induced field on the membrane and highest induced field in the cytoplasm of
the cell are observed. It should be noted that all the following analysis will be
based on Ncell = 0, which is assumed to be the most appropriate model for the
real living biological tissues or cells in this microdosimetry study. The
simulations were performed at the transformed intermediate frequency of 10
GHz and the overall model was then transformed to the intended lower
frequencies. Table 8.6 reports the transformation factors at 900 MHz, 1800
MHz, 2000 MHz and 2450 MHz used in the analysis.
Figure 8.23 illustrates the 10 GHz field distribution on the xz-plane of the
simulated structure. The distributions of the electric field through the centre of
the simulated structure, along the incident wave propagation direction, at 900
MHz, 1800 MHz, 2000 MHz and 2450 MHz are given in Figures 8.24 and
8.25, where Figure 8.25 is an enlarged version of Figure 8.24. The magnified
version of Figure 8.25 is also plotted in Figure 8.26 and 8.27 as well. From
inspection of Figures 8.25, 8.26 and 8.27, the field inside the cells is not
constant and the induced field intensity is directly proportional to the
frequency. In other words, the higher the operating frequency that is used to
excite the model, the higher the electric field intensity that will be induced
within the analyzed structures.
To complete the simulation, the Hodgkin-Huxley models were embedded
in the surface of the spherical cells, in a direction normal to the surface, to
represent the membrane effect of the tissue model. Versions including this
were studied at frequencies of 900 MHz and 2450 MHz. As can be seen in
Figures 8.28 and 8.29, there is a difference of approximately 15% in the field
strength due to the contribution of the membrane effect from the Hodgkin-
Huxley model. These variations were in well agreement with expectations [16,
23, 38,46,47,48].
Since living cells, when compacted into connected tissue, are not perfectly
spheres, a cluster of cubical cells was also chosen for study on the foundation
of the previous spherical cells analysis. Figure 8.30 depicts the proposed
cluster of cubical cells in a three dimensional view of the FDTD computational
domain. In order to compare the results obtained from the previous model with
290 C. H. See, R. A. Abd-Alhameed, M. J. Ngala et al.
this analysis, an FDTD simulation was executed, keeping the same parameter
values unchanged. The 2-D view of the electric field inside the cubical cell
tissue is shown in Figure 8.31. The field distributions along the propagation
direction of the incident wave, through the centre of the simulated structure at
various frequencies, are illustrated in Figures 8.32 and 8.33, while the
magnified version of Figure 8.33 are depicted in Figures 8.34 and 8.35. The
contribution of the Hodgkin-Huxley model to the cubical tissue model has also
been investigated, as shown in Figures 8.36 and 8.37. The effects of adding the
Hodgkin-Huxley model are about 15% difference in field magnitude, as can be
easily observed in the figures.
The peak field values on the membrane of the cubical structure are
observed to be about three times higher than in the cytoplasm, which agrees
well with the results from the structure based on spherical cells. However, the
absolute field strength is approximately doubled in the spherical cell case,
presumably because of the curvature at the points studied. It is to be expected
that much higher fields would be observed at the corners of the cubical cells,
but it might be argued that, as a localised matter, these points do not
correspond well with biological reality.
Figure 8.18. The 2-D view of the simulated structures in FDTD computational domain.
Computation of Electromagnetic Field Inside a Tissue … 291
The shape of the living cells can be so diverse. In order to have better
understanding on the subject of EM field interaction with different geometry
of biological tissue, a cluster of cylindrical cells model of the tissue is
proposed, as depicted in Figure 8.38. This analysis is performed in which the
material properties unchanged as in the case of spherical and cubical cell
structures. Figure 8.39 describes the 2-D view of electric field distribution of
the proposed cylindrical cell tissue at 10 GHz, while the electric field
distribution along the centre of the analysed structure is shown in Figure 8.40
to 8.43, where the Figures 8.41-8.43 are the amplified version of the Figure
8.40.
The loading effect of the HH model into the cylindrical cell tissue has also
been studied. Figure 8.44 and 8.45 demonstrate the difference of 15% in field
magnitude with and without the presence of HH model in the proposed
simulated structure. The results show consistent difference with the previous
spherical and cubical cells tissue simulated structures.
Figure 8.19. The 3-D view of the simulated spherical structures in FDTD
computational domain.
292 C. H. See, R. A. Abd-Alhameed, M. J. Ngala et al.
Figure 8.20. Penetration of electric field along z-axis, through the centre of the
simulated structure with different location of the Floquet boundary condition.
structure is found to be about 1.7 times higher than in the cytoplasm, which is
distinct with the previous two models. As can be noticed, the peak field value
of this cylindrical model is higher than the cubical model and lower than the
spherical model, whereas the peak field value on the cytoplasm is the about the
same as in the spherical model and double the value found in the cubical
structure.
Figure 8.22. Ncell = 20, modulus of the electric field on xz-plane at intermediate
frequency 10 GHz (logarithmic scale).
Figure 8.24. Penetration of electric field along z-axis, through the centre of the
simulated structure.
Figure 8.28. Electric field distribution along z-axis, through the centre of the simulated
spherical structure in Figure 8.19, incorporating Hodgkin-Huxley model and driven at
900 MHz.
Figure 8.30. The 3-D view of the simulated cubical structures in FDTD computational
domain.
Figure 8.32. Penetration of electric field along z-axis, through the centre of the
simulated structure.
298 C. H. See, R. A. Abd-Alhameed, M. J. Ngala et al.
Figure 8.36. Electric field distribution along z-axis, through the centre of the simulated
cubical cell structure in Figure 8.30, incorporating Hodgkin-Huxley model and driven
at 900 MHz.
300 C. H. See, R. A. Abd-Alhameed, M. J. Ngala et al.
Figure 8.38. The 3-D view of the simulated cylindrical structures in FDTD
computational domain.
Computation of Electromagnetic Field Inside a Tissue … 301
Figure 8.40. Penetration of electric field along z-axis, through the centre of the
simulated structure.
Figure 8.44. Electric field distribution along z-axis, through the centre of the simulated
cubical cell structure in Figure 8.38, incorporating Hodgkin-Huxley model and driven
at 900 MHz.
Computation of Electromagnetic Field Inside a Tissue … 303
CONCLUSION
An approach to microdosimetric modeling of bioelectromagnetic
interactions at the cellular level has been presented. This uses the FDTD
method, combined with an arbitrarily-oriented implementation of the
304 C. H. See, R. A. Abd-Alhameed, M. J. Ngala et al.
REFERENCES
[1] Xiong, R; Chen, B; Yi, Y; Chen, H-L; Zhao, H-F. “Two-step method for
the FDTD analysis of long apertures having depth,” IEEE Antennas and
Wireless Propagation Letters, vol. 12, 39-42, 2013.
[2] Kamepally, SK; Kumar, BP; Paidimarry, CS. “FDTD estimation for
accurate specific absorption rate in a tumor,” Annual International
Conference on Emerging Research Areas and International Conference
on Microelectronics, Communications and Renewable Energy
(AICERA/ICMiCR), 1-5, 2013.
[3] Whittow, WG; Njoku, CC; Vardaxoglou, YC. “Fine scale simulations of
patch antennas with heterogeneous substrates,” USNC-URSI Radio
Science Meeting (Joint with AP-S Symposium), 223, 2013.
Computation of Electromagnetic Field Inside a Tissue … 305
Chapter 9
M. A. Mangoud1, R. A. Abd-Alhameed1
and P. S. Excell2
1
Mobile and Satellite Communications Research Centre,
Bradford University, Bradford, UK
2
Institute for Arts, Science and Technology,
Glyndwr University, Wrexham, UK
ABSTRACT
A heterogeneous hybrid computational electromagnetic method is
presented in this chapter. Different parts of an antenna simulation
problem are treated by using different methods enabling the most
appropriate method to be used for each part. The method employs the
method of moments and finite-difference time-domain techniques to
compute the fields in two regions. The two regions are interfaced by
surfaces on which effective sources are defined by the application of the
equivalence principle. An extension to this permits the conduction
currents to cross the boundary between the different computational
domains. Several examples are studied and the results are compared and
analysed. The method paves the way for accurate simulation of the
behaviour of an antenna that is closely coupled with lossy dielectric
volumes.
310 M. A. Mangoud, R. A. Abd-Alhameed and P. S. Excell
9.1. INTRODUCTION
Electromagnetic computation algorithms are widely used for the
prediction of the behaviour of antennas of all types. In general, these are
broadly classified into two types, integral equation methods, often known as
the method of moments (MoM) [1-10] and differential equation methods (such
as finite-difference time-domain) [11-20], which have significantly different
strengths and weaknesses. Integral equation methods are widely used for
simulation of broadcasting antennas since they can simulate the behaviour of
arbitrarily shaped metallic structures with relative ease and a good level of
fidelity. However, problems arise if the antenna is associated with substantial
dielectric volumes, particularly lossy dielectrics such as earth, sea water and
the human body. The integral equation methods inherently compute the
coupling between each discretization element in the model and every other
element. This does not cause a problem for metal because it is considered to be
effectively impenetrable and hence only the surface has to be discretized;
further, wires and thin rods can be represented as one-dimensional elements.
Penetrable dielectrics, on the other hand, have to be represented with a 3-D
assembly of volumetric elements. Leaving aside the mathematical complexity
of modelling such elements, the key difficulty that arises is that, in being a 3-D
structure, the number of elements required in the model rapidly becomes very
large and the computational task can become unmanageable. An alternative
approach for penetrable dielectrics that is computationally much more
tractable is to use a differential equation method, such as the FDTD method.
This is a well-known and highly developed algorithm that directly implements
a discretized form of Maxwell’s equations. It is very efficient and relatively
easy to implement, but it has one major disadvantage in that it requires the
problem space to be discretized into three-dimensional rectangular elements.
This results in a significant deficiency in the modelling of structures with
surfaces that are not parallel to the principal axes, since the surfaces then have
a stepped shape (known as ‘staircasing’) as a result of conformance to the
rectangular quantization grid. This poses particular problems when attempts
are made to model linear (wire-type) antennas with components that are not
parallel to the axes. These two methods have unique advantages in their own
domain of application but it is not possible to apply either to the other domain
without significant disadvantages (e.g., using integral equation methods to
model the ground or using FDTD to model wires). A method was thus
developed to enable both of them to be used for a single problem, combining
them via a hybridization algorithm in which a bounding surface was
Simulation of Antennas Coupled to Lossy Dielectric Volumes 311
Figure 9.1(a) shows two regions, one enclosing the source, the other
enclosing the scatterer. The source region is further subdivided into two sub-
regions, A and B as shown in Figure 9.1(b). In general there is no physical
attachment between source region B and the scatterer region, but the scatterer
312 M. A. Mangoud, R. A. Abd-Alhameed and P. S. Excell
can be very close to, or attached to, the source region A. In Figure 9.1(b) the
surface Sc encloses the entire source region. The fields due to the induced
currents from both source sub-regions A and B can be computed in order to
evaluate the surface currents Jsi and Msi over the modified closed surface Sc1:
M si ( E ( J A , M A ) E ( J B , M B )) nˆ (9.1)
J si ( H ( J A , M A ) H ( J B , M B )) nˆ (9.2)
where JA, JB, MA and MB are the electric and magnetic surface currents of
source regions A and B; n̂ is the unit vector normal to the surface Sc1.
JA
E E FDTD (9.3)
( / t / 2)
Simulation of Antennas Coupled to Lossy Dielectric Volumes 313
J si Sc1
E EFDTD (9.4)
( / t / 2)
t
H H FDTD M si Sc1 (9.5)
where Jsi and Msi are the surface electric and magnetic currents, and with the
suffix Sc1 are only those defined at the surface Sc1. is the permeability of the
medium. The voltage induced on the source region A due to the back-scattered
fields from the scatterer is:
where Efree and Escat are the induced electric fields in the region A, calculated
in the FDTD domain, without and with the scatterer, and s is the FDTD cell
size. The voltage induced in region B (the scattered field region) is found using
MoM by employing the Reaction Theorem:
VB E ts J ib H ts M ib dsc1 (9.8)
S c1
where:
1
EiB jA(r ) V (r ) F (r ) (9.9)
314 M. A. Mangoud, R. A. Abd-Alhameed and P. S. Excell
j
V (r )
J
s ib g ( r , r )dsc1 (9.11)
S c1
F (r ) M ib g (r , r )dsc1 (9.12)
S c1
jk r r
e
g (r , r ) (9.13)
4 r r
where rn is the position vector of the centre of the cell surface and an is the
surface area of the cell. Therefore Jibn and Mibn are considered to be the
equivalent surface currents at the centre of the surface cell n. At the boundary
of the hard source and the equivalent surface (Sc1) FDTD requires at least two
or three cells of additional margin to recalculate the voltage induced in region
Simulation of Antennas Coupled to Lossy Dielectric Volumes 315
B, as seen from Figure 9.1(b). When the voltages have been obtained, the new
currents (JA, Jsi and Msi) can be computed and the method can be repeated until
the steady state is reached.
the hybrid MoM/FDTD [25] method (with a Huygens surface surrounding the
full dipole) and the pure MoM NEC programs [27].
Example 9.2: A wire scatterer was introduced adjacent to the half-
wavelength dipole of the previous example, in the form of a PEC rod of length
0.3528 and radius 0.001, parallel to the source dipole at a distance 0.0196,
equivalent to two FDTD cells, see Figure 9.2(a). The scatterer was located in
the region of the impressed current source. The dipole here was modelled
using the current crossing technique, using Huygens surface to enclose 19.6%
(0.098) of the length of the dipole. The impressed currents were used for the
rest of the dipole inside the FDTD part of the program. The input impedance
of the source dipole versus the number of iterations is shown in Figure 9.4.
The results are in good agreement with those from a uniform MoM treatment
[27]. The electric field contours in the y-z plane, which contains the source and
the scatterer, are shown in Figure 9.5 for the cases with and without the
scatterer. This figure clarifies the methodology used and the technique of the
current crossing the boundary using the equivalent surface principle.
Example 9.3: In this example, illustrated in Figure 9.2(b), the source
region contains a dipole with total conductor length 1.0 and radius 0.001.
The working frequency was doubled to become 1923 MHz which means
wavelength of 0.156 m. Similar to the previous example, cell size of 0.00306
m was used. One of dipole arms is coiled into a helix with three turns of radius
0.12 and pitch 0.004. This simulates a canonical form of a transmitter
cabinet directly feeding a helical antenna adjacent to a scatterer, in a case
where the cabinet is intended to be simulated by FDTD. The scatterer is a wire
of length 0.473 (24 cells) and radius 0.001. The scatterer and the source are
separated by 0.01182 (6 cells). The equivalent surface encloses the helix
region (treated by MoM) while the impressed current treatment (FDTD region)
is used for the other arm of the dipole. The input impedance versus the number
of iterations is shown in Figure 9.6. The results converge to Z = 420 - j190
which is also in good agreement with a non-hybrid NEC treatment Z = 462 -
j220 [27]. The accuracy of this proposed modified hybrid method is
demonstrated in Figure 9.7. This shows the computed near-field magnitude
distribution contours for the example given in Figure 9.2, when conduction
currents are crossing domain boundaries for both the dipole and the helix-
dipole in free space and with a wire scatterer. The plane displayed contains the
scatterer and the second arm of the source dipole that is represented by the
straight wire. The near-field distribution shows the important advantage of the
Simulation of Antennas Coupled to Lossy Dielectric Volumes 317
method over the conventional hybrid technique [24, 25, 28-30] which requires
a more substantial separation distance between the source and the scatterer.
In the free space simulation of Figure 9.7(a) the symmetric patterns of the
half-wavelength dipole in the total field region outside the Huygens surface
can be observed. Figure 9.7(a) shows a reasonable pattern with expected
asymmetry for the helix-dipole combination of length 1.0, with a minimum
near-field value around the middle of the straight half of the dipole. In both
Figures 9.5(a) and 9.7(a), the near-field values inside the Huygens surface are
close to zero. Also, it is noticed that the boundary area of connection between
the wire and the virtual surface shows very little disruption to the fields in the
scattered region, hence indicating that the proposed current crossing treatment
is implemented successfully. Figures 9.6(b) and 9.7(b) illustrate near-field
contours in the presence of scatterers. It can be noted that more back-scattered
field is observed inside the Huygens surface for the helix-dipole as compared
with the dipole. This is because of the closer proximity of the scatterer to the
crossing boundary in the first case.
Figure 9.3. Comparison between the computed near electric fields of a dipole in free
space (Example 9.1).
318 M. A. Mangoud, R. A. Abd-Alhameed and P. S. Excell
Figure 9.4. Input impedance of the source dipole with scatterer of Example 9.2, versus
number of iterations, for current-crossing hybrid treatment.
(a) (b)
Figure 9.5. Vertical slice cut for computed near electric field magnitude distribution
contours (in dB) for the Example 9.2: (a) without the scatterer; (b) with the scatterer.
Simulation of Antennas Coupled to Lossy Dielectric Volumes 319
Figure 9.6. Input impedance of the semi-helical dipole with scatterer of Example 9.3,
versus number of iterations, for hybrid current-crossing treatments.
(a) (b)
Figure 9.7. Vertical slice cut for computed near-field magnitude distribution contours
(in dB) for the Example 9.3: (a) without the scatterer; (b) with the scatterer.
320 M. A. Mangoud, R. A. Abd-Alhameed and P. S. Excell
(a)
(b)
Figure 9.9. The input resistance and reactance of the source dipole in Example 9.4, for
three different distances from the dielectric.
322 M. A. Mangoud, R. A. Abd-Alhameed and P. S. Excell
(a)
(b)
Figure 9.10. The far-field E component for Example 9.4, at two different cuts: (a)
vertical cut at = 90; (b) horizontal cut at = 90, for three different separation
distances from the dielectric. (— — 0.1, — — 0.15, — — 0.2).
Simulation of Antennas Coupled to Lossy Dielectric Volumes 323
Figure 9.11. Different views of duo-conical monopole antenna on perfect ground plane
(infinite).
324 M. A. Mangoud, R. A. Abd-Alhameed and P. S. Excell
Figure 9.13. Electric field magnitude (Ez) versus distance parallel to z-axis at 27 m
radial distance from the origin.
CONCLUSION
A novel hybrid method has been presented which permits conduction
currents to cross the boundary between different computational domains. The
Simulation of Antennas Coupled to Lossy Dielectric Volumes 325
objects can be physically connected or separated. The method uses the FDTD
method in the scatterer and part of the source region and uses the standard
frequency-domain MoM program in the remaining part of the source region to
give the complete solution of the fields in the two regions. Several validation
cases are examined and the results compared with available data. Also a short
wave broadcasting monopole over ground plane has been successfully
implemented using this modified technique. The hybrid method will have
great advantage over alternative methods in the following electromagnetics
applications such as cellular telephone dosimetry investigations close to real
biological tissues, complex satellite quadrifilar antennas, SAR reduction using
phased array antennas, base station safety assessment, subsurface radar, and in
general, interaction between complex source structures and inhomogeneous
dielectric materials.
REFERENCES
[1] K. Tap, P. Pathak and R. Burkholder, “Complex source beam - moment
method procedure for accelerating numerical integral equation solutions
of radiation and scattering problems,” IEEE Transactions on Antennas
and Propagation, vol. 62, no. 4, pp. 2052 - 2062, 2014.
[2] D. J. Ludick, R. Maaskant, D. B. Davidson, U. Jakobus, R. Mittra and D.
De Villiers, “Efficient analysis of large aperiodic antenna arrays using
the domain Green's function method,” IEEE Transactions on Antennas
and Propagation, vol. 62, no. 4, part 1, pp. 1579 - 1588, 2014.
[3] S. R. Zang and J. R. Bergmann, “Analysis of omnidirectional dual-
reflector antenna and feeding horn using method of moments,” IEEE
Transactions on Antennas and Propagation, vol. 62, no. 3, pp. 1534-
1538, 2014.
[4] X. -M. Pan and X. -Q. Sheng, “Improved algebraic preconditioning for
MoM solutions of large-scale electromagnetic problems,” IEEE
Antennas and Wireless Propagation Letters, vol. 13, pp. 106-109, 2014.
[5] M. A. Echeverri Bautista, M. A. Francavilla, F. Vipiana and G. Vecchi,
“A hierarchical fast solver for EFIE-MoM analysis of multiscale
structures at very low frequencies,” IEEE Transactions on Antennas and
Propagation, vol. 62, no. 3, pp. 1523-1528, 2014.
[6] M. Amir, S. Bedra, S. Benkouda and T. Fortaki, “Bacterial foraging
optimisation and method of moments for modelling and optimisation of
326 M. A. Mangoud, R. A. Abd-Alhameed and P. S. Excell
Raed A. Abd-Alhameed
Mobile and Satellite Communications Research Centre
Bradford University
Bradford, UK
Email: r.a.a.abd@bradford.ac.uk
Peter S. Excell
Institute for Arts, Science and Technology
Glyndwr University
Wrexham, UK
INDEX
B C
bandwidth, 7, 143, 144, 150, 153, 154, 155, capacitance, 46, 76, 272
156, 157, 160, 170, 171, 172, 176, 178, capacitor, 278
179, 180, 181, 182, 183, 185, 186, 187, Cartesian, 6, 19, 38, 61, 107, 209
189, 190, 191, 192, 195, 199, 203 cell, viii, 6, 8, 23, 25, 28, 33, 40, 41, 117,
basis function, vii, 5, 7, 14, 46, 51, 52, 54, 124, 138, 140, 213, 217, 220, 244, 252,
55, 62, 63, 67, 69, 77, 78, 81, 82, 87, 99, 255, 263, 266, 267, 268, 272, 275, 276,
100, 103, 149, 151, 167, 168, 169, 175, 277, 278, 279, 281, 283, 284, 285, 288,
193, 314 289, 290, 291, 299, 302, 304, 306, 313,
Berenger, 4, 7, 11, 12, 17, 38, 41, 43, 44, 314, 315, 316
135, 235, 242, 244, 247, 251, 252, 262, CFL criterion, 28
267, 304, 307 charge density, 48, 49
Bessel function, 236, 240 chiral, vii, 6, 7, 143, 145, 146, 147, 150,
bianisotropic, vii, 6, 7, 143, 144, 145, 154, 151, 156, 157, 158, 160, 162, 166, 167,
167, 171, 172, 173, 185, 186, 187, 188, 168, 169, 170, 171, 172, 175, 185, 186,
189, 192, 195, 196, 198, 200, 204, 205 187, 188, 191, 192, 196, 197, 200, 201,
biaxial, 148, 171, 185 204, 205
bi-isotropic, 168, 169, 186, 191, 193, 196 chirality, 7, 143, 145, 150, 154, 155, 162,
bioelectromagnetic, 8, 234, 235, 303 164, 165, 166, 167, 169, 170, 171, 175,
biological, viii, 3, 6, 8, 140, 208, 231, 234, 183, 184, 185, 187, 188, 191, 192, 193,
244, 252, 255, 263, 264, 265, 266, 267, 197
272, 275, 278, 282, 285, 286, 287, 288, chirostrip, 157, 186, 197, 199, 201
290, 291, 304, 308, 320, 325 chiro-waveguide, 158, 200
biological cell, viii, 6, 140, 252, 255, 265, chloride, 273
266, 272, 275, 278, 285, 286, 287, 288, circumferential, vii, 7, 46, 76, 80, 81, 82,
304, 308 84, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96,
biological tissue, 3, 208, 231, 266, 267, 282, 97, 98, 102
289, 291, 304, 320, 325 coarse grid, 4, 114, 115, 117, 128, 217
bi-static, 56, 58, 63 components, 3, 4, 17, 21, 23, 24, 25, 28, 29,
Bodies of Revolution, 66 30, 31, 32, 33, 34, 35, 36, 37, 38, 40, 41,
bone, 247 67, 69, 74, 76, 78, 82, 88, 89, 90, 91, 92,
boundary, viii, 2, 4, 5, 6, 7, 8, 9, 10, 11, 12, 93, 94, 95, 96, 97, 118, 119, 120, 145,
23, 32, 38, 42, 43, 46, 48, 50, 56, 67, 68, 146, 147, 153, 167, 170, 171, 173, 177,
115, 116, 117, 118, 119, 136, 137, 148, 180, 186, 188, 191, 192, 214, 215, 216,
173, 195, 210, 217, 218, 227, 234, 235, 218, 236, 239, 244, 268, 269, 270, 271,
236, 242, 244, 247, 251, 252, 261, 262, 273, 307, 310
264, 265, 266, 267, 268, 271, 275, 282, computational domain, 3, 4, 18, 115, 119,
288, 292, 304, 305, 307, 309, 311, 314, 124, 128, 129, 214, 242, 244, 245, 246,
316, 317, 323, 324, 327 247, 251, 255, 268, 277, 289, 290, 291,
Boundary Element Method, 67 297, 300, 309, 324
computational electromagnetic, 2, 17, 230,
309
conductance, 272, 273, 307
Index 333
conductivity, 19, 40, 41, 77, 118, 119, 126, 201, 202, 203, 208, 214, 217, 227, 230,
127, 208, 217, 235, 237, 243, 247, 251, 247, 260, 261, 263, 264, 267, 285, 286,
285, 313, 320 305, 307, 309, 310, 320, 321, 322, 323,
converge, 124, 316 325, 327
cosine, 268, 288 dipole, 68, 86, 87, 88, 90, 98, 99, 214, 217,
Courant stability, 4, 115, 117, 128, 130, 288 263, 315, 316, 317, 318, 319, 321
cross-section, 76, 78, 80, 88, 93, 98, 196, Dirac, 53, 54, 62
203 dispersion, 4, 12, 27, 186, 195, 205, 210,
cubical, 21, 216, 267, 282, 289, 290, 291, 226
292, 297, 299, 302 divergence, 48, 50, 82
curl, 2, 48, 211 duo-conical, 323
curl equations, 2, 211 dyadic, 148, 149, 151, 175, 196
current, vii, 5, 6, 7, 13, 15, 19, 46, 48, 50,
53, 56, 61, 63, 67, 68, 69, 74, 76, 77, 78,
82, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, E
97, 98, 101, 102, 103, 114, 116, 125,
EBG, 160, 185, 203
133, 140, 145, 149, 151, 168, 175, 176,
efficiency, viii, 2, 4, 6, 17, 28, 77, 97, 99,
193, 196, 260, 272, 273, 277, 279, 280,
144, 157, 171, 186, 199, 207, 222, 225,
306, 308, 311, 315, 316, 317, 318, 319,
226
323
eigenfunction, 208
current density, 5, 13, 19, 48, 53, 63, 67, 69,
electric charge density, 48
74, 77, 82, 90
electric conductivity, 19, 243
current distribution, vii, 5, 6, 7, 15, 46, 77,
electric conductor, 116, 234
78, 87, 90, 98, 103, 140, 151
electric current density, 19
curvilinear, 80, 82, 83, 87
electric field, 3, 7, 8, 9, 19, 21, 22, 23, 45,
cut-off frequency, 159, 160, 164, 166, 167
47, 48, 49, 50, 53, 57, 67, 78, 104, 117,
cycle, 117, 266
118, 119, 120, 121, 122, 124, 131, 132,
cylinder, 59, 60, 61, 63, 64, 65, 66, 98, 116,
133, 134, 140, 145, 149, 167, 168, 173,
130, 204
175, 194, 211, 214, 215, 222, 233, 236,
cylindrical, 46, 101, 198, 199, 267, 282,
247, 252, 255, 258, 259, 267, 268, 281,
291, 292, 300
282, 283, 289, 290, 291, 292, 293, 294,
cytoplasm, 255, 267, 281, 285, 286, 287,
295, 297, 298, 299, 301, 302, 308, 313,
289, 290, 293, 295, 298, 302
315, 316, 317, 318, 323, 326
Electric Field Integral Equation, 48
D electric flux density, 19
electric permittivity, 19, 38
DC, 29, 199, 278, 288, 308 electric potential scalar, 49
DC content, 29 electromagnetic, vii, viii, 1, 2, 3, 5, 6, 7, 9,
delta, 53, 54, 62, 87, 92 10, 11, 13, 14, 15, 18, 28, 38, 41, 43, 45,
delta-function, 87 46, 47, 101, 102, 103, 105, 113, 114,
deterministic, 281 115, 116, 117, 136, 137, 139, 141, 160,
dielectric, vii, viii, 5, 6, 8, 10, 12, 14, 15, 23, 168, 185, 186, 193, 197, 199, 201, 203,
42, 100, 114, 115, 116, 117, 119, 136, 208, 209, 210, 220, 227, 228, 229, 231,
137, 144, 151, 154, 158, 159, 163, 166, 233, 234, 235, 242, 259, 260, 262, 263,
171, 172, 179, 180, 185, 195, 198, 200,
334 Index
264, 265, 266, 267, 282, 306, 307, 309, 198, 207, 213, 222, 225, 226, 236, 241,
325, 327, 328 285, 286, 287, 314, 323, 325
electromagnetic band gap, 160 fundamental, 17, 41, 154, 167, 186
Electromagnetic Compatibility, 10, 11, 13,
43, 102, 104, 137, 139, 194, 200, 228,
229, 259, 261, 264, 307, 308, 326 G
EM wave, viii, 6, 130, 209, 243
Galerkin, vii, 6, 7, 9, 53, 69, 82, 104, 143,
energy, 114, 260, 263, 306
149, 158, 167, 170, 175, 176, 194, 210,
equivalence principle, 6, 33, 67, 68, 220,
259, 260, 326
309, 311
Gaussian pulse, 29, 32, 120
equivalent sources, viii, 311
geometry, 2, 4, 5, 46, 57, 76, 80, 87, 88, 89,
equivalent surface, viii, 7, 209, 211, 212,
90, 95, 96, 115, 172, 185, 203, 208, 210,
214, 235, 244, 247, 311, 314, 315, 316
215, 216, 218, 244, 268, 281, 291, 312,
equivalent surface currents, 211, 212, 214,
315, 320
311, 314
grading factor, 41, 243, 244, 247, 252, 255,
excitation, 28, 29, 31, 32, 41, 46, 51, 52, 56,
288
57, 62, 69, 70, 73, 74, 79, 87, 92, 114,
Green, 8, 49, 56, 71, 79, 99, 100, 145, 148,
144, 157, 159, 160, 166, 171, 196, 197,
149, 151, 158, 161, 167, 168, 172, 173,
209, 212, 214, 226, 276, 278, 279, 281,
175, 191, 193, 195, 196, 200, 213, 314,
283, 284, 288, 306, 314
325
extracellular, 285, 286, 287
grid, 4, 9, 10, 14, 20, 27, 28, 29, 31, 42, 114,
117, 119, 120, 121, 122, 124, 128, 132,
F 133, 136, 137, 209, 213, 217, 227, 242,
244, 245, 261, 267, 268, 276, 279, 281,
far field, viii, 7, 57, 214, 216, 320 288, 305, 310, 326, 327
fat, 247 grid dispersion, 27
FDTD lattice, 23, 24, 28, 29, 31, 38, 41, 220 gyro-electric, 145
FEM, 45, 47, 209, 230, 234, 260 gyrotropy, 154
ferrite, 157, 196, 197, 198, 199
fine grid, 4, 114, 115, 117, 119, 120, 128,
217 H
Finite-Difference Time-Domain, v, vii, 228
Haar, 74
Floquet, viii, 6, 8, 265, 266, 267, 268, 269,
half-power, 7, 143, 144, 150, 154, 170, 176
275, 282, 288, 292, 304, 307, 308
half-wavelength, 214, 217, 315, 316, 317
Fourier spectrum, 29
Hankel function, 236, 241
Fourier transform, 152, 168, 172, 176
harmonic, 197, 282
free space, 6, 8, 30, 49, 87, 90, 115, 116,
helical, 5, 6, 7, 15, 77, 80, 98, 104, 316, 319
118, 128, 147, 213, 214, 217, 220, 226,
helix, 46, 77, 79, 80, 81, 82, 84, 85, 86, 87,
233, 235, 237, 242, 243, 244, 246, 247,
94, 95, 96, 97, 99, 100, 108, 316, 317
248, 249, 250, 258, 314, 315, 316, 317
heterogeneous, 8, 10, 41, 136, 208, 227,
function, vii, viii, 5, 7, 8, 14, 20, 29, 32, 41,
234, 261, 304, 308, 309, 326
47, 49, 51, 52, 53, 55, 56, 71, 78, 79, 81,
high voltage, vii, 114, 116, 124, 125, 128,
83, 84, 87, 99, 100, 115, 126, 127, 145,
129, 130, 141
148, 151, 152, 158, 167, 169, 195, 196,
Index 335
Hodgkin, 8, 266, 272, 288, 289, 290, 296, 98, 99, 100, 101, 102, 103, 104, 115,
299, 302, 304, 306, 307 137, 143, 145, 149, 152, 158, 167, 172,
Hodgkin-Huxley, 8, 266, 288, 289, 290, 175, 176, 193, 194, 213, 229, 259, 310,
296, 299, 302, 304 314, 325, 326
homogeneous, 14, 28, 128, 130, 196, 207, integral equation, 5, 7, 8, 9, 15, 45, 47, 50,
209, 210, 216, 234, 235, 243 51, 53, 54, 56, 59, 61, 66, 67, 98, 99,
horizontal, 217, 220, 221, 222, 226, 322 101, 102, 103, 104, 115, 137, 143, 145,
human body, viii, 6, 7, 13, 140, 207, 208, 149, 152, 167, 172, 175, 176, 193, 194,
209, 221, 222, 226, 228, 230, 308, 310 229, 259, 310, 325, 326
human head, 216, 217 interaction, vii, viii, 2, 6, 7, 8, 14, 38, 44,
Huygens, 210, 214, 217, 220, 226, 244, 246, 79, 115, 117, 124, 130, 135, 140, 186,
247, 251, 252, 255, 316, 317, 323 207, 209, 210, 226, 229, 233, 234, 262,
Huygens surface, 210, 214, 217, 220, 226, 264, 265, 266, 291, 304, 307, 325, 327
244, 246, 247, 251, 252, 255, 316, 317, interference, 114, 116, 139
323 intermediate frequency, 289, 293, 297, 301
hybrid, viii, 1, 4, 6, 7, 8, 12, 13, 41, 46, 47, ions, 272, 273, 277, 306
99, 100, 158, 207, 208, 209, 210, 213, isotropic, 10, 42, 47, 48, 59, 137, 148, 150,
220, 226, 228, 229, 230, 264, 306, 309, 161, 163, 165, 167, 169, 171, 175, 176,
316, 318, 319, 323, 324, 327, 328 178, 179, 180, 186, 187, 191, 193, 196,
227, 305
iteration, 252
I iterative, 116, 204, 208
loop, 5, 15, 76, 86, 87, 90, 92, 93, 99, 102, 137, 145, 146, 173, 186, 194, 227, 259,
103, 263 262, 305, 310, 326
lossless, viii, 6, 8, 233, 235, 243, 244, 251, membrane, 8, 255, 263, 265, 266, 267, 272,
253, 254, 258, 268 273, 275, 277, 278, 279, 280, 281, 285,
lossy, vii, viii, 6, 8, 38, 116, 136, 233, 235, 286, 287, 289, 290, 292, 295, 299, 302,
242, 243, 244, 252, 255, 256, 257, 258, 304, 306, 307
267, 275, 278, 281, 283, 304, 309, 310, metal, 3, 115, 117, 126, 131, 310
320, 323 Method of Moment, v, vii, 45, 79, 230
lumped element, 266, 281, 304, 306 MIC, 156, 157
microdosimetry, 289
microstrip, vii, 3, 6, 7, 10, 99, 143, 144,
M 150, 151, 156, 157, 158, 160, 167, 170,
171, 172, 179, 185, 186, 187, 188, 189,
macroscopic, 145
192, 194, 195, 196, 197, 198, 199, 200,
magnetic, 3, 7, 8, 17, 18, 19, 21, 23, 24, 25,
201, 202, 203, 204, 205, 259, 326
28, 29, 30, 31, 32, 33, 36, 38, 41, 43, 45,
microwave, 2, 5, 11, 15, 43, 145, 156, 191,
47, 48, 49, 50, 53, 56, 67, 68, 69, 98,
194, 199, 204, 260, 263, 304, 306
101, 115, 117, 121, 122, 124, 128, 131,
Microwave Integrated Circuits, 156
132, 133, 134, 135, 140, 147, 148, 150,
mobile, viii, 6, 8, 78, 103, 140, 222, 228,
154, 157, 171, 172, 173, 185, 196, 199,
229, 244, 252, 258, 262, 264, 287, 293,
202, 203,209, 211, 212, 213, 220, 226,
306, 327
236, 239, 243, 244, 247, 252, 255, 278,
MoM, v, vii, viii, 1, 2, 3, 5, 6, 7, 9, 13, 41,
312, 313, 314
45, 47, 51, 77, 78, 100, 103, 104, 144,
magnetic current density, 19
193, 203, 207, 208, 209, 210, 211, 212,
magnetic field, 3, 8, 17, 18, 19, 23, 24, 25,
214, 217, 219, 226, 229, 230, 234, 259,
28, 29, 30, 33, 36, 41, 45, 48, 50, 56, 67,
310, 311, 313, 314, 315, 316, 323, 325,
115, 117, 121, 122, 124, 128, 131, 132,
327
133, 134, 135, 140, 147, 157, 172, 173,
monolayer, 160, 171, 177, 178, 179, 183,
209, 213, 220, 226, 236, 239, 244, 247,
187, 188, 192
252, 255, 278, 314
monopole, 323, 325
Magnetic Field Integral Equation, 50
multilayer, 149, 175, 195, 198
magnetic flux density, 19
magnetic permeability, 19, 38
magnetic potential vector, 49 N
magnetic resistivity, 19
magneto-chirality, 197 near field, 7, 216
magneto-electric, 145, 151, 166, 171, 172, Neumann function, 236
183, 186, 187, 188, 191, 192, 193 noise, 281, 284
magnitude, 30, 57, 83, 214, 222, 247, 285, non-uniform, 4, 12, 15, 46, 76, 77, 90, 97,
290, 291, 316, 318, 319, 323, 324 98, 116, 208
matrix, 5, 47, 51, 52, 56, 62, 63, 66, 69, 72, normal incidence, 244, 247, 252, 268, 271
73, 74, 79, 83, 98, 100, 106, 146, 149, numerical dispersion, 4, 12, 28, 210
150, 152, 158, 173, 176, 193, 240, 281 numerical grid, 27
Maxwell, v, 1, 2, 3, 9, 10, 13, 17, 18, 38, 41,
42, 45, 47, 49, 104, 113, 114, 115, 118,
Index 337
185, 187, 189, 193, 195, 198, 199, 200, 244, 246, 247, 251, 252, 266, 267, 278,
201, 202, 203, 236, 268, 310 279, 289, 290, 306, 309, 310, 315, 317
rectangular grid, 20 sine wave, 288
refinement, 128 singularity, 56, 84, 163
reflection, 6, 8, 12, 38, 43, 115, 117, 136, sinusoidal, 3, 11, 29, 30, 32, 47, 103, 120,
209, 235, 242, 243, 247, 251, 252, 258, 130, 193, 214, 227, 278, 288, 314
304 sodium, 273, 306, 307
relative permittivity, 126, 127, 144, 157, soil, 7, 126, 127, 128, 130, 136, 260
217, 237, 243, 320 solution, vii, 1, 4, 5, 6, 7, 10, 11, 15, 17, 18,
relaxation, 285, 286 41, 42, 44, 45, 46, 50, 52, 67, 76, 86, 88,
resistance, 77, 101, 199, 321 98, 101, 102, 103, 113, 114, 137, 144,
resonant, 7, 77, 99, 143, 144, 150, 151, 152, 147, 150, 163, 167, 176, 201, 214, 227,
153, 154, 155, 156, 157, 170, 171, 172, 234, 235, 236, 240, 243, 244, 247, 261,
176, 177, 178, 179, 180, 181, 182, 183, 262, 267, 276, 305, 307, 320, 325
184, 185, 186, 187, 188, 189, 190, 192, source, 7, 8, 9, 28, 29, 32, 40, 41, 42, 50, 51,
199, 202 58, 79, 80, 87, 92, 93, 103, 113, 114,
resonant frequency, 7, 143, 144, 150, 151, 115, 124, 125, 136, 160, 173, 209, 210,
152, 153, 154, 155, 156, 170, 171, 172, 211, 212, 213, 214, 217, 220, 226, 227,
176, 177, 178, 179, 180, 181, 182, 183, 228, 244, 261, 278, 288, 305, 311, 312,
184, 185, 186, 187, 188, 189, 190, 192, 313, 314, 316, 318, 321, 325, 326
199, 202 space resolution, 27, 28
resonator, vii, 7, 143, 144, 150, 155, 156, spatial resolution, 25, 117, 128, 129, 136,
157, 167, 176, 185, 186, 193, 194, 197, 208, 214
200 spectral domain, vii, 6, 7, 63, 143, 152, 167,
168, 172, 196
Spectral Domain Approach, 144
S sphere, 72, 106, 235, 236, 241, 243, 244,
247, 248, 249, 250, 251, 253, 254, 255,
sag angle, 323
256, 257, 258, 263
scalar, 19, 49, 52, 55, 62, 78, 191, 193, 208,
spherical, 107, 203, 235, 236, 237, 239, 244,
228, 260
255, 263, 266, 267, 281, 282, 288, 289,
scaling, 8, 69, 118, 119, 136, 233, 246, 258,
290, 291, 292, 296, 304
293, 304
stability, 4, 7, 12, 17, 27, 28, 41, 115, 117,
scatterer, 38, 41, 115, 210, 211, 214, 246,
124, 128, 130, 138, 216, 252, 288
251, 252, 311, 313, 315, 316, 317, 318,
standard deviation, 222
319, 320, 325
steady state, 25, 214, 244, 252, 272, 275,
scattering, 1, 2, 3, 5, 7, 8, 10, 11, 13, 14, 15,
276, 277, 278, 315, 320
18, 28, 30, 31, 32, 41, 43, 45, 46, 47, 99,
step, 4, 5, 9, 21, 25, 28, 29, 30, 32, 40, 41,
100, 101, 102, 103, 105, 116, 195, 196,
115, 117, 129, 130, 136, 217, 220, 226,
197, 201, 203, 204, 205, 208, 209, 210,
234, 246, 261, 266, 278, 281, 285, 286,
229, 235, 242, 246, 263, 264, 307, 308,
288, 304, 326
325, 327, 328
subgrid, 7, 119, 120, 121, 122, 124, 128,
sensor, 214, 216
129, 130, 131, 133, 134, 214, 222, 224
shock hazard, 132
simulation, 3, 8, 11, 12, 18, 27, 28, 38, 117,
124, 129, 130, 140, 145, 217, 226, 242,
Index 339
subgridding, vii, 4, 6, 7, 12, 13, 41, 115, TM, 29, 31, 40, 53, 54, 56, 57, 58, 59, 61,
116, 118, 119, 135, 138, 139, 208, 209, 62, 63, 64, 65, 66, 124, 147, 149, 158,
214, 217, 218, 219, 226 161, 164, 166, 171, 175, 186, 192, 197
substrate, vii, 6, 7, 15, 137, 143, 144, 150, TM grid, 29
151, 153, 154, 156, 157, 159, 160, 163, transmission lines, vii, 6, 7, 14, 113, 114,
165, 166, 167, 168, 170, 171, 172, 173, 116, 117, 124, 132, 135, 136, 139, 195,
175, 176, 185, 187, 188, 189, 190, 191, 202
192, 195, 196, 197, 199, 200, 201, 203, transverse, 31, 78, 146, 173, 191, 192, 193,
205 204, 236, 315
superconductor, 77 triangle, 5, 167
superstrate, 144, 171, 172, 173, 175, 176, triangular, 5, 53, 54, 55, 82, 185, 202
177, 178, 179, 180, 181, 182, 183, 184,
185, 186, 187, 188, 189, 190, 191, 192,
195, 197, 200, 202, 203 U
surface current, vii, 5, 6, 7, 46, 76, 77, 78,
82, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, underground, vii, 6, 7, 114, 116, 117, 135,
97, 98, 101, 103, 149, 169, 175, 176, 139
211, 212, 214, 311, 312, 313, 314, 323 uniaxial, 144, 145, 149, 150, 155, 170, 171,
surface current distribution, vii, 5, 6, 7, 46, 175, 176, 177, 178, 179, 180, 181, 182,
76, 77, 78, 87, 90, 98, 103 183, 193, 195, 196, 197, 200, 201
surface kernel, vii, 6, 7, 46 utility, vii, 114, 117, 124, 125, 135
T V
tangential components, 147, 236, 239, 270, vector, 3, 23, 25, 28, 30, 31, 32, 38, 40, 41,
271 45, 49, 50, 58, 61, 67, 78, 79, 105, 106,
tangential electric, 33, 48, 145, 167, 173, 118, 147, 172, 176, 201, 211, 213, 234,
268, 270, 271 236, 237, 238, 260, 312, 314
Taylor series, 21 vector dot product, 118
TE, 53, 54, 56, 57, 58, 59, 60, 61, 62, 63, vertical, 209, 220, 221, 222, 223, 226, 229,
64, 65, 66, 101, 147, 149, 157, 158, 161, 322
164, 166, 171, 175, 186, 192, 197 voltage, vii, 80, 87, 114, 116, 124, 125, 128,
three-phase, 114, 116, 117, 130 129, 130, 132, 136, 139, 141, 212, 213,
time constant, 274, 285 214, 272, 273, 276, 277, 278, 279, 280,
time domain, 5, 12, 13, 18, 42, 77, 114, 115, 281, 306, 313, 314
116, 140, 208, 228, 229, 262, 264, 305 volume density, 217
time resolution, 27, 28
time step, viii, 4, 6, 8, 21, 23, 25, 27, 28, 29, W
30, 32, 40, 113, 115, 117, 129, 130, 136,
217, 220, 226, 233, 234, 246, 258, 266, water, 127, 263, 267, 306, 310
278, 279, 281, 288, 304, 313 wave, viii, 2, 3, 6, 7, 11, 13, 15, 27, 28, 29,
tissue, 140, 208, 217, 226, 228, 231, 262, 30, 31, 32, 38, 43, 44, 45, 53, 105, 113,
264, 266, 267, 272, 282, 285, 288, 289, 116, 120, 130, 139, 143, 147, 157, 158,
291, 292, 304, 306, 320 159, 160, 167, 186, 196, 197, 198, 199,
200, 208, 209, 211, 217, 228, 233, 234,