Computation of electromagnetic fields in inhomogeneous media: Scattering and guiding properties

FIELDS IN INHOMOGENEOUS MEDIA:

SCATTERING AND GUIDING PROPERTIES

G.

MUR

COMPUTATION OF ELECTROMAGNETIC

FIELDS IN INHOMOGENEOUS MEDIA:

SCATTERING AND GUIDING PROPERTIES

f1|ill^»il!i!ra!!il!ll*>ll

iiilliiii

j S S i u ^

SCATTERING AND GUIDING PROPERTIES

PROEFSCHRIFT

TER VERKRIJGING VAN DE GRAAD VAN DOCTOR

IN DE TECHNISCHE WETENSCHAPPEN AAN DE

TECHNISCHE HOGESCHOOL DELFT, OP GEZAG VAN

DE RECTOR MAGNIFICUS PROF. IR. L HUISMAN,

VOOR EEN COMMISSIE AANGEWEZEN DOOR HET

COLLEGE VAN DEKANEN TE VERDEDIGEN OP

WOENSDAG 24 MEI 1978 TE 16.00 UUR

door

GERRIT MUR

elektrotechnisch ingenieur

geboren te Kockengen

DIT PROEFSCHRIFT IS GOEDGEKEURD DOOR DE PROMOTOR PROF. DR. IR. A.T. DE HOOP

Katja.

CONTENTS 7

GENERAL INTRODUCTION II

CHAPTER I

CALCULATION OF REFLECTION AND TRANSMISSION COEFFICIENTS IN ONE-DIMENSIONAL WAVE

PROPAGATION PROBLEMS * t 15 1. Introduction 15 2. Description of the configuration and

formulation of the problem 16 3. Reduction of the integral equation to

one of the Volterra type 20

4. The general case 21 5. Method of computation 24 6. Numerical results 25

References 28

CHAPTER II

COMPUTATIONAL ASPECTS OF THE SCATTERING OF ELECTROMAGNETIC WAVES BY A DIELECTRIC OBSTACLE

IN A WAVEGUIDE OF RECTANGULAR CROSS-SECTION ** t 29 1. Introduction and formulation of the problem 29 2. Fredholm integral equation (F.I.E.) method 32 3. Volterra integral equation (V.I.E.) method 35

4. Differential equation (D.E.) method 39 5. Numerical results and conclusions 43

CHAPTER III

A DIFFERENTIAL-EQUATION METHOD FOR THE

COMPUTATION OF THE ELECTROMAGNETIC SCATTERING

BY ANINHOMOGENEITY IN A CYLINDRICAL WAVEGUIDE *** 49

1. Introduction 49 2. Description of the configuration 50

3. Derivation of the "Generalized Telegraphist's

Equations" 51 4. Transformation of the system of differential

equations into a stable one 56 5. Discussions of the numerical techniques 57

6. Numerical results 59

7. Conclusion 65 Appendix A: The method of invariant imbedding 65

References 67

CHAPTER IV

COMPUTATION OF ELECTROMAGNETIC FIELDS SCATTERED BY A CYLINDRICAL INHOMOGENEITY IN A HOMOGENEOUS

MEDIUM OF INFINITE EXTENT 69 1. Introduction 70 2. Description of the configuration 71

3. The electromagnetic-field equations and the

related system of differential equations 72 4. Transformation of the system of differential

equations into one that is stable 79 5. Discussion of numerical techniques 83

6. Numerical results 85

7. Conclusion 92 References 95

CONTENTS

CHAPTER V

COMPUTATION OF ELECTROMAGNETIC FIELDS GUIDED BY A CYLINDRICAL INHOMOGENEITY IN A HOMOGENEOUS

MEDIUM OF INFINITE EXTENT 97 1. Introduction 98 2. Description of the configuration 99

3. The electromagnetic field, equations and

the related system of differential equations 100 4. Transformation of the system of differential

equations into one that is stable 105 5. Discussion of the numerical techniques 109

6. Numerical results 112

7. Conclusion 115 References 1 16

SAMENVATTING

* This part has been published by G. Mur and A.J.A. Nicia in: Journal of Applied Physics, Vol. 47, 1976, pp. 5218-5221.

** This part has been published by G. Mur, D. Quak and G.J. van Dijk as a conference paper for the International Conference on Numerical Methods in Electrical and Magnetic Field Problems,

S. Margherita Ligure (Italy), June 1-4, 1976.

*** This part has been published by G. Mur in:Journal of Engineering Mathematics, Vol. 12, 1978, pp. 1 5 7 - 175.

t The author of this thesis wishes to thank the coauthors of Chapters I and II for their kind permission to include these papers in the present thesis.

A.J.A. Nicia carried out the numerical computations for the one-dimensional wave propagation problem in Chapter I when being a research student. D. Quak formulated and programmed the Fredholm integral-equation method that is reported in Chapter II, he is a staff member of the Laboratory of Electromagnetic Research, Department of Electrical Engineering of the Delft University of Technology. G.J. van Dijk carried out the numerical computations pertaining to the Volterra integral-equation method and the differential-equation method when being a research student.

GENERAL INTRODUCTION 11

GENERAL INTRODUCTION

In telecommunication engineering, microwave power engineering, plasma physics as well as in other areas of application for electromag-netic fields, the computation of time-harmonic electromagelectromag-netic fields in inhomogeneous media is a subject of considerable importance. In some applications, the inhomogeneity can be considered as an obstacle that either scatters or absorbs the incident electromagnetic-field energy or both. In other applications, the inhomogeneity serves to guide elec-tromagnetic energy from one point in the configuration to another. The methods used for the computation of the relevant fields are either of an approximate or of a rigorous nature. In the present thesis we primarily investigate rigorous methods for computing time-harmonic electromagnetic fields in inhomogeneous media. Of the rigorous methods employed for this purpose, integral-equation methods, using integral equations of the Fredholm type, are a standard tool. By using the method of moments, these integral equations are replaced by a system of algebraic equations, that is then solved numerically.

In the present thesis, two alternatives to this method are described and investigated. The first alternative also starts with a formulation of the field problem in terms of a Fredholm integral equation. Rather than being solved, this integral equation is transformed into a related system of coupled Volterra integral equations of the second kind. The system is then solved by a direct numerical quadrature. The major part of the thesis deals with a second alternative which avoids the use of integral equations and starts from the electromagnetic-field equations transforming them into a system of coupled first-order ordinary dif-ferential equations, by introducing a certain direction of propagation of the waves and an appropriate transverse modal expansion of the fields. The latter system can then be solved using standard techniques. The system of equations thus obtained, however, turns out to be numerically unstable and it is therefore transformed into a stable system through a specific transformation scheme.

Having described, in general terms, the subject of the present thesis, we will now give a resume of its five chapters. The chapters are each self-contained and they are arranged in the chronological order of their composition.

In Chapter I, we investigate one-dimensional wave reflection and transmission in an inhomogeneous medium. The wave propagation problem is first formulated, rigorously, in terms of an integrodifferential equation of the Fredholm type. Since the numerical solution of this equation is complicated and time consuming, it is transformed into an integral equation of the Volterra type. The solution of the latter can be carried out easily and rapidly. Reflection coefficients are presented for a layered troposphere and for a metal film that is deposited on a glass substrate.

In Chapter II, three different methods for the computation of the scattering properties of a two-dimensional cylindrical obstacle in a waveguide of rectangular cross-section are compared. These methods are (1) an integral-equation method of the Fredholm type, (2) an inte-gral-equation method of the Volterra type and (3) a method using a system of coupled, first-order ordinary differential equations. When applied to the same obstacle, these methods yield, within a high degree of accuracy, the same result. From the point of view of the amount of computation time required, the Volterra integral-equation method and the differential-equation method proved to be superior to the Fredholm integral-equation method. As compared with the Volterra integral-equation method, we note that the advantages of the differential-equation method are that it is less complicated and somewhat faster and that it can easily be generalized to more complicated configurations. A disadvantage is, however, its instability and, consequently, it cannot be used to compute the scattering properties of relatively "large-size" obstacles.

