#036307-1-T Focused applications software for Ferrite Antenna Analysis SBIR Phase II MRC/UM/EMA review held at the University of Michigan Ann Arbor December 17, 1997 36307-1-T = RL-2503

Agenda for Ferrite Antenna Review Project with MRC/UM December 17, 1997 Univ of Michigan 3400 EECS Building, Ann Arbor, MI 8-12am It is assumed Questions will be asked only at the end of each presentation. 8:15-8:40 Updates on PRISM (Volakis/Eibert/Ozdemir) * Spiral modeling and design using combined distorted PRISMs/BRICKs * Periodic version of PRISM (FSSPRISM) * Fast integral methods (AIM/FMM) performance in PRISM 8:40-9:10 Measurements and MRC Project Activity (Kempel/Trott) 9:10-9:20 Overall Project Status at UM (Volakis) 9:20-9:40 PRISM/TETRA Geometry Driver (Zhifang Li) * Description and Status Break 10:15-10:45 GMRES Solvers (Arik Brown) Eigenvalue Solver Issues (Lars Andersen/Volakis) 10:45-11:15 Tetrahedral Ferrite Code Implementation (A. Brown) * Formulation * Code Status

UPGRADES TO FEMA-PRISM FOR ARRAY AND BROADBAND ANTENNA SIMULATIONS WITH FAST INTEGRAL METHIDS T. Eibert, K. Sertel, T. Ozdemir, Y. Erdemli, J. L. Volakis Radiation Laboratory Dept. of Electrical Engin. and Computer Science University of Michigan Ann Arbor, MI 48109-2122 December 1997

H.PRISM FOR ANTENNA/ARRAY ANALYSIS _ _ _ _ _ L t X FEATURES/CAPABILITIES *Finite Element Antenna Analysis Code *Prismatic elements and Hexahedrals for mesh generation (retains geometrical generality in one dimension) *Boundary Integral or artificial absorber for mesh truncation *Automatic meshing for rectangular/circular patches, log-periodics, archimedean spirals, dipole/patch arrays. *Visualization of mesh and fields using interfaces with Matlab and UM's Meshview

PRISV CUTRENT FEATURES * Flexible (Materials and Geometry), Rigorous and Fast. * Incorporates fast solvers (AIM and FMM) * Planar and Curved Elements/Arrays * Periodic and non-periodic broadband element Frequency selective surface planar and non-planar inhomogeneities * Antenna fe Periodic Arrays and Frequency Selective Surfaces.. S::A 4 j i. S; y::;,3 l,, I:111& ' ' ' f 'i T.i Bandgap using periodic array features -Multilayered powerc tombining networks 2

