UM Report #025921-6-T TECHNICAL REPORT FOR NASA Grant NAG-2-541 NASA Technical Monitor: Alex Woo Grant Title: A NEW TECHNIQUE FOR SIMULATING COMPOSITE MATERIAL Task 1: A Finite Element Conjugate Gradient FF7l Method for Scattering. Institution: The Radiation Laboratory Department of Electrical Engineering and Computer Science The University of Michigan Ann Arbor, MI 48109-2122 Period Covered: February 1989 - September 1989 Report Tide: A Combined Finite Element and Boundary Integral Formulation for Solution via CGFFI of 2-Dimensional Scattering Problems Report Authors: Jeffery D. Collins and John L. Volakis Principal Investigator: John L. Volakis Telephone: (313)764-0500

en 7 gNI t-LRl-' 6,~" We (7 2 _. )q.fi

Abstract A new technique is presented for computing the scattering by two-dimensional structures of arbitrary composition. The proposed solution approach combines the usual finite element method with the boundary integral equation to formulate a discrete system. This subsequently solved via the conjugate gradient (CG) algorithm. A particular characteristic of the method is the use of rectangular boundaries to enclose the scatterer. Several of the resulting boundary integrals are therefore convolutions and may be evaluated via the fast Fourier transform (FFT)in the implementation of the CG algorithm. The solution approach offers the principle advantage of having O(N) memory demand and employs a one-dimensional FFT versus a two-dimensional FFT as required with a traditional implementation of the CGFFT algorithm The speed of the proposed solution method is compared with that of the traditional CGFFT algorithm, and results for rectangular bodies are given and shown to be in excellent agreement with the moment method.

Contents 1 Introduction 4 2 Analysis 7 2.1 Discretization of the Object and Field Quantities............... 9 2.2 Derivation of the Finite Element Matrix................... 11 2.3 Evaluation of Boundary Integral....................... 15 2.4 A CGFFT Algorithm............................. 25 3 Computational Considerations 29 3.1 Storage Efficiency................................ 29 3.2 Computational Efficiency........................... 35 4 Far Field Computation 37 5 Code Validation 40 6 Conclusions and Future Work 50 7 Program Manual 52 7.1 Description of FECGFFT..................... 52 7.2 Mesh Generator for Curved Bodies........................ 55 7.3 Mesh Generator for Rectangular Bodies................... 66 O

List of Figures 2.1 Geometry of the scatterer........................... 8 2.2 Partially discretized body........................... 9 3.1 Example of the mesh of a penetrable structure............... 33 3.2 Example of the mesh of an impenetrable structure............. 33 5.1 E, backscatter from a.25 x 2A body...................... 41 5.2 Hz backscatter from a.25 x 2A body.. 42 5.3 E, backscatter from a.25x 1A perfect conductor with a A/20 thick material coating containing the properties E, = 5. - j.5,4r, = 1........... 43 5.4 H,, backscatter from a.25 x 1A perfect conductor with a A/20 thick material coating containing the properties E, = 5. - j.5,,r = 1.. 44 5.5 E, backscatter from a.25x 1A perfect conductor with a A/20 thick material coating containing the properties Er = 5. - j.5, r, = 1.5 - j.5... 45 5.6 H, backscatter from a.25 x 1A perfect conductor with a A/20 thick material coating containing the properties Er = 5. - j.5,p,r = 1.5 - j.5....46 1

5.7 Ez backscatter from a.25 x 1, perfect conductor with two A/20 thick top material coatings. The lower layer has the properties E7 = 2. - j.5, /, = 1.5 - j.2, and the upper layer has properties E7 = 4. - j.5, tir = 1.5 - j.2.. 47 5.8 Hz backscatter from a.25 x 1A perfect conductor with two A/20 thick top material coatings. The lower layer has the properties Er = 2. - j.5, /r 1.5 - j.2, and the upper layer has properties Ec = 4. - j.5, Ar = 1.5 - j.2.. 48 5.9 HI scattering from a composite body. Both the perfect conductor and the dielectric body are A/2 in each dimension. The material properties are E, = 5.- j.5,1i, = 1.5 - j.5.......................... 49

List of Tables 2.1 Definitions for the finite element mesh................... 10 2.2 Definitions of the field vectors......................... 11 3.1 List of major memory-consuming variables.............. 31 3.2 List of major memory-consuming variables (continued)........... 32 3.3 Summary of major memory consumption.................. 34 3.4 Summary of mnajor memory consumption: filled mesh............ 34 3.5 Summary of major memory consumption: open mesh............ 35 3.6 A comparision of computation efficiency of the FE-CGFFT with the CGFFT method.................................... 36 3

Chapter 1 Introduction Two-dimensional problems have been studied extensively from both analytical and numerical standpoints for many years. The demand for three-dimensional (3-D) methods has increased in recent years, and as a result two-dimensional (2-D) methods are finding their niche as testing grounds for 3-D applications. The step required to generalize a 2-D method to 3-D is almost always large in analytical and geometrical complexity. Most importantly, though, the demands in computation time and storage are often prohibitive for electrically large 3-D bodies. Vector and concurrent (i.e., hypercube, connection, etc) computers are beginning to eleviate the first of these demands ([1]-[7] to mention a few of the papers addressing this), but storage demands remain problematic. A reduction in storage requirements is therefore essential. The traditional Conjugate Gradient Fast Fourier Transform (CGFFT) method [8], [10] is one such frequency domain solution approach which requires O(N) storage. This method involves the use of FFT whose dimension equals that of the structure under consideration and therefore demands excessive computation time when used in an iterative algorithm. With this in mind, a new solution approach is proposed for solving scattering 4

problems that address this difficulty. The proposed method will be referred to as the Finite Element-Conjugate Gradient Fast Fourier Transform (FE-CGFFT) method and was inspired by Peters [10]. The FE-CGFFT method requires that the scatterer be surrounded by a double rectangular box. Inside the box boundaries, the Helmholtz equation is solved via the finite element method. The boundary constraint is satisfied by an appropriate integral equation, which implicitly satisfies the radiation condition. Along the parallel sides of the box, this integral becomes a convolution and is, therefore, amenable to evaluation via the FFT. The dimension of the required FFT in this method is one less than the dimensionality of the stucture leading to an O(N) memory demand making it attractive for 3-D simulations. The proposed method. is similar to the moment-method version developed by Jin [11]. Jin's method was in turn based upon work published in the early 70's by Silvester and Hsieh [12] and McDonald and Wexler [13], who introduced an approach to solve unbounded field problems. The proposed method is also similar to other methods, a few of which will be mentioned here. The unimoment method [14] uses- finite. elements inside a ficticious circular boundary and an-eigenfunction expansion to represent the fields in the external region. The coefficients of the expansion are then determined by enforcing continuity at the finite element (FE) mesh boundary. The coupled finiteboundary element method [15] uses the finite element method within the boundary and the boundary element method to provide the additional constraint at the termination of the mesh. Unlike the proposed method, the solution was employed by direct matrix inversion as in [11], and the outer mesh boundary was not rectangular to take advantage 5

of the FFT for the evaluation of the resulting convolution integrals. Further, only one boundary was employed, and therefore an analytical evaluation of the Green's function was required. In this work, we consider only rectangular structures and results derived from this formulation are compared to those based on traditional moment method techniques. Nevertheless, the proposed method is equally applicable to more complex geometries by using available sophisticated finite element modelling packages.

Chapter 2 Analysis Consider a cylindrical body of arbitrary cross-section and composition illuminated by the plane wave o e(Ta)= hilC(p) = sierko Pcom( o) (2.1) as indicated in fig. 2.1 (A time dependence of eiwt has been assumed and suppressed.). To evaluate the fields scattered from this object, two boundaries are placed tightly around the body as shown in fig. 2.2. Inside the outer boundary, the Finite Element Method is applied to solve the Helmholtz equation given by V [v(p)VO(p)] + ku(-p)q() = 0 (2.2) where =p) = E,(p) (2.3) v(7) = LP (2.4) u(7) =,r(p) (2.5) for E-polarization and S(p) = Hz(P) (2.6)

