032578-1-T ISOMOM:A SPARSE MOMENT METHOD TECHNIQUE FOR A WIDE CLASS OF 2D AND 3D SCATTERING PROBLEMS J. L. Volakis S. Bindiganavale A. Chatterjee Wright Wright OHIO Laboratories/AARA Patterson AFB 45433-7001 Veda, Inc Suite 200 5200 Springfield Pike Dayton, OH 45431 January 1995 32578-1-T = RL-2462

A sparse method of moments technique for a wide class of 2D and 3D scattering problems John L. Volakis, Sunil S. Bindiganavale and Arindam Chatterjee Radiation Laboratory Department of Electrical Engineering and Computer Science The University of Michigan Ann Arbor, MI 48109-2122 Abstract This report introduces a new approach for generating and solving sparse moment method matrices. The technique is based on an iterative refinement procedure, which guarantees convergence to the exact solution, provided certain sparsity criteria are satisfied. The technique has been applied to scattering by both 2D and 3D geometries. The full moment method matrix was made sparse by zeroing out as a first approximation, the matrix entries corresponding to elements lying physically far away from each other. As expected, it was demonstrated that this technique works very well for physically large bodies. Scattering parameters, computed with different degrees of matrix sparsity and at different stages of the refinement procedure, are reported. 1