unite ELIilLD~l~r neral versus F inite Eleiment~ A-1rtificial Absobers;e Matrix. {L1 - {tnc PM4L Deteriorate Covergence Dense *Numb eY oflInriiowns is reducedb a factor of or mnore leFor Narrowv slots Bimatrix is smlal *FF can be used to speed-up) comp utation of matrix vector prod ucts 3

Ct VI l' A E LO~ CD 0 t- Tcor~ (Dt -oVI >D. VQ H (D 4N U Jy1 VI - 1O (D L. -~ II iiI8 0)8 0 00 0 Sw E U)

jF-IELD BIEHA VIOR INSIDE TILE SLaOT 80 Zin 71-[j 0 C ARG=-94dB 0.8 GHz -80 L -100 0 Distance (cm) 100 80 cm 0 -5 Zin= 46+j2 Q =5 dB AR. =2 dB -80OL -1I0 1 0 0 Distance (cm) 100 80 CT5 -5p Zin 37-j4 Q G -8dB AR =1 dB3 3.0 GHz -80o L. -100 0 Distance (cm) 100 5

ORIGINAL AND IMPR O VED DESIGNS 10 8 6 4 2 c — cr c_5 0 -2 _FROM FEMA-PRISM -*^- COMPUTATION.. 4 turns - Measured - 4 turns - Simulated - 5 turns - Simulated - 6 turns - Simulated 2.5 3 HHz) r =3.38. 05 1 cm '1 1 O i 0.58cm -4 \* 7 m la* 4 7 -8v -? Original — 8 -10 o I 1 1.5 Fr e 2 u ency (G Lumped resistors (Dissipates same power) 'Resistive card Simulated (4,5,6 turns) Actual (4 tuLrn s) 6

FS round Plan or Radom -i j.1 I 1,,., iij, I. m T:!, 5 'm.... TI i. - - i....!;,.i Pzj c r~~N~. kI*I pw, - I~2 * 6- -diilII cro a.1-1 Ikailai ui tUt.:s *C I-, ",Ack.I. 'i t I I,i~b u 9']~ k s"I -k- *t s's I 'Ii i, 7

FSS PRISM Periodic Boundary Condition (FE) Unit Cell: Infinite Periodic Array DYP cot Y / 22 7 Periodic Free-Space Green's Function (BI) Unit Cell with Image Sources: 8 9 lo 2 7 1 14 3 19 20 II. DXP 8

*! Hybrid FEM-BI Formulation for Arrays * Periodic Boundary Condition (PBC) for FE-part Metallic backing (Antenna) BI PGF + Periodic Green's Function (PGF) for BIpart I PBC " PBC Metal * Metallic patches in all layers possible (Antenna, FSS) BI on top and bottom surfaces (Transmission, FSS) BI/PQF * Probe current feeds and plane wave / excitation __ _ PBC '...' PBC BI/PGF * Lumped impedances / resistive sheets _ 9

Array Sinmuiation Results Vt'''''_w k r 1 'k - "L 1'. If U x '' 'C'" Z, Fss "Vc' v- | 1 7....... Power Reflection Coefficient 1if - FSS-Erick: * 2x2 Array 3x3 Airray EIGER: Infinite Array - Dipole O 0.4O cm -. Antenna hFSS- 0.4:0cmt l<-af-^ FSS-i I,S' Frequency, GHz h2- 0.25cmT (:a FSS-2~t Dipole-2:: (0.25x0.6cm2) 1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 n.........................................,......................................... 4,........................ '.........i.,f........... ~ ~ ~~~~~~~ii~~......... i..........i.,.......... - - ----! /{ I., h1= 0.l5cm* 4 4 0-9 i514 K' - cm L= 0.95 cm v........ 0 10 20 30 40 50 60 SCAN ANGLE (DEGREES).. I, 70 80 90 10

FSS-PRIS NonComesuat VaiaIo in m inm- m - r --- —--— ' I I _______ _______________________________________________ I _______________________________________________ _______________________________________________ I __________ I I I 1 I -I - I - I ~ I I Unil cell 9/ 9 9 9 g 27 'm -rn —rn * gtljm 9 1 8 1 8 1 8 12 g 6 14 1 4 14 floi U I c1 4.-J C', 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 r r=2.3b j.... I r w 0 1000 2000 3000 4000:5000 6000 7000 8000 Frequency InGQHz per frequency nt Workstation CPU abo ut 1/ hour 11

CCl Coo) 0, PO0 S,o 0 00 00

* PUand Memory Requirement Reductions for O'nL5) MemorySorage *_MMAJM vs. Conventional Metho 0n2) v s. On15) FLP 0n2) vs Onl.5) Memory Locations 13

, ADAPTIV INTEGAL MEHDAIM) 0 riigina1 Expansion J Ir' Nv I~f~(rI ~kn= I8( -Am4[2g An~ 14

ts bCIA CA 16.0 4c) ci An, n ci r~l"~i mLn d) dd Q1 cncr I ag ~~ I I I %"* **.*1 Ij I I Qci I 0 0 0 0 - - - - - - - - -- - - - - - - - - -: - - - - - - - - - - - - - - -- - - - -

CPU Times 20x20 Grid Run Data DOFs - 7435(vol) + 2x1240(Surface) Memory: with AIM: 588k elements (7MB) w/o AIM: CPU(IBM RS6000/590): 3000k elements (36MB) AIM additi FILL=3 sec(FEM)+685 sec(BI)+5 sec(AIM n BCG Solver= 227 Iterations (228sec) BCG Solver=234 Iterations (166 sec) I Conventional 13 sec(FEM)+ 1 400sec(BI) | AIM reduced memory by a factor of. CPU by a factor of 50% W f I 16

CPU Times 32x32 GrdRnDt DO s = 542vl + 2x3 12(urface) Memory: wit AIM 213 elements(36B wlvo AIM: 20,000 elements(20B CPU(66Mz PC): AIM FIL sec(FM+ 1 0sec(BI+ 1Osec(I BCG Solver=- 41terations (01see) Total CPU=3028sec (50mm.) AI redu ced mlemo ry byI afctr o 6.5 Souto wl AI i toIxpnive 17

DO s = 6700vol + 2x586(urface) Mem-oryr: with AIM 4672k elemnents (5Mlbtes pere element) w/vvo AIM: 71,000 elements(82B AI IL3sec (E )30sec(I +22sec(i finll BICG Solver-67 Iter. (573sec=1.4 hours) mtemoyDya faIcto of 15 SoutIon w/o AI is to expensive 18

PRISM with Fast Multipole Method FE-BI Revisited Iterative Solution of the Resulting System Sparse Matrix Dense At cross int A1 o ols E' nt1 Tcross bound bound ound t \^ g X XndA Ao E UJ fl GL LE J [. -- -Z~.= —~F~/1 -~ Plain FE-BI Strategy 19

FMM's Ideare for (/ IvAatrix-Vector Product - 0 Plain FE-BI Strategy 20

Ci) 0 rm IT0,M z ml 104r CA E -x q -J IA (H T VI 4- 7 -t-4 V"J\, \1 4" ',,I ITT, TN "111,41) VIII) Tq \J T-I 1" \l" 'I) \P TIP, '11 d\ 1\ W l\ T" 'ITT '\I\ W E E Q 'Tt Itt Wdfit — MD- -ml4F- 0 -I 0 CY) (0. - CC) LO 000 -0 0 LNAN C\J, ICKo E. ON ei ------ '1*11 I i i II IN I C-) I 0

0 ------------- -------------------------......................................... 0 0 0 I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I E I I VI CIA 0-0...... l~ 0I o C14 0 0 0 0 0 0 ~r (I) (C4 T U.PI -- I d — - - ---- - - 0 qbC. t0~ C. 0.w C. 0 0 40 40 0 No 0 0 0 so to 0 N 0l V # i N

MRC Presentation Leo Kempel Keith Trott

PRISM/TETRA Geometry Driver Zhifang Li and John L. Volakis Radiation Laboratory Department of Electrical Engineering and Computer Science University of Michigan Ann Arbor, MI 48109-2122, USA

Outline * Features of the Driver * Mesh Generating Steps * Mesh Formats and Examples o Surface Mesh o Volume Mesh

Features of the Driver * Capable of meshing rectangular/circular/log-periodic patches/slots * Non-uniform surface mesh available * Volume elements could be right/distorted prisms/tetrahedrons * Mesh truncation wih boundary integral (BI) or artificial absorber (AA) * Mesh generated by FORTRAN code and can be viewed in MatLab or I-DEAS

Mesh Generating Steps a Surface Mesh I rI Surface Mesh Grows Above and Below I Volume Mesh Using Tetras I Prisms Divided into Tetras i| Volume Mesh Using Prisms

Mesh Termination Options * Boundary integral (BI) * Artificial Absorber (AA) Antenna A Cavity-backed: Microstrip: Ground plane (PEC) Cavity Substrate Ground plane (PEC) Modeling of the above two configurations with Bi and AA terminations: (Computation space is circled with dashed lines) Termination echnique ch e Boundary Integral (BI) Artificial Absorber (AA) Configuration Absrober Cavity-backed -.. Absrober Microstrip Does not apply gm iii

Surface Mesh * Structured mesher o Input file format for structured mesher Line #1: I I I t 1 = Boundary Integral (BI), 0 = Artificial Absorber (AA), 1 = Log-periodic, 2 = Circular, 3 = Rectangular E 1 = Printed, 0 = Slot R: Real, I: Integer, - - - - - -1 Lines #2-4 depend on the entries on Line #1 (see the following pages). o Output files o SurfMesh - connectivity table of surface nodes and triangles o MatLab files for mesh viewing on all platforms o Could also viewed by MeshView on PC * Unstructured mesher: will not be discussed in details here

Surface Mesh Examples 1 - Format BI termination - Rectangular Meshin Line#1: 1 31 or 1 30 Line #2: R R I I Ax1 Ayl NxA NyA Ax 2 Ay2 NxC NyC (sampling cell size in x-direction) (sampling cell size in y-direction) (number of antenna cells in x-direction) (number of antenna cells in y-direction) (sampling cell size in x-direction) (sampling cell size in y-direction) (# of cells between the antenna and the cavity wall in x-direction) (# of cells between the antenna and the cavity wall in y-direction) Line #3: R R I I./.. I/ / / /,/1I,//i/I / './. z /I/ I\ Ayl Ay2 Ax2 Axl For this example: NxA=10, NyA=2, NxC=2, NyC=2

Surface Mesh Examples 1 - View with MatLab MeshIn - 9 BI termination / Rectangular - 2 Example run - I Line #1: 1 3 1 (printed) Line #2: 1 1 12 2 Line #3: 3 4 After running the mesher, one can view the mesh, number the nodes and assess the mesh quality using MatLab as shown below: MATLAB INTERFACE 1. Display the mesh >> etup??? Undefined funct ion or variable etup. 2. Number the nodes and zoom in

Surface Mesh Examples 2- Format Bl termination - Circular Mesh Format Line#1: 1 21 or 1 20 Line #2: R I AR1 (radial thickness of antenna rings) - Na (# of antenna rings) Line #3: R I - AR2 (radial thickness of rings between the antenna and cavity ) Nc (# of rings between the antenna and the cavity wall). -- ' " ""-......'.' ' -.................... -—: —,. -.....;. - -...... -—. *....-. " -............ '................-.' j: / ' ' / 5 \ _ / j r...,...,...-.'.:. ~..........",..... \ ";:;-.,0 \.,...,........i/i, \.. \....' "'.~'.. '....-.::.........:......:..:.... A.. For this example, Na=5, Nc=3

Surface Mesh Examples 2 - View with MatLab BI termination - Circular Mesh Viewing Line #1: 1 2 0 (slot antenna) Line #2: 1 5 Line #3: 2 After running the mesher, one can view the mesh, number the nodes and assess the mesh quality using MatLab as shown below: MATLAB INTERFACE 1. Display the mesh 2. Number the nodes and zoom in commanas 0or more lnrormai ion: help, whatsnew, info, subscribe

Mesh Grows to Whole Volume Using Prisms * A single prism element 3 r 2 b6 5 4 * Mesh grows up- and down-ward for planar/non-planar substrates/superstrates Non-planar substrate Important! Works with planar platforms with Boundary Integral termination and circular cavities. Platform is asssumed to be located at z=O. Substrate accupies the z<O space. z Patch z=O r / z>O z<O - cut 0o - Locations for which the coordinates (z,r) are provided in "MainIn" file. Example Geometry: This example's result is provided on the distribution disc under the directory called "Demo-BI-NP". It is the same circular antenna modeled in Demo-AA-1 but with a non-planar substrate. (-0.45, 2.6) (-0.4, 1.2) (-0.4, 0) I / / = --- e=1.6 ------ g=1.2 (-.8, ) (-0.8, 0) \ (-0.8, 1.2) (-0.95, 2.6)

Prisms Divided to Tetrahedrons 2 5 4 A single prism is divided into 3 tetrahedrons 2 2 1 7 II / 6\ 6 -- 5 4 * Keep all the nodes generated in prism mesh * Add three diagonal edges in each prism * Each prism is divided into 3 tetrahedrons * Construct new connectivity tables for tetra mesh

Volume Mesh Input File Format MainIn (Part 1) 0 D 1 = Boundary Integral (BI) termination, 0 = Artificial Absorber (AA) termination * I I = Adaptive frequency sweep, 2 = Uniform frequency sweep I * R Desired amount of increment in the input impedance (in Ohms) the previous line has I. the entrv "1 * I I 1 = Printed, 0 = Slot, > I = Compute memory allocation, 0 = Analysis * I I I I I I I I I I, \ I. Io R R I? - # of substrate layers I = all substrate layers are identical and planar, 0 = non-uniform but planar, 2 = Non-uniform and non-planar i # of superstrate layers ( enter zero for no superstrate) 1 I = all superstrate layers have the same thickness and material parameters, 0 = otherwise C C C C } Ordered from the bottom of the cavity up, each row corresponds to a substrate layer. Only one row: J is needed if all layers are identical ( row has the info for a single layer ). C C I w'I I Thickness of the layer I- - Relative permittivity of the layer " I Relative permeability of the layer I C --- —-------------------------------------------------------------------- ---------- ---------------------------------------------- / * I C C \: p 'ubstrate glatform layers conforms to I I I I I I I I I I I I I I I I I I I I I I I 0 0 0 -- -- I + Number of linear segments - Relative permittivity of the layer -- Relative permeability of the layer R R R R..R R z- r coordinates of the segment junctions (number of pairs equals to the first entry on the previous line) I C C R R R R..R R I C C R R R R..R R Ordered from bottom up (first entry for the bottom ground plane. No entry for the aperture, which is always planar) \ substrate layers do not necessarily conform to platform, which is planar Works with BI option only! For this section, enter "2 "for the second entry on the previous line. / I Geometric Parameters

Volume Mesh Input File Format (Continued) MainIn (Part 2) -- - - - - - - - - - - - -- - - - - - - - - - - ------------------- -- / * R C C, * R C C Same as above but for the superstrate. Ordered from the antenna surface up (first row corresponds to the layer just above the antenna surface). ', R C C./ - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - ---— __ _ -_ -_- - - - - - - - _ —.- - - - - - - - -^ I - - - # of probe feeds present only if there is superstrate I I I C * I I I C and works with the AA option only 1 II C Each row corresponds to a probe feed I I I C. Surface node number #1 - Surface node number #2 } Probe current flows from node #1 to node #2. -> Layer # (layer within which the normally oriented probe is located, or the layer at top of which the laterally oriented probe is located). Entry can be positive or negative and increase away from the surface of the antenna with zero corresponding to the layer immediately below the antenna. ~ Complex amplitude of the probe current * R R R I I R I I -. Starting frequency in GHz - Final frequency in GHz - Increment frequency in GHz - - Frequency run to save (1 = save the first freq. run, 2 = next frequency, etc.) - -- Tolerance( - 0.01 ) - 1 = monitor convergence (dump residual error at each iter.), 0 = otherwise - 1 = compute element matrices assuming distorted prism (must for doubly-cureved or non-planar substrates), 2 = compute assuming right prism (saves on CPU time for planar platform with planar cavity-bottom) - Maximum number of iterations / - - 1 = Read in user specified termination parameters (given in the following row), 0 = code will figure out the optimum parameters (this is the safe course if one is not familiar with the artificial absorber termination). R I I C 0 RI C. - Thickness of one layer (all layers have the same thickness) For AA only! —. Total number of layers from the top of the outer-most superstrate layer to the termination boundary (Ignored by BI) - - Number of absorber layers X, --- Relative permittivity of the absorbing layers. -- -- -- -- -- -- - ---- - -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- ------ Electrical Parameters

Volume Mesh View with I-DEAS * Prism Mesh * Tetra Mesh

Format of Major Output Files from Volume Mesh * Edge Table 'EDGY' edge number x-coordinate y-coordinate z-coordinate * Edge Connectivity Table 'EGLOB' ele no. local edge no. node 1 node 2 material no.| * Element Connectivity Table 'ELNO' element no. node 1 node 2 node 3 node 4 material no. * Node Table 'ENODDY' node number x-coordinate y-coordinate z-coordinate * PEC Surface Table 'ESURFC' global edge number on PEC * Aperture Edge Table 'ESURFD' ele no. on aperture edge 1 and two nodes...

Ferrite Patch Antennas- Phase II December 1997 Update Students: Arik Brown, Zhifang Li, Lars Andersen PI: John Volakis Radiation Laboratory University of Michigan Ann Arbor, MI 48109-2122 volakis@umich.edu Phone: 313-764-0500 Fax:-313-647-2106 HomePage: http://wwwv-personal.engin.umich.edu/-volakis/ University of Michigan 1

(Ni — Co Co Co am CU LU) 0=j LL cmU) a, LCU U-C C,) a) C,) a0) 0> a/) c,)I LV >LM >0 4a, CUV 00 (1) CU) am) CU CU a)~ C,) ENO;,No FM F=;MO t TI 1900 CM) ct Co) 0: Q) 0 0

() il r% WIO ' I,i _If p C U - -. 'C i II II II II Ii ' i ii I ~7 C~ d i ~ I,,~ O~~~~ zz! I I ' ~ I' C ~ ~ Q

0 *0 I I I 0 11 V4 1 0 1 0 1 0 2. Backg round, r Agrees to wvithin 1.8 data -ror = 6 AMmwb -a dw m m Wm - - -.L. m m*m bmi,. w i I. University of Michi-gan Univrsit of ichian4

w 2 m 0 - a w a -M - - W- a 5 a 0 w 9 v 0 v a (I) a C') ~2i aca n icsI AM0O 0i Un o 16 - IL W Cl)l Uin A WI WI PMPL=Oqmg.i

tIe ' " (:D I~ o 03 a? - "CT 0j -7 0 0 ) a sc -:~ O T OD 4-aS 1)_ > _ I 4.) ^~m~- ' |2 ( >_,- o a) i D C)0 C" c-) 4-' 4-' ( -,_ C a3 r <D (o) o C) 4 a, A (I) >- m in E 0 1- 0 _ E: --- c- co ciO (:o 0)-~ (D 2 U) (1) >~~Q J) W O0'-... ~~ wE. In 0 I c7cOz

Phase II Tasks - Year 1 (cont'd) As Planned at Kickoff 4. Validation of tetrahedral code using Phase I Brick code. - solver choices and convergence studies - array simulation - biasing studies - non-uniform bias field effects 5. PML implementation with tetrahedral elements for non-planar printed antennas (by June 1998). University of Michigan 7

Summary of Actual Progress to Date * Geometry Tasks Completed (by Zhifang Li) - MeshDriver (patches, log-periodics....) - Prism-to-Tetra Mesher - Mesh Formats, Interface with FEM Solver - Mesh and Field Visualization Using Matlab and Meshview * Solver Tasks (A. Brown and Y. Botros) - Ported and Tested GMRES in Ferrite Brick - Solver Parameter Investigation - Completed Ferrite Antennas Paper University of Michigan 8

Summary of Actual Progress to Date-Cont * Ferrite Tetra Solver (A. Brown and Lars Andersen) -Ferrite Formulation - FE-BI Formulation and integration - Eigenvalue Solver for Non-Symmetric Complex Matrices (two tetra codes have been completed for ferrite modeling) University of Michigan 9

Solver updates GMRES and Eigenvalue Solver University of Michigan 10

LIST OF SOLVERS * L INPACK I NASA angley MultIpat CIorp — CG BCGj, etc * iia Resdua - -QMR, GMRES — DGRES University of Michigan 11

Iterative Solvers * CG Family - CG, CGS, BCG * Minimal Residual Family - QMR (several Versions) - GMRES - DQMRES - Flexible GMRES(FGMRES) - Housholder GMRES University of Michigan * Minimizes r= Ax-bl or related residuals by expanding x along direction vector = Xo + aOl + a2X2 +...... *minimize same residual *maintain orthogonality of solution projection vectors *Robust Solvers *Multiple right hand side *Strong preconditioners *Matrix Compression 12

GMRES* Features *Robust ---quarantees convergence *Smallest error for fixed No. of Iterations *Predictable error history *Adaptable to several preconditioners and variable preconditioners *Can incorporate matrix compression (using QR, SVD) *Can Incorporate FMM/AIM for CPU/Memory Reduction * Multiple right hand sides University of Michigan 13

GMRES and FGMRES-Algorithms Algorithm Select M as PC, Say diagonal PC. Initialize I- Initialize x, r=Ax-b, beta=norm(r). Starting vector v(., ) ) r b e ta; m: no. of restarts. H= zeros (m +, m). M: Preconditioning for -=I:m Do;, Mtrin |ziJ)-=M _*v':,tL Difference betw en the GMRES and FGMRE S w (',i)AZ fr i=:i Do FMM/AIM (, = (w (j) v,:,i)) Construct Hessenb( can be Added w(:,)=w(.j)-()v; atrix Using the 'r Here /+, arm._ w m I A Orthogonal + I =w f // - n ) Spanning/Restart EndDo Vectors Compute y that minimizes norm (betae - H *y) where, e=v(:, ]) *v set x=x+y. S olve Hessenberg if converge, stop, end. Matrix for Solution University of Michigan erg.I, 14

GMRES/FGMRES-Parameters * CPU time in all iterative algorithms is due to matrix vector product [A]{x} * GMRES has two crucial parameters that control its performance (CPU and memory): - Number of restarts (m), - m vectors are used at each iteration step to approximate solution * Larger m increase CPU and memory, but provide robustness - Preconditioning (FGMRES allows for flexible preconditioners) University of Michigan 15

', *= * a- - ---- ---....... L''''''''' ''..- l --- — r -- aI- I '. -'.., ' '..................... '......-;.. 4 --- II Cl -'. " -I-I w 0 --- -' i --- -..... ~ ----I — -'..........-.........-......... — ---- ------- ~ 4 0 ci..o PEA /Z....E-4. IRV....... ~' ~- ~i~ ' i ~c 00 Q.Q CL I! IIl r~ /' ~E~...... i......!.. i> i o '.......:~..... i~...... i...........-.o.......,,,.... ~~ os~sa o ~Ei

'I 6 -6) ('5 0 ~ 0) 6: 6D 6 — OD 6) 6D CD 6= 6D F-2 6l 6 6: 6D 6 6) 6) 6) 6) 6j CZ) Ci) rA W tL w;00 0 C/ C40 0 PENO POMO W;MO 0 C40-;, w rA M P04 0 QPM 0 * - 9 O *0 0 0 * - 0

