Triangular prisms for edge-based vector finite element analysis T. Ozdemir and J. L. Volakis Radiation Laboratory Department of Electrical Engineering and Computer Science University of Michigan Ann Arbor, Michigan 48109-2122 March 22, 1995 Abstract This paper deals with the derivation and validation of edge-based shape functions for the distorted triangular-prism. Although. the tetrahedron is often the element of choice for volume tessellation, mesh generation using tetrahedra is cumbersome and CPU intensive. On the other hand, the distorted triangular prism allows for meshes which are unstructured in two dimensions and structured in the third dimension. This leads to substantial simplifications in the meshing algorithm and many printed antenna and microwave circuit geometries can be easily tessellated using such a mesh. The new edgebased shape functions presented in this paper are validated by computing the eigenvalues of three different cavities (rectangular, cylindrical and pie-shell). 1

1 Introduction The brick and tetrahedron are popular elements in finite element analysis of electromagnetic problems. The first is indeed attractive because of its simplicity in constructing volume meshes whereas the tetrahedron is a highly adaptable, fail-safe element. It is often the element of choice for three dimensional (3D) meshing but requires sophisticated and CPU-intensive meshing packages. The distorted prism (see Figure 1) is another volume element which provides a compromise between the adaptibility of the tetrahedron and the simplicity of the brick. Basically, the distorted prism allows for unstructured meshing (free-meshing) on a surface and structured meshing in the third dimension. An approach for growing prismatic meshes is illustrated in Figure 2 and most volumetric regions in antenna and microwave circuit analysis can be readily tessellated using such a mesh. As seen, once the surface mesh at the different surface levels is constructed, the prismatic mesh can be generated by simply connecting the nodes between the adjacent mesh surfaces. This avoids use of CPU-intensive volumetric meshing packages and in many cases the user/analyst can construct the mesh without even resorting to a surface meshing package. Examples include printed circuits and antennas on planar platforms. Moreover, because of their triangular cross-section, the prisms overcome modeling difficulties associated with bricks at corners formed by planes or edges intersecting at small angles. A special case of the distorted prism is the right prism which is characterized by the right angles formed between the vertical arms and the triangular faces [1]. The top and bottom faces of the right prism are necessarily parallel and equal, restricting them to a limited range of applications, namely, geometries with planar surfaces. In contrast, the distorted triangular prism is almost as adaptable as the tetrahedron with the exception of cone-tips which are not likely to occur in printed antenna and microwave circuit configurations. In this paper, we introduce edge-based basis functions for the most general distorted prisms. These prisms have non-parallel triangular faces and each of their three vertical edges can be arbitrarily oriented. In the following, we first present the edge-based shape functions and then proceed with the derivation of the finite element matrix. Eigenvalue, radiation and scattering computations are included and used to validate the new edge-based basis functions for prismatic elements. 2

zv 3 Y 7 x Figure 1: The distorted triangular prism shown with the directions of the edge vectors. Surface normal urface geometric fidelity Distorted triangular prisms stacked up Figure 2: Illustration of an approach for constructing prismatic meshes (unstructured surface mesh and structured mesh along the third dimension)

Perspective view Top view Interior point (x,y,z) Triangular cross-section containing the interior point Figure 3: Variation of the triangular cross-section along the length of the prism 2 Vector Edge-based Basis functions Consider the distorted prism shown in Figure 3. The prism's top and bottom triangular faces are not necessarily parallel to each other and the three vertical arms are not perpendicular to the triangular faces. On our way to specifying a set of shape functions for the prism, we proceed with the identification of a triangular cross-section of the prism which can be uniquely defined given a point (x, y, z) interior to the prism. The cross-section contains this point as illustrated in Figure 3. A way to specify the nodes of such a triangular cross-section is to use the parametric representation ri = rib + s(rit - rib), i = 1,2, 3, (1) where r2 are the nodal position vectors of the triangular cross-section (see Figure 4(a)) whereas rit and rib are associated with the top and bottom triangular faces, respectively, and 0 < s < 1. Thus, for s = 0 ri specify the nodes of the bottom triangle and for s = 1 ri reduce to the nodes of the top triangle. Since s assumes values between these two limits, the cross-section sweeps the entire volume of the prism. Clearly, each (x, y, z) point belongs to one of the triangles specified by ri and consequently there is a unique correspondence between s and a given point in the prism's volume. 4

(X. z 3 (x!' 1 Z d A d ~((x z Y Z- d.' * (xyz) Ac1 " hi (x,.y,z), * - -- (x,.y,.z - h h d -- cross-section g- * oC l (x,y lb 1 —.a h, (X,y bZ --- (X"Y-b'Z 2Y Z' 2 (a) (b) Figure 4: (a) Nodal coordinates, (b) triangular cross-section with the local coordinates 4 and rI. To determine s in terms of (x,y,z), we first note that (1) represents nine equations with ten unknowns (s and the coordinates of the three nodal points). Thus, additional equations are needed and these are the equations defining the plane coinciding with the triangle, viz. Xi Y~ Zi i + Yi + i +1 0, i= 1,2,3 (2) + - + - +1 = (3) a b c in which the constants a, b, and c correspond to the locations at which the plane intersects the x,y, and z axes, respectively. Substituting in (2) the values of xi, yi and zi obtained from (1) yields three equations involving four unknowns a, b, c and s. Solving these three linear equations for 1, I and c in terms of s and substituting them in (3), we end up with the third order polynomial a3s3 + a2s2 + als + ao = 0 (4) whose real root is the desired value of s. The expressions for the coefficients ai can be given in closed form in terms of x, y, and z, but are rather lengthy and bear no significance on the rest of the analysis. Therefore, they have been omitted but can easily be derived as outlined above. The appropriate solution of (4) in terms of ao,1,2,3 is [2] S = (S1 + 2) a2 (5) 3a3

where S1 = [r+ (q3+r2')~] (6) s2 = [r-(q3+r2)~]2 (7) r = 6 a(a2 - 3aoa3) - -(a2/a3)3 (8) 1 1 2 q -al -a2 (9) and this completely specifies s in terms of x, y, and z. We next proceed with the derivation of the basis functions. We choose to represent the field variation across the triangular cross-section (defined by s) using the two-dimensional Whitney form [3]. A simple linear variation will be assumed along the length of the prism. Specifically, the vector basis functions for the top triangle edges can be expressed as N1 = di (L2VL3-L3VL2)s N2 = d2 (L3VL1L1 VL3 )s (10) N3 = d3 (L1VL2- L2VL1 )s and correspondingly those for the bottom triangle edges will be M1 = di (L2VL3-L3VL2)(1 - s) M2 - d2 ( L3VL - LVL3 )(1 -s) (11) M3 = d3 (LVL2- L2VL )(1 - s). The subscripts in these expressions identify the edge numbers as shown in Figure 1 and the distance parameters di are equal to the side lengths of the triangular cross-section containing the observation point (see Figure 4(b)). Also, hL L(2, 77) cosa3 s 7n3 (12) h2.h2 cos a2 sin a2 L3(, ) C =+ 77. h3 h3

2X2 as a function of x,, and z. node heights from the opposite side. The variables, and r/ represent the 1'' U ith edge whereas Ni * j -Mi * ej - 0 for i 5 j. Their vector character 2 " (a) (b) (C) Figure 5: (a) Vector map of nd2 or M2 (b) Vector map of K (c) Variation of unode as a function of x, y, and z. are the usual two dimensional scalar node-based basis functions [4] for the same triangle with shap denoting the interior vertex angles and h being the node heights from the opposite side. The variables l and r7 represent the local coordinates and are illustrated in Figure 4(b). As required with all edge-based shape functions, N defi and Mi * i have unity amplitude on the ith edge whereas N * = * of = 0 for i j. Their vector character is depicted in Figure 5(a) and as seen they simply "curl" around the node opposite to the edge on which their tangential components become unity. It remains to define'the shape functions for the three vertical edges and we chose to express these by the linear representations K2 V,77) = n(~,r)L2((,,) (13) As before, Li are the node-based shape functions defined in (12) and a pictorial description of K1 is found in Figure 5(b). Of particular importance in (13) is the unit vector V'((, 77). It is a linear weighting of the unit vectors v1, v2 and V3 associated with the vertical arms (see Figure 5(c)), and is given 7

by 3'(=, Z) = E L~(~, r)ii.- (14) i=l This particular choice of vi minimizes tangential field discontinuity across the faces. Contrary to the rectangular brick, tetrahedron and right prism. the edge-based vector basis functions for the distorted prism do not ensure tangential field continuity across the faces, nor are they divergenceless. The same is true for the curvilinear bricks [5] and, in the case of the distorted prism, the field discontinuity and divergence increase with surface curvature (see Figure 2). This is certainly an undesirable feature but in most cases (particularly when sampling at 15 or so points per linear wavelength and the radii of curvature are greater than a wavelength), the angular deviation of the vertical arms is quite small. Consequently, for all practical purposes the field discontinuity and divergence are negligible. 3 Eigenvalue Computation In this section, we examine the validity of the presented edge-based functions. Specifically, we consider the eigenvalues of three different cavities using the edge-based distorted prism as the tessellation element. We begin by first deriving the matrix elements following Galerkin's testing. The weighted residuals of the vector wave equation are Ni * (V x V x E - kE)dV i= = 1,2,-3, (15) I IMi (V x V x E-kE)dv, i= = 1,2,3, (16) I II/|Ki * (V x V x E - k2 E)dV = 0, i = 1,2, 3, (17) in which Ni, Mi and Ki comprise the nine edge-based vector basis functions defined in the previous section and E is the electric field vector. The matrix equations are generated by introducing the representation 3 E E [EiNN(r) + EiMMi(r) + EiKKi(r) ] (18) i=l 8

swhere LN. Ei.Ei and EitK are the expansion coefficients, and correspond to the average amplitudes of the field vector at the ith edge. Substituting (18) into (15)-(17), and invoking the divergence theorem yields the element equations 3 3 E EjN[NNCij - k2NN;Dij] + Z Ejm,[N IC, j - koNa Diij] + j=l j=1 3 Z EjK[NKC0j - k2NNKDtij] = 0 (19) j=l 3 3 Z EJN[MNC j - k 2MNDlOj] + E EjM[MMC j - ko2MMDZj] + j=1 j=1 3 EjK[MKCi, - k2MKDij] = 0 (20) j=l 3 3 E EjN[KNCi, - k2KNDij] + E Ej,[KMC,j - k2KMDj] + j=1 j=1 3 EEjK[KKCij - ko2KKDij] = 0 i = 1,2,3. (21) j=l where NNCie = JI (V X Ni) (V x Ne)dV (22) NMCit = / (V x Ni) * (V x Me)dV (23) NKCit = III(v x Ni) * (V x Ke)dV (24) MMCit = I JI(V X Mi) * (V X Me)dV (25) MKCit = I II(V X Mi) * (V x Ke)dV (26) KKCie = I (V x K ).(V x Ke)dV (27) NNDie = IvNi N.dV (28) NMDie = JJ Ni MedV (29)

Region I (O<z<z,) Region2(z,<z<z,) Region3(z <z<z,) Triangular cross-section Quadrilateral cross-section Triangular cross-section'z / Ig / i% NKD/ = | Nt KtdV (30) MMDit = || Mi MjdV (31) MKDite = || Mi eK dV (32) KKD = |||Ki KedV (33) The above integrals must be evaluated numerically and to integrate over the volume, the prism is divided into three regions (see Figure 6) of triangular or quadrilateral cross-section. In the case of the right prism, the integrals have closed form expressions and are given in the Appendix. Upon assembly of (19)-(21) and boundary condition enforcement, we obtain the generalized eigenvalue system 2[A]x} - = -k[B]x (34) 10 10

lcm Mode k. cm % Error 0.5cm (Exact) Prism Brick Tetra. - - as TE l0ol 5.236 -0.99 -1.36 0.44 -'. _-' 0TM 0o 7.025 -4.44 -2.23 0.70 0.75cm TEol1 7.53 1 0.07 -2.58 1.00 v~, -..TE 201 7.531 -0.25 -3.13 -0.56 - (Actual mesh) TM,,, 8.179 -0.31 -2.09 2.29 (a) (b) Figure 7: (a) Rectangular cavity discretized using right triangular prisms, (b) eigenvalues for the air-filled cavity. Example 1: Rectangular Cavity The first example is the rectangular cavity shown in Figure 7(a). Results based on brick, tetrahedron [6] and prism discretizations are given in Figure 7(b), where they are compared to the exact values. We observe that, overall, the data based on the prism are better than those based on the brick and, at least, as good as those associated with the tetrahedron. Example 2: Circular Cavity The tessellation of the circular cavity (drum) is shown in Figure 8(a). The mesh was constructed by first creating a surface grid of triangles at the base of the drum. The prisms were then generated by growing the mesh along the axis of the cylinder. In this case, the resulting volume element is the right prism and the eigenvalues are given in Figure 8(b) along with the exact. Except for the higher order mode TE211, the remaining eigenvalues were computed to within three percent of the exact. Example 3: Pie-shell Cavity The third example is a pie-shell sector as shown in Figure 9(a). It is obtained by bending the rectangular cavity considered earlier and the resulting volume element is the distorted prism with a vertical arm angular deviation of five degrees. The computed and exact eigenvalues for the first five domi

Radius = lcm (Actual mesh) Mode k, c-l % Error (Exact) (Computed): *, 1 -- — ~~-____~ TMolo 2.405 1.29 II _ -- = - TEiii 3.640 2.17 4 i _TMl0o 3.830 -2.90 t 1t. r --. — -- r TMoll 3.955 0.81 - - ~ Y TE211 4.380 -8.97 (a) (b) Figure 8: (a) Circular cavity discretized using right triangular prisms, (b) eigenvalues for the air-filled cavity (Actual mesh) I Mode k, crm- % Error <' ] (Exact) (Computed) TMi,, 4.693 1.52 0,.75cm TM210 6.009 i 0.95 T rE Ill 6.640 4.96 01cm / 21 TE211 7.513 -1.11 O(700) 1 TE Oil 7.579 j 1.71 (a) (b) Figure 9: (a) Pie-shell cavity discretized using distorted triangular prisms, (b) eigenvalues for the air-filled cavity 12

Side view: A Z Top view: A x x FEM mesh ABC surface ABC surface ~,,1 1,, 0 _-_____________ PEC 0.3 o, t _ _ _ _ _ _ _ _ _ _ _ 10.3 X Figure 10: The geometrical set up for finite element-ABC formulation of radiation and scattering from cavities in infinite ground planes. nant modes are given in Figure 9(b), and these testify to the accuracy of the distorted prisms in modeling curved geometries. 4 Radiation and scattering from cavity backed structures As the last set of tests, we have looked at radiation and scattering from cavity backed radiators recessed in an infinite planar or cylindrical metallic platform. The general set up geometry is shown in Figure 10. As it is seen, a three dimensional cavity is discretized using the triangular prisms. The finite element mesh is extended outside the cavity in all directions about a fraction of the wave length and truncated using the socond order ABC reported in [7]. Example 1: Radiation from dipoles in a cavity on planar platform We have first looked at the radiation from a pair of current elements residing inside the cavity; one horizontally and one vertically oriented as shown in Figure 11. Radiation pattern is shown in Figures 12a and b. The solid and dotted line results were obtained using brick elements and closing the finite element domain at the aperture of the cavity via boundary integral (BI) method, which is an exact formulation and is explained in [8]. The 13

0.3 xk (Side view) y-directed element Actual gridding z-directed element 0.525 7 rI (Top view) 0.75 Xo Figure 1 1: Locations of current elements with respect to the gridding inside the cavity triangular prism results were obtained in conjunction with the absorbing boundary condition (ABC) to truncate the mesh 0.3A0 away from the cavity aperture in all directions. As the plots clearly demonstrate, the triangular prisms together with the ABC perform well. Example 2: Plane wave scattering from a cavity on planar platform Next, we have looked at the plane wave scattering from a larger cavity whose dimensions in x, y, and z directions are 1.2Ao, 0.75A0, and 0.3A0, respectively. The geometrical set up for the finite element-ABC formulation is the same as in the previous radiation problem which is illustrated in Figure 10. Figures 13 and 14 show bistatic and backscattering radar cross sections (RCS), respectively, of the cavity. It is seen that the FEM-ABC method using triangular prism elements exhibit sufficient accuracy in computing plane wave scattering from cavities in a ground plane. Example 3: Radiation from a patch antenna on cylindrical platform As the last test, we have looked at the radiation from a cavity backed patch antenna on a circular-cylindrical platform shown in Figure 15(a). The distorted triangular prisms used to discretize the cavity in Figure 9 have been used here also to discretize the computational domain. Finite element mesh has been terminated 0.3 wavelengths away from the cavity aperture by the first order cylindrical ABC derived in [7]. Figure 15(b) shows the radiation pattern computed as compared to the finite element-boundary integral 14

(a) (b) -10 -10 -20 -20 Brick-B0 (-poll -- Brick- Bric- -pol) Trang~ular prism-ABC * Tnranlula prism-ABC -30 -30 0 30 60 90 0 30 60 90 Observation angle 0 (deg) Observation angle (deg) Figure 12: Radiation pattern of two current elements located inside the cavity; (a) -=45~; (b) 0=45~ (a) (b) 10.'.,.,..10.. N rs N 1 ~ L - Brick-BI (B-pol) -30 Brick-BI (-pol) |30 - Brick-B! (+-pol) -3 - - B30 - - Brick-B (s-pol)) e prism-ABC Triangular prism-ABC (S-pol) * Triangle prism-ABC (S-pol) -40 T-40 - ~ Triangular prism-ABC (t-pol) * Trinalie prism-ABC (s-pol) -50 -50 0 30 60 90 0 30 60 90 Observation angle e (deg) Observation angle t (deg) Figure 13: Bistatic RCS of a cavity of dimensions 1.2 O x 0.75 o x 0.3 3o for a plane wave with its electric field vector making an angle of 45 witPthe unit vector in 0 direction incident on the cavity aperutre at oinC=45 o 4inc=30 (a) 4=60, (b) 0=45~. 15

-20 - - -20 -20, — Bnck-EBI (-pol, -30 - - Bnck-BI &-pol) -30 B -po30 I Bric-Bl i-pot) Triangle prism-ABC (0-pol) Triangle prism-ABC (s-poa) -40 * Triangle pnsm-ABC (*-pol) 40, ~. \,*Triangle pnsm-ABC (s-pol) -50 -50 0 30 60 90 0 30 60 90 Observation angle o (deg) Observation angle 0 (deg) (a) (b) Figure 14: Backscatter RCS of a cavity of dimensions 1.2ko x 0.75 Xo x 0.3Xo for an incident plane wave with its electric field vector making an angle of 300 with the unit vector in 0 direction; (a) 4=30~; (b) 0=600. solutions [9]. Again the results are excellent. 5 Conclusions The distorted prism is an indispensable tool for the analysis of many doubly curved antenna and microwave geometries. It provides simple volume meshing without compromising accuracy (see the comparison with the tetrahedron in Figure 7(b)). In this paper, we introduced a set of edge-based shape functions for the distorted prism and used them to generate the element equations. For validation purposes, the eigenvalues of three different cavities (rectangular, cylindrical and pie-shell) were computed and compared to the exact solutions. In almost all cases, the error was less than three percent. Radiation and scattering results are also excellent when compared to the finite element-boundary integral solutions. 16

r No,~ —.... -PEC infinite cylinder (Actual mesh) 0 o.05o _- Feed location- (a) Cavity wall o (o y Cr,;; 0'23'o_ 1' I, Patch II 0.2X0___ ABC at 0.3 oX 0=30 o 0=600 - 40.-, —,-,-,-,., —.. -,., - - 30. -4030 F4 30 20 slo % 20 1 0 e —-.I. I'i,,~l 10 -o( 0X0 -20 ~~~~~~~0-20 ~~~~~0 ~~0 -30 ~n -30 -40 -40 0 30 60 90 120 150 180 0 90 180 270 360 Observation angle 0 (deg) Observation angle ) (deg) FE-BI (9-Pol) * FE-ABC (0-Pol) -~- FE-BI ($-Pol) * FE-ABC ( —Pol) Figure 15: Radiation from a rectangular patch on circular-cylindrical platform: (a) geometry, (b) radiation pattern. 17

Appendix For the right prism, closed form expressions of the integrals in ('22)-(33) are possible. Referring to Figure Al, we have ddt( cos 3kn COS 3jm COS 3km El~C~e = - ( Xjm + - Xkn - jn, - c hkh, hjhm hkh, COS /jn 2 c2h1dl Xkm + hhkhmhsin3jk sin3mn ) hjhn 3 hjhkhmhn dide cos fkn cos Ojm Cos 3km ENMCe = (- XJm h-hm Xkn + Xhkhm + c hkh, hjhm hkhm cos/ j, 1 c2h1dl Xkm + hhkhhn sinjk sin3mn ) hj hn 3 hjhkhmhn ENKCie = hhld di (cos je cos 3ke 6 hjhe hkhe EMMCie = ENNCie EMKCie = -ENKCie hEldl COS Oil 2 hih1 dide COS /kn COS/ jm COS 3km COS_ Xkm 3 hk + hjh Xk- hkhm Xjn - Xkm) ENMDie = -ENNDit 2 ENKDie = 0 EMMDie = ENNDie EMKDie = 0 18

EKK IDi; = ci uw here,- 0 ifr=s wah + a, otherwise Xrs = 2 {wrs + hi [(Cota3 - Cotoa2)(rlsWr + 7qrWs) + 2(~sWr'+ &rWs)] + 2 3 h2 -1[3(cota3 - cota2)(77,sr +,r.s) + 2r71Sr(Cot2a2 - COta2COta3 + 12 cot2cr3) + 6,Ur]},a r, s = 1, 2,3, the set {i, j, k} takes on the value {1,2, 3}, {2, 3, 1} or {3, 1, 2}, and the same is true for the set {1, m, n}, W1 = 1 il -h- 1 -~0 sin ~a3 w2.=-O _2= h2 T= h2 =Co02 sin a2 w3=0 3 h3 ~3 h3 References [1] Sacks, Z., S. Mohan, N. Buris and J. F. Lee, "A prism finite element time domain method with automatic mesh generation for solving microwave cavities," IEEE APS Int. Symp. Digest, Vol. 3, pp. 2084-87, Seattle, Washington, June 19-24, 1994. [2] Abramowitz, M and I. A. Stegun, Handbook of Mathematical Functions, Dover Publications, Inc., New York, 1972. [3] Whitney, H., Geometric Integration Theory, Princeton University Press, Princeton, New Jersey, 1957. [4] Zienkiewics, O. C. and R. L. Taylor, The Finite Element Method, 4th ed., Vol. 1, McGraw-Hill, New York, 1989. 19

UNIVERSITY OF MICHIGAN 11111111111111111111111111111111 111111111111111111111111111 3 9015 03466 4121 3 N2 *M N ~~~M,~~~~~~~~~3 K 3A h3 K K di K2 (a) (b)d3 (a) (b) Figure Al: Right triangular prism: (a) Indexing of nodes and edges, (b) parameters defining the triangular cross-section. [5] Antilla, G. E. and N. G. Alexopoulos, "Scattering from complex threedimensional geometries by a curvilinear hybrid finite-element-integral equation approach," J. Opt. Soc. Ame. A, Vol. 11, No. 4, pp. 1445-57, April 1994. [6] Chatterjee, A., J. M. Jin and J. L. Volakis, "Computation of Cavity Resonances Using Edge-Based Finite Elements," IEEE Trans. Microwave Theory Tech., Vol. 40, No. 11, pp. 2106-2108, Nov. 1992. [7] Chatterjee, A., "Investigation of finite element-ABC methods for electromagnetic field simulation", Ph.D. Thesis, Chap. VI, EECS Dept., Univ. of Michigan, Ann Arbor, MI 48109-2122, September 1994. [8] Jin, J. -M. and J. L. Volakis, "A finite element-boundary integral formulation for scattering by three-dimensional cavity-backed apertures," IEEE Trans. Antenn. Prop., vol. 39, No. 1, pp. 97-104, January 1991. [9] Kempel, L. C., J. L. Volakis and R.J. Silva, "Radiation by cavity-backed antennas on a circular cylinder," IEE Proceedings, 1994. 20