025921-24-T Application of. Edge-Based Finite Elements and Vector ABCs in 3D Scattering A. Chatterjee, J. M. Jin and J. L. Volakis National Aeronautics and Pacific Missile Space Administration Test Center Ames Research Center Pt. Mugu CA Moffett Field CA 94035 03042-5000 January 1992

t* st: yOti rjC~C 025921-24-T TECHNICAL REPORT for NASA Grant NAG-2-541 NASA Technical Monitor: Alex Woo Grant Title: Development of 3D Electromagnetic Modeling Tools for Airborne Vehicles Report Title: Application of Edge-Based Finite Elements and Vector ABCs in 3D Scattering Institution: Radiation Laboratory Department of Electrical Engineering and Computer Science The University of Michigan Ann Arbor MI 48109-2122 Period Covered: September 1991 - February 1992 Report Authors: A. Chatterjee, J.M. Jin and J.L. Volakis Principal Investigator: John L. Volakis Telephone: (313) 764-0500

APPLICATION OF EDGE-BASED FINITE ELEMENTS AND VECTOR ABCs IN 3-D SCATTERING A. Chatterjee, J.M. Jin and J.L. Volakis Radiation Laboratory Department of Electrical Engineering and Computer Science University of Michigan Ann Arbor, MI 48109-2122

Table of Contents 1. Introduction 1 2. Derivation of the ABCs 2 2.1 Two-Dimensional ABCs 2 2.2 Three-Dimensional ABCs 3 3. Formulation 7 3.1 Two-Dimensional Case 7 3.2 Three-Dimensional Case 7 4. Results 9 5. References 10 Figures follow page 10

1 Introduction Open region scattering problems are frequently modelled using an integral equation formulation[1] and solved using the method of moments. For an inhomogeneous dielectric scatterer, the integral equation must be discretized over the entire volume of the object. Since this leads to a full matrix system, the integral equation technique proves to be expensive in terms of storage and solution time. Differential equation methods, such as finite elements, have therefore been tried as an alternative solution approach. Although these techniques give rise to matrix systems that are considerably larger than those generated using the integral equation formulation, the finite element matrices are highly sparse leading to smaller storage and lower solution time. These methods, however, can only solve bounded field problems whereas most electromagnetic problems are open boundary-infinite domain types. To solve for unbounded field problems, the finite element mesh needs to be truncated artificially at some distance from the scatterer. A wide variety of methods like ballooning, harmonic series expansions and infinitesimal scaling have been developed to accomplish this. However, all of the above methods are'global' techniques which produce fully populated submatrices that are computationally expensive and spoil the sparse, banded structure of the finite element system. The most promising method of mesh termination developed so far has been the application of an absorbing boundary condition (ABC) on an artificial outer boundary to minimise the non-physical reflections from the boundary. The essence of the method is to force the field components at the mesh termination plane to satisfy the differential equation (absorbing boundary condition) for an outward travelling wave. The degree of accuracy of the finite approximation to the differential equation determines the quality of the mesh termination. Stability considerations limit the degree of approximation. The advantage of the ABCs over global methods is that they are local boundary conditions and hence retain the sparse structure of the finite element formulation. Moreover, the additional computational effort when using ABCs is small when compared to a bounded field problem. One disadvantage, however, is that ABCs are approximate and do not model the exterior field exactly. The local boundary conditions should also give rise to inaccuracies when travelling waves are excited on the scatterer and substantial global coupling exists between widely separated parts of the scatterer. In this report, we consider a FE-ABC solution of the scattering by arbitrary three-dimensional structures. The computational domain is discretized using edge-based tetrahedral elements. In contrast to the node-based elements, edge elements can treat geometries with sharp edges, are divergenceless and easily satisfy the field continuity condition across dielectric interfaces [1]. They do, however, lead to a higher unknown count but this is balanced by the greater sparsity of the resulting finite element matrix. Thus the computation time required to solve such a system iteratively with a given degree of accuracy is less than the traditional node-based approach [2].

The purpose of this report is to examine the derivation and performance of the ABCs when applied to two and three dimensional problems and to discuss the specifics of our FE-ABC implementation. The two dimensional ones investigated here are the well-known Bayliss-Turkel[3] and Engquist-Majda[4] absorbing boundary conditions. The three dimensional ABCs presented here are those derived by Mittra[5] and Kanellopoulos and Webb[6]. All the ABCs discussed here are derived from the Wilcox expansion and can only be applied on a spherical outer boundary. Some results are then presented which demonstrate that remarkably accurate solutions can be obtained by enforcing the ABC a small fraction of a wavelength from the scatterer. This is in contrast to our experience with two-dimensional applications and is probably due to the faster (1/r rather than 1/r/F) decay of the propagating fields. 2 Derivation of the ABCs 2.1 Two-dimensional ABCs Absorbing boundary conditions for two-dimensional problems have been extensively studied over the years. The most notable ABC to come out of this research has been the one proposed by Bayliss, Gunzburger and Turkel[3] who employed an asymptotic analysis to derive a series of local operators. Using the pseudo-differential operator theory, Engquist and Majda[4] have generated a set of different operators for minimisation of the reflections from the mesh termination boundary. However, since the Engquist-Majda ABCs are not as accurate as the Bayliss-Gunzburger-Turkel(BGT) ABCs, only the BGT operators will be discussed here. Consider a perfectly conducting cylindrical scatterer, shown in Figure 1, whose cross-section is defined by the contour Pl. Let the exterior region of the scatterer lie in the domain Qf. For a TM-polarized incident wave, we need to solve the wave equation V2u + k2U = 0 (1) in Q. The wave function u is proportional to the z-component of the scattered electric field and satisfies the boundary condition ui+u = 0 on rl (2) where ui is the incident wave function. Following Wilcox, an asymptotic expression for u can be written as follows: u(p, a) = a F ) 2() JP n=o p

Defining up as the derivative of u with respect to p, we have from (3) u+k = -____pF aa(q) a2(4) 23/2 a() + 3 + 5 2 + (4) 2p3 / p Therefore, from (4), we obtain u, + jku [/2] Thus if we neglect terms of the order O(p-3/2) and smaller, we obtain up + jku = 0, which is equivalent to the Sommerfeld radiation condition for the wave function u in two dimensions. It is shown next that a higher-order boundary operator B1 can be obtained that yield terms of the order O(p-5/2) when applied to u. The operator B1 is given by B1 = +jk+ (5) Op 2p It is readily found that, for a given p, B1 introduces a higher order error in p-' than does the Sommerfeld radiation condition such that terms of the order O(p-5/2) and smaller are neglected. Continuing along similar lines, the next higher order operator B2 can be derived by first defining v = Blu, and then showing that + jk++- v = O(p-/2) The second-order absorbing boundary operator B2 is thus given by B2 = -(0 jk~+- +j )(+ ) (6) In fact, it is shown in [3] that a generalised operator Bm can be constructed by repeating the above procedure such that _ a 2p - 3/2 Bm = i [+jk+ P ] (7) and Bm u = [2i/2] p2m+1/2 Since the error between the exact field and the approximated value decreases as we employ higher order boundary conditions, the mesh termination boundary can be brought closer to the scatterer resulting in a smaller number of unknowns. However, the higher order conditions make the system of equations progressively more ill-conditioned and the choice of the order of the

ABC is often limited by this constraint. The second order ABC is usually employed in practice. An alternate derivation[5] of the boundary operators can also be carried out by imposing the requirement that u(p, 0) in (3) satisfy the wave equation and then deriving a recursion relationship in terms of the angular derivatives. The second order ABC thus reduces to ( (-jk 2p 8jkp2+ 8k2p3) u + 2jkp2 +2k2p3) U + (p9/2 c(p)u + 3(p)uoo (8) Higher order boundary operators in the above formulation involve higher order angular derivatives and are thus simple to implement numerically. 2.2 Three dimensional ABCs Three dimensional absorbing boundary conditions can be classified into two categories: scalar ABCs and vector ABCs. Most of these boundary conditions are untested till date; it is, therefore, difficult to comment on the performance of these boundary operators. However, in tests carried out using the vector ABC derived in [6], reliable results have been obtained on placing the mesh truncation boundary a fraction of a wavelength distant from the scatterer. The first section presents a scalar ABC proposed in [5] and the second section is devoted to the formulation of vector ABCs derived in [5] and [6]. 2.2.1 Scalar ABC Let u(r, 0, () satisfy the scalar wave equation (1) and be expressed asymptotically as — jkr oo an(sl ) u = (, r (9) rn n=O It can then be shown that the coefficients in the series expansion, an, satisfy the following recursion relationship -2jknan = n(n- 1)anl + Dan-1 (10) where D is Beltrami's operator and is given by 1 ( O) 1 02 D = sn&(sin0) + si20$ (11) Differentiating (9) with respect to r and incorporating the the recursion relationship repeatedly in the resulting expression yields ur = -jk {c(r)u +/3(r)Du} (12) up to and including terms of order r-4. In (12), ca(r) 1 + j-r /3(r) 2(kr)2cQ(r) ( 1) 1~

2.2.2 Vector ABC As for the scalar case, we begin with the Wilcox representation of the electric field which has an expansion of E(r) -jr r A=(0,) (13) r - n where r, 0, q are spherical coordinates. From (13), we get 1 D x+D1 e-jk ~ n~hA,0 V7xE-= {- jk - x + E- 2 E kn (14) V E yk rx+ E r2 r where Ant = j x An is the transverse component of An and, for a vector F, D1F is given by 1 [a aF9] D1F = sin (sinOF[) - + s in sinOF + as+ (15) Using the recursion relation -2jknAnt = n(n- 1)An-l,t + D4An-l where D4An = n sDA~D An) + (DAX +D, An)c D A A 1 A 2cos0 A$n D A, = 2A 1 AO n -50 sin20 n sin20 ao 2 0A - 1i 2cos0 0A' D, - sinO 0q sin20 n + sin20 (oq and D is Beltrami's operator defined in (11), we can derive the representation correct to r-4. Applying the recursion relation in ( 14) yields the desired relationship for the vector ABC: V x E = a(r)E + (r)D4E (16) where I 1 at (r) -j k 1+ r:2jkr2 (1 - 1/jkr) An alternate derivation of the three-dimensional vector ABC is given in [6] which has a higher order of accuracy and is easier to implement than the 5

previous one. The basic building block is the differential operator LN, defined as LN(U) = x' x u-(jk+-Lu, N = 0,1,2,... (17) where r is a unit vector in the radial direction. It can be shown that for N > 0 and n > 0, rn+l) = (n- rn+2 LN(v {gr(9 ) - (n +l- N)Vt {gAnr( q)} (19).n+. rn+2 where g is e-jkr and the subscripts t and r denote the transverse and radial components of a vector, respectively. In both cases, LN has the effect of multiplying by 1/r while leaving the transverse dependence unchanged. The operators BN, N = 1, 2,... can now be defined such that BN(E) ( 2N+i) where BN (g ) (n + 1 - N)(n + 2 - N).. (n)g ) rn+l =rn+l+N +s(n + -N)(n + 2-N)...(n-1)Vt g nr(+N (20) Since AOr is zero, it can be seen that the RHS of (20) vanishes for n - 0, 1,..., N - 1, i.e., BN annihilates the first N terms of the vector expansion (13). Thus, BN = 0 is an approximate absorbing boundary condition on the surface S of a sphere of radius r. After a certain amount of algebra and making use of the fact that E satisfies the vector wave equation, the operators B1 and B2 can be written as B1(E) = x V x E-aEt + (s-1)VtEr (21) B2(E) = (r x V x E)+ Et + V x [r(V x Er] +(s - 1)Vt(V Et) + (2 - s)aVtEr (22) where a = jk and,3 = 1/(2jk + 2/r). The operator B2 can be incorporated into a weighted residual formulation for s = 1 but leads to an unsymmetric matrix problem. However, for s = 2, the functional yields a sparse, symmetric matrix. Thus, the second order boundary operator can be expressed as B2(E) = -(i x V x E)+P2(E)

where P2(E) = aEt + fV x [ri(V x E)r] + (s - 1)Vt(V * Et) + (2 - s)acVtEr(23) which can be easily incorporated into the surface integral of the electric field functional. 3 Formulation 3.1 Two-dimensional case The scattered field must satisfy the wave equation (1) within the domain of interest. The first step in the finite element formulation is to multiply (1) with a test function v and integrate the product over the entire problem domain. This gives j (vV2u + k2vu) dS = O (24) where u consists of a set of known basis functions weighted by its corresponding unknown coefficients. The second step involves transferring the derivative from the basis function u to the testing function v via Green's identity as follows VV2U = -VU* VvdS +] v-dl (25) Substituting the above identity into (24), we arrive at the weak form of the Helmholtz equation | (Vu Vv-k uv)dS = a dl (26) where rl and r2 describe the boundary of the solution domain as shown in Figure 1. The absorbing boundary condition for two dimensions is to be applied at the contour boundary ('2) and must be expressed in terms of the normal derivative of the wave function. For a Galerkin formulation of (26), we replace the testing function v with the basis function u. 3.2 Three-dimensional case For simplicity, we consider the problem of 3-D scattering by a conducting body. To solve this problem using the finite element method, it is necessary to enclose the scatterer within a fictitious surface denoted by S. Since the ABCs are derived only for spherical surfaces, the fictitious outer surface in our case is usually a sphere. The scattered field, denoted by E, satisfies the Helmholtz vector equation interior to S and the absorbing boundary condition at S. The equivalent variational problem for this is given by SF(E) = 0O

where F denotes the functional given by [6] F(E) = [( x E) (V x E)- k2E. E] dV + j [aE, + (V x E)2- (V Et)2] dS (27) in which V denotes the volume enclosed by S. To discretize this functional, the volume V is subdivided into a number of small tetrahedra, each occupying the volume Ve (e = 1,2,...,M), where M denotes the total number of elements. Within each element, the electric field is expressed as m E= e E;W- = {We}T{E)} = {Ee}T{We} (28) j=l where Wq are the edge-based vector basis functions given in [2], Ee denote the expansion coefficients of the basis, m represents the number of edges in the element and the superscript stands for the element number. On substituting (28) into (27), we obtain M Ms F = {ET }[Ae]{Ee} + Z]{E }T[Bs]{Es} (29) e=l s=l where M, denotes the number of triangular surface elements on S. Also, the elements of the matrices [Ae] and [Bs] are given by Aj= [(V x W?). (V x we)- kW. W;]dV BAi=J aRs -.W I (N7-. We''V C W t) d B,~J=j i [aw t W~ + (V x Wt)r. (V x W;)r - (V W, )(V. W;t)]dS where Ve denotes the volume of the eth tetrahedron and S, denotes the surface area of the sth triangular surface element on S. On assembling all M elements, (29) can also be written as F = E}T[KI]{E} where [K] is a N x N sparse, symmetric matrix with N being the total number of edges in the geometry and {E} is a N x 1 column vector denoting the edge fields. The system of equations is then obtained by applying the Rayleigh-Ritz procedure which amounts to differentiating F with respect to each edge field and then setting the differential to zero. This yields [K]{E} = {0} (30) Upon imposing the boundary condition at the surface of the scatterer which states that the total tangential electric field vanishes at the conducting surface, we obtain the final system whose solution yields the field components over the entire domain. In our implementation, only the non-zero elements of the [IK] matrix were stored and the resultant system was solved using the biconjugate gradient algorithm.

4 Results Finite element codes have been extensively tested using two-dimensional absorbing boundary conditions and have been known to perform quite well for circular mesh termination boundaries. In figure 2(a), we consider a cylindrical problem[7] to compare the BGT and Engquist-Majda ABCs. This problem is interesting since strong resonant phenomena are observed for small bodies. It consists of a conducting cylinder with a radius of one free-space wavelength covered by a half free-space wavelength thick dielectric with a relative permittivity of 2 and a relative permeability of 2. Figure 2(b) shows the comparison of the bistatic pattern (0i,, = 1800) of two solutions of the BGT ABCs, with the outer boundary at a distance of 1.5A and.5 A away from the scatterer. As expected, the 1.5A distant boundary produces better results than the.5A distant one. Figure 2(c) plots the performance of the EngquistMajda ABCs for the same geometry and angle of incidence. The BGT ABC with its outer boundary a mere.5A away outperforms the Engquist-Majda ABC when its outer boundary is 2.5A distant from the scatterer. Figure 3(a) analyses two conducting wedges, 5Ao long, with the area between them filled with a dielectric having er = 2 and u, = 2. The incident angle of the incoming plane wave is 0~. Again, the results obtained using the BGT absorbing boundary condition in figure 3(b) matches the data obtained from a hybrid finite element / boundary element analysis. Figure 3(c) shows the corresponding results obtained with the Engquist-Majda ABCs. In figure 4, we present the backscatter data from a 0.3A x 0.3O plate. The Sommerfeld radiation condition was employed at the mesh termination boundary placed 0.3A away from the edge. The comparison with a standard integral equation code is seen to be excellent. The backscatter pattern of the same plate is shown in figure 5 with a conformal outer boundary. The zero thickness plate lying at z = 0 on the x - y plane was enclosed within a boundary resembling a cut-off sphere. This was achieved by making a plane cut at z = +.3A of a.45A radius sphere and setting r in (23) to oo over the planar surfaces. Storage was reduced by about 10% on using the conformal outer boundary but the reduction should be more striking for larger geometries. Figure 6 compares the measured bistatic cross-section (Oi = 180, i - 900) of a metallic cube having an edge length of 0.755A with the corresponding pattern computed by the three-dimensional FE-ABC code. The secondorder vector ABC was employed at the mesh truncation boundary which was placed only 0.1A from the edge of the cube. About 33000 unknowns were used for the discretization of the computational domain and the [A] matrix contained a total of 264000 distinct non-zero entries. The storage requirement of this matrix is consequently much smaller than the 2000 unknown system generated by the boundary integral approach which has 4 million non-zero entries. Figure 7 presents backscatter data for a cylinder of radius 0.3A and height 0.6A. The data from the three-dimensional finite element code again compare

well with that obtained from a moment method-body of revolution code. The mesh was terminated at a distance of 0.3A from the edge of the scatterer and the system consisted of about 33000 unknowns and converged to the solution in about 350 iterations when using the Sommerfeld radiation condition. Each iteration took approximately 0.3 seconds on a Cray YMP and on the average it was found that for N > 25000, the number of iterations required was of the order of N/100. Figure 8 shows the backscatter pattern on enclosing the pec cylinder in figure 7 with a conformal cylindrical boundary capped at the ends by a spherical surface. The storage was reduced by about 20% but the computed values are off by about 1dB for most incidence angles. The poor results are probably due to the fact that the ABCs were derived only for spherical surfaces and therefore fail on cylindrical boundaries. The above results point out that the mesh truncation boundary could be placed a few tenths of a wavelength away from the scatterer to obtain reliable results. This is probably because in three dimensions, the scattered fields decrease as 1/r whereas in two dimensions, the rate of decrease is proportional to 1/v'T. 5 References 1. A. Bossavit, "A rationale for'edge-elements' in 3-D fields computations", IEEE Trans. Magn., vol. 24, no. 1, pp. 74-9, Jan. 1988. 2. M.L. Barton and Z.J. Cendes, "New vector finite elements for threedimensional magnetic field computation", J. Appl. Phys., vol. 61, no. 8, pp. 3919-21, April 1987. 3. A. Bayliss, M. Gunzburger and E. Turkel, "Boundary conditions for the numerical solutions of elliptic equations in exterior regions", SIAM J. Appl. Math., vol. 1, no. 3, pp. 371-85, September 1980. 4. B. Engquist and A. Majda, "Radiation boundary condition for the numerical solution of waves", Math. Comp., vol. 31, no. 139, pp.629-51, July 1977. 5. R. Mittra, O. Ramahi, A. Khebir, R. Gordon and A. Kouki, "A review of absorbing boundary conditions for two and three-dimensional electromagnetic scattering problems",IEEE Trans. on Magnetics, vol. 25, no. 4, pp. 3034-9, July 1989. 6. J.P. Webb and V.N. Kanellopoulos, "Absorbing boundary conditions for finite element solution of the vector wave equation",Microwave and Opt. Tech. Letters, vol. 2, no. 10, pp. 370-2, October 1989. 7. J. D'Angelo and I.D. Mayergoyz, "On the use of local absorbing boundary conditions for RF scattering problems", IEEE Trans. on Magnetics, vol. 25, pp. 3040-2, July 1989.

Incident Wave e.c. scatte rer Figure 1: Geometry of p.e.c scatterer enclosed in an outer circular boundary ~~~~~~~~~~~~- I

v=2,; =2 incident field5 incident field(a) 20 -Closcd Form B aylissl-Turcl ABC. 1.5 ) B 8ayliss-Turkcl ABC,O.5. 10 0 5- GEquit-Majdh ABC. z5 X fto. 0 - - O 30 60 90_ ~ 150 180 (C) Figure 2: (a)Cylinder geometry (b)Bistatic pattern using BOT absorbing boundary conditions (c)Bistatic pattern using Engquist-Majda absorbing boundary conditions

conducting wedge conducting wedg I I incident field 5Xo (a) 20 10 e' 0 30 60 120 150 180 (b) 10CL~ ~() boundary conditions (c)Bistaticpattern using Engquist-Majda absorbing boundary conditions

Backscatter pattern of a.3Xx.3X plate 0 I -.:...... C)o it |. | BiCG,=0 ) |. | 0....................................-j..F.:..0)....... *c FEM( 0) -50....-................. -.30 - 0 15 30 45 60 75 90 Observation Angle 0, deg. Figure 4: Backscatter pattern of plate using Sommerfeld radiation condition

~~~Z Y~~~~~~~z y x Top view Side view Backscatter pattern of a.3Xx.3X plate 2nd order ABC; cut-off sphere'0 -0 Q) - 20........................................................................................ 0 z_30't s 0.......C... )."... (...i............ 3 40 i EMH <.....;............, FEM( =O)' BCC.(H.,=O). 0 30 60 90 Angle(0) in degrees Figure 5: Backscatter pattern of the plate in figure 4 enclosed by a conformal outer boundary and using the 2nd order absorbing boundary condition.

Bistatic pattern of a.755 A metallic cube O,=180o -15 I.cu'-'. 30 -. a i Observation Angle 8x, deg. Figure 6: Bistatic pattern (#i = 1800) of a metallic cube having an edge length of 0.755A using the 2nd order vector absorbing boundary condition. length of 0.755A us-ing the 2nd order vector absorbing boundary condition.

Backscatter pattern of a metallic cylinder base radius=.3>' height:.6A 5DO A 0 o -,~'' Cl~cIe(E =0) Q4~~~~~~~~~) *~~ ~: sABC(E= —o) -10 Ci9(H~9~)........................................' g.,.0' ABC(H =0) 0 jO 60 90 Observation Angle 8~, deg. Figure 7: Backscatter pattern of a metallic cylinder of radius 0.3A and height 0.6A using the Sommerfeld radiation condition. The axis of the cylinder coin ~: with the z-axis of the coordinate system.

I. - v Backscatter pattern of pec cylinder (r.3h,h=.6X) 0 30 60 90 2nd ormal outer boundary and using the second order absorbing boundary reduction3 9015 02514 8472 #<), 5............:to'P ~ FE:(H =0) 0 30 60 90 Angle(O) in degrees Figure 8: Backscatter pattern of the cylinder in figure 7 enclosed by a con