International Journal of Infrared and Millimeter Waves, Vol. 27, No. 4, April 2006 (© 2006) DOI: 10.1007/s10762-006-9101-z

NUMERICAL SIMULATION OF OPEN WAVEGUIDE CONVERTERS USING FDTD METHOD Maxim L. Kulygin, Gregory G. Denisov, Alexey V. Chirkov and Sergey V. Kuzikov Institute of Applied Physics Russian Academy of Sciences 46 Ulyanov Street, Nizhny Novgorod 603950, Russia Received 28 November 2005 Abstract We

study

3-dimensional

asymmetric

diffraction

problems

for

waveguide-based electro-dynamic systems, radiating to infinite free space. For calculations we utilize the FDTD (Finite Difference Time Domain) numerical simulation method with the UPML (Unsplit Perfectly Matched Layer) absorbing boundary conditions. This paper states that the FDTD method, in spite of its relatively low calculation speed, has an approved ability of solving certain problems that cannot be solved by the other traditional numerical simulation methods (the method of integral equation, the method of scattering matrix). Keywords: FDTD, PML, UPML, open, waveguide, converter.

1. Introduction The article demonstrates a method of solving 3-dimensional asymmetric diffraction problems for waveguide-based electro-dynamic systems, radiating to infinite free space. The area of interest covers waveguide and mirror mode converters, open transmission lines and antennas used in radars, powerful generators and amplifiers of millimeter band radiation. An example of an open 591

592

Kulygin et al.

waveguide mode converter from an eigenmode to a Gaussian beam (without a phase corrector) is displayed at Fig.1. It represents a segment of a circular waveguide with a specifically cut end (we call it a cap).

Fig.1. Open waveguide mode converter from an asymmetric eigenmode to a Gaussian beam

Traditionally for calculation of open waveguide systems we use stationary numerical methods – the method of integral equations (IE), the method of scattering matrix and approximate analytical solution methods. Unfortunately these methods have shortages and are unable to work for a definite area of problems. We can see that the IE method works only with a quasi-optical approximation case, i.e. operates only high waveguide eigenmodes with a fixed polarization; the method of scattering matrix has the standard limitations of linear algebra, restricting the maximum number of eigenmodes and propagation directions in the whole system, also it has a problem of series convergence. These drawbacks are off for the FDTD (Finite Difference Time Domain) method [1, 2]. Till recently an insufficient level of computers did not let this method to work for electro-dynamic problems in a volume of about 1000 cubic wavelengths. Now the limitation has been taken off, and effective technologies of infinite free space modeling for the FDTD method have been created – perfectly [3, 4]. The FDTD method, in spite of its own drawbacks including relatively low calculation speed and precision and high consumption of

Numerical Simulation of Open Waveguide Converters Using FDTD Method

593

computer resources, significantly extends the area of problems solved by the traditional methods. The following problems are covered in the present article: •

test problem of open waveguide system for the FDTD method, that has the

known analytical solution [5] (see the Appendix); •

problems of the open waveguide converters from the symmetrical

eigenmodes of a circular waveguide to a Gaussian beam [6] at the limitation boundary of the IE method and outside the boundary; •

problem of the open waveguide converter from the asymmetric eigenmode

with low mode indexes (available for calculation using the FDTD method) to a Gaussian beam; •

more accurate definition of approximate analytical solution for the

converter from the asymmetric eigenmode to a Gaussian beam on a “strong deflection”. The goals of the presented study are: •

approbation of the software developed by the authors, using the FDTD

method with the UPML absorbing boundary conditions (see the Appendix); •

exact definition of limitations of the traditional numerical simulation

methods. 2. Circular waveguide converter of symmetrical modes The circular waveguide converter of symmetrical modes (Fig.2) consists of a segment of regular waveguide with a cap. Transversal section of the cap is half of a circle [6]. The length of the cap L is selected the way that outgoing radiation from the waveguide, according to the Brillouine conception, is unable to reach the outer surface of the cap. For real applications a close-to-Brillouine conception situation is obtained only for high waveguide eigenmodes. For lower eigenmodes a “diffraction leakage” is significant, to reduce it the cap is extended by several tens of percent (by 30% in this study). For given Brillouine angle T