/ ino sBs Y Oo8 Ut...............~.... Figure 2.1: Geometry of the scatterer v(p) = E() (2.7) u(P) = Ar(P) (2.8) for H-polarization. Also, ko, = w(bjy; is the wave number, and Ar and E, are the relative permeablility and permittivity, respectively. The appropriate boundary condition is enforced on the surface of the impenetrable boundary, while the radiation condition is satisfied implicitly by evaluating the integral equation (.)=4'~ifC(-)-f {G(Tph)[~ a )] - (Qn)[ G('P)]} d' (2.9) on the boundary ra, where G(~,') = -l/H(2)(kolp-'I) (2.10)

1 r0 ^ + 3! I Figure 2.2: Partially discretized body is the 2-D Green's function in which ( denotes the zeroth order Hankel function of the second kind. Furthermore, p and' are the usual observation and source position vectors, respectively, and I2P'h= (x_')2+(y y)2 (2.11) The normal derivatives are taken in the direction of the outward normal of rid 2.1 Discretization of the Object and Field Quantities In fig. 2.2, rTp is the field/observation point boundary, and r, is the integration contour, which is placed midway between r. and rb. Also, rd denotes the contour contour, which is placed midway between F4 and rb. Also, Fd denotes the contour

Definitions for Finite Element Mesh Ng = total number of nodes in the finite element mesh N, = total number of elements in the finite element mesh N. = number of nodes on ra or rb along the x-direction Nu = number of nodes on ra or rb along the y-direction Na = total number of nodes on ra Nb = total number of nodes on rb Nab = Na + Nb La = Lrai Lb = ~?Irbi (za, Ya, )'- coordinates of a point on contour ra, (zb, Ybi) - coordinates of a point on contour rb, Table 2.1: Definitions for the finite element mesh enclosing the impenetrable portion of the scatterer. Herewith, each side of r,, rb aind r, will be numbered counterclockwise starting from the top side, as indicated in fig. 2.2. The fields in the region between ra and rd satisfy (2.2) in conjunction with the required boundary condition on rd. The boundary integral equation (2.9) will be enforced on ra. To numerically solve (2.2), it is required that the region within ra be discretized. This is done in a traditional manner when employing the finite element method. The 10

Definitions of Field Vectors (in terms of field unknowns at nodal points) Oi = fields corresponding to the nodes on Fa, iNi = fields corresponding to the first N, or N. (whichever is appropriate) nodes on Lb, ib, = fields corresponding to the nodes on rb, OI = fields corresponding to region I enclosed by rb and rd ld = fields corresponding to the nodes on the rd Table 2.2: Definitions of the field vectors global node numbering starts from the right endpoint of contour Lal and proceeds counterclockwise. The numbering continues beginning at the right endpoint of contour rb1 and proceeds counterclockwise. Within rb, the nodes are numbered from the lower left corner and proceed column -by column. The definitions pertaining to the FE mesh are given in table 2.1. Each node corresponds to an unknown field value to be determined. It is important to associate the unknown field values corresponding to the various node groups on contours ra and rL by using different variables The labeling scheme is given in table 2.2, and this discrimination of the nodal fields is required, since they are treated differently in the analysis. 2.2 Derivation of the Finite Element Matrix One of several approaches may be used to generate the finite element matrix, such as the variational approach or the method of weighted residuals. In this development, we will utilize Galerkin's method, which is a special case of the method of weighted residuals

with the distinction that the testing functions are the same as the basis functions. Proceeding with the finite element analysis, we may rewrite (2.2) as.dz( [v ( y) 5(X + )] + kau(x, y)(x, ) = 0 (2.12) the residual of which is given by R -T (Y) ( - k'u(x,y)k(xy) (2.13) az 1' ax I ay OY J The field within rF may be represented as a summation of piecewise continuous functions and, thus, may be written as Ne +(x, y) = Z 4,(X Y) (2.14) e=1 where 4e(x, y) is the field within element e. It is expanded as j —— 13 N(XYJ(2.15) j=1 where Ne's are the standard shape functions (found in any standard finite element book), S0's are the fields at the nodes, and n is the number of nodes per element. The weighted residual equation for the eth element is defined by Js NRdxdy = 0 i =...,n (2.16) where Se denotes the surface area of the eth element. Inserting (2.15) and (2.13) into (2.16) yields zJ]N [ A 7) (-V AZ ) } - A y -) ku],dzxdy = 0 12

