T I 034702-1-T PERFECTLY MATCHED ABSORBERS (PML) THEORY, ANALYSIS AND IMPLEMENTATION Y. Botros, S. Legault, J. Gong J. Volakis and T. Senior Quarterly Report Compact Software, Inc. 201 McLean Boulevard Paterson, N.J. 07504 October 1996 THE UNIVERSITY OF MICHIGAN Radiation Laboratory Department of Electrical Engineering and Computer Science Ann Arbor, Michigan 48109-2122 USA 34702-1 -T = RL-2473

PROJECT INFORMATION PROJECT TITLE:: Accelerated 3D Arbitrary EM Simulation Methods REPORT TITLE: Perfectly Matched Absorbers(PML): Theory, Analysis and Implementation U-M REPORT No.: 034702- 1 CONTRACT START DATE: END DATE: 1 August 1996 31 July 1997 DATE: SPONSOR: SPONSOR CONTRACT No.: U-M PRINCIPAL INVESTIGATOR: November 1996 (1st Quarterly Report) Jason Gerber Compact Software Inc. 201 McLean Ave. Paterson, N.J 07504 96-0502 John L. Volakis EECS Dept. University of Michigan 1301 Beal Ave Ann Arbor, MI 48109-2122 Phone: (313) 764-0500 FAX: (313) 747-2106 volakis @umich.edu http://www-personal.engin.umich.edu/-volakis/ CONTRIBUTORS TO THIS REPORT: Y. Botros, S. Legault, J. Gong, J. Volakis, T. Senior

Forward This document represents our first quarterly report on this project. The report is rather extensive (over 40 pages) and covers all aspects relating to the theory and performance of the perfectly matched absorbers (PML) for mesh truncation. It contains data and research background on the PML when used in conjunction with the finite element method. During this quarter we considered the implementation of the PML for truncating microwave structures and Figures 16, 17, 19 and 20 are examples of such applications. The optimization of the PML (i.e. maximum performance for minimum thickness) is a continuation of work which began before the start of this contract. This study is very significant since we concluded that a simple formula can be used for obtaining such information. Furthermore, we were able to confidently conclude that although the PML is highly effective as an absorber, its improved performance came at price. The convergence of the FEM system deteriorated and the higher sampling is need within the PML region because of the fields rapidly decaying behavior. The issues of convergence and increased sampling requirements were of primary concern to. us over the past month. It is important that the convergence be improved and the PML is implemented using minimal resources. At this point we are investigating two possible approaches to address these issues. Our next report will address these issues and related progress in some detail.

Perfectly Matched Absorbers (PML): Theory, Analysis and Implementation Youssry Y. Botros, Stephane R. Legault, Jian Gong John L. Volakis and Thomas B.A. Senior Radiation Laboratory Department of Electrical Engineering and Computer Science, The University of Michigan Ann Arbor, MI 48109-2212 Abstract The recently introduced perfectly matched uniaxial absorber layer (PML)for frequency domain finite element simulations has several advantages. In this report we present the application of PML for microwave circuit simulations along with design guidelines to obtain a desired level of absorption. In addition to providing accuracy control and conformality, the PML absorbers are also easy to implement and result in codes with no parallelization bottlenecks. Several example applications to three dimensional structures are included to demonstrate their performance. 1

Contents 1 Introduction 2 Formulation 3 5 3 Optimization Study 6 3.1 Numerical Model............................... 7 3.2 Dependence on a and i......................... 8 3.3 Dependence on Q, N andt................................ 9 3.4 Design Curves............................... 11 4 Applications 12 4.1 Rectangular Waveguide....................... 13 4.2 Microstrip line............................. 14 4.3 Band Eliminator............................ 15 4.4 Sphere Scattering........................... 16 4.5 Meander Line.............................. 17 5 Conclusions 5.1 Appendix: FEM Formulation...................... 18 22 2

1 Introduction One of the most important aspects of finite difference and finite element implementations is the truncation of the computational domain. An ideal truncation scheme must ensure that outgoing waves are not reflected backwards at the mesh termination surface, i.e. the mesh truncation scheme must simulate a surface which actually does not exist. To date, a variety of non- reflecting or absorbing boundary conditions (ABCs) have been employed for truncating the computational volume at some distance from the radiating or scattering surface, and applications to microwave circuits and devices have also been reported. The ABCs are typically second or higher order boundary conditions and are applied at the mesh termination surface to truncate the computational volume as required by any PDE solution. Among them, a class of ABCs is based on the one-way wave equation method [1, 2] and another is derived starting with the Wilcox Expansion [4, 5]. Also, higher order ABCs using Higdon's [6, 9]formulation and problem specific numerical ABCs have been successfully used, particularly for truncating meshes in guided structures [14]. There are several difficulties with traditional ABCs. Among them is accuracy control, conformality, ease of parallelization and implementation difficulties when dealing with higher order ABCs. Also, the applications of ABCs in microwave circuit modeling requires a priori knowledge of the propagation constants which are typically not available for high density packages. An alternative to traditional ABCs is to employ an artificial absorber for mesh truncation. Basically, instead of an ABC, a thin layer of absorbing material is used to truncate the mesh, and the performance for a variety of such uniform isotropic absorbers has already been considered [13] and [19]. Nevertheless, these lossy artificial absorbers (homogeneous or inhomogeneous) still exhibit a non-zero reflection at inci 3