In Chapter III, the differential-equation method is applied to the computation of the scattering properties of an arbitrary inhomogeneity in a cylindrical waveguide. The system of differential equations to be used for the computation of the fields in the region where the inhomo-geneity is present, is the well-known system of "Generalized Telegra-phist's Equations". In this chapter, a method is introduced for coping

GENERAL INTRODUCTION 13

with the numerical instabilities mentioned earlier. Reflection and transmission factors are presented for a number of obstacles of varying complexity, present in a waveguide of rectangular cross-section. We mention the so-called E-plane applicator, used to dry a wet leather sheat with the aid of microwave power, reflection by a radially inho-mogeneous cylindrical plasma column as used in plasma diagnostics and reflection by a leaning pyramidal waveguide termination in a WR 90 waveguide. The latter problem is a three-dimenisonal one.

In Chapter IV, the differential-equation method is applied to the computation of electromagnetic fields scattered by a cylindrical inhomogeneity in a homogeneous medium of infinite extent. Again a

method is given to cope with numerical instabilities. Far-field scattering patterns are presented for an inhomogeneous dielectric cylinder and

for a number of homogeneous cylinders of various cross-sections. The

incident fields are either E- or ö-polarized, or have an arbitrary

angle of indicence.

Finally, in Chapter V, the differential-equation method is applied to the computation of surface-wave modes guided by a cylindrical inhomo-geneity in a homogenous medium of infinite extent. The analysis of this configuration runs, for the major part, parallel to the analysis presen-ted in Chapter IV, the main difference being that in Chapter V we are interested in the eigenvalues, and their corresponding solutions, which entails an iterative procedure. Propagation coefficients are presented for surface-wave modes guided by dielectric waveguides having a rectangular and an equilaterally triangular cross-section.

CALCULATION OF REFLECTION AND TRANSMISSION COEFFICIENTS IN ONE-DIMENSIONAL WAVE PROPAGATION PROBLEMS

G. MUR and A.J.A. NICIA

Department of Electrical Engineering, Laboratory of Electromagnetic Research,

Delft University of Technology, P.O. Box 5031,2600 GA Delft,

The Netherlands.

One-dimensional electromagnetic wave propagation problems can be formulated exactly in terms of an integrodifferential equation of the Fredholm type. The numerical solution of this equation is often complicated and time consuming. In this paper a method is described of transforming this equation into an integral equation of the

Volterra type. The solution of the latter equation can be carried out easily and within a short time.

1. INTRODUCTION

In several domains of physics the computation of reflection and transmission coefficients associated with one-dimensional wave propa-gation problems in inhomogeneous media is a subject of interest. The methods that are known often use fairly complicated mathematics, some of them yielding only first-order approximations (in some sense) [1,2,3], others using asymptotic methods such as the WKB approximation 14]. The present approach to the problem is first formulated, exactly, in terms of a certain Volterra integral equation of the second kind. Then this integral equation is solved numerically by an efficient procedure that is rather easily carried out. A somewhat different

16 ONE-DIMENSIONAL WAVE PROPAGATION

Volterra integral equation approach has been presented by Wang [5]; his integral equation is of such a structure that it can only be solved by a time consuming iterative procedure. We believe that our method is to be preferred.

2. DESCRIPTION OF THE CONFIGURATION AND FORMULATION OF THE PROBLEM The theory will be presented for the case of a plane electro-magnetic wave normally incident upon an inhomogeneous planar layer of finite thickness (Fig. 1). A Cartesian coordinate system is intro-duced, with the X and y axes parallel to the layer and the z axis

normal to it. We subdivide the configuration into three domains D.,

with i = 1,2,3, each of them being infinite in extent in the x and y direction and such that

öj: - » < 2 < 0, with e(2) = e^ , vi(3) = p^;

D^: 0 < z < I, with c(z) = c^(,z) , p(3) = U2(3); D : I < z < «>, with e(3) = e , u(3) = y,.

The layer under consideration is uniform in the transverse direction, but inhomogeneous in the 2 direction. The inhomogeneity is assumed to be present both in the permittivity e(2) and in the permeability viz); these quantities can be either real or complex. Considering a

linearly polarized electric field, we can take E = E i- as the

funda-— x-x

mental unkown quantity, which entails no restriction as regards the generality of the solution. Assuming a complex time dependence through

the factor exp(—toot), where i is the imaginary unit, o) is the circular

frequency, and t denotes the time, the source-free electromagnetic

field equations for the problem reduce to

3 fl - iiLEE = 0 , ( l a )

z y X ' d E - iüjyff = 0 . ( l b )

First we shall describe the method for the case where the medium in

0_ has the same electromagnetic properties as the medium in D (i.e.,

E- = e, , VT = p,). We return to the general case in Sec. 4. Introducing the wave numbers

k^ = co(e^y^)^ Im(fe^) > 0 (i = 1,2,3), (2)

we can write the electric field of the incident plane wave, which is defined everywhere in space, as

EI = exp(tfe,2). (3)

The incident wave satisfies the free-space electromagnetic field equations 3 H'^ - iu^e.E'^ = 0, (4a) z y \ X ' ^ ' 3 È"^ - imv.H'^ = 0. (4b) z X I y ' s s

The scattered field {E ,H } is now introduced as

x' y

{E^,lf-] = {E - E'^, H - H^] (5)

x' y' X x' y y ^ ' everywhere -.n space. Subtracting (4) from (1) and rearranging the result we arrive at

\^y - ^"^^l^x = ^'^(^ - ^O^x' (6^)

From (6) we observe that the scattered field is generated by two source distributions, viz., a source distribution of the magnetic type with volume density ia)()j - \i )H and a source distribution of the electric type with volume density ztü(e - £.)£• . These source distributions are

18 ONE-DIMENSIONAL WAVE PROPAGATION

located in D.. We shall determine the contribution from these two source

distributions to the scattered field separately, and afterwards super-impose the results,

2

X,.

y

z = 0 Z--1

Fig. 1. One-dimensional inhomogeneous configuration.

A. Inhomogeneity in the p e r m i t t i v i t y .

I yields

In this case we have e * e,, u = u,. Elimination of H from (6)

1 ' 1 y

(^^^X = --'^l(--^l)V

If the right-hand side of (7) were known, the solution of (7) would,

using the Green's function technique, be obtained as

•I 2

(7)

El(z ) = / ^ a ) % [e (2) - £ ]£ (s)G(2 2)d2,

X f

_{r x'}

in which the one-dimensional Green's function is given by

Giz^z) = (•i/2fej) exp(i;fe||2^ - 3|), (9)

B. Inhomogeneity in the permeability. I

yields

In this case we have e = e., y * y.. Elimination of E from (6)

Again using the Green's function technique we obtain

ff^(2p) = I Q cü^e,[y2(2) - V^'\B^{z)G{Zp,z)dz. (11)

In this case, using (lb) and (6a), the representation for the scattered electric field is obtained as

rl y (3) - u

^'(V = - ' z j 0 p^(2) i\E^iz))Giz^,z)dz. (12)

Adding the two expressions (8) and (12), and using (5), we obtain for the total electric field

E^iz^) = EliZp) + J Q Oi M^U^iz) - E,]i?^(2)G(2p,2)d2

v^{z) - V (13)

% y I \ ( 2 ) ' 0/,(3))G(3p.2)d2.

As an integral representation for the total field (13) holds everywhere in space; when 2p e D„, it is a Fredholm integrodifferential equation of the second kind. Next we introduce the electric-field reflection factor R and the electric-field transmission factor T by writing

Ê'^(2p) = exp(£fc|2^) + R exp(- ik^z^) for z^ e D^, (14)

and £'^(3p) = T exp(i/Cj2p) f o r z^ e Dy (15) From (13) we o b t a i n , by t a k i n g z e D , u s i n g ( 9 ) , R = ^ | - JQ Ü) y j [ e 2 ( 2 ) " e^ ]Ê'^(2)exp(ffe|2)<i2 ; y , ( s ) - y. - 'i !Q — ^^z^x^""^^ e x p ( i f e , 2 ) d 2 . y2(2) ( 1 6 )

