Symmetric Galerkin boundary integral formulation for interface and multi-zone problems

INTERNATIONAL JOURNAL FOR NUMERICAL METHODS IN ENGINEERING, VOL. 40, 3085–3101 (1997) SYMMETRIC GALERKIN BOUNDARY INTEGRAL FORMULATION FOR INTERFACE AND MULTI-ZONE PROBLEMS L. J. GRAY1∗ AND GLAUCIO H. PAULINO2 1 Computer Science and Mathematics Division; Oak Ridge National Laboratory; P.O. Box 2008; Bldg. 6012; MS-6367; Oak Ridge; TN 37831-6367; U.S.A. 2 Department of Civil and Environmental Engineering; University of California; Davis; CA 95616-5294; U.S.A. ABSTRACT Domains containing an ‘internal boundary’, such as a bi-material interface, arise in many applications, e.g. composite materials and geophysical simulations. This paper presents a symmetric Galerkin boundary integral method for this important class of problems. In this situation, the physical quantities are known to satisfy continuity conditions across the interface, but no boundary conditions are speci ed. The algorithm described herein achieves a symmetric matrix of reduced size. Moreover, the symmetry can also be invoked to lessen the numerical work involved in constructing the system of equations, and thus the method is computationally very ecient. A prototype numerical example, with several variations in the boundary conditions and material properties, is employed to validate the formulation and corresponding numerical procedure. The boundary element results are compared with analytical solutions and with numerical results obtained with the nite element method. ? 1997 by John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng., 40, 3085–3101 (1997) No. of Figures: 15. No. of Tables: 5. No. of References: 28. KEY WORDS: boundary integral method; symmetric-Galerkin approximation; interfaces; multi-zone boundary element analysis; composite materials 1. INTRODUCTION The interface problems of the title refers to any domain containing an ‘internal boundary’, one for which no boundary data is speci ed. Typical geometries are illustrated by Figure 1. These problems arise in many important applications, with composite materials1 being an area of much current interest. A particular example is provided by the analysis of composite beams, as shown in Figure 1(b). In this gure, the composite consists, for example, of bre reinforced layers joined at the interfaces. As the bre orientation in the layers is not the same, the material properties are Correspondence to: L. J. Gray, Oak Ridge National Laboratory, P.O. Box 2008, Bldg. 6012, MS-6367, Oak Ridge, TN 37831-6367, U.S.A. ∗ Contract grant sponsor: U.S. Department of Energy; Contract grant number: DE-AC∅5-96OR22464 Contract grant sponsor: Lockheed Martin Energy Research Corp. This article is a US Government work and, as such, is in the public domain in the United States of America CCC 0029–5981/97/163085–17$17.50 ? 1997 by John Wiley & Sons, Ltd. Received 29 July 1996 Revised 26 November 1996 3086 L. J. GRAY AND G. H. PAULINO Figure 1. Examples of interface problems: (a) Multilayered elastic rock mass; Ei and i , i = 1; : : : ; 4 are Young’s modulus and Poisson’s ratio of each layer. (b) Cantilever ‘sandwich’ beam with a concentrated load at one endpoint. (c) Heat transfer in a cross-section of a hollow shaft consisting of two concentric materials; T0 and T1 denote temperatures in the interior and exterior of the shaft, respectively Figure 2. Division of a thin geometry into zones for multi-zone boundary element analysis di erent. Assuming perfect bonding, the constraints on the interface are continuity of displacements and equilibrium conditions (equal and opposite tractions), but there are no prescribed values. In general, boundary value problems with interfaces exhibit relatively large surface to volume ratio, and in this case the boundary element method (BEM) is expensive compared to domain techniques such as the nite element method (FEM). Nevertheless, the BEM approach is very attractive for this class of problems, the overriding reason being that it provides a very natural treatment of the interface. As noted above, the physical conditions at the interface are usually continuity of the principal function and its derivative, e.g. displacement and traction in elasticity, or potential and ux in potential theory. For the displacement-based FEM, enforcing the continuity of traction (or ux) is a very dicult task. For boundary integral equations however, the derivative quantity appears directly in the formulation, and thus the interface conditions can be incorporated simply and accurately. This issue is quite important in practical applications. For instance, for a bi-material interface crack problem, the interface lies in the critical region ahead of the crack tip where an accurate solution is imperative. For the boundary integral method, there is another important class of problems for which an ecient interface algorithm is required, namely for long thin domains.2 These problems can be solved by means of a multi-zone boundary element analysis. Although adding surface area is generally a bad idea, in this situation it is computationally ecient to decompose the domain into small subdomains by means of arti cial interfaces (Figure 2). The reason is that constructing the boundary integral equations for each subdomain only requires integrating over a small piece of the boundary, and this more than compensates for the presence of the additional boundary surface. Note too that the resulting coecient matrix has a block-sparse structure, adding to the computational eciency. Int. J. Numer. Meth. Engng., 40, 3085–3101 (1997) ? 1997 by John Wiley & Sons, Ltd. SYMMETRIC GALERKIN BOUNDARY INTEGRAL FORMULATION 3087 The symmetric Galerkin (SG) method is a highly robust and ecient boundary integral approximation. It was introduced by Hartmann et al.3 in 1985, and further developed by Maier and co-workers4 and many others. A key advantage of the Galerkin formulation is the ability to treat hypersingular equations with standard continuous elements (e.g. linear, quadratic). The more commonly used collocation approximation5 requires a di erentiable boundary interpolation,6; 7 which is inherently dicult and computationally expensive (e.g. Hermite8; 9 or Overhauser10; 11 elements). Considering the essential role of hypersingular equations in the ecient solution of fracture and other problems12; 6; 13 this is a signi cant advantage of Galerkin-type formulations. Compared to collocation methods with standard continuous elements, Galerkin is also quite expensive, but the added symmetric aspect makes this approach computationally ecient. Although non-symmetric Galerkin is roughly an order of magnitude slower than collocation, SG is as fast or faster.14; 15 As a consequence, SG is an ideal algorithm for crack geometries.16; 17 In addition, the FEM will produce a symmetric coecient matrix, and thus SG is especially advantageous for coupling boundary integral analysis with nite elements.18 As will be explained in more detail below, for standard boundary value problems, a symmetric matrix is obtained by suitably combining the standard boundary integral equation (BIE) and its derivative, the hypersingular boundary integral equation (HBIE). The choice of equation enforced at a node depends upon the boundary condition at this point. For the Dirichlet-type surface (i.e. speci ed potential, displacement, etc.) the BIE is employed, whereas the HBIE is written on the Neumann-type surface (i.e. speci ed ux, traction, etc.). For the general class of interface problems, however, no data is supplied on the ‘internal boundaries’ (see Figures 1 and 2). A natural question therefore arises—is it possible to develop a fully symmetric Galerkin formulation for interface problems? This paper will demonstrate that the SG approach does have a natural extension to interface geometries. Moreover, this formulation is computationally very ecient: the coecient matrix is of reduced size, solving for values only on one side of the interface, and the symmetry can be invoked to lessen the computational work required to set up the equations. 2. SYMMETRIC GALERKIN The boundary integral equations in this work are considered in a limit to the boundary sense.19; 20 This approach allows a uni ed and direct treatment of the singular integrals which appear in both the BIE and the HBIE. Moreover, this technique permits writing the same equation for points either inside the domain or on the boundary, bypassing potential diculties with the geometry dependent coecient of the ‘free term’ outside the integral.20; 21 The discussion below will consider the Laplace equation for the potential , ∇2  = 0, but it will be clear that the method is generally valid. This choice is not only for simplicity of notation, but also because an issue (to be discussed below) involving the material parameters arises in the potential formulation that does not appear in elasticity. This will require a slight reformulation of the BIEs for potential problems. The singular BIE can be written in the form5 Z Z @G @ (1) (P) + (Q) (P; Q) dQ = G(P; Q) (Q) dQ @n @n where n = n(Q) is the unit normal at a point Q on the domain boundary , and @(·)=@n denotes the normal derivative with respect to Q. The fundamental solution G(P; Q) depends upon the dimension of the problem, and is not unique. For the following discussion, we need to remark ? 1997 by John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng., 40, 3085–3101 (1997) 3088 L. J. GRAY AND G. H. PAULINO Figure 3. Bi-material interface problem that this function, and its derivatives, satisfy certain symmetry properties, i.e. G(P; Q) = G(Q; P) @G(P; Q)=@n = − @G(Q; P)=@n @2 G(P; Q)=@N@n = @2 G(Q; P)=@N@n (2) @G(P; Q)=@N = − @G(Q; P)=@N @G(P; Q)=@N = − @G(P; Q)=@n where @(·)=@N indicates the normal derivative with respect to P. The corresponding hypersingular equation is obtained by di erentiating equation (1) with respect to P in the direction N = n(P), the normal to the boundary at P. This results in Z Z @ @2 G @ @G (P) + (Q) (P; Q) dQ = (P; Q) (Q) dQ (3) @N @N@n @N @n The above HBIE is an equation for the normal derivative, but it is the ux F(P) ≡  @ (P) @N (4) which satis es the continuity conditions across a bi-material interface. The constant , which goes by various names (e.g. thermal conductivity) depending upon the application (e.g. heat transfer, electrostatics, potential ow), is di erent in each material. The continuity equations on the interface are then A = B (5) and A @B @A = −B @n @n (6) where the subscripts A and B refer to the top and bottom layer in Figure 3. By contrast, the material properties in an elasticity formulation are already embedded in the fundamental solution, and the hypersingular equation is most often written directly for the traction.12; 17 Thus, rather than equation (6), the interface equation would simply be cA = −cB , where c denotes the traction vector. Int. J. Numer. Meth. Engng., 40, 3085–3101 (1997) ? 1997 by John Wiley & Sons, Ltd. SYMMETRIC GALERKIN BOUNDARY INTEGRAL FORMULATION 3089 For the SG interface method, it will be convenient to rewrite the BIE, equation (1), and the HBIE, equation (3), in terms of the ux, i.e.  Z  Z @G 1 G(P; Q) F(Q) dQ (7) (P) + (Q) (P; Q) dQ = @n  and F(P) + Z   2 Z @G @G (P; Q) dQ = (P; Q) F(Q) dQ (Q)  @N@n @N (8) respectively. Note that including the constant  with G and @G 2 =@N@n does not alter the basic properties required for the symmetric Galerkin method, equations (2). These kernel functions are still symmetric with respect to P and Q, and the integrals containing a single derivative of G (which also play a role in achieving symmetry) remain unchanged. While this is a trivial rewriting of the boundary integral equations, it is key for directly embedding the interface continuity equations into the formulation and for obtaining a symmetric matrix. In a Galerkin approximation, equations (7) and (8) are satis ed in an averaged sense. Speci cally, the weighting functions are chosen to be the basis shape functions l (e.g. linear, quadratic) employed in the approximation of  and F on the boundary. Thus, Z Z Z @G (Q) (P; Q) dQ dP l (P) l (P)(P) dP + @n  Z  Z 1 G(P; Q) F(Q) dQ dP (9) = l (P)  and   2 Z Z Z @G (P; Q) dQ dP (Q)  l (P) l (P)F(P) dP + @N@n Z Z @G = (P; Q) F(Q) dQ dP (10) l (P) @N The additional boundary integration, with respect to P, is the last ingredient required to obtain a symmetric matrix. After discretization, the set of equations in matrix form can be written as [H ]{} = [G]{F}, and in block-matrix form these equations become #( ) " #( ) " bv G11 G12 Fu H11 H12 = (11) H21 H22 u G21 G22 Fbv The rst row represents the BIE written on the Dirichlet surface, and the second represents the HBIE on the Neumann surface. Similarly, the rst and second columns arise from integrating over Dirichlet and Neumann surfaces. The subscripts in the vectors therefore denote known boundary values (bv) and unknown (u) quantities. Rearranging equation (11) into the form [A]{x} = {b}, and multiplying the hypersingular equations by −1, one obtains ) #( ) ( " H12 −H11 bv + G12 Fbv Fu −G11 (12) = H21 bv − G22 Fbv u G21 −H22 T T , H22 = H22 , and H12 = G21 , now follows from The symmetry of the coecient matrix, G11 = G11 the properties of the kernel functions (equations (2)). ? 1997 by John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng., 40, 3085–3101 (1997) 3090 L. J. GRAY AND G. H. PAULINO 3. INTERFACE AND SYMMETRY To describe the interface algorithm, it suces to deal with a geometry having a single common boundary, as in Figure 3. The extension to more complicated interface geometries, if not immediately obvious, is nevertheless seen to be relatively straightforward after a little thought. This topic will be further discussed in Section 4. For convenience, we will refer to the two subdomains A and B as the top and bottom regions, respectively, as shown in Figure 3. The basic idea is to write the usual SG equations on the non-interface boundaries in each subdomain, together with an appropriate combination of equations on the common boundary, as indicated by Table I. In addition, only one set of variables is employed on the interface, i.e. potential and ux on the top side. Whenever the interface ux on the bottom side appears in an equation, it is replaced by the negative of the top ux. It will again be convenient to use block-matrix notation, and only the coecient matrix (the analogue of the left-hand side in equation (12)) will be considered. For the interface problem, the matrix is partitioned into a 4 × 4 block structure and takes the form   0 SAFI SAI SAA    0 SBB −SBFI SBI    (13)  SF A SF B SF F SF   I I I I I   I SI A SI B SI FI SI I The rst row corresponds to the SG equations written for the non-interface boundary in the top (A) domain, equation (9) or equation (10), depending upon the boundary conditions. The rst column therefore indicates integrations involving unknown quantities on the top (non-interface) boundary. The second row and column (B) are the analogous entries for the bottom material. Thus, SAA and SBB are (square) symmetric matrices, a consequence of the SG procedure. The o -diagonal (1; 2) and (2; 1) blocks are, as shown, equal to zero, as the top equations do not involve the bottom geometry, and vice versa. The third and fourth columns refer to the unknown interface ux (FI ) and interface potential (I ), respectively. The top and bottom equations (i.e. the rst and second rows) integrate over the interface and thus the (1; 3), (1; 4), (2; 3) and (2; 4), entries are in general non-zero. At the risk of being repetitive, these rows do not include equations for the interface, i.e. equations (9) and (10) with the weighting functions l centred on an interface node. The minus sign multiplying SBFI is due to the change in sign for the interface ux, as mentioned above. The key to achieving a symmetric matrix is to ll in rows three and four with appropriate interface equations (Galerkin weighting functions l centred on interface nodes) for determining I and FI . The successful procedure is as follows (see Table I). Table I. Symmetric Galerkin integral equations for an interface or multi-zone problem Surface Dirichlet Neumann Interface (FI ) (I ) Int. J. Numer. Meth. Engng., 40, 3085–3101 (1997) Integral equations BIE HBIE BIE(top) – BIE(bottom) HBIE(top) + HBIE(bottom) ? 1997 by John Wiley & Sons, Ltd. SYMMETRIC GALERKIN BOUNDARY INTEGRAL FORMULATION 3091 For the interface ux, row three, write the BIE (9) for the top domain, minus this equation constructed for the bottom domain. This is reasonable, as equations for ux in the usual SG procedure indicate that potential is speci ed, and therefore the BIE is employed. Thus, it follows T , the T superscript indicating transpose, as this is the SG procedure for the top that SFI A = SAF I domain. One can convince oneself of this simply by assuming that the top boundary is either all Dirichlet or all Neumann, and noting which integration kernel comes into play in forming this matrix. Symmetry also holds for the SFI B contribution, as taking the negative BIE over the bottom domain compensates for the negative sign in the (2; 3) position. It is required that the diagonal block SFI FI be symmetric on its own, and this easily follows from the observation that this entry comes from integrating the symmetric kernel G(P; Q)= over the interface. For the interface potential, row four, write the HBIE (10) for the top domain, plus this equation constructed for the bottom domain. Again, as these equations are thought of as the equations to determine the interface potential, the use of the hypersingular equation is consistent with the usual T T follows again from the basic SG algorithm. The and SI B = SB SG procedure. That SI A = SA I I symmetry of the (4; 4) block follows immediately by noting that this matrix arises from integrating the symmetric kernel @2G(P; Q)=@N@n over the interface. It therefore remains to show that SFI I = STI FI . Note rst that SFI I originates from the left hand side of equation (9). However, as the single integral term is the same for both top and bottom equations, this term drops out, leaving just the @G=@n integration. On the other hand, SI FI comes from the ux terms in equation (10). The single integral term once again drops out, this time because of the change in sign of the ux across the interface, and thus both matrices come from the rst derivative kernel. Moreover, the di erence of the two equations used for the third row is once again matched by the change in sign in the ux, and so the symmetry follows. It is interesting to note that, for the Laplace equation and a at interface, SI FI = 0. The normal derivative of the Green’s function is given by 1 n·R @G =− @n 2 r 2 (14) where r = kRk = kQ − Pk. Due to the n · R factor, the only surviving term in integrating over a at surface is the singular contribution. However, as in the single integral term, the singular part is the same on both sides of the interface (note that the material constant  is not present), and therefore cancels. In some applications, this observation could be invoked to reduce the computational e ort. 4. REMARKS Remark 1 (Multiple Interfaces). As mentioned above, the extension of the interface algorithm to more complicated geometries is not dicult, but nevertheless warrants some further discussion. Consider the problem with three zones, labeled A, B and C, illustrated in Figure 4. The noninterface boundary segments for each subdomain are denoted by a; b; c, respectively, and d; e; f denote the interfaces. As before, the usual SG procedure is followed on the non-interface segments, and combinations of BIE and HBIE (see Table I) are used on the interfaces. Thus, for example, the equations written for the interface segment d are potential and ux equations for regions A and C, while the equations for e involve regions A and B. ? 1997 by John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng., 40, 3085–3101 (1997) 3092 L. J. GRAY AND G. H. PAULINO Figure 4. A more complicated interface problem The coecient matrix for this problem takes the form   Saa 0 0 Sad Sae 0    0 Sbb 0 0 Sbe Sbf       0 0 Scc Scd 0 Scf     Sda 0 Sdc Sdd Sde Sdf       Sea Seb 0 Sed See Sef    0 Sfb Sfc Sfd Sfe Sff (15) Here, to simplify notation, the unknown potential and ux on the interface segments have been combined. Thus, the fourth block column corresponds to potential and ux on the d interface. The symmetry of the upper-left 3 × 3 block again follows from the SG procedure. Veri cation of the symmetry for the rest of the matrix requires more thought, but follows along the same lines as the reference problem discussed in Section 3. Remark 2 (Corners). As illustrated by the above example (Figure 4), multiple interface con gurations can produce dicult corner problems. For three-dimensional problems, even more complicated corner/edge con gurations can occur. With respect to the A and B subregions of Figure 4, there are two unknown uxes to be determined at the corner point where all three subdomains meet. This situation is analogous to a boundary value problem in which Dirichlet data is supplied at a boundary corner. For a collocation approximation, this situation requires special treatment,22 as additional equations must be written to determine the unknown values. Another distinct advantage of the Galerkin formulation is that the extra equations are easily handled by proper de nition of the Galerkin weighting functions at the corner point. All that is required is appropriate use of ‘double nodes’ in the discretization process. Thus, the SG provides and e ective treatment of corners with conforming elements, which includes the type of corners produced by interface con gurations. Remark 3 (Free interface). Another special situation worth mentioning is illustrated in Figure 5. This con guration arises in applications such as the transmission problem in wave scattering,21 and in this terminology the scattering object is completely embedded in the host material. In this case the interface forms the entire boundary surface for the inclusion subregion. The SG algorithm in this case is straightforward, the boundary integral equations for the object consisting solely of Int. J. Numer. Meth. Engng., 40, 3085–3101 (1997) ? 1997 by John Wiley & Sons, Ltd. SYMMETRIC GALERKIN BOUNDARY INTEGRAL FORMULATION 3093 Figure 5. The interface boundary can comprise the entire boundary of a subdomain, as in this wave scattering / transmission problem interface equations. The coecient matrix will take the form of equation (13), with the rst row and column removed, and hence, is clearly symmetric. Remark 4 (Computational savings). As in the SG formulation for a non-interface problem, the primary computational advantage of the symmetric formulation is in the reduced time required to solve the system of equations. For a direct solution, this is of the order of N 3 =6 versus N 3 =3 for a non-symmetric matrix. Thus, for suciently large problems, this factor of two can translate into signi cantly less computational time.14 Some savings can be realized in constructing the matrix, but in general this is unfortunately limited: although the coecient matrix is symmetric, the contributions to the right-hand side vector must still be computed for each equation. Thus, the number of calculations that can be bypassed due to symmetry is, speaking somewhat loosely (as this can depend upon the equation being solved), not large relative to the total operations required. For interface problems, however, the symmetry of the coecient matrix can be used to good advantage. Note that in (13) the upper right 2 × 2 block need not be calculated, due to symmetry. These matrix elements come from the boundary integral equations written for the non-interface boundaries, integrating over the interface. As indicated above, the only reason for now integrating over the interface in these equations is for constructing the right-hand side vector. However, as there are no boundary conditions on the interface, there is no contribution to the right-hand side. Thus, in forming the equations for the non-interface boundary, integration over the interface is not required. Depending upon the geometry of the problem, this can result in considerable computational savings. For instance, for the prototype problem considered in the next section, invoking this symmetry, i.e. bypassing the interface integrations for the non-interface equations, produced roughly a 13 per cent reduction in total computation time—1.05 versus 1.21 s. For this relatively small problem (102 nodes), the matrix solution time is minimal, and thus the percent reduction (measured with total time) will be less when larger matrices are involved. Nevertheless, it is clear that a signi cant amount of computation has been avoided. The 13 per cent number is consistent with related symmetric Galerkin fracture calculations. In this case, integrations over the crack surface can be skipped when forming the non-crack equations.17 Remark 5 (Related work). At the nal stage of this work, we became aware of related work by Layton et al.23 They have used the symmetric Galerkin concept to develop a multi-zone BEM formulation which results in a partially symmetric coecient matrix. The symmetry is restricted to the non-interface degrees of freedom of each zone, and the degrees of freedom associated with the ? 1997 by John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng., 40, 3085–3101 (1997) 3094 L. J. GRAY AND G. H. PAULINO interfaces lead to non-symmetric subblocks. They demonstrate that condensing out the symmetric portion of the system and solving the resulting reduced non-symmetric system leads to an e ective algorithm. However, these authors note that a fully symmetric formulation, as presented herein, would be highly advantageous. For large-scale calculations, the condensation techniques discussed in Reference 23 would still prove useful when applied to the present fully symmetric system. In this case, the condensation process is similar to the substructuring techniques in the FEM.24 5. CALCULATIONS A prototype numerical example modelling heat transfer in an eccentric annulus geometry is considered. To validate the interface algorithm, several variations in the boundary conditions and material properties are examined. This problem has characteristics which make it very suitable for testing purposes, e.g. curved boundaries, a non-convex region, and an interesting (i.e. non-trivial) corner situation at the interface junctions. The basic interpolation approximation employed in these tests consists of standard isoparametric, conforming, linear elements. Figure 6 shows the geometry, boundary conditions, and material properties for the rst calculation. The thermal conductivities for the top and bottom parts are denoted by kA and kB , respectively. If only one material is present, i.e. kA = kB , then there is a closed-form solution for this problem. This solution can be obtained by means of conformal mapping and complex variable techniques, and is given by (see, for example, the book by Greenberg25 ) "
# √ ln u2 + v2 ; a=2+ 3 (16)  = 100 1 − 2 ln a where u(x; y) = (ax − 1)(x − a) + ay2 ; (ax − 1)2 + a2 y2 v(x; y) = (a2 − 1) y (a x − 1)2 + a2 y2 (17) Figure 6. Geometry, boundary conditions and material properties for heat conduction in an eccentric annulus with two types of material Int. J. Numer. Meth. Engng., 40, 3085–3101 (1997) ? 1997 by John Wiley & Sons, Ltd. SYMMETRIC GALERKIN BOUNDARY INTEGRAL FORMULATION 3095 To verify the computer code, the material constants kA and kB were set equal (kA = kB = 1·0), and the problem solved using interfaces as shown in Figure 6.∗ The boundary element mesh employed for the top-half of the domain is shown in Figure 7, the discretization for the bottom-half is identical. A comparison of the theoretical solution (equations (16) and (17)) with the numerical solution is provided in Tables II–V for some representative points along the boundary. The actual matching interface uxes are zero, and the SG-BEM values were of the order of 10−12 or less. These results indicate that the SG interface algorithm is correct. The more interesting situation, kA 6= kB , has also been investigated. For example, if kA = 1/2 and kB = 1·0, then the numerical solution on the bottom semicircles is essentially the same as before, but now the uxes on the top semicircles are half of what they were before. This happens because the actual matching ux across the interface is zero. The boundary element solution in this case was consistent with this observation, and thus provides additional validation of the present numerical solution procedure for interface problems. Figure 7. Boundary element mesh (58 nodes and 54 elements) considering the symmetric part of the eccentric annulus. The corresponding mesh for the interface calculations, obtained by symmetry with respect to the horizontal line (y = 0), has 116 nodes and 108 elements Table II. Flux at the top semi-circle with radius = 1·0 @=@n Angle (rad) 0 =4 =2 3=4  Theory SG-BEM Error (%) 131·5187 101·7244 65·7595 48·5829 43·8397 133·7729 101·2764 65·4230 48·4524 44·7766 −1·71 0·44 0·51 0·27 −2·14 Table III. Flux at the top semi-circle with radius = 0·25 @=@n Angle (rad) 0 =4 =2 3=4  Theory SG-BEM Error (%) −350·7162 −334·3935 −300·6147 −273·0345 −263·0383 −346·1688 −329·7187 −300·2498 −269·2404 −259·4839 1·30 1·40 0·12 1·40 1·35 Note that for the single material case, one could take advantage of symmetry and model only the top (or bottom) part of the problem. However, the goal here is to test the interface algorithm ∗ ? 1997 by John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng., 40, 3085–3101 (1997) 3096 L. J. GRAY AND G. H. PAULINO Table IV. Potential at the LHS interface (y = 0)  x −1·0000 −0·8333 −0·6667 −0·5000 −0·3333 −0·1667 −0·0833 0·0000 Theory SG-BEM Error (%) 100·0000 92·0218 82·3856 70·4045 54·8778 33·4085 18·8852 0·0000 100·0000 92·2483 82·6837 70·8167 55·5480 34·0673 19·4347 0·0000 0·00 −0·25 −0·37 −0·59 −1·22 −1·97 −2·91 0·00 Table V. Potential at the RHS interface (y = 0)  x 0·5000 0·5625 0·6250 0·7500 0·8750 1·0000 Theory SG-BEM Error (%) 0·0000 19·5925 35·7162 61·6269 82·3856 100·0000 0·0000 20·1178 36·3434 62·2651 82·8032 100·0000 0·00 −2·68 −1·76 −1·04 −0·51 0·00 Figure 8. Geometry, new boundary conditions and material properties for heat conduction in an eccentric annulus A new boundary value problem for the eccentric annulus con guration is shown in Figure 8. Since there is no analytical solution available for this problem, the SG-BEM results will be compared with those obtained by a FEM calculation. The nite element mesh for the problem of Figure 8 is shown in Figure 9. The graphs in Figures 10–15 compare the numerical results obtained by the two methods. Figures 10 and 11 compare the nodal temperature on the top Int. J. Numer. Meth. Engng., 40, 3085–3101 (1997) ? 1997 by John Wiley & Sons, Ltd. SYMMETRIC GALERKIN BOUNDARY INTEGRAL FORMULATION 3097 Figure 9. Finite element mesh (440 nodes and 400 elements) for the problem of Figure 8 Figure 10. Comparison of the nodal temperatures at the top semi-circle of radius = 1·0 obtained by the FEM and the SG-BEM Figure 11. Comparison of the nodal temperatures at the top semi-circle of radius = 0·25 obtained by the FEM and the SG-BEM ? 1997 by John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng., 40, 3085–3101 (1997) 3098 L. J. GRAY AND G. H. PAULINO semi-circles of radius 1·0 and 0·25, respectively. Note that the results obtained by the two methods are in good agreement. Figures 12 and 13 compare the nodal ux at the bottom semi-circles of radius 1·0 and 0·25, respectively. In both graphs, the FEM nodal thermal uxes are average nodal values, except at the endpoints where the corresponding element nodal thermal ux was used. The bullets at the endpoints of these graphs denote the average nodal ux. In Figure 13, there is a relatively small gap between the SG-BEM and the FEM curves. This is probably due to the coarse boundary element discretization of the bottom semi-circle of radius 0·25 (see Figure 7). Figure 14 shows a comparison of the computed nodal temperature along the interface, the two curves stemming from the left and right segments of the interface. It is interesting to observe that the results obtained by the two methods are within plotting accuracy. Figure 15 is the corresponding comparison for the interface ux, and three types of results obtained with the FEM are displayed in this gure. ‘FEM(average)’ means nodal averaging, which should not be used since there are di erent materials at each side of the interface. ‘FEM(average top)’ refers to the average nodal ux considering the elements above the interface. These results compare quite well with the SG-BEM. ‘FEM(average bottom)’ denotes the average nodal ux considering the elements Figure 12. Comparison of the nodal ux at the bottom semi-circle of radius = 1·0 obtained by the FEM and the SG-BEM. The bullets at the endpoints of the graph denote the average FEM nodal ux Figure 13. Comparison of the nodal ux at the bottom semi-circle of radius = 0·25 obtained by the FEM and the SG-BEM. The bullets at the endpoints of the graph denote the average FEM nodal ux Int. J. Numer. Meth. Engng., 40, 3085–3101 (1997) ? 1997 by John Wiley & Sons, Ltd. SYMMETRIC GALERKIN BOUNDARY INTEGRAL FORMULATION 3099 Figure 14. Comparison of the nodal temperatures along the LHS and RHS interfaces obtained by the FEM and the SG-BEM Figure 15. Comparison of the nodal ux along the LHS and RHS interfaces obtained by the FEM and the SG-BEM below the interface. This graph illustrates the danger of using average derivative values when di erent materials are present. In this case, care should be taken to consider appropriate averaging of nodal quantities. 6. CONCLUSIONS AND EXTENSIONS This paper has demonstrated that a fully symmetric Galerkin boundary integral equation formulation can be achieved for interface and multi-zone problems. The standard (i.e. single domain) SG formulation is based on use of the dual equations (i.e. the BIE and the HBIE) according to the type of boundary conditions. The double integration over the boundary, together with the symmetry properties of the kernels guarantees that the system matrices are symmetric. This approach on the non-interface boundaries, together with a strategic use of both equations on the interface, yields a symmetric system, while automatically incorporating the continuity conditions on the interface. The present interface formulation is therefore a natural extension of the SG procedure. The two most important bene ts provided by this boundary integral approach to interface problems are (a) a symmetric coupling of boundary and nite elements, and (b) a computationally ecient algorithm. In a single domain SG algorithm, the eciency arises from the linear algebra ? 1997 by John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng., 40, 3085–3101 (1997) 3100 L. J. GRAY AND G. H. PAULINO phase of the calculation, solving a symmetric versus a non-symmetric system. However, for an interface problem, the symmetry can be pro tably invoked in the matrix construction phase. In constructing the outer boundary equations, integrations over the interface can be bypassed, clearly a considerable bonus for this approach. Detailed timing comparisons with equivalent FEM calculations have not been performed in the work, but should be carried out. Numerical examples involving heat transfer in an eccentric annulus geometry, with several variations in the boundary conditions and material properties, were considered. These results were in good agreement with analytical solutions, when available, and with numerical solutions obtained with nite elements. Although a two-dimensional scalar eld formulation was discussed herein, there is no diculty in extending this method to three-dimensional and vector- eld problems. Fracture problems were also not considered, but this algorithm clearly provides an e ective framework for studying interfacial cracks. Coupling this approach with special elements for capturing crack tip singularities, or interface corner singularities (as occur in multi-material thermoelastic problems), should be possible, but requires further investigation. The interface formulation should also be directly applicable to non-linear analysis, in particular the symmetric Galerkin formulation for plasticity developed by Maier and co-workers.26; 4; 27 Finally, it is expected that this technique can be naturally incorporated into a general SG error estimation and adaptivity procedure.28 Investigations in these areas are currently being pursued. ACKNOWLEDGEMENT This research was partially supported by the Applied Mathematical Sciences Research Program of the Oce of Mathematical, Information, and Computational Sciences, U.S. Department of Energy under contract DE-AC05-96OR22464 with Lockheed Martin Energy Research Corp. LJG would like to thank Prof. M. R. Ramey for hospitality at UC-Davis while this work was performed. GHP acknowledges the “Junior Faculty Research Fellowship” from the UC-Davis Academic Senate. REFERENCES 1. J. N. Reddy and A. Miravete, Practical Analysis of Composite Laminates, CRC Press, Boca Raton, FL, 1995. 2. J. H. Kane, Boundary Element Analysis in Engineering Continuum Mechanics, Prentice-Hall, Englewood Cli s, N.J., 1994. 3. F. Hartmann, C. Katz and B. Protopsaltis, ‘Boundary elements and symmetry’, Ing.-Arch., 55, 440–449 (1985). 4. G. Maier, G. Novati and S. Sirtori, ‘On symmetrization in boundary element elastic and elastoplastic analysis’, in G. Kuhn and H. Mang (eds.), Discretization Methods in Structural Mechanics, Springer, Berlin, 1990, pp. 191–200. 5. C. A. Brebbia, J. C. F. Telles and L. C. Wrobel, Boundary Element Techniques, Springer, Berlin, 1984. 6. G. Krishnasamy, F. J. Rizzo and T. J. Rudolphi, ‘Hypersingular boundary integral equations: their occurence, interpretation, regularization and computation’, in P. K. Banerjee and S. Kobayashi (eds.), Developments in Boundary Element Methods—Advanced Dynamic Analysis, Vol. 7, Chapter 7, Elsevier, Amsterdam, 1991. 7. P. A. Martin and F. J. Rizzo, ‘Hypersingular integrals: how smooth must the density be?’, Int. J. Numer. Meth. Engng., 39, 687–704 (1996). 8. T. J. Rudolphi, ‘Higher order elements and element enhancement by combined regular and hypersingular boundary integral equations’, in B. S. Annigeri and K. Tseng (eds.), Boundary Element Methods in Engineering, Springer, Berlin, 1990, pp. 448–455. 9. K. Tomlinson, C. Bradley and A. Pullan, ‘On the choice of a derivative boundary element formulation using Hermite interpolation’, Int. J. Numer. Meth. Engng., 39, 451–468 (1996). 10. H. G. Walters, J. C. Ortiz, G. S. Gipson and J. A. Brewer, ‘Overhauser boundary elements in potential theory and linear elastostatics’, in T. A. Cruse (ed.), Advanced Boundary Element Methods, Vol. 16, Springer, Berlin, 1988, pp. 459–464. Int. J. Numer. Meth. Engng., 40, 3085–3101 (1997) ? 1997 by John Wiley & Sons, Ltd. SYMMETRIC GALERKIN BOUNDARY INTEGRAL FORMULATION 3101 11. W. S. Hall and T. T. Hibbs, ‘The treatment of singularities and the application of the Overhauser C 1 continuous quadrilateral boundary element to three dimensional elastostatics’, in T. A. Cruse (ed.), Advanced Boundary Element Methods, Vol. 16, Springer, Berlin, 1988, pp. 135–144. 12. L. J. Gray, L. F. Martha and A. R. Ingra ea, ‘Hypersingular integrals in boundary element fracture analysis’, Int. J. Numer. Meth. Engng., 29, 1135–1158 (1990). 13. G. Krishnasamy, L. W. Schmerr, T. J. Rudolphi and F. J. Rizzo, ‘Hypersingular boundary integral equations: Some applications in acoustic and elastic wave scattering’, J. Appl. Mech. ASME, 57, 404–414 (1990). 14. C. Balakrishna, L. J. Gray and J. H. Kane, ‘Ecient analytical integration of symmetric Galerkin boundary integrals over curved elements; thermal conduction formulation’, Comput. Methods Appl. Mech. Eng., 111, 335–355 (1994). 15. C. Balakrishna, L. J. Gray and J. H. Kane, ‘Ecient analytical integration of symmetric Galerkin boundary integrals over curved elements; elasticity formulation’, Comput. Methods Appl. Mech. Eng., 117, 157–179 (1994). 16. L. J. Gray, C. Balakrishna and J. H. Kane, ‘Symmetric Galerkin boundary integral fracture analysis’, Eng. Anal. Boundary Elements, 15, 103–109 (1995). 17. L. J. Gray and G. H. Paulino, ‘Symmetric Galerkin boundary integral fracture analysis for plane orthotropic elasticity’, Computational Mechanics, in press. 18. S. M. Holzer, ‘The symmetric Galerkin BEM in elasticity: FEM/BEM coupling, integration and adaptivity’, Comput. Mech., submitted. 19. L. J. Gray, ‘Boundary element method for regions with thin internal cavities’, Eng. Analy. Boundary Elements, 6, 180–184 (1989). 20. E. D. Lutz and L. J. Gray, ‘Exact evaluation of singular boundary integrals without CPV’, Commun. Numer. Meth. Engng., 9, 909–915 (1993). 21. F. J. Rizzo, D. J. Shippy and M. Rezayat, ‘A boundary integral equation method for radiation and scattering of elastic waves in three dimensions’, Int. J. Numer. Meth. Engng., 21, 115–129 (1985). 22. L. J. Gray and E. D. Lutz, ‘On the treatment of corners in the boundary element method’, J. Comput. Appl. Math. 32, 369–386 (1990). 23. J. B. Layton, S. Ganguly, C. Balakrishna and J. H. Kane, ‘A symmetric Galerkin multi-zone boundary element formulation’, Int. J. Numer. Meth. Engng., submitted. 24. S.-H. Hsieh, G. H. Paulino and J. F. Abel, ‘Recursive spectral algorithms for automatic domain partitioning in parallel nite element analysis’, Comput. Methods Appl. Mech. Eng., 121, 137–162 (1996). 25. M. D. Greenberg, Foundations of Applied Mathematics, Prentice-Hall, Englewood Cli s, N.J., 1978. 26. G. Maier and C. Polizzotto, ‘A Galerkin approach to boundary element elastoplastic analysis’, Comput. Methods Appl. Mech. Eng., 60, 175–194 (1987). 27. G. Maier, S. Miccoli, G. Novati and S. Sirtori, ‘A Galerkin symmetric boundary element methods in plasticity: formulation and implementation’, in J. H. Kane, G. Maier, N. Tosaka and S. N. Atluri (eds.), Advances in Boundary Element Techniques, Springer, Berlin 1993, pp. 288–328. 28. G. H. Paulino and L. J. Gray, ‘Error estimation and adaptivity for the symmetric Galerkin boundary element method’, in preparation. ? 1997 by John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng., 40, 3085–3101 (1997)