45q the length of the cap is L 1.3d ˜ ctgT .

594

Kulygin et al.

Fig.2. Scheme of circular waveguide converter of symmetrical modes

For convenience the source frequency is selected equal to 30 GHz. 6 lowest symmetrical eigenmodes of the waveguide have been calculated – from TE01 to TE06. For the TE06 case the sizes of lattice are 360u360u870 nodes along all three spatial coordinates requiring about 450 Megabytes of RAM. The step value of the lattice ' (a distance between a couple of adjacent coordinate nodes, or a half of FDTD cell – the lattice period) is about 0.035 cm. So the vacuum wavelength of the source is equal to 29' or 14.5 cells. The waveguide radius is

R

d / 2 128'

4.415 cm,

the

cap

length

L

334' 11.5 cm,

the

waveguide thickness is always 2'. For the other eigenmodes the sizes of the system are decreased so the number of cells per wavelength is approximately constant (14.5…15). This condition defines an order of calculation precision.

(a)

(b)

(c) Fig.3. Diagrams of the field components in longitudinal (a), (b) and transversal (c) central sections of the system; (a) and (c) – transversal electric field, (b) – transversal magnetic field

Numerical Simulation of Open Waveguide Converters Using FDTD Method

595

Fig.3 displays the results of calculation for TE06 mode represented on three diagrams of different field components. A hard source [2] is located at the left end of the waveguide. The source radiates an eigenmode transversal structure propagating to the right. Fig.3 demonstrates an ordinary closed-waveguide solution nearby the source looking like “circular” areas near field maxima (white color) between field minima (black color). Coming away from the source the field structure transforms to a free-space solution looking like a spherical wave. At Fig.3 (a), (b) we can see a field minimum nearby the upper metallic wall. This is the place where the geometrical optics rays (shown at Fig.2) are finished according to the Brillouine conception. But some power still leaks from the cap. The 5-cell UPML absorbent (see the Appendix) is located at the boundary walls except the source area at the left. A specific UPML effect is clearly visible at Fig.3 (b), (c) close to the vertical metallic walls – the absorbent amplifies normal-to-boundary field components. This is a normal behavior of UPML allowing electromagnetic field components to satisfy Maxwell’s equations with no reflection from absorbent layers. For convenience of interpretation of the results and comparing them to ones obtained by the integral equation (IE) method let us consider field diagrams on unrolled surfaces of circular waveguide (not on a plane surface like above). These two numerical methods are quite different – the FDTD method operates all 6 time dependent field components and the IE method calculates only one stationary field component, here this is a longitudinal magnetic field (the most informative for TE-modes). Also the IE method divides the field component to a slow amplitude A and phase \ that correspond to a regular electric field component E of the FDTD method:

>

@