20 ONE-DIMENSIONAL WAVE PROPAGATION

In the same way, by taking 2- e D- and using (9), we have

i tl 2

T = 1 + 2 ^ JQ 0) y|[e2(2) - e.j]£^(2) exp(- ik^z^dz 0 - y,

,(2) 2 X

r^ y,(2) - y, ^'^^

+ W Q ^ ,,,—^ O^ff^Cs)) exp(- ik,z)d2.

'2

3. REDUCTION OF THE INTEGRAL EQUATION TO ONE OF THE VOLTERRA TYPE Upon adding and subtracting terms in the right-hand side of (13), we rewrite this equation, using (2), (9), and (17), as

E^iZp) = T exp(ife,2p . 2 ^ \l^ .\,l^^iz) - e,]^^(2)

X {expLife (2 - 2 p ) ] - expCife.(2 - z)^]dz tl Vj{z) - y , - i / g „ r„s OjJz)^ {exp[ife,(2 - 2„)] (18) Zp U^iz) Z X ' ^ 1 + exp[i?C.(2_ - z)]}d2, for Zp e D .

The solution of (18) is constructed with the aid of the solution of the following Volterra integrodifferential equation [8]:

i>(Zp) = expiik^Zp) + 2F" ^2 '^^^^^^2'-^^ ~ ei]'J'(2)

X {exp[ife|(2 - 2 ^ ] - explïk^iZp - 2 ) ] } d s

fl V^^z) - y , (19) " V ^ P >^2<^> (\HB)){explik^iz - Zp)3

+ explïk^iZp - z)l}dz, for Zp e D^.

Let us assume that (19) has been solved; then

^^(3p) = TH^p) (20)

is obtained upon substituting (20) into (17). This results into the following expression:

= M - 2 ^ IQ '^ ^\^^2 ~ '^i^'''^^) exp(- ik^z)dz

(1^2 " ^"1 \-l

ij Q-^ Og^Cs)) exp(- ik^z)dz) .

(21)

Having solved (19), the integrals in (21) have a known value and hence the value of T is obtained.

4. THE GENERAL CASE

In the case where the electromagnetic properties of the medium in

D differ from those in D , slight modifications have to be made to the

method of solution. The incident field E is now chosen to be the field

X

that would exist if the properties of the medium in ö_ were the same as those in D . Hence, we have

E (2) = expiik z) + p explik ill -s)], for z e D u Ö

T explik I + -ik iz - 1)1, for z e D ,

(22)

where

P = • i r (23)

(£,/y,)^ + ic^/vj)'

is the Fresnel reflection factor of a plane boundary and 2(E,/y,)^

(e,/y,)^ +

ieJu.)'

(24)

is the corresponding Fresnel transmission factor. The influence of the domain 0 , where the actual configuration differs from the configuration

in which the incident field has been defined, is again described with the aid of a source distribution that generates the scattered field.

22 ONE-DIMENSIONAL WAVE PROPAGATION

Again a representation for this scattered field is obtained, using a Green's function technique. The relevant Green's function is now defined in the same configuration as the incident field, and is obtained as

GiZp^z) = (i/2fej){exp(tfcj]2p - z\) + p expHk^ill - z - Zp)l}

(25) for z e D u D , z e Dy,

= iil2k,)i expiikAl - 2) + ikAZp - Z)]

for Zp c D^, z e Dj.

Proceeding in the same way as in Sec. 2 , we now obtain an integral

't •

equation of the same form as (13), E and G being now defined as in

(22) and (25), respectively. The reflection factor R and the transmission

factor T for the inhomogeneous configuration are introduced through

^X^^P^ = exp(-z;fej2p) + R exp(- ik^Zp) for Zp e D^, (26)

E iZp) = T expCife Z + ik^izp- l)~\ for Zp e Dy ill)

Taking 2- e D. we now obtain from (13), using (22) and (25),

i I 2

R = p expilik^l) + 2 ^ JQ (O yj[e2(2) - z^lE^iz)

X {exp(t?c 2) + p expCife (2Z - z)~\}dz rl v^iz) - p, ^

X {expiik z) + p expiik ill - z)l}dz.

Taking 2- e C,, we arrive at

T = T(I + ^ JQ co^y,[e2(2) - ^^lE^iz) exp(- ik^z)dz

+ J^ t^Z^^^ " ^r (3/^(2)) exp(- i/c,2)d2). y O y (2)

(28)

Again, the integrodifferential equation can be reduced to one of the Volterra type. Upon adding and subtracting terms in the right-hand side of (13) and using (22), (25), and (29), we now arrive at

\

^w^ Lp'^^'^i^^i^^^ - ^i^^x(^)

-^Jzp u,iz) (^V-)>

The solution of (30) is constructed with the aid of the solution of the following integrodifferential equation of the second kind:

\piZp) = expiik^Zp) + p exp[ifej(2i. - Zp)l

^W-jlp '^^^\^^2^^'>

- ^i]^(«)

-

^ P ~ ^

Let u s assume t h a t (31) h a s been s o l v e d ; then

E_^iZp) = iTh)i,iZp) (32)

24 ONE-DIMENSIONAL WAVE PROPAGATION

(29). This results in the following expression:

r l 1

^''^ ' (' ~ W' ^0 " f^jtez^^) " ^i^'''^^) exp(- ik^z)dz

fl y2(2) - y, N_,

- V o y,(2) (V(^)> ^^P(- ^^«^^^) •

(33)

5. METHOD OF COMPUTATION

We shall briefly outline the numerical method employed for solving either (13) or (31). First, the equation to be solved is written in a

form in which the dependence on z„ is removed from the integrands.

The resulting equation contains four integrals. The evaluation of each

of these integrals is carried out stepwise, with N steps of step size

h = l/N, from 2 = Z to 3 = 0. The integration over the wth step

in = 1 N) is carried out as follows: First we notice that the

value of the integrals and of i|'(2p) and 3 iiis^) is known at the upper

P z p

boundary of the interval Iz^ = I - in - 1)«]. We now make an

inter-mediate approximation for iiil - nh) and 3 iiil - nh) at the lower

Zp

boundary of the nth interval, for example, by using a linear extra-polation formula. Using this approximation, we compute the integrals over the interval under consideration with the aid of the Euler-McLaurin summation formula [61. The value of the integrals being known, we

determine a new approximation for ipil - nh) and 3 i|)(l - nh) , and update the numerical approximation for these functions. This iterative process is continued until a sufficiently accurate result has been obtained. In practice, only one or two steps in the iteration are required. Now the process of integration requires the evaluation of 3 ip and it is

^P

noted, that this is done without applying numerical differentiation. Consequently, in each step of the integration process, the error in the

4 integrals will be of the order Oih ) .

In the analysis of our problem we have chosen the electric field £ = £• i as the fundamental unknown quantity. As an alternative, the

magnetic field H = H i could also have been chosen as such. The

y~y

given for the electric field. The fact that now two equivalent, but different, formulations of our problem are available can be used as a check on computer programs. An additional check on the numerical results is obtained using the law of energy conservation.

6. NUMERICAL RESULTS

Numerical results are presented for two different configurations of physical interest. The first configuration to be considered is that of a slowly varying dielectric medium which serves as a model for a layered troposphere whose refractive index profile follows from the

-70-, -80- -90- ?.100-u V "

-120-JiPhase dtflerence" solution (n si)