1 Introduction Integral equation methods and in particular the moment method using subdomain basis functions for solving open domain problems were introduced first in the 1960s [1]. In the ensuing years the moment method achieved spectacular success and was applied to model a wide variety of complicated electromagnetic phenomena. Among its chief attractions were the exactness of the solution technique and the ability to model complicated geometries with ease. However, these advantages came for a price. The price to be paid for the exactness of the solution was the full matrix systems that were a necessary part of the technique. Since full matrices need a storage of O(N2), MoM suffers severe problems of scalability, especially for 3D problems and for large 2D problems. Recently, researchers have been trying to develop IE techniques which make the matrix system sparse by using special basis functions [2] or clumping the matrix elements in the far-field [3]. The first method has been applied solely to 2D problems and its extension to 3D problems may not be possible. The second method, also known as the Fast Multipole Method (FMM), is yet to gain widespread acceptance as a valid method for solving 2D & 3D problems. In this report, the problem of making the full matrix sparse is attacked in a different manner. It was observed that the interaction between elements lying physically far away from each other is small. If they are sufficiently far away, the interaction should be almost negligible, i.e., the matrix entry can be zeroed out as a first approximation. Thus, the full matrix reduces to a sparse matrix, although the number of elements to be zeroed out depends on the size of the body and accuracy of the desired results. For solving large body scattering, this technique should work even better since the percentage of sparsity increases as a function of the physical size of the body. We have examined the idea both for 2D as well as 3D geometries. The performance of the technique for 2D has been worse since the magnitude of far-off elements taper off as 1/./F opposed to the 1/r taper of 3D. However, the problem was corrected by using the Physical Optics (PO) currents as the first approximation and iteratively refining the solution. Reasonably accurate results were obtained with a 20% full system. In 3D, the method performed comparatively better; introduction of the PO currents, however, produced no significant improvement. An important component of this technique (to be referred to as Isomom) is the capability to improve upon the approximation in an iterative manner without using additional storage resources. One or two iterations beyond the first yields the correct result provided certain sparsity criteria are satisfied. 2 Two dimensional problems The governing integral equations for computing two dimensional scattering from metallic bodies are the well-known MFIE and EFIE. The TM (or E-polarization) integral equation is given by [4] k0Z0 jX Jz(J)H2) (koIT - dl - E'(p) (1) 2

Expansion of the current using pulse basis functions and point matching to obtain a set of linear equations, was employed in the moment method solution of these integral equations, resulting in [Z]{J} = {V} (2) where An { (2)(koRn) (m ~ n) Zmn a- [1- O (ln (ka) -1 )] (m = n) (3) with Aim A In 2 Rmn = [m- Xn + - COS m - 2 COS On 2 2sn /lm Alsin 2 + Ym - n + in m sin (4) and V - ejk~[(xm+ Alm cos m) cosqo+(ym+am-sinAm)sinqo] (5) Alm and Ain indicate the segment lengths. The full system in (2) is inverted to obtain the unknown current values for different excitation vectors. The TE (or H-polarization) integral equation is given by 1 & Jl() + 1 (')- 2 (k - )1 p-4 =) dl' =-H ( y) (6) 2 4 On' Discretization of (6) using pulse-basis point-matching results in a system of linear equations as in (2) where jko oaA.I Al )H(2) Zmn = S(m n) + j (n )I )(koRn dl (7) 2 4 (koR~) d (7 with I 1 -mn = Rp (8) Rmn= (xm - n) + (lm cOm - 21' cos On)] + Y [(Ym -y n)+ 2(Alm sinOm 21 sin n)] (9) Vm - ejko[(xm+ lm COS0m)COSqo+(ym+Alm-sin0m)sinqo] (1 Vm e, 2 2 1 (11) 3

Owing to the operation of the derivatives on the Green's function, the discretization of the integral equation for H-polarization (TE incidence) results in a more diagonally dominant matrix than that for E-polarization (TM incidence). A single row of the impedance matrix for the cylinder whose cross-section is depicted in Figure 1 is shown in Figure 2. This row, being representative of the character of the matrices encountered in the pulse-basis pointmatching solutions of 2D scattering problems, hints that any attempt of making full matrices sparse should address both polarizations separately, owing to their innate different degrees of sparsity. y + Segment #50 - - - - - - x Segment# 275 Segment # 1 L Segment # 325 15X Figure 1: Test body #1 - Rectangular cylinder with circular end- caps In this report we first consider the H-polarization, with its strongly diagonal impedance matrix, and next the E-polarization and its weakly diagonal matrix. 2.1 H-polarization Owing to the strongly diagonal nature of the H-pol impedance matrix, the H-pol currents can be solved for by directly applying a tolerancing criterion. This implies that the matrix entries corresponding to the interaction between elements lying physically far away from each other is zeroed out, once they drop below a certain threshold value. Tolerancing criteria were applied using the distance criterion, where elements which lay physically farther than a specified "window" distance were discarded, and the magnitude criterion, where elements of each row with magnitudes less than a specified fraction of the self-cell were discarded. For 2D applications, it was found that in most cases the above criteria were equivalent, while in 3D the magnitude criterion was found to be more favorable. Hence the essence of this technique is, as a first iteration, to solve the sparse system [ZS]{J(l)} = {V} (12) where {J(l)} is a first approximation of the current, {V} is the excitation vector, and the original impedance matrix [Z] is split as [Z] = [ZS] + [ZR] (13) 4

I. - -, - - H-pol -- E-pol 1I I 0.5 i - -,h ---hh- - -- L — --. 0 100 200 300 400 500 Segment number n (a) 600 (b) (c) Figure 2: (a) Impedance matrix of test body #1 (row=150) (b) and (c) Regions of interest with [ZS] being a predominantly diagonal sparse matrix and [ZR] being the residual matrix. The strength of this new technique is an iterative refinement procedure which can be employed to improve on the first approximation. The iterative refinement technique performs the following mathematical operation, corresponding to a second iteration [ZS]{J() = {V} - [ZR]{J(l)} (14) In general, the mth iteration is described by [ZS]{J(m)} = {V} - [ZR]{J(m-) (15) The convergence of this iteration process can be investigated by considering (14) and multiplying each side by [ZS]-1 to get {J(l)} = [ZS]- ([I]- [ZR][Zs]-1) {V} (16) 5

After a subsequent iteration {j(2)} = [ZS]- ([I]- [ZR][ZS]-1 + ([ZR][ZS]-1)2) {} and at the mth iteration {J(m)} = [ZS]-{[I] - [ZR][ZS]- + ([ S) ([)2 [ZR][ZS]-1)3 +... ([ZR][ZS] -l)m}{V} (17) (18) This geometric series will converge if every eigenvalue of [ZR][ZS]-l has a magnitude less than one [5]. The result of this convergence after an infinite number of iterations would be {j(oo)} = [ZS]- ([I] + [ZR][ZS]-1)1 {V} (19) and would imply convergence to the right answer if [Z]- = [Z]- ([I] + [ZR][Z]-1)-1 (20) which is indeed the case, and can be proved by multiplying each side of (20) by ([Zs] + [ZR]). Application of (12) to the test body shown in Figure 1 results in the backscatter RCS computation shown in Figure 3. The sparse matrices were generated by employing a thresh m,C$ cn u 04 CO U CQ m Ad 40 35 30 25 20 15 10 5 0 -5 0 15 30 45 60 75 90 Observation angle (deg.) Figure 3: H-pol Backscatter RCS of a 15A by 2A rectangular cylinder with circular end-caps olding operation. Mathematically, an element in the ith row is discarded if IA(i,j)l < {min(lA(i,j)l,Vj)+ T.F* IA(i, i)l} (21) 6

40 35 - I m-...... 11% Full 30 -25 7 -' 20 A X (1 Y f Q 15 U c-6 10 - 5 0 15 30 45 60 75 90 Observation angle (deg.) Figure 4: H-pol backscatter RCS of a 30A x 2A rectangular cylinder with circular end-caps where T.F. is a threshold factor and A(i, i) is the diagonal element. The threshold factor was set to 0.35% to generate the 23% full matrix and when set to 0.5% it generated the 19% full matrix. The backscatter computation indicates that the maximum error in RCS prediction was less than 3dB with a 19% full matrix and this was at angles close to grazing. As expected, with the 23% full matrix, the maximum error was even less (2dB or so) at angles close to grazing. The degree of sparsity can be increased for physically larger bodies and this is illustrated in Figures 4 and 5 where bodies of identical cross-section as in Figure 1 are considered but with broadside dimensions doubled and tripled, i.e. 30A and 45A respectively. With the threshold factor set to 0.5%, the sparse matrix is 11% full for the 30A body and 7% full for the 45A body but the predicted RCS is again computed very accurately everywhere except near grazing where a maximum error of 3dB is observed. We next consider the RCS of another 2D scatterer which in contrast to the previous target, contains a sharp edge. The subject scatterer is an isosceles triangular cylinder with a circular end-cap as shown in Figure 6. A thresholding factor of 0.3% results in a 21% full matrix and backscatter RCS computed using this sparse matrix and the iterative refinement procedure of (14) is shown in Figure 7. While the first sparse matrix iteration yields considerable error at the nulls, the refinement technique corrects this in the second iteration, displaying an error of only -0.5 dB at ( = 180~ (normal to the circular end-cap) and an error of only 1 dB at ( = 0~ (edge-on). A third body considered was the circular cylinder with a 10A diameter. In this case, a threshold factor of 0.5% results in a 23% full matrix and from the RCS computation in Figure 8 it is seen that the refinement procedure corrects the slight oscillations seen in the first iteration. 7

40 - - l ' ' | * l| *- l| * l- - 35 -- Exact --------- 7% Full 30 1 -C/ 0 - 20 U u~10 0 15 30 45 60 75 90 Observation angle (deg.) Figure 5: H-pol backscatter RCS of a 45A x 2A rectangular cylinder with circular end-caps 2.2 E-polarization In the case of TM incidence, the lack of diagonal dominance in the impedance matrix necessitates an alternative to (12). Since the currents computed using (12) had large discrepancies for this polarization, the physical optics currents were employed as an initial estimate and the iterative refinement technique was employed to improve on this estimate. The iterative refinement technique performs the following mathematical operation, corresponding to a first iteration [ZS]{J()} = {V} -[ZR]{J(PO)} (22) where {J(l)} is a first approximation of the current and {J(PO)} corresponds to the physical optics currents. Subsequent iterations are implemented according to (15). From the RCS computations for the rectangular and triangular cylinder in Figures 9 and 10, it is seen that the second iteration is almost coincident with the exact solution. This appreciable improvement is attributed to the fact that the residual matrix [ZR] contains entries of orders larger than the corresponding matrix for the H-pol. The currents being iterated in this manner play a more significant role in the correction process. For the circular cylinder it is of more interest to examine the currents, and these are depicted in Figure 11 for an 8A diameter circular cylinder. Here the exact and PO currents as well as the iteration currents are plotted as a function of the segment number on the cylinder. The numbering of the segments on the cylinder is done in a counter-clockwise direction with segment number 1 and 300 being at the extremities of the face exposed to the illumination. It is seen that the second iteration predicts currents which are nearly identical to the exact. 8

y + I I I I I- - - - - - - ---- x 12 t Figure 6: Geometry of the isosceles triangular cylinder with a circular end-cap 20 10 r0 0 X -10,a O it -20 -30 -40 0 30 60 90 120 150 180 Observation angle (deg.) Figure 7: H-pol backscatter RCS of a 12A isosceles triangular cylinder closed with a circular end-cap 2A in diameter 9

15 ~ I ~ I ~ w ~ ~ I.... I I. I ~ ~ ~ I.. I I I i e 10 CU C) 0 fP ca Exact -- 23% Full- Iter #1 23% Full - Iter #2.. I.... I.... I.. I.... I...... 0 15 30 45 60 75 90 Observation angle (deg.) Figure 8: H-pol backscatter RCS of a 10A diameter circular cylinder C-, co CU Cl 40 35 30 25 20 15 10 5 0 -5 0 15 30 45 60 75 90 Observation angle (deg.) Figure 9: E-pol backscatter RCS of a 15A x 2A rectangular cylinder with circular end-caps 10

30 PI-9 c) rd 0 P3 20 10 0 -10 0 30 60 90 120 150 180 Observation angle (deg.) Figure 10: E-pol backscatter RCS of a end-cap 2A in diameter 2.50.... I,..... 2.00 - 12A isosceles triangular cylinder closed by a circular I to c) 1.50 1.00.50.00 0 100 200 300 400 500 Segment number 600 Figure 11: E-pol currents on a 8A diameter circular cylinder 11

3 Three dimensional problems The governing integral equation for solving 3D scattering problems is given in [6]. The approach combines the advantages of triangular patch modeling and the EFIE formulation resulting in a simple and efficient algorithm. The electric current is discretized in the form of special basis functions associated with each nonboundary edge of the triangulation. Galerkin's method is then used as the testing procedure to obtain the final system of linear equations. The resulting dense system is inverted to obtain the unknown current values for different excitation vectors. As mentioned in the section on the 2D implementation, Isomom zeroes out the nonsignificant values of the impedance matrix using a tolerancing technique. The tolerancing is much simpler to implement in 2D rather than in 3D owing to reduced geometric complexity. In 3D, the matrix elements are normalized with respect to the diagonal entry of the corresponding row and zeroed out if their value falls below the user specified tolerance. We refer to this type of tolerancing as the magnitude criterion. The distance criterion which was successful in 2D only works for the simplest of planar 3D geometries, principally because the physical nearness of two points on a closed surface does not necessarily imply strong field interaction. 3.1 Results 15 -I,- 5 - 0 -5 - 1< b -10 --15 Exact 64% full 40%fuIlI -20 -? t71% -l 56% full - -25 -....,.....,....-. 0 30 60 90 Observation Angle O., deg. Figure 12: aso backscatter pattern of a 1A square plate (lying on the x-y plane) The perfectly conducting plate was the first 3D geometry Isomom was tested on. Geometrically, it is also the simplest since the structure is planar and the physical separation between edge-pairs is more meaningful. We computed the backscatter patterns from a 1A square metallic plate for the 00 polarization. In Figure 12, the exact solution is compared with progressively sparser approximations of the impedance matrix. It is found that even 12

with 44% sparsity, the deviation from the exact values are extremely small. We employed a rectangular filter for the truncation and observed that a trapezoidal filter of the same width degraded the final result. The tolerancing was carried out by the distance criterion since the geometry is confined to a single plane. 25 20 < tm 10 -5 - \...... b 0 Exact 51% full 25% ful \ -, ~~~........................... -10 80% full 43%full -15 -0 30 60 90 Observation Angle B., deg. Figure 13: cra backscatter pattern of a 2A square plate (lying on the x-y plane). A rectangular filter was employed for truncation. We next consider the backscatter RCS of a 2A square to test our assumption that tolerable sparsity would increase as the problem gets larger. In Figure 13, we plot the backscatter pattern for the 2A plate for the more difficult q( polarization. A rectangular filter was again employed. It is observed that the backscatter pattern is accurately predicted, even when the impedance matrix is 43% full. In the instance when 25% of the impedance matrix is non-zero, the far-field values deteriorate near edge-on and the dips in the backscatter pattern are not picked up very well. The next geometry in our testbed was a perfectly conducting cylinder. A cylinder was chosen since broadside incidence excites creeping waves; the accurate prediction of such phenomena accurately would be a big test for Isomom. For a cylinder, the distance criterion has obvious limitations. Calculating the distance between two points on the cylinder and accepting or rejecting elements on that basis, may lead to non-physical coupling and incorrect answers. Since the main idea of Isomom rests on the observation that the coupling between matrix elements decrease with increasing distance, we normalized off-diagonal elements with respect to the diagonal element and zero out those elements which fall below a user-specified tolerance. The higher the value of the user-specified tolerance, the sparser the matrix. This tolerancing is thus done by the magnitude criterion. In the case of the cylinder, we tried two schemes of tolerancing: combined magnitude/distance criteria and solely magnitude criterion. Figure 14 plots the backscatter pattern from a perfectly conducting circular cylinder of radius 0.5A and height 1.5A. The tolerancing is done by the combined magnitude/distance scheme. The next figure, Figure 15, shows the backscatter pattern for the same geometry 13

15 -20 im u -5 -0 30 60 90 -15 -Exact 37% full -20 -31% full -25 -0 30 60 90 Observation Angle 6O, deg. Figure 14: ass backscatter pattern of a PEC cylinder (radius=0.5A; height= 1.5A). Tolerancing done by magnitude/distance criteria. and for the same polarization, but the tolerancing is carried out solely by the magnitude criterion. It is observed that both schemes provide the same degree of accuracy. However, since it is difficult to implement the distance criterion for general 3D geometries, we will henceforth consider only the magnitude criterion for tolerancing. In Figure 15, it is observed that a system with 64% sparsity yields decent far-field values. Next, we use the iterative refinement algorithm to see whether the result converges to the correct value. In Figure 16, we plot the backscatter pattern for the same cylinder for different iterations. The matrix is now 37% full and it is observed that the error oscillates around the reference result between iterations - as predicted by the derivation in Section 2 - but converges to the correct value. However, if the error in the first iteration is already very high, i.e., the matrix is too sparse, then the convergence to the correct result is too slow to be justified computationally. In the 2D case, drastic improvements were seen on using the Physical Optics (PO) currents as a initial guess. In our case, the result does not change at all. This is probably because the remainder matrix in the 3D case contains elements which are about three orders of magnitude lower than the diagonal entries: thus, the effect of PO currents is negligible at best in the far field. The last geometry considered in this report is a metallic almond of length 3A. It was noticed that the characteristics of Isomom and the conclusions derived from them for the PEC cylinder hold true for the almond. In Figures 17 and 18, the exact solution computed for the q)q polarization is compared with the Isomom solution for varying degrees of matrix fill. It is observed that the solution obtained through Isomom is reasonably accurate for almost all incidence angles. Even a 66% empty matrix predicts excellent far-field values for a wide range of incidence angles. However, the solution starts to degrade at near grazing incidence, an expected phenomenon in Isomom since the traveling wave information is lost when small magnitude matrix elements are zeroed out for achieving sparsity. If near normal 14

15 -20 - 52f 26%ful 5 0/ " - 5" -15\ Exact 36%fii -20 - 52% full- 26% full_ -25 -0 30 60 90 Observation Angle 6., deg. Figure 15: ass backscatter pattern of a PEC cylinder (radius=0.5A; height= 1.5A). Tolerancing done by magnitude criterion. backscatter values are desired, the Isomom method may be able to tolerate very low sparsity levels without a significant degradation in the computed result. 4 Conclusion Isomom is a technique which solves a sparse version of the moment method matrix. This leads to an approximate solution of the scattering problem but the solution can be improved by successive iterations which account for the far-zone interactions of the moment matrix elements. The iterative refinement algorithm is an attractive feature of Isomom and it was already shown in Section 2, using simple matrix concepts, that the iterative procedure converges to the exact value in the limit. Typically, though, only one more iteration is necessary to yield a solution of good accuracy. Isomom was tested on three types of 3D geometries and a wide variety of 2D structures till now. Clearly, the methodology has obvious merit for solving large scattering problems. Basically, Isomom performed as expected without any major surprises. Most importantly, unlike other techniques, such as the fast multipole method and wavelets, Isomom is easily extended to 3D analysis. Moreover, tolerable sparsity will only increase as the problem size gets larger. This was clearly demonstrated in the 2D applications of Isomom where it was possible to compare with the exact solution. It was found that, generally, only elements within a radius of 0.75A need be kept for tolerable accuracy during the first iteration. Unfortunately, because of the large resources required for inverting large full matrices, it was not possible within the course of this short study to examine the performance of Isomom on large 3D structures. Nevertheless, the large 2D body simulations provide sufficient proof that Isomom's accuracy will improve with body size and greater than 95% sparsity can be achieved for targets 3A to 4A long (The projected savings in memory and CPU time are in shown in Figures 19 and 20). Such sparsities are, of course, 15

15 5 - a Cltr t I' A' /~ '^" Ii~.~. lid /: '/ U -10 _ _ _ Exact -20 -20 -- - Iter#2 2 Iter#3 -25.,.........,,.. 0 30 60 90 Observation Angle 00, deg. Figure 16: oera iteratively refined backscatter pattern of a PEC cylinder (radius=0.5A; height= 1.5A). Impedance matrix is 36% full. necessary to allow for fast calculations of large 3D structures (30-50A long). 4.1 Future Plan The following tasks are suggested for implementation in the near future * Application of Isomom to large 3D structures. To apply Isomom to large 3D structures, it is necessary to develop a new code which retains in core only the non-zero elements of the sparse moment method matrix. Also, the iteration need be performed without a need to bring in core any of the far-zone elements. Actually, only the multiplication of the current vector with the matrix row is needed, and this can be dlone efficiently on-the-fly without calculating the explicit far-zone matrix elements. We anticipate that a 20A long 3D structure can be calculated using Isomom over the next year. * Material modeling. Since it employs standard moment method matrices, existing codes such as CARLOS3D (which incorporates material modeling) can be used to construct the sparse Isomom matrix. This is indeed one of the attractive features of Isomom. * Combined Isomom-xpatch implementation. Although, high frequency methods can compute the contribution from large smooth section of the structure, they lack accuracy in modeling surface details, irregularities and junctions as well as material. As was already demonstrated, Isomom can be used to improve upon the physical optics result. Since the physical optics results is typically a good approximation over large smooth surfaces, this kind of hybridization removes a major part of the computational stress. 16

* Combined Isomom-finite element implementation. There are several advantages in combining Isomom with the finite element method. Basically, Isomom is used for truncating the finite element mesh very close to the target, thus, minimizing the volume of the computational domain. In contrast to the usual boundary-integral implementation which results in a partly full, partly sparse matrix, Isomom is associated with sparse matrices and the entire finite element-Isomom matrix is therefore sparse. On the other hand, the finite element method has the attractive features of geometrical adaptability, robustness and ease of material handling. References [1] R.F. Harrington. Field computation by moment methods. Macmillan: New York, 1968. [2] F.X. Canning. Fast integral equation solutions using GTD-like matrices. Radio Sc., 29(4):993-1008, Jul-Aug 1994. [3] V. Rokhlin. Rapid solution of integral equations of scattering theory in two dimensions. J. Comput. Phys., 86(2):414-439, 1990. [4] J.L. Volakis. EECS 633 course pack. University of Michigan, 1993. [5] Capt. Paul Skinner. Personal communication. Air Force Institute of Technology, 1994. [6] S.M. Rao, D.R. Wilton, and A.W. Glisson. Electromagnetic scattering by surfaces of arbitrary shape. IEEE Trans. Antennas Propagat., 30(3):409-418, May 1982. 17

mQ il.b 10 5 0 -5 -10 -15 -20 -25 -30 Observation Angle So, deg. Figure 17: ao0 backscatter pattern of a 3A long almond (azimuth cut). im.O.l< b 20 10 0 -10 -20 -30 -40 -90 -60 -30 0 30 60 90 Observation Angle 8., deg. Figure 18: ar, backscatter pattern of a 3A long almond (elevation cut). 18

/-N 0 'S 30 25 20 15 10 5 0 0 500 1000 1500 2000 2500 3000 3500 4000 4500 5000 Number of unknowns Figure 19: Comparison of memory requirement for the moment method and the new technique 20000 C-:4 F&. 15000 10000 5000 0 0 2000 4000 6000 8000 10000 Number of unknowns Figure 20: Comparison of solution times for the moment method and the new technique 19