034891-1-T Scattering From Plates Containing Small Features Using the Adaptive Integral Methods (AIM) Sunil S. Bindiganavale, John L. Volakis, and Hristos Anastassiu May 1997 34891-1-T = RL-2480

Scattering from Plates Containing Small Features Using the Adaptive Integral Method (AIM) Sunil S. Bindiganavale, John L. Volakis and Hristos Anastassiu Radiation Laboratory Department of Electrical Engineering and Computer Science The University of Michigan Ann Arbor, MI 48109-2122 Abstract Fast integral equation algorithms provide for a reduction of the computational complexity of moment method integral equation solutions. The adaptive integral method (AIM) is one such algorithm which has been demonstrated to reduce memory down to 0(N1"5) and computational complexity down to O(N15 log N) for surface problems consisting of N unknowns. The reduction is even further for planar bodies and in this paper we demonstrate the suitability of AIM to model planar scatterers with intricate geometrical details. It is shown that AIM is extremely accurate while saving a significant amount of memory, especially for bodies requiring high tessellation due to the inclusion of very small features. AIM is particularly attractive for modeling planar three-dimensional scatterers even on workstation platforms, reducing the need for elaborate parallelization and domain decomposition procedures. The capability of AIM to model small and intricate perturbations on otherwise smooth bodies and to predict near fields accurately (as reported in this paper), even for bodies which are not electrically large, makes it attractive for many scattering and radiation computations. 1

1 Introduction Fast integral equation methods were introduced in the early 1990s and have been shown effective in accelerating the computation of matrix-vector products in iterative solvers. The Adaptive Integral Method (AIM) [1] and the Fast Multipole Method (FMM) [2] belong to this class of fast integral solution techniques. Both AIM and FMM reduce the solution time and memory requirement of the moment method (MM) solutions and their initial applications focussed on electromagnetic scattering from large conducting bodies. More recently, they have also been used successfully in hybrid methods [3],[4] to evaluate the scattering from composite structures. FMM achieves its CPU reduction by grouping the farzone unknowns and interacting their weighted contributions. In the case of AIM, the CPU reduction is achieved by mapping the original MM discretization onto a rectangular grid and exploiting the Toeplitz property of the Green's function on this grid. That is, the Fast Fourier Transform (FFT) is invoked to compute the matrix-vector products in the iterative solver. For an arbitrary three dimensional body, a three dimensional FFT is required and as can be understood, this calculation is very time consuming. For planar scatterers the dimensionality of the FFT is reduced by one, thereby significantly accelerating the solution. In this paper, we examine the benefits of AIM when the body is not electrically large, but is highly tessellated owing to its intricate construction, thus leading to a large unknown count. We show that significant savings in CPU and memory can be achieved by AIM and examine its accuracy for near field and far field computations. 2

2 AIM for Planar Scatterers The AIM algorithm for arbitrary three dimensional bodies has been given in [1]. In this section, we describe its specialization to planar scatterers, giving only the essential details reqluired for its implementation. Consider a planar, arbitrary, conducting body whose surface S is illuminated by an incident plane wave E. The boundary condition enforced on S is (E' + Es) t =0 (1) where the scattered field Es is given by Es = -jA - V, (2) in which e-jkR A(r)= / Js(r') RdS' (3) is the vector potential and f s e-Jk O(r) = 4rwe j Vs J(r ) R dS (4) is the scalar potential. Substituting (2)- (4) into (1) gives the necessary integral equation which must be discretized for the solution of J,. To do so, the equivalent electric current density J. is expanded using linear basis functions [5] as J,(r') = E Infn(r') (5) N where In are the unknown coefficients. Application of Galerkin's technique leads to the linear system [Z]{I} = {V} (6) 3

with [Z] being the elements interaction matrix, whereas {I} is the vector of the unknown coefficients and {V/} is the excitation vector. The matrix [Z] is fully populated, demanding 0(N2) storage, and each [Z]{I} matrix-vector product requires 0({N2) multiplications. Fast algorithms such as FMM and AIM are used to reduce the operation count from NV2 down to NA, where a < 1.5. Both algorithms work on approximating the far zone interactions. In the case of AIM, the CPU reduction is achieved by first splitting the matrix as [Z] = [Z ]n + [Z '] (7) based on a threshold distance referred to as the near-zone radius. The matrix [Z,"ar] contains the interactions between elements separated less than the threshold distance, whereas [Zfar] contains the remaining interactions. The elements of [Znear] are evaluated with the exact MM while those of [Zfar] and the product [Zfar]{I} are evaluated in an approximate manner as prescribed by the AIM procedure [1]. Application of AIM requires that the whole geometry be enclosed in a regular rectangular grid. Basically, the fields of each interior edge is re-expressed using a new expansion based on delta sources located at the nodes of the uniform AIM grid as depicted in Figure 1. For the mth edge, this new expansion has the form M2 m = E ((X - Xmq)5(y - Ymq)[Axnq + A\ny] (8) q=1 where rmq are the position vectors of M2 points on the square surrounding the center of the edge and 6(x) is the usual Dirac delta function. The coefficients AX,! are suitably chosen so that the new expansion is equivalent to the original representation using triangular elements. 4

A similar expansion is used for the divergence of the basis functions M2 = 6 ( - Xm,)6(y - ymq)An (9) q=1 To find a relation between the AxY and In coefficients, we equate moments of the two expansions up to order M. Specifically, we set Mqlq2 = Fq,q2 (10) where Mqq2 - / )qI (y - y)q2dxdy for 0 < ql, q < M = 5(Xmq - X)1(yq- y)2(y[Axq + AM q] with q = ql + q2 (11) q=l F00 f y00 Fq = / f - (x-x)q (y-a )q2dxdy (12) ql-2 o -- -00 J - 00 Similarly, by equating moments of V, Ja with the new expansion (9), we establish a relation between Ad and I,. That is, we set D = HM (13) Dqlq2 = ql,q2 (13) where r0 ro M2 DMq = / f(x - x~)q(y - y)q2dxdy = "(xmq - xa)ql(ymq - y)q2Adq (14) qo ql rco roo Hqlq2 = J J fm(X - - xa)ql(y - ya)2dxdy (15) (10) and (13) give three M2 x M2 systems yielding the equivalence coefficients as the solution. This process is depicted pictorially in Figure 1. 5

Y ZL AIM grid -* * *, * * I Original NMM discretization (section) * e Y Y Am2 A ml * m9A mm7 I AIM' Representation for edge "m" (Note: Each of the * b original basis functions is * ^ _ defined on triangle pairs) 2 MN =3.. AIM coefficients = MN =9 Original MM discretization Figure 1: The process of transformation from the original MM grid onto the AIM grid Were we to use the equivalent expansions to represent the currents everywhere, the resulting impedance matrix will be of the form 3 [Z]AM= [A[G] [A] T i=l (16) In this, [A]i are the sparse matrices containing the coefficients of the expansion (11 ) and (9) whereas [G] is the Toeplitz matrix whose elements are the free space Green's function evaluated at the grid points. It has been shown in [1] that [Ztfj-] is not of sufficient accuracy ~AIMJ for modeling the interactions between the nearby current elements. To take advantage of the Toeplitz structure of [G] and sparsity of [A] we can still use [Z1Jo1a] to represent the far element interactions. However, we will retain the exact interaction matrix elements for the near element interactions. That is, we rewrite [ZTtoj] as _I AIMJ aLs [ztotal = [ near lJAIM = JAAIM + [z JAIM 6 (17)

Comparing this to (7) and setting [Z]'far [Z]faM we can rewrite the original [Z] matrix as [z] ([r - [Z]ea) + [Z]tota (18) or [Z] [S] + Z[A]i[G][A]T (19) i=1 where [S] = [Z] - [Z],a5 is a sparse matrix corresponding to the difference between the near field interactions computed by MM and AIM. The Toeplitz property of the Green's function, defined on the regular grid, enables use of the FFT to accelerate the computation of the matrix-vector product. The sequence of operations involved in the construction of the coefficient and Green's function matrices are indicated in Figure 2(a); those for the matrix-vector product execution are outlined in Figure 2(b). In the computation of the matrix-vector product, the initial step of transforming the currents from the original MM grid onto the uniform AIM grid is comparable to the grouping operation of the FMM. While the FMM relies on grouping to reduce the number of scattering centers, the sequence of operations in AIM can be interpreted as a realignment of scattering centers onto a regular grid. Although, this process does not reduce the number of scattering centers, the regularity of their location enables use of the FFT for fast computation of matrix-vector products. 3 Results When examining the merits of a fast integral algorithm such as AIM, of importance is the memory and CPU requirements, both contrasted to the delivered accuracy. Although approximate analytical expressions have been derived in [1] for some of these parameters, 7

Construct cartesian grid which extends over the object Calculate moments of the original basis functions (RWGs).....j d.??,,............................................................................................................................................. Solve for equivalence coefficients................................................................................................. I Calculate Toeplitz Green's function over the cartesian grid................................................................. Compute the near-field matrix (interactions between elements closer than threshold) (a) Compute I = A I ( Performs transformation of current from MMn grid to AIM grid) Compute I= ' (Compute the Fourier transform of the discrete current distribution)!............................................................................................... Compute V = (G r Far-field Fourier Transform - ( is the Green's function),.......................................,,. I Compute V= ' V I (Inverse transform to AIM grid) l.............................................................................................. Compute V AV (Far field back to MM representation)................................................................................................... Com ut ear near Compute V = Z I (Near field with exact MNM)..................................................................................................... Add to form total field far near V(+ V (b) Figure 2: (a) Matrix build operations and (b) Matrix vector product computation in AIM 8

these refer to implementations involving cubical grids and the three-dimensional FFT. Our goal in this paper is to assess the accuracy of AIM in treating small details within an aperture/surface and to provide the reader with quantitative measures on the performance of AIM when implemented with the two dimensional FFT. The near-zone radius or threshold distance has a dramatic impact on the CPU requirements since it controls the non-zero element population of the system matrix. In the case of AIM, because of the inherent mapping to a uniform grid, we are highly interested in examining its suitability to model small and fine details embedded in much larger scale structures. The calculations for the plate configurations given next are intended to address this issue by examining the method's performance for a number of representative and practical situations. All of the included results were generated using single precision arithmetic on an HP9000/C-110 workstation with a rated peak speed of 47 Mflops (the level 4 optimization option was also used). In all cases, a third order (M=3) multipole expansion was used with a grid spacing of 0.05A. Figures 3-7 depict the 00 and + polarization radar cross section patterns (-0O cut) as calculated by AIM for the different threshold distances indicated on the figures. The first circular plate has no holes and was used to validate the method. From the pattern comparisons, it is clear that AIM recovers the exact result very well. As given in Table I and 2, AIM achieves this with at least a factor of five less memory than the traditional MM, even though the geometries are still rather small to demonstrate the full impact of AIM. Also, Table 2 shows that a near zone radius of 0.3A is sufficient to maintain good accuracy (below one dB in RMS error [6]). The advantage of AIM is more pronounced when gaps are inserted into the plate's surfaces 9

and this is the primary reason that one may prefer AIM over other fast integral methods for planar structures. As depicted in Figures 4 and.5, AIM maintains its accuracy for the same threshold criterion even though the gaps/slots have a dominant effect on the RCS pattern as shown in Figure 4. In the case of narrow slots (or thin ridges in the plates) of width 0.03A, the memory requirements of the traditional MM increase quickly due to the higher element density. For the geometry in Figure 6, AIM yields memory saving of 79% and the CPU time is reduced by a factor of 12 while retaining the monostatic pattern accuracy to within a tenth of a dB. This is achieved by using a uniform AIM grid density of 20 points per linear wavelength even though the cell density of the original plate mesh is much greater due to the narrow slot. One may assume that this change in grid density will affect the near zone field. However, our observations indicate that the surface current is equally accurate. For the configuration in Figure 6 the average current density error is 7.3% for a threshold distance of 0.2A and 6% for a threshold distance of 0.4A. The currents for the geometry in Figure 6 along the center narrow strip are plotted and compared in Figure 8. These results demonstrate the important feature that the near zone threshold criterion is not affected by the specific geometrical details, leading to tremendous memory savings. Moreover, the accuracy of the results provide a convincing argument that AIM can efficiently handle highly irregular and resonant (i.e. antenna) geometries as well as smooth scatterers. At the same time, the convergence rate of the AIM system is unaffected indicating that the system condition is unchanged. This is of critical importance for fast iterative solutions, since an increase in the iteration count would annul the faster computation of the mnatrix-vector product. 10

Figure 7 shows the monostatic RCS pattern for a grating structure which acts as a "polarization filter". The thin ridges in the grating cause a strong specular return for the Ac polarization (almost 10 dB above the return in the absence of the gratings) as is evident from the results in Figure 7(d). Of importance is that the MM triangular mesh in Figure 7 required a cell size of 0.02A per linear dimension because of the narrow grating. However, the overlaid rectangular AIM grid could be selected to have a much coarser discretization. More specifically, we chose grid spacings of 0.05A and 0.1A for the AIM grid and, thus, computational requirements of AIM were much lower. For the 0.1A grid spacing the solution time was reduced from 2.75 minutes down to only 12 secs at the expense of some accuracy (fraction of a dB). To further increase in accuracy, we employed a 0.05A grid spacing and as shown in Figure 7(b) the AIM curve is now indistinguishable from the reference MM result (within 0.1 dB). From Tables 1 and 2, the AIM computational and memory requirements are 8 times and 9 times less, respectively, without loss of accuracy. This is a significant observation and we have found that both the convergence rate and condition of the AIM system remains essentially unchanged from the original moment method system. 4 Summary The performance of AIM is much improved when applied to scattering from flat complex scatterers. A memory reduction of 5 to 10 times over traditional MM was observed without compromise in accuracy when using a threshold radius of 0.2A. This CPU reduction is achieved without resorting to parallelization or optimization techniques (as is known AIM 11

is particularly amenable to such improvements). More importantly, the AIM algorithm is capable of modeling very small details in large bodies with a high degree of accuracy, while simultaneously saving considerable memory. This is of importance when modeling broadband antennas (spirals or log-periodics) and gratings which are both large in overall size but can contain features as small as A/100 in size. Discretization Geometry Facets Edges Unknowns MM memory (MB) MM solution time 00 pol (0 == 0~ inc.) Figure 3 586 908 850 5.51 32 sees Figure 4 554 890 772 4.54 29 secs Figure 5 1130 1806 1584 19.14 4 mins 50 secs Figure 6 1036 1667 1441 15.84 4 mins Figure 7 1038 1957 1157 10.21 2 mins 45 secs Table 1: Solution CPU time and memory requirement of the moment method References [1] E. Bleszynski, M/. Bleszynski, and T. Jaroszewicz. AIM: Adaptive integral method for solving large-scale electromagnetic scattering and radiation problems. Radio Science, 31(5):1225-1251, 1996. 12

AIM Data Geometry Threshold Non-Zeros Memory Solution time RMS Error(dB) (A) in Near Z (MB) 00 pol (O = 0~ inc.) 00 pol < pol 0.3 59928 0.68 23 sees 0.1718 0.0755 Figure 3 0.4 100182 1.14 25 sees 0.1490 0.0693 0.7 257390 2.94 28 sees 0.0728 0.0490 Figure 4 0.4 79030 0.9 21 sees 0.0728 0.0583 0.6 157994 1.8 27 sees 0.0721 0.0520 Figure 5 0.7 283774 3.24 3 mins 32 sees 0.8017 0.5185 Figure 6 0.2 296250 3.39 20 sees 0.1063 0.0949 0.4 649556 7.43 31 sees 0.0548 0.0632 Figure 7 0.2 120220 1.37 18 sees 0.0469 0.0469 Table 2: Solution CPU time, memory requirement and RMS error of AIM (all entries in this table were computed with an AIM grid spacing of 0.05A) [2] V. Rokhlin. Rapid solution of integral equations of scattering theory in two dimensions. Journal of Computational Physics, 86(2):414-439, 1990. [3] S. S. Bindiganavale and J. L. Volakis. A hybrid FEM-FMM technique for electromagnetic scattering. IEEE Transactions on Antennas and Propagation, 45(1):180-181, 1997. [4] N. Lu and J. M. Jin. Application of the fast multipole method to finite element-boundary integral solution of scattering problems. IEEE Transactions on Antennas and Propaga

tion, 44(6):781 —786, 1996. [5] S.M. Rao, D.R. Wilton, and A.W. Glisson. Electromagnetic scattering by surfaces of arbitrary shape. IEEE Transactions on Antennas and Propagation, 30(3):409-418, 1982. [6] S. S. Bindiganavale and J. L. Volakis. Guidelines for using the fast multipole method to calculate the RCS of large objects. Microwave and Optical Technology Letters, 11(4):190 -194, 1996. 14

Plate orientation - reference YAL z x Polarization 00 Polarization Z --------- Polarization reference y oQ Polarization X 25,15 CZ -15 25 ' 5;J.15 -25.............. r-... I.... w. I.. I - I — ' 1 I... ' I..., MM El AIM - Threshold = O.3 A AIM - Threshold = 0.43 AI * - Threshold = 0.7/% 25 20 - 3 15 15 cn 10 Q Cu -5 -in....... I.. l-................ 0 15 30 45 60 Observation Angle e (deg.) 75 90 0 15 30 45 60 Observation Angle e (deg.) 75 90 Figure 3: Monostatic RCS for a circular plate of diameter 2A; Comparison of the standard MM & AIM 15

Plate orientation reference Y z j 25 j 15 -2, 5 -15 -25 1.6X (a) 00 Polarization 25 x 15 -5 -15 -25 15 30 45 60 7' Observation Angle 0 (deg.) (b) ~~ Polarization '... I I I.. I ' " -- MM a AIM - Threshold = 0.4i. AIM - Threshold = 0.6,6 5 90 i 0 15 30 45 60 75 90 0 15 30 45 60 75 90 Observation Angle 8 (deg.) z Observation Angle e (deg.) (c) i, (d) Polarization reference y X Figure 4: Monostatic RCS for a circular plate of diameter 2A with three slots computed with standard MM & AIM (a) Geometry (b) Effect of the slots on the RCS (c) 00 polarization backscatter RCS for the plate with slots (d) Hi polarization backscatter RCS for the plate with slots 16

Plate orientation reference Y X Z X 3.2 I 4 X 1.9 X 00 Polarization Z Polarization reference y P0 Polarization X 30 25 320 r, ~15 y mo -5 -10 0 15 30 45 60 75 90 0 Observation Angle 0 (deg.) 15 30 45 60 Observation Angle 0 (deg.) 75 90 Figure 5: Monostatic RCS for a circular plate with standard MM & AIM of diameter 3A with three holes computed 17

Plate orientation reference Z X Z A, / Y0, -10 7c F v r0 n -20 L I 1 - UU t0;m A c -30:Q MM A LM - Threshold = 0.2. A LN, - Threshold = 0.4, I I I I 0 15 30 45 60 Observation Angle e (deg.) 75 90 0 15 30 45 60 75 90 Observation Angle e (deg.) Figure 6: Monostatic RCS for a circular plate of diameter 1A sampled at 0.03 A (smaller than the AIM grid spacing) due to the narrow center ridge 18

..l._mMc^^ - WoooooeMuW^m 0 O.04A t 0.05 k 0.02) lA 1 k 0.1). QRtiyjwcxxxx Plate orientation reference YAL z x oZ X Z ~~ ---~ 'A0, 0.04) 0.02), (b) ce / Y X (~~ Polarization XL ^ 1a) (a) Polarization referen 00 Polarization.......... I..... I.......... -o -s -10.-15 -20 c' -25 -30 - - - AM - Grid spacing = 0.05k (For Figoar ' ' \\ _ - - AI1 - Grid spacing = 0.1 (ForFig (a) \ -- MM (For Fig (b)) m la "c C, a t. C. 2 v 1i ~. a e Lj 0:0 15 10 1 5 i1 CS -5 -10 21C,....,I.. I..... I. I..... I..... _ MIIMM (For Fig (a)) - - - -AIMI - Grid spacing = 0.05. (For Fig (a)e - - - AI - Grid spacing = 0.1k (For Fig (a)) --- MM (ForFig(b)) '\ A -N d\./ \.</ \ /... I... I... I.... I..... 0 15 30 45 60 Observation Angle 9 (deg.) (c).-X 75 90 0 15 30 45 60 Observation Angle 0 (deg.) (d) 75 90 Figure 7: (a) Geometry and mesh of the grated plate (b) Geometry and mesh of the "Groove" plate without gratings (c) 90 polarization backscatter RCS computed by AIM and MM (d) dx polarization backscatter RCS computed by AIM and MM 19

Plate orientation reference YA z x 0.6 k Currents at the centroid of the shaded triangular patches on the narrow ridge -- MM * AIM - Threshold = 0.4 it u 0 *a I E wo 0 o c.S,s 0i 4.5 4.3 4.1 3.9 3.7 3.5 33 3.1 2.9 2.7 2.5 x 1-4 K 10 ---— h ------— ~ ----~ --- —-- ---- — T — ---- I —r- I --- — I I — -— r ----r --- —-- I — - - - I - - - I I I - - - I — - - - I- - I I- - - - -- - - - - - - - - --- - I I I Il I I I..... -0.2 -0.1 0 0.1 y in wavelengths 0.2 0.3 Figure 8: Electric currents (Solution coefficients) on the narrow ridge for the geometry of Figure 6 20