t,(z)=(j(1.00001 5in(rm2//l)

a 2-!rll\ ,0

Fig. 2. Reflection by a tropospheric inhomogeneity of the type e = e [ 1 -i- 0.0001 sin(mr2/Z)].

situation in Fig. 1, where the value of the permittivity is taken as one also considered by Tang [3]. The value of the permittivity in region D„ is taken as

£2 = EQCI + 0.0001 sin(mi2/Z)] n = 1,2,3,4, (34)

and y_ = y_; E_ and \i^ denote the permittivity and the permeability

•^0' 0 ^0

of vacuum, respectively. In the regions D and 0 the medium is assumed

to be vacuum, whence E = e , E_ = E „ , and y

26 ONE-DIMENSIONAL WAVE PROPAGATION

I 12

the reflectance |P| of the inhomogeneity m the troposphere is plotted as a function of 1T\1/X^ for n = 1,2,3,4; A„ denotes the wavelength of the incident wave in free space. In Fig. 2 we have also plotted the results for n = 1 that have been obtained by Tang using the "phase difference method" [3]. A considerable difference between our and

his results is noted, especially for small values of I. Since the

relative permittivity is D. in very close to unity, the Rayleigh-Gans or Born approximation [7] will lead to relatively good results provided

that I is not too large. Hence the correct solution should be close to

the one obtained with the Born approximation. Our results show this behavior whereas those of Ref. [3] do not. Further, in the limit of zero width of D„ the reflectance should go to zero. Our results are in accordance with this property, some of those in Ref. [3] are not. Finally, we have successfully applied the methods for checking our computer programs that are given in Sec. 5.

The second configuration to be considered consists of a metal film of thickness J>i_ (A_ being the wavelength of the incident field) which is deposited on an infinitely thick glass substrate. In Fig. 3 we have inserted a picture of the configuration together with a

description of the properties of the media in D., D., and P-, and a

plot of the assumed density Niz) of the free electrons whose charge

is q and whose mass is m. The electron density Niz) continuously fits the neighboring values with zero slope through the use of a third-degree polynomial. For e we have

E2 = Eg -H to/o), (35)

with

2

a = EQCO^ /(V^ - tio), (36)

where v denotes the collision frequency of the electrons and oj denotes their plasma frequency

2.,.^., ___,^

j)p = iq Niz)lz^my (37)

J25%,„,

-1 1 1 1 1 1—

2 ^ ,^iu. 6

Fig. 3. Reflection of a plane wave incident upon a glass substrate with metal film. For 0 < 2 < Z/4

we have chosen Niz) = [48(2/^)^ - \18iz/l)^W^.

In Fig. 3, the modulus of the reflection coefficient R of the

config-uration has been plotted as a function of üj/ojp, for several values of (.p/v^.

ACKNOWLEDGMENT

The authors wish to thank Professor A.T. de Hoop for his suggestions and remarks concerning the research presented in this paper.

28 ONE-DIMENSIONAL WAVE PROPAGATION

REFERENCES

[1] A. Banos Jr., J. Math. Phys. j_4, 963 (1973). [2] D.D. Evans, Radio Sci. 8, 1083 (1973). [3] C.C.H. Tang, J. Appl. Phys. 45, 1115 (1974).

[4] K.C. Yeh and C.H. Liu, "Theory of Ionospheric Waves", (Academic,

New York, 1972), Chap. 5, p. 252.

[5] T.N.C. Wang, IEEE Trans. Antennas Propag. AP-13, 986, (1965). [6] E. Isaacson and H.B. Keller, "Analysis of Numerical Methods",

(Wiley, New York, 1966), Chap. 5, p. 340.

[7] M. Born and E. Wolf, "Principles of Optics" (Pergamon, London,

1959), Chap. 8. p. 452.

[8] G.F. Drukarev,"Theory of Electron-Atom Collisions", (Academic, New

COMPUTATIONAL ASPECTS OF THE SCATTERING OF ELECTROMAGNETIC WAVES BY A DIELECTRIC OBSTACLE IN A WAVEGUIDE OF RECTANGULAR CROSS-SECTION

by G. MUR, D. QUAK and G.J. van ÜIJK Department of Electrical Engineering,

Delft University of Technology, P.O. Box 5031,2600 GA Delft

The Netherlands.

SUMMARY

In this paper a comparison is made between three methods for computing the scattering properties of a dielectric obstacle in a waveguide. The first method employs the two-dimensional Green's theorem and leads to two coupled integral equations of the Fredholm type. For the second method, the obstacle is replaced by a polarization current distribution; it leads to a system of coupled integral equations of the Volterra type. In the third method, which is new and appears to be the most promising one, the electromagnetic field equations are transformed into a system of coupled first-order ordinary differential equations.

1. INTRODUCTION AND FORMULATION OF THE PROBLEM

In this paper a comparison is made between three different methods for computing numerically the scattering properties of a cylindrical, dielectric obstacle in a waveguide of rectangular cross-section. After formulating the problem we recall the formulas from which we start the discussion about the methods of solution. A Cartesian coordinate system

30

obstacle is uniform; the 2-axis is chosen parallel to the axis of the waveguide. The waveguide walls are assumed to be perfectly conducting. The medium in the guide is taken to be homogeneous, linear, isotropic and lossless with permittivity e . The medium of the obsta'sile is taken to be linear, isotropic and uniform in the i/-direction. It may be inhomogeneous (in the x,2--plane) as well as lossy; its complex permit-tivity is Zj = e2(^'^^* '^'^^ permeability is taken to be y.

(= 4 ÏÏ X 10 H/m) everywhere. All fields are taken to vary sinusoidally in time with angular frequency u. The complex time factor exp(j(i)t;) is omitted throughout. A longitudinal section of the configuration is shown

1 1 Region 1 I Region 2 i Region 3

-l^

Z Z=Z, Z=Zj

Fig.1. Longitudinal section of the configuration,

in Fig.1. Depending on the method of solution or the fields to be dis-cussed, we divide, in two ways, the longitudinal section into different regions. The first method distinguishes the Region S inside the obstacle

and the Region S' outside it. In the second one we divide the waveguide

into Region 1 , which is uniform and occupies - ' ^ < z < z . Region 2 which contains the obstacle and occupies 2 < 2 < 2„ and Region 3 which again is uniform and occupies z^ < z < «>. The choice of z and z. is such that 2, - 2, is minimized. For simplicity, we take as incident field the LSM --dominant mode (Collin, 1960a). In this case, no coupling arises

1 ,u

between LSE- and LSM-fields, and E = E (x,2)i . The source-free

electro-- y electro--y

magnetic field equations reduce to

^2^X -

\ h - ^-^^^'^^%

= 0.

'z^y - J'^^O^x =

'^' (•)

where E ( X , 2 ) denotes the local value of the permittivity. Elimination

of H and H from (1) leads to the equation for the only nonvanishing

X 2

component of the electric field

(3^ -H 3^ + k^ix,z))E ix,z) = 0 » (2)

X 3 y

with

k = uiie.]!-.)^ when (x,2) e S',

kix,z) = (3)

k ix,z) = oiie ix,z)v ) ^ when (x.s) e S,

In the uniform regions, the field can be expressed in terms of LSM „--modes whose modal functions follow from

+ — * —1

(() (x,s) = sin(mTrx/a)M iz), (f (x,2) = sin(mTrx/a)M (2), (4)

m

m

m = 1, 2, ..., with

u iz) = exp(- T z), u (2) = exp(r 2 ) , (5)

m '^ m ' m ^ m ^ ^ where

r^ = iim-n/a)^ - fej)^ Re(r^) > 0, Im(r^) > 0, (6)

m = 1, 2 As incident field E we take the LSM -mode travelling

y ' »o in the positive 2-direction

i

= *;• (7)

The scattered field E is now defined in the entire waveguide as

E® = £ - E^. (8)

32 COMPUTATIONAL A S P E C T S OF SCATTERING THEORY

In Region 1 the scattered field can be w r i t t e n as

E^ = 1° R f , (9)

y m=\ m ^m '

w h e r e R denotes the reflection factor for the m - t h m o d e . Using ( 7 ) , (8) and (9) the total field in Region 1 is obtained as

E = (f)! -t- s " , /? (fi~. (10)

y ^l m=\ m m

In Region 3 the transmitted or total field can be written as

