375458-4-T MODELING OF PERIODIC ARRAYS USING ACCELERATED FINITE ELEMENT-BOUNDARY INTEGRAL METHODS AND HIGHER ORDER ELEMENTS D. Filipovic T. F. Eibert L. S. Andersen J. L. Volakis Rockwell International Science Center 1049 Camino Dos Rios P.O. Box 1085 Thousand Oaks, CA 91360 January 1999 375458-4-T = RL-2519

PROJECT INFORMATION PROJECT TITLE: REPORT TITLE: U-M REPORT No.: CONTRACT START DATE: END DATE: DATE: Hybrid Finite Element-AIM Methods for Antennas Modeling Of Periodic Arrays Using Accelerated Finite ElementBoundary Integral Methods And Higher Order Elements 375458-4-T October 1998 March 1999 February 1, 1999 SPONSOR: Marek Bleszynski Rockwell International Science Center 1049 Camino Dos Rios P.O. Box 1085 Thousand Oaks, CA 91360 SPONSOR CONTRACT No.: U-M PRINCIPAL INVESTIGATOR: CONTRIBUTORS TO THIS REPORT: 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/ D. Filipovic, T. F. Eibert, L. S. Andersen and J. L. Volakis

375458-4-T MODELING OF PERIODIC ARRAYS USING ACCELERATED FINITE ELEMENTBOUNDARY INTEGRAL METHODS AND HIGHER ORDER ELEMENTS D. Filipovic T. F. Eibert L.S. Andersen J. L. Volakis January 1999

PROJECT INFORMATION PROJECT TITLE: REPORT TITLE: U-M REPORT No.: CONTRACT START DATE: END DATE: DATE: SPONSOR: SPONSOR CONTRACT No.: U-M PRINCIPAL INVESTIGATOR: Algorithms for Phased Arrays Modeling Of Periodic Arrays Using Accelerated Finite ElementBoundary Integral Methods And Higher Order Elements 375458-4-T June 1, 1996 May 30, 1999 February 1, 1999 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/ D. Filipovic, T. Eibert, L.S. Andersen and J. Volakis CONTRIBUTORS TO THIS REPORT.

Contents 1 Introduction 1 2 Infinite Arrays 3 3 Formulation 6 3.1 Finite Element Boundary Integral............................ 6 3.2 Periodic Radiation Conditions......................... 9 3.3 Periodic Boundary Conditions........................... 12 4 Numerical Results 15 5 Conclusions and Future Work 21 6 Appendix I 25 6.1 Extension of the matrix transformation approach for nonequal sideface m eshes...................................... 25 6.2 Generalized PBCs based on a finite-sized FE domain............ 27 6.3 Modified generalized PBCs based on a finite-sized FE domain....... 30 7 Appendix II 31 7.1 On Hierarchical Tangential Vector Finite Elements for Tetrahedra..... 31

List of Figures 1 Important parameters of an infinite periodic array.............. 3 2 Plane wave (a) incident on an array, (b) radiate from an array....... 4 3 Geometry of an infinite double periodic antenna array.............. 6 4 Edge based tetrahedral finite element...................... 8 5 FE mesh consisting of tetrahedral elements.................. 13 6 Infinite periodic (a) slot array, (b) microstrip patch array..............15 7 Case 1: Infinite periodic slot array (E-plane)(a) Input resistance, (b) Input reactance...................................... 16 8 Case 2: Input scanning (a) resistance, (b) reactance............... 17 9 Case 2: Infinite periodic microstrip patch array (a) E-plane, (b) H-plane.. 18 10 Case 3: Infinite periodic microstrip patch array, E-plane scan............. 19 11 Case 3: Mesh (a) HFEM,(b) FE-BI-TETRA.................. 20 12 Waveguide simulator measurement (a) Case 4, (b) Case 5.......... 20

1 Introduction The concept of infinite arrays has been accepted as the most promising for accurate evaluation of phasing properties of large periodic arrays, frequency selective surfaces, diffraction gratings, artificial dielectrics, etc. It is found that the boundary elements in large arrays, although having a different environment than the interior ones, do not have a significant influence on the overall array performance. Full wave analysis of phase arrays is computationally very expensive (memory and time) so an alternative approach must be used. By utilizing Floquet theory [1], an infinite periodic geometry can be accurately modeled and analyzed with a single periodic array cell (to be refered as unit cell) and the contributions of other elements can be included through the appropriate infinite periodic Green's function [2]. Infinite periodic arrays have been successfully modeled in the past using integral equations [3] and finite element methods [4, 5, 6, 7]. In this report, we implement the finite element-boundary integral (FE-BI) method for the analysis of infinite antenna arrays. We particularly emphasize the use of tetrahedral hierarchical mixed-order tangential vector finite elements (TVFEs) [8]. Also, the infinite periodic Green's function accelerated by the Ewald transform [7] is used as the kernel of the radiation integral. Periodic boundary conditions (PBCs) for lower (0.5) and higher (1.5) order tetrahedral finite elements are implemented through the element matrix transformation algorithm [4]. It is found that by properly placing higher order elements in regions where high field variations are expected and lower order elements elsewhere, accurate values for the scanning reflection coefficient can be obtained with relatively coarse discretizations. This is particularly true for scanning angles away from blindness. To improve field modeling around scan blindness, higher order elements should be placed throughout the unit cell. This way, power transfer between neighboring cells is more accurately modeled. Due to the hierarchical nature of the employed vector basis functions, the extension of the PBC algorithm from lower to higher order elements is fairly straightforward, requiring reordering of the nodes along the faces on opposite unit cell boundaries so that continuity of the tangential fields is preserved. The important parameters and properties necessary for understanding an infinite array 1

analysis are briefly reviewed in section 2. Decomposition of the infinite domain down to a unit cell and the appropriate FE-BI formulation are discussed in section 3. Tetrahedral finite elements are used for the unit cell modeling inside the volume and the boundary integral along with appropriate boundary conditions for enclosing the computational domain (section 3.1). Periodic radiation conditions (section 3.2) and periodic boundary conditions (section 3.3) are also formulated. Five different examples necessary for the validation are given in section 4. Conclusions and suggestions for future work are given in section 5. Several ways of implementing periodic boundary conditions in the case of unequal periodic boundary meshes are discussed in appendix I. Finally, several papers important for understanding the developed tetrahedral elements are given in appendix II. 2

2 Infinite Arrays In an infinite array with certain periodicity between the antenna elements, every element is surrounded by the same environment. Thus, in a large array, all elements but those on the boundaries have approximately the same properties. Consequently nearly all elements in a large regular array can be represented by a single periodic element residing in an infinite array. The following parameters are important in evaluating the properties of an infinite array (see Fig. 1): * Scan angle - the angle from broadside at which an array is phased for maximal radiation; * Plane of the scan - defines the plane where steering takes place; * Polarization - polarization of the plane wave that radiates or scatters. intercardinal plane (IC) " El ED I | El E orH - cardinal _1 X l m plane (C) X E D L 1m 1. r --- — i_ broadside radiation direction direction (IC scan) Figure 1: Important parameters of an infinite periodic array. The effect of phase difference used for array steering is shown in Fig. 2(a) for plane wave incidence and scattering and in Fig. 2(b) for radiation. The phase difference on opposite sides of a rectangular unit cell can be calculated from ATX = koDxcos(O)sin(q); A = -koDysin(0)sini(), (1) 3

where D, and Dy are spatial periods in x and y directions and (0, f) are the steering angles. ko / p,I) -- ikO coS() X * sinf) /,I /,' Dy I AT / / T,, Dx Dx Dy k (a) ko,((pe) /I 'D - O S 7 Dy- I / I --- Dx Dx Dy (b) Figure 2: Plane wave (a) incident on an array, (b) radiate from an array. One of the principle problems in designing arrays is impedance matching between the radiation space in front of the array and the internal circuit behind the array. The objective is to reduce the reflection at the interface. A parameter often used to characterize the behavior at the array element input port is the scanning (active) reflection coefficient Psc Zin(f 0) - Zin(f, 0 o. 0) Lscann = (2) Zin (f, 0) + Zien (f, 0 0~)( Here, Zin(f,0) is scanning input impedance for scanning angle 0 in the plane X = const, and Zin(f,0=0~) is the scanning input impedance at broadside. When rscann O then no wave is radiated and the angle(s) at which this occurs is referred to as "scan blindness angle(s)". This corresponds to the physical interpretation that the excited surface wave inside the substrate captures all energy. Another important parameter in the infinite array analysis is the scanning (active) element pattern which can be easily obtained from the scanning reflection coefficient and is given by Pscann = (1 - rscann (f,0) 12)COS0. (3) 4

This parameter can be referred to as the radiating pattern of the infinite array when only one element (unit cell) is fed and all others are terminated with broadside matched loads. 5

3 Formulation 3.1 Finite Element Boundary Integral Consider an infinite double periodic antenna array as depicted in Figure 3. The most rigorous of the implementations is to decompose the infinite periodic array down to a single unit cell and employ the finite element method (with appropriate periodic boundary conditions PBCs for the side walls) to model the unit cell volume. For truncating the finite element mesh on the unit cell aperture (top and/or bottom) the boundary integral provides for a rigorous implementation. The contribution of other cells is included through the Floquet phasing conditions (see (1)) incorporated into the infinite periodic Green's function and the PBCs. / Y 4.....4'....4.... - - - - - - - - - - L — — r / -- -- D - PBC BI or PEC Figure 3: Geometry of an infinite double periodic antenna array. To develop the necessary linear equations, one approach is to begin with the vector wave equation for the electric field V x (ir-'. V x E) - W2~ E -juJ. (4) 6

After applying Galerkin's testing in conjunction with the weighting function T, introducing vector identity V * (A x T) = T. (V x A)- A. (V x T), (5) where A=/ ir-lV x E, and the divergence theorem, we obtain the weak form of the vector wave equation (weak form since the differentiability over E is "weakened" by the derivative over the weighting function T): J {(V x T) (r-1 V x E) - 2 T. = ~E} dV = -|f| T * J dV + i f T (n x H) dS, (6) As usual, E and H denote the electric and magnetic fields, r and =J are the relative tensor permittivity and permeability of the unit cell filling (possibly inhomogeneous), Sa represents the non-metallic portions of the aperture. The volume V, refers to the volume occupied by the impressed sources J. Note that the latter integral refers to H on the antenna aperture Sa and over the feed opening. Also, the unit normal h is directed outward from the boundary surfaces Sa. To solve (6) for E, we require knowledge of H over Sa. In the context of the FE-BI method, the relation between H and E is determined by the boundary integral equation H = -2ikoYo Jj G(r, r') (z x E(r')) dS' (7) where G is the half space dyadic Green's function G -I - VV Gp(r,) (8) with I is the unit dyad. Here Gp(r, r') is a free space infinite periodic Green's function (see 3.2). To construct a linear set of equations from (6) and (7), we must first tessellate the volume and introduce expansions for each of the tessellation elements. We chose to work with hierarchical mixed order TVFE for tetrahedra [8]. Beside high flexibility in modeling (owing to tetrahedral shape) the hierarchical nature of these elements allow the introduction of higher order expansions over a chosen set of elements. 7

Figure 4: Edge based tetrahedral finite element. Choosing tetrahedrals as the tessellation elements (see Fig. 4), the field is expanded within the unit cell as 20 E EJeW = [W]{Ee} j=1 (9) where [W]e = [{W}, {Wy}, {WW}] and {Ee} = {Ee, Eq,... E0})T. 20 basis functions are needed for higher order expansion: twelve are associated with edges (six provide constant tangential components at edges and six model the corrections from constant values of the tangential fields), and eight are face based expansions. These and other higher order basis functions available in literature are given in Tab. 1 [18]. TVFE Vector basis functions Order Label Edge-based Face-based 0.5 Whitney (iV( - (jV( 1.5 Peterson and (iVCj ck(iiVj - jVi) Savage (jV(i j ((kV - (iVck) 1.5 Graglia (3Ci - 1)((iV(j - CjV() (k((iV(j - (jV(i) et.al. (3j - 1)(iCVj - (jV(i ) (j((kV(i - (ioVk) 1.5 Andersen and ( -(jV(i (k((iV(j -(jV(i) Volakis ((i-( ((iV(j- -(jV(i) j((kV(i - (iVk) 1.5 Webb and iV(o - jV(i (k((iV - j ) Forghani (iV(j + jV(i (j((kVi - (iV(k) 1.5 Lee and (iV(j (k(iV( Cendes _jVi (kCjV(i Table 1: Vector basis functions for six different TVFEs; (i, j, k) follow cyclic order, with each index ranging form 1 to 4. 8

On the unit cell aperture, the electric field is reduced to 3 ES(r) E>ES'(r)= [S]'{Es} (10) i=i where [S], = [S, Sy] and Si = 2Ae x (r - ri) (11) Here, r and ri refer to the position vector within the triangle and at the ith node of the triangular face on the aperture. The parameters li and Ae denote the length of the ith edge of the triangle and its area, respectively. To generate a linear system for E, (9) and (10) are substituted into (6) and (7) and Galerkin's method (setting T = W) is employed to yield the assembled system [A] { {EV [0 ] 1[ {{E} = { {} J (12) In this system, {EV} denotes the field unknowns within the unit cell volume enclosed by SO, whereas {ES} represents the corresponding unknowns on the boundary Sa. The excitation column {bv} are due to internal sources. 3.2 Periodic Radiation Conditions The edge-based basis functions for tetrahedra reduce to Rao-Wilton-Glisson basis functions [9] on the unit cell aperture. Therefore, these basis functions are used to represent the field in the boundary integral in (3). In the spatial domain, the periodic Green's function Gp(r, r') has the form Gp(r, r')= E e-ktPmn e Rmn (13) Z ~-j~~m-Rmn m=-oo n=-oo 43mn where Rmn-= r-r -Pmn. (14) In the spectral domain, Gp(r, r') becomes 1 0 o 1, Gp(r, r') = 2 jk tmn.(p-p )e-j - (15) m=-oon=-oo J zmn 9

where A = IPa x Pbl is the cross sectional area of the unit cell, r = p+zz, (16) 27 ktmn = ktoo +. [m(pb X z) + n(z x Pa)] (17) is the so-called reciprocal lattice vector, and kzmn - - 2 - ktmn * ktmn, (18) where Re(kzmn) > 0,Im(kzmn) < 0. In many cases, the spectral domain representation (15) has satisfactory convergence behavior if applied in a spectral-domain formulation of the integral equation. However, for arbitrary array configurations analyzed in the space domain, having strongly as well as weakly coupled array elements, it is necessary to have a representation that converges faster than either (13) or (15). This can be achieved by employing the so-called Ewald transformation [12, 13]. The Ewald transformation starts from the spatial domain representation of the periodic Green's function (13) and makes use of the identity e-jkoRmn 2 00 o2 k2 =- _-| eRTn ds, (19) Rmn Va where s is a complex variable. In order that the integrand converges as s -+ 0 for a wavenumber ko with an arbitrary amount of loss, the path is chosen so that arg(s) = 4 as s - O. To have convergence as s -+ oo, the path is chosen so that -7r/4 < arg(s) < 7r/4. Next, (19) is substituted into (13) and the parameter E is introduced to split the integral into two terms as Gp(r, r') = Gp1(r, r') + Gp2(r, r') (20) where 1 00 00 2 2 k2 Gp(1jkto 2 fE Rm S2+T- ds (21) 4rm:-oo n=-oo V/ Jo 1 00 00 2 00 2 2 Gp2(r,r') - 4 e-jktoo-P,, / -RmnS ds (22) 4,F m=-oo n=-oo V7 JE Using the identity [14] (eq. 7.4.34) 2 / e-Rm'n2+ ds = [eomn erfc (RmnE- 2 VT 2R mn 2E( + ejkoRmn erfc RmnE + j-, (23) \27 10

where erfc is the complementary error function, Gp2(r, r') can be written as 00 0o e-jktoo- pmn k \ Gp2(r, rl) m=-o _ ____ [e- erfc - jkm m=-oo n=-oo 87rRm e mn 2E + ejkoRmn erfc RmhE +, (24) ( 2E (24) which is essentially a "modified" spatial domain portion of the periodic Green's function. Making use of the Poisson transformation, or alternatively following the procedure in [13, 15], employing a transformation formula for the series expansion of the t9-function, (21) is finally transformed to 1 km 1 Gp (r,r') - A e-jktmn (Pkzmn m=-oo n=-oo 4jzmn [e-jkzmnZ-' erfc (jk _ Iz - zE) ( \2E) + eikzmnlz-zI erfc (J + z - z'IE)], (25) where A is defined immediately after (15)). Equation (25) can be identified as a "modified" spectral domain portion of the periodic Green's function. For planar BI surfaces, we can select z = z, = 0, giving the simplified form 1 00 00 1 ) erfc kzmn Gpl(r,r') = E ES -2jk p-) jkzmne- P (26) m=-oo n=-oo 2jkzmn The two expressions (24) and (25) or (26) both converge exponentially (Gaussian convergence) and their computation is therefore very efficient requiring only a few terms of the series. The parameter E controls the convergence rate. As E becomes larger, the spatial series (24) converges faster, but (26) converges slower. The optimum parameter is that which makes the two series converge at the same rate, so that equal numbers of terms are required in the calculation of both series (assuming the same calculation time for corresponding terms in each of the two series). By analysis of the asymptotic behavior of the series terms, the optimum parameter Eopt is found to be [15] Eopt A (27) Choosing this value for E and adjusting the summation limits so that the most dominant terms are kept, in almost all practical cases it is sufficient to include only 9 summation 11

terms in (24) and (26) (i.e., the summation limits are from -1 to +1). With these few terms, the error level is usually less than 0.1%. For the implementation of the BI portion of the method, we apply the same phase transformations to the matrix elements associated with edges on Fri and Fup (see Fig. 5) as done in the FE portion of the implementation (see section 3.3). For our approach, source and test triangles are always inside the unit cell, therefore guaranteeing that the singularities of the neighboring array cells are never inside the test triangle. However, it is still necessary to carefully deal with the singularities of the neighboring array elements that are close to the test triangles. 3.3 Periodic Boundary Conditions The derivation of PBCs for unit cell meshes with equal surface meshes on opposite side faces, starts from a FE mesh of the whole infinite periodic structure. Thus, the bounding surface S of the FE solution domain only consists of the top and bottom termination surfaces of the mesh. The PBCs for the vertical side faces of the finite mesh of one unit cell are derived in a way that each unknown of the unit cell mesh sees an environment as it would see in the infinite mesh. If we consider for example a unit cell mesh as illustrated in Fig. 5, typically all unknowns on the right and upper side faces, Fri and Fup, are eliminated and matrix elements involving these unknowns are transformed to the corresponding unknowns on the left and lower side faces, Fle and Fo. Thus, the unit cell mesh is wrapped around and connected at opposite side faces so that no vertical boundaries are present anymore. The relation between opposite unknowns is found from the periodicity condition E(r + m Qa = n Qb) E(r) e-jktoo(m Qa+n Qb) H(r + m Qa + b) = H(r) e-iktoo(mQa+nQb) (28) with ktoo - kxoo x + kyoo y ko sin 0o cos po x + ko sin 0o sin po Y. (29) 12

r m l/e '[x m ri r nlo lo Figure 5: FE mesh consisting of tetrahedral elements. Here, ko is the wave number of free-space and 00, Po describe the scan direction of the array. If ek is the unknown field at an edge on one of the vertical side faces, the value ek, of the field at the corresponding edge on the opposite side face is given by ek, = eke-jkto' (30) where Ar = (~Pa or i Pb) is the vector joining the two edges (see Fig. 3 for definitions of Pa and Pb). Special treatment is needed for the unknowns on the edges at the upper right corner of the unit cell. These unknowns need to be transformed two times so that they finally relate to the unknowns corresponding to the edges in the left lower corner of the unit cell mesh. To obtain the correct phase relation, valid for this two stage transformation process, Ar must be set to ~(pa + Pb)The necessary transformation algorithm for the matrix elements Akl encountered in the matrix generation process can be summarized as Aki - Ak'l' = Akle-jktoo~Arkk e+jktoo-Arll (31) and means that a matrix element Akl is multiplied by a phase factor and creates a new element Akl/I. Basically, this equation can be applied to all matrix elements. However, it is only necessary if either edge k or edge I is on the right or upper vertical side walls. Otherwise the image edge is the same edge as the original one and no transformation 13

needs to be carried out. Important to note is that Ar is oriented from the original edge (on right or upper side wall) to the image edge (primed). The phase factor for the column index I in (31) is directly obtained from (30) by writing the original unknown in terms of the image unknown. However, for the row index transformation, it must be kept in mind that within the weak formulation of the FE problem, the basis function related to the row index is the test function. Therefore, in order to perform testing with the image unknown instead of the original unknown, the phase transformation must be carried out the other way round. This means that the original unknown is multiplied with the corresponding phase factor and replaced by the image unknown, as directly implied by (30). Due to the nature of the implemented higher order elements, three types of matrix elements are encountered for the matrix element transformation algorithm. If the global directions of parallel edges on opposite boundaries are opposite, then to preserve the continuity of the fields (in infinite environment), the transformed matrix element (associated with the Whitney basis function (iV(j-(jV(- ) should be multiplied with -1. This is not the case for the other edge-based basis function (((i - (j)((V(j -(jV(j)) since the tangential component of the field is always directed from the edge end to the edge center. To assure proper transformation of the matrix elements associated with a face-based basis function ((k((iV(j-(jV(-)) it is necessary that global numbering for the parallel faces on opposite sides be the same. In that case, the continuity of the electric fields is also preserved (see appendix II). 14

4 Numerical Results To examine the use of tetrahedral hierarchical mixed-order TVFE for modeling infinite periodic arrays and validate presented formulation, three test cases are analyzed: * IPSA - infinite periodic slot array; * IPMPA - infinite periodic microstrip patch array; * WGSMP - waveguide simulator for microstrip patch array. The geometries are given in Fig. 6 and appropriate design parameters are given in Tab. 2. Dx Px,w Dx Px I Dy PY4 Dy Fx d | (a) Figure 6: Infinite periodic (a) ] (b) slot array, (b) microstrip patch array. First, consider the infinite periodic, horizontally fed, slot array given in Fig. 6(a) (operating wavelength Ao=lOcm). The reference results are found using the FE-BI method with prismatic finite elements (a code developed at the University of Michigan) [7]. These results will be referred to as FE-BI-PRISM. To achieve convergence, Ao/40 discretization has been applied. The mesh discretization for the FE-BI method with hierarchical mixedorder tetrahedral elements (referred to as FE-BI-TETRA) is Ao/20. Since the slot is very 15

Example Parameters defined in Fig. 6(a),(b) Case Problem Dx Dy Px Py Fx Fy d ei1 e2 1 IPSA 5.5 5.5 5 0.5 2.5 - 0.5 - 1 2 IPMPA 5 5 3 3 0.75 1.5 0.2 1 2.55 3 IPMPA 5 5 1.5 0.98 0 0.75 0.6 - 12.8 4 WGSMP 2.37 2.22 1.8 1.8 0.9 1.32 0.159 - 2.33 5 WGSMP 2.378 2.215 1.8 1.8 0.9 1.32 0.242 - 2.58 Table 2: Parameters for testing cases. (Unit for length is cm.) narrow, higher field values can be expected in the proximity of the slot. Higher order elements are therefore placed in the light gray area in Fig. 6(a). The results for the scanning input impedance in the E-plane (0-900) are given in Fig. 7. While there are differences in results when only lower order elements are used (HO), introduction of higher order elements (HO/HI) leads to significant improvements. A similar trend is obtained in the H-plane. z C: Z z O IU) () Iz Z 0 Z SCANNING ANGLE IN DEG (a) 20 30 40 50 60 70 80 SCANNING ANGLE IN DEG (b) Figure 7: Case 1: Infinite periodic slot array (E-plane)(a) Input resistance, (b) Input reactance. Next, consider the vertically fed microstrip patch array on a thin substrate shown in Fig. 6(b) (Ao=10cm). In addition to the FE-BI-PRISM results, method of moments (MoM) reference results [3] are also used for comparison. The computational cost, array broadside resonance and broadside input impedance (Zbroad) at the analysis frequency 3GHz for the examined cases (HO-0.5 elements, HO/H1-0.5/1.5 elements in the light gray 16

area with linear PBCs, and HI-1.5 elements with higher order PBCs) are given in Tab. 4. The unit cell mesh discretization is Ao/40 for FE-BI-PRISM and Ao/28 for FE-BI-TETRA (HO, HO/HI, HI). For better field modeling on the substrate surface, an extra layer of air is added (see Fig. 6(b)). The expected higher field variations at the surface are now modeled with higher order vector basis functions. Results for broadside scanning impedance and scanning reflection coefficient in the E-plane (q$=0~) and H-plane (0=90~) are given in Fig. 8 and Fig. 9 respectively. Method Broadside resonance Zbroad Number of Time per angle [GHz] [Q] unknowns [min] FE-BI-PRISM 2.903 43.8-i18.7 2744 - FE-BI-TETRA (HO) 2.725 3.93+i3.7 2542 22.1 FE-BI-TETRA (HO/HI) 2.915 24.1-izl.03 10360 34.1 FE-BI-TETRA (HI) 2.913 24.1-i11.4 13127 56.1 Table 3: Infinite periodic microstrip patch array. 2.5 2.6 2.7 2.8 2.9 3 3.1 -. -,. -",.4 -.5 2. 2.7 z. FREQUENCY IN GHz FREQUENCY IN GHz (a) (b) Figure 8: Case 2: Input scanning (a) resistance, (b) reactance. The curves in Fig. 8 display a clear shift in a resonance frequency predected by the solutions based on the HO and HI elements. This clearly shows that lower order elements can not accurately predict scanning properties of an array. This is particularly revealed 17

in Fig. 9. When the array is phased for scanning in the H-plane, the FE-BI-TETRA results converge to the FE-BI-PRISM values provided higher order elements are introduced throughout the unit cell. In the E-plane, scan blindness cannot be perfectly reached, but the maximum scanning reflection coefficient is shifted to 83.40 (85.6~ in [3]) and the resulting curve closely follows the reference. Since the broadside resonant frequency is shifted toward lower frequencies, i.e. the broadside resonant wavelength is larger than the operating wavelength (Ares > Ao as compared to AresXo in [3]), scan blindness occurs at lower scanning angles, as observed. To obtain better agreement with MoM, a denser mesh is needed with higher order elements throughout the unit cell at scanning angles near grazing. As expected (see Fig. 9), a very small difference for broadside resonance and Zbroad between the HO/HI and HI cases is observed (area around radiated edges is modeled with higher order elements). Also, the total CPU time for the H1 case is increased as compared to the HO/HI case, not only due to the larger number of unknowns but also due to the deteriorated convergence when higher order PBCs are applied., '-' —'-'' - ' 1 ' ': ' - ' ' '-'.. o 0o 0.9 MoM (Pozar) 0.9 MoM (Pozar) - MoM (Pozar) ~- ~9-! — FE-BI-PRISM /' ---- FE-BI-PRISM ---- FE-BI-TETRA (HO) m 0.8 FE-BI.TETRA z. FE-BI-TETRA (HO) —. E-BI-TETRA (HO) '/ - ru i x 0 x A /0.8 FE-BI-TETRA (HO/H1) / 0 o o FE-BI-TETRA (HO/H1) FE-BI-TETRA (H1) - e:. S N ALG 0S.7 ( 0.6 - '/ 0.6 o, 0 diguri, o subsea (casepr 3 min Ta.2 paretchnarred (a) Fr, a coplario with our method, referenced MoM results [3] and hybrid FE-BI results from [4] obtained 0. W LL 0.4 <3 0.3 (/3 Z 0.2 Z 0.2 0.1 0.1 0 10 20 30 40 50 60 70 80 90 0 10 20 30 40 50 60 70 80 SCANNING ANGLE IN DEG SCANNING ANGLE IN DEG (a) (b) Figure 9: Case 2: Infinite periodic microstrip patch array (a) E-plane, (b) H-plane. Infinite periodic arrays, of vertically fed microstrip patch antennas on a thick high dielectric constant substrate (case 3 in Tab. 2) are considered next. For a comparison with our method, referenced MoM results [3] and hybrid FE-BI results from [4] obtained with linear tetrahedral elements (refered as the HFEM in Fig. 10) are also included. The 18

computational costs for the two finite element based methods are given in Tab. 4, and scanning reflection coefficient results are given in Fig. 10 (b). Also, the 3D meshes used for both methods are displayed in Fig. 11. Due to the nature of the problem (thick dielectric with high dielectric constant) the H1 simulation did not converge. Method Number of Number of Solve time / Matrix size mesh cells mesh edges fill time [MBy] HFEM 4291 4862 4.45 6.2 FE-BI-TETRA (HO/H1) 2352 2838 2.4 12.8 Table 4: Computional costs for HFEM [4] and FE-BI-TETRA. 1 --- ~7! — z o.9- - 3 0.7 U. o 0.6 z O 0.5 cc 0. FE-BI-TETRA(HH1 z z < 0.2 - o MoM (Pozar) a HFEM (McGrath) 0.1 0 o FE-BI-TETRA (HO/Hi] 0 10 20 30 40 50 60 70 80 90 SCANNING ANGLE IN DEG Figure 10: Case 3: Infinite periodic microstrip patch array, E-plane scan. Finally, two examples of a waveguide simulator (case 4 and case 5 in Tab. 2) are studied and results are given in Fig. 12. The results are in terms of the clasically defined reflection coeficient (Zo=50Q). To achieve better agreement with the referenced results ( [3] in Fig. 12(a) and [11] in Fig. 12(b)) better feed modeling is required. However, the accuracy improvement when H1 elements are used is obvious. 19

/':' 4Zc~ ~ /i:i.ii: e i - *:. (a) (b) Figure 11: Case 3: Mesh (a) HFEM,(b) FE-BI-TETRA. I. 0.9 0.8 0.7 z 0 M 0.6 o. 0 0.5 Z 0 F m 0.4 le 0.3 0.2 0.1 - MoM (Pozar) o Meas. (Pozar)..... FE-BI-TETRA (H1) FE-BI-TETRA (HO) 0 0," O \ /,O '>s? 0 SA I u 4.4 4.6 4.8 5 5.2 FREQUENCY IN GHz 5.4 5.6 (a) (b) Figure 12: Waveguide simulator measurement (a) Case 4, (b) Case 5. 20

5 Conclusions and Future Work A hybrid FE-BI method for infinite periodic array analysis is developed. The method utilizes tetrahedral hierarchical mixed-order TVFEs for modeling the unit cell volume. Radiation boundary conditions are implemented using the standard boundary integral with the appropriate free space periodic Green's function. Due to the slow convergence of the double infinite sum, Ewald's transform is used for accelerating the computation of the periodic Green's function. Periodic boundary condition based on the matrix transformation algorithm are used for bounding side walls of the unit cell, and its implementation is carried out in both the finite element and boundary integral part. These PBCs are developed for linear and higher order tetrahedral TVFEs. For validation, several examples are given and the general agreement is very good. It is shown that by using hierarchical basis functions and combining lowest and higher order elements, the field modeling is dramatically improved and the array scanning properties are therefore more accurately evaluated. To predict the scanning angle at which blindness occurs, higher order elements must be used throughout the unit cell. This is very important for phased arrays where the energy is trapped into the excited substrate surface waves. Finally, a detailed discussion on various aspects of the periodic boundary conditions and observations on their robustness is included. Several methods were implemented and obtained results were not satisfactory. However, improved ways for attacking the problem are proposed and these should be thoroughly investigated and implemented in the near future. Beside that, the following topics should be investigated in the future: * Implementation of the PBCs discussed in Appendix I for unequal meshes on periodic boundaries. This type of implementation has been carried out for simpler prismatic elements and is quite promising. Its implementation for tetrahedral elements is more complex and carries more risks. Nevertheless, the success of this implementation will lift the serious restriction of having identical meshes (image-like meshes) on the opposite periodic boundaries. In any case, it should be noted that drastically unequal meshes between periodic boundaries are not easy to handle and would lead 21

much worse system condition. * Inclusion of the option of LU and iterative solution and comparison the two in terms of time and memory. * Inclusion of the option of using Lee-Cendes HI elements vs. the Andersen-Volakis elements. This involves a straightforward transformation and will allow for a simpler inclusion of the periodic array code as a module into the HFSS. * Implementation of the second order ABC truncation scheme with PBCs to conform with the current methodologies used in HFSS. Note that ABCs will introduce inaccuracies and given our experience in computing the scanning reflection coefficient, the resulting numerical system may be inadequate for predicting this necessary array parameter. * Introduction of fast integral algorithms into the present FE-BI-TETRA version for low memory and much faster computation. The fast spectral domain algorithm (FSDA), adaptive integral method (AIM) and fast multipole method (FMM) are already available in other University of Michigan codes. This code will be a faster and more rigorous module for reliable array analysis. 22

References [1] N. Amitay, V. Galindo, C. P. Wu, Theory and Analysis of Phased Array Antennas, Bell Telephone Lab., 1972. [2] A. Ishimary,Electromagnetic Wave Propagation, Radiation and Scattering, Prentice Hall, 1991. [3] D. M. Pozar, 'Analysis of an infinite array of rectangular microstrip patches with idealized probe feeds', IEEE Trans. Antenna Propagat., pp. 1101-1107, October 1984. [4] D. T. McGrath, Phased Array Antenna Analysis Using Hybrid Finite Element Methods, PhD Dissertation, Air Force Institute of Technology, Wright Patterson Air Force Base, Ohio 1993. [5] E. W. Lucas and T. W. Fontana, 'A 3-D hybrid finite element/boundary element method for the unified radiation and scattering analysis of general infinite periodic arrays', IEEE Trans. Antennas Propagat., pp. 145-153, February 1995. [6] J. D'Angelo and I. Mayergoyz, 'Phased array antenna analysis', in Finite Element Software for Microwave Engineering, Edited by T. Itoh, G. Pelosi, and P. Silvester, John Wiley and Sons, Inc., 1996. [7] T. F. Eibert, J. L. Volakis, D. R. Wilton, and D. R. Jackson, 'Hybrid FE/BI modeling of 3D doubly periodic structures utilizing triangular prismatic elements and MPIE formulation accelerated by the Ewald transformation', IEEE Trans. Antennas Propagat., accepted for publication. [8] L. S. Andersen and J. L. Volakis, 'Hierarchical tangential vector finite elements for tetrahedra', IEEE Microwave and Guided Wave Letters, vol. 8, pp. 127-129, March 1998. [9] S.M. Rao, D.R. Wilton, A.W. Glisson, "Electromagnetic Scattering by Surfaces of Arbitrary Shape," IEEE Trans. Antennas Propagat., vol. 30, no. 3, pp. 409-418, May 1982. 23

[10] D.S. Filipovic, J.L. Volakis and L.S. Andersen, 'Efficient modeling and analysis of infinite periodic antenna arrays by tetrahedral finite elements', paper submited to IEEE AP Conf., Orlando, FL, 1999. [11] M.A. Gonzales, at all.,'Full-wave analysis of cavity-backed and probe-fed microstrip patch arrays by a hybrid mode-matching generalized scattering matrix and finiteelement method', IEEE Trans. Antennas Propagat., pp. 234-242, February 1998. [12] P.P. Ewald, Dispersion und Doppelbrechung von Elektronengittern (Kristallen), Dissertation, Miinchen, 1912, (also Ann. Phys. 49, p. 1, 1916). [13] P.P. Ewald, 'Die Berechnung optischer und elektrostatischer Gitterpotentiale,' Ann. Phys. 64, pp. 253-287, 1921. [14] M.Abramowitz, I.A.Stegun, Handbook of Mathematical Functions, Dover Publications, New York, 1965. [15] K. E. Jordan, G. R. Richter, P. Sheng, "An Efficient Numerical Evaluation of the Green's Function for the Helmholtz Operator on Periodic Structures," J. of Comp. Phy., 63, pp. 222-235, 1986. [16] L.S. Andersen and J.L. Volakis, 'Development and application of a novel class of hierarchical tangential vector finite elements for electromagnetics', to appear in IEEE Trans. Antennas Propagat., January 1999. [17] L.S. Andersen and J.L. Volakis, 'Accurate and efficient simulation of antennas using hierarchical mixed-order tangential vector finite elements for tetrahedra', submited to IEEE Trans. Antennas Propagat. [18] L.S. Andersen and J.L. Volakis, 'Condition numbers for various FEM matrices', submited to IEEE AP Conf., Orlando, FL, 1999. 24

6 Appendix I 6.1 Extension of the matrix transformation approach for nonequal sideface meshes In the case of unequal surface meshes an opposite sidefaces, the individual unknowns corresponding to the edges in the right and upper vertical sidewalls cannot be directly related to one unique image unknown. However, the matrix transformation algorithm can formally be extended to account for this case, if it is assumed that each unknown is now related to several image unknowns. Thus, by applying the transformation, each transformed matrix element results in several new matrix elements. This generalized transformation can be written as Ak - [AkI] Ak [f(k, 1, k, 1)], (32) where i - 1,..., Ik, j 1,..., J1. Therefore, the unknowns k and I are related to Ik and Jk image unknowns, respectively and [Aij] denotes a matrix with elements aij. f (k, 1, k, 1) is a complex coefficient which represents the mapping of the original unknowns on the image unknowns and also contains the phase informationa in accordance to the scan direction of the array. At this point, the formal transformation of the matrix elements is defined. However, the mapping of the unknowns on the right and upper vertical walls on their image unknowns on the left and lower sidewalls still needs to be specified. This is equivalent to the determination of the coefficients f(k,, k, l,). In principle, there should be several ways of mapping the unknowns. However, it is obvious that each mapping procedure must resemble to the one by one relation for equal sideface meshes. Also, it cannot be expected that the resulting algorithm will give good results for strongly unequal meshes since the derivation is kind of heuristic and there is obviously no theoretical foundation for the whole procedure. In the following, a mapping algorithm is presented which is based on a weak enforcement of the periodicity conditions for the electric field intensity. The mapping relation must allow to calculate all unknowns related to the right and upper vertical side walls from the unknowns related to the left and lower vertical sidewalls. Unknowns in the cor 25

ner of the unit cell mesh belonging to a right/upper sidewall and a left/lower sidewall are considered to belong to the right/upper sidewall only and need to be eliminated. Using (28), the fields in the opposite vertical sidewalls can be related by (as discussed earlier) Eri(U, v) E= e(, v) e-ktooPa Eup(s, t) = Elo(s, t) e-joo-Pb, (33) where (u, v) and (s, t) are appropriate parametrizations of the left/right and lower/upper sidewalls, respectively. Introducing the FE discretization of the unit cell, the relations become Nri Nle E Eni fni (U, v) = 5 Ele fnle (, v) e-jktooPa, nri-1l ne=l1 Nup Nlo Enpfnup ( t) Enlo f nl (s, t) e-jktooPb, (34) nup=1l no~l where f denotes the tangential components of the basis functions refered to the vertical sidefaces. Multiplying the equations with weighting functions fmri and fmup and integrating over the right/left and upper/lower sidefaces, respectively, gives Nri Nte Enri dmmi E(f, Vf v)f ( U, v) dudv e-jktoo~~' nri=l (uv) nle=l (UV) mri 1,..., Nri I Nup N10o Z t)ft)dsdt = E jk))o- Pb nup-1 (s,t) nlO-l (st) mup = 1,..., Nup (35) as a weak enforcement of the periodicity condtions. The two equation systems in (35) are coupled via the unknowns belonging to both the right and upper sidefaces. Therefore, the two equation systems are combined to one equation system. This is achieved by introducing new numberings for the unknowns in the sidefaces. All unknowns in the right and upper sidefaces get a unique index nriup and the remaining unknowns in the left and lower sidefaces, which do not belong to the right or upper sideface at the same time, get a new index nielo. Now, the equations in 26

(35) can be rewritten in a way that the Nriup unknowns with indices nriup and the Nlelo unknowns with indices nlelo are collected on the left-hand and right-hand sides of the resulting equation system, respectively. This equation system can be written in matrix form as [Amriupnriup] {Enriup = [Amriupnlelo] Enlelo} (36) Solving this equation system for the Enriup according to {EJnriup} = [Amriupnriup] 1 [Amriupnlelo] {Enieio}, (37) gives the necessary relations for the generation of (31). As already mentioned, the above extended matrix transformation algorithm was implemented in a test code utilizing prismatic elements. Initial results imply that the algorithm works for slightly different surface meshes on opposite sidefaces, however, strong instabilities were encountered for larger differences. 6.2 Generalized PBCs based on a finite-sized FE domain The matrix transformation algorithm was derived under the assumption that the finite sized unit cell mesh behaves like an equivalent infinite mesh of the whole periodic structure. This strategy is obviously strongly dependent on the mesh structure and its extension to unequal surface meshes on opposite sidefaces of the unit cell is not straightforward. A theoretically correct statement of the problem can be obtained by deriving the FE formualation for a finite-sized FE domain, which is the unit cell under consideration. This requires that the vertical sidewalls of the unit cell are considered to be part of the bounding surface S in (6). In the FE implementation, this is achieved by evaluating the surface integral term in (6) also for the vertical sidewalls. For this purpose, an expansion of the tangential components of the magnetic field intensity H in the vertical sidewalls is introduced. In the following, it is assumed that H is discretized using the same set of basis functions as used for E and that the numbers of basis functions with nonzero expansion coefficients in the individual sidefaces are the same for H and E. 27

Based on this assumption, the H fields in the four sidefaces can be written as Ne1 Hle(u,V) = Hnlefnl(uU), nle=l Nro Hri(u, ) = Hnifnri(u) nri=l Nlo=l H10(s, t) E HnlO fno(sI t) n10=1 Nup Hup(S,t) =- Hnupfnup(S,t), (38) nup=l where (u, v) and (s, t) are defined as in (33). Now, the surface integral terms in (6) can be evaluated and assembled into the global FE equaiton system. For this purpose, the global numbering scheme, which typically only involves the E field unknowns in the FE solution domain, must be extended so that it also accounts for the H field unknowns on the vertical sidewalls. Nevertheless, the resulting FE equation system has less equations than unknowns and can, therefore, not be solved. However, the remaining equations necessary to form a uniquely defined system can be obtained from the periodicity conditions for the the E and H fields in the vertical sidewalls. The necessary relations for the fields in the lef/right and lower/upper sidewalls are derived from (28) and result in Eri(u, v) = Ele(U, v)e-jktooPa~ Eup(s, t) = El1(s, t)e-jkt0~Pb, (39) Hri(u, v) = Hle(U, ) e-jkt0ooa, Hup(S, t) = Hlo(s, t) e-jkt~Pb. (40) Introducing the expansions for the fields gives Nri NIe Enri f (, v) E nfn v (u, ) e-kta, nri=l1 nle=l Nup Nlo Enupf nup (S) t) Enlofnl (s, t) e -jktoo'Pb (41) nup=1 n1o=l Nri NMe E Hnri fnri (u ) = 5 Hnlefnl (u, v) e -kto0Pa nri=l nle=l 28

Nup Nlo E Hnupfnp (s, t) = HnOfnl (s, t) e-jktooPb. (42) nup=l nlo=l For identical surface meshes, these equations can be fulfilled in a one-by-one relation for the individual expansion coefficients, however, for the general case of arbitrary surface meshes on the opposite sidefaces, a weak enforcement of the equations is carried out. Therefore, the equations are multiplied by appropriate weighting functions and integrated over the domains of the parametrizations (u, v) and (s, t). The resulting linear algebraic equations are Nri E Eni JJfmi (u,)fnri(u ) dudv = nrT=1 (u,v) Nup E Enufp //f no (s, t) f nup(s, t) dsdt nup=1 (s,t) Nri E HnriJf fmi(U V)fnri(Uv) dudv nri=l (u,v) Nup E Hnup f Mup (s, t)f (st) (sdsdt nup=l (s,t) Nie E En, f l,,f (u, v) f (u, v) dudv e-jktoo-P nle=1 (u,v) mle =,..., Nle Nlo E Eno10 JfM, O f1 n (s, t) dsdt e-jkto~Pb, n1i=1 (s,t) mIo = 1,..., No, (43) Nie > Hnle ]f,ri (, V)f ( ) dudv e-jktoPa, nle=1 (u,v) mri - 1..., Nri Nlo E Hlo /fmup (s, t)f n, (s, t) dsdt e-jkto~~Pb nlo=l (s,t) mup l,..., Nup. (44) In contrast to the matrix transformation method described in section 6.1, these equation can be assembled into the global FE equation system to obtain a final discretized formulation of the periodic FE problem. The global numbering scheme, which typically only involves the E field unknowns in the FE solution domain, needed already to be extended to account for the H fields in the surface integral contributions of (6). All equations related to global row indices of the E field expansion functions are obtained from the discretization of the functional (6). Therefore, the equations in (43) and (44) are considered to be global equations corresponding to the global row indices of the H field expansion functions. Since the unknowns related to the vertical edges in the four corners of the unit cell mesh belong to two vertical sidewalls, the corresponding two equations in (43) 29

and (44) add up to one equation in the global system. Also, it should be noted that the left/lower basis functions are used for testing in (43) whereas the right/upper basis functions are used in (44). This guarantees that the correct number of equations equal to the number of H field unknowns is generated. The presented generalized PBC algorithm was implemented in a FE code based on prismatic elements. First results are very promising, however, severe problems with the condition of the final equation system were encountered. It is necessary to solve the equation system by singular value decomposition since several eigenvalues are extremely small. It is assumed that these problems arise from an unproper handling of the connection between the BI surfaces and the vertical sidewalls. The explicit introduction of the H fields on the vertical sidewalls is obviously not compatible with the H field free formulation of the BI termination. Therefore, it is suggested to investigate this problem further and try to find a consistent formulation. Another issue may arise with the handling of metallic sections in the vertical sidewalls for which the electric field unknowns must be set to zero. At the moment, it is not completely clear what this would mean for the magnetic field unknowns and how a consistent system of equations can be obtained. This points needs also further investigation. 6.3 Modified generalized PBCs based on a finite-sized FE domain The algorithm presented in section 6.2 contains redundancy with regard to the H field expansion. Basically, there is no need to discretize the H fields in both sidewalls which are opposite to each other. The evaluation of the surface integral term in (6) can also be carried out if only one discretization is introduced. In this case, the periodicity conditions for the H fields are enforced in a strong way according to (40). This implementation is, of course, more efficient than the one presented in section 6.2, however, the discussed issues, which need to be solved, still remain. 30

7 Appendix II 7.1 On Hierarchical Tangential Vector Finite Elements for Tetrahedra This appendix contains four papers written by Lars S. Andersen and John L. Volakis on a new class of hierarchical mixed order TVFEs for tetrahedra. These TVFEs are also used in developing the formulation for the periodic array analysis. These are: * Development and application of a novel class of hierarchical tangential vector finite elements for electromagnetics - this paper proposes a new class of hierarchical mixed-order TVFE for solving various boundary value problems [16]. * Hierarchical tangential vector finite elements for tetrahedra - higher order basis function for tetrahedra are presented in this paper as a natural extension of 2D basis function for a a triangle [8]. * Accurate and efficient simulation of antennas using hierarchical mixedorder tangential vector finite elements for tetrahedra - the efficiency of antenna modeling when both lower and higher order elements are used is demonstrated in this paper [17]. Conditioitn numbers for various FEM matrices - the analysis of the system condition arising when different interpolatory and hierarchical TVFEs are used within the contex of finite element is provided in this paper [18]. 31

Development and Application of a Novel Class of Hierarchical Tangential Vector Finite Elements for Electromagnetics Lars S. Andersen and John L. Volakis Radiation Laboratory, Department of Electrical Engineering and Computer Science, University of Michigan, Ann Arbor, MI 48109-2122, USA Abstract Tangential vector finite elements (TVFEs) overcome most of the shortcomings of node based finite elements for electromagnetic simulations. For a triangular element, this paper proposes a class of hierarchical TVFEs that differ from traditional TVFEs. The hierarchical nature of the proposed TVFEs makes them ideally suited for employing an efficient selective field expansion (the lowest order TVFE employed within part of the computational domain and a higher order TVFE employed within the remaining part of the computational domain). This is an attractive feature not shared by non-hierarchical TVFEs for which a more traditional approach (the same TVFE employed throughout the computational domain) must be applied. For determining the scattering by composite cylinders, this paper argues that the performance (in terms of accuracy, memory and, in most cases, CPU time) of the proposed class of hierarchical TVFEs when applying selective field expansion is superior to that of the lowest order TVFE and a traditional mixed-order TVFE. This is the case when an artificial absorber as well as a boundary integral is used for truncating the computational domain. A guideline is given as to how lowest and higher order TVFEs shall be combined for optimal performance of the proposed class of hierarchical TVFEs. 1

1 Introduction Node based expansions in finite element method (FEM) solutions are suitable for modeling scalar quantities but typically not so for simulating vector electromagnetic fields. When assigning vector field values to element nodes, values may need to be specified at locations where the true field is undefined (corners, edges), spurious modes can be generated and the enforcement of the boundary conditions occurring in electromagnetics can be a challenging task. Tangential vector finite elements (TVFEs) based on expanding a vector field in terms of values associated with element edges have been shown to be free of these shortcomings [1] and therefore TVFEs are of considerable practical interest. A TVFE is referred to as polynomial-complete to a given order, say n, if all possible polynomial variations up to and including order n are captured within the element and on the element boundary. If a TVFE captures polynomial variations of order n interior to the element and polynomial variations of order n - 1 on the element boundary, it is referred to as polynomial-complete to order n - 1. For a triangular element, polynomial-complete expansion to order I, 1 and 2 requires 3, 6 and 12 vector basis functions, respectively. Nedelec pointed out [2] [3] that it is not necessarily advantageous to employ polynomialcomplete TVFEs when applying the FEM. It was proven that a polynomial-complete expansion of a vector field A can be decomposed into a part representing the range space of the curl operator (V x A O0, A Z V+) and a part representing the null space of the curl operator (V x A = 0, A = V+). For representation of electromagnetic fields, it suffices to employ a TVFE that is complete in the range space of the curl operator. Such a TVFE is referred to as a mixed-order TVFE or (less commonly) a reduced-order TVFE. For rigorous criteria for spurious mode elimination and extensive discussions of mixed-order TVFEs versus polynomial-complete TVFEs, see [2] [3]. For discussions of element completeness and spurious modes for TVFEs, see also [4] [5] [6] [7]. A class of TVFEs is referred to as hierarchical if the vector basis functions forming the nth order TVFE are a subset of the vector basis functions forming the (n + l)th order TVFE. This desirable property allows for selective field expansion using different order TVFEs in different regions of the computational domain. Hence, lowest order TVFEs can be employed in regions where the field is expected to vary slowly whereas higher order TVFEs can be employed in regions where rapid field variation is anticipated. This selective choice of TVFEs 2

over the computational domain can lead to a memory and CPU time reduction as well as improved accuracy. Scalar FEM analysis using hierarchical finite elements is a well-known approach, see for example [8] [9]. For a triangular element, the lowest order TVFE was originally introduced by Whitney [10] and is referred to as the 'Whitney TVFE' or the 'Whitney element'. It provides a constant field value along element edges and a linear field variation inside the element. Hence, the Whitney TVFE is complete to order -. Several non-hierarchical TVFEs of higher order 2 than the Whitney TVFE have been introduced by Mur and de Hoop [11], Lee et al. [12] and Peterson [13]. Recently, Carrie and Webb [14] presented a hierarchical TVFE derived directly from a set of scalar finite element basis functions [9]. A hierarchical TVFE for a tetrahedral element (the special case of a triangular element follows trivially from the case of a tetrahedral element) up to and including second order was introduced by Webb and Forghani [15]. Recent reviews of TVFEs have been given by Peterson and Wilton [7] and by Graglia et al. [16]. The discussion above presents the concepts of polynomial-complete, mixed-order and hierarchical TVFEs and summarizes different TVFEs that have been proposed for a triangular element. However, hierarchality has never been exploited for TVFEs for a triangular element to numerically clarify whether selective field expansion using hierarchical TVFEs is an advantageous approach. The purpose of this paper is to introduce a class of hierarchical TVFEs for a triangular element and to demonstrate its performance when a selective field expansion is employed. The hierarchical TVFEs differ from those previously presented for a triangular element. The derivation is based on an attractive set of higher order vector basis functions recently presented by Popovic and Kolundzija [17] [18] for expanding surface currents in conjunction with method of moments (MoM) solutions. We demonstrate the performance of the proposed class of hierarchical TVFEs by comparing FEM results obtained using the Whitney TVFE and the proposed hierarchical TVFEs (applying a selective field expansion) to MoM solutions. The structure of this paper is as follows. Section 2 presents the derivation of the proposed class of TVFEs based on the expansion introduced by Popovic and Kolundzija. Section 3 discusses the merits of the proposed TVFEs and compares them to existing TVFEs. Expressions in terms of simplex coordinates are given and vector plots are added to provide a physical understanding of the TVFEs. Section 4 presents a set of numerical results that 3

demonstrate the performance of the proposed class of hierarchical TVFEs. Section 5 summarizes and concludes the paper and outlines future work to be carried out. Parts of this work were presented earlier [19]. 2 Derivation of a class of hierarchical TVFEs The proposed class of hierarchical TVFEs is based on an expansion introduced by Popovic and Kolundzija [17] [18] for the surface current on a perfectly electrically conducting (PEC) generalized quadrilateral. For this expansion, it is demonstrated in [17] [18] that the surface current can be expanded using approximately ten unknowns per square wavelength as opposed to approximately one hundred unknowns per square wavelength for traditional subdomain pulse basis functions. This suggests that the expansion introduced by Popovic and Kolundzija is very efficient and below we convert it to a class of TVFEs for a triangular element. As a degenerate case of the generalized quadrilateral considered in [17] [18], we consider a triangular element with nodes 1, 2 and 3 described by position vectors ri, r2 and r3 with respect to the origin 0 of a rectangular coordinate system, see Fig. 1. The edges from node 1 to node 2, node 2 to node 3 and node 3 to node 1 are referred to as edge #1, edge #2 and edge #3, respectively. The area of the triangle is denoted by A. Simplex (or area) coordinates (i, (2 and (3 at a point P described uniquely by a position vector r are defined in the usual manner, viz. n = An/A where An denotes the area of the triangle formed by P and the endpoints of the edge opposite to node n. We let n denote a unit normal vector to the surface of the triangle. Popovic and Kolundzija expands the surface current J, over the triangle as [17] 3 3 Js= J= Z Tnv (1) n=l n=l where r - rn Vn- 2A (2) is a vector whose direction is from node n to P and On is a polynomial function of position that provides the amplitude variation of the vector current component Jsn = TnVn. The polynomial II contains a number of unknown expansion coefficients. Its specific form is 4

irrelevant at this point and will be given later. As in the Rao-Wilton-Glisson expansion [20], Js, has no normal component along the two edges sharing node n and Jsn has both a normal and a tangential component along the edge opposite to node n. Thus, the quantity 2A Fn - li x sn = J-n A x Vn - O 2A = -nWn (3) with ft x (r —rn) n 2A (4) has no tangential component along the two edges sharing node n and has both a tangential and a normal component along the edge opposite to node n. This suggests that the vector basis functions multiplying the expansion coefficients in the expansion of Fn can be employed as vector basis functions for the edge opposite to node n when applying the FEM. Considering all three edges, the FEM expansion of an unknown vector quantity F becomes 3 3 3 F n x Js Jsn - Fn nWn (5) n=l n=l n=l where expressions for On (depending upon the order of the expansion) and Wn (independent of the order of the expansion) are to be presented. Introducing coordinates over the triangle and using relations from [17], it is shown in Appendix A that Wl = -2VC3 - C3VC2 (6) W2 -= 3VC1 - ClVC3 (7) W3 = C1V(2 - C2V(1 (8) The polynomial An is a function of position that in terms of a number of unknown expansion coefficients provides the amplitude variation of the vector component Wn. It can be defined using normalized coordinates Un and Vn over the triangle. Specifically, we choose Un = 0 at node n, Un = 1 along the edge opposite to node n and Vn = ~1 along the two edges sharing node n. From [17], we have n =2 bn n -2-) - (9) j=1 i=3 5

where rn and nu are integer constants determining the order of the approximation and b3n and an are the expansion coefficients. Also, ul = 1 2 + (3, UlV = (2 - 3 U2 -= 3 + -1l U2'V2 = - (1, U3 -= 1 + f2 and U3V3 = - - 2. The expansion (5) for F along with the expressions (6)-(8) for W1, W2 and W3 and the expression (9) for on describes the proposed class of TVFEs. However, a certain simplification provides a more familiar form. By regrouping the terms in (9) for A, the expansion (5) for F can be cast into 3 Nnvnu F = n c^Uw U (10) k=l m=l where NnTv = nv(nu - 1) denotes the total number of vector basis functions per edge for the given values of nu and nv. Also, cnukm are expansion coefficients corresponding to edge #k while Wnvnu are vector basis functions associated with edge #k. Wvu is given (in terms of the simplex coordinates (i, 2 and (3) as a function of (1, (2 and 43 times a direction vector (1V(2 - (2V(1 for k = 1, 2V(3 - (3V2 for k = 2 and 3V(1 - 1V(3 for k = 3). Except for normalization constants, the vector basis functions Wnvu form the proposed class of TVFEs. Each choice of nu and rn gives a different set of vector basis functions Wnvnu and thereby a different TVFE. 3 Discussion and comparison to traditional TVFEs In this section, we examine the properties of the TVFEs introduced in the previous section and contrast them to those of traditional TVFEs. From (10), it is seen that the total number of vector basis functions employed for expanding a vector field over a triangular element is 3Nnvnu = 3nv (nu - 1). For (nv, nu) = (1,2), we have 3Nnvnu = 3. This is the number of vector basis functions required -for the lowest order polynomial-complete expansion (complete to order -) over a triangular element. For other values of rn and n,, the expansion (10) can contain vector basis function of higher order than that to which the element is complete. Note, however, that this does not imply that the expansion (10) is not a polynomial-complete expansion (10) is not a polynomial variations to a given order are represented for properly chosen nv and nu. In addition, polynomial variations of higher orders may be represented and can be excluded if so desired. From (10), we recover for (n, nu) = (1, 2) the lowest order TVFE introduced by Whitney 6

and characterized by three vector basis functions. For larger values of rin and nu. the proposed TVFEs include additional vector basis functions all of which maintain the same fundamental direction vectors ((1V(2 - C2V(1 for edge #1, (2V(3 - 43VC(2 for edge #2 and (C3V( - (1V(3 for edge #3). Thus, the proposed higher order vector basis functions differ from the lowest order vector basis functions only in magnitude and hence in a given point of the triangle, the field is represented as a linear combination of vector basis functions having only three fundamental directions. This is one of the major differences between the proposed and traditional TVFEs. For the latter, the higher order vector basis functions differ from the lowest order vector basis functions in both magnitude and direction. The field in a given point of the triangle is again represented as a linear combination of vector basis functions but in this case the number of fundamental vector directions used for representing the field grows with the order of the TVFE. An important property of the proposed TVFEs is that the vector basis functions W4nvu, k=1,2,3, m=l,,Nnvnu, are a subset of the vector basis functions Wk +l), k=1,2,3, m=l,..,N(nv+l)(nu+l). This shows that the proposed TVFEs are hierarchical, a very desirable property. Hierarchical TVFEs are ideally suited for employing an efficient selective field expansion where different order TVFEs (in this case different values of n, and nu) are employed in different regions of the computational domain. Hence, for a uniform mesh, the lowest order TVFE can be employed in regions where the field is expected to experience smooth variation (regions where the relative material parameters are (nearly) unity, away from edges, etc.) whereas a higher order TVFE can be employed in regions where the field is expected to vary rapidly (near edges, close to material boundaries, in dense materials, etc.). Similarly, for a non-uniform mesh (for example where geometric complexity requires detailed meshing in a given region), the lowest order TVFE can be employed where the mesh is dense while a higher order TVFE can be employed where the mesh is coarse. Regions where higher order TVFEs are employed can be fixed a priori or an adaptive scheme can be developed where lowest order TVFEs are initially employed throughout the computational domain and higher order TVFEs are subsequently employed in regions where the error is estimated to be large. Below, we examine the TVFEs given by the expansion (10) for different values of nv, and nu and compare them to traditional TVFEs. Explicit expressions for some of the vector basis functions discussed here are given in Appendix B. 7

For the special case of (n, rnu) = (1,2), we obtain a set of three vector basis functions W1, k=1,2,3, forming a TVFE identical to the Whitney TVFE [10], see Appendix B. This result was expected since the lowest order expansion adopted by Popovic and Kolundzija is identical to the Rao-Wilton-Glisson expansion [20] whose vector basis functions reduce to the Whitney vector basis functions when converted using the procedure applied above. This lowest order TVFE provides a constant tangential field along element edges, a linear variation of the normal field component along element edges (referred to as CT/LN field variation along element edges) and a linear field variation interior to the element. A traditional second order polynomial-complete TVFE is described by twelve vector basis functions, see [5]. However, by excluding the four vector basis functions associated with the null space of the curl operator, this second order polynomial-complete TVFE can be reduced to a mixed-order TVFE [2] [3] formed by eight vector basis functions as presented by Peterson [13], see Appendix B. The latter provides a linear variation of the tangential field component along element edges, a quadratic variation of the normal field component along element edges (referred to as LT/QN field variation along element edges) and quadratic field variation interior to the element. To compare Peterson's mixed-order TVFE to one of the proposed TVFEs, the TVFE corresponding to (nv nu) = (2,3) is also reduced to eight vector basis functions providing LT/QN field variation along element edges and quadratic field variation interior to the element, see Appendix B. Peterson's mixed-order TVFE [13] has the desirable property of being complete to second order in the range space of the curl operator. This property ensures a complete second order expansion of a field with non-zero curl and guarantees eigenvalue solutions free of spurious non-zero eigenvalues. Due to the existence of a linear transformation (see Appendix B) from Peterson's eight vector basis functions to the proposed eight vector basis functions, the proposed eight vector basis functions form a mixed-order TVFE that has the same desirable property. However, the two TVFEs are not identical. Their vector basis functions span the same space but have different properties and may not be equally efficient numnerically. For both mixed-order TVFEs, six of the vector basis functions provide a linearly varying tangential component along element edges while the remaining two vector basis functions (identical for the two different TVFEs) provide a quadratic variation of the normal field component along element edges. However, the linear variation of the tangential component along element edges is obtained in two different ways. For Peterson's mixed-order TVFE, 8

the two unknowns per edge represent the magnitude of the field at edge endpoints. For the proposed mixed-order TVFE, the two unknowns per edge represent the average field value along the edge and the deviation from this average value at edge endpoints. To pictorically illustrate the behavior of the proposed vector basis functions, we consider a triangular element where the nodes 1, 2 and 3 have the coordinates (0, 0), (1, 0) and (0,1) in a rectangular (x, y) coordinate system. The lowest order vector basis function (2V(3 - 3v(2 associated with edge #2 is plotted in Fig. 2. For the proposed mixed-order TVFE, the linear variation of the tangential field along edge #2 is provided by C2V(3 - C3V(2 (due to the hierarchical nature of the proposed class of TVFEs) and by (2 - (3) (2'C3 - (3 2) plotted in Fig. 3. For Peterson's mixed-order TVFE, the linear variation of the tangential field along edge #2 is provided by (2V(3 and by (3V(2 individually, see Fig. 4 and Fig. 5, and this mixed-order TVFE therefore does not compare to the lowest order TVFE in a hierarchical fashion. The two vector basis functions providing quadratic normal field variation along element edges are the same for the two mixed-order TVFEs. The vector basis function 1(C2VC3 - 3V(2) associated with edge #2 (zero field along edge #2) is plotted in Fig 6. All vector basis functions are seen to provide the desired variation. For TVFE orders higher than (n, nu) = (2, 3), a similar procedure can be followed. A proposed TVFE can again be reduced to form a hierarchical mixed-order TVFE differing from (but spanning the same space as) a traditional mixed-order TVFE. The principle is similar but the procedure becomes more cumbersome and such TVFEs shall not be explicitly examined in this paper. 4 Numerical results In the previous section, properties of the proposed TVFEs were examined. Specifically, the lowest order TVFE (corresponding to (n,r nu) = (1,2), see also [10]), Peterson's mixedorder TVFE [13] and the proposed mixed-order TVFE (corresponding to (n11, nu) = (2,3) and reduced from twelve to eight vector basis functions) were compared. It is the aim of this section to numerically demonstrate the performance of the proposed class of hierarchical TVFEs when the field is selectively expanded using the lowest order TVFE in part of the computational domain and the proposed mixed-order TVFE in the remaining part of the computational domain. A guideline will be given as to how lowest and higher order TVFEs shall be combined for optimal (with respect to accuracy, memory and CPU-time) 9

performance. A FEM computer code was developed to evaluate the scattering of a transverse electric (TE) or transverse magnetic (TM) polarized plane wave by an arbitrary infinite cylinder composed of PECs and isotropic dielectric and/or magnetic materials. The code is based on a standard FEM formulation for two-dimensional problems where the transverse field component is expanded using a TVFE and the solution domain is truncated using a homogeneous isotropic artificial absorber (AA) (a fictitious material of relative permittivity and permeability 1 - j2.7 backed by a PEC structure) of thickness 0.25Ao (Ao denotes the free space wavelength) placed a distance 0.5Ao from the scatterer [21]. The resulting sparse FEM matrix equation system is solved using a quasi minimal residual solver [22]. For validation, MoM results were successfully compared to FEM results using each of the three TVFEs individually as well as the two different TVFEs selectively over the computational domain. Let us consider a square PEC cylinder of side length Ao situated in a free space region characterized by the permittivity Eo and the permeability luo. Centered on the upper side of the cylinder is a rectangular groove of length Ao/2 and of height Ao/4. The groove is filled with a material characterized by the relative permittivity r = 2 - j2 and the relative permeability r, = 2 - j2. The cylinder is illuminated by a TE (with respect to the cylinder axis) polarized homogeneous plane wave whose propagation vector forms a 45~ angle with all sides of the cylinder, as illustrated in Fig. 7. In the following, we compare the scattering by the cylinder using different TVFE options and different uniform discretizations to demonstrate the merits of the proposed TVFEs when the field is selectively expanded over the computational domain. In Fig. 8, we compare results for the two-dimensional radar cross section (RCS) (2-D normalized to Ao as a function of the observation angle q 1. The MoM result is denoted 'MoM'. For a mesh where the generic element edge size is 0.15Ao, the FEM result using the lowest order TVFE is denoted 'FEM - 1 TVFE - Coarse' and the FEM result using selective field expansion (with the groove and a layer surrounding the scatterer as the region in which the mixed-order TVFE is employed) is denoted 'FEM - 2 TVFEs - Coarse'. For a mesh where the generic element edge size is O.lo, the FEM result using the lowest order TVFE is denoted 'FEM - 1 TVFE - Dense'. The 'FEM - 1 TVFE - Coarse' result is seen to compare reasonably well with the exact 10 = 450 corresponds to backscatter and q = 2250 corresponds to forward scattering, see Fig. 7. 10

MoM result. However, discrepancies can be seen and this is not surprising since the mesh is relatively coarse. For the denser mesh, the 'FEM - 1 TVFE - Dense' result shows a slight improvement. However, by keeping the original mesh and employing the proposed mixed-order TVFE close to the scatterer where the field can be expected to vary rapidly and accurate modeling is therefore necessary, the 'FEM - 2 TVFEs - Coarse' result shows a significant improvement. It matches the MoM result exactly except in regions surrounding nulls and it was obtained using less computational resources (less unknowns, less non-zero matrix entries and less matrix solution time) than the 'FEM - 1 TVFE - Dense' result. In conclusion, we observe selective field expansion to be superior to the more traditional approach of using a denser mesh and the same TVFE throughout the computational domain. We now consider a slightly different cylinder geometry by introducing a slab of length Ao and height Ao/4 on top of the cylinder. As depicted in Fig. 9, the groove is filled with free space and the slab has the relative permittivity Er = 2 - j2 and the relative permeability Pr = 2 - j2. For the same illumination as before, results similar to those in Fig. 8 are given in Fig. 10 and they reinforce the conclusions from the previous case: The 'FEM - 1 TVFE - Coarse' result compares reasonably well with the exact MoM result and the 'FEM - 2 TVFEs - Coarse' result is, though found using less computational resources than the 'FEM - 1 TVFE- Dense' result, significantly more accurate than the 'FEM - 1 TVFE - Dense' result. Explicit parameter values quantifying the computational savings for the results in Fig. 8 and Fig. 10 are given in Tab. 1 and Tab. 2, respectively. In both cases, improved accuracy is obtained for less non-zero matrix entries (i.e., less memory) and less solution time. To test the validity of the reported observations for an alternative mesh truncation scheme, the FEM code was modified to use a boundary integral (BI) for truncating the finite element mesh. Where the AA mesh truncation scheme is approximate, the BI is (at least until discretized and coupled with a FEM system) exact and hence it is attractive for truncating finite element meshes. For our test, the integral contour is situated a slight distance away from the scatterer so that a piecewise constant (lowest order) expansion can be employed for discretizing the BI. As illustrated in Fig. 11, we consider a rectangular PEC cylinder of width 3A0 and height 0.25Ao covered by a dielectric material of width 3A0 and height Ao whose relative permittivity is Er = 2 - j'0.5. On top of the dielectric is a grating structure of height 0.25Ao consisting of three PEC strips of lengths 0.75A0, 0.5Ao and 0.75Ao, respectively, separated by dielectric inserts of length 0.5A0 having the relative permittivity 11

er = 10. A structure of this type (but of different size and different material composition) is of practical interest for guiding electromagnetic waves and below we demonstrate how a selective field expansion can lead to accurate modeling of the fields in and near the grating structure and hereby accurate prediction of the scattered field. The structure is situated in free space and illuminated as the previous two cylinders. Results similar to those in Fig. 8 and Tab. 1 are given in Fig. 12 and Tab. 3. The results again reinforce the conclusions reported above, except that the matrix solution time for the 'FEM - 2 TVFEs - Coarse' result is larger than that for the 'FEM - 1 TVFE - Dense' result. This fact is due to the condition number of the resulting FEM-BI equation system and might change if a different iterative solver had been applied. Further, we note that the pre-processing time is significantly larger for the 'FEM - 1 TVFE - Dense' result than the 'FEM - 2 TVFEs - Coarse' result due to the larger BI system. This must be kept in mind when interpreting Tab. 3. We note that at least two other approaches could be utilized for improving the accuracy of FEM results. For example, higher order TVFEs (either Peterson's mixed-order TVFE or the proposed mixed-order TVFE) could be applied throughout the computational domain. This approach was tested and the two mixed-order TVFEs gave similar and accurate results but could not measure up with the selective approach in terms of computational resources. Alternatively, non-uniform meshing could be utilized. However, this could be employed for all the TVFE options described in this paper and was therefore not tested. Moreover, mesh regeneration for improved solution accuracy is not an attractive option. Nevertheless, it; is reasonable to assume that this approach would lead to accurate results with a denser mesh close to the scatterer where the field is expected to vary rapidly. 5 Conclusions and future work We introduced a class of hierarchical TVFEs for FEM discretization. The properties of the proposed class of TVFEs were discussed and a comparison to those of traditional TVFEs was given. A set of numerical results were presented that demonstrate the effectiveness of the proposed class of hierarchical TVFEs when the computational domain is selectively discretized using the lowest order TVFE in part of the domain and a mixed-order TVFE in the remaining part of the domain. Hence, the computational domain can initially be discretized using lowest order TVFEs and the accuracy of the solution can then be improved by selectively superimposing more vector basis functions where rapid field variation is anticipated, i.e. in regions near edges, near material boundaries, in dense dielectrics etc. 12

Although the class of hierarchical TVFEs was presented for a triangular element, the approach has the potential to be more general. The derivation of a class of hierarchical TVFEs for a generalized quadrilateral and, as a special case, a curved triangle would again begin with a suitable polynomial expansion for a surface current as presented by Popovic and Kolundzija. Such elements conform well to almost all geometries and are thus attractive for FEM discretization. Further, corresponding three-dimensional elements could be developed. References [1] J.P. Webb, 'Edge elements and what they can do for you', IEEE Trans. Magn., vol. MAG-29, pp. 1460-1465, March 1993. [2] J.C. Nedelec, 'Mixed finite elements in I13', Num. Math., vol. 35, pp. 315-341, 1980. [3] J.C. Nedelec, 'A new family of mixed finite elements in 13', Num. Math., vol. 50, pp. 57-81, 1986. [4] Z.J. Cendes, 'Vector finite elements for electromagnetic field computation', IEEE Trans. Magn., vol. MAG-27, pp. 3958-3966, September 1991. [5] A.F. Peterson and D.R. Wilton, 'A rationale for the use of mixed-order basis functions within finite element solutions of the vector Helmholtz equation', Proc. of the 11th Annual Review of Progress in Applied Computational Electromagnetics, Monterey, CA, USA, pp. 1077-1084, March 1995. [6] D. Sun, J. Manges, X. Yuan and Z.J. Cendes, 'Spurious modes in finite-element methods', IEEE Antennas Propagat. Mag., vol 37, pp. 12-24, October 1995. [7] A.F. Peterson and D.R. Wilton, 'Curl-conforming mixed-order edge elements for discretizing the 2D and 3D vector Helmholtz equation', Chapter 5 of T. Itoh, G. Pelosi and P.P. Silvester (Editors), Finite element software for microwave engineering, John Wiley and Sons, Inc., 1996, pp. 101-124. [8] J.P. Webb and S. McFee, 'The use of hierarchal triangles in finite-element analysis of microwave and optical devices', IEEE Trans. Magn., vol. MAG-27, pp. 4040-4043, September 1991. [9] J.P. Webb and R. Abouchacra, 'Hierarchal triangular elements using orthogonal polynomials', Int. J. for Num. Meth. in Eng., vol 38, pp. 245-257, 1995. 13

[10] H. Whitney, Geometric integration theory, Princeton University Press, 1957. [11] G. Mur and A.T. de Hoop, 'A finite-element method for computing three-dimensional electromagnetic fields in inhomogeneous media', IEEE Trans. Magn., vol. MAG-21, pp. 2188-2191, November 1985. [12] J.-F. Lee, D.-K. Sun and Z.J. Cendes, 'Full-wave analysis of dielectric waveguides using tangential vector finite elements', IEEE Trans. Microwave Theory Tech., vol. MTT-39, pp. 1262-1271, August 1991. [13] A.F. Peterson, 'Vector finite element formulation for scattering from two-dimensional heterogeneous bodies', Trans. Antennas Propagat., vol. AP-42, pp. 357-365, March 1994. [14] C. Carrie and J.P. Webb, 'Hierarchal triangular edge elements using orthogonal polynomials', Proc. of the IEEE Antennas and Propagation Society International Symposium 1997, Montreal, Quebec, Canada, Volume 2, pp. 1310-1313, July 1997. [15] J.P. Webb and B. Forghani, 'Hierarchal scalar and vector tetrahedra', IEEE Trans. Magn., vol. MAG-29, pp. 1495-1498, March 1993. [16] R. Graglia, D.R. Wilton and A.F. Peterson, 'Higher order interpolary vector bases for computational electromagnetics', Trans. Antennas Propagat., vol. AP-45, pp. 329-342, March 1997. [17] B.D. Popovic and B.M. Kolund2ija, Analysis of metallic antennas and scatterers, IEE Electromagnetic Waves Series, vol. 38, 1994. [18] B.M. Kolund2ija and B.D. Popovic, 'Entire-domain Galerkin method for analysis of metallic antennas and scatterers', Proc. IEE H, vol. 140, pp. 1-10, 1993. [19] L.S. Andersen and J.L. Volakis, 'A novel class of hierarchical higher order tangential vector finite elements for electromagnetics', Proc. of the IEEE Antennas and Propagation Society International Symposium 1997, Montreal, Quebec, Canada, Volume 2, pp. 648-651, July 1997. [20] S.M. Rao, D.R. Wilton and A.W. Glisson, 'Electromagnetic scattering by surfaces of arbitrary shape', Trans. Antennas Propagat., vol. AP-30, pp. 409-418, May 1982. 14

[21] L.S. Andersen, 'Summary of finite element method formulations for solution of twodimensional bounded- and unbounded-domain electromagnetic problems', Internal report, Radiation Laboratory, Department of Electrical Engineering and Computer Science, University of Michigan, Fall 1996. [22] Y. Saad, Iterative methods for sparse linear systems, PWS Publishing Company, 1996. 15

Appendix A In this appendix, explicit expressions for W1, W2 and W3 are derived. To derive an expression for W1, we introduce two coordinates (ui1, v) over the triangle. These are degenerates of similar coordinates for a generalized quadrilateral [17]. Their ranges are Umin < u1 < Umax and Vmin < Cl < Vmax where ul = umin, u1 = Umax, Vl = Vmin and vl =,max are constant coordinate lines describing the boundary of the generalized quadrilateral. For a triangular element (which is a degenerate generalized quadrilateral), we let v, = -1 along edge #3, vi = 1 along edge #1, ul = 0 at node 1 and u = 1 along edge #2. Using these coordinates, the position vector r defining P can be expressed as [17] r = r1 + ul ru1 + ut1v ru1l1 (11) where 1 ru = [(r3- r) + (r2- r)] (12) 2 1 rul = -(r2- r3) (13) 2 Further, ui and vl can be shown to be related to the simplex coordinates (1, (2 and (3 via U1 = 2 + (3 (14) vl (15) 52 + 3(15) From (4) for n = 1, trivial algebra then leads to W1 = C2V(3 - C3V(2 (16) To derive expressions for W2 and W3, we can similarly introduce coordinates (u2, v2) and (U3, v3) where v2,3 = f1 along the two edges shared by node 2, 3, U2,3 = 0 at node 2, 3 and 2,3 = 1 at the edge opposite to node 2, 3. The algebra is similar and we arrive at u2= C3 + -1 (17) C3 - C1 V2 (18) 2 3 + (1 W2 = C3V(l - ClVC3 (19) 16

U3 = '1 + (2 (20) V3 = (21) (1 + (2 W3 = l1V(2 - C2V(1 (22) Appendix B In this appendix, explicit expressions for the vector basis functions discussed in section 3 are presented. The basis functions are not normalized. The Whitney TVFE The Whitney TVFE is characterized by the three vector basis functions W1 = 1V2 - (2V1 (23) W2 = C2VC3 - C3VC2 (24) W3- 3V(1 - '1V(3 (25) Peterson's mixed-order TVFE Peterson's mixed-order TVFE is characterized by the eight vector basis functions W1 = 1VC2 (26) W2 = -2V 1 (27) W3 = C2V(3 (28) W42 = C3VC2 (29) W52= V3Vl (30) W2 = 1iV(3 (31) W7= C3(C1VC2 - C2V(1) (32) W8 = 1(C2VC3- C3VC2) (33) 17

The proposed mixed-order TVFE The proposed mixed-order TVFE is characterized by the eight vector basis functions W' = ClV(2 - C2V(C (34) W2 = C2V(3 - C3VC2 (35) W33 = C3V1 - C1VC3 (36) = (1 - (2)(1V(2 - )2V(1) (37) W3 = (C2 - C3)(C2V(3 - C3V(2) (38) W3 = (C3 - C1)(C3V(C - C1VC3) (39) W3 = 3(C1VC2 - C2VC1) (40) W8- =Cl(C2V3 - C3V(2) (41) Transformation between mixed-order TVFEs The traditional mixed-order TVFE presented by Peterson [13] and the proposed mixed-order TVFE presented in this paper are related through a simple linear transformation. Let [W2] be a column vector containing Peterson's eight vector basis functions W2, j = 1,,8, and [W3] be a column vector containing the proposed eight vector basis functions W3, i = 1,, 8. The relationship between the two sets of vector basis functions is then given by the matrix equation [W3] = [Aij][Wj] (42) where [Aij] is the sparse 8 x 8 transformation matrix 1 -1 0 0 0 0 0 0 0 0 1 -1 0 0 0 0 0 0 0 0 1 -1 0 0 [. 1 1 0 0 0 0 1 2 [Aj] - (43) 0 0 1 1 0 0 -2 -1 0 0 0 0 1 1 1 -1 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 1 18

List of Figures 1 Geometry of a triangular element and illustration of the vectors n x (r - r), n=1,2,3, describing the directions of the vector basis functions at the point P. 20 2 Plot of the proposed vector basis function C2V(3- 3V'(2 for a triangular element. 20 3 Plot of the proposed vector basis function ((3 - -3)(C2V(3 - 3V(2) for a triangular element............................... 21 4 Plot of the traditional vector basis function C2VC3 for a triangular element.. 21 5 Plot of the traditional vector basis function 3VX(2 for a triangular element. 22 6 Plot of the vector basis function (1((2Vi3 - (3V(2) for a triangular element. 22 7 Coated square cylinder with crack and illuminated by TE polarized plane wave. 23 8 Bistatic RCS of the cylinder in Fig. 7......................... 23 9 Coated square cylinder with crack loaded by a dielectric slab and illuminated by TE polarized plane wave............................ 24 10 Bistatic RCS of the cylinder in Fig. 9.......................... 24 11 Grating structure on top of a grounded dielectric and illuminated by TE polarized plane wave................................ 25 12 Bistatic RCS of the cylinder in Fig. 11.......................... 25 19

#2 Figure 1: Geometry of a triangular element and illustration of the vectors n' x r - r,, n=1,2,3A describing the directions of the vector basis functions at the point P. _., i::: ~~~:: -i-iirii~iiiiiiii.-............ - ~ i iiiii~ii!ii~ iii~~~........ -............ -i:... i.i:...!!i......................... Plto 0 rp sdvco iiiii........... E 7.. E — -....... i.............t f...................... t;-......., —... Figue 2:Plotof he poposd vetorbasi funtion42V3-43S2 fr a.riagula.eleent 20

Figure 3: Plot of the proposed vector basis function ((3 - 3)(2V(3 -3V(2) for a triangular element. Figure 4: Plot of the traditional vector basis function (2V(3 for a triangular element. 21

0 CD+ cj~ 0. A* CD CD (ID

co - o 2o 2 4o 4 F H1 r =2-j2 gr =2-j2 I - koo Figure 7: Coated square cylinder with crack and illuminated by TE polarized plane wave. '.. I -V 10 0 IL....... I I I I I I I MoM FEM - 1 TVFE - Coarse mesh FEM - I TN FEM - 2 TN - 0 0 1: /FE - Denser mesh /FEs - Coarse mesh " f v'V -10 -20 - -30 -L I. — tv 45 65 85 105 125 145 165 Observation angle q [degrees] 185 205 225 Figure 8: Bistatic RCS of the cylinder in Fig. 7. 23

o0 to Hi 4 o;ko 2 PE-C --- PEC \ko p 4 r =2-j2 Lr =2-j2 Odl Figure 9: Coated square cylinder with crack loaded by a dielectric slab and illuminated by TE polarized plane wave. 20 MoM............ FEM - 1 TVFE - Coarse mesh 10 ------ FEM - 1 TVFE - Denser mesh - -- - - FEM - 2 TVFEs - Coarse mesh 0 -Q ~ ~~~ — 2,,'. 225 Observation angle b [degrees] Figure 10: Bistatic RCS of the cylinder in Fig. 9. 24

H1 S 3Xo 4 to xo 2 2. 2o 2 3,o 4. o 4 PEC PEl PEC l PEC M- -?.Jim: | | / PEC P _ - 4 I / Er=10 ~r=2-j0.5 Figure 11: Grating structure on top of a grounded dielectric and illuminated by TE polarized plane wave. 30 MoM FEM - 1 TVFE - Coarse mesh 20 - - FEM - 1 TVFE - Denser mesh - FEM - 2 TVFEs - Coarse mesh -.". -, I/nA\./ A 0 -~. Ia ~:1'I' 225 Observation angle p [degrees] Figure 12: Bistatic RCS of the cylinder in Fig. 11. 25

List of Tables 1 Comparison of relevant parameters for the three FEM results in Fig. 8.... 27 2 Comparison of relevant parameters for the three FEM results in Fig. 10... 27 3 Comparison of relevant parameters for the three FEM results in Fig. 12... 27 26

FEM Code Unknowns Non-zero mat. entries Mat. solution time RMS error 1 TVFE - Coarse mesh 811 3955 9 seconds 3.8657 dB 1 TVFE - Denser mesh 3977 19664 111 seconds T 2.1059 dB 2 TVFEs - Coarse mesh 1070 7292 21 seconds 1.5614 dB Table 1: Comparison of relevant parameters for the three FEM results in Fig. 8. FEM Code Unknowns Non-zero mat. entries Mat. solution time RMS error 1 TVFE - Coarse mesh 900 4398 12 seconds 2.7353 dB 1 TVFE - Denser mesh 4469 22118 153 seconds 1.2883 dB 2 TVFEs - Coarse mesh 1280 9466 35 seconds 0.6290 dB Table 2: Comparison of relevant parameters for the three FEM results in Fig. 10. FEM Code Unknowns Non-zero mat. entries Mat. solution time RMS error 1 TVFE - Coarse mesh 1230 19662 172 seconds 5.3980 dB 1 TVFE - Denser mesh 2421 43961 301 seconds 0.9807 dB 2 TVFEs - Coarse mesh 1716 25884 534 seconds 0.7705 dB Table 3: Comparison of relevant parameters for the three FEM results in Fig. 12. 27

Hierarchical Tangential Vector Finite Elements for Tetrahedra Lars S. Andersen and John L. Volakis Radiation Laboratory, Department of Electrical Engineering and Computer Science, University of Michigan, Ann Arbor, MI 48109-2122, USA Abstract Tangential vector finite elements (TVFEs) overcome most of the shortcomings of node based finite elements for electromagnetic simulations. Hierarchical TVFEs are of considerable practical interest since they allow use of effective selective field expansions where different order TVFEs are combined within a computational domain. For a tetrahedral element, this paper proposes a set of hierarchical mixed-order TVFEs up to and including order 2.5 that differ from previously presented TVFEs. The hierarchical mixed-order TVFEs are constructed as the three-dimensional equivalent of hierarchical mixed-order TVFEs for a triangular element. They can be formulated for higher orders than 2.5 and the generalization to curved tetrahedral elements is straightforward. 1

1 Introduction Tangential vector finite elements (TVFEs) based on expanding a vector field in terms of values associated with element edges have been shown to be free of the shortcomings of node based finite elements [1]. TVFEs are therefore of considerable practical interest. Nedelec pointed out [2] [3] that it may not necessarily be advantageous to employ polynomialcomplete TVFEs when applying the finite element method (FEM). This lead to the introduction of attractive mixed-order TVFEs. A set of TVFEs is referred to as hierarchical if the vector basis functions forming the nth order TVFE are a subset of the vector basis functions forming the (n + l)th order TVFE and this desirable property allows for effective selective field expansions combining different order TVFEs in different regions of the computational domain. For a large class of electromagnetic problems, hierarchical mixed-order TVFEs are therefore attractive for FEM discretization. For a tetrahedral element, the lowest order TVFE was originally introduced by Whitney [4]. It provides a constant tangential / linear normal (CT/LN) field along element edges and a linear field at element faces and inside the element (complete to order 0.5). Mixed-order TVFEs providing a linear tangential / quadratic normal (LT/QN) field along element edges and a quadratic field at element faces and inside the element (complete to order 1.5) were presented by Lee et al. [5], Webb and Forghani [6], Savage and Peterson [7] and Graglia et al. [8]. Only the TVFE presented by Webb and Forghani compares to the Whitney TVFE in a hierarchical fashion. Non-hierarchical mixed-order TVFEs providing a quadratic tangential / cubic normal (QT/CuN) field along element edges and a cubic field at element faces and inside the element (complete to order 2.5) were presented by Savage and Peterson [7] (a correction to this TVFE was recently given by Peterson [9]) and Graglia et al. [8]. Hierarchical mixed-order TVFEs for a tetrahedral element have only been proposed up to and including order 1.5 [6] and these were written up by inspection. The purpose of this paper is to propose a set of hierarchical mixed-order TVFEs for a tetrahedral element beyond order 1.5. Specifically, hierarchical mixed-order TVFEs are presented up to and including order 2.5 where the mixed-order TVFE of order 1.5 differs from the one presented by Webb and Forghani [6]. We derive the hierarchical mixed-order TVFEs from existing non-hierarchical mixed-order TVFEs for a tetrahedral element [7] [9] and existing hierarchical mixed-order TVFEs for a triangular element [10] [11] in a systematic fashion that makes the proposed set of hierarchical mixed-order TVFEs for a tetrahedral element the direct three-dimensional 2

equivalent of the set of hierarchical mixed-order TVFEs for a triangular element [10] [11]. Hierarchical mixed-order TVFEs for higher orders than 2.5 can be derived by modifying the TVFEs proposed by Graglia et al. [8] and their extension to curved tetrahedral elements is straightforward via a simple mapping, see for instance [8]. 2 Formulation We consider a tetrahedral element with nodes 1, 2, 3 and 4 as shown in Fig. 1. The volume of the tetrahedron is denoted by V. Simplex (or volume) coordinates (1, (2, (3 and (4 at a point P are defined in the usual manner, viz. n = 1/I/V where V, denotes the volume of the tetrahedron formed by P and the nodes of the triangular face opposite to node n. Below, vector basis functions will be formulated in terms of these coordinates. Vector basis functions associated with an edge or a face of the tetrahedron will be referred to as edgebased or face-based vector basis functions, respectively. All other vector basis functions will be referred to as cell-based vector basis functions. A mixed-order TVFE of order 0.5 providing CT/LN variation along element edges and linear variation at element faces and inside the element is characterized by 6 linearly independent vector basis functions. Whitney initially presented 6 such vector basis functions [4]. The three-dimensional equivalent of the two-dimensional CT/LN vector basis functions presented in [10] [11] is identical to the vector basis functions presented by Whitney [4]. The 6 edge-based vector basis functions are 1 CV(j - Cjv, i < j. (1) A mixed-order TVFE of order 1.5 providing LT/QN variation along element edges and quadratic variation at element faces and inside the element is characterized by 20 linearly independent vector basis functions. Savage and Peterson [7] proposed the 12 edge-based vector basis functions iVCj, i j, (2) 1The vector basis functions presented in this paper are not normalized. Furthermore, the indices i, j and k in (1)-(12) are implicitly assumed to belong to the set {1,2,3,4}. 3

and the 8 face-based vector basis functions Ck(Cv-CV() }V,) (k - /( i) i < j < k. (3) (j ((kV(i -(i V(k) ' The 20 linearly independent vector basis functions (2)-(3) do not compare to the Whitney vector basis functions (1) in a hierarchical fashion. We propose to replace the 12 edge-based basis functions (2) by (CI- ( - )(CV( - CV) J -I - < (4) The 20 linearly independent vector basis functions (3)-(4) form a mixed-order TVFE of order 1.5 that compares hierarchically to the proposed mixed-order TVFE of order 0.5. A mixed-order TVFE of order 2.5 providing QT/CuN variation along element edges and cubic variation at element faces and inside the element is characterized by 45 linearly independent vector basis functions. Savage and Peterson 2 [7] [9] proposed the 18 edge-based vector basis functions C(2Ci - 1)VCj, ij, (5) Cj(Vc- j), i <j, (6) the 24 face-based vector basis functions Ck(2Ck - 1)(C, VO3 - CjVC) '1 k(2(k - 1)(ij - (V(i), i < j < k (7) (j(2(j - 1)((kv< - (iVk) J v(cck), i < j <, (8) (C2(CVC3 - (CVh), i k j ki, (9) and the 3 cell-based vector basis functions C2C3(ClVC4 - 4V(1C) C2C4(C1V(3 - C3V1) (10) C3C4(ClVC2 - C2VC1). The 45 linearly independent vector basis functions (5)-(10) do not compare to the Whitney vector basis functions (1) in a hierarchical fashion. We propose to replace the 18 edge-based 2A correction of the QT/CuN vector basis functions initially proposed by Savage and Peterson [7] was given by Peterson [9]. This corrected set is the one presented here. 4

basis functions (5)-(6) by (CV(C - (jvc( (- Cj)(CvCj - CjvC), i<j. (11) ((- (,)2(Co,- (CjV) <o) Further, we propose to replace the 8 face-based vector basis functions (7) by (kVCj- jV), i < j <. (12) (j('kV(i ( k) (12) The 45 linearly independent vector basis functions (8)-(12) form a mixed-order TVFE of order 2.5 that compares hierarchically to the proposed mixed-order TVFEs of order 0.5 and 1.5. The vector basis functions (1), (3)-(4) and (8)-(12) form a set of hierarchical mixed-order TVFEs of orders 0.5, 1.5 and 2.5, respectively. Such a set offers several advantages over nonhierarchical mixed-order TVFEs, especially for FEM solution of electromagnetic problems where the field varies non-uniformly over the computational domain. In such cases, a lower order TVFE can be employed in regions where the field varies smoothly whereas a higher order TVFE can be employed in regions where the field varies rapidly thus leading to an effective discretization of the unknown electromagnetic field. 3 Conclusions For a tetrahedral element, we proposed a set of hierarchical mixed-order TVFEs up to and including order 2.5. These differ from previously presented TVFEs and were constructed as the three-dimensional equivalent of hierarchical mixed-order TVFEs for a triangular element. TVFEs for higher orders than 2.5 can be formulated in a similar manner and the generalization to curved tetrahedral elements is straightforward. 5

References [1] J.P. Webb, 'Edge elements and what they can do for you', IEEE Trans. Magn., vol. MAG-29, pp. 1460-1465, March 1993. [2] J.C. Nedelec, 'Mixed finite elements in 313', Num. Math., vol. 35, pp. 315-341, 1980. [3] J.C. Nedelec, 'A new family of mixed finite elements in Rt3', Num. Math., vol. 50, pp. 57-81, 1986. [4] H. Whitney, Geometric integration theory, Princeton University Press, 1957. [5] J.F. Lee, D.K. Sun and Z.J. Cendes, 'Tangential vector finite elements for electromagnetic field computation', IEEE Trans. Magn., vol. MAG-27, pp. 4032-4035, September 1991. [6] J.P. Webb and B. Forghani, 'Hierarchal scalar and vector tetrahedra', IEEE Trans. Magn., vol. MAG-29, pp. 1495-1498, March 1993. [7] J.S. Savage and A.F. Peterson, 'Higher-order vector finite elements for tetrahedral cells', IEEE Trans. Microwave Theory Tech., vol. MTT-44, pp. 874-879, June 1996. [8] R. Graglia, D.R. Wilton and A.F. Peterson, 'Higher order interpolary vector bases for computational electromagnetics', IEEE Trans. Antennas Propagat., vol. AP-45, pp. 329-342, March 1997. [9] A.F. Peterson, Private communication, e-mail, August 1997. [10] L.S. Andersen and J.L. Volakis, 'A novel class of hierarchical higher order tangential vector finite elements for electromagnetics', Proc. of the IEEE Antennas and Propagation Society International Symposium 1997, Montreal, Quebec, Canada, vol. 2, pp. 648-651, July 1997. [11] L.S. Andersen and J.L. Volakis, 'Development and application of a novel class of hierarchical tangential vector finite elements for electromagnetics', Submitted to IEEE Trans. Antennas Propagat., September 1997. 6

List of Figures 1 Illustration of tetrahedral element and the numbering of nodes and edges... 8 7

4 Figure 1: 8

Accurate and efficient simulation of antennas using hierarchical mixed-order tangential vector finite elements for tetrahedra Lars S. Andersen and John L. Volakis Radiation Laboratory Department of Electrical Engineering and Computer Science University of Michigan Ann Arbor, MI 48109-2122, USA 1

Abstract Hierarchical mixed-order tangential vector finite elements (TVFEs) for tetrahedral elements are attractive for accurate and efficient finite element method simulation of complicated electromagnetic problems. They provide versatility in the geometric modeling of physical structures, guarantee solutions free of spurious modes and allow local increase of resolution by selective expansion of the unknown electromagnetic field, i.e. by combination of mixed-order TVFEs of different orders within a computational domain. For a realistic antenna radiation problem, this paper demonstrates that field expansion using lowest and higher order hierarchical mixed-order TVFEs selectively is vastly superior (in terms of accuracy, memory as well as CPU-time) to field expansion using a lowest order mixed-order TVFE only. 2

1 Introduction The finite element method (FEM) has proven attractive for simulation of complicated electromagnetic problems. Use of tetrahedral meshes provides versatility in the geometric modeling of physical structures and field expansions based on mixed-order tangential vector finite elements (TVFEs) guarantee solutions free of spurious modes and provide for easy enforcement of boundary conditions. Hierarchical mixed-order TVFEs (where the vector basis functions forming a mixed-order TVFE of a given order are a subset of the vector basis functions forming mixed-order TVFEs of higher orders) have the additional advantage of permitting selective field expansion. That is, they allow for combination of mixed-order TVFEs of different orders within a computational domain for efficient expansion of the unknown electromagnetic field. This is relevant for simulating a large class of realistic electromagnetic problems characterized by disjoint regions with high and low field variation. Thus, FEM analysis with fields expanded using hierarchical mixed-order TVFEs for tetrahedral elements is attractive for accurate and efficient solution of certain classes of electromagnetic problems. Hierarchical mixed-order TVFEs for tetrahedral elements have been proposed up to and including order 1.5 by Webb and Forghani [1] and up to and including order 2.5 by Andersen and Volakis [2, 3]. However, no test has yet been carried out that numerically demonstrates the potential of selective field expansion for realistic electromagnetic problems. The purpose of this paper is to do so using the hierarchical mixed-order TVFEs of order 0.5 and 1.5 for tetrahedral elements proposed in [2, 3]. For eigenvalue computation, the convergence rate of the hierarchical mixed-order TVFE of order 1.5 is shown to be comparable to that of a nonhierarchical mixed-order TVFE of order 1.5 [4]. For a realistic antenna radiation problem, 3

field expansion using the hierarchical mixed-order TVFEs of order 0.5 and 1.5 selectively is shown to be vastly superior (in terms of accuracy, memory as well as CPU-time) to field expansion using the mixed-order TVFE of order 0.5 only. This paper is organized as follows. Section 2 presents the hierarchical mixed-order TVFEs of order 0.5 and 1.5 for tetrahedral elements that will be used for field expansion. Section 3 presents the numerical results. Section 4 concludes the paper. 2 Presentation of TVFEs We consider a tetrahedral element with nodes 1, 2, 3 and 4. The volume of the tetrahedron is denoted by V. Simplex (or volume) coordinates 41, (2, C3 and (4 at a point P are defined in the usual manner, viz. n = Vn/V where V, denotes the volume of the tetrahedron formed by P and the nodes of the triangular face opposite to node n. A mixed-order TVFE of order 0.5 providing constant tangential / linear normal variation along element edges and linear variation at element faces and inside the element is characterized by 6 linearly independent vector basis functions. Whitney [5] initially presented such a TVFE. It is characterized by the 6 edge-based vector basis functions 1 VCj - jV, i < j. (1) A mixed-order TVFE of order 1.5 providing linear tangential / quadratic normal variation along element edges and quadratic variation at element faces and inside the element is characterized by 20 linearly independent vector basis functions. Andersen and Volakis [2, 3] 1The vector basis functions presented in this paper are not normalized. Furthermore, the indices i, j and k in (1)-(3) are implicitly assumed to belong to the set {1, 2, 3, 4}. 4

presented a mixed-order TVFE of order 1.5 that compares hierarchically to the mixed-order TVFE of order 0.5 presented by Whitney [5]. In addition to the 6 edge-based vector basis functions (1), it is characterized by the 6 edge-based vector basis functions (C< - )((Vj - jV(), <j (2) and the 8 face-based vector basis functions C(k((iV(- (jV(c), i < j < k. (3) Cj((kvC - CVk) 3 Numerical results The objective of this section is to numerically demonstrate the potential of the hierarchical mixed-order TVFEs of order 0.5 and 1.5 proposed in [2, 3] and summarized above. The eigenvalues of a homogeneous and isotropic rectangular cavity are determined numerically for different uniform tetrahedral meshes to show that the convergence rate of the hierarchical mixed-order TVFE of order 1.5 is comparable to that of a non-hierarchical mixed-order TVFE of order 1.5 [4]. The input impedance (and hereby the resonant frequency) of a probe-fed square metallic patch antenna backed by a dielectric-filled cavity recessed in an infinite metallic ground plane is determined numerically for different uniform tetrahedral meshes to show that field expansion using the hierarchical mixed-order TVFEs of order 0.5 and 1.5 selectively is vastly superior (in terms of accuracy, memory as well as CPU-time) to field expansion using the mixed-order TVFE of order 0.5 only. Consider a homogeneous, isotropic rectangular cavity of normalized dimensions 1 x 0.75 x 0.5. The exact eigenvalues for this geometry are well-known [6]. A FEM solution for the eigenvalues of the cavity is carried out for various uniform tetrahedral meshes of different 5

average edge length with the hierarchical mixed-order TVFEs of order 0.5 and 1.5 used for field expansion (for the formulation, see for instance [7]). The convergence rate for the two cases is illustrated in Fig. 1 where the average error of the first eight eigenvalues is plotted in percent as a function of the average edge length in the mesh (log-log plot). The approximate distribution around a straight line suggests that the average error decreases as xP for a decreasing average edge length. For the mixed-order TVFE of order 0.5, the exponent is p = 2.37 which is slightly larger than the expected value of 2 [4]. This is due to the very low average error 0.56% for the average edge length 0.175. Similarly, for the hierarchical mixed-order TVFE of order 1.5, the exponent is p = 4.66 which is again larger than the expected value of 4 [4] and the exponent p= 3.86 found in [4] for a, different and non-hierarchical mixed-order TVFE of order 1.5. This demonstrates that the hierarchical mixed-order TVFE of order 1.5 in [2, 3] has slightly better convergence properties than the non-hierarchical one in [4] for this particular geometry and for the employed meshes. However, a relatively large uncertainty range can be expected for such numerically obtained exponents and thus no general statement can be made regarding the rigor of results based on the two different mixed-order TVFEs of order 1.5. Consider a square metallic patch antenna backed by a rectangular cavity recessed in an infinite metallic ground plane, as illustrated in Fig. 2 (side view) and Fig. 3 (top view). The cavity-backed patch antenna is situated in free space characterized by the permittivity s0 and the permeability [o. The cavity is of dimensions 1.85 cm x 1.85 cm x 0.15 cm and filled with a dielectric material of permittivity 10 eo and conductivity 0.0003 S/cm. The patch is of side length 0.925 cm and centered in the cavity aperture. It is fed by a vertical coaxial line whose outer conductor is attached to the ground plane and whose inner conductor is 6

attached to the patch at the mid point of an edge, as illustrated in Fig. 2 and Fig. 3. The coaxial feed will be modeled as a vertical probe of constant current. An almost identical antenna was considered by Schuster and Luebbers [8]. In [8], the cavity walls and the ground plane was removed and a similar patch on a similar but finite grounded dielectric substrate was analyzed using the finite difference time domain method. In spite of these geometrical differences, the two antennas are expected to have nearly the same input impedance and, consequently, nearly the same resonant frequency since the dominant fields are confined to a volume under and in the near vicinity of the patch. The resonant frequency was estimated in [8] to be 4.43 GHz. The resistance at resonance was found to be 400 Q while the reactance was in the range of 230 Q to -170 Q close to resonance. The patch antenna is analyzed using the finite element - boundary integral (FE-BI) method (for the formulation, see for instance [7]). We discretize the cavity into tetrahedral elements and consequently discretize the surface forming the boundary between the cavity and free space into triangular faces. Two different TVFE options are applied. The first TVFE option is to use the mixed-order TVFE of order 0.5 throughout the mesh. For a mesh of average edge length 0.260 cm (Case 1), the input impedance is determined as a function of frequency and the resonant frequency of the patch is predicted. The coarse discretization of Case 1 means that this resonant frequency is most likely not accurate. For meshes of average edge lengths of 0.188 cm (Case 2), 0.153 cm (Case 3) and 0.133 cm (Case 4), more accurate resonant frequencies but also higher computational costs can be expected. The second TVFE option is to use the mixed-order TVFE of order 1.5 close to the radiating edges (where we expect high field variation) and the mixed-order TVFE of order 0.5 elsewhere (where we expect little field variation). For the meshes of average edge length 0.260 cm (Case 5) and 7

0.188 cm (Case 6), the input impedance is again determined and the resonant frequency is again predicted. The effectiveness of this approach (Case 5-6) in terms of accuracy / CPUtime / memory requirements is compared to the previous one (Case 1-4). The six cases are summarized in Tab. 1. Real and imaginary parts of the input impedance as a function frequency are given in Fig. 4 - 5 for Case 1-6 and corresponding resonant frequencies are provided in Tab. 1. For Case 1-4, a larger and larger resonant frequency is observed as the mesh becomes denser and denser. However, even for Case 4, the error as compared to the result obtained by Schuster and Luebbers is quite large (2.98 %) for resonant frequency computation. Use of selective field expansion (Case 5-6) leads to a significant accuracy improvement. Case 5 (error 2.42 %) gives a more accurate result than Case 1-4 and Case 6 (error 0.16 %) matches the result by Schuster and Luebbers almost exactly. The computational cost (number of unknowns, number of BI unknowns, number of non-zero matrix entries (memory usage) and CPU-time per frequency point) to obtain these results are also given in Tab. 1. It is evident that the second TVFE option corresponding to Case 5-6 is significantly more attractive than the first TVFE option corresponding to Case 1-4. Case 5 gives a more accurate result than Case 4 but uses only 4.22 % of the memory and 2.15 % of the CPU-time that Case 4 does. The accuracy of Case 6 is vastly superior to that of Case 4 and yet Case 6 uses only 14.88 % of the memory and 10.02 % of the CPU-time that Case 4 does. 8

4 Conclusions The potential of the hierarchical mixed-order TVFEs of order 0.5 and 1.5 for tetrahedral elements proposed by Andersen and Volakis [2, 3] was demonstrated. For eigenvalue computation, the convergence rate of the hierarchical mixed-order TVFE of order 1.5 was shown to be comparable to that of a non-hierarchical mixed-order TVFE of order 1.5 [4]. For a realistic antenna radiation problem, field expansion using the hierarchical mixed-order TVFEs of order 0.5 and 1.5 selectively was shown to be vastly superior (in terms of accuracy, memory as well as CPU-time) to field expansion using the mixed-order TVFE of order 0.5 only. 9

References [1] J.P. Webb and B. Forghani. 'Hierarchal scalar and vector tetrahedra'. IEEE Trans. Magn., vol. MAG-29, pp. 1495-1498, March 1993. [2] L.S. Andersen and J.L. Volakis. 'Hierarchical tangential vector finite elements for tetrahedra'. IEEE Microwave and Guided Wave Letters, vol. 8, pp. 127-129, March 1998. [3] L.S. Andersen and J.L. Volakis. 'Hierarchical tangential vector finite elements for tetrahedra'. In Proc. of the IEEE Antennas and Propagation Society International Symposium 1998, Atlanta, Georgia, USA, vol. 1, pp. 240-243, June 1998. [4] J.S. Savage and A.F. Peterson. 'Higher-order vector finite elements for tetrahedral cells'. IEEE Trans. Microwave Theory Tech., vol. MTT-44, pp. 874-879, June 1996. [5] H. Whitney. Geometric integration theory. Princeton University Press, USA, 1957. [6] R.F. Harrington. Time harmonic electromagnetic fields. McGraw-Hill Book Company, USA, 1961. [7] J.L. Volakis, A. Chatterjee, and L.C. Kempel. Finite element method for electromagnetics. IEEE Press, USA, 1998. [8] J.W. Schuster and R. Luebbers. Private communication. 10

List of Figures 1 Convergence rate for expansion of field within homogeneous and isotropic rectangular cavity using hierarchical mixed-order TVFEs of order 0.5 and 1.5. 12 2 Side view of square metallic patch antenna backed by a dielectric-filled rectangular cavity recessed in an infinite metallic ground plane............ 13 3 Top view of square metallic patch antenna backed by a dielectric-filled rectangular cavity recessed in an infinite metallic ground plane........... 14 4 Real part of the input impedance of the antenna in Fig. 2 and Fig. 3 for Case 1-6....................................... 15 5 Imaginary part of the input impedance of the antenna in Fig. 2 and Fig. 3 for Case 1-6....................................... 16 11

At) I l. 10l L0 a) 0) a) 3 1 v v TVFE of order 0.5 TVFE of order 0.5 - curve fit o o TVFE of order 1.5 TVFE of order 1.5 - curve fit v v v V v ^6 ~,ao / / 0 / / 0.3 0.1 O /,6 / 0 0 ( 11:( I 0.2 0.3 Average edge length 0.4 0.5 0.6 Figure 1: Convergence rate for expansion of field within homogeneous and isotropic rectangular cavity using hierarchical mixed-order TVFEs of order 0.5 and 1.5. 12

co go 0.15 cm 10 Co go 0.0003 S/cm Figure 2: Side view of square metallic patch antenna backed by a dielectric-filled rectangular cavity recessed in an infinite metallic ground plane. 13

Figure 3: Top view of square metallic patch antenna backed by a dielectric-filled rectangular cavity recessed in an infinite metallic ground plane. 14

CaseC 2 ase 2 2300 i 200 100.7 3.8 3.9 4 4.1 4.2 4.3 4.4 4.5 4.6 4.7 Frequency [GHz] Figure 4: Real part of the input impedance of the antenna in Fig. 2 and Fig. 3 for Case 1-6. 15

ro" N E * 4.1 4.2 4.3 4.4 4.5 4.6 4.7 Frequency [GHz] Figure 5: Imaginary part of the input impedance of the antenna in Fig. 2 and Fig. 3 for Case 1-6. 16

List of Tables 1 Computational effort for Case 1-6 for antenna in Fig. 2 and Fig. 3......... 18 17

Average Time per edge Resonant # of frequency TVFE length frequency # of # of BI matrix point Case order(s) [cm] [GHz] unknowns unknowns entries [sec] 1 0.5 0.260 3.974 345 120 17119 7.52 2 0.5 0.188 4.147 817 288 89695 44.78 3 0.5 0.153 4.258 1489 528 291359 222.92 4 0.5 0.133 4.302 2361 840 725791 771.59 5 0.5/1.5 0.260 4.323 827 120 30675 17.33 6 0.5/1.5 0.188 4.437 1467 288 107963 77.28 Table 1: Computational effort for Case 1-6 for antenna in Fig. 2 and Fig. 3. 18

Condition numbers for various FEM matrices Lars S. Andersen* and John L. Volakis Radiation Laboratory Department of Electrical Engineering and Computer Science University of Michigan Ann Arbor, MI 48109-2122, USA E-mail: andersen@umich.edu, volakis@umich.edu Abstract A study is presented that examines thoroughly the inter-relationships between condition numbers of finite element method matrices based on various interpolatory and hierarchical mixed-order tangential vector finite elements (TVFEs). The validity of the generally accepted premise that interpolatory higher order TVFEs lead to better conditioned matrices than hierarchical higher TVFEs is found to be questionable. Introduction Tangential vector finite elements (TVFEs) whose degrees of freedom are associated with element edges, faces and cells have been shown to be free of the shortcomings of node based expansions [1]. A TVFE is referred to as polynomial-complete to a given order, say n, if all possible polynomial variations up to and including order n are captured within the element and on the element boundary. Nedelec pointed out [2] that for representation of electromagnetic fields it suffices to employ TVFEs that are complete only in the range space of the curl operator, i.e. the space of fields having nonzero curl. Since such a TVFE captures polynomial variations of order n interior to the element and polynomial variations of order n - 1 for the tangential component along element edges, we call it complete to order n - 0.5 and refer to it as a mixed-order TVFE. A TVFE is referred to as interpolatory if the vector basis functions forming the TVFE interpolate to (tangential) field values at discrete points within or on the boundary of the element. A class of TVFEs is referred to as hierarchical if the vector basis functions forming the TVFE of a given order are a subset of the vector basis functions forming the TVFEs of higher orders. This desirable property allows for selective field expansion, i.e. combination of different order TVFEs within a computational domain for effective discretization of the unknown field. Linear equation systems resulting from application of the finite element method (FEM) are often solved using one of several available iterative solvers which all perform relatively poorly (require many iterations) when the FEM matrices are badly conditioned. This drawback makes it desirable to construct FEM matrices having small condition numbers. It is generally accepted that higher order TVFEs lead to larger condition numbers than the lowest order TVFE and that hierarchical highler order TVFEs lead to larger condition numbers than interpolatory higher order TVFEs [3]. These are two premises that make selective field expansion using hierarchical higher order rVFEs significantly less attractive.

TE formulation TM formulation TVFE Unnormalized Normalized Unnormalized Normalized Wh.3 3 5 6 Pe 428 90 479 149 Gr 144 37 144 61 An 402 43 451 72 We 827 71 926 101 Table 2: Condition numbers for global FEM matrices for rectangular waveguide. E formulation H formulation TVFE Unnormalized Normalized Unnormalized Normalized Wh 8 12 23 31 Pe 2173 287 5939 958 Gr 684 99 1387 239 An 1834 195 5341 656 We 4238 472 11564 1414 Table 3: Condition numbers for global FEM matrices for rectangular cavity. References [1] J.P. Webb. 'Edge elements and what they can do for you'. IEEE Trans. Magn., vol. MAG-29, pp. 1460-1465, March 1993. [2] J.C. Nedelec. 'Mixed finite elements in R3'. Num. Math., vol. 35, pp. 315-341, 1980. [3] J.S. Savage. 'Comparing high order vector basis functions'. In Proc. of the 14th Annual Review of Progress in Applied Computational Electromagnetics, Monterey, CA, USA, vol. 2, pp. 742-749, March 1998. [4] R. Graglia, D.R. Wilton, and A.F. Peterson. 'Higher order interpolary vector bases for computational electromagnetics'. IEEE Trans. Antennas Propagat., vol. AP-45, pp. 329-342, March 1997. [5] J.P. Webb and B. Forghani. 'Hierarchal scalar and vector tetrahedra'. IEEE Trans. Magn., vol. MAG-29, pp. 1495-1498, March 1993. [6] A.F. Peterson. 'Vector finite element formulation for scattering from two-dimensional heterogeneous bodies'. IEEE Trans. Antennas Propagat., vol. AP-42, pp. 357-365, March 1994. [7] J.S. Savage and A.F. Peterson. 'Higher-order vector finite elements for tetrahedral cells'. IEEE Trans. Microwave Theory Tech., vol. MTT-44, pp. 874-879, June 1996. [8] L.S. Andersen and J.L. Volakis. 'Development and application of a novel class of hierarchical tangential vector finite elements for electromagnetics'. Accepted for publication in IEEE Trans. Antennas Propagat. [9] L.S. Andersen and J.L. Volakis. 'Hierarchical tangential vector finite elements for tetrahedra'. IEEE MAicrowave and Guided Wave Letters, vol. 8, pp. 127-129, March 1998. [10] II. Whitney. Geometric integration theory. Princeton University Press, USA, 1957.

Tlhe condition numbers of element matrices based on the mixed-order TVFEs of order 1.5 developed by Graglia et al. [4] (interpolatory) and Webb and Forghani [5] (hierarchical) have been stlludied before [3]. However, no study has beeni presented that examines more thoroughly the interrelationships between condition numbers of global FEM matrices based on various interpolatory and hierarchical mixed-order TVFEs. It is the aim of this paper to do so for the interpolatory mixed-order TVFEs of order 1.5 developed by Peterson [6] (two-dimensional problems), Savage and Peterson [7] (three-dimensional problems) and Graglia et al. [4] and the hierarchical mixedorder TVFEs of order 1.5 developed by Andersen and Volakis [8, 9] and Webb and Forghani [5]. The study includes results for two- as well as three-dimensional problems, transverse electric (TE) as well as transverse magnetic (TM) field formulation (two-dimensional problems), electric (E) as well as magnetic (H) field formulation (three-dimensional problems) and unnormalized as well as normalized vector basis functions. This paper is organized as follows. First, vector basis functions for mixed-order TVFEs of order 0.5 and 1.5 for triangular and tetrahedral elements are reviewed. Subsequently, condition numbers for FEM matrices based on various TVFE options are presented. Finally, conclusions are given. Presentation of TVFEs We consider a triangular element of area A with nodes 1, 2 and 3. Simplex coordinates (1, (2 and (3 at a point P are defined in the usual manner, viz. (n = An/A where An denotes the area of the triangle formed by P and the nodes of the edge opposite to node n. Similarly, we consider a tetrahedral element of volume V with nodes 1, 2, 3 and 4. Simplex coordinates (1, (2, (3 and (4 at a point P are again defined in the usual manner, viz. (n = Vn/V where Vn denotes the volume of the tetrahedron formed by P and the nodes of the triangular face opposite to node n. Table 1 presents the vector basis functions characterizing the mixed-order TVFE of order 0.5 developed by Whitney [10], the interpolatory mixed-order TVFEs of order 1.5 developed by Peterson [6] (two-dimensional problems), Savage and Peterson [7] (three-dimensional problems) and Graglia et al. [4] and the hierarchical mixed-order TVFEs of order 1.5 developed by Andersen and Volakis [8, 9] and Webb and Forghani [5]. The four mixed-order TVFEs of order 1.5 span the same space, their difference is the form of the vector basis functions. The indices i, j and k in Table 1 are implicitly assumed to fulfil i < j < k and belong to the set {1,2,3} for triangular elements and the set {1,2,3,4} for tetrahedral elements. We note that the vector basis functions are not normalized. To normalize them, each edge-based vector basis function must be multiplied by the length of the edge it is associated with. TVFE Vector basis functions Order Label Edge-based Face-based 0.5 Wh (iV(j - iV(1.5 Pe (iV(j (k((iV( - (jV(i) jV( f. j((kV(, - (2v(k) 1.5 Gr (3(i - 1)((iV(o - (jV(i) (k((Vo - (jV(,) (3 - 1)((V( - (kV - )(vi - V) (v -V(k) 1.5 Anl ( - (v(, k((, vj - jV(f) -(C(i% ( ( - j)(i) - V() 1.5 \ev - (,jV(i (k(( --, jV~) __ __ V +, V(( G ((kvk( - (,v-k) Table l: Definlition of vector basis functions for five different TVFEs.

Condition numbers for various FEM matrices Linear equation systems resulting from application of the FEM are often solved using one of several available iterative solvers which all perform relatively poorly (require many iterations) when the FEM matrices are badly conditioned. This is the case for (generalized) eigenvalue problems as well as excitation problems. The drawback makes it desirable to use TVFEs that lead to global FEM matrices having small condition numbers. However, given a TVFE, different applications lead to structurally different global matrix systems and hence it is impossible to uniquely define a global matrix whose condition number characterizes the TVFE for all applications. If we consider a waveguide or a cavity with metallic walls, FEM analysis based on a given TVFE leads to element matrix equations of the form [Ae]{x}e = A2[Be]{xe} (1) where each entry of [AC] is the integration of the dot product of the curl of two basis functions over the element and each entry of [Be] is the integration of the dot product of two basis functions over the element. Assembly of element equations leads to a global matrix equation of the form [A]{x} = A2[B]{} (2) The matrices [AC] and [A] are singular while [B2] and [B] are regular. In this work, we will perform FEM analysis of waveguides and cavities with metallic walls based on different TVFEs and use the condition number of [B] as an indicator of the matrix contitioning that the different TVFEs lead to. We note that the condition number used in this work is the ratio of the maximum to the minimum eigenvalue. Consider a rectangular waveguide of normalized dimensions 1 x 0.5 discretized into 64 triangular elements and a rectangular cavity of normalized dimensions 1 x 0.75 x 0.5 discretized into 130 tetrahedral elements. For TE as well as TM formulation (waveguide) and E as well as H formulation (cavity), Table 2 and Table 3 give the condition number of the global FEM matrix [B] for unnormalized as well as normalized vector basis functions. The higher order TVFEs are seen to lead to much larger condition numbers than the lowest order TVFE. For interpolatory TVFEs, the one by Graglia et al. leads to better conditioned matrices than the one by Peterson (waveguide) / Savage and Peterson (cavity). For hierarchical TVFEs, the one by Andersen and Volakis leads to better conditioned matrices than the one by Webb and Forghani. The hierarchical TVFE by Andersen and Volakis leads to better conditioned matrices than the interpolatory TVFE by Peterson (waveguide) / Savage and Peterson (cavity), especially when the vector basis functions are normalized. The interpolatory TVFE by Graglia et al. leads to better conditioned matrices than the hierarchical TVFE by Andersen and Volakis. We note that all conclusions for the global matrix [B] in (2) also hold for the individual element matrices [Be] in (1). Conclusions A study was presented that examined thoroughly the inter-relationships between condition numbers of FEM matrices based on various interpolatory and hierarchical mixed-order TVFEs. The validity of the generally accepted premise that interpolatory higher order TVFEs lead to better conditioned matrices than hierarchical higher order TVFEs was found to be questionable. This is of importance wlheu performing selective field expansion using hlierarchical higher order TVFEs.