-Approx. Inverse Preconditioner Preconditioning Matrix [M] is Evaluated by Minimizing the Residual R=[I]-[A][M] [A]: System Matrix *Robust and efficient if applied on a percentage of the entire system *Works even for indefinite systems *Column by Column mimimization is employed to reduce CPU *Special error control has been devised to automate percentage of matrix preconditioning as needed *0 to 30% are typical for matrix preconditioning *No of restarts (m) is 1-2 % of unknowns for large systems University of Michigan 18

Preconditioning in FGMRES *Preconditioning is a Crucial Step for Fast Convergence *Should be Done with Minimal User Intervention *Preconditioning Increases CPU(per Iter.) and Memory Comparison Between FGMRES with Diag. PC and FGMRES with AIPO 4500 ' '.; DPC: Diagonal\. 0 FGMRES with.PC *....FGM.with DPC Preconditioner 4 3500 00 AIPC: Approx. \ p Inverse/ A Preconditioner / I 40 50 No. of restarts "m" University of Michigan 19

--- 1 GMRE m -AIPC Solver o r. I aZ1T3 LNk. C7So0~ EXk-w) ft -10 -20 -30 0 L -40 -50 -60 -70 -80 0 10 20 30 40 50 60 70 80 90 100 The effect of the nuimber of NORMI minimization steps per column: for 1=-:N, N is the FEMf size for k= 1; —K K/s the No. of Norm m"innimiztio tp/oun M(;2=*IT~)/Ifs the identy matrx Mis th~e PC. r=I,;i)- A *M)W.1, Ar= A *r; 'Is thecojgttrnps a 1=Ar"'*r; a2=Ar*Ar-fend,' end; 1% 1 --- --- -- -- -- — dff

p- MMWMl Data for te used matrix: No. of No zero elment =639 Generaitin geometry: Spira Anenna. Mesh tye: Non unfom Mesh elemens Mied elements (ricks an pri~sms) SytmQaur:Bal 'condItioned BI par ofa themarxSie70 z z v

# 2. 0 500 0 500 1000 1500 200 50 2000 2500

