S) (. . .I' 7 1-]; If”. 5" .2 .1 J. 4 IllllllllllllllllllllllllllllllllllllllllllHlllllllllllll 1 31293 00788 58 52 LIBRARY "uh“!!! State l University L __ * This is to certify that the dissertation entitled On W Effie/(lam MQH/Lool with l/«chfJf {36198 filming presented by Honj 2 hanfi 34% has been accepted towards fulfillment of the requirements for PE» ~ D degree“, Apnhécgh Matl lwnah CS Date 2’3 (1 5 7 MSU i: an Affirmative Action/Equal Opportunity Institution O~12771 PLACE IN RETURN BOX to remove this checkout from your record. TO AVOID FINES return on or before date due. DATE DUE DATE DUE DATE DUE |C____________ fl(_ #%= MSU Is An Affirmative Action/Equal Opportunity Institution emu“: ON THE GALERKIN METHOD WITH VECTOR BASIS FUNCTIONS By Hong Zhang Sun A DISSERTATION Submitted to Michigan State University in partial fulfillment of the requirements for the degree of DOCTOR OF PHILOSOPHY Department of Mathematics 1989 C»? u (984 3601 ABSTRACT ON THE GALERKIN METHOD WITH VECTOR BASIS FUNCTIONS By Hang Zhang Sun The Galerkin method with new vector basis functions has been used by researchers in electrical engineering. Many experiments show that the new vector basis functions are ideally suited for representing surface electric current and surface magnetic current on the triangulated surface of the given dielectric scatterer. The efficient and simple numerical algorithms have been developed. In this work, a theoretical analysis on the convergence, error estimate and numerical stability of this method is given for a certain types of electric field integral equations. Finally, the formulation of coupled surface integral equations and corresponding matrix equation is demonstrated and numerical results are presented for the case of a homogeneous dielectric sphere. ACKNOWLEDGEMENTS I would like to thank Professor Tien Yen Li, my thesis advisor, for his invaluable inspiration and guidance throughout my graduate study. I would like to acknowledge all the faculty members and students who gave me help and assistance during my studying at Michigan State University. In particular, I would like to express my appreciation to Ms. Xiao-Yi Min, for providing me the paper related to this work and her insightful comments. Also I would like to thank Professor Richard 0. Hill and Clifford E. Weil who took time to read my drafts and give helpful suggestions. I am very grateful to my husband, Xian-He Sun, my parents and my sister for their love and years support; and to my wonderful son Alan for being so lovable. TABLE OF CONTENTS Introduction .............................................................................................................. 1 Chapter 1 Basis Functions ..................................................................................... 4 Chapter 2 Assumptions and Projection Method .................................................. 10 Chapter 3 Convergence and Error Bound ........................................................... 23 Chapter 4 Computational Stability ....................................................................... 26 Chapter 5 Numerical Results ................................................................................. 31 5.1. The formulation of coupled surface integral equations ................................ 31 5.2. Matrix equation ............................................................................................. 36 5.3. Numerical results .......................................................................................... 37 Summary .................................................................................................................. 41 Bibliography ............................................................................................................. 43 ii LIST OF FIGURES 1. Finite circular cylinder dielectric scatterer ............................................................ 5 2. Local coordinates associated with an edge ............................................................ 6 3. Three edge currents associated with a triangle ...................................................... 7 4. Geometry of vectors to centroids of triangles associated with an edge ................ 9 5. Design of triangular patches .................................................................................. 13 6. Triangle T and T’ ................................................................................................... 17 7. Geometry of a homogeneous lossy dielectric scatterer in an isotropic free space medium ......................................................................... 32 8. Triangular surface patching for top half of sphere ................................................ 37 iii LIST OF TABLES Table l ....................................................................................................................... 38 Table 2 ....................................................................................................................... 39 Table 3 ....................................................................................................................... 40 iv INTRODUCTION There is a very broad class of practical problems in engineering electromagnetics which are of the boundary value type for perfect electrically conducting bodies. Such problems generally require the determination of various electrical quantities for a con- ducting body and are called electrodynamic problems. In such problems, the current dis- tribution is often the fundamental quantity. Knowledge of the current distribution on the body is sufficient to determine other electrical quantities such as the charge distribution, scattering pattern or radar cross-section. We usually call this type of problem a scatter- ing problem. A rapid expansion, both in computer capabilities and in the availability of efficient and highly stable numerical techniques, has resulted in an increasing interest in develop- ing computer codes for heating electrodynamic problems. The most notable step in this direction is using the electric field integral equation formulations, both in the time and frequency domains, in which the body is modeled as a surface patch model. The surface patch model is easy to describe to the computer and a number of methods have been developed for two and three-dimensional scattering problems. Among them a simple efficient Galerkin procedure with two dimensional vector basis functions defined on tri- angular surface patches or three-dimensional vector basis functions on tetrahedral volume elements has been used by many researchers. In this work a theoretical analysis on the convergence, error estimate, and numerical stability of this procedure is given for l a certain types of electric field integral equations. The analysis can be applied to other types of electric field integral equations in both two—dimensional and three-dimensional cases similarly. Usually an electric field integral equation is formulated as [TJ(r) llm=E‘(r)|.... re s. where T is a linear integral operator, J is an equivalent electric surface current, E" is an incident electric field, and S is a regular two-dimensional region. The subscript "tan" refers to the S tangential component. Under appropriate assumptions this integral equation is well-posed. The equation is then solved by the Galerkin method as follows: 1) Triangular elements are used to model a scattering two-dimensional body in which the electrical parameters are assumed constant in each triangular patch. 2) Special vector basis functions { f" 191:1 are defined on the triangular elements to insure that the normal electric field satisfies the correct jump condition at interfaces between different dielectric media. 3) A test function i: z 1,, r, is substituted into the integral equation to obtain Tj=E‘+ERR. Setting < ERR, r", > =0 (m =1, .. .,N), (here <,> denotes the inner product for the Hilbert space under consideration), the integral equation is converted to a matrix equation: ZN I = V , where ZN is an NxN matrix with < Tr, , r, > as its (m, a)“ entry, I = (11, . . . , IN )T is an unknown vector of coefficients of the test function 3, and v: (< E‘,f1>,...,< E‘ , fN >)T is aknown vector with < E‘ , r, > as its nth com- ponent. Solving this system numerically, an approximate solution JN of the elec- tric field integral equation is obtained. This work is organized as follows. Chapter 1 describes vector basis functions and their properties. In Chapter 2, the projection method (the Galerkin method is a special case) is introduced. The definition of Galerkin projection used in this work and its pro- perties are given. Chapter 3 shows the nonsingularity of the matrix ZN , the convergence and the error bound of the approximate solution to the exact solution. In Chapter 4, the bound of the condition number of the matrix ZN is derived, which shows the stability of the numerical computation. In Chapter 5, the formulation of coupled surface integral equations and the corresponding matrix equation is demonstrated. Then the numerical results are presented for the case of a homogeneous dielectric sphere. CHAPTER 1 Basis Functions In this Chapter we describe a set of vector basis functions, originally proposed by Glisson[1], which are suitable for use with the integral equation under consideration and triangular patch modeling. To begin with, an arbitrary two-dimensional body Sin R3 is modeled by triangular patches. In this context, S is sometime called the scatterer. The body is assumed to be connected, orientable, of finite extent, and composed of non-intersecting surfaces. In general, a triangulated surface modeling an arbitrary body consists of planar triangular faces, vertices, and edges. These geometrical elements are illustrated in Fig.1. Note that real world problems are often defined over bodies with curved boundary or curved sur- face. They are not amenable to planar triangular modeling. In practice, we usually assume that the geometrical discretization error can be decreased by reducing the size of the elements and ignoring misfits. Eventually, we shall restrict our attention to a plate S, that can be modeled by planar triangular patches perfectly. Once the two-dimensional scatterer has been appropriately modeled by uiangular patches, the edges of the triangles are of primary importance for the development of the basis functions. Fig.2 shows two triangles, T}; and T7,, associated with nth edge of a tri- angulated surface modeling a scatterer. Points in T; may be designated either by the position vector r with respect to the origin 0, or by the position vector p: = p: (r) defined with respect to the free vertex of TI. Similar remarks apply to the position vector p; 4 x it A Triangular urfaee Patch Vertex 1 Edge / % S ’ / I 3 / TRIANGULATED SURFACE z . E k Hi 0 ; Fig.1. Finite circular cylinder dielectric scatterer. nth edge :1 Figure 2. Local coordinates associated with an edge. ' except that it is directed toward the free vertex of T;. The plus or minus designation of the triangles is determined by the choice of a positive current reference direction for the nth edge which is assumed to be from T; to T;. The electric current flows along the radial direction p: in the triangle T: and similarly flows along the radial direction p; in the triangle T;. We define a function associated with the nth edge as Apt/2A}, forallrinT; f,(r)=< 1,. p;/2A;, forallrinT; butnotonthe nth edge (1.1) 0, otherwise where I,I is the length of the edge and Afi is the area of the triangle Ti. (Note that we use the convention, followed throughout this work, that subscripts refer to edges while superscripts refer to faces). It is easy to verify that { f,- } are linearly independent. The unknown electric current can be approximated by a N ‘10-) = Z In qu‘) - (12) 3-1 The summation is over the N edges that make up the triangular model of S. The basis functions f,I (r) have several properties that make them useful for representing 3 (r). Constant Vectors 12f141lello® l3fj-llf3 II t0® 0 Fig.3. Three edge currents associated with a triangle. 1) Within each triangle 3(r) is the sum of three basis functions that are associated with the three edges. Fig.3 shows that the superposition of the basis functions with a tri- angle conveniently represents a constant current flowing in an arbitrary direction 2) 3) 4) 5) 6) within the triangle. At each edge, f,I (r) has no component normal to that edge except at the nth edge. The component of f, (r) normal to the nth edge is constant and continuous across the edge because the normal component of pi along the nth edge is just the height of T2: with nth edge as the base and the height expressed as (2.43 )/l,.. This latter fac- tor normalizes f,,(r) in (1.1) such that its flux density normal to the nth edge is unity, hence ensuring continuity of the component of f,(r) normal to the edge. The surface divergence of the basis function is f InlAt, rintheinteriorofT; V-f,(r)=+—l,,/A;, rintheinteriorofT; . (1.3) 0 , otherwise where the surface divergence in T3 is (it/pi) 3(pf f,,)/api’.‘ (f. is the component of f, in the direction of p"). The charge density is thus constant in the interior of each triangle, and the basis functions for the charge evidently have the form of pulse doublets. The moment of f, is given by (A: + A;) if." where I. (Animist: I f.ctv=-2—(p:++pf:) (1.4) =In(rn--rft+) and of is the vector from the free vertex to the centroid of Ti with p? directed away from the vertex and p.‘ directed toward the vertex, as shown in Fig.4, and rfft is the vector from the origin 0 to the centroid of 7?. If edge n is on the boundary of S, then only one of the triangles, T; or T;, is interior to S. In this case it is assumed that f, is defined over only the triangle in the interior :1 :1 ii" + 92‘) Fig.4. Geometry of vectors to centroids of triangles associated with an edge. of S. An important interpretation of (1.2) that follows from 1) and 3) is that the expan- sion coefficient 1,. represents the normal component of j (r) at the nth edge. CHAPTER 2 Assumptions and Projection Method Projection methods are a class of techniques for obtaining an approximate solution to an operator equation T x = y , (2.1) where T: X -) X is a bounded linear operator with bounded inverse, X is a normed linear space. The approximate solution it” of (2.1) is an element of a finite-dimensional subspace XN of X satisfying PN T XN = PN y , (22) where PM is a bounded projection operator (a linear idempotent operator) on X with range X”. From this general scheme we obtain the Galean method. Now we want to develop a process yielding the approximate solution of our electric field integral equation by this method. First, we need the following preliminary assump- tions: 1) In equation (2.1), X = L2(S), where S is a plate that can be modeled by triangular patches perfectly, i.e., 10 2) 3) 11 L2(S)= x: S 4C2, II x(r)|2 dr=jx(r)-'y'(r)dr, x, ye L2(S). s The operator T is an operator of form T = E + K , where E is the identity operator and K is a completely continuous integral operator. Thus equation (2.1) becomes (E+K)x=y, x,yeL2(S). (2.3) It is a Fredholm equation of the second kind. Here and throughout, || M denotes the L2 norm. There is no solution of the equation Tx=0 other than x = 0. It follows from the Fredholm Alternative that '1"1 exists and is bounded. Therefore equation (2.3) is well-posed. For all the triangulations of S that we consider, there are constants co and 90 with 0 < 90 < 1t/2, such that Ig/ljSCo and 9059? Sit-90, where I; is the length of the ith edge, and 6? is the angle in the triangle T? opposite 4) 12 the ith edge. For each triangulation I, let I 1 = max 1;. It follows that I -L '11 Sli S11 , Co sine 2° 4% SA?E 51% . (2.4) 260 where A? is the area of the triangle T?. First, let GN be an N xN matrix with the inner product of the basis functions < f; , fj > as its (j, i)"' entry. It is easy to verify that On =12 I'N . (25) where I =11, TN is a positive definite matrix which depends only on the shape of triangular elements. Throughout the approximation process we attempt to keep the elements in the same shape. This means that we only change the size of the patches. The shape of elements is designed to satisfy 5 min l-S max is" forallN, 2.6 60 1.60mi) ‘ 1.6007») p ' 80 ( ) where 80, 80 are some positive constants, and 0(1),) is the spectrum of the matrix F”. Thus conda‘N) s 80 / so for all N , (2.7) where cond (TN) denotes the condition number of 1"”. For example, cond (1),) S 7/3 and cond (FN) < 4 if we choose triangles as shown in Fig.5 (a) and (b) respectively. We shall see in Chapter 4 that this is one of the conditions that guarantees the stability of the numerical computation. 13 (a) 0)) Fig.5. Design of triangular patches. 5) For a fixed triangulation 1:, define N XN = {z I,l f" | 1,, (n =1,..., N)arecomp1ex numbers}, i=1 where N is the number of edges of the triangular model. XN is an N-dimensional 14 linear subspace of L2(S). From previous assumptions, the projected equation(2.2) becomes (E+P~K)x~=PNy, (2.8) where y e L2(S), X” e X”. For the Galean method a specific projection operator PN is defined as follows: N Definition. For each x in L2(S), let PN x = z 1,, f,,, where the coefficients n=l I = (11, . . . , IN )T are determined by the solution of the linear system GNI=(,...,)T. (2.9) Since GN is a positive definite matrix, equation (2.9) has a unique solution for any x e L2(S ). Thus PM is a well defined operator from L2(S) onto the N-dimensional sub- space XN. Also it is easy to verify that IPA; is a projection with the properties: (1) =(m=1,..,N)forallxe L2(S). (2) || PN||S1forallN. Let x' denote the exact solution of equation (2.3). A projection method to approxi- mate x' proceeds as follows: N Choose an approximation X” = Z I,,f,, e XN, and determine its coefficients [1,} n=1 by requiring that T xN — y vanishes under the projection PN. Let ERR = T xN — y. From the definition of PN, PN ERR = 0 if and only if =O, (m=1,...,N). Therefore the coefficients of xN are determined by =, (m=1,...,N), 15 that is ZN I = V , (2.10) where Z” isaanNmatrix with < 'r r,- , fj > as its (1,1)“ entry, I=(I,, . . ., IN )T are coefficients of the approximation xN, and V = ( < y , f1 >,...,< y , fN > )T is a vector with < y , f,- > as the ith component. In order to analyze this process and give an error estimate of the approximate solu- tion, we need the following theorem. Theorem 2.1. PN x —) x as N —> no for all x e L2(S). In particular, IIP~x—x||=0(1). ifxeczm. (2.11) where C 2(S ) = {x : S —) C2, x has continuous second derivatives. } Proof. Using the property (1) of P”, for any y e XN, =— = o , that is, PN x - x is perpendicular to the subspace XN. Therefore lanx-XI|=dis(X.X~). (2.12) where dis( x , XN )=mi)r{1|| x — X” II. It is sufficient to show that XNE N dis(x,XN)—>O aSN—)°°. 16 Case 1, xe C2(S)(x¢0): N Choose X” = 2 I,- f,-, where i=1 1: =X(I‘i) ' ni . (2.13) r,- is an arbitrary point on the ith edge and n,- is a unit vector normal to the ith edge. Let the error function e(r) be defined as e(r) = W) - xN(r) . Inside a triangle T, we draw a small triangle 7" with sides parallel to the sides of T and with vertices F,- ( i = 1, 2, 3). 8 is the uniform width between T and T' which satisfies 2mm:- 0<5< —-"——— (2.14) 3“ XII.» and the Lebesgue measure of S - UgT’i satisfies meas(s - U112) <12 . (2.15) Here II x II... = magtl x(r) I. Without loss of generality, T is labeled as shown in Fig.6. re Atr=F1, 3 901012 = x(i"1) ' n2 - Z 1i I’i(|‘:1) ' n2 is] - 11 12 13 =X(l’1)'"2- [13:514’12 E47012 +52)+13 E53 =x(i‘,)-n2—12—8, 9 4'; n; n 3 I‘12 lg . T, 13 F3 i”'2 I 1 Fig.6. Triangle T and T’. 3 wherel8,-|=8,AistheareaoftriangleT, and 8:33-21,- I-,-l8,-.From(2.13) e(i‘1) ‘ u2 = x(i"1) ' n2 - X(l’2) ' n2 - 5 =IVx(r2+tz)zdt-n2-8, WhCYCZ=F1-l'2. 0 Using (2.14), yields |e(F1)-n2| Sll Vx||.. 4+1 = const - l , (2.16) 18 where const = (II Vx||.. +1 ) is independent of t. For simplicity, we shall use cons: to denote constant numbers that are independent of I. Here and throughout, const’s may not have the same value. Similarly, we get I e(F1)-n3| Scans: - l . (2.17) Let i,j be an orthonormal basis on S. Then there exist real numbers c1 and (:2, such that 1:61 112 +62 I13 . Applying n2, n3 on both sides of this equation, we have i-n2=ct+n3-n2cz, i-n3=n2-n3 €1+Cz . The above two equations can be written as a matrix equation 1 "2'03 Cr i'nz nz-ng 1 (:2 i - 113 Let B be the angle between 112 and 113. From assumption 3) and Fig.6, (n2 -n3 )2 = 005213 = c0320, 5 005260. Thus the determinant of the above coefficient matrix 1 n2'113 det =1---(“2'“3)2 n2°n3 1 21—cos260>0. 19 It follows from Cramer’s Rule that |c,-|Sconst i=1,2. Here, cons: is independent of triangular elements. Therefore |e(F,)-i|=| c1e(F,)-n2+c2e(F1)-n3|Sconst-l by (2.16) and (2.17). Similarly, we obtain |e(F1)-j| Sconst-l. Thus I C(F1)I=VIC(F1)'i|2+I C(Fl)'j|2 SCOIISI '1. Repeating this process, we get |e(i",-)| Scans: - l , i=1, 2, 3, (2.18) where each cons: depends only on H Vx H. and is proportional to 00- Next we are going to show that (2.18) holds for all points on the edges of 7": PM = {fl T=Fg +£(Fj -F;), I E [0, 11]}, (i, j) E {(1, 2), (1, 3), (2, 3)}. Once this is proved, (2.18) will be true for all the points of the segments which have end points on those edges. Therefore |e(r)| Sconst-l, forallre 7,, i=1,2,..., n1, (2.19) here n1 is the number of triangles. We shall see that cons: here depends only on IID“XII..(I al 52)andCo. 20 It is sufficient to prove (2.18) on edge 1‘12. On 1",; , there is a point r0 = r", + to(i'-2 — F1) (to e [0, 1]), such that I e(ro)| = maxl e(r)| . ref}, Without loss of generality, we assume that I e(ro) I > 0. We let = e(ro) no |e(I'o)| and g(t)=e(F1+t(F2—i"1))-no . 1610.1]- Since g e C 2[0, 1] and g(to) is an extrema, g(0) = g(to) + 332$) (—t0)2 , for some § 6 [0, 1] . Consequently |e(ro)l=lg(to)15|8(0)l+%i gv<§>l .<_| e(F1)-no| +const 'M -| F2 -i"1|2 Sconst-l-+-const'M'l2 , where M = ”mags; D“x(r)| . Thus |e(r)| Sconst-l, forallre [‘13. It is easy to see that (2.19) can be proved by using the same technique. Now, Hence, 21 || e||2=I| e(r)!2 a.- S = I |e(r)|2dr-t- I |e(r)|2dr Ur, S—UT', Swim-[2+2 I |x(r)|2dr+ I |x,,,(r)|2 dr by (2.19) S-Ur, S‘Uri s. const - 12 +2 n xuz. +( 3” it". max" 11".. )2 meas(S —Ur.-) I 5 const - 12 by (2.4) and (2.15). dis(x, XN)S|| e|| Sconst'l, where const depends only on H Bax ll... (I al 5 2) and is proportional to Co- Case 2, x e L2(S): Since C2(S) is dense in L2(S), there exists a sequence [ x9") } in C2(S), such that x‘“’ —> x in L2 norm as m —->oo. For any 8 > 0, choose x‘") satisfying || x-x‘") H < 812. Then dis(x , XN)S|| x-x("') || +dis(x("') , XN) OasN—)oo. El CHAPTER 3 Convergence and Error Bound In this chapter, we show that the projection method defined by a subspace XN and an operator PN determines a unique solution of the linear system (2.10), which then gives an approximate solution X” of the integral equation (2.1). The sequence { xN } converges to the exact solution x. of (2.1) as the diameters of elements tend to zero uni- formly. Moreover, the error bound can be found and the method gives a first order approximation if x' is smooth enough. Let the operator K and projections [PM be defined as in chapter 2. We denote HN =E +PN K and RN = HNI x, . Here HN is the restriction of fly on subspace X”. Then equation (2.8) can be written as min, = PNy . (3.1) This equation is uniquely solvable for any y e L2(S) if and only if HIT/l exists. Since {RN} converges to T, which is invertible and its inverse is bounded, the existence of HIT/1 and HE} is guaranteed by the following lemma and theorem: Lemma3.l. II PNK-Kll —-)0 as N —>oo. Proof. The assertion follows from theorem 2.1 and lemma 2 in [2, pg.53]. CI Theorem 3.2. Let X be a Banach space and L (X) the normed vector space of bounded linear operators from X into X. Suppose L, L0, L61 6 L (X). If 23 24 A=||L—Lo||||L51|| <1,thenL-1e L(X) and L"=[ 2(L5‘1Lo-L»"1La‘. n=0 Proof. See [3. 138.87]. CI Lemma 3.1 enables us to choose N 0, such that q~=llP~K-KIIII'1“‘||<1. ifN>No- (3.2) Using theorem 3.2, H17} exists when N > N o and Hi1=i i (1"1 (K-P~K»"1T‘. n=0 Then we have the inequality II HIT/1 II s 1 || T1”, whenN>No. (3.3) l—qN Since equation (3.1) is uniquely solvable for any y e L2(S) if and only if the linear sys- tem (2.10) is uniquely solvable, the matrix ZN is nonsingular if N > N 0- This proves the feasibility of the approximation procedure. In order to prove the convergence and give an error analysis, we apply PN on both sides of (2.3): rN(E+K)x‘=PNy. Using (2.8), we find that HN(x’ —x~)=x' —PNx‘. From(3.3), we get 25 H X' - XN II 5 II “N1 II H X‘ - PNX’ || 1 s l-qN u T" || n x‘ —PNx‘ || ifN >No. Applying theorem 2.1 and lemma 3.1, yields II X*—XNII -)0, (LYN—900, where X” is the unique solution of (3.1) and its coefficients are determined by the solu- tion of matrix equation (2.10). Moreover, from theorem 2.1 we have the error bound II x. —xNII Scans: ° I , ifx’ e C2(S). Here the constant number cans: depends only on H Dax' II... (I Gil S 2) and is propor- tional to co defined in chapter 2. CHAPTER 4 Computational Stability The approximated solution X” is determined in terms of the basis functions {f,,} of subspace XN used in the computational scheme. To examine the influence of these basis functions on the stability of the linear system ZN I = V , (4.1) we have the following theorem: Theorem 4.1. The condition number of the matrix ZN of (4.1) satisfies the inequality cand( ZN ) S cons: - cond ( GN ) when N is large enough. (4.2) Here cans: is a constant number independent of N, and GM is an N XN matrix defined in assumption 4) of chapter 2. Proof. Let C be the set of all complex numbers, and (D : XN -9 C” be defined by ,...,)T, forallxe X”, N i.e., if X” = Z I,- f,- e XN, then i=1 I2 i=1 N 52(llftllz I Ix~|2dr) i=1 TiUTT N 2 smalls-"2x2: i Ix~l dr) ' i=17‘1‘u17 s(x~) = GN I. (4.5) Using assumption 4) of chapter 2, we obtain IND-1“)“2 =11 xN”2 N 2 =Il 2 lift“ '=1 =I‘GNI = ( Gi’i‘c )‘ G~ ( Gn‘c) by (4.5). That is II We)“2 =c‘ GIT; c 1 n _. =l—2c 1"” c sin FN‘IIIICIIZ 12 1 1 2 S—— c b (2.6) :2 50 || H y cans: =7" c”2- Here cans:( = —1- ) is independent of N. Thus 50 29 II (VI II S cans: /l . (4.6) Next, applying (D to (3.1) and noting that ZNI=¢(PN Y). we find Z" r = a ti” xN =’1GN I) foraner C". This implies that zN=fiN o4 G”. Thus we obtain cont“ ZN )=|| ZN || °|| Z11/1" S(|| 4’“ || “71 || )2“ “N" H H1711 ll c0nd(GN) Scanst°cand(GN) ifN >N0 , where N o is the positive integer used in the inequality (3.2). [3 Using assumption 4) of chapter 2, we have cand(GN)S80/80 forallN. Therefore the numerical scheme is stable as N —> oo. The bound (4.2) suggests that we should design triangles with considerable caution to make cond( GN) as small as possible. Also, the modulus of operators E + K and 30 (E + K )'1 will influence the stability of the scheme. With small modulus, we may expect more accurate approximation. CHAPTER 5 Numerical Results In this chapter we demonstrate the formulation of coupled surface integral equations(CSIE) and corresponding matrix equation. Then the numerical results are presented for the case of a homogeneous dielectric sphere in free space and excited by a plane wave. 5.1 The formulation of coupled surface integral equations The detailed derivation of the coupled surface integral equations can be found in [4] and [5], but for completeness and further numerical development only a summary of the CSIE equations is given below. Referring to Fig.7, S denotes the surface of a homo- geneous, lossy dielectric scatterer having a volume V contained in region 2 and bounded by the surface S. The scatterer is located in region 1 representing an isotropic, lossless free space medium. Let ( E’ , Hi )= scattered electric and magnetic fields in region 1 ( E’ , H5 )= scattered electric and magnetic fields in region 2 . Then, referring to the electromagnetic equivalent principle, various scattered electric 31 32 E‘ 11 2 REGION 1 k (81,111.01) (E2. Hz) (82 9'12 ’62) "8 Surface S / ’ > 0 Volume Y V REGION 2 / X / J (El 1 HI) M Homogeneous Lossy Dielectric Scatterer Fig.7. Geometry of a homogeneous lossy dielectric scatterer in an isotropic free space medium. and magnetic fields in regions 1 and 2 are given by 33 {(r) = -j(nA1(r)- vv,(r) — imer) (5.1.1) {(r) = -ij1 (r) - VU1(r) + ItLVXAlm’ for r on or outside S (5.1.2) 1 E§(r) = jcoA2(r) + VV2(r) + EI—erqm (5.1.3) 2’ 5(r) = jorF2(r) + VU2(r) - Il-VXAN'). for r on or inside S (5.1.4) 2 where the various vector potentials A,- and F,- and the scalar potentials V,, and U,, for i=1,2aregivenby A.(r> = f—n II J(rt)G.-(r. r.) dS(r') (5.1.5) mm = i,- ISI M(rI)G,-(r, rr) dS(rr) (5.1.6) mm = 4—;— II p‘(rt)G.-(r. r.) (ism) (5.1.7) mm = 4—1”; I] p"(r')G.-(r. r') dS(r') (5.1.8) e,r=e, [14%] (5.1.9) and ‘ I --—l I' I P(r)-J.m (V. 10)] (5.1.10) 34 m I -:—1- I- I P (r)- 1.0) [Vs M(r)1- (5.1.11) In obtaining the above expressions, 81“" time dependence is assumed for various field quantities and or is the fiequency in radians per second. The Green’s function defined in (5.1.5)-(5.1.8) fori = 1, 2 is given by elk-R 010'. r') = (5.1.12) R: l r-rvl (5.1.13) and the propagation constant is k,- = [ 8211,1531” . (5.1.14) In (5.1.5) and (5.1.6), J is the equivalent electric current and M is the equivalent mag- netic current on the surface of the dielectric scatterer. The equivalent electric and mag- netic currents are, in fact, related to the surface total magnetic and electric fields tangen- tial to the surface S: J(rr) = n,xH(rr) (5.1.15) M(rr) = E(rr)xn,, rr on the surface S (5.1.16) where n, is an outward unit normal on S shown in Fig.7. Further, in the above expres- sions, (81, 111. 0, =0) and (82, 112. 02) are the permittivity, permeability, and conductivity for the regions 1 and 2. On enforcing the boundary condition that the total tangential electric field and the total tangential magnetic field should be continuous across the surface of the arbitrary dielectric scatterer, the following coupled surface integral equations are obtained in terms of the unknown surface equivalent electric and magnetic currents: 35 5‘81... = 1 1°01 A1(l‘) + A209 1+ 1 Win-1+ “’20) 1 61’ 82’ +Vx[ Fm + F2“) I )|,,,, (5.1.17) H‘s]... = 1 jth. (r) + F2011 +1 VU1(|') + VU2(r) 1 mm + A2(r) ill #2 —Vx Illa... ronsurfaceS (5.1.18) where E‘ and Hi are the incident electric and magnetic fields in region 1 and the sub- script "tan" refers to tangential component. Thus, coupled surface integral equations in general can be expressed as TnJ(r) + r,,M(r) = E‘(r) (5.1.19) r,,.1(r) + rumor) = H‘ (r) (5. 1.20) for r e S, where T,- are linear integral operators, J, M are equivalent electric and mag- netic surface current, E‘, H‘ are incident electric and magnetic fields respectively, and S is a regular two dimensional region of interest. Set T11 T12 T21 T22 3(r) = [ J(r) . mo 1’ fi(r)=15‘(r) . Mr) 1’ . Equations (5.1.19) and (5.1.20) then become 36 Tj(r)=E(r) re S. Under appropriate assumptions, the theoretical analysis demonstrated in the previous chapters can be extended to certain types of coupled surface integral equations. 5.2 Matrix equation Referring to the dielectric scatterer shown in Fig.1, the surface electric and mag- netic current, J and M, distributions are expanded in terms of vector basis functions defined in chapter 1. Let N represents the total number of edges. Then the testing func- tions are written as N J(rt) = 2 1,1,(1-2) (5.2.1) as] M(rr) = ilMJArI) (5.2.2) where 1,, and M, are coefficients yet to be determined. Since the normal component of f, at the nth common edge connecting T; and T; is unity, each coefficient of 1,, and M. can be interpreted as the normal components of the electric and magnetic current density flowing past the nth common edge. In order to find the current coefficients, the coupled integral equations (5.1.19) and (5.1.20) are tested with respect to testing functions. The electric current and magnetic current expansion terms defined in (5.2.1) and (5.2.2) are now substituted into (5.1.19)- (5.1.20) and test the equations based on L2 inner product. Hence, we obtain the matrix equation 37 +== (5.2.3) += (5.2.4) for m = 1... N. For the detailed numerical evaluation of matrix elements, the reader may refer to [6]. 5.3 Numerical results Although the coupled surface integral equations shown in section 5.1 are not second kind Fredholm equations as we assumed in previous chapters, the numerical / Fig.8. Triangular surface patching for top half of sphere (plane view). 38 Table 1 Mi) 9 Numerical Solution Exact Solution 0.261799 0.771489 0.867912 0.785398 0.625174 0.710197 1.308997 0.455707 0.503564 1.832596 0.426590 0.461907 2.356194 0.530623 0.578552 2.879793 0.613626 0.670465 Error(e“)) 6.53093e-02 No. of faces = 60; No. of edges = 90; Matrix size = 180x180; Mesh size(l“)) = 2112/6. results show that the convergence rate is very close to our estimation. We shall present the experimental results for the case of a homogeneous dielectric sphere in free space and excited by a plane wave. As we know, the sphere is not amenable to planar triangular modeling. Assuming that the geometrical discretization error can be decreased by reducing the size of the ele- ments and such misfits are ignorable if the size of the elements is small enough, the sur- face of the sphere is modeled in terms of triangles having arbitrary edges and vertices arranged to depict the shape of a sphere. Fig.8 shows the plan view (top hemisphere 39 Table 2 M1 9 Numerical Solution Exact Solution 0.196350 0.818644 0.877422 0.589049 0.731751 0.783832 0.981748 0.584576 0.628048 1.374447 0.458582 0.485250 1.767146 0.433632 0.454231 2.159845 0.504070 0.530279 2.552544 0.592780 0.621703 2.945243 0.639877 0.675944 Error(e (2)) 3.87571e-02 No. of faces = 112; No. of edges = 168; Matrix size = 336x336; Mesh size(lm) = 211/8. only) of the triangular scheme adopted. In our experiments, the electrical size of the sphere is kla =1 where the free space propagation constant k, = 211/10 and the radius of the sphere is a. The relative dielectric constant of the sphere is e, = 4. The sphere is located in free space and is excited by an axial incident plane wave. Table 1 and Table 2 show the numerical results of the induced magnetic current on the surface of the sphere along a circumferential are in xz plane. Along the are, there are two components of the 40 Table 3 Error(e (2) / 8‘”) 0.593439 Ratio Mesh Size(l (2) / l (1)) 0.750000 magnetic current M g and M 4,. The results of the induced magnetic current are shown nor- malized with respect to incident electric field. The ratio of the errors (e 2/e1) and the corresponding ratio of mesh sizes are shown in Table 3. The experiments show that the convergence rate is at least linear. It is a lit- tle better than our estimation. This may be caused by the geometrical discretization error, the numerical integration of matrix elements and linear system solver. We have not taken into account these errors in our analysis. Since the basis functions are linear inside the patches and discontinuous on the edges, it is difficult to expect that the Galer- kin method on such basis functions can achieve higher a convergence rate than linear convergence. SUMMARY The Galean method with new vector basis functions has been used by researchers in electrical engineering. Many experiments show that the new vector basis functions are ideally suited for representing surface electric current J and surface magnetic current M on the triangulated surface of the given dielectric scatterer. Efficient and simple numeri- cal algorithms have been developed. In this work, a theoretical analysis on the conver- gence, error estimate and numerical stability of this method is given for certain type of electric field integral equations. More specifically, we are concerned with Fredholm integral equations of the second kind: (E+K)J(r)=E"(r). J. E‘ e L2(S). where E is the identity operator and K is a completely continuous integral operator. This restriction is sufficient but not necessary. If the integral equation T J (r) = Ei (r) is well posed on L2(S) or on a subspace of L2(S) which contains U XN, that is, the integral N=1 operator T is invertible and its inverse operator is bounded on the image of T, and in addition PN( T — E) -) T - E as N —> oo, then all the results of this work still hold. In practice the physical properties of electric field integral equations may provide the first condition, while the second condition is not easy or practically is impossible to verify. In chapter 5, The numerical results of the coupled surface integral equations on a homo- 41 42 geneous dielectric sphere suggest that the convergence rate is at least linear, which confirms our estimation. The analysis of the Galean method with vector basis functions for general integral equations. especially for integral equations of the first kind remains an open question. BIBLIOGRAPHY BIBLIOGRAPHY [1] [2] [3] [4] [5] [6] A. W. GLISSON, 0n the Development of Numerical Techniques for Treating Arbitrarily-Shaped Surfaces. Ph.D. dissertation, Univ. of Mississippi, 1978 K. E. ATKINSON, A Survey of Numerical Method for the Solution of Fredholm Integral Equations of the Second Kind. Society for Industrial and Applied Mathematics, 1976. V. HUTSON AND J. S. PYM, Applications of Functional Analysis and Operator Theory. Academic Press Inc. (London) LTD. 1980. K. R. UMASHANKAR AND A. TAFLOVE, Analytical Models for Elec- tromagnetic Scattering. Part-I, Final Rep., Contract F19628-82-C-0140 to RADC/ESD, Hanscom Air Force Base, MA, June 1984. l. R. MAUTZ AND R. F. HARRINGTON, Electromagnetic Scattering From a Homogeneous Material Body of Revolution. Arch. Elek. Ubertragung, vol. 33, no.4, pp.7l-80, Apr. 1979. K. R. UMASHANKAR, A. TAFLOVE AND S. M. RAO, Electromagnetic Scattering by Arbitrary Shaped Three-Dimensional Homogeneous Lossy Dielectric Objects. IEEE Transactions on Antennas and Propagation, vol. AP-34, no.6, June 1986. 43 I MICHIGAN STATE U 111 1111111111111111111115 31293007885852