Further, by using the identities O ( ON 0 ONN ON (21 ay v O y ( y y (y N (2.19) and the divergence theorem I(se -aa + -y dxdy = (u+ vy) iidl (2.20) where Ce is the contour enclosing the eth element, (2.17) becomes 1 il (ON? + N N j OqU C jdNdy Ox Ox ay ay Oy N v+, + y)'ndl i = 1,2,..., n (2.21) This may be written in mnatrix form as AC4 = b= (2.22) where ONe ONf ONe 0N ([Ac],j = || (v + v, i _ kCulsCNN) ddy (2.23) and {i= f Ne ( Z'- + v' i1,2,... n (2.24) For linear triangular elements, N', are given by N, 2= j(a? + bex + cy) (2.25) 13

and 11 1 -C 1det | 1 y y (b -c b-jc) (2.26) 1 x3 Y3 a! = xye - xkj (2.27) = y -e (2.28) c= X- (2.29) where (xz, yi) is the coordinate of the ith node of the eth element. From (2.25) Nr i (2.30) e - (2.31) Oy 20C Substituting these and the formula ] (NF)p(N2)qdxdy= 2'( p!q! (2.32) into (2.23), we find [Ae]ij = 4(be ( + ccj) - k i2Ue i (2.33) where 2 if i=j aifj = (2.34) 1 otherwise In (2.33) we have assumed that u and v (the material constitutive parameters) are constant within each element and are given by uc and vC, respectively. By summing over 14

all elements as implied by (2.14), we may write the overall system more explicitly as Aaa Aab 0 0 ba Aba Abb AbI 0 ||b 0 (2.35) O AIb AII AId' 0 0 0 AdI Add J d 0 In this, the values of the elements in the submatrix Apq are the contributions associated with the nodes in group p which are connected directly to the nodes in group q. One can easily show that the line integral contribution (see 2.24) of those elements vanishes everywhere, unless the element has a side in common with r.. As a result, be contributes only from the boundary ra of the finite element region, as indicated by the presence of the vector ba in (2.35). Without a priori knowledge of the total field on that boundary, ba cannot be determine. We may, however, provide the appropriate condition on this boundary by utilizing the integral equation (2.9). This amounts to replacing the first block-row of the matrix (associated with Oa on ra) with a discrete form of this integral equation. 2.3 Evaluation of Boundary Integral The boundary integral in (2.9) may be written as a summation of four integrals, one for each side of the contour r, as S(j0) = OknC(p) {Jrc [G(P,). ) (T') - O' P-G('Pl) dl, + j [G(P').n') + (7')G(,')] dl, + f| [GP,7 —. (p)+' + G(),-)1 dlC3 15n, O 15

d+ O [pP3) zG(p7p) dlC4 f (2.36) where the derivatives along the z and y directions, denoted by;- and -, respectively, have been left in this form for the later convenience of determining them numerically. More explicitly, we have (x, ya 1) = kinc(x, ya ) 3 3 d 6 -'jrc1G(x -Ic x,y )() -an' YW )- (x',Yc,) y-yG(x - x', y, y yc) dx' + I G(, xxc2, Y' (' ) + (x2y)G(x"xc2ya') dy C2 + [G(x - X' Ya i Yc3 (XYC3) + q(X', Yc3) -G(x - X',y,YC3)1 dx' F a + j [G(x, Zc4, 1,) y') (X, y') - Y (X4' y') G( c,,y4,')Yj y')]dy' C4 3 n.,, 3 (2.37) and O(x2, Y) oi= nTC(Xa2,IY) 4 4 wh-rethef G( t fnxa2 x' Y y r(', bc,) - c(x', yc) G(xs a2 x', re) dx' Co 4o Ofany 42' j G(Xa2,x,2 Y - Y')!) C(2 Y'x) + O(XC2 y')~/G(x,x)]1 C2 4 + j G[Tz 2 1, XY - Y') -O(XC4,Y') - 0(z IyY')-G(Xa2IXC47Y- y') dy e3 On, OX, 4 (2.38) where the first subscript on x or y refers to the contour (a, b or c), and the second refers to the contour number (see fig. 2.2 and table 2.1). It is noted that the arguments of the 16

Green's functions have been modified to imply a convolution when appropriate, and this representation will be used throughout. The normal derivatives of q may be evaluated via the central difference formulas n'(xcY') = [((aY - b, )] + O(A) (2.39) (Oz Yc = [(x', Ya) - (z', Yb)] + O(A2) (2.40) where A is the displacement of rl from rb (A is usually less than one tenth of a wavelength). Substituting (2.39) and (2.40) into (2.37) and (2.38), we obtain (x, ya1 ) = inc(zx, Yal ) 3 3 - {Jr [K; G(x - X' Ya1 Yc)~(X', Ya) - KIG(xz'Yai Yc,)4(X',Yb)] dx' + KI [KG(x Zxc2 Ya,l Y')(Z2,y) - K; G(x zcI x Ya, y')(xb.Y')] dy'' - c3 3 + j + [K (x - X'G, Ya,,YC3)k(X,ya3) - KG(x - X'Ya1 *YC3)(.Yb3)] dx' +'r [K- G(xZ Ic4 Ya, ) Y)O(xa4, y')- K.+ G(x, c4, Ya l Y')4(xb4* Y)] dy'} 3 3 (2.41) and O(Xa, y) = Oi'nc(xa, y) 4 4 - { [Ky G(Xa 2,X', Y. Y)(Z', Yal) - K+G(xa2 lX Y.Ycl)D(Z Yb)] dx' + G( [ a2,z,y - y')(zl2,y') - K G(za 2, 2,, Y - y')(zb2')] dy' + j [K (X',y, Yg X',ya3) -KG(xz.a2,'YYc3 )0('Yb3) dx' + c4 [K; G(xa2,cx4Y - y')t(Z4, y') - K+G(Xa2, c4 Y - Y')O(xb4,')] d'} (2.42) 17

in which K~ = zv + 28' (2.43) KV = -, -' (2.44) Assuming a pulse basis expansion for the nodal fields (centered at the nodal positions along contours ra and rb), a midpoint integration may be performed for the evaluation of the integrals in (2.41) and (2.42), to obtain O(Xi, Ya ) = "nc(Xi, Y ) 3 3 -iA E [KyG(xi -Xi, Ya,, Yc,)O(xj Ya -KV G(xi -ZXi, Ya1 Yc)O(XiYbi)},j2 jv,.=1 q 3 3 +[K+ —' — 3 A K-3 + A l K+cG(,~ s xc,,YiyI)~(a21yj )- K.- G(X- a I ty Yi-Yi)(Xb2 Yj ) j=1 18 (2.45) and 4 4

+ A K;-G(Xa,xc4, Xc, Y-yj)q(Xa),,j) x G - Y Y) (2.46) In these xi and yj denote the ith matching/testing points corresponding to the nodal locations on ra, while xj and yj denote locations on rb. We recognize some of the terms in (2.45) and (2.46) as discrete convolutions amenable to numerical evaluation via FFT. The subsystem (2.45) and (2.46) may be written more concisely as a) 1,iin' Sa T+ S+ R Tb4 1 al1 a al bi al b2 al b3 i, + S+ T+ S R i a2 b a2b a2b3 a2b q~ \q~,inc +b R T S+ Sa-3b I a3b2 a3b3 Ta3b4 3 qoh4 gfTaic S+ R1 T+ S+ 6, ra4 T ab a4b2 a4b3 a4b4 S + R, T+ arbi ab61 Q$b3 aR 4 ~ T+ S+ e i STc + a2b abd a2b3 oab4 p S+ RX T Sb3 T+ a3b2 -s abl b3 aab4 3ba T+:' S+ a4bl, Sasb2RY 77,b3 Sb, (2.47) with the various parameters to-be given explicitly later. The matrices RZ,,, simply reverse the order of the unknown vector so that the convolutions may be performed properly. This is required solely because of the employed counterclockwise nodal numbering scheme. Since i = 1,2,3, 4 ('i')last = (b)i first = (2.48) j = 2,3,4,1 19