Poor Convergence Using BiCG Layer 1. '0li 0..O a) 'c: d=4.445 cm Layerl d = 1.905 cm, r = 1.0. a = 1.0 Lay r-2 d = 0 535 cm. ea = 1 3.9. 4rt==BOOGi. Ho=10(Xe. =I=10le Laer3 d = 1.905 cm, r = 1.0. r= 1.0 University of Michigan 21

Improved Convergence for Biased SubstratesUsing GMRES Lawer I d=41445 cm Layer3 1 11 Layerl d=O 105cm, r =1D., [ = 1.0 Layer d = 0.35 cm, Er = 132, 4nMs=OIOG,,H=l OOe Layer3 d =1305cm, r = 1, = 1.0 1.6 1.8 Fie uierny CHZ4 University of Michigan 22

Requirements for general-purpose eigenvalue solver * Eigenvalue analysis is important for verification of computer codes * For the most general ferrite applications, we need an eigenvalue solver that o Solves the generalized eigenvalue problem o Deals with complex-valued matrices o Allows non-symmetric matrices o Works well for large eigenvalue problems (limited CPU-time growth as number of unknowns increases) o Takes advantage of matrix sparsity (limited memory growth as number of unknowns increases) o Allows non-positive definite matrices o Can handle ill-conditioned matrices o Can handle eigenvalue multiplicity o Can find eigenvalues in any user-specified range o Is robust and accurate * This is not a trivial task * Several methodologies/packages were investigated * ARPACK was finally opted for

Outline of Tetrahedral Code Development * Variational Formulation * FEM Formulation * Validation

Variational Formulation. Radiation F(() = IJJ[-1 (Vx) (VxE) - ko2(r ) ]dV V +Il [jkZ int. E-r1 int VxE]dV JJ [ikoZo i E- r ' '' V fiE. [2k2 x [z x ')] Go(r, r')dS' dS Sa Sa * Scattering F(E) = jJJf[r 1(Vx)~ (Vx)-k2(er~ ) E]dV V 2|JE [2ko x JJ [ xr Go(r, r')dS']dS Sa Sa +2jkoZo [2 x J(r)] tinc(r)dS Sa Go(r, r') = I + k2VV )Go(, r') = Free Space Dyadic Green's Functior * Surface integral containing the dyadic is where the Boundary-Integral(BI) is performed. * In the case of a PML/Artificial Absorber, the surface integral term is omitted * Focus will be on the term: lJJl [- (Vx) (Vxt)- k2(r~ E) ]dV V

FEM Formulation 4 3 1 2 Tetrahedral Element Node Connectivity Table edge i il i2 1 1 2 2 1 3 3 1 4 4 2 3 5 4 2 6 3 4 * Use edge elements of the form: 6 E =,Ei]i * After substitution into the functional, need to solve: Ae. ii +Be. = ii JJ [(r1 Vx7i) ~ (Vx.j) V -k2(,r i). Rj]dV 0 i

FEM Formulation (cont.) * What is the form of the basis function? * 3 different but equivalent forms for the edge basis functions. 1) i = i.(L.i I 11 VL -. L. VL. ) 12 12 1l i= 1,...,6 2) i = fi+ gi x r i= 1,...,6 -) (7- i) - (7- i), i 6Vi X i 7-i ie7-i ^ 6V ' e7-i r(7 17-i 1. 3) i = 9V2 [Ai ( - k)]Ai2 -[A i2 k A;- 2 12L -r)]Ai Al =sbqxbr S = sgn[b1 x b2 - Yq2 + Z q2 b3] q I -x )+ (y q2 q2 q2 -q2) * 3) was used because of it's simple implementation - After substitution 3) yields a simple result for the element matrices. - Makes the programming implementation easier.

FEM Formulation(cont.) * To solve for the element matrices the following formula is used: 21. Vx xi -(iAi2) Vi 9v(Ai1 2 * Final closed form expressions for the element matrices: -— e Aii le le 9Ve 1 f 180Ve m )) b(7- j)( 1) m = 1,2,3,4 fi = 2[(r An) An+(~r t) A] -(r ') At-(r At) An} f2 = 2( * A,) * At + (~r * Ap) * P- (er.) * AP)- (A r * Ap) At f3 = - 2(Er * An) t - (Er. A"P) A XP + (E r Arp) t + (PE ' n) r np 4 ( r An) At-(Er An) A- (~rA At + (Er 'Aq) Ap * Element matrix entries are simple and compact!

Geometry of a Circular Patch Used for Validation * BiCG solver * 12,014 total edges, 2,168 PEC edges, 9,846 unknowns - Er=2 ~r=1 - ~r=1- j2.718 lLr=1 for all 3 layers Circular Patch

Input Impedance(Zin=Rin+jXin) vs. Frequency 100 50 I r I I. I ol t- Ii I -50 100 - c (."o 0a: C c -150 -200 H -250 Rin, isotropic form. - - Xin, isotropic form. + Rin, anisotropic form. ' Xin, isotropic form. ". -300 F -350 - -400 L 0 I I I I I I I 0.5 1 1.5 2 Frequency(GHz) 2.5 3 3.5 4

Eigenvalue computation using ARPACK * Development: o Lehoucq, Sorensen and Yang at Rice * Availability: o Public domain - Netlib * Methodology: o Implicitly restarted Arnoldi method (Krylov method) * Main advantages: o Fulfils all requirements from the previous slide c Does not assume anything about the matrix storage format. Via a reverse communication loop, the user specifies the action of a matrix on a vector. The user can do this using any desired storage format (full, banded, sparse,...) o Fixed pre-determined storage requirements o User-specified numerical accuracy o Numerous sample drivers available for a variety of problems => often requires very little implementation to solve a given problem

Eigenvalue computation using ARPACK * Consider the generalized eigenvalue problem Ax AMx * This is equivalent with the standard eigenvalue problem (A - M)1Mx - x A* Eigenvalues A closest to C of the generalized eigenvalue problem correspond to the numerically largest eigenvalues (A - o)-1 of the standard eigenvalue problem * To solve this eigenvalue problem, ARPACK requires the computation of w = (A - M)-1Mv for many vectors v * This is not trivial since (A - oM)-1 is not explicitly known * However, the computation of w corresponds to solution of (A - (M)w = Mv for w * This can be done effectively using an LU-decomposition scheme that takes advantage of matrix sparsity * Presently, an algorithm based on banded matrices is used

Eigenvalue computation using ARPACK * The flow of the code then becomes the following o Compute A and M o Convert A and M to banded matrices (REDUCE) o Compute A - aM (banded) o LU-decompose A - aM (banded) o Enter reverse communication loop. When ARPACK requires the computation of w (A - oM)-1Mv, simply solve (A - rM)w = Mv for w using substitution. Note that the same LU-decomposition is used every time. O Upon exit of the reverse communication loop, extract eigenvalues * The LU-decomposition is currently based on a LAPACK routine requiring storage of all elements in the band - zeroes as well as non-zeroes. This is (though significantly better than storing full matrices) a complete vaste of memory and should be changed. Note, however, that this storage issue is independent of ARPACKs eigenvalue computation.

Hierarchical tangential vector finite elements * Tangential vector finite elements (TVFEs) are superior to node based finite elements * Nedelec pointed out the attractiveness of mixed-order TVFEs * Hierarchality is of interest for certain applications * Development of hierarchical mixed-order TVFEs for triangular element * Thorough test of these for two-dimensional scattering computations * Conclusion: When hierachality is exploited, we obtain more accurate results using less CPU time and less memory as compared to the case where the same TVFE is used throughout the computational domain * Hope/expectation: Conclusion holds for three-dimensional applications * Development of hierarchical mixed-order TVFEs for tetrahedral element * Thorough test of these for computation of eigenvalues of metallic cavities (partly) filled with isotropic or anisotropic dielectric and/or magnetic materials is in progress * If hierarchality offers advantages for modeling ferrite materials, more advanced ferrite structures will be analyzed

80. 70. Fig. i Cylindnrical anisotropic resonator in i F i. II 3 -4 \1 f 0 0 '= 0 f 0 0 0 c, 0 02 q 06 08 0 12 1 Fi, 2 Normalized resonant frequencies for:-dep dent nmot, f t!, cylindrical cavity completely filled with gyromagnetic medium cr,, K/p alus. 1. - = R = p.. C=.xact Ivale approximate values obtained bv means of the c(-R-R mciWo( IL +-4 11= JK 0 L - jK 0] /t. 0. 0 /J. -I (-J K 0 0.2 0 I K. - I. - j 06 08 t0 p - gio 3. Normalized reo)nant frequencies fu(r:-indepedcr the cylindrical cavi\ts o lv filled ith gyromagn- ', ver-sus |K/| valucs. -. u = - - exact v alues: inatl values obtainedb 3m-t-r-o the C( - R method!I I(.d( II

IEEE TRANSACTIONS ON ANTENNAS AND PROPAGATION, VOL. 45, NO. 11, NOVEMBER 1997 1689 Designing Broad-Band Patch Antennas Using the Sequential Quadratic Programming Method Zhifang Li, Panos Y. Papalambros, and John L. Volakis, Fellow, IEEE Abstract- The utility of numerical codes is greatly enhanced if they can be used in design, a situation that typically involves iterative optimization algorithms. An attractive way is to use gradient-based algorithms developed for solving nonlinear programming (NLP) problems. In this letter, we examine the performance of a general sequential quadratic programming (SQP) optimization algorithm for designing patch antennas in conjunction with a finite-element boundary-integral code. Index Terms- Microstrip antennas. I. INTRODUCTION NTENNA design involves the selection of the physical antenna parameters to achieve optimal gain, pattern performance, VSWR, bandwidth, and so on, subject to specified constraints. Over the past ten years, sophisticated computer codes have been developed for antenna analysis [1]-[3] based on a variety of popular methods. By and large, these codes have not been extended to include design capabilities primarily because of their complexity and nonlinearity with respect to the physical properties of the antenna (material constants, dimensions, feed location, and type, etc.). Some design algorithms have been proposed but these are applicable to specialized antenna shapes and do not address the general antenna optimization problem [4]. Recently, genetic algorithms (GA's) have been examined for array design and absorber optimization [5]-[7]. However, GA's, although robust, require large number of function evaluations to complete the optimization study. Also, GA's are more suitable for discrete variable problems. In contrast, antenna simulations rely on complex computationally intensive codes, which generate continuous functions. It may, therefore, be impractical to generate a sufficiently large sample space for carrying out an optimization study using GA's. An alternative optimization algorithm is the sequential quadratic programming (SQP) method, suitable for continuous nonlinear objective functions such as the input impedance, gain, pattern shape, etc. with both equality and inequality constraints. Convergence is typically achieved in a few iterations and, therefore, their interface with rigorous (but expensive) numerical antenna analysis codes is much more practical. SQP and other similar algorithms are routinely used for large structural design problems involving finite-element Manuscript received April 18, 1997; revised July 10, 1997. Z. Li and J. L. Volakis are with the Department of Electrical Engineering and Computer Science, University of Michigan, Ann Arbor, MI 48109 USA. P Y. Papalambros is with the Department of Mechanical Engineering and Applied Mechanics, University of Michigan, Ann Arbor, MI 48109 USA. Publisher Item Identifier S 0018-9'6X(97)07992-1. analysis [8] and, thus, we can benefit from the extensive experience available in other disciplines. In this letter, we examine the performance of a general SQP code [9] for designing patch antennas in conjunction with a finite-element boundary-integral code [10]. Both are rigorous general-purpose codes. The main point of the paper is to examine the suitability of SQP for antenna parameter optimization to achieve the design objectives subject to constraints. We will illustrate the performance of the optimizer using a few illustrative examples from simple to more complex. II. SEQUENTIAL QUADRATIC PROGRAMMING ALGORITHM DESCRIPTION SQP is a gradient-based class of methods that became prominent in the late 1970's [11]. They are considered the most efficient general-purpose nonlinear programming algorithms today. The basic principle of sequential approximations is to replace the given nolinear problem by a sequence of quadratic subproblems that are easier to solve. Consider the equality constrained problem min f(x) subject to h(x) = 0 (1) where x is the design variable vector, f(x) is the objective function, and h(x) is the vector of equality constraints. Using a Lagrange-Newton method (see, for example, [11]), at the kth iteration, we have Wk Ak 1 - vfk1 LAk 0 [Ak+l L -hk j (2) where W = V2f + ATV2h. A = Vh, and A is the vector of Lagrange multipliers. Solving the above equations iteratively, we obtain the iterates Xk+l = Xk + Sk and Ak+l which should eventually approach x and A, the optimal values. We observe that the above equation can be viewed as the first-order optimality (Karush-Kuhn-Tucker) conditions for the quadratic model min q(sk) = fk + VxLkSk + 5 kWkSk subject to Aksk + hk = 0 (3) where VxLk - Vfk + A TVhk. Solving the quadratic programming subproblem (3) gives the same Sk and Ak+i as solving (2) and thus the two formulations are equivalent. The values of x and A can be obtained from solving a sequence 0018-926X/97$10.00 ~ 1997 IEEE

1690 IEEE TRANSACTIONS ON ANTENNAS AND PROPAGATION, VOL. 45, NO. II, NOVEMBER 1997 ANALYSIS PROCESS DESIGN PROCESS Cavity (14.325 cm x 12.375 cm) Patches (5.73 cm x 6.60 cm) a,. PEC - I di=1.59 mm o Gd2 =? oa dI=1.59 mm Coax Feed Ground Plane Fig. 2. Geometry of the dual-patch antenna. Fig. 1. Flow chart of SQP combined with FEM. 1.42 of quadratic programming (QP) subproblems, hence the name SQP methods for the relevant algorithms. Proper convergence proeprties are achieved with some modifications on this basic SQP algorithm. We may view Sk as a search direction and define the iterate as Zk+l = Xk + akSk where the step size ak is introduced and computed by minimiznig an appropriate merit function along the search direction. The QP subproblem can be solved efficiently by well-developed QP solvers, based, for example, on projection or augmented Lagrangian methods. Using an active set strategy, problems with both equality and inequality constraints can be solved. III. COMBINING SQP WITH THE FINITE-ELEMENT METHOD In the next section, we give two examples of microstrip antenna optimization using the SQP algorithms. For the calculation of the objective function, a hybrid finite-element algorithm is used to compute electromagnetic scattering and radiation by an open three-dimensional rectangular cavity recessed in an infinite ground plane [10]. The cavity may support microstrip patch or slot antennas and may be filled with layered dielectric material. The problem is formulated using the finite-element boundary-integral method and the resultinng system of equations is solved via the biconjugate gradient method. Besides different layers of dielectrics, the cavity may also have lumped loads, probe feeds, and short circuit pins. The flow chart of the whole process is shown below in Fig. 1. Accordingly, the finite-element method (FEM) code first computes the objective function using an intial set of antenna parameters which is used by the optimizer to determine the new search direction and step size. The process is repeated until convergence within the given tolerance is achieved. IV. EXAMPLE APPLICATIONS A. Probe-Fed Dual Patch A number of techniques have been suggested and implemented to improve the bandwidth of the microstrip patch antenna. One of them is stacking patches horizontally or vertically [12]. Several parameters (size of patches, substrate 0 -o a.0 Ce3 CZ, u, m, 0 I0 o o CD r_ fo 0. 0 3" 1 2 3 4 5 6 7 8 Iteration Number Fig. 3. Iteration history for the optimization of the dual-patch antenna. thicknesses, feed locations, etc.) can be used to maintain a low voltage standing wave ratio (VSWR) over a given frequency range. For our purposes, we choose the VSWR to be around two, which corresponds to a return loss of about 3 dB. Thus, we can write the problem statement as Ro/Rin <2 where Ro refers to the input resistance at resonance and Rir is the corresponding resistance at nearby frequencies. Fig. 2 illustrates the configuration of the stacked antenna under investigation. The top and lower substrates have a dielectric constant of E, = 2.2 and thickness h = 1.59 mm and the driven patch is designed to operate at 1.53 GHz. The middle substrate has a relative permittivity of,E = 1.1. We wish to find the optimum length, width, and separation of the patches to achieve a 15% bandwidth. As a starting point in the optimizer, we use the values from the cavity model [13], i.e., L = 5.73 cm and W = 6.60 cm, respectively. The top patch was chosen to have the same dimensions. After nine iterations of the optimizer, a value of d2 = 14.25 mm was determined which delivers a bandwidth of 15%. The VSWR is less than 2.09 within the entire bandwidth. The center frequency is 1.53 GHz, if the patch size is 5.73 cm x 6.6 cm and Fig. 3 gives the iteration history for the dual-patch optimization. B. Slot-Fed Dual Patch There are many ways to improve impedance bandwidth such as impedance matching and multiple resonances [14].

LI et al.: BROAD-BAND PATCH ANTENNAS USING SEQUENTIAL QUADRATIC PROGRAMMING METHOD 1691 Top Patch Substrate; Bottom Patch Substrate Slot Substrate o 3 C4. r 0 -00;o CD 0 f-9 9 (p CD 0111 Ground Plane Fig. 4. Three-dimensional view of the stacked patches with aperture. The latter (borrowed from tuned electronic amplifier design) is a popular approach and introduces additional resonant patches to provide two or more closely spaced resonances. The design principle of the previous example falls into this category. However, that example dealt with only a single-variable optimization. To improve the bandwidth for a dual-patch configuration, we now consider a multivariable optimization. The structure of the stacked-patch antenna is shown in Fig. 4. The size of the top patch is slightly larger than that of the bottom patch to get a better bandwidth and a slot feed is put below the bottom patch. Among all geometric parameters which are sensitive to the design, we choose three for optimization. Specifically, the slot length, and the two substrate thicknesses were chosen for optimization with the goal of achieving a bandwidth more than 15% with VSWR < 2. The dimensions of the other parameters are given below: TOP PATCH: length = 3.89 mm/ width = 5 mm thickness to be decided; rE = 2.33; BOTTOM PATCH: length = 3.5 mm width = 4.5 mm thickness to be decided;,r = 2.2; SLOT FEED: slot-width = 0.5 mm slot length to be decided; GROUND SUBSTRATE: thickness = 0.508 mm; sr = 2.2. The stacked patches are residing in a cavity 1.4 cm x 1.8 cm in size. The problem statement is the same as in the last example, i.e. min Ro/Rin plus some size constraints, and the antenna analysis simulation is also based on the same hybrid finiteelement code. After six iterations, the SQP optimizer found the following optimal values for the three unknown parameters: thickness of the substrate supporting top patch = 0.85 mm; thickness of the substrate supporting the bottom patch = 0.55 mm; and slot length = 4 mm. The performance of this antenna is shown in Fig. 5, with a VSWR equal to only 1.414 over an 18% bandwidth. V. CONCLUSION The advantages of the SQP algorithm are fast convergence and reduced number of function evaluations. This is attractive 17 17.5 18 18.5 19 Frequency (GHz) Fig. 5. Performance of the optimized dual patch-slot antenna. when expensive computations are needed within the optimization loop. SQP achieves its speed by restricting its search in a more narrow range (local optimization) than genetic algorithms and by making use of gradient information. This is suitable for a large class of antennas where a limited range of parameter values is sufficient for good performance. REFERENCES [1] D. M. Pozar and D. H. Schaubert, Microstrip Antennas. New York: IEEE Press, 1995 (reprints volume). [2] J. L. Volakis, J. Gong, and T. Ozdemir, "Application of the FEM to conformed antennas," in Finite Element Software for Microwave Engineering, T. Itoh, G. Pelosi, and P. Silvester, Ed. New York: Wiley, 1996, pp. 313-345. [3] C. Penny and R. J. Luebbers, "Radiation and scatteirng of a square archimedean spiral antenna using FDTD," Electromagn., vol. 14, no. 1, pp. 87-98, 1994. [4] R. L. Haupt, "Introduction to genetic algorithms for electromagnetics," IEEE Antennas Propagat. Mag., vol. 37, pp. 7-15, Apr. 1995. [5] R. L. Haupt, "Thinned arrays using genetic algorithms," IEEE Trans. Antennas Propagat., vol. 42, pp. 993-999, July 1994. [6] D. S. Weile, E. Michielssen, and D. E. Goldberg, "Genetic algorithm design of Pareto optimal broadband microwave absorbers," IEEE Trans. Electromagn. Compat., vol. 38, pp. 518-525, Aug. 1996. [7] E. Michielssen, M. J. Sajer, S. Ranjithant, and R. Mittra, "Design of lightweight, broad-band microwavew absorbers using genetic algorithms," IEEE Trans. Microwave Theory Tech., vol. 41, pp. 1024-1030, June/July 1993. [8] N. Olhoff and G. I. N. Rozvany, Eds., Structural and Multidisciplinary Optimization. Oxford, U.K.: Pergamon-Elsevier, 1995. [9] R. L. Crane, K. E. Hillstrom, and M. Minkoff, "Solution of the general nonlinear programming problems with subroutine VMCON," Argonne National Labs., 1990. [10] J. L. Volakis, J. Gong, and A. Alexanian, "A finite elemen t-boundary integral method for antenna RCS analysis," Electromagn., vol. 14, no. 1, pp. 63-85, 1994. [ 11] P. Y. Papalambros and D. J. Wilde, Principles of Optimal Design. New York: Cambridge Univ. Press, 1988. [12] J. S. Dahele, S. H. Tung, and K. F. Lee, "Normal and inverted configurations of the broadband electromagnetically coupled microstrip antennas," in IEEE Antennas Propagat. Soc. Int. Symp. Dig., Philadelphia, PA, June 1986, pp. 841-844. [13] J. R. James and P. S. Hall, "Handbook of microstrip antenna," Inst. Elect. Eng. Electromagn. Waves. London, U.K.: Peter Peregrunus, 1989, ser. 28. [14] D. M. Pozar, "A review of bandwidth enhancement techniques for microstrip antennas," Microstrip Antennas: The Analysis and Design of Microstrip Antennas, D. M. Pozar and D. H. Schaubert, Eds. New York: IEEE Press, 1995, pp. 157-166. [151 F. Croq and D. M. Pozar, "Millimeter-wave design of wide-band aperture-coupled stacked microstrip antennas," IEEE 7rans. Antennas Propagat., vol. 39, pp. 1770-1776, Dec. 1991.

1692 IEEE TRANSACTIONS ON ANTENNAS AND PROPAGATION, VOL. 45, NO. 11, NOVEMBER 1997 Zhifang Li was born in China in June 1972. She received the B.S. degree in mathematics from Peking University, China, in 1995, and the M.S. degree in electrical engineering from the University of Michigan, Ann Arbor, in 1996. She is currently working toward the Ph.D. degree in electrical engineering and computer science at the University of Michigan. Her research interests include antenna optimization and numerical methods. John L. Volakis (S'77-A'79-S'80-M'82-SM'88-F'96), for photograph and biography, see p. 507 of the March 1997 issue of this TRANSACTIONS. Panos Y. Papalambros received the Diploma degree in mechanical and electrical engineering i from the National Technical University of Athens, Greece, in 1974, and the M.S. and Ph.D. degrees from the Design Division of Mechanical Engineering, Stanford University, Stanford, CA, in 1976 and 1979, respectively. He joined the University of Michigan faculty in 1979 where he currently serves as the Chair of the Department. He has published over 100 articles in the field of design automation and has co-authored (with D. J. Wilde) the textbook, Principles of Optimal Design: Modeling and Computation. He serves on the editorial boards for the following journals: Artificial Intelligence in Engineering Design and Manufacturing, Engineering Design, Engineering Optimization, Integrated Computer-Aided Engineering, and Structural Optimization. He is also past editor of the Japan Society of Mechanical Engineers International Journal, Global Optimization and the ASME Journal of Mechanical Design. Dr. Papalambros is the Founding Director of the multi-university Automotive Research Center in Modeling and Simulation of Ground Vehicles, funded by the United States Army. He is a member of ASME, INFORS, MPS, SIAM, SME, SAE, and ASEE and a Fellow of ASME.

Patch Antennas on Ferromagnetic Substrates Arik D. Brown, John Volakis, Leo C. Kempelt, and Youssry Botros Radiation Laboratory Department of Electrical Engineering and Computer Science University of Michigan Ann Arbor, MI 48109-2212 tMission Research Corporation 147 John Sims Parkway Valparaiso, FL 32580 Abstract Patch antennas on ferrite substrates allow for pattern control, frequency shifting, and scattering reduction. This is achieved by external magnetic field biasing coupled with the inherent magnetization of the ferrite substrate. Measurements and analytical studies based on the Moment Method (MoM) have verified these attractive properties of ferrite substrates. However, verification of the analysis is difficult, and furthermore previous models have relied on uniform biasing across the substrate. In this paper, we present a hybrid finite element-boundary integral (FE-BI) method which permits modeling of the true non-uniform bias fields within the substrate for a more accurate prediction of the ferrite patch performance. After validation of the proposed simulation, and a demonstration of the inherent properties of the ferrite patch, it is shown that non-uniform biasing is responsible for additional frequency shifts. We also identify the poor condition of the resulting matrix systems and relate this situation to the predictable occurrence of non-propagating substrate modes. A more robust iterative solver with preconditioning is therefore proposed and applied to handle these situations. 1 Introduction Patch antennas on ferrite substrates are attractive because they offer greater agility in controlling the radiation characteristics of the antenna. Their inherent anisotropy and non-reciprocal properties [1], permit variable frequency tuning [2], [3], [4], and antenna polarization diversity [5]. External biasing of the ferrite substrate also allows for beam steering [6], [7], [8], [9], pattern shape control, and radar cross section control [10], [11] by forcing the ferrite into a cut-off state [12]. Several papers have already considered the performance of ferrite patch antennas. These works [13], [14] employed the Moment Method (MoM) technique in conjunction with the sub

strate Green's function. Validation of the results given in [13], [14] have so far been difficult to achieve. Also, MoM formulations do not permit modeling of non-uniform biasing and inhomogeneous constitutive parameters, a situation which inherently occurs when the ferrite is biased. To provide for greater flexibility in modeling the ferrite substrate and the substrate cavity (see Figure 1), we performed an analysis of the ferrite patch using the Finite Element-Boundary Integral (FE-BI) method. As usual, the substrate housed within the cavity is modeled by the Finite Element Method (FEM) using an edge-based formulation [15], [16]. Consequently, multiple substrate (and superstrate) layers can be handled easily including lateral material inhomogeneities within each layer. In our formulation, the FE mesh is truncated at the surface of the cavity using the rigorous BI method. Thus, the proposed FE-BI implementation is equally rigorous to the traditional MoM (employing the substrate Green's function) and allows modeling of finite and inhomogeneous substrates. In the following sections, we first give a brief description of the formulation and it's implementation for ferrite materials. We then proceed with the presentation of some simple ferrite patch antenna calculations involving scattering and radiation examples. These serve to validate the implementation and preciseness of the results. They also reveal that serious solution convergence difficulties can arise for certain bias states of the substrate. It is shown that these difficulties can be predicted a priori. For our implementation we resorted to a more robust iterative solver, the Generalized Minimal Residual Method(GMRES) with preconditioning. Calculations showing the effects of biasing on the antenna resonance and scattering characteristics are given in the latter part of the paper. Our final example is a calculation simulating a non-uniform substrate (due to natural biasing). This example demonstrates the importance of modeling these non-uniformities accurately for evaluating the performance of the radiating element. 2 Formulation The geometry to be considered, (see Figure 1), consists of a patch situated in a rectangular cavity. The presence of the cavity eliminates radiation loss via surface waves, and does not affect the radiation pattern provided the cavity perimeter is placed at some small distance from the patch edges. To obtain the unknown field, in the context of FEM, the variational equation F(E)=0 (1) is solved [16], where F(E) = J [(V X E) (V x E) - kcrE.E]dV 2 r 1 + I IjkZJint. E- 1 Mint. (V x E)]dV +jkoZo J (E x H).zdS (2) 2

z x x- y E > H Patch Antenna Ferrite Substrate Figure 1: Geometry for a patch antenna on an anisotropic substrate. In this equation, V denotes the cavity volume, S is the cavity aperture, cr and /ur are the relative permittivity and permeability of the ferrite substrate, Jint and Mint are internal electric and magnetic sources due to the antenna feeds, and the last term in (2) is the BI term. Discretization of (1) using Galerkin's method leads to the linear system, [A]{E} = {B} (3) where [A] is an NxN matrix and {B} is an Nxl column vector given by [15]. When modeling gyromagnetic substrates, the functional must be modified to incorporate the inherent anisotropy of the ferrite material. Specifically, for general anisotropic media, we have (E) = Jv[l x E) (V x E) - k M- E. E]dV + J [jkZJint E - Mint (V x E)]dV v r +jkoZo ( S(E x H) zdS (4) where e and Lr are the relative permittivity and Polder permeability tensors. The element matrices in the FE assembly process, resulting from this functional, are given in Appendix A. For a z-biased ferrite, uZr is given by 3

I jK 0 - -jK m o. (5) 0 0 O ] and for other biasing directions(x and y), the tensor entries are simply rotated accordingly [1]. Here the parameters, [ and K, are functions of frequency given by (1 + ) (6) 0 - Jo( WWm ) (7) where W = 7(Pou0) (8) and Wm = 7(yoM). (9) Also, Ms is the saturation magnetization, Ho is the DC bias field, 7 is the gyromagnetic ratio, and wo and,m are the precession and forced precession frequencies, respectively. When dealing with ferrite materials, the field behavior is determined by the propagation direction and its orientation with the applied magnetic bias field direction. There are two separate cases which determine the effective permeability([teff) within the ferrite [1], [19] - the longitudinal case where propagation is parallel to the applied bias field and the transverse case where propagation is perpendicular to the applied bias field. In the longitudinal case /Pff = 1 ~ (10) whereas in the transverse case H2 2 I1eff = (11) For both propagation modes, the propagation constant within the ferrite is calculated as - = jVEoErtoPeff (12) = a +j3 (13) The modes due to the propagation constant play a major role in the FE solution. Because of their tensor properties, ferrites introduce a great deal of complexity into the formulation when solving radiation problems. When using the FE method, it is observed that the system matrix becomes asymmetric and can be poorly conditioned at certain values of the ferrite 4

E 0 250 - x 200 - E 150 -C 100 - II I 50 II -50l 0 1 2 3 4 5 6 Frequency (GHz) Figure 2: Real Input Impedance parameters. Initially, the Bi-Conjugate gradient(BiCG) method was used for solving the matrix system. To improve performance, a preconditioned BiCG algorithm was also examined. However, the BiCG was not robust under certain bias conditions. In these cases we resorted to the GMRES method as described later. 3 Applications and Validation 3.1 Probe-Fed Patch Antenna Consider the probe fed patch antenna geometry given in Figure 1. For this example the ferrite substrate parameters were 47rMs = 650 G, Ho = 600 Oe, and cr = 10. The calculated input impedance and radiation pattern are given in Figures 2 and 3. As expected, biasing caused a shift in resonance and this is clearly seen in Figure 2. Specifically, the ferrite substrate decreased the lowest resonance of the patch from 4.44 GHz to 2.24 GHz, thus reducing the overall size of the patch for operation at the same frequency. From Figure 3, we also observe that the biased patch exhibits a null along the horizontal direction. This patch was also considered by Schuster and Luebbers [17] 5

Figure 3: Radiation pattern in the yz plane on an x-biased ferrite substrate(Frequency = 2.2 GHz, - Isotropic Substrate, - - Anisotropic Substrate) using the Finite Difference Time Domain (FDTD) method. Our computed resonance shift was within 40 MHz of their values (1.8%). Although this type of agreement is considered very good for patch antennas, the small difference may be attributed to possible numerical implementation inaccuracies. Our simulation used a cavity size of 4.085cm x 4.085cm x.015cm and the FE-BI system consisted of 3766 unknowns. 3.2 Ferrite Filled Cavity 3.2.1 Biased Substrate Consider now a cavity with several magnetized layers as shown in Figure 4. Layers 2 and 4 are magnetized in the y direction, e.g. K 0 jK \ r = 0 o 0 (14) — jKi 0 u This is a particular example considered by Kokotoff [18]. The RCS of the layered ferrite cavity for different biasing values(Ho) is given in Figure 5, and our calculations are seen to be in agreement to those of Kokotoff [18] for all cases. This example again demonstrates the frequency shifting property of ferrite materials with biasing and validates the employed FEM formulation. The number of unknowns for this example was 6,776(BI unknowns = 420). 3.2.2 Unbiased Substrate We next consider a 3 layer cavity consisting of a ferrite layer between two free space layers. The ferrite layer is magnetized with parameters Er = 13.9 and 47rM, = 800 G. However, no biasing is 6

z x VZ-Do d=5.158 cm Layerl Layer2 Layer3 Layer4 Layer5 d =0.726 cm, ~r =2.2, trl=.O d = 1.790 Cm, er = 13.9, 4itMs=800G, AH=50e d =0.737 cm, Er =2.2, gir =1.0 d = 0.762 Cm, Er = 13.9, 4nrMs=800G, AH=5 Oe d =1.143 Cm, Er =1.0, grl=.0 Figure 4: Geometry of a cavity with ferrite layers. 0 -5 Ho=4450e --- Ho=5OO0e -10.- - Ho=6000e..... Ho=7000e + Kokotoff Data 15' 20 030/ 0.85 Frequency(GHz) Figure 5: Effect of biasing on the RCS of the cavity in Figure 4. All computations were carried out using the FE-BI method except as noted. 7

z x * —^ y ] __ Layer 1 d=4.445 cm Layer 3 Layeri d = 1.905 cm, Er = 1.0, gir = 1.0 Layer2 d = 0.635 cm, Er = 13.9, 47Ms=800G, AH=100e Layer3 d = 1.905 cm, Er = 1.0, ur = 1.0 Figure 6: Ferrite Cavity Geometry applied. As shown in Figure 7, calculations using the BiCG solver for Ms = 0 were in complete agreement with results given by Kokotoff. Figure 8 shows results for the same geometry with the ferrite layer magnetized. About 780 unknowns (BI unknowns c- 180) were used to simulate this cavity. This example presented us with convergence difficulties when the BiCG solver was used. An investigation of several other cases demonstrated that, in general, convergence difficulties were encountered when the propagation constant /, as given in (13), was zero for one of the modes. For this example, /3 vanished for one of the longitudinal modes corresponding to Pjeff = pt + K and the transverse modes. The actual values of/3 for all three modes are given in Figure 9 and we observe that for the aforementioned two modes, /3 vanishes from about 1 GHz to 2.5 GHz (see Figure 9. In concert, the BiCG solver failed to converge within this frequency range. Results based on a direct solver were also inaccurate due to the poor system condition. To overcome convergence difficulties for those frequencies where /3 = 0 for one or more of the modes, we resorted to a more robust iterative solver such as the preconditioned flexible GMRES (FGMRES) [20]. Features that made the FGMRES algorithm attractive were its guaranteed convergence, ability to adapt variable preconditioners, and a predictable error history (i.e. a smooth and monotonic convergence pattern as compared to the erratic convergence pattern of the BiCG algorithm). An important parameter for the GMRES solver is the number of interior iterations(m) before restarting the solver. These initial iterations control the number of spanning basis vectors used for an initial approximation of the solution. For our examples, the minimum m used was 70 while the maximum 8

I!.I...... -10 0 ca -20 ~ -30 cr -5 -50 I +/ I I/ $' + +, I I I + I / I + Kokotoff Data I ' - - BiCG Validation II + -60 ' I I I I' 0.8 1 1.2 1.4 1.6 1.8 2 Frequency (GHz) 2.2 2.4 2.6 2.8 Figure 7: RCS(normal incidence), using the FE-BI, of the loaded cavity in (i.e. no magnetization). Figure 6 with Ms = 0 Figure 8: RCS(normal incidence), using the FE-BI, of the loaded cavity shown in Figure 6 with the values of Ms and AH as given there. 9

Frequency(GHz) 2 2 Figure 9: P Normalized to k, - IPeff = - - leff -= + /, -. [Ueff = - - c for a ferrite medium having 47rM, = 800G, AH = 10 Oe, and c = 13.9. m was 280. For frequencies where the system is ill-conditioned a higher value for m is required along with preconditioning. From our analysis, these points occur near resonance, which is approximately 1.98 GHz (Figure 8). Using the GMRES solver, with the Approximate Inverse Preconditioner(AIPC) [20], convergence was obtained at all points for the geometry in Figure 6. Figure 8 shows the results, and we observed that the GMRES solution tracks the data in [18] quite well. Given the poor condition of the system, it is not clear as to which of the curves in Figure 8 is not accurate. 4 Nonuniform Biasing When building a ferrite antenna a permanent magnet is required to produce the applied magnetic bias field. Due to the finite nature of the magnet, the field is no longer uniform and thus the electrical material properties become inhomogeneous. Since many analysis methods assume a uniform bias field, this produces a solution which is no longer accurate. In contrast, the FEM allows for arbitrary specification of the material within the volume which is an inherent advantage of FEM over other numerical methods. To observe the effect of non-uniform magnetization, let us consider the modeling of the measured bias field in a ferrite cavity as given in Figure 10 [18]. Indeed, Figure 10 reveals a difference of more than 1000G at different locations within the cavity, showing the necessity of the FEM technique to handle this inhomogeneous behavior. RCS calculations for this non-uniform biasing are provided in Figure 11 for a 6x6x1 cm cavity filled with this material. It is clear, that the resonance of the cavity is substantially affected by the non-uniformity of the bias field. 10

5 Conclusion In this paper we presented several results and validations demonstrating the attractive properties of ferrite patch antennas. The high dielectric constant of the ferrite, inherent magnetization, and external biasing all serve to minimize the size of the patch, in addition to providing pattern control and lower radar cross section over a given band. The employed hybrid FE-BI method also permitted an investigation on the effects of the typical non-uniform bias fields which occur across the substrate volume. These non-uniform bias fields cause inhomogeneities which affect the operation frequency and overall response of the antenna and may be a cause of discrepancies between measurements and calculations. Our study also showed that poor matrix conditions and solution convergence (lifficulties may be traced to band regions where one or more ferrite modes are non-propagating. This situation prompted the use of more robust iterative solvers, and, to achieve convergence, a preconditioned version of the GMRES method was used. GMRES proved effective in cases where the usual conjugate and biconjugate gradient algorithms failed. 11

250(:, ----.-15. - ~ ~ I ".I 1..' UV:I:'l........~ ~~ -........ I... I...I..... I....: 2 magnets with 4 steel spac.s 2 rneWIS wh 4 Al spac.s 2 metws 2 2 Al s..ewI pa.2 m gne........... 6 magnets - -.................... - ^> **. o.I.I o - - -.........,............................................................................... 01 -0:75 -.05. 0. 25.25 O.5.,75 1 Measurement of the nonuniform magnetic field within a cavity. Figure 10: FE co CO1 z~ 0.0 -10.0 -20.0 -30.0 -40.0 -50.0 -60.0 0.0 1.0 2.0 3.0 Frequency [GHz] Figure 11: RCS due to a nonuniform magnetic field (see Figure 10) across a 6x6x1 cm cavity. 12

References [1] D. M. Pozar, Microwave Engineering, Addison Wesley Publishing Company, Reading, MA, 1996. [2] P.J. Rainville and F. J. Harackiewiz,"Magnetic tuning of a microstrip patch antenna fabricated on a ferrite film," IEEE Microwave Guided Wave Letters, vol. 2, p. 483-485, 1992. [3] H. How and C. Vittoria,"Radiation frequencies of ferrite patch antennas," IEEE Electronics Letters, vol. 28, no. 15, pp. 1405-1406, 1992. [4] D.M. Pozar and V. Sanchez,"Magnetic tuning of a microstrip antenna on a ferrite substrate," Electronic Letters,vol. 24 no.12, pp.729-731, 1988. [5] J.S. Roy, P. Vaudon, A. Reineix, F. Jecko, and B. Jecko,"Circularly polarized far fields of an axially magnetized circular ferrite microstrip antenna,"Microwave Opt. Tech. Lett., vol.5, pp. 228-230, 1992. [6] N. Okamoto and S. Ikeda,"An experimental study of electronic scanning by an antenna loaded with a circular array of ferrite rods,"IEEE Trans. Antennas Propagat., vol. AP-27, no. 12, pp. 426-430, 1979. [7] H. How, T. Fang, D. Guan, and C. Vittoria,"Magnetic steerable ferrite patch antenna array," IEEE Trans. on Magnetics, vol. 30, no. 6, pp. 4551-4553, 1994. [8] N. Buris, T. B. Funk, and R. S. Silverstein,"Dipole arrays printed on ferrite substrates," IEEE Trans. Antennas Propagat., vol.41, no. 2, pp. 165-175, 1993. [9] H. Maheri, M. Tsutsumi, and N. Kumagai,"Experimental studies of magntically scannable leaky-wave antennas having a corrugated ferrite slab/dielectric layer structure,"IEEE Trans. Antennas Propagat., vol. 36, no. 7, pp.911-917, 1988. [10] A. Henderson and J. R. James,"Magnetized microstrip antenna with pattern control,"Electronics Letters, vol. 24, no. 1, pp. 45-47, 1988. [11] D.M. Pozar,"Radar cross-section of microstrip antenna on normally biased ferrite substrate,"Electronics Letters, vol. 25, no. 16, pp. 1079-1080, 1989. [12] D.M. Pozar,"RCS reduction for a microstrip antenna using a normally biased ferrrite substrate," IEEE Microwave Guided Wave Letter, vol. 2, no. 5, pp. 196-198, 1992. [13] H.Y. Yang, J.A. Castaneda, and N.G. Alexopoulos,"The RCS of a microstrip patch on an arbitrarily biased ferrite substrate,"IEEE Trans. Antennas Propagat., vol. AP-41, no. 12, pp. 1610-1614, 1993. 13

[14] D. M. Pozar,"Radiation and scattering characteristics of microstrip antennas on normally biased ferrite substrates,"IEEE Trans. Antennas Propagat., vol. AP-40, no. 9, pp. 1084-1092, 1992. [15] J. Jin and J.L. Volakis," A hybrid finite element method for scattering and radiation by microstrip patch antennas and arrays residing in a cavity," IEEE Trans. Antennas Propagat., vol. 39, no. 11, pp.1598-1604, 1991. [16] J.L. Volakis, T. Ozdemir, and J. Gong,"Hybrid finite element methodologies for antennas and scattering,"IEEE Trans. Antennas Propagat., vol. 45, pp. 493-507, 1994. [17] J. Schuster and R. Luebbers, "FDTD for three-dimensional propagation in a magnetized ferrite," IEEE Trans. Antennas Propagat., pp. 1648-1651, 1996. [18] D. M. Kokotoff,Full wave znalysis of a ferrite-tuned cavity-backed slot antenna. PhD thesis, Arizona State University, 1995. [19] H. How, T. Fang, and C. Vittoria,"Intrinsic modes of radiation in ferrite patch antennas," IEEE Trans. on Magnetics, vol. 42, no. 6, pp. 988- 1994. [20] Y. Saad, Iterative Methods for Sparse Linear Systems, PWS Publishing Company, Boston, 1996. 14

A Anisotropic Formulation In the FEM formulation, the relevant integrals to be computed in the volume domain are Ej = j / /V v x Ni r ' v x Nj)dVe Fj = JJJ Ni. (r Nj)dVe (15) (16) where =-I -xx Pr = Iyx izx E(xx ( (zx -xy z \ Lyy Cyz Etzy Ezz ) Exy ~xz %Yy (yz I czy czz (17) (18) (19) Exx Exy Exz Ee = Eyx Eyy Eyz Ezx Ezy Ezz and Fxx XX Fe - F Fr - Fyx Fzx Fxy Fxz Fyy Fyz. Fzy Fzz (20) The values for the brick element matrices, in a general anisotropic medium are 2 -2 1 -1 -2 2 -1 1 i -1 2 -2 -1 1 -2 2 (21) 15

K2 z- i)(22) K3( iIj{ (23) K4~( H (24) K5E2 A1 I 1) (25) K6~Hj i (26) - lxzz K1 + lxlYIIYYK ix/Izy ixyzKi (27) 61Y 6lz 4 4YK 6lz 61X 4 4 3(8 E _ Izy K lxlzIIxxK2lZIIYK~ lIXKT (29) 6lx61 4 4 3 -E-l ftlx K5 + K4 + K6 + lYfYx KT (30) 6 4 4 41Z~ Ex 'Y K xzfz 3+ K4 + K6 (31) 16

1-1, 0 - c 4 N N z::- -,t N 4 -1) N z:- "t,-Q CY') 4;zl) N z::- N N — 4. LO 4 N H ::- C.0 N I I I I N 1) w 1I 11. Q r --- Cq CIIA -!t CIA r --- "Itt CIA cq It r — Cl-I llt cq CIA r — H I.-L) N ('0 m;:Y,-Q H 11. 1lz P;4 C-q -1 E-l -- r-4 r-q E-~r -4 -4 N4 ~N ~ CIIIA CIA II cII cq CN

Hierarchical Tangential Vector Finite Elements for Tetrahedra Lars S. Andersen and John L. Volakis Department Radiation Laboratory, 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. I

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

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. Cn = V,/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 (iVQ-(jV(i, z<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 VC;J, 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 k(Cv(j- Cik), i< < k. (3) 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 -i - <j. (4) ((i- Cj)(~QVCj - (jVC,) } ' 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 (i(2i - 1)V(j, (5) ((V(i- V(j), < j, (6) the 24 face-based vector basis functions (k(2k - 1)(iV(j - (jV(i) 1( <i <j < k, (7) (j(2(j 1-)(kVc - cV(k) J V((ijk), i < j < k, (8) 5(Ck(VCj - CjVC), i j # k # i, (9) and the 3 cell-based vector basis functions C2C3(C1V(4 - 4V(l) 2C4((1 V(3 -3V(1 ) (10) (3C4(1VC(2 - (2V(1)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 2 A correction of the QT/CuN vector basis functions initially proposed by Savage ad Peterson [7] was given by Peterson [9]. This corrected set is the one presented here. 4

basis functions (5)-(6) by (.V( - 1V ((i - (j)((ivj - jv,), i < j. (11) (( - j)CvC - (jV,) Further, we propose to replace the 8 face-based vector basis functions (7) by 'k(C,vc,- jvcn) }, i < j < k.( < i<j k. (12) j ((k VC - ~( v) J 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.

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.C. Nedelec, 'A new family of mixed finite elements in R3', 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

Figure 1: 8