E FDTD (M , z ) { Re A IE (M , z ) ˜ exp(i\ IE (M , z )) ,

(1)

Where M is the azimuth, counted along the waveguide surface, z is the longitudinal coordinate. The values corresponding to amplitude and phase of the IE method, from time-dependent field component of the FDTD method at

596

Kulygin et al.

close-to-stationary state, are obtained the following way: A(M , z )

1 T

\ (M , z )

ª « arccot « « ¬

t 0 o f,

³

t0 T

t0

T

E (M , z , t ) ˜ dt , t0 T

³ ³

t0 t0 T t0

E (M , z , t ) ˜ cos(Z c t ) ˜ dt º» »  const E (M , z , t ) ˜ sin(Z c t ) ˜ dt » ¼

(2)

2S . Zc

Here T is a source oscillation period corresponding to the central frequency of the wave source Zc , t0 is an arbitrary moment of time for close-to-stationary state. All the integrals in (2) mean averaging by a period of the central frequency. The results similar to Fig.3 are displayed at Fig.4 at an unrolled surface of the waveguide. The upper pair of diagrams corresponds to the FDTD method. The lower one – to the IE method, amplitudes on the left and phases on the right. Here for convenience the longitudinal coordinate is originated from the beginning of the cap (before it all field components at the surface are zero). TE06 FDTD Amplitude

TE06 IE Amplitude

TE06 FDTD Phase

TE06 IE Phase

Numerical Simulation of Open Waveguide Converters Using FDTD Method

597

Fig.4. TE06 amplitude and phase diagrams for longitudinal magnetic field component at the unrolled surface of the circular waveguide calculated by the FDTD and IE methods

The color scale of amplitudes corresponds to maximum in black color and zero in white (a square area of constant white color to the upper left corresponds to the metallic half-cylinder cap). The color scale of phases is from -S to S, at the cap the phase is set to a constant equal to zero. The left diagrams look almost equal; to measure quantitatively the similarity we calculate the normalized integral of field component coupling: 2

S IE _ vs _ FDTD _ no _ phase

FDTD § (M , z ) ˜ A IE (M , z ) dzdM ·¸ ¨ M ,z A © ¹ . FDTD 2 IE (A (M , z )) dzdM ˜ ( A (M , z )) 2 dzdM

³

³M

,z

³M

(3)

,z

For TE06 S IE _ vs _ FDTD _ no _ phase = 99.34% demonstrating a good agreement between the results calculated by two different numerical methods. With phases taken into account we have (coordinate dependences are now omitted): 2

S IE _ vs _ FDTD _ with _ phase

§ A FDTD ˜ A IE ˜ exp(i (\ FDTD  \ IE )dzdM · ¨ M ,z ¸ © ¹ . FDTD 2 IE 2 (A ) dzdM ˜ ( A ) dzdM

³

³M

,z

³M

(4)

,z

For TE06 S IE _ vs _ FDTD _ with _ phase = 95.32% decreasing the methods’ agreement, this can be explained by a presence of two effects. The basic effect is caused by a staircase model of circular waveguide boundary approximation in Cartesian FDTD lattice. It can be indirectly confirmed by barely visible longitudinal stripes at the phase FDTD diagram. Also the decrease of agreement, especially for the lower modes, is explained by a specific error of the IE method caused by a deviation from the approximation of geometrical optics. The TE06 mode is high enough to the approximation of geometrical optics. At the same time the FDTD method can calculate this mode by consuming a reasonable amount of computational resources. These two reasons make the TE06 mode to be a “reference point” for two methods. A typical disagreement between two methods is shown at Fig.5 for TE01 mode.

598

Kulygin et al.

Fig.5. TE01 amplitude and phase diagrams for longitudinal magnetic field component at the unrolled surface of the circular waveguide calculated by the FDTD and IE methods

Several diffraction indexes have been calculated for the 6 lowest symmetrical TE modes: x

two coupling integrals given by (3) and (4);

x

relative power (Pointing vector integral when the source power is 100%) through the square under the cap, the square has a longitudinal size of 2 Brillouine lengths and azimuthal size of 180q;

x

coupling integral for the longitudinal magnetic field component with the field pattern of constant amplitude at the square meant above, without phase and with a linear propagating phase corresponding to the regular waveguide.

The numerical results for all 6 indexes are displayed at Fig.6. The best agreement is reached for amplitude without phase, the agreement decreases for lower modes as explained above. Taking phase into account we face another decrease of agreement due to the FDTD approximation error. All indexes demonstrate similar tendencies – the diffraction effects increase for lower

Numerical Simulation of Open Waveguide Converters Using FDTD Method

599

modes.

Fig.6. Six quantitative diffraction indexes for the six lowest eigenmodes

3. Circular waveguide converter of non-symmetrical modes

The non-symmetrical mode converter differs from the symmetrical one by a form of the cap. The cap is now formed by two helical lines [6]. The generic non-symmetrical mode converter looks like that shown at Fig.1. The Brillouine area to the cut can be selected arbitrary. Here we consider a simplified one where one of the helical lines is a straight line parallel to the waveguide’s axis of symmetry. The working mode is now TE53, this mode is frequently used in microwave applications and it is currently available for calculation by the FDTD method. Fig.7 displays an unrolled surface of non-symmetrical

600

Kulygin et al.

waveguide converter.

Fig.7. Unrolled surface of non-symmetrical waveguide converter

The source frequency is 30 GHz, the Brillouine angle is 45q. So the waveguide radius is R

3.1482236 cm 128' , '

0.0246 cm. According to [6] for the

selected values of parameters the cap length is L 15.328 cm, the rise angle of the second helical line is D = 0,911556 = 52,228q. For FDTD calculations we use a lattice of 337 u 337 u 1123 cells requiring about 500 Megabytes of RAM. The results for the longitudinal magnetic field component at the waveguide’s unrolled surface are shown at Fig.8. The diagrams demonstrate the situation close to geometrical optics when almost all outgoing power is located under the cap and the phase corresponds to a propagating wave with almost constant phase velocity. TE53 H|| Amplitude

TE53 H|| Phase

Fig.8. TE53 amplitude and phase diagrams for longitudinal magnetic field component at the unrolled surface of the circular waveguide calculated by the FDTD method

For the same converter Fig.9 shows the diagrams for longitudinal electric

Numerical Simulation of Open Waveguide Converters Using FDTD Method

601

field component which describes the effect of radiation cross-polarization in the system due to symmetry disturbance caused by the cap. A significant part of cross-polarization power is located near the edge of the cap and it has a structure resembling spherical wave. TE53 E|| Amplitude

TE53 E|| Phase

Fig.9. TE53 amplitude and phase diagrams for longitudinal electric field component (cross-polarization) at the unrolled surface of the circular waveguide calculated by the FDTD method

The relative coefficient of cross-polarization is calculated as follows: CTE 53

³M ³M

2

,z

,z

E|| dM ˜ dz 2

H || dM ˜ dz

.

(5)

Here the integration is performed at the whole unrolled surface. For our converter CTE 53 = 0.5% which is in a good agreement with experimental results for typical converters. 4. Precise calculation of the “strong deflection” waveguide converter

A waveguide converter on a “strong deflection” [8] (Fig.10) represents one period of sinusoidal corrugated circular waveguide placed between two segments of regular circular waveguide. The input of the converter is fed by an eigenmode of regular circular waveguide. The converter is intended to “focus” a transversal field structure inside the output waveguide at a distance from the period of corrugation. The focused multi-mode structure, where field values at metallic walls are close to zero, can be used to make systems with an open gap or waveguide corners [8]. In the presented study the input mode is TE53, the

602

Kulygin et al.

source frequency is 30 GHz, the Brillouine angle is T

40q . Therefore the

radius of the circular waveguide is R = 3.4632 cm.

Fig.10. Longitudinal section of the “strong deflection” waveguide converter

From [8] for corrugated segment we have an approximate analytical solution:

d

h54  h52 2S h , , 2 h 2ka cos T | 0.9.

(6)

Here d is a period of corrugation, a is a corrugation depth, h54 and h52 are longitudinal wave numbers for TE54 and TE52 modes of circular waveguide, k is a full wave number. For these values of parameters the period of corrugation is

d

7.419 cm, the corrugation depth according to the analytical solution is

a0

0.0935 cm. The FDTD is required to define the optimal corrugation depth

(or “deflection value”) for which the maximum focusing can be obtained inside the output waveguide (with the minimum field value at the walls). We have made several calculations for different corrugation depths. The length of the input waveguide is chosen equal to 50 ' for FDTD convenience reasons. The output waveguide length should exceed the distance of focusing. According to [8] its approximate value is: L|

1 2 m2 kR cos 3 T ˜ (1  2 ) , 2 Q

(7)

where m is an azimuthal index of the input eigenmode, Q is its Bessel’s root. For this converter we have L | 15 cm. The considered problem is actually 2-dimensional due to an axial symmetry

603

Numerical Simulation of Open Waveguide Converters Using FDTD Method

of its geometry and all existing modes have the same azimuthal index value. However for approbation of the method authors use generalized 3-dimensional FDTD procedures similar to the mentioned above. The results of calculation for the value of corrugation depth given by the analytical solution (6) are displayed at Fig.11.

(ɚ)

(b)

(c)

Fig.11. Diagrams of a longitudinal magnetic field component at longitudinal (a) and transversal sections: (b) – “defocusing”, (c) – “focusing”, a = 0.0935 cm

Fig.12. Absolute field value along the straight line close to the surface, a = 0.0935 cm

Fig.11 (a) shows a diagram of a longitudinal magnetic field component at the central longitudinal section of the system. Small white arrows placed beyond the waveguide walls mark position of transversal field structure “defocusing” and “focusing”, the corresponding diagrams are shown at Fig.11 (b), (c). These positions are also marked at Fig.12 where the plot of the longitudinal field component is shown in dependence of longitudinal coordinate, close to the

604

Kulygin et al.

waveguide surface. The vertical axis means an instant field value in dimensionless units at close-to-stationary state. The dashed line shows the place where the corrugated segment is located. The area with zero field value corresponds to the concave of corrugation. It is obvious that the “focusing” of the transversal field structure is not enough in this case, the amplitude value in “focused” section is only 5 times less than the one in the regular input waveguide.

(ɚ)

(b)

(c)

Fig.13. Diagrams of a longitudinal magnetic field component at longitudinal (a) and transversal sections: (b) – “defocusing”, (c) – “focusing”, a = 0.14 cm

Fig.14. Absolute field value along the straight line close to the surface, a = 0.14 cm

Fig.13, 14 display similar results for the case a = 0.14 cm (50% increase from the theoretical value). For the corrugation depth value of 0.14 cm we can see “overfocusing”; it is obvious that the optimal corrugation value is between

Numerical Simulation of Open Waveguide Converters Using FDTD Method

605

0.0935 cm and 0.14 cm. For the value of 0.117 cm the “focusing” is optimal (Fig.15, 16).

(ɚ)

(b)

(c)

Fig.15. Diagrams of a longitudinal magnetic field component at longitudinal (a) and transversal sections: (b) – “defocusing”, (c) – “focusing” a = 0.117 cm

Fig.16. Absolute field value along the straight line close to the surface, a = 0.117 cm

Fig.16 also displays an envelope to a solution to the same problem obtained by D.A. Lukovnikov using the method of coupled waves. Both solutions are almost equal except the right maximum of the envelope. This is because the slight difference of the problem formulation for two methods – the method of coupled waves solves the problem with a longer segment of the output waveguide. The noticed solution agreements prove the ability of the FDTD method to solve

606

Kulygin et al.

open problems. 5. Conclusion

A new software technology of solving open waveguide problems has been successfully approved in this study. For numerical simulations we use the FDTD method with the UPML absorbing boundary conditions. The high efficient UPML technique allows a 5-cell layer to achieve reflection coefficient less than -60 dB in broad band of parameters of incident radiation. At the same time the relative calculation performance decrease caused by the UPML is less than 1% for typical problems. The typical calculation times are about a few hours for a regular PC. The FDTD method is able to solve open waveguide problems with working volumes of about 1000 cubic wavelengths. It can normally operate either symmetrical or non-symmetrical waveguide eigenmodes with coordinate indexes less than 10. The results of calculations are in good agreement with other traditional numerical methods (if possible), the difference is less than 1% for the IE method at its found limitation boundary. Besides the FDTD method obtains much more information since it operates all 6 components of electromagnetic field. In particular the FDTD method calculates a cross-polarization. The FDTD method allows solving problems with reasonable precision at the areas where traditional numerical methods are unable to work. Acknowledgments

The authors thank D. A. Lukovnikov for the method of coupled waves solution, A. M. Gorbachev and V. A. Koldanov for valuable discussions.

Numerical Simulation of Open Waveguide Converters Using FDTD Method

607

Appendix: UPML absorbing boundary conditions for the FDTD method The FDTD method consumes lots of computational resources – random access memory and processor time. According to the condition of physical approximation meaning, one should use not less than 10 nodes of discrete lattice per a typical wavelength of the system [1, 2 and 7]. The modern level of computer development makes it possible to calculate 3-dimensional electro-dynamic problems with typical sizes of about several 10s of wavelengths per each of three spatial dimensions. Therefore the whole computational domain is rather limited and lattice truncation surfaces (usually perfect electric conductor - PEC boundary conditions) cannot be located far enough from objects of calculation. So reflected radiation is returned back making significant interference at rather small time values. Therefore, to model free space the FDTD method requires matched absorbing boundary conditions. The most effective of them are perfect matched layers (PML) [3, 4] providing small longitudinal sizes of several layers and low reflection coefficients for a wide spectrum of frequencies, incident angles and transversal structures. This article uses the UPML (Unsplit Perfectly Matched Layer) [4] absorbing boundary conditions. One of the advantages of the UPML (i.e. over standard PML [3]) is its satisfaction to the Maxwell’s set of equations: G G G G 1 wB 4S Zˆ G rotE  rotH  V m H i P H , c wt c c

G 4S G 1 wD VE  c c wt

i

Zˆ G H E.

Here V and V m mean tensors of electric and magnetic conductivity, iZˆ

c

(8)

w is an wt

operator of time derivative. The UPML layer located between vacuum area and truncation surface should not reflect anything, so the tensors of permittivity and permeability are:

P

H

§ s y sz ¨ ¨ sx ¨ ¨ 0 ¨ ¨ ¨ 0 ¨ ©

0 sx sz sy 0

· 0 ¸ ¸ ¸ 0 ¸˜ ¸ sx s y ¸ ¸ s z ¸¹

(9)

608

Kulygin et al.

The dimensionless parameters sD are responsible for wave propagation along 1

corresponding Cartesian coordinates, so sD

4SV D (D ) . An operator standing at iZˆ

the denominator means a backward operator. For convenience we omit the operator symbol “^” at further expressions.

sx

Let

us

sy

1, s z

consider 1

a

simplified

4SV , where V iZ

case

of

1-dimensional

UPML:

V z (z ) . This UPML can be useful as a load for

an electro-dynamic system with a waveguide of constant transversal section parallel to z-axis at its output. The tensor (9) is the following now:

P

§ ¨ ¨1  4SV ¨ iZ ¨ ¨ 0 ¨ ¨ 0 ¨ ¨ ©

H

0 1

4SV iZ 0

· ¸ ¸ 0 ¸ ¸ 0 ¸˜ ¸ 1 ¸ 4SV ¸ 1 ¸ iZ ¹

(10)

Placing (10) to (8) for transversal field components we have similar equations like 1§ w · ¨  4SV ¸ E x c © wt ¹

rotHG

x

, that transforms to an ordinary FDTD scheme with

conductivity: 1 E xt 1 / 2  E xt 1 / 2 c 't E xt 1/ 2



G 4S VE xt  rotH x . c





1  2SV't t 1/ 2 1 c't H zy 1  H zy 1  H yz 1  H yz 1  Ex 1  2SV't 1  2SV't 2'

Considering average E xt |



E

t 1 / 2 x

 E xt 1 / 2 2



(11)



(12)

we have the final equation (12). In

expressions (11) and (12) upper subscripts mean coordinate and time values that differ from current 4-vector values of (x, y, z, t) in a part of lattice period, dimensionless for convenience; ' is a distance between adjacent lattice nodes (lattice period is 2'), 't is a time step. For all of 4 transversal field components the FDTD schemes are similar.

Numerical Simulation of Open Waveguide Converters Using FDTD Method

609

However, due to nonlinear dependence of longitudinal tensor component on conductivity, the longitudinal field components expressions are obtained another way. Let us consider a longitudinal component of electric flux density: 1

§ 4SV · ¨1  ¸ Ez . iZ ¹ ©

Dz

(13)

Placing (13) to (8) we obtain a finite difference equation corresponding to vacuum for Dz component:

D zt 1 / 2

D zt 1 / 2 

c't x 1 H y  H yx 1  H xy 1  H xy 1 . 2'





(14)

These longitudinal electric and magnetic flux density components corresponding to the UPML area should be retained in computer’s random access memory like the ordinary 6 field components. Using (13) we have:

wDz  4SVDz , wt

wE z wt E zt 1 / 2

(15)

E zt 1 / 2  D zt 1 / 2 (1  2SV't )  D zt 1 / 2 (1  2SV't ) .

(16)

Equations (14) and (16) let us calculate the longitudinal field components in 1-dimensional UPML. The general case of 3-dimensional UPML the tensor (9) is equal to a commutative multiplication of three “1-dimensional” tensors similar to (10). It means that operators containing Z may be applied in any order. The way of obtaining finite difference equations is similar for all field components now. Here we describe only the x-component from (8):

§ 4SV y iZ ¨¨1  iZ © Considering D1x

§ 4SV y ¨1  ¨ iZ ©

an

1

rotHG .

·§ 4SV z ·§ 4SV x · ¸¨1  ¸ Ex ¸¨1  ¸© iZ ¹© iZ ¹ ¹ electric

flux

x

density

(17)

component

1

·§ 4SV z ·§ 4SV x · ¸¨1  ¸ Ex ¸¨1  ¸© iZ ¹© iZ ¹ ¹

we obtain an ordinary “vacuum”

equation:

D1t x1 / 2

G D1t x1 / 2  c't ( rotH ) x .

(18)

610

Kulygin et al.

4SV y § ¨1  ¨ iZ ©

Considering another auxiliary vector component D2 x

·§ 4SV z · ¸¨1  ¸ E x we ¸© iZ ¹ ¹

have: D2t x1 / 2

D2t x1 / 2  D1tx1 / 2 (1  2SV x 't )  D1tx1 / 2 (1  2SV x 't ) .

Further we have to consider another vector component D3 x D3t x1 / 2

§ 4SV y ¨1  ¨ iZ ©

(19)

· ¸ E x , so: ¸ ¹

1  2SV z 't t 1 / 2 1 D3 x  D1tx1 / 2 (1  2SV x 't )  D1tx1 / 2 (1  2SV x 't ) (20) 1  2SV z 't 1  2SV z 't

>

@

From definition of D3x we have: E xt 1 / 2

1  2SV y 't 1  2SV y 't

E xt 1 / 2 

>

@

1 D3t x1 / 2  D3t x1 / 2 . 1  2SV y 't

(21)

The system of equations (18), (20) and (21) allows us to calculate any of 6 electromagnetic field components in 3-dimensional UPML layer by a cyclic change of the coordinate subscripts. All D1 and D3 components need to be stored in RAM (D2 doesn’t need storing). Usually the problem of storing is simplified using a single FIFO storage array for all auxiliary components because their sequence order is invariable in time. The conductivity components VD depend on their spatial coordinates and responsible for fading along corresponding coordinate axes. They usually have the m

following dependence V a (D )

§D · ¸ , 0  D  d , where d is the layer thickness, ©d¹

V max ¨

usually it is constant for all coordinates counted from a UPML-vacuum boundary. The dependence can be selected freely, the authors use power dependence with UPML order m = 4. The maximum conductivity constant Vmax corresponds to a conductivity value at a UPML-PEC boundary (at the FDTD lattice truncation plane). It can be chosen the way when the minimum thickness of the layer provides no reflection from a PEC. In [4] estimation for the maximum conductivity value is provided, we show it in the CGS system of units:

V max

2˜

1 c't m  1 , ˜ ˜ 2S't 2' N UPML

(22)

Numerical Simulation of Open Waveguide Converters Using FDTD Method

611

where NUPML is a double number of layers (the layer thickness in units of '). This estimation has been experimentally defined more precisely by the authors. The numerical experiment assumes a generic problem of a waveguide excitation by different eigenmode sources. The waveguide is loaded with a 1-dimensional UPML layer. A reflection coefficient is calculated using a standing wave coefficient. Using a half-deflection method the optimal value of Vmax is acquired:

V max

3.5 ˜

1 c't m  1 ˜ ˜ . 2S't 2' N UPML

(23)

With a thickness layer of 5 cells the reflection coefficient is less than -60dB in a wide band of frequencies, incident angles and transversal structures of incoming radiation. The optimized UPML algorithm developed by the authors has a fast, compact and effective source code based on polymorphism and inheritance links to a standard vacuum FDTD procedure. The 5-cell UPML algorithm causes a relative performance decrease of less than 1% for typical 3-dimensional problems with a volume of about 1000O3. The developed UPML algorithm has been tested on a well-known problem of eigenmode diffraction at an open end of circular waveguide. The problem has an accurate analytical solution by L. A. Weinstein in [5]. A corresponding test problem for the FDTD method with UPML is formulated by the following: we consider a circular waveguide radiating to a parallelepiped limited by PEC planes. Each PEC plane is covered by a UPML layer. In [5] the solution to the problem is given by a directional radiation pattern. The graphs for by-power normalized “electric” V (- ) and “magnetic” V~ (- ) radiation characteristics, proportional to square values of some field components, are displayed at Fig.17.

612

Kulygin et al.

Fig.17. Directional radiation patterns for TM11 eigenmode diffraction at an open end of circular waveguide [5] Here the direction -

0q corresponds to backward scattering, - 180q - to forward

scattering. The diagrams (instant photographs) of field components corresponding to “electric” and “magnetic” field characteristics, calculated by the FDTD method, are displayed at Fig.18. On the left of both diagrams a segment of radiating waveguide is displayed at the central longitudinal section passing through the waveguide’s axis of symmetry. All truncation PEC planes are covered by UPML layers. 90q 0q

V (- )

180q

V~ (- )

Fig.18. Radiation characteristic diagrams calculated by the FDTD method There’s no visual reflection from PEC boundaries on both diagrams. Performing a qualitative comparison with the results from Fig.17 one can see an obvious coincidence of all significant extrema of the field characteristics. A quantitative comparison is difficult because Fig.17 displays a pattern at an infinite distance from the end of waveguide but the FDTD method can operate only limited distances. So this comparison approves the developed UPML-FDTD algorithm for solving open

613

Numerical Simulation of Open Waveguide Converters Using FDTD Method

problems.

References [1] D. V. Vinogradov and G. G. Denisov, “Mode conversion in curved waveguide with variable curvature”, Izv. VUZov Radiofizika, v.33, No.6, p.726, 1990. [2] M. Thumm. “Modes and mode conversion in microwave devices”, Proc. of the 48th Scottish universities summer school in physics: Generation and application of high power microwaves, St. Andrews, a NATO Advanced study institute, p.121, August 1996. [3] D. V. Vinogradov and G. G. Denisov, “Waveguide mode converter”, USSR patent SU 1566424 A1. [4] K. S. Yee, “Numerical Solution of Initial Boundary Value Problems Involving Maxwell's Equations in Isotropic Media”, IEEE Transactions on antennas and propagation, Vol.AP-14, No.8, p.302, 1966. [5] A.

Taflove

and

S.C.

Hagness,

“Computational

Electrodynamics:

The

Finite-Difference Time-Domain Method”, Artech House, Boston and London, 854 pages, 2000. [6] M. L. Kulygin, “Calculation of dispersion characteristics of circular waveguide with deep helical corrugation using FDTD method”, Izv. VUZov Radiofizika, v.47, No.1, p.69, 2004. [7] J. P. Berenger, “Three-dimensional perfectly matched layer for the absorption of electromagnetic waves”, Journal of Computational Physics, v.127, p.363, 1996. [8] L. A. Weinstein, “Electromagnetic waves”, Radio and communications, Moscow, 340 pages, 1988.