the vector u[ bi 4 lb2 ib 4 (2.49) can be related to the actual unknowns on the contour rb via a transformation Db as 4b = D0b4 (2.50) and (2.47) may then be written as (I + Laa)ka - LabDbOb = nc (2.51) or [ABI[:] = n1fc (2.52) Lb where I+ Sa TSaI Sb3 R, T'a I + L,, = T2b + (2.53) L bl TSa Rvb Taa4b3 I+Sb4 and s+lbl T32 Sib3RX Tb2 b3 = S4T Lab = A Salb2 Ta2b3 a2b4RV T+ s _~ T+ S+3b61RR T,3b Sa3b3 T+3b4 T+ S+ carbl b S`L RV TT b3 Sa4b4 (2.54) 20

After replacing the first block row of (2.35) with (2.51), the complete system may be finally written as I + Laa -LabDb 0 O Oa 0 Sanc Aba Abb AbI 0 | b 0 (2.55) 0 AIb AII Aid OI 0 0 0 AdI Add J qd 0 to be solved via the CG algorithm. The Elements of ABI The elements of ABI defined above may be evaluated via the discrete Fourier transform. Specifically, we have Sa,bb = DFT- { DFT[G(. yal. Yb ) ~ G,(x, yY. a, )]DFT[,b ] (2.56) Sa3bi ~b = DFT-1 DFT[G(x, ya3 7 yb ) Gy(z, y,,3 Yb )]DFT[qb 1] (2.57) -3 3 3 3 3) Sab~ qb2 = DFT-1 DFT[G(Za2,Zb2, Y) G,(za2,zb2 Iy)]DFT[kb ]} (2.58) 4 4 4 Sab bi =DFT- 1DFT[G(aZb2,Y) ~ G(Za4,.Zb2,y)]DFT[b2l (2.59) 44 4 4 in which DFT denotes the discrete Fourier transform operator. Also G(ZZ, yyl) = 4 H(2)(kolp - p'l) = 4H2) (k z - )2 + (y - y2 (2.60) 21

G,(z,', y,y') = - -G(x,',y, y') j R}~Il) (ko. I-')2 + (y- y,) =- - Ak X /(z x')2 + (y y) -') (2.61) G(x x', y,y') =G(x,',y,y') 2 ay' - 1k( 1 (2) + ( y _ y/)2 8ko +(Y - Y') (2.62) Special cases of the convolution operators for the chosen mesh are given as Ga(x2,Xb2, y - y') = ~(-.) k Ib + 2 - Xb2 4 4(8o - b2 )2 + (y _- y')2 Za2 H(2) ko /( _ Z)2 + (ya, - Yb, )2 G,(z Z', Ya.Yb,)= A(- k)^ko 3 IYa, Yb, I 33 3 -- Yb2 3 3 { (2.64) Ya3 and the corresponding expressions for G are implied by the arguments in the previous two equations. The cross-term element submatrices are given by [T4 ] = G(xi, Z2, Yal,yj) Gx(zi, b2 z, Yaly) (2.65) 3 4 ij 4 3 4 3 [T b ] = G(Za2,zjyi, Yb)+ Gy(Za2,jYiYbl) (2.66) 4 ~3 i~ 43 4 3 with G.(xi, Zb2,yai, Yi) 4 3 (2)2(k = )2 + (Ya -yj)2x = ~( )ko (2 )2 zi - b2 (2.67) j(zi - Zb2)2 + (Yi - yj)2 4 ZXb 22

and Gy(Xa2,j, Yi, ybi ) 4 3 = +- 11 ko ~ /(a2 -j)2 + (yi Yb )2), - 3 -~(-2.)k ( 4+3)IN Yb 1A r Yb3 (2.68) 8 /a2 - zj)2 + (Yi - Yb1 )2 3 Yb{ where again the corresponding expressions for G are implied by the arguments of the previous two equations. Making the substitutions (Xi- Xb2 )2 = (i - _5)2z2 (2.69) 4 (ya _ -y)2 = j22 (2.70) 3 /(XZ- Zb )2 + (Ya -yj)l2 = i/(2 (2.71) and (xa2 - j)2 = j2A2 (2.72) 4 (Yi - Yb, )2 = (i- 5)2A2 (2.73) 3 X~za2 X /( -Z yjb)2 = j (i )2 (2.74) 3 we may write G. and G, as GG(~,,:, ybl, Y~yj) =:t:( )iko (2) + o, G-(zi, z2j, yi, y ba ) = (-( )k (2.75) 4 3V8' i -.5)2 + (-.75) GV(Z.2 Z 1i, Yi1 b) = i - i (k0 ( /2 + (i.51)i -.5 { (2.76) to be used in the actual implementation of the system. Since each of the above relations are similar, we are required to store only one of them and alter the signs accordingly. It should be noted that, however, this is not the most efficient method of storage. Storing 23

only a few of the cross term values and using an interpolation scheme will reduce the storage considerably. Of course, an interpolation table of (2.75) and (2.76) will lead to a substantial reduction in memory at the expense of some computational efficiency. Assuming that the positive sign is chosen in equations (2.75) and (2.76), we have T+1 = b 3 3 T+ -T+ a2b3 a2b3 4 4 T+ -fal b4 a1 b 3 3 (2.77) Choosing the positive sign for the (2.63) and (2.64), we also find s+ = S+ albl albl 3 3 s+ =s "361 a3bI 4 4 S+b =S 4 4 (2.78) Thus, the elements of ABI may be written as Sa- T+ T+ I + S aal, %lb b3 fb IT+ I = 5+I+ T2 S1 3b1R 7a3b I + Sa3b33 Ta-3b4 8Taibl Sarb T+ 3 I+ Sab 24

+ S-b R aib1 ~zab2 Sjb3 T +- S T R Laabbl =Sa2b2 a2 b3 a2 b4R, Lab = Sba3R;a'3b2 sa363 a3b4 a-b S R Tab3 S+b (2.79) The elements of the adjoint of ABI required in the implementation of the CG algorithm are I + ( )a' (T )a RT(S+ )a (T+ )a (I+ Laa)=)a I+(S ) (a32) R(Sa+b2)| RT)(S+ )a (ab3)a I +(Sa33 ) ( )a (' alb3 3 I (Sa3b a4 b) (LabDb)* = (|) (T )a (r)G R)a(S g)a D(7 T), (S) * + )a (;b )) R(S b2)a Initialize the residual and search vectors s = ) (2.2 )) a nT )"So -b) )b3 (S b )a," R (sT )a (sa )a+ ) (2-.0) 2.4 A CGFFT Algorithm The CG algorithm to be employed for solving the system (2.55) is as follows Initialize the residual and search vectors?~ = Il[~e1[ in" 0o IT 112=11 b11 (2.81) s = AO(~) (2.82) 25

r(l) = b-s s = Aar(l) (2.83) 7a = 11 112 (2.84) P(0) = r,1 (2.85) p() = O(O)S (2.86) Iterate for k 1,..., Ng s = Ap(k) (2.87) 7. = 119112 (2.88) or(k) = r1 (2.89) 4O(k+l) = (k) + (k)p(k) (2.90) r(k+l) = r(k) a(k)p(k) (2.91) 7r -rII P II~ (2.92) + = t11 (k+l) 112 (2.92) s = Aar(k+l) (2.93) 7. = 1I S112 (2.94) p(k) = a-l (2.95) p(k+l) p(k) + (k)s(k) (2.96) Terminate when k = Ng or ~ < tolerance. The individual operations associated with the ABI matrix-vector products are quite numerous and, therefore, will not be shown explicitly. However, it can be shown that the total system may be decomposed into a summation of two matrices; one involving 26

operators associated with the boundary integral and another involving the elements of the finite element matrix. The system matrix may then be written as {S} = (SBI} + {SFE}. (2.97) where I+ La LabDb 0 0 Z1 0 0 0 0 Z2 {sBI} = (2.98) 0 0 0 0 Z3 0 0 0 0 z4 and o o 0 0 z Aba Abb AbI 0 Z2 {SFE} = (2.99) 0 AIb AII Aid z3 o0 0 AdI Add z4 For the adjoint operations, we have I+ La 0 0 0 Z DbTLa 0 00 0 z2 (BI} = b a (2.100) 0 0 0 0 Z3 0 0 0 0 Z4 and 0 Aa 0 0 zI 0 Aa Aa 0 z2 {SFE} = 0 Ab Azb 0 (2.101) 0 AaI Aa1 At z3 0 0 A~ Ad d z4 27

In each case, the operation is performed such that only the resulting vector {s} need be stored. 28

Chapter 3 Computational Considerations This method is efficient in terms of memory usage and computation time. We discuss each of these aspects in detail below. 3.1 Storage Efficiency The fundamental advantage of this method is the reduction of storage requirements, so that the scattering by electrically large bodies may be evaluated. To show that the low storage requirement of O(Ng) is assured, we refer to tables 3.1 and 3.2. These contain a list of all major memory consuming variables. A summarized list is also given in table 3.3. Specifically, table 3.3 includes the memory requirements pertaining to the finite element matrix (FE), fast Fourier transfafm-s (FT), boundary integral cross terms (Cross) and conjugate gradient variables (CG). NC is one less than the number of elements connected to a particular node, and a typical value of 5 is used here. To put the quantities of table 3.3 in terms of N,, the total number of nodes, we consider two special geometries. The mesh in fig. 3.1 corresponds to a penetrable body, while that of fig. 3.2 corresponds to an impenetrable structure tightly enclosed 29

by the picture frame. Within each special case two extremes are considered; a mesh corresponding to a square object (N, = N.) and a long strip (N. >> Nv). In each case, Ng is assumed to be large. Alluding to table 3.4 the total storage is O(N9) for the square region, but is somewhere between O(Ng) and O(N2) for the (N. >> N,) case. This is based on the assumption that all cross terms are individually stored, but by using an interpolation table, the O(Ng) memory requirement can be assured regardless of the value of N, with respect to N,. In table 3.5, more dramatic results for the storage of the cross term are listed. In this case, all of the unknowns are on the outer two boundaries, so it appears that the storage is O(N2) for the square case. One must note, however, that the factor multiplying the Ng term may be quite small. The strip case, on the other hand, requires an O(N2) storage. This case would be an unlikely candidate for the use of this method, since that structure would be handled much more efficiently via a direct implementation of the CGFFT method. As noted above, the storage of the cross terms may be brought down to O(Ng) for all cases by using an interpolation table, and this will certainly be necessary in a 3-D implementation. 30

a Memory Consumption variable type count comment Mesh Xg R Ng X coordinate of global nodes Yg R Ng Y coordinate of global nodes Nglob I 3Ne Needed for FE matrix formation Pointers ABint I Nab Points to observation and integration points Pnodes I Pnum Points to nodes on conducting bodies dmatl I Ng - Nab Element material properties Finite Element Matrix (FE) Ar C (Nt2 )(NS - Na) Non-zero values of FE matrix col I (N )(Ng Na) Column pointer of FE matrix row I Ng - Na Pointer to rows of FE matrix Conjugate Gradient (CG) Phiz C Ng Unknown vector CG1 C Ng Residual vector CG2 C N, Search vector CG3 C N, Temporary vector q C MAX(N, N,) Temporary vector phiinc C Na Incident field vector Table 3.1: List of major memory-consuming variables 31

Memory Consumption (continued) code variable type count comment Fourier Transforms (FT) FTHxl C 2Nx Fourier transform along x-direction FTHx2 C 2N, FTHx3 C 2N, FTHx4 C 2Nx FTHy1 C 2Ny Fourier transform along y-direction FTHy2 C 2Ny FTHy3 C 2N, FTHy4 C 2Ny FT C 2MAX(Nx, Ny) FT of unknown sub-vector WR R 2MAX(N., Ny) Temporary array WI R 2MAX(N:, Ny) Temporary array Cross-Term Matrices (Cross) QP C ~ MAX(N,, N, ) PQm C MAX(N:, Ny) Legend R = REAL*4 C = COMPLEX I = INTEGER*4 32 Table 3.2: List of major memory-consuming variables (continued)

Figure 3.1: Example of the mesh of a penetrable structure Figure 3.2: Example of the mesh of an impenetrable structure 33

Major Memory Consumption (N., Nv) Item Type Count FE COMPLEX (N-2 )[Ng - 2(Nx + Ny)] FT COMPLEX 12N= + 8Ny Cross COMPLEX 2N2 CG COMPLEX 4Ng Table 3.3: Summary of major memory consumption Major Memory Consumption: Penetrable Item N = Nt J Nx >> NI FE ( 2 )( - 4V/g). (2)Ng(l - N+2 FT 208? 12N,/(Nv + 2) Cross 2N 2(N/Ny)2 CG 4N, 4Ng total ~ 9N I 2(N)2+ +7Ng Table 3.4: Summary of major memory consumption: filled mesh 34

Major Memory Consumption: Impenetrable Item N.=N NI Nx >> N FE ( )Ng/2 (I)Ng/2 FFT 5Ng/2 3Ng Cross N2 /32 N2/8 CG 4Ng 4Ng total N2N/32 + 8Ng j N22/8 + 17N2/ 2 Table 3.5: Summary of major memory consumption: open mesh 3.2 Computational Efficiency Since the primary importance of this method is storage reduction, a comparable level of efficiency with alternative methods is a bonus. A method for which a fair comparison may be made is the standard CGFFT. This requires a 2-D FFT, which is slower than using multiple 1-D FFTs for large bodies. We compared the two methods for a specific penetrable scatterer using an Apollo 3500 without code optimization. The pertinent CPU times are compared in table 3.-6. The comparison provides only a relative measure of the speed difference. 35

Il Body Properties F-CGFFT CGFFT II Composition Dimensions T/I (s) I Total T/I (s) I Total dielectric J 2A x 2A 8.6 155 1333 170 33 5610 E~ =4- j.1 Legend T/I = time/iteration I = number of iterations Table 3.6: A comparision of computation efficiency of the F-CGFFT with the CGFFT method 36

Chapter 4 Far Field Computation The scattered fields may be computed as qs(~) = -. {G(;,3') [-7Vn,(' )] -( ) [ } (4.1) Using the discretization scheme developed earlier, we have q(X, y) = - (Il [K;G(x, x'., yYC, )(' (xYal, )- K+G(xZ x', yYC, )(x' Yb)] dx' + | [K+G(x Zxc2, Y, Y')M(Zxa2,'Y') - K G(x, xc2s Y, y')(Zb2, Y')] dy' + j /KGz(x, x, Y 1 Yc)(', y - KG(x, x', Y Yc3 )(X', yb3)] dx' YKG(xxc Y) = |. fY-+J G(x, c, y, y')dy' (4.4) 37 21

7x(xyyc) y, y -)dy. (4.5) x(ZXcY7 ) = z+2 d 0G(x,x, y,y')dy' (4.6) (4.2) becomes Oq(x,y)= E 1 ([l3X(xy Yyc)_ (XyZ )YYcl)]{a)j - [Px(Z YYcl) + (Z. Y, Ycl)] {bl )j) Ny + E ([3Y(X, xc2 ) Y) + Y(x, xC2? Y)] {'Aa2 }j-[3y(x, c2, y)- y(x xC2, Y)] {fb2 }j) j=1 7 C + Z([13X(x, Y) + 7X(X, y, Yc3)].{a j - [3X(Xyc Yc3) - 7X(Zx,Y Y3 ) {fbb3 }j) Ny j=1+ E ( c4Y) /(C4 y()X a4 } ()j [ (X, XC4, Y) + (X XC4 Y)],O4 }j) (4.7) valid for all observation points (x, y). The specialize (4.7) to far zone computations, we must introduce the appropriate asymptotic expansion for the Hankel functions implied in (4.3)-(4.6). In doing so, we have x(x, Y, Yc ) = jfo(P)f1 (, Yc, )eikojxisO (4.8) 3 3 py(X xcCy) = jfo(P)f2(O, Xc2 )ejkoyjsin (4.9) 7'(x, y,y ) = -fo(p)fi(O,yc, )koA sin Oeiko-xosO (4.10) 3 3 yY(x,XC2,Y) = -fo(P)f2(,zc2: )koAcosOeikoysin6 (4.11) 4 4 where fo(p) = eJkp (4.12) 38

jkloyc1 sin8 fi(O,Ycl) = -e 3 (4.13) 3 j kox 2 cOs 0 f2(o, x2) = -e 4 (4.14) 4 in which (p, 9) imply the usual cylindrical coordinates of the observation point. Substuting expressions (4.8)-(4.11) into (4.7), we obtain b(x, y) = -fo(P) + Z([j + koA sincO] { }l - iij koA sinO] {bl }j) fi(O,yxc)eijkO jc j=1 + E ( - k0o sin 9] {qa'}, - Li + ko~ Csin 9] {(kb2 }i) f2( c2)eJkoyj sin j=1 Nx + ([j-koA sin 0] {fa3}j - [ + koA sin 0] {b3}j) fl (O, Yc3)ejik~i Cos~ Ny + ([j + ko4 cos ] { a4} -,j - k-ko A cos ] {]b4 }i) f2 (, x,)eJkoyj sin } j=l (4.15) The echowidth is defined by a = lim 2irp (4.16) - ir Ik'p incl2 and from (4.15) we have =- Z'I ([i + koA sin 9] {,}j- L[- koA sin ] {cbl }j) fl(,ycl )ejk'i J~ Ny + E (Li - koa COs 0]{ 2}j - [ + ko cos0] {4b2}j) f2(, Xc2)ejkyjsin j-=1 Nx + E ([ - koA sin o] {1a~} - U' + koA sin o] { b3}b) fl(o> YC3 ) eko~j co~I j=l Ny 2 + E (U' + ko Cos 10] {a4 }j - [- kos co~O] {1b44 }j) MfO XC4()e,)i si (4.17) 39