ff = Z°° , T (t>* , (11) y m=l m m ^

w h e r e T denotes the transmission factor for the m-th m o d e . In the next m

three sections w e shall successively discuss three different m e t h o d s for computing the reflection and transmission factors. In Section 5 w e shall present numerical results that have b e e n obtained using these three methods and m a k e a comparison between the results in the three cases and the computational effort involved.

2. F R E D H O L M INTEGRAL EQUATION (F.I.E.) METHOD

T h e first method of solution employs as the m a i n tool the two-dimen-sional Green's theorem. Since in this method Green's functions are involved, which can only be determined easily in homogeneous r e g i o n s , the applicability of this method is restricted to the case E_, = constant (i.e. a homogeneous o b s t a c l e ) . W r i t i n g E^^\x,z) w h e n (x,2) e S' , E ix,z) = (12) y (2) E^ ix,z) when (x,2) e S, (2) results into

(3^ + ^l + fe^)£'^'^(x,2) = 0 when (x,2) e S', (13)

? ? ? (2)

'•x ^ z '^ kpE^ ' ix,z) = 0 when (x,2) e S. (14)

Applying Green's theorem to S', the scattered field E and the Green's

function (Collin, 1960b)

oo - 1 1 1

G = E _ (ar ) sin(mTrx/a) sin(wTTXp/a) exp(- V \z - 2p|), (15)

which satisfies the boundary conditions at the waveguide walls, we obtain

^^(G,3 E^ - £^3 G,)ds = {1 , ^ ,0}£'^(x„j2j,)

^C \ n y y n ^ '^' y^ P' P

when iXp,Zp) e {S',C,S}. (16)

The principal-value integral in (16) applies to the case (Xp,2p) e C, Applying Green's theorem to S, the Green's function G and the incident field, we obtain

^C^^l V ^ "

''I'n^^^^' =-{OAA)EliXp,Zp)

when iXp,Zp) £ {S',C,S}. (17)

Upon adding (16) and (17) we obtain, using (8) and (12),

4^(G,3 E^'^ - ff^'^3 G,)ds

+ E^ix^,z^) =

^C \ n y y n \ y P' P

= {l,i,0}Ê'^'^(Xp,2p) when iXp,Zp) e {S',C,S}.

Applying Green's theorem to S, the Green's function

.2 . ,„ . .2J

(18)

34 COMPUTATIONAL ASPECTS OF SCATTERING THEORY

where Y. denotes the Bessel function of the second kind and of order

(2)

zero, and the total field E , we obtain

y

h^^l'rf^r - 'y'^'n^2^^'

=

{OA,l}El'Uxp,Zp)

(20) when iXp,Zp) e {S',C,S}.

At the boundary C between S and S'the tangential components of both the electric and the magnetic field are to be continuous. This results into the conditions

(a^£^')(x,2))p=

iY^^'hx,z))p

when P £ C, (21)

In order to determine the distribution of E and 3 2? on C,

y n y

integral equations are derived by choosing in (18) and (20) the point of observation P on C. Using (21) we arrive at

4^(G,3 ff^'^ - £^'^3 G,)ds •^ £'^(Xp,2„) = |P^'^ (x„,2p),

'^C I ny y n \ y P' P ^ y P' P '

when iXpfZp) e C

(22)

when iXp,Zp) e C.

From (22), which is a coupled system of two Fredholm integral equations

of the second kind, the unknown distributions of E and 3 £" on

y ny C can be obtained. Substituting (15) in (18) and using (4),(5) and (10) the reflection factors are obtained as

R = iaV ) " ' i^(3 E^^\* - ff^'^3 i})'^)ds (23)

for m = 1, 2, .. . Using (11) the transmission factors are obtained as

T = 6 , + iaV )•"' ^^(3 E^^\~ - S^'^3 f)d3 (24) m m,l ^ m ^C ny m y rrm

for m = 1, 2, ... . In (24), 6 , denotes the Kronecker delta. The

coupled integral equations (22) are solved numerically using the method

of moments (Harrington, 1968) with N pulse functions as a set of

expansion functions for E and 3 E . By subsequent use of the

y n y ^

one-dimensional Dirac delta function as a testing function, the system

of integral equations is replaced by an approximating set of IN linear

r

algebraic equations, which are solved using matrix inversion. In order to save computation time, a technique for accelerating the convergence of G and its derivatives has been applied (De Jong, 1972).

3. VOLTERRA INTEGRAL EQUATION (V.I.E,) METHOD

The second method of solution starts with the derivation of an integral equation for the field distribution inside the obstacle instead of equations for the field and its normal derivative at the boundary of the obstacle. As to this method, the medium of the obstacle may be inhomogeneous in the x,2-plane, since no Green's function for the

Region S will be used. The solution of the scattering problem now starts from (2) which is rewritten as

(3^ -1-3^-1- k^E = - ik^ - khs . (25)

X z \' y ^ \ y

The incident field satisfies the two-dimensional Helmholtz equation

(3^ -1-3^-1- khs^ = 0. (26)

X z \ y ^ S u b t r a c t i n g (26) from (25) we o b t a i n

(3^ + d^ + khiE - E^) = - ik^ - khs . ( 2 7 ) X z ] ^ y y \ y

36 COMPUTATIONAL ASPECTS OF SCATTERING THEORY

Now, the solution of (27) is identical with the solution of

E iXp,Zp) = E\xp,Zp) + jj^ik (x,2) - k^)E ix,z)G^ix,z;Xp,Zp)dA, (28) where G denotes the Green's function for the uniform waveguide given

by (15). Equation (28), when iXp,z ) e S, constitutes a Fredholm

integral equation of the second kind for the unknown function E . This

equation can, in principle, be solved (Richmond, 1965), for instance using the method of moments. We shall however discard this possibility and transform (28) into a system of coupled integral equations of the Volterra type. The latter system has the advantage that it can be solved by an integration procedure. This starts with writing

E (x,2) = I _, t (2)sin(mTx/a) 0 < x < a ; - ' » < s < < » . (29)

Using (4), we now substitute (15) and (29) in (28), multiply the right-and left-hright-and side of (28) by sin(mïïXp/a), m = 1, 2, ..., and integrate the resulting equation over the interval 0 < Xp < a. As a result we obtain 2 I I °° exp(- r 2 - z„\) Z ,V (3)ijj (2)d2 3=3 ^ m' P' n=\ mn '^n (30) for m = 1, 2, ..,, - <» < 3p < °>, with

V iz) = iaT ) ! Ak ix.z) - fe,) s in (miTx/a) sin (mix/a) dx, (31)

— °° < 2 < °°,

We note that

V iz) = 0 when - < » < 3 < 3 , o r 2 ^ < 2 < ' » . (32)

mn 1 2 ^ '

Using (4) and splitting up the interval of integration, (30) can be rewritten as

f^P -1

ip (2„) = uAzA& , + u iz-A / u iz) Z ,V iz)é iz)Az -^ ^m P 1 P m,l m P J z=z. m n=\ mn n

(33) — 1 I 7 uo

•I- u izA / u iz) l" J iz)^ ( 2 ) d s , m P J z=z m^ ' n=\ mn ^n '

m = 1, 2, .,. . On account of (32), the representation (33) holds

every-where in the waveguide. Using matrices, (33) can more succinctly be written as l(3p) =£(3p)(6^'^ +13^3 £"'(3)Z(2)*(3)d3)

+ £~'(V 13=3 £(2)Z(2)i(2)d3,

(1) T — 1 2 and &_ = ( 1 , 0 , 0, ...) and the square matrices £(3) = (w iz)) with u iz) = u iz)& ,_{mn m mn =^ '} and Viz) = iV iz)). _{mn = ^ '} The matrix U iz) denotes the inverse of ILiz) . By writing in the first integral on the right-hand side

of (34)/ ^ =/ ^ - / ^, (34) can be written as (Kouri, 1973)

J 3 . J ^\ ' ^p

-1 r^2

i(3p) = £(3p)r -I- £ (3p)y g^3 £(2)Z(2)lJ;(2)d2

+

^ (35) 2.. ^ Ë ( V / ' = - = '(2)£(2)l^(3)d3 Z = Zr with

:(1)

J'2

y-1,

We note that r is an as yet unknown constant column vector. In order to obtain the solution of (35) WE

is defined as the solution of

obtain the solution of (35) we introduce the matrix f(3) = (é iz)), which

38 COMPUTATIONAL ASPECTS OF SCATTERING THEORY

3

i(3p) = £(2p)(i -ƒ gfg £''(2)£(2)l(2)d3)-H

where J^ denotes the unit matrix. Equation (37) constitutes a system of coupled Volterra integral equations of the second kind which can be solved by a direct numerical quadrature. Once (37) has been solved, the solution of (35) can be expressed as

i(3p) = i(3p)r, (38)

Substituting (38) into (36) the remaining unknown constant T is obtained

as

3

z ='l

T=

(J^ -ƒ

^l^ g\z)Viz)giz)dz)~^6}^\

(39)

To sum up, the scattering problem is solved by first integrating

(37), subsequently computing T from (39) and finally substituting these

results into (38), From the latter equation, the electromagnetic field

is finally obtained by using (29), When z^ < z < <», (35) reduces to

i(3p) = giZp)T, in Region 3, (40)

Using (29) and (4) we can now identify T as the column matrix containing

the transmission factors T that have been defined in (11), When

m

- oo < 3 < 2 (34) can, using (38) be written as

i(3p) = £(3p)6('^ + £~'(2p)P in Region 1 (41)

with

r^2

Using (41), (29), (10) and (4), P can be identified as the column matrix

containing the reflection factors R that have been defined in (9). For

the numerical solution, the expansion of E in (29) is truncated after

M terms. M„ is chosen such that a sufficiently accurate solution is

obtained. As a consequence of this truncation, the vectors in the above

analysis contain M elements and the square matrices have the dimension

M X M , As to the computation of V^, we note that, using the product formulas for goniometric functions, it can be reduced to the computation

of a vector containing 2M + 1 elements. For the numerical quadrature of

(37) the interval of integration is divided into Ny equal intervals,

Over each of these intervals (37) is integrated iteratively, using the Euler-McLaurin summation formula (Isaacson and Keller, 1966) as

integration rule,

4, DIFFERENTIAL EQUATION (D,E,) METHOD

The third method of solution does not make use of integral equations but is based entirely on the use of a suitable system of ordinary

differential equations. The medium of the obstacle may be inhomogeneous in the x,2-plane. In the uniform Regions 1 and 3 the field is expressed in terms of LSM -modes; the expressions for it are given in (10) and (11). In Region 2 the total field satisfies the source-free electromag-netic field equations (1). The first step towards the solution of the

problem consists of eliminating the longitudinal field component H

from (1), which yields

d H = - (jü)y„)~'3^£' + jtoEE" ,

3 X 0 X y "^ y '

in Region 2. (43) 3 £• = jüiVnH •

z y ^ 0 X

We shall now transform this system of partial differential equations into a system of coupled first-order ordinary differential equations. To this aim we write

40 COMPUTATIONAL ASPECTS OF SCATTERING THEORY

E (x,3) = Ï. ,a iz) sinin-nx/a) in Region 2, (44)

y n=\ n '^ U^ix,z) = (jojy^) Z^^jB^(3)sin(mTx/a) in Region 2. (45)

We now substitute (44) and (45) in (43) and perform the differentiations with respect to x. Multiplying the right- and left-hand sides of (43) by sin(mTTx/a), m = 1, 2, ..., and integrating over the interval 0 < X < a we obtain 3 6 (2) = E .W (3)a (2), z m n=\ mn n ' 3 a (s) = 6 iz), z m m in Region 2, (46) m = \ , 2 with

W iz)_{mn m,n} = 6 im-n/a) - _{v < / •'x=0 ^ ' ' ^ ^ / / >} (2/a) f .k (x,3)sin(mTrx/a)sin(wTrx/a)dx,

(47)

m,n = 1, 2, .... Using matrices (46) can be written as

3gB = V a,

in Region 2, (48) 3^a = 1.

T where we have introduced the column vectors a(3) = (a (2), 0.(2),...)

T — 1 2

and £(3) = (B (3), 6,(3),..,) and the square matrix £(3) = iW iz)).

The system of differential equations (48) must be supplemented by

boundary conditions at 3 = 3 and 2 = 2 , which follow from the continuity of_{•' y X} E and H at the planes 3 = 3 , and 2 = 2., Using (10),

"^ 1 2 6 V / >

(4) and (44) and the orthogonality of the sin(?TiTx/a) functions over

the interval 0 < x < a, the continuity of P at 3 = 3 yields

«m.l"l(^l) ^^m\ (^1> =-.(^l)- ^''^

Similarly, using (1), (10), (4) and (45) the continuity of H at 3 = 2 yields

r (- 6 ,uAz.)+ R M~'(3,)) = 6 (2,), (50)

m m, 1 1 1 m m \ m l

Elimination of R from (49) and (50) yields as boundary condition at

3 = 3

r a(3j) - i(3|) = 2£ £(2 J ) 5^*^' ^ (51)

where r = (r ) with r = r 6 , J/ and 6 have been defined in Section

= mn mn m m,n = —

3, As boundary condition at 2 = 3. we obtain in a similar way

£0(22) + 6(22) = 0. (52)

Now, (48) together with the boundary conditions (51) and (52) consti-tutes a two-point boundary-value problem. The solution {«(3), 6(2)} of this boundary-value problem is constructed by considering a related set of initial-value problems. The latter are obtained from the two--point boundary-value problem by removing the boundary condition (51) at 3 = 3 and replacing it by

J'Phz^) = &}P^

p = 1, 2 (53)

where 6 denotes a column matrix whose elements are zero except for the p-th element which equals one. The initial-value problem defined by (48), (52) and (53) can, for p = 1, 2 be solved yielding a set

of solutions {a ^ (3), 6^^^ (3)} p = 1, 2, ... , In general, none of

these solutions satisfies (51). However, since the problem is a linear one, a superposition of the solutions can be obtained such that (51) is satisfied. In order to obtain the latter solution we define the square

42 COMPUTATIONAL ASPECTS OF SCATTERING THEORY

matrices

A(2) = (a('^(2), a^^hz),..A; giz) = ( B ( ' ^ ( 3 ) , 6(^^(2),,.,), (54)

It follows from (53) that A^izA = J^ and consequently we can write the general solution {0^(3), _B(3)} of (48) as

a(3) = 4(3)0(22) (55)

and

6(3) = l(3)a(32) (56)

where aizA remains to be determined such that (51) is satisfied,

Substituting (55) and (56) into (51) we obtain an equation from which

aizA can be solved as

a(32) = 2(r4(2j) - I ( 3 , ) ) " ' L £ ( 3 J ) 6 _ ( ' \ (57)

Finally, R = (Pj, P2, ...)''' and T = (T^, T^,...)"^ can easily be obtained from (49) and the continuity of P at 3 = 2 as

R = £(2,)(_a(2j) - £(3,)6^('^) (58)

and

T = £"'(32)0(22). (59)

For the numerical solution, the expansions (44) and (45) are

trun-cated after M terms. As to the computation of the square matrix W^ we

note that it can be reduced to the computation of a vector containing 2M_ + 1 elements. For the numerical solution of the initial-value problems a fourth-order Runge-Kutta method (Keller, 1968) is employed,

5, NUMERICAL RESULTS AND CONCLUSIONS

The three different methods for the numerical solution of the scattering problem under consideration are applied to the configuration depicted in Fig,2. The angular frequency oj is chosen such that

0) = 1 ,5 X 0) , c

D

E, = E.= 8854xlO''V/ a E,=15E.

Fig,2, Longitudinal section of the configuration for which numerical results are presented,

where M denotes the cut-off frequency of the dominant LSM _-mode of the waveguide. As a consequence, there is only a single propagating mode in the uniform waveguide. Three computer programs for the solution of the scattering problem have been written by the authors, using the three different methods of solution outlined in the Sections 2, 3 and 4. In Figs. 3 and 4, the modulus of the reflection factor ff and of the

transmission factor T , respectively,are depicted as a function of D/a

(see Fig.2), In these figures we also have inserted results that have been obtained using a perturbation method that consists of substituting

the Rayleigh-Gans or Born approximation E = E in the right-hand

side of (28), The results obtained using the three methods differ relatively by less than 1 per cent, the relative difference between the results obtained using the V,I,E, method and the D,E. method are less than 0,1 per cent. In Table 1 the methods of solution are compared with regard to the computation time and storage requirements for the

compu-tation of the scattering by an obstacle with Dfa = 0,5. We observe that

the Fredholm integral equation method (F.I.E,) turns out to be compu-tationally very expensive. The comparison of the V,I.E.-method and the

44 COMPUTATIONAL ASPECTS OF SCATTERING THEORY

D,E,-method regarding time and storage requirements does not give a clear indication which method should be preferred. Regarding the problems encountered in formulating the problem and in writing the computer program, however, the D,E,-method is much easier to apply than

03-

02--^ <]

Born approx

Fig,3. Reflection factor |p | as a function of D/c

Method Computation time (seconds) Storage req. (kbyte) R F.I.E. iN„ = 36) r V.I.E. (My = 7) D.E. (Mp = 7) 200 448 160 160 0.20700 0.20682 0.20680

Table 1. Comparison of results for obstacle with D/a = 0.5,

the V,I,E, method. Due to its greater simplicity the D,E, method has the advantage that it can more easily be generalized to problems

involving coupled LSM- and LSE-modes, For large obstables, z„ - z = a

say, the D.E. method tends to show instabilities. It is expected that these difficulties can be coped with using other methods for solving the system of differential equations. The F,I,E. method does not offer an attractive alternative since the accompanying computation time would be in the order of 1 hour,

Finally, in order to demonstrate the applicability of the D,E, method to inhomogeneous, lossy obstacles we have computed the scattering by a cylindrical electron-plasma column, present in a waveguide (Cupini, Molinari and Poli, 1973). The configuration is depicted in Fig.5. In

C^M»

••^J.U'^

-(l-4000r*)) F/m

üXüJ-jV)

a=00229mjp=001m

Ü) =5O.5x10'racl/s;v=1.6.10'' s"*

46 COMPUTATIONAL ASPECTS OF SCATTERING THEORY

Fig,6. our results are compared with those given by Cupini et al. (1973), A good agreement is observed.

ACKNOWLEDGEMENT

The authors wish to thank Professor A.T, de Hoop for stimulating discussions and critical remarks concerning the research presented in this paper, 0.fr

t

Fig,6. Reflection factor |P.| as a function of the frequency ƒ = a)/2iT,

REFERENCES

COLLIN, R.E, (1960): Field Theory of Guided Waves, McGraw-Hill Book Company I n c , a: p.215; b: p, 198.

CUPINI, E., MOLINARI, V.G., and POLI, P. (1973: "Reflection coefficient of an electromagnetic wave by a plasma column of variable electron density in a waveguide", Alta Frequenza, vol. XLII, p,62-2E.

HARRINGTON, R,F, (1968): Field Computation by Moment Methods. The Macmillan Company.

ISAACSON, E., and KELLER, H,B, (1966): Analysis of Numerical Methods, John Wiley and Sons, I n c , p.287.

JONG, G. de (1972): "Scattering by a perfectly conducting cylindrical obstacle in a rectangular waveguide". Int.J,Elec-tronics, vol,32, p,153.

KELLER, H,B, (1968): Numerical Methods for Two-Point Boundary-Value problems, Blaisdell Publishing Company,

KOURI, D,J. (1973): "On the noniterative solution of integral equations for scattering of electromagnetic waves", J.Math. Phys,, vol,JI4_, p,l 1 16.

RICHMOND, J.H. (1965): "Scattering by a dielectric cylinder of arbitrary cross-section shape", IEEE Trans, on Antennas and Propagation, vol.AP-13, p.334,

A differential-equation method for the computation of the

electromagnetic scattering by an inhomogeneity in a cylindrical

waveguide

G. MUR

Department of Electrical Engineering, Laboratory oj Electromagnetic Research, Delft University of Technology, Mekelweg 4, Delft-IlOi, The Netherlands.

SUMMARY

In this paper an exact method is descrit)ed for compuiing numerically the scattering by an inhomogeneity in a cylindrical waveguide. The "Generalized Telegraphist's Equations" are used to transform the electromagnetic-field equations into a system of ordinary differential equations. The latter system behaves numerically unstable. A method is given to cope with this difficulty. Numerical results are presented for two- and three-dimensional obstacles in a waveguide of rectangular cross-section and they are compared with those obtained by other methods. Our method requires, in general, a relatively small amount of computation time and storage capacity. Another advantage of the method is its flexibility.

1. Introduction

In this paper an exact method is described for computing numerically the electromagnetic scattering properties of an inhomogeneity in the dielectric and/or the magnetic properties of the medium in an otherwise uniform, cylindrical waveguide (Fig. 1). In the waveguide, three different regions are distinguished, viz. the Regions I and 3 in which the medium is homogeneous and Region 2 in which the inhomogeneity in the properties of the medium is located. In the homogeneous regions of the waveguide, the fields can be expressed in terms of an uncoupled system of T£-, TM- and, if present, T£M-modes. In Region 2, where the inhomogeneity is located the "Generalized Telegraphist's Equations" [1] are used to transform Maxwell's partial differential equations into an infinite system of ordinary

•amaaam U n i t e r m w a v e g u i d e Region 1 | I n h o r ^ o g e n e o u s l / ' U n i f o r m tillecy w a v e g u i d e . w a v e g u i d e t I e 2 ( u , v , z ) , M , ( u ^ I Ej.uj ^ i^ I Region 2^:~ { Region 3

_,

^é^ 1

.—__

50 G. Mur

differential equations in z, the longitudinal coordinate. The presence of the inhomogeneity causes these differential equations to be coupled through z-dependent coupling coefficients. At either boundary between a homogeneous and the inhomogeneous region, continuity conditions apply. The relevant system of differential equations appears to be numerically unstable. In Section 4 we present a method to transform the system of differential equations into a related one with better stability properties as regards its numerical solution. To illustrate the applicability of the method, the scattering properties of obstacles of varying complexity, present in a waveguide with rectangular cross-section, are computed. The results are presented in Section 6.

The hitherto-known methods that deal with scattering problems of the type under consideration, are either approximate in nature, such as quasi-static approximations [2], or exact, using an integral-equation formulation of the problem [3,4, 5]. Although being exact, integral-equation methods have the disadvantage that they require the (often very time consuming) evaluation of dyadic Green's functions, and, in general, a large amount of storage capacity. Our method of solution is exact, too, but requires only a relatively small amount of computer time and storage capacity [6].

2. Description of the configuration

The configuration under consideration (Fig. I) consists of the three different regions that are defined in Table 1. In this table, also the permittivity e and the permeability ^ applying to these regions are indicated, as well as the relevant wavenumbers. We shall consider time-harmonic fields. The complex time factor expijwt), where j = imaginary unit, co denotes the angular frequency and t the time, is omitted in the formulas. Since the media in the waveguide may be lossy, the constitutive coefficients are in general complex with a negative imaginary part.

TABLE 1

Regions in the waveguide. Region Region 1 Region 2 Region 3 Location (u,v)eD — 00 < z < z Zj < Z < 00 Constitutive coeflficients £, /i of the medium C„£j(u, V, Z), ^„//j(u, 1!, z)

e„ and /!„ denote the permittivity and the permeability respectively, denotes the wavenumber in the reference medium.

Wavenumt)er Re(/()>0 k, =*„(£,/;,)»

of the reference medium; k„ •• = C0(£^„)» > 0

The electromagnetic fields in the configuration satisfy the source-free electromagnetic-field equations

V X H = jcoeE, V X £ = -jtu////, (1) where e = e(r) and /j = fi{r) denote the local value of the permittivity and the permeability,

« X £ = 0 (2) holds (n denotes the unit vector along the normal to the waveguide wall, as shown in Fig. I).

3. Derivation of the "Generalized Telegraphist's Equations"

As has been stated in the introduction, we shall in Region 2, where the inhomogeneity is present, transform the electromagnetic-field equations into a system of "Generalized Telegraphist's Equations". In the derivation of these equations we shall closely follow the line of thought due to Schelkunoff [1]. Thus, we first consider the modes of propagation in a guide filled with a homogeneous reference medium, not necessarily being vacuum. Each mode in this waveguide is described by the transverse distribution pattern T of either a potential or a stream function. To indicate the distinction between TE- and TM-modes we follow the convention adopted in [1], where the ordinal numbers for T£-modes are enclosed in brackets and those for TM-modes in parentheses. The T£Af-modes will be treated as a special case of the TM-modes. For convenience, we use a single-subscript notation, while arranging the modes in the order of non-decreasing cut-off frequencies. The function T{u, v), where u and v are suitable coordinates in a cross-section D with boundary C (see Fig. 1), is a solution of the two-dimensional Helmholtz equation in the cross-sectional domain. For T£-modes we have

(V^ + <i)^in, = 0 in D, (3) together with the boundary condition

a.7;„, = o on C. (4)

In (3), /Cj^i denotes the (real) cut-off wavenumber of the n-th T£-mode, while V^ = V —i^5^. For TM-modes we have

(V^ + <,)7;„, = 0 in D, (5) together with the boundary condition

?;„, = 0 on C. (6) In (5), ftTj^, denotes the (real) cut-off wavenumber of the n-th TM-mode. The T-functions, that

are chosen to be real, satisfy a number of orthogonality relations that are given in [1]; they are normalized such that

7;^,d-4

jjv,7;„,)(v,7;„,)d^ = K-(^

JD J JD

5 2 G. Mur The transverse components Ej and Hj of the electric- and the magnetic-field strength, respectively can be expressed in terms of the potential and stream functions, 'F and U for T£-modes, V and 77 for TM-modes, as

£T. = -VjV+i^ X VT-V, HT = -V-^U ~ i^ x VT.77 (8)

where i^ denotes the unit vector in the direction of increasing z. In general, the potential and stream functions can be written as

nu, <;, z) = - I J'[„,(z)7;„,("^, vl U{u, t), z) = - 17[„,(z)7;„,(t<. v),

(9) n « , V,Z)=-J: V^„,(z)7;„(u, i;), 77(u, v,z)=-X IJz)T,Ju, v),

n n

where the minus signs have been inserted to avoid a preponderance of minus signs in later equations. The Ks and the 7's of a mode represent the "voltages" and the "currents" in the equivalent transmission lines. Substituting (9) in (8) we obtain

£ T = X (ya')'^TTJu, V) - K X I^„,(z)V^7;„,(«, t>)), «T = S (/,„,(Z)VT7;„I(", ^) + K X /,„,(Z)VT7;„,(I<, f)).

n

We note that in obtaining (10), the medium in the reference waveguide, apart from being homogeneous, has not been specified. Furthermore we note that the expansion (10) of the transverse field components £y and H-^ is in terms of a complete sequence of transverse field distributions.

Having derived general expressions for the transverse field components in the reference waveguide, we now consider the fields in Region 2. To obtain expressions for the transverse field components in this region, we separate the transverse field components in the electromagnetic-field equations (1) from the axial ones and write the result as

5jH^ = V.p77. - jtuEi. X Ej, d.Ej = VjE, + yofii. x H-^, £ , = (jcu£)- ' V T • (H^ X i.), H. = {jcofir'V^ii.x E-,).

Asa consequence of the completeness of the sequence of expansion functions, the expression (10) for Ej and Hj can also be used in Region 2. Substituting (10) in (11) and eliminating the longitudinal field components £ . and H. we obtain

K ^ . - V I ' ^ T T ; , , + z('V,„,)'.-X v.r7;„,

= - j ' ^ ' ^ K K n , ^ X ^T7;„, + ^i„,^TÏ;„l) - I < , » ' , . , V T ( ( J ' " / ' ) ^ '7;„|), (12)

n n

i('\-K-.,)^r7;„,-i('\K.ii''--xVT7;„i

by the vector functions i^ x V.p7j^j and V^Tj^,, respectively. Secondly, we carry out a scalar multiplication of (12) by the vector functions Vj7j„, and i^ x V^-TJ^,, respectively. Subsequently, we integrate the four resulting equations over the cross-section D of the waveguide. As a result we obtain, using the orthogonality relations of the T-functions and the normalization conditions (7), a system of differential equations that can be written as [1]

n -ö,K,„,(z) = I (Z,„„.,(z)7,.,(z) + Z,„„.,(2)/,.,(z)), 11 n -8j^„p) = I (i^„,[„(z)^^„,(z) + >;„„„,(z)v;„,(z)), in Region 2 (14) with ^ M n l ( ^ ) = J " ' ' ^ 0 2[„„„,(2) = -i<^Mo 2,„,[„|(2) = -jCÜ/io 2,„,(„,(Z) = JG>/io A ' A x V r 7 ; „ , ) ( i , xV,7;„,)d/l, M2^TT,„,-(KxVTTi„,)dA, t'2^TTim)'^jTMdA + '(IK (15) and

W i ( ^ ) = •i'"^o J I «2VT7J„I • VT^i.id^ + K^fljJ

(ja)£oÊ,)-'7;„,7;„,dJ, D (ja)//o//,)-'7;„,7;„,d^, 51»1,„|(^) = j'^'^o

IJ

^m)(™)(z) = jÖ>£o '^zO'. xv^7;„,)-(i, xv^7;„,)dA

Now, (14) constitutes the system of generalized telegraphist's equations for the in-homogeneity in Region 2. These first-order ordinary differential equations are to be supplemented by boundary conditions at the planes z = z^ and z = Z2 bounding the waveguide section in which the inhomogeneity is located. Across these planes the tangential field components Ej and .77.^ are continuous. In order to formulate the relevant boundary conditions we must first obtain expressions for the transverse field components in Region 1 and Region 3.

54 G. Mur In Region 1 we assume that the field consists of incident modes travelling in the positive z-direction, having a voltage V" and a current 7', and scattered modes travelling in the negative z-direction having a voltage V' and a current 7'. The transverse field in Region 1 can therefore be written as

in Region 1 (17)

+ <> E (/;.,„,+ /!,,jv^7;, imy

The voltages and the currents of these (uncoupled) modes are related through

V' = 7 /' V' = — 7 I' V' =7 I' V = — 7 / '

' l . d n ) ^ l . l m ) ' 1 . ( 1 1 1 ) ' l . d n l ' ' 1 . ( m r 1, (m)'

(18)

In (18), Z j j^j and Z , („^ denote the wave impedances of the m-th TE- and TM-mode in Region 1, respectively. We have

^ i . M = J'^/^o^i/A.M' ^1.1».) = •^i.l-./iö'Vi' (19) where 7", ,^, and 7", ,„| denote the propagation coeflficients of the m-th TE- and TM-mode in

Region 1, respectively, i.e.

/".,,„, = « , - k\)' with Im(r,,,„,) > 0, ^^^^ A.,™, = « , - fci)* with I m ( T , , , „ , ) > 0 .

The propagation factor of a mode is exp( —Tz). In Region 3, the field consists of incident waves travelling in the negative direction and scattered waves travelUng in the positive z-direction. The transverse fields in Region 3 are written as

m -'•zXl(K',„,)+V-.,„,)V,7;„,, ffl m + ' z X K n , , ™ , + ^5.,™,)VT7;™,. m

where the voltages and the currents are related through

in Region 3 (21) V' = —7 I' V' = 7 P 'i.lm] -^3,(mi's.[m)' '3.|m) ^ 3 , [ m r 3 , [ m ] ' V' = — 7 ƒ' V' = 7 1' ' 3 , ( m ) ^ 3 . ( i > i r 3 . ( i t i ) ' ' 3 , ( m ) ^ 3 . ( m r 3 . ( m ) -(22)