dence angles away from normal. Recently, Berenger [16] introduced a new approach for modeling an artificial absorber that is reflectionless at its interface for all incidence angles. In two dimensions, his approach requires the splitting of the field components involving normal (to the boundary) derivatives and assigning to each component a different conductivity. In this manner an additional degree of freedom is introduced that is chosen to simulate a reflectionless medium for all incidence angles. Provided the medium is lossy, this property is maintained for a finite thickness layer. Berenger refers to the latter as a perfectly matched layer (PML) and generalizations of his idea to three dimensions have already been considered [17] and [18]. Implementations of the absorber for truncating finite difference-time domain(FDTD) solutions have so far been found highly successful. Nevertheless, it should be noted that Berenger's PML does not satisfy Maxwell's equations and cannot be easily implemented in finite element (FEM) solution. In this report we examine a new uniaxial anisotropic artificial absorber [15] for truncating FEM meshes. This artificial absorber is also reflectionless at all incidence angles. Basically, by making appropriate choices for the constitutive parameter tensors, the medium impedance can be made independent of frequency, polarization, and wave incidence angle. A PML layer can then be constructed by introducing sufficient loss in the material properties. The implementation of this artificial absorber for truncating finite element meshes is straightforward and, moreover, the absorber is Maxwellian. Below, we begin with a brief presentation of the proposed artificial absorber, and this is followed by an examination of the absorber's performance in terminating guided structures and open domain volume meshes. Results are presented which show the absorber's performance as a function of thickness/frequency and for different loss factors. 4

2 Formulation Consider the waveguide, shielded microstrip line and scatterer shown in Figures 2 and 3. We are interested in modeling the wave propagation in these structures using the finite element method. For a general anisotropic medium, the functional to be minimized is = i v x E (7 V xE) - k a E EdV E x V x E) dS, (1),in +Sout in which r and r denote the permeability and permittivity tensors whereas E is the total electric field in the medium. The surface integrals over Sin, and Sout must be evaluated by introducing an independent boundary condition and the ABC serves this purpose. Alternatively, an absorbing layer may be used. In this study we consider the performance of a thin uniaxial layer for terminating the FE mesh in a rectangular waveguide, a microstrip line and for open domain scattering problems. Such a uniaxial layer was proposed by Sacks et.al. [15] who considered the plane wave reflection from an anisotropic interface (see Figure 1 ). If r and Er are the relative constitutive parameter tensors of the form a2 0 0 r = r = 0 b2 0 (2) 0 0 C2 the TE and TM reflection coefficients at the interface (assuming free space as the background material) are TE cos0 - cos Ot RTE t cosOi + a cos(3 TM cost - cosOi RTM a cosOi +. cosOt 5

9 By choosing a2 = b2 and c2 = - it follows that RTE = RTf = 0 for all incidence angles. implying a perfectly matched material interface. If wve set a2 = Q - j3. the reflected field for such a metal-backed uniaxial layer is I R(o,)= - (4) where t is the thickness of the layer and 0. is the plane wave incidence angle. The parameter ak is simply the wavenumber in the absorber. Basically, the proposed metal backed uniaxial layer has a reflectivity of -30 dB if 3tcosO, = 0.275A or -55dB if n3tcosi =' 0.5A, where A is the wavelength of the background material. The reflection coefficient (4) can be reduced further by backing the layer with an ABC rather than a PEC. However, the PEC backing is more attractive because it eliminates altogether the integrals over the surfaces. Clearly, although the interface is reflectionless, the finite thickness layer is not and this is also true for Barenger's PML absorber. Below we present a number of results which show the performance of the proposed uniaxial absorbing layer as a function of the parameter,3, the layer thickness t and frequency for the guided and scattering structures shown in Figures 2 and 3. We remark that for the microstrip line example it is necessary to let a2 = ~rb({ - j3) for the permittivity tensor and a2 =,Urb( - j3) for the permeability tensor, where Erb and Lrb are the relative constitutive parameters of the background material (i.e. the substrate). 3 Optimization Study As implied from the previous section, we can use different types of layers for mesh terminations, namely, the isotropic, the anisotropic backed with PEC and that backed with ABC. For these three types of layer the theoretical behavior of \RI is relatively simple. In the case of the isotropic material, an increase in f3 and/or t decreases 6

IR(O)]. The uniaxial material has this behavior for all real angles of incidence, and while a plays little or no role, large values of a do produce higher absorption for complex angles. It follows that for a uniaxial layer of given thickness t a and;3 can be chosen sufficiently large to produce high absorption over a broad angular spectrum with angles near =- ir/2 (see Figure 4) providing the only exception. Unfortunately, the analytical results do not immediately translate into numerical performance. Because of the discretization inherent in an FEM implementation, the fields inside the layer are reproduced only approximately, and this is particularly true for a rapidly decaying field, (see Figure 11). To design a good absorber it is necessary to understand the impact of the sampling rate on the choice of a, fi and t, and our objective is to find the minimum number of sampling elements (or discrete layers) to achieve a specified IRI. It is anticipated that the errors introduced by the discretization will have a number of consequences. In particular, for a given number N of discrete layers and given t, increasing,i will ultimately lead to an increase in IRI because of the inability to model the increasing attenuation. An increase in c will likely produce a similar effect. To obtain some insight into the roles played by N, a, p and t, we now consider a simple FEM model of the layers. 3.1 Numerical Model A one-dimensional FEM code was used to examine the numerical performance of the absorbing layers. The computational domain was limited to the discretized layer structure shown in Figure 4, with the appropriate boundary conditions applied at the interface x = 0 and the PEC backing x = t. We consider first a homogeneous isotropic layer at normal incidence for which the theoretical reflection coefficients are given in (4). In spite of the fact that the magnitudes are the same for both polarizations, a polarization dependence shows 7

up in the FEM implementation. This is illustrated in 5. and Nwe note that as N increases, the FEM values of |R(0)j converge to the common theoretical value for bothi polarizations. 3.2 Dependence on a and P For a layer of constant thickness the theoretical value of IR(O)I is independent of a and polarization, but in the numerical implementation the behavior is much more complicated. Figure 7 shows IR(O)I plotted versus a and d3 for a layer of thickness t = 0.25Ao made up of 5 (=N) elements, where the darker tones indicate lower values. For small 3 the results are in close agreement with theory. As evident from the level lines, JR(O)J is almost independent of a and decreases exponentially with /3, leading to a linear decrease on a dB scale. For large /3, however, the behavior is quite different, and the most striking feature is the series of deep minima whose spacing in a increases with increasing a and decreasing /. These are numerical artifacts which are common to both polarizations and may depend on the particular numerical code employed. The minima for the two polarization are interlaced, and for H polarization the first minimum occurs at a = 0, p = 1.6. Their locations also depend on t and N. If N is fixed, the spatial sampling is inversely proportional to t. Decreasing t results in better sampling, pushing the minima to higher values of /3 and producing agreement with the theoretical values for larger /3 than before. Increasing t has the opposite effect. On the other hand, if t is fixed, increasing N improves the accuracy, and shifts the minima to higher /. Apart from the minima, the reflection coefficients for fixed 0 increase slightly with increasing a, and it is therefore sufficient to confine attention to the lower values of a. In Figure 6 the reflection coefficients are plotted as functions of P for the same 8

layer with a = 0 and a = 0.75. The curves correspond to vertical cuts through the patterns in figure 4, and we also show the theoretical value obtained from (4). \Ne observe that as /3 increases the reflection coefficients decrease initially at almost the same rate implied by (4), but beyond a certain point they begin to increase. The deep minimum at a = 0 and /3 = 1.6 in Figure 7 is clearly seen, but for design purposes it is logical to focus on the worst case, i.e. the polarization for which the reflection coefficient is larger. The upper curves in Figure 6 are almost identical and constitute this case. Since they correspond to two different values of a, either of them would suffice, but for reasons that will become evident later, we choose a = 0. 3.3 Dependence on 3, N and t We now seek a connection between the values of /, N and t for which JR(0)j is minimized. To this end, we first examine JR(0)j as a function of f3 and N for constant t, and the resulting plot is shown in Figure 8 for E polarization with t = 0.25Ao and a = 0 as before. For fixed /3 the reflection coefficient tends to its exact values as N increases. This is evident from the level curves and, as expected, the convergence is better for the smaller /. Consider now the behavior of |R(0)j for fixed N. As /3 increases from zero, the reflection coefficient decreases to a minimum and then increases. The location of the minima is indicated by the solid line. This is consistent with the behavior shown in Figure 9 and the upper curve is, in fact, just a vertical cut through Figure 9 at N = 5. The solid line in Figure 9 therefore gives the value of f at which \R(0)I is a minimum as a function of the number of elements. If the process is repeated for other layer thicknesses, it is found that for minimum IR(0)I the curve of /3t/Ao versus N is virtually the same for all thin layers. The observation that ft/Ao is a scalable parameter is an important conclusion of our study, 9

and by choosing a constant layer thickness we can produce a universal curve for the optimal choice of N and,3 in FEM simulations. Such a curve is shown in Figure 8 and can be interpreted as giving the value of 3t/Ao for a specified V to minimize the reflection coefficient {.R(0)|. For example, if t = 0.2Ao and N = 3, then 3 = 2.13. If a smaller value of fi is chosen, |R(0)| will be larger (see Figure 5), and this can be attributed to the fact that the field reflected from the metal backing has not been attenuated sufficiently. Iff3 is set to a value larger than 2.13, IR(0)I will still be larger because the chosen N is too small to reproduce the rapid field decay within the layer. So far we have considered only normal incidence, but for the anisotropic layer it is a simple matter to extend the results to all real angles of incidence k < 90~. As evident from the exponent in (4), the absorption at normal incidence is reproduced at any angle if the layer thickness is inversely proportional to cos/. This can be achieved by specifying the layer thickness t as a fraction 6 of the wavelength A. along the normal (or x axis) to the material interface. Since t = 8Ao/ cos = 8AX, t/At is now independent of t and the scalable parameter O3t/Ao (at normal incidence) becomes A3t/Ah. For the anisotropic layer, all the results obtained at normal incidence are made applicable for arbitrary ( by substituting A, for A0. For example, plotting the optimum,3t/A, as a function of N duplicates the curves shown in Figure 8. This notion can also be used to account for problems where the outer medium is not free space. In such cases, we have A, = Ao0/\' (at normal incidence) where cr is the relative permittivity of the outer medium. Although the scaling property of ft/AT has been established only for a = 0, it holds to a reasonable degree for small a $ 0, but as a increases, the f3t/AT versus N curves become increasingly dependent on a. The scalability also extends to the associated values of IR(0)|, and this enables us to provide a simple design prescription 10

for an absorbing layer. 3.4 Design Curves Since the quantities f3t/A\ and |R()I| are the same for layer thicknesses up to about 0.5Ar at least, design curves can be obtained by plotting JR()j\ and N versus 3t/Ar on the same figure as shown in Figure 10. To see how to use the figure, suppose that the desired reflection coefficient at normal incidence is -50 dB. At - = 0 we now have A, = oA. In Figure 10 we observe that the \R(0)j curve intersects the -50 dB line at ft/Ao - 0.58, and referring now to the N curve, the number of elements required is N = 10. The value of 3 can then be found by specifying either the element size or the layer thickness. Thus, for elements 0.025Ao thick, we have t = 0.25Ao and 3 = 2.32. By increasing N we can improve the performance up to the limit provided by the theoretical value of IR(0)I which has been included in Figure 10. A good approximation to the short-dashed curves in Figure 10 obtained by linear regression is t = -0.01061R| + 0.0433 (5) N = 0.147exp[7.353/t/Al] (6) where A, = A0/ cos X, RI is measured in dB and N is the smallest integer equal to or exceeding the right hand side of (10). These equations hold for <, 0~ < ~ < 900, for the anisotropic layer, and 4 = 0~ for the isotropic one. The design criterion provided above applies to specific angles of incidence. In the case where a specific absorption level is required over a range of angles of incidence the layer must be designed for the largest angle occurring. Doing so ensures that the absorption is superior at all smaller angles. 11

The performance can be improved by making the anisotropic material inhomogeneous, and to illustrate this we consider the case ^(x) = -j3(x/t)' for which the theoretical reflection coefficient is the same as before. The scalability is still preserved and the resulting curve is shown in Figure 10. The fact that the curve for the quadratically tapered layer lies below that of the homogeneous material confirms the improvement in performance, and we can now achieve a reflection coefficient of -50 dB by choosing 3t/A, - 0.64 corresponding to N = 9. Approximations to the long-dashed curves in Figure 10 are t = -0.01191RI + 0.0451 (7) Ax N = 0.298exp [5.263ft/A^] (8) where IRI and N are as before. Compared with the homogeneous material the decrease in the number of elements required becomes more pronounced as \RI is reduced. As previously mentioned, these equations hold for all real angles of incidence in the case of the anisotropic layer and for normal incidence if the layer is isotropic. 4 Applications As noted earlier, a PML is particularly attractive for terminating a finite element mesh in the simulation of microwave circuits. For these applications a PML has an advantage over a traditional ABC because it does not require a priori knowledge of the guided wave propagation constant. To demonstrate the applicability of the design equations in three dimensions, we consider a shielded microstrip line and a rectangular waveguide as well as some other examples. 12

4.1 Rectangular Waveguide section has dimensions 4.755 cm x 2.215cm and is chosen to propagate only the TE0o mode. It is excited by an electric probe at the left, and Figure 11 shows the mode field strength inside the waveguide which has been terminated by a perfectly matched uniaxial layer. As expected, the field decay inside the absorber is exponential and for 3 values less than unity the wave does not have sufficient decay to suppress reflections from the metal backing of this 5cm layer. Consequently, a VSWR of about 1.1 is observed for /3 = 0.5. However, as 3 is increased to unity, the VSWR is nearly 1.0 and the wave decay is precisely given by e-OktCOS9i = evrCcosGp where t is the wave travel distance measured from the absorber interface, P = 2/3t/A9 and here 6 = 44.5~. It is noted that when /3 is increased to larger values, the rapid decay is seen to cause unacceptable VSWR's. For optimum VSWR, /3 must be chosen (for a given absorber thickness) to provide the slowest decay without causing reflections from the absorber backing. That is, the lowest reflection may be achieved when the entire absorber width is used to reduce the wave amplitude before it reaches the absorber's backing. Not surprisingly (see equation (4)), for this example, the value of a does not play an important role in the performance of the absorber and this is demonstrated in Figure 13. As seen, setting always a = 3 gives the same performance as the case of a = 1 shown in Figure 12. Our tests also show that other choices of a give the similar absorber performance. 13

4.2 Microstrip line The performance of the perfectly matched uniaxial laver in absorbing the shielded microstrip line mode is illustrated in Figure 14 where the reflection coefficient is plotted as a function of 23t/Ag, where Ag = Ao/V/Te/f and oeff is the effective dielectric constant. In this case, the microstrip line is terminated with a 1.87 cm thick, 5-layered absorber and the line is extended up to 4 layers inside the absorber to avoid an electric contact with the metallic wall. Similarly to the waveguide, we again observe that an optimum: value exists and it was verified that in the absorber the wave exhibits the same attenuation behavior as shown in Figure 11. The reflection coefficient at the optimum,3 = 1 is now -42dB and if better performance is required, a thicker absorbing layer may be required. Again as in the case of the waveguide example, the value of a plays little role in the performance of the absorber and this is illustrated in Figure 15. However, of importance is the behavior of the reflection coefficient as a function of 23t/Ag. For the waveguide and microstrip examples, we observe that the absorption is maximized for approximately the same value of 23t/Ag (about 0.8). Thus these curves can be used for other applications as well, although it should be noted that the discretization rate plays an equally important role and this needs further investigation. 14

4.3 Band Eliminator Another guided structure that was considered is the stripline band eliminator shown in Figure 16. The band eliminator consists of a circular disk 1.64 cm in radius, placed on a substrate having a relative dielectric constant of 2.4. The disk is directly coupled to 50 Q striplines and to take advantage of symmetry, only half of the problem's domain is used in the FEM computation. The thickness of the absorber used in this calculation is 1.5cm, and in Figure 17 we show the insertion loss results as computed by the FEM using a mesh density of about 10 edges per wavelength. The values of a and 3 for the absorber were both 3 and, as seen, our calculations compare very well with measured data [20], thus, demonstrating the accuracy of the absorber termination. 15

4.4 Sphere Scattering As another example we consider the scattering by a metallic sphere illuminated by a plane wave (see Figure 3). The unknown quantity in the FEMNI implementation was the scattered electric field, and the excitation due to the incident plane wave was incorporated by enforcing the boundary condition Et = -Et on the surface of the sphere. The FEM mesh was terminated by cubical metal backed PML absorbers (with a = d 4.2) of thickness O.lAo and 0.2Ao, corresponding to roughly 1 and 2 PML layers, respectively, whose interface was placed only 0.05Ao away from the sphere's surface. The numerical results for both O.1Ao and 0.2Ao thick PML absorbers are presented in Figure 18 and we note that the sampling rate was approximately 14 elements per wavelength. As can be seen from Figure 18, the value of f3t/Ag closer to 0.8 results in better accuracy. The computed data are within IdB of the reference curve and this difference may be due to several reasons including surface modeling fidelity as well as mesh truncation accuracy. 16

4.5 Meander Line Another example is the meander line shown in Fig. 19. For the FENI simulation. the structure was placed in a rectangular cavity of size 5.8mm x 18.0mm x 3.1757mm. The cavity was tessellated using 29 x 150 x 5 edges and only 150 edges were used along the y-axis. The domain was terminated with a 10 layer PML, each layer being of thickness t = 0.12mm. The S11 results are shown in Figure 20 and are in good agreement with the measured data. 17

5 Conclusions A new method of mesh truncation based on an artificial absorbing lay er with anisotropic material properties was implemented for terminating the FENI mesh for propagation and scattering examples. It was demonstrated that reasonable performance can be achieved with relatively thin absorbing layers (a small fraction of a wavelength) placed close to the computational domain. Also, because there is no surface integral or ABC at the mesh truncation boundary, the method is easily implemented in FEM solutions and is well suited for parallel computations. 18

References [1] B. Engquist and A. Majda, "Absorbing boundary conditions for the numerical simulation of waves," Math. Comput. Vol. 31, pp. 629-651. 1977 [2] L.Halpern and L.N. Trefethen, "Wide-angle one-way wave equations," J. Acoust. Soc. Amer. Vol. 84, pp. 1397-1404, 1988 [3] W. Sun and C.A. Balanis, "Vecter one-way wave absorbing boundary conditions for FEM applications," IEEE Trans. Antennas Propagat. Vol. AP-42, pp. 872 -878, 1994 [4] Webb and Kanellopoulos, " Absorbing boundary conditions for the finite element solution of the vector wave equation," Microwave Opt. Tech. Lett. Vol. 2, pp. 370 -372, 1989 [5] A. Chatterjee and J.L.Volakis, "Conformal absorbing boundary conditions for the vector wave equation," Microwave Opt. Tech. Lett. Vol. 6, pp. 886-889,1993 [6] R.L. Higdon, "Absorbing boundary conditions for acoustic and elastic waves in stratified media," J. Comp. Phys. Vol. 101, pp. 386-418, 1992 [7] Z.P. Liao H.L. Wong, B.P. Yang and Y.F. Yuan, "A transmitting boundary for transient wave analysis," Sicentia Sinica Vol. 28, pp. 1063-1076, 1984 [8] M. Moghaddam and W.C. Chew, "Stabilizing Liao's absorbing boundary conditions using single-precision arithmetic," it IEEE AP-S Int. Symp., London, -Canada, pp. 430-433, 1991 [9] J.Fang, "ABCs applied to model wave propagation in Microwave integratedcircuits," IEEE Trans. MTT, Vol. 42, No. 8, pp. 1506-1513, Aug. 1994 19

[10] R. Luebbers and C. Penney. "Scattering from apertures in infinite ground planes using FDTD," IEEE Trans. Antennas Propagat. Vol. 42, pp. 731-735, May 1994 [11] B. Stupfel and R. Mittra, "Efficiency of numerical absorbing boundary conditions for finite element applications," URSI Radio Science Meeting, pp. 165, 1994 [12] K.K. Mei, R. Pous, Z. Chen, Y.W. Liu and M.D. Prouty, " Measured equation of Invariance: A new concept in field computations," IEEE Trans. Antennas Propagat. Vol. 42, pp. 320-328, March 1994 [13] T. Ozdemir and J.L.Volakis, "A comparative study of an absorbing boundary condition and an artificial absorber for truncating finite element meshes," Radio Science Vol. 29, No. 5, pp. 1255-1263, Sept. 1994 [14] J.-S. Wang and R. Mittra, "Finite element analysis of MMIC structures and electronic packages using absorbing boundary conditions," IEEE Trans. Microwave Theory and Techn., Vol. 42, pp. 441-449, March 1994. [15] Z.S. Sacks, D.M. Kingsland, R. Lee and J.F. Lee, "A perfectly matched anisotropic absorber for use as an absorbing boundary condition," IEEE Trans..Antennas and Propagation, December 1994. [16] Berenger, J.P. "A Perfectly Matched Layer for the Absorption of Electromagnetic Waves," J. Comp. Physics, 114: 185-200, 1994. [17] Chew, W.C. and Weedon, W.H. "A 3-D Perfectly Matched Medium from Modified Maxwell's Equations with Stretched Coordinates," Microwave and Optical Technology Letters, p.599-604. September, 1994. 20

[18] Katz, D.S, Thiele. E.T., and Taflove. A. "Validation and Extension to Three Dimensions of the Berenger PML Absorbing Boundary Condition for FD-TD Meshes," IEEE Microwave and Guided Wlave Letters. August, 1994, p.268-270. [19] C. Rappaport and L. Bahrmasel, "An absorbing boundary condition based on anechoic absorber for EM scattering computation," J. Electromagn. Waves Appl., Vol. 6, No. 12, pp. 1621-1634, Dec. 1992. [20] R. R. Bonetti and P. Tissi, "Analysis of Planar Disk Networks", IEEE MTT-26, July 1978, pp. 471-477. 21

5.1 Appendix: FEM Formulation To generate an appropriate linear system from equation (1), it is necessary to first discretize the computational volume into smaller volume elements. The field in each element can be approximated by introducing some linear or higher order expansion. Subsequently, by using different testing or weighting functions Wj, (j = 1, 2 3,.....N) a linear set of equations can be generated for the solution of the field expansion coefficients. The fields within each element are expressed as N, E= EeW e; (9) i=1 where Nv is the number of edges forming each element, W1e is the element basis or expansion functions. For our analysis, they have been chosen such that the expansion coefficents Ee represent the value of the field component along the ith edge of the eth element when evaluated at the location of the same edge. For the brick element shown in figure 21, it is convenient to first introduce the definitions. {Et} ={~i, y z} j. %. A T {W } = {ixNjyNY, JNz} (10) where j = 1,2,3,4. Thus, for example, Ee =,2, E = Oy3, W = xN2, W = yNy3, etc. The field expansions can then be written as 4 4 4 E= E Nxj(y,z);j; E E Ny(zx); z,=; (11) j=1 j=1 j=l N"j (x, y)X (12) 22

where (b - y')(c- ') y'(c-') (b- y'):' y'z' T = (c - - x') N '(a - x') = (c ) (13) ~ be a -r ~ - 3 ca be - be re (c - x')(b - y') e (b y) (a - x')Y' - x' ) y2 ~3 ~ -,,. (13) ab I z2 ab ab 3 - ab In these, (x', y', z') denote the local coordinates specifying a point within the eth element and from an examination of the expansion functions we observe that O44 represents an average of the field component Ex along the edge segment (1, 2). Likewise, 442 is associated with the E: component along the edge (3,4) and so on. The elements of the 12 x 12 [Ae] matrix defined in (10) can now be readily evaluated in a straightforward manner since Npj (p = x, y, z) are simple linear functions. To do so, for notational and programming convenience, we shall denote these elements as A(pj),(qk), where p = x, 1y, z; q = x, y, z; j = 1,2, 3, 4; and k = 1,2, 3, 4. We next define the matrices 4 2 2 1 abc 2 4 1 2 [Ko] = -36 36 2 1 4 2 1 2 2 4 1224 2 -2 1 -1 -2 2 -1 1 [K1] - 1 -1 2 -2 -1 1 -2 2 23

2 1 -2 -1 1 2 -1 -2 [K2] = -2 -1 2 1 -1 -2 1 2 2 1 -2 -1 -2 -1 2 1 [K3] = 1 2 -1 -2 -1 -2 1 2 where a, b and c are the side lengths of the brick as defined in Figure 21. On using this notation we have [Asj,xk] [Aj,yk] [Arj,zk] [A] = [Ayj,xk] [Ayj,yk] [Ayj,Zk] (14) [Azjrk] [Azj,yk] [Azjzk] If we define the material uniaxial tensors as follows; erxx 0 0 r r 0 eryy 0 (15) 0 0 rzz and; /rxx i 0 0 = 0 ryy 0 (I6) 0 0 Yirzz > 24

then. the element m~atrices [A4.rj]z = [Axiyk]= [Av3,x~k] [Av 3,Vk]= [Azj',xk] [.Az3,zk] wIll be on the form-, { ac [h a b 2 - I [IK3], [Ax 2,zk] b I-T 6/Irzz 6/dryy 3 { c ab bc[K2] + - b[Kit] -k2 6alrx 6cu,uzz O EryyjA b [AK3], [4 a Ki, yk] =-6iz bc[KI] ~ b [K2] -k'rzKO 6 a/uryy 6aprxx 0 (17-) (18) (19) (20) (21)I (22) 25

List of Figures 1 Plane wave incidence on an interface between two diagonally anisotropic half-spaces............................................... 28 2 A rectangular waveguide (a) and a microstrip line (b) truncated using the perfectly matched uniaxial absorbing layer.............. 29 3 Geometry of sphere scattering example.................. 29 4 Geometry of the metal-backed absorber layer (a) and its FEM implem entation (b)................................. 30 5 Numerical results for a homogeneous isotropic layer with t = 0.15Ao and b = -j2.5. The top four curves are for E pol., the bottom four for H pol.: ( ) exact, (- - -) N=3, ( ---) N=6 and (- -) N=12... 30 6 IR(0)| with t = 0.25Ao and N = 5: ( ) Exact, (- -) a = 0 E pol., (- - -)a = 0 H pol., (- — ) a =0.75 E pol. and (...) c= 0.75 H pol. 31 7 Plot of JR(0)I in dB for (a) an H polarized and (b) an E polarized wave incident on a homogeneous isotropic layer with t = 0.25Ao and N = 5. The solid curves are level lines.......................... 32 8 3t/Ao computed from the E pol. case with a = 0: ( ) t = 0.15Ao, (- -) t = 0.25Ao and (- -) t = 0.5A0...................... 33 9 IRI as a function of / and N for E pol. with t = 0.25Ao and a = 0. The dashed curves are level lines and the solid curve gives the location of the minimum IRI........................... 33 10 Absorber design curves. The straight lines give \RI in dB, and the curved ones give N: ( ) exact IR[, ( — -) homogeneous case and (- ) inhomogeneous case.......................... 34 26

11 Field values of the TE0o mode inside a waveguide terminated by a perfectly matched uniaxial layer. The absorber is 10 elements thick and each element was 0.5 cm which translates to about 13 samples per wavelength at 4.5 GHz.................................. 35 12 Reflection coefficient vs 23t/Ag (a = 1) for the perfectly matched uniaxial layer used to terminate the waveguide shown in Figure 2..... 36 13 Reflection coefficient vs 2/t/Ag, with a = /, for the perfectly matched uniaxial layer used to terminate the waveguide shown in Figure 2... 37 14 Reflection coefficient vs 2/t/Ag with a=l, for the shielded microstrip line terminated by the perfectly matched uniaxial layer.......... 38 15 Reflection coefficient vs 203t/Ag with a = /, for the shielded microstrip line terminated by the perfectly matched uniaxial layer.......... 39 16 Geometry of a band eliminator...................... 40 17 Insertion loss for the band eliminator shown in Figure 9........ 41 18 Bistatic scattering for a 0.4A diameter sphere for different values of 3t/A. The FEM mesh density was approximately 14 elements per A and was terminated by a cubical thin absorber whose inner and out diameter was 0.5A and 0.7A, respectively................ 42 19 Illustration of a meander line geometry used for comparison with measurement........................................... 43 20 Comparison of calculated and measured results for the meander line shown in Fig.19............................. 44 21 Local Geometry of the Rectangular Brick................45 27

Region 1 Reflected wave I- - - - - - - - - -- - Oi /...,. O _ ai 0 0 Incident wave 0 b1, 0 O0 C O Region 2 Transmttted wave ao o 0 b2 0 00 C2 Figure 1: Plane wave incidence on an interface between two diagonally anisotropic half-spaces. 28

Electric Probe \_ I'I ~=~= -=aijp =x = =y a-j d+t=40cm t= 5 cm:ross-section: 4.755x2.215 cm i —.w I dc (a). waveguide Microstripline Absorbing layer rotri l = Ar, L W::: L T.:::; jh-*K t1 d + t = 12.0 cm L = 2.38 cm h= 0.21 cm H = 1.06 cm w = 0.548 cm T -w tM i ti i (b). Microstrip Line, = 3.2 Figure 2: A rectangular waveguide (a) and a microstrip line (b) truncated using the perfectly matched uniaxial absorbing layer. -- Anisotropic Absorber....... Solution domain 0.1 A...... -.... ----: ~f::: Figure 3: Geometry of sphere scattering example. 29

y y.. N -- it 1 x (a) (b) Figure 4: Geometry of the metal-backed absorber layer (a) and its FEM implementation (b). 0.6 sin n 1.2 Figure 5: Numerical results for a homogeneous isotropic layer with t = 0.15Ao and b = -j2.5. The top four curves are for E pol., the bottom four for H pol.: (- ) exact, (- - ) N=3, (- - -) N=6 and (- -) N=12. 30

0 — 10 -20 -30 -40 -50 -60 0 1 2 3 4 5 Figure 6: JR(O)j with t = 0.25Ao and N = 5: ( ) Exact, ( ---) a- 0 E pol., ( ---) a= 0 H pol., ( — -) a 0.75 E pol. and (...)a=0.75 H pol. 31

0 -20 -30 IRI 2 [d8] -50 — 60 0 1 2 3 4 5 a (a) 0 -10 -20 -30 IRI 2 [dB] -40 -50 0-60 0o 1 2 3 4 5 a (b) Figure 7: Plot of IR(O)J in dB for (a) an H polarized and (b) an E polarized wave incident on a homogeneous isotropic layer with t = 0.25Ao and N = 5. The solid curves are level lines. 32

0.7 0.6 0.5 - 0.4 0.3 0.2 Figure 8: 3t/Ao computed -) t = 0.25Ao and ( ---) t 7-I P J N from the E pol. case with a = 0: ( ) t = 0.15Ao, (= 0.5Ao. 0 -10 -20 — 30 IRI [dBe -40 ~....4o -50 _ — 60 n -- 2 4 6 8 10 12 N 14 16 18 20 Figure 9: IRI as a function of f3 and N for E pol. with t = 0.25Ao and a = 0. The dashed curves are level lines and the solid curve gives the location of the minimum \R\. 33

7 -20 30 -30 '"25 -40 - 20 -m - ^.. \ / -50 2 15 -60,- '"";., 10 -70 5 -80, 0.3 0.4 0.5 0.6 0.7 0.8 Pt/Xx Figure 10: Absorber design curves. The straight lines give RI in dB, and the curved ones give N: ( — ) exact IRJ, (- - -) homogeneous case and (- -) inhomogeneous case. 34

64 66 68 70 72 74 76 78 80 Segments number along waveguide (alpha=1.0, f=4.5GHz) Figure 11: Field values of the TE1o mode inside a waveguide terminated by a perfectly matched uniaxial layer. The absorber is 10 elements thick and each element was 0.5 cm which translates to about 13 samples per wavelength at 4.5 GHz. 35

-5 -5 --- f-=4.0GHz --- f-4.5GHz 10 -- f=5.0GHz -5- - -20 - -35- - i //' -40 - -45 -50 0 2 4 6 8 10 12 14 2pt in Xg (a=1.0) Figure 12: Reflection coefficient vs 2Pt/Ag (a = 1) for the perfectly matched uniaxial layer used to terminate the waveguide shown in Figure 2. 36

O f I I II I I I -5 -__f4Q~ --- f=4.OGHz -- f —4.5GHz -10 f=5.0GHz..., ^" o -15 -m _25i i -,/ -30 // -35 -40 - / -45 --50, --- —--- I -I, I I.... I L 0 2 4 6 8 10 12 23t In kg (a=p) Figure 13: Reflection coefficient vs 2/3t/Ag, with a = w/, for the 2 uniaxial layer used to terminate the waveguide shown in Figure 2. 14 perfectly matched 37

f 2 2.5 2tt in kg (a=1.0) Figure 14: Reflection coefficient vs 23t/A9 with a=l, for the shielded microstrip line terminated by the perfectly matched uniaxial layer. 38

0. -15 -20 0 0.5 1 1.5 2 2.5 3 3.5 4 2't In g (a-p) Figure 15: Reflection coefficient vs 2,3t/Ag with a -, for the shielded microstrip line terminated by the perfectly matched uniaxial layer. — 30 -35 -40 -45 -50-'- - 0 0.5 1 1.5 2 2.5 3 3.5 4 213t in Xg ((a=3) line terminated by the perfectly matched uniaxial layer. 39

6.0 cm output port Top View 5 mm 6.0 cm - - input port 6.0 cm "I — a Cross Section View i I -------— __ -----------:6.4 mm 5mm 5.36.81.30.6 Figure 16: Geometry of a band eliminator. 40

Band Eliminator 0.0 -10.0 -20.0 -30.0 V 'a Cl) C,) 0 -J Q,) c c -40.0 -50.0 L 2.0 2.5 3.0 3.5 Frecuencv (GHz) 4.0 Figure 17: Insertion loss for the band eliminator shown in Figure 9. 41

I 0.0.................... — ---------------........ cO -10.0 o —o Mie Series A A PML, t = 0.2 wavelength + --- PML, t = 0.1 wavelength -20.0 '....... 0.0 30.0 60.0 90.0 120.0 150.0 180.0 Figure 18: Bistatic scattering for a 0.4A diameter sphere for different values of 3t/A. The FEM mesh density was approximately 14 elements per A and was terminated by a cubical thin absorber whose inner and out diameter was 0.5A and 0.7A, respectively. 42

I II 1.525 mm.61 m ic x Figure 19: Illustration of a meander line geometry used for comparison with measurement. 43

frsqecy (GHz) Figure 20: Comparison of calculated and measured results for the meander line shown in Fig.19. 44

12 z i y / x Figure 21: Local Geometry of the Rectangular Brick. 45