Chapter 5 Code Validation The scattering patterns from several rectangular structures are presented. The echo width is computed for each structure and compared to the results of the moment method. The bodies are discretized at a sampling rate of 20 samples/free-space wavelength. Results are presented for the following cases: * perfectly conducting bodies (figs. 5.1 and 5.2) * partially and fully coated perfectly conducting cylinders (figs. 5.3 - 5.8) * composite body (fig. 5.9) In each figure, the discretized geometry is shown along with the corresponding results. As seen in all cases, the generated patterns using the FE-CGFFT formulation are in excellent agreement with the moment method data.

.25 x 2. X Perfect Conductor (E-pol) 30.0 FE-GCFFTI' 20.0 - 4 MOM 10.0 - 0.0 -10.0 M -20.0 -30.0, 0.0 15.0 30.0 45.0 60.0 75.0 90.0 Angle (deg) Figure 5.1: Ez backscatter from a.25 x 2A body. 41

.25 x 2. X Perfect Conductor (H-pol) 30.0 FE-GCFFT 20.0 MOM 10.0 -10.0 mp -20.0 -30.0 0.0 15.0 30.0 45.0 60.0 75.0 90.0 Angle (deg) Figure 5.2: H,z backscatter from a.25 x 2A body. 42

.25 x 1. X Coated Perfect Conductor (E-pol) 30.0 FE-GCFFT 20.0 10.0 - 0.0 -10.0 mP -20.0 -30.0, 0.0 15.0 30.0 45.0 60.0 75.0 90.0 Angle (deg) Figure 5.3: E, backscatter from a.25 x 1A perfect conductor with a A/20 thick material coating containing the properties E, = 5. - j.5, -,. = 1. 43

.25 x 1. X Coated Perfect Conductor (H-pol) 30.0 FE-GCFFT 20.0 MOM 10.0 0.0 -10.0 go -20.0 -30.0 0.0 15.0 30.0 45.0 60.0 75.0 90.0 Angle (deg) Figure 5.4: Hz backscatter from a.25 x 1A perfect conductor with a A/20 thick material coating containing the properties E, = 5. - j.5,ls = 1. 44

.25 x 1. X Coated Perfect Conductor (E-pol) 30.0 FE-GCFFT 20.0 -------- OM 10.0 0.0 Cu -10.0 mQ -20.0 -30.0......... 0.0 15.0 30.0 45.0 60.0 75.0 90.0 Angle (deg) Figure 5.5: E, backscatter from a.25 x 1A perfect conductor with a A/20 thick material coating containing the properties Er = 5. - j.5, lr = 1.5- j.5. 45

.25 x 1. X Coated Perfect Conductor (H-pol) 30.0 FE-GCFFT 20.0,~~ ~ MOM 10.0 0.0 -10.0 m -20.0 -30.0 0.0 15.0 30.0 45.0 60.0 75.0 90.0 Angle (deg) Figure 5.6: H, backscatter from a.25 x 1A perfect conductor with a A/20 thick material coating containing the properties E, = 5. - j.5, /,r = 1.5 - j.5. 46

.25 x 1. X Double Topped Conductor (E-pol) 30.0.... FE-GCFFT 20.0 o MOM 10.0 0.0 =3 -10.0 m2 -20.0 -30.0............................ -90.0 -60.0 -30.0 0.0 30.0 60.0 90.0 Angle (deg) Figure 5.7: E, backscatter from a.25 x 1A perfect conductor with two A/20 thick top material coatings. The lower layer has the properties E, = 2. - j.5, it = 1.5 - j.2, and the upper layer has properties er = 4. - j.5,l = 1.5 - j.2. 47

.25 x 1. X Double Topped Conductor (H-pol) 30.0.....'.. i I FE-GCFFT'-~ 20.0 - 10.0 0.0 -~ -10.0 CP -20.0 -30.0 -90.0 -60.0 -30.0 0.0 30.0 60.0 90.0 Angle (deg) Figure 5.8: H, backscatter from a.25 x 1A perfect conductor with two A/20 thick top material coatings. The lower layer has the properties Er = 2. - j.5,ir = 1.5 - j.2, and the upper layer has properties Er = 4. - j.5, /l, = 1.5 - j.2. 48

.5 x 1. X Composite Body (H-pol) 30.0... FE-GCFFT mL, 20.0 OM 10.0 0.0 -10.0 m -20.0 -30.0 0.0 30.0 60.0 90.0 120.0 150.0 180.0 Angle (deg) Figure 5.9: XH scattering from a composite body. Both the perfect conductor and the dielectric body are A/2 in each dimension. The material properties are cE = 5. -j.5, Plr 1.5-j5.5. 49

Chapter 6 Conclusions and Future Work A procedure was developed for computing the scattering by 2-D structures. This procedure combined the boundary integral and finite element methods, and the resulting system was solved via CGFFT. The main advantage of the new approach was a reduction in memory demand to O(N) compared to O(N2) required with traditional solution techniques. A detailed map of the storage requirements was presented, and the principle memory consuming arrays were discussed. Also, the computational efficiency of the technique was examined and found favorable. To validate the proposed solution approach, several backscatter patterns were presented and compared with results based on traditional solution methods. A goal is to extend the technique to 3-D applications. In this case, the cross terms must be efficiently stored using an interpolation table to ensure an O(N) storage requirement. Also, the use of a simple boundary (as in [15]) in the application of the boundary integral equation would be desirable for additional storage reduction. Higher order elements are further of interest to increase the CG convergence rate. Second order elements are also within the solution domain of Maxwell's equations and would allow 50

a more accurate evaluation of the normal field derivatives. In addition, there are other numerical difficulties that must be addressed in 3-D applications. The modeling of the fields near corners of the scattererrequires some care (an obvious approach is to avoid placing a node at the corner location). Also, the field discontinuity at material transitions must be handled properly. The standard basis functions ensures continuity across a boundary, but this will require some modification.

Chapter 7 Program Manual In this section, brief descriptions of the pre-processing programs (mesh generators) and the main processing program are provided. They have been executed on an Apollo workstation an IBM 3090-600E and a Cray Y/MP. 7.1 Description of FECGFFT The main processing program, FECGFFT, is a menu driven program which allows the user to load the desired pre-generated mesh file, choose the type of computation (E- or H-polarization,backscatter or bistatic echo width), generate the desired data and store the resulting output. Some initial post-processing is also performed. For instance, if the near-field values on the grid are stored (this option is only available for bistatic computation), an additional file may be generated for a contour plot. The following menu is produced by FECGFFT during its execution. ****** Main Menu ****** Pre-processing (1) Load Finite Element Mesh File (2) Set new CG parameters 52

Analysis (3) E-polarization (Backscatter) (4) H-polarization (Backscatter) (5) E-polarization (Bistatic) (6) H-polarization (Bistatic) _ _ —— __- --— _ ----— _______ —-— ____......_ __-. — — _ —-- -- --— _________ Post-processing (7) Generate 3-D plot file......._ —--......._________________ Test Routines (10) Element node ordering (11) Test integral matrix: scattered fields (12) Free-space field comparision (20) Quit Item (1) allows the user to load a mesh data file generated from one the mesh generators to be discussed later. Actually, any mesh generator may be used, but the file must contain the correct information and format. This information can be found by examining the module MLOAD. Item (2) allows the user to change the CG residual error tolerance value and interval for printing the iteration number and the associated residual error. Items (3) and (4) are selected for the generation of backscatter data for E- and Hpolarization, respectively. When either of these is selected, the starting angle, stopping angle and angle increment will be prompted. The file name for the far-field data is also requested. A response of "none" will produce no file. A prompt for the pad size will then be requested. The suggested response of "1" will automatically determine the proper pad size. When the program has finished generating the desired data, it returns to the main menu. Items (5) and (6) are selected when bistatic data for E- or H-polarization are desired. 53

When either of these is selected, the incident angle will be prompted followed by the starting angle, stopping angle and angle increment. The file name for the far-field data is requested, followed by a file name for storing the nod'al field values. A prompt for the pad size follows as before. Item (7) allows the user to generate data in MPLOT format for contour plots. At the present time, only rectangular bodies will work for this option. Items (10)-(12) direct the user to test routines, not used for normal operation. The pertinent files which contain groups of subroutines associated with the accompanying description are as follows: file name brief description fecgfft.ftn main program fevec-sub.ftn - vector operation subroutines feio-sub.ftn file i/o routines fe-test3_sub.ftn near-field test routines fecross-sub.ftn cross-term subroutines fematrix.sub.ftn FE matrix routines fetest5_sub.ftn node order test routine fe fft-sub.ftn FFT routines fethree-sub.ftn three-dimensional plot data generation fefield-sub.ftn near/far field computation feftestsub.ftn free-space test routine fe-summary.sub.ftn generates session summary 54

These files must be compiled and linked prior to execution. 7.2 Mesh Generator for Curved Bodies This program is under development for various specific types of bodies. It may, however, be used to generate a mesh for virtually any desired body. The mesh generation is accomplished by first dividing the region between the impenetrable surface (if any) and the rectangular enclosure into first, second or third order serendipity elements. These are subsequently mapped to a square domain, subdivided and mapped back. Examples of these are shown in fig. 7.1. It works well for modelling regions with curved boundaries, but generally produces more unknowns than necessary for the solution method. An input file to this program may be generated either with option (2) from the main menu, or manually. Selection of option (3) from the main menu processes this file and places the results in a specified output file. The output file is then used as input to the program FECGFFT. Running the program produces the following menu: Main Menu (1) Preferences (2) Create an input file (3) Process an existing input file (5) Plot routine (10) Quit Item (1) has not been incorporated as of yet. The selection of Item (2) produces the menu: 55

4 1 \ (D 2 4 0i 2 4 10 12 73 Figure 7.1: Typical serendipity elements used in the region descretization process 56

Input File Creation Menu (1) Conducting Strip (2) Rectangular Coated Slabs (3) Coated Ogive (4) Circular Cylinder (10) Return to Main Menu Only Items (3) and (10) are operational at this time. Selection of (3) yields the menu: Ogive Menu (1) Enter geometry (2) Modify geometry (3) RETURN Choosing Item (1) results in a series of prompts outlined as follows: 1. a,b for the coating, where a = height of the arc and b = half-length as indicated in fig. 7.2 2. the relative permittivity and permeability of the coating 3. a,b for the conductor 4. sampling interval (in wavelenths) at the integration boundary 5. number of circumferential samples in the free-space region 6. number of circumferential samples in the material 7. arc distance along coating 8. arc distance along conductor 9. comment to appear in the input file 57

a H- 2b " Figure 7.2: Arc used for part of an ogival structure. 10. input file name Upon completion of the session, a file is generated to be used at the input to FECGFFT. Selecting Item (3) results in a prompt for the input and output filenames. Manual generation of-the input file requires that the scattering body be surrounded by a rectangular boundary displaced approximately one element from the body. The region within the rectangular boundary and the impenetrable body surface (if present) is then subdivided into either linear, quadratic or cubic elements, examples of which are given in fig. 7.1. Note that every node and side of each element is numbered as indicated. The output file contents are listed as indicated in tables 7.1 and 7.2 with the variable and descriptions in table 7.3. The first four lines are self explanatory. The next group of lines contains the coordinates of the nodes, and the order of these pairs determines the global node numbering scheme. Two real numbers followed by a "C" are assumed to be in cylindrical (r,6) coordinates centered at the previously specified value. An "N" following the "C" will change the center coordinates to xc, yc. Immediately after the node coordinates definition, the elements or "blocks" are defined. The local/global node 58

relationship defines the block. This format is repeated for each of the Ne elements. The elements must then be connected by specifying the sides of adjacent elements. This avoids the time-consuming task of comparing the coordinates of every node to the others for spatial commonality. The impenetrable and integration boundaries designation are present for a similar reason. Finally, the material properties are requested. The order of these determines the number to be used in determining the element material property number. The free-space value is always present. 59

Line number Contents of Line 1 Comment 2 Ne, 1, 0 3 fi, f2, f3, f4, f5 4 Nn 5 x,y (or r,.,C or r,,CNx, yc) Nn+4 x,y (or r,.,C or r,.,CN,x, yc) N,+5 Comment (1-st Block) Nn+6 - 0, 0, M1 Nn+7 N,I, Na1, 01 Nn+8 1., 1. Nn+9 L G(1), L G(2),, L G(O1) Nn + Ne+3 Comment (N,-th Block) N, + Ne+4 0, O, MN8 Nn + Ne+5 NaNe, N.Ne) ONC Nn + Ne+6 1., 1. Nn + Ne+7 L * G(1), L 4 G(2),..., L * G(ON,) Table 7.1: Beginning portion of the input file format 60

Line number Contents of Line Nn + Ne+8 Comment -(Connection) Nn + Ne+9 Nc Nn + Ne+10 ei, sj, ek, sl Nn + Ne + N0+9 ej, sj, ek, sl Nn + Ne + Nc+10 Comment (Conducting Boundary) Nn + Ne + N+11 Nbc Nn + Ne + Nc+12 ei, sj Nn + Ne + Nc + Nb+11 ej, sj Nn + Ne + Nc + Nb+ 12 Comment (Integration Boundary) Nn + Ne + Nc + Nbc+13 Nbi Nn + Ne + Nc + Nbc+14 ei, sj Nn + Ne + Nc + Nbc + Nbi+13 ei, sj Nn + Ne + Nc + Nbc + Nbi+14 Comment (Material Property Table) Nn + Ne + Nc + Nbc + Nbi+15 Np N, + Ne + Nc + Nbc + Nbi+16 Re{ E}, Im{Er}, 7Ze{lrp}, rZm{pu} N, + + Ne + N + Nbc + Nbi + Np+15 |Z e{fEr} m{E}, ZeIm{,}, e{Im{L,} Table 7.2: Remaining portion of the file input format 61

variable type description Ne I total number of modelling elements Af I put element numbers on mesh (l=yes, 2=no) f2 I put node numbers on mesh (1=yes, 2=no) f3 I put material numbers on mesh flag (1=yes, 2=no) f4 I surround mesh with scale (1=yes, 2=no) f5 I generate PostScript file (1=yes, 2=no) Nn I total number of nodes (x, y) R,R cartesian coordinate pair (r, 9) R,R cylindrical coordinate pair (xc, yc) R,R coordinates of arc center Mi I material number of ith element Nsi I number of samples on side i of element j Oi I order of ith element L - G(i) I global node number of ith local node Nc I total number of element sides in contact ei I element number St I side i of element j Nbc I number of elements adjacent to conducting boundary Nbi I number of elements adjacent to integration boundary Np I number of entries in material table Table 7.3: Variable description and FORTRAN declaration type. 62

Figure 7.3: The mesh of a conducting circular cylinder. Example The following file is an example of a perfectly conducting cylinder of radius.5 A in free space. A figure of the resulting mesh is shown in fig. 7.3. Circular cylinder 8 Aug 1989 8,1,0, 2,2,2,2,2 40,.5,90.,CN,O.,0..5,112.5,C.5,135.,C.5,157.5,C.5,180.,C.5,202.5,C.5,225.,C.5,247.5,C.5,270.,C.5,292.5,C.5,315.,C.5,337.5,C.5,0.,C.5,22.5,C.5,45.,C.5,67.5,C 63

.525,90.,C.64,135.,C.525,180.,C.64,225.,C.525,270.,C.64,315.,C.525,0. C.64,45.,C.55,.55.275,.55 0.,.55 -.275,.55 -.55,.55 -.55,.275 -.55,0. -.55,-.275 -.55, -.55 -.275, -.55 0.,-.55.275,-.55.55,-.55.55,-.275.55,0..55,.275 1-ST BLOCK 0,0,1, 11,3,8, 1. 11., 1,15,25,27, 16,24,26,17 2-ND BLOCK 0,0,1, 11,3,8, 1.,1., 3,1,27,29,2,17,28,18 3-rd BLOCK 0,0,1, 3,11,8, 1.,1., 31,5,3,29,19,4,18,30 4th BLOCK 0,0,1, 3,11,8, 1.,1., 33,7,5,31,20,6,19,32 5th BLOCK 0,0,1, 11,3,8, 1.,1., 33,35,9,7,34,21,8,20 6th BLOCK 0,0,1, 11,3,8, 1.,1., 35,37,11,9,36,22,10,21 7th BLOCK 0,0,1, 3,11,8, 64

1.,1., 11,37,39,13,22,38,23,12 8th BLOCK 0,0,1, 3,11,8, 1.,1., 13,39,25, 15,23,40,24,14 CONNECTION 8 1,1,2,3 2,1,3,4 3,2,4,4 5,1,4,2 6,1,5,3 7,2,6,3 7,4,8,2 8,4,1,3 Perfectly Conducting Boundary 8, 1,2 2,2 3,3 4,3, 5,4, 6,4 7,1 8,1 Integration Boundary 8, 1,4, 2,4 3,1 4,1 5,2 6,2 7,3 8,3 Dielectric property table Epsr, Epsi, Mur, Mui 1 The program contains the following files: 65

file name Jbrief description gen.ftn main program gensub.ftn contains associated subroutines fe.grid.ttz.sub.ftn for plotting a mesh with rectangular elements on the Apollo screen using graphics primatives fegrid sub.ftn for plotting a mesh with triangular elements on the Apollo screen using graphics primatives fe-post-sub.ftn for generating a postscript version fegridsub.ftn These programs should be compiled with the SAVE option and linked before execution. 7.3 Mesh Generator for Rectangular Bodies This program is useful for Fenerating the mesh associated with coated rectangular bodies. Executing the program produces the following menu: **** Mesh Generation Menu **** (1) Conducting Strip (2) Composite Bodies (3) View an existing file (10) Quit Item (1) should be ignored, since the data file it generates for the strip does not distinguish between the nodes above and below the strip. It is thus invalid for H66

polarization computations. Item (2) allows the user to generate a rectangular composite body. Upon its selection, the user is promped for the following items: 1. the size of the square building block cell (in wavelengths) 2. the permittivity and permeabilities of the various material layers for the structure 3. the length and width of the main scattering structure in integer multiples of the initially specified building block size in 1. 4. the type of scattering body (conductor or material) 5. the number of layers for each side plus the material number from the material table generated in 2. 6. the number of cells between the scattering body and each of the four boundaries (usually 0, unless a larger grid is desired) 7. the name of the output file to be used by the FECGFFT program

Figure 7.4: The mesh of a rectangular partially coated cylinder. Example To generate a.5A x 2A conducting body with a.1A material coating of Er = 5. - j.5 and eu, = 1.5 - j.1 over sides (1) and (2) (see fig. 7.4, the following output is displayed: **** Mesh Generation Menu **** (1) Conducting Strip (2) Composite Bodies (3) View an existing file (10) Quit 2 Enter del (size of building block) in wavelengths.05 Enter dielectric materials to be used [(-1.,0.) to quit] (Remember, Imaginary parts <=0.) Epsilon 1 = (1.000000,0.0000000) Mu 1 = (1.000000,0.0000000) Enter Epsilon 2 (5.,-.5) Enter Mu 2 (1.5,-.1) Enter Epsilon 3 (-1.,o.) Enter length and width of main body in units of del 40,10 Main body composition: (0) Conductor (1) Dielectric 0 68

21 14 I I 3 Enter Onumber of dielectric layers sides 1, 2, 3, 4: 2,2,0,0 Index Epsilon Mu 1 1.000 0.000 1.000 0.000 2 5.000 -0.500 1.500 -0.100 Material property number for: Side 1 layer 1 2 Material property number for: Side 1 layer 2 2 Index Epsilon Mu 1 1.000 0.000 1.000 0.000 2 5.000 -0.500 1.500 -0.100 Material property number for: Side 2 layer 1 2 Material property number for: Side 2 layer 2 2 Number of rows and columns of blank cells surrounding the body 0,0 Generate PostScript file? (l=yes, 2=no) 2 Enter file name for data storage testout **** Mesh Generation Menu **** (1) Conducting Strip (2) Composite Bodies (3) View an existing file (10) Quit The output file testout contains the mesh information required by FECGFFT. Item (3) provides for viewing the plot on an Apollo screen. Upon its selection, the user will be prompted for a file name. Entering the interactive mode results in the following menu: Interactive Mode Menu 69

(1) Max and Mins (2) Picture orientation (3) Picture size (4) Picture offset (5) Tick spacing (6) Legend contents (7) Legend offset (8) Legend label size (9) Label contents (10) Label size (11) Number size (12) Number format (13) Print option flags (14) View on screen (15) Get hard copy (16) PsPreview hard copy (17) Reset default values (20) Return to main menu Currently, options (6)-(10) have not yet been incorporated. The remaining items are self-explanatory. The programs contains the following files: file name brief description mgenlin nc-new2.ftn main program fegrid.sub.ftn for plotting a mesh with triangular elements on the Apollo screen using graphics primatives fepostsub.ftn for generating a postscript version fegrid-sub.ftn which should be compiled with the SAVE option and linked before execution. 70

Bibliography [1] M Swaminathan, E. Arvas and T.K. Sarkar, "A survey of the various computer architectures for solving large electromagnetic field problems," 1989 APS/URSI Symposium, San Jose, CA. [2] T. Cwik and J.E. Patterson, "Solutions of large or time consuming problems on the hypercube computer,"1989 APS/URSI Symposium, San Jose, CA. [3] L.A. Takacs and V.P. Cable, "Low and High frequency electromagnetic methods and calculations on a connection machine," 1989 APS/URSI Symposium, San Jose, CA. [4] K.D. Tatalias, "ParalleL computing for the electromagnetic analyst,"1989 APS/URSI Symposium, San Jose, CA. [5] J.E. Patterson, T. Cwik, R.D. Ferraro, W.A Imbriale, N. Jacobi, P.C. Liewer, T.G. Lockhart and G.A. Lyzenga, "Computation of electromagnetic scattering on a Hypercube,"1989 APS/URSI Symposium, San Jose, CA. [6] A.T. Perlik, "Computational electromagnetics using a connection machine," 1989 APS/URSI Symposium, San Jose, CA. 71

[7] S.D. Gedney and R. Mittra,"Solving electromagnetic scattering problems via the method of moments on the connection machine," 1989 APS/URSI Symposium, San Jose, CA. [8] T.K. Sarkar, E. Arvas and S.M. Rao, "Application of FFT and the conjugate gradient method for the solution of electromagnetic radiation from electrically large and small conducting bodies," IEEE Trans. Antennas Propagat., vol. AP-30, pp. 409-418, May 1986. [9] T.J. Peters and J.L. Volakis, "On the Formulation and Implementation of a Conjugate Gradient FFT Method," submitted to J. Electromagnetic Waves and Applications, 1989. [10] T.J. Peters, "Computation of the Scattering by Planar and Non-planar Plates using a Conjugate Gradient FFT Method," Ph.D. dissertation, Radiation Laboratory, The University of Michigan, 1988. [11] Jian-Ming Jin and Valdis V. Liepa, "Application of hybrid finite element method to electromagnetic scattering from coated cylinders," IEEE Trans. Antennas Propagat., vol. AP-36, pp. 50-54, Jan. 1988 [12] P. Silvester and M.S. Hsieh, "Finite-element solution of 2-dimensional exterior field problems," Proc. Inst. Elec. Eng., vol. 118, pp. 1743-1747, Dec. 1971. [13] B.H. McDonald and A. Wexler, "Finite-element solution of unbounded field problems," IEEE Trans. Microwave Theory Tech., vol. MTT-20, pp. 841-847, Dec. 1972. 72

[14] K.K. Mei, "Unimoment method of solving antenna and scattering problems," IEEE Trans. Antennas Prop., vol AP-22, pp. 760-766, Nov. 1974. [15] K.L. Wu, G.Y. Delisle, D.G. Fang, and M. Lecours, " Waveguide discontinuity analysis with a coupled finite-boundary element method," IEEE Trans. Microwave Theory Tech., vol. MTT-37, pp. 993-998, Jun. 1989.

UNIVERSITY OF MICHIGAN 3IIIl l I1II5lI III4 111111161111110' 3 9015 02947 5160