INTERNATIONAL JOURNAL FOR NUMERICAL METHODS IN ENGINEERING, VOL. 40, 2913—2931 (1997) A SYMMETRIC GALERKIN MULTI-ZONE BOUNDARY ELEMENT FORMULATION J. B. LAYTON,*s S. GANGULY,t C. BALAKRISHNA° AND J. H. KANE} Mechanical and Aeronautical Engineering Department, Clarkson ºniversity, Box 5725, Potsdam, N½, 13699, º.S.A. ABSTRACT The recent development of the symmetric Galerkin approach to boundary element analysis (BEA) has been demonstrated to be superior to the collocation method for medium to large problems. This fact has been shown in both heat conduction and elasticity. Accounts of collocation multi-zone analysis techniques have also been prevalent in the literature, dealing with multiple boundary integral relations associated with portions of overall objects. This technique results in overall system matrices with a blocked, sparse, but unsymmetric character. It has been shown that multi-zone techniques can produce smaller solution times than a single zone analysis for large problems. These techniques are useful for multi-material problems as well. This paper presents an approach for combining the benefits of both techniques resulting in a Galerkin multi-zone method, that is overall unsymmetric but contains a significant amount of block symmetry. A condensation technique in the multi-zone solver is shown to exploit the symmetry of the Galerkin formulation as well as the blocked sparsity of the multi-zone technique. This method is compared to collocation multi-zone on two elasticity problems from the literature. It is concluded that an appropriate implementation of the symmetric Galerkin multi-zone BEA indeed has the potential to be superior to the collocation based multi-zone BEA, for medium to large-scale elasticity problems. ( 1997 by John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng., 40, 2913—2931 (1997) No. of Figures: 7. KEY WORDS: No. of Tables: 2. No. of References: 22. boundary elements; symmetric Galerkin; multi-zone; condensation INTRODUCTION Some potentially useful alternative boundary element approaches have recently emerged in the technical literature. Some of these new approaches involve concepts including the employment of hypersingular boundary integral equations (HBIEs) and the Galerkin approach to BEA. The collocation approach to BEA, a well-known one, involves the satisfaction of boundary integral equations (BIEs) at a discrete number of source points. The Galerkin BEA (GBEA) approach is an alternative one that involves the satisfaction of the governing BIE in an integral or weighted * Correspondence to: Jeff Layton, Mechanical and Aeronautical Engineering Department, Clarkson University, Box 5725, Potsdam, NY, 13699, U.S.A. s Assistant Professor t Ph.D. Candidate ° Director, Computer Simulation, Center for the Advancement of Instruction in Science and Engineering (CAISE) Currently, Robust Engineer, GM Truck Group, 2000 Centerpoint Pkwy, Pontiac MI, 48431-3147, U.S.A. } Associate Professor, Director, Center for the Advancement of Instruction in Science and Engineering (CAISE) CCC 0029—5981/97/162913—19$17.50 ( 1997 by John Wiley & Sons, Ltd. Received 20 May 1996 Revised 17 December 1996 2914 J. B. LAYTON E¹ A¸. residual sense, leading to double surface integrations. If used strategically in concert with both singular and hypersingular BIEs, the Galerkin approach can generate symmetric sets of algebraic equations. Investigation of the potential for computational efficiency in the double integration and equation solving stages of the overall symmetric GBEA algorithm has been described in recent papers. Balakrishna et al.1,2 have employed curved continuous boundary elements for elasticity and heat conduction formulations, using a direct analytical treatment of the singular and hypersingular double integrations and employing a limit to the boundary definition to gain computational efficiency. They have shown that symmetric Galerkin with the analytical evaluation of integrals can produce superior results in terms of speed compared to the conventional collocation method for medium-to-large-scale problems, both in heat conduction and elasticity. Multi-zone BEA techniques break up the object of concern into separate zones. A multi-zone collocation BEA strategy can then be applied to each zone to produce the overall system matrix. This matrix has a blocked, sparse, but unsymmetric character, which makes efficient solution somewhat of a challenging task. Multi-zone BEA techniques are introduced in a number of texts.3—5 Several papers deal with the assembly and solution of these types of matrix equations.6—10 These references have shown that the multi-zone BEA strategy, when applied to collocation techniques, significantly extends the range of model shapes that can be successfully treated. At the same time they showed that substantive computational economy in the numerical integration, equation solving and response recovery phases of the overall analysis process is obtained. All these approaches describe a direct matrix block triangular factorization process that exploits the substantial block sparsity present in multi-zone analyses. Static condensation of degrees of freedom in BEA has also been investigated.11—13 The references described procedures that also couple boundary elements and finite elements in the same analysis. The boundary element matrices are condensed to just the degrees of freedom on the boundary element—finite element interface. Kane and his co-workers14,15 developed a collocation multi-zone analysis approach that optionally condenses the degrees of freedom exclusively in a single zone. Kane and Saigal16 and Kane et al.17 also describe a sparse blocked assembly and BEA equation solving procedure. This strategy incorporates boundary element sub-structuring in a completely arbitrary fashion. It allows for both condensed and uncondensed boundary element zones to consistently coexist in the same multi-zone problem. For problems with a large number of degrees of freedom, multi-zone solution techniques seem to be generally superior to the solution of an equivalent single zone problem.18,19 Such an approach takes advantage of the block sparsity of the multi-zone approach to more efficiently solve the system of equations. Moreover, multi-zone techniques can improve the conditioning of the resulting system matrices, improving the stability of the solution process.18 This paper presents an approach that allows for the melding of these two techniques: symmetric Galerkin boundary element analysis and multi-zone solution techniques. The resulting symmetric Galerkin multi-zone approach needs less storage than conventional collocation multi-zone methods. Also, as will be shown, it can result in substantial speed improvements for medium-to-large-scale problems. This paper will begin by briefly reviewing the elasticity Galerkin boundary formulation. It then continues on to the symmetric Galerkin boundary integral, hybrid analytical/numerical integration approach used in this work. A brief review of multi-zone techniques is then discussed, focusing on the use of this technique with conventional collocation approaches. Then the integration of the two approaches is discussed, paying particular attention to the insights that are needed for efficient coupling of the symmetric Galerkin and multi-zone techniques. Finally the Int. J. Numer. Meth. Engng., 40, 2913—2931 (1997) ( 1997 by John Wiley & Sons, Ltd. SYMMETRIC GALERKIN MULTI-ZONE BOUNDARY ELEMENT FORMULATION 2915 approach is demonstrated on two elasticity problems involving two and three zones respectively to investigate the efficacy of the method. The general efficiency of the symmetric Galerkin multi-zone technique is also discussed by increasing the number of zones with a fixed number of elements for a plane elasticity problem. ELASTICITY GALERKIN BOUNDARY FORMULATION For isotropic elasticity problems, the governing BIE is Somigliana’s identity, P! ti (x)uik (x, d) d!(x)"P! tik (x, d)ui (x) d!(x)#uk (d) (1) In this expression, u (x, d) and t (x, d) are the ith displacement and traction fundamental (Kelvin) ik ik solution components at a sample point x due to a point load at d in the direction of the kth co-ordinate axis. It should be noted that d is assumed to be located in the interior of the domain ), and ! is the boundary and the summation convention is implied for repeated subscripts. The spatial gradient components of u (d) can be obtained by differentiating this expression with k respect to d . These components can be employed with the strain displacement relations and i utilized in the stress strain law to obtain a stress component BIE,18 P! Dijr (x, d)tr (x) d!(x)!P! Sijr (x, d)ur (x) d!(x) !p (d)" ij (2) The Cauchy stress transformation is used to convert this second rank tensor expression into a traction vector BIE, P! G(1)ik (x, d)ti (x) d!(x)!P! G(2)ik (x, d)ui (x) d!(x) t (d)" k (3) where G(1) (x, d)"!D (x, d)n (d) ik kji j (4) G(2) (x, d)"!S (x, d)n (d) ik kji j The fundamental solutions u , t , D , and S are discussed in Reference 18. In a 2-D plain stress ik ik ijk ijk or strain context, the strongly singular function G(1) behaves essentially as Dx!dD~1, while ik the hypersingular kernel G(2) behaves as Dx!dD~2, for small Dx!dD. An important observation ik to be made here is that the roles of x and d are symmetric in the hypersingular fundamental solution G(2) . ik SYMMETRIC GALERKIN BOUNDARY INTEGRAL APPROACH An approach that can produce symmetry in the BEA coefficient matrices requires that the conventional governing BIE must be multiplied by a weight function and integrated with respect to d around ! for a second time. In the Galerkin approach, one chooses this weight to be the same as the boundary element interpolation function associated with a particular node, so that x and d are treated symmetrically. Equation (3) is a hypersingular(traction) BIE that can be utilized in ( 1997 by John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng., 40, 2913—2931 (1997) 2916 J. B. LAYTON E¹ A¸. a GBEA approach. In this case, the matrix with the symmetric coefficients becomes the one that premultiplies the vector of nodal displacements. This follows from the fact that the fundamental solution G(2) present in the hypersingular BIE (HBIE) is symmetric with respect to x and d. For ik mixed boundary value problems, BEA matrix column exchanges are required to move all unknowns to the left-hand side of the overall system of equations. This maneuver will destroy the symmetric character of the coefficient matrix suggested above. However, the singular BIE shown in equation (1) contains a term involving the (unsymmetric) fundamental solution t that allows ik the symmetry to be restored in a mixed boundary value problem. The details of restoring symmetry by making use of the unsymmetric kernel functions present in the displacement and traction BIEs are described in References 1 and 2. The surface response can be approximated via standard interpolation functions, starting with the traction BIE given in equation (3). As the HBIE is now multiplied by one of these interpolation functions h(r) and integrated a second time with respect to d, it will be necessary to interpolate the surface response with respect to both x and d as u (x)"h(N)(a)u(N) i i (5) t (x)"h(N)(a)t(N) i i u (d)"h(M)(b)u(M) i i (6) t (d)"h(M)(b)t(M) i i The summation convention is again implied over the repeated superscripts, where the range of the superscripts is from one to the number of nodes on the element. The symbol a is used to denote the intrinsic variable of integration with respect to x on element E, while b is the intrinsic variable on an element ¸ for which h(r) is non-zero. The same interpolation functions are employed for both displacements and tractions. Tractions at shared nodes in adjacent elements are allowed to take on independent values, while displacements at nodes shared by adjacent elements are constrained to be the same. Following these steps, the following expression for equation (3) is developed: Term (1)"Term (2)!Term (3) kr kr kr (7) where Term (1)" + kr L3j Term (2)" + kr L3 j Term (3)" + kr L3 j GP GP GP `1 ~1 `1 NEL h(r)(b) + ~1 E/1 `1 H h(r)(b)h(M)(b)t(M) J(b) db k NEL h(r)(b) + ~1 E/1 AP AP (L) (8) B B H H `1 (E) (L) G(1) (x, d)h(N)(a)t(N) J(a) da J(b) db ik i ~1 (9) `1 (E) (L) G(2) (x, d)h(N)(a)u(N) J(a) da J(b) db ik i ~1 (10) In this expression, ¸ symbolizes the element being integrated with respect to d, and j is the set of these elements for which h(r) is non-zero. For continuous elements, j will contain more than one Int. J. Numer. Meth. Engng., 40, 2913—2931 (1997) ( 1997 by John Wiley & Sons, Ltd. SYMMETRIC GALERKIN MULTI-ZONE BOUNDARY ELEMENT FORMULATION 2917 element for interpolation functions associated with edge nodes. Moreover, for continuous elements, element contributions to the overall matrix equations associated with edge nodes on adjacent elements get added. Thus, the columns associated with such edge nodes are assembled from contributions from adjacent elements. Symmetric treatment of the rows in the Galerkin approach requires that the definition be generalized so that the interpolation functions mean the union of all interpolation functions associated with a particular node. Thus, the rows associated with an edge node on a continuous element are now also assembled from element contributions associated with this edge node. The subscript k and r remain free indices, representing the co-ordinate direction of the unit load and the node number associated with the interpolation function used as the weighting function of the outer integration, respectively. A direct analytical treatment of the singular and hypersingular double integrations, employing a limit to the boundary definition, was developed in References 1 and 2 to reduce the computational burden associated with the numerical accurate evaluation of the double singular integrations present in this approach. The analytical regularization emplıoyed in these references separates the potentially singular Galerkin integrands into an essentially singular but simple part, plus a remainder that can be integrated numerically. The references go on to illustrate that such a hybrid analytical/numerical scheme can result in faster solutions for symmetric Galerkin BEA than collocation techniques when applied to moderate-to-large-scale problems. UNSYMMETRIC MULTI-ZONE BOUNDARY ELEMENT ANALYSIS Multizone BEA3—10 is accomplished by breaking up an entire boundary element model into portions, referred to as zones. The surface of each zone totally encloses a subvolume of the entire model and may generally contain boundary elements and nodes that belong exclusively to a single zone, and boundary elements and nodes that are on an interface with other zones. The governing boundary integral relationships can then be written for each zone. By evaluating the discretized BIE at a set of locations in the BIE corresponding to the node-point locations (in a collocation context) for the zone in question, a matrix system of equations can be generated for each zone, [F(i)]Mu(i)N"[G(i)]Mt(i)N (11) where [F(i)] and [G(i)] are square and rectangular coefficient matrices, respectively, and Mu(i)N and Mt(i)N are column vectors of node-point displacement and traction components, respectively, of conformable size. The superscript (i) means no sum and that the quantities are associated with the ith zone. The rectangular nature of [G(i)] is due to jumps in specified traction that can exist when continuous boundary elements are employed. The matrix isolated relations for individual zones are combined for use in an overall analysis by employing the conditions of displacement continuity and force equilibrium at the zone interfaces. The continuity condition requires that nodal displacements, calculated for zone-i at an interface between zone-i and zone-j, must equal the nodal displacement components calculated in zone-j on the same interface. A similar relationship exists due to equilibrium considerations for components of the traction vector at interface nodes between the two zones, except that a negative sign must be present to account for the opposite directions of the outward surface ( 1997 by John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng., 40, 2913—2931 (1997) 2918 J. B. LAYTON E¹ A¸. normal in the two zones in question as shown in equations (12) and (13): Mu(i) N"Mu(j) N ij ij (12) Mt(i) N"!Mt(j) N ij ij (13) The double subscript notation is used to convey that the vector in question is a collection of components entirely on the interface between zone-i and zone-j. For collocation-based multi-zone methods, there is a particular ordering approach to assembling the zones, which has been shown to yield good results.18,19 Consider the small BEA model as shown in Figure 1. This two-zone model is composed of three-node discontinuous elements. The first zone in this model has four elements, 12 nodes, and 24 displacement degrees of freedom, while the second zone has five elements, 15 nodes, and 30 displacement degrees of freedom. Both the displacement and tractions are unknown at the interface. The matrix system of equations for the first zone can be formed with an ordering so that the six degrees of freedom on the interface are put last in the matrices and vectors for that zone. This can be done for both the left side [F]MuN and the right side [G]MtN of the discretized BIE. The matrix—vector product can be expanded to include the 24 degrees of freedom associated exclusively with zone-2. Note that although zone-2 has 30 degrees of freedom associated with it, six of these are shared with zone-1. Note further that, to retain the correct relationship among the quantities expressed in the original zone-1 discretized BIE, the coefficients associated with the degrees of freedom exclusively in zone-2 in the expanded relationship must be zero. Thus, the zone-1 matrix system of equations has been partitioned into blocks. The [F ] block contains the columns of [F] with coefficients that multiply degrees of 11 freedom exclusively in zone-1. The [F ] block contains columns with coefficients that multiply 12 degrees of freedom on the zone-1, zone-2 interface. An analogous type of reordered and expanded matrix equations can be written for zone-2. The integration of the two zones is accomplished by first writing the degrees of freedom exclusively in zone-1, followed by the interface degrees of freedom, and then the degrees of freedom exclusively in zone-2. The resulting partitioned matrix equations for the two zones are [F1]Mu1N"[G1]Mt1N (14a) Mt1 N Mu1 N 11 11 [[F1 ] [F1 ] [0]] Mu1 N "[[G1 ] [G1 ] [0]] Mt1 N 12 11 12 11 12 12 Mt222N Mu222N (14b) Figure 1. Two-zone BEA model Int. J. Numer. Meth. Engng., 40, 2913—2931 (1997) ( 1997 by John Wiley & Sons, Ltd. SYMMETRIC GALERKIN MULTI-ZONE BOUNDARY ELEMENT FORMULATION 2919 [F2]Mu2N"[G2]Mt2N (15a) Mt1 N Mu1 N 11 11 [[0] [F2 ] [F2 ]] Mu2 N "[[0] [G2 ] [G2 ]] Mt2 N 22 12 12 12 12 11 N Mu2 Mt222N 22 (15b) where the superscript indicates the zone. The resulting two blocked matrix equations can be rewritten with the unknown Mt N brought to the left-hand side: 12 G H G H GH GH G H GH Mu1 N Mt1 N 11 11 Mu1 N M0N 12 "[[G1 ] [0] [0] [0]] [[F1 ] [F1 ] ![G1 ] [0]] 12 12 11 11 Mt1 N M0N 12 Mu2 N M0N 22 (16) Mu1 N M0N 11 Mu2 N M0N 12 "[[0] [0] [0] [G2 ]] [[0] [F2 ] ![G2 ] [F2 ]] 22 12 12 22 Mt2 N M0N 12 Mu2 N Mt2 N 22 22 (17) It is important to notice that the vectors on the left-hand side are not the same in equations (16) and (17). The continuity and equilibrium relationships shown in equations (12) and (13) can be used to make the vector on the left-hand side of equation (17) the same as that contained in equation (16): Mu1 N M0N 11 Mu1 N M0N 12 "[[0] [0] [0] [G2 ]] [[0] [F2 ] #[G2 ] [F2 ]] 12 12 22 22 Mt1 N M0N 12 Mu2 N Mt2 N 22 22 (18) equations (16) and (18) can now be combined to form a single matrix equation by appending the rows of equation (18) to the bottom of equation (16): C [0] [F1 ] [F1 ] ![G1 ] 12 12 11 [0] [F2 ] #[G2 ] [F2 ] 22 12 12 D G H Mu1 N 11 Mu1 N [0] [G1 ] [0] [0] 12 " 11 Mt1 N [0] [0] [0] [G2 ] 12 22 Mu2 N 22 C D GH Mt1 N 11 M0N M0N (19) Mt2 N 22 There are several techniques for solving the above blocked, sparse, unsymmetric equations such as block-based-solvers,9,18,19 lumping techniques,16,20 and iterative solvers.21 However, all of them have to deal with the efficient solution of unsymmetric systems of equations. ( 1997 by John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng., 40, 2913—2931 (1997) 2920 J. B. LAYTON E¹ A¸. SYMMETRIC GALERKIN MULTI-ZONE FORMULATIONS For collocation multi-zone problems, a general goal is to solve blocked, sparse, unsymmetric equations as efficiently as possible. One of the key problems in achieving this goal is the efficient handling of the unsymmetric blocks in the solver. Since they are unsymmetric, the entire block must be stored and many more operations are required for block factorization or solution as compared to symmetric systems of equations. This section of the paper will discuss how symmetric Galerkin techniques can be integrated with multi-zone methods to create symmetric blocks that deal exclusively with a single zone. Within the final system of equations, such as those in equation (19), the ordering of the blocks is arbitrary. For any ordering of the blocks, it is always possible to formulate the relevant [F] and [G] matrices. However, certain orderings are better than others. A typical ordering for collocation problems was discussed in the previous section. That particular ordering was chosen such that the non-empty blocks are as close as possible to the main diagonal. This minimizes the block fill-in effect that occurs during the solution of these sparse blocked equations.18,19 The order is determined by first simply listing all permutations of the two zones as shown below: 11 12 21* 22 (20) For three zones the permutations would be the following: 11 12 13 21* 22 23 31* 32* 33 (21) For permutations where the first digit is less than or equal to the second digit, blocks of displacements are generated; otherwise, blocks of traction components are generated. The permutations associated with blocks of traction components are shown with an asterisk in the above list. Using symmetric Galerkin BEA, as described in References 1 and 2 in a multi-zone sense, symmetry is restricted to the degrees of freedom associated with a particular zone only. The blocks associated with the interfacial degrees of freedom will not be symmetric. Therefore, some of the matrices in equation (19) will be symmetric (i.e. F1 , F2 ) while the remaining matrices, 11 22 which are associated with the interfacial degrees of freedom, are unsymmetric. A natural way to combine substructuring with multi-zone BEA for symmetric Galerkin is to allow for the possible condensation of degrees of freedom that appear exclusively in any particular zone.14,15 In this case, the partitions to be eliminated by the condensation process coincide exactly with certain partitions already present in the multi-zone analysis procedure. When zones are selected for condensation, a second level assembly procedure for the formation of overall sparse blocked system of equations is used to assemble condensed or uncondensed zone contributions to the overall matrix. For symmetric Galerkin multi-zone, it becomes very logical to condense out the symmetric portions of the blocked sparse system of equations, which are associated with degrees of freedom that appear exclusively in any particular zone. This is because the computations present in such condensation is dominated by the factorization of the left-hand side block matrix. A symmetric factorization is roughly twice as fast as the factorization of an unsymmetric block of the same size. For collocation multi-zone BEA these matrices are unsymmetric, and condensation, while improving the overall solution process, adds a large burden to the computational task. Int. J. Numer. Meth. Engng., 40, 2913—2931 (1997) ( 1997 by John Wiley & Sons, Ltd. SYMMETRIC GALERKIN MULTI-ZONE BOUNDARY ELEMENT FORMULATION 2921 However, the degree of freedom ordering as discussed previously for the collocation-based multi-zone does not take full advantage of all the symmetry present in the symmetric Galerkin multi-zone BEA. For the two zone problem defined in equation (20), the zone noted as 11 is symmetric while the zone noted as 12 is unsymmetric. Any condensation beyond condensing zone 11 will not take full advantage of the symmetry present from the symmetric Galerkin BEA, because the resulting matrices from condensing zone 11 are unsymmetric and no further benefits are gained from the symmetry. Consequently, it is necessary to change the ordering of the degrees of freedom so that all the symmetry can be exploited. In this paper, a new approach is described that reorders the degrees of freedom so that the degrees of freedom that appear exclusively in any particular zone are listed first for all the zones, one after the other. The degrees of freedom associated with the interfaces come at the end. For example, the new ordering for a two-zone problem would be 11 22 12 (22) and the new order for a three-zone problem is as follows: 11 22 33 12 13 23 (23) This new ordering can have a great deal of impact on the solution of the problem. This is due to the fact that the 11 and 22 zones are both symmetric, which allows the condensation process to proceed much faster than if they were unsymmetric. The block left after condensation is unsymmetric. The effect of the re-ordering can be seen in Figure 2, which is a plot of the uncondensed matrix population associated with the model used in the first design study. Notice the symmetric blocks associated with the two zones exclusively. To elaborate on the procedure, consider a generic two-zone problem as illustrated in Figure 1. The resulting equations can be collected in the form shown in equation (24) with the ordering Figure 2. Plot of the matrix population ( 1997 by John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng., 40, 2913—2931 (1997) 2922 J. B. LAYTON E¹ A¸. appropriate for symmetric Galerkin multi-zone, [A ] [0] [A ] Mx N Mb N 11 13 1 1 [0] [A ] [A ] Mx N " Mb N 22 23 2 2 [A ] [A ] [A ] Mx N Mb N 31 32 33 3 3 (24) where Mx N are the degrees of freedom exclusively associated with zone 1, Mx N are the degrees of 1 2 freedom exclusively associated with zone 2, and Mx N are the degrees of freedom along the 3 interface. The degrees of freedom associated with Mx N can be condensed out of the global 1 problem by considering the matrix equation associated with the first row: [A ]Mx N#[A ]Mx N"Mb N 11 1 13 3 1 (25) Solving for the degrees of freedom exclusively associated with zone 1 results in an expression for Mx N: 1 Mx N"[A~1]Mb N![A~1][A ]Mx N 11 13 3 11 1 1 (26) Notice that the submatrix [A ] is symmetric through the symmetric Galerkin formulation. It is 11 important to note that the matrix inverse need never be computed nor used in the condensation. Its use in the equations is purely symbolic. When implementing the equations, the matrix inverse always premultiplies either a column vector or a rectangular matrix. Consequently, the result of the multiplication can be determined without ever having to compute the inverse.18,22 Now consider the third matrix equation in equation (24) and substitute equation (26) into the expression [A ]([A ]~1Mb N![A ]~1[A ]Mx N)#[A ]Mx N#[A ]Mx N"Mb N 31 11 1 11 13 3 32 2 33 3 3 (27) Collecting terms, [A ]Mx N#([A ]![A ][A ]~1[A ])Mx N"Mb N![A ][A ]~1Mb N 32 2 33 31 11 13 3 3 31 11 1 (28) Symbolically this can be written as [A ]Mx N#[A ]*Mx N"Mb N* 32 2 33 3 3 (29) where the superscript * indicates those matrices or vectors that have been modified due to the condensation. The condensed matrix system is now C DG H G H [A ] [A ] 22 23 [A ] [A ]* 32 33 Mx N Mb N 2 " 2 Mx N Mb N* 3 3 (30) From the first row of equation (30) Mx N can be found, 2 Mx N"[A ]~1Mb N![A ]~1[A ]Mx N 2 22 2 22 23 3 Int. J. Numer. Meth. Engng., 40, 2913—2931 (1997) (31) ( 1997 by John Wiley & Sons, Ltd. SYMMETRIC GALERKIN MULTI-ZONE BOUNDARY ELEMENT FORMULATION 2923 where [A ] is symmetric due to the symmetric Galerkin formulation. Substituting equation (31) 22 into the last row of equation (30) results in the following: [A ]([A ]~1Mb N![A ]~1[A ]Mx N)#[A ]*Mx N"Mb N* 32 22 2 22 23 3 33 3 3 (32) Collecting terms, ([A ]*![A ][A ]~1[A ])Mx N"Mb N*![A ][A ]~1Mb N 33 32 22 23 3 3 32 22 2 (33) Symbolically this can be written as [A ]**Mx N"Mb N** 33 3 3 (34) where [A ]**"[A ]*![A ][A ]~1[A ] 33 33 32 22 23 (35) Mb N**"Mb N*![A ][A ]~1Mb N 3 3 32 22 2 The double ** indicates the matrices and vectors that have been modified twice by condensation. The matrix [A ]** is unsymmetric. Solving for Mx N from equation (34) one can substitute back 33 3 into equation (31) to compute Mx N, then one can substitute Mx N back in equation (26) to get Mx N. 2 3 1 It is important to note that the block factorization process described above can be done in place. The computer memory used to store the block entries of the original matrix can be reused to store the matrix blocks as they evolve.18 The only extra memory needed in this blocked factorization process is for the storage of the largest matrix that is encountered in the factorization process. The factorization and condensation process for the symmetric Galerkin multi-zone BEA requires much less computer time and storage than the collocation approach to multi-zone since only half of the storage of collocation is required. Furthermore, the condensation strategy just described is useful for problems with non-linear boundary conditions. In this case, the problem would be zoned about the non-linearity. Then the degrees of freedom not associated with the non-linearity would be condensed out. Then a solution strategy for the non-linear portion of the problem can be used without carrying along the baggage associated with the other, linear degrees of freedom. DESIGN STUDIES This section presents two design studies that are used to examine the efficacy of the approach. Two different problems are used, a two-zone problem, and a three-zone problem. The first example problem involves a rectangular plate in plane stress, subjected to a uniform traction on one end and clamped at the other end as shown in Figure 3. A two-zone model is used for this case. The second problem is a rectangular plate with a centrally located hole in plane stress. The plate is subjected to a uniform traction pulling on the left- and right-hand sides of the plate. Only one quarter of the plate is modeled since the problem has symmetry in two directions. The non-traction boundary conditions are formulated to ensure symmetry of the results (i.e. simply supported in the x-direction along the vertical face, and simply supported in the y-direction along the horizontal face). ( 1997 by John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng., 40, 2913—2931 (1997) 2924 J. B. LAYTON E¹ A¸. Figure 3. Plane surface under uniform traction Figure 4. Quarter symmetry model of plate with a hole under uniform traction Both the problems are test problems from the multi-zone literature.22 Both specimens use a common material with a Young’s modulus of elasticity of E"0·3]105 units and a Poison’s ratio l"0·3. A standard traction value of 1000 is used. To gain confidence in the ability of the present method to solve accurately the problems in question, the response using the symmetric Galerkin multi-zone BEA method is compared to the response computed using a collocation multi-zone method. Two-zone plane stress problem Table I presents the resulting x displacement and local principal stress p for selected 22 locations for the collocation multi-zone method using 50 elements and the present method using 50 elements. The agreement is quite good for this model. The next step is to examine the benefits of the method for efficient solution of multi-zone problems. A series of BEA meshes were used in this study. These meshes begin with a coarse mesh of 8 elements and increases to 300 elements. The examination is done by calculating the times needed to compute the solution for the multi-zone collocation BEA method and the symmetric Galerkin BEA multi-zone. The times are broken down into several components (Note: the term NM indicates the time was not measurable using the clock resolution of the compiler). There are separate numbers for the integration and solution phase for the different solution techniques. In addition, there are separate numbers for the condensation phase and the solution phase for the symmetric Galerkin multi-zone technique. Note that the collocation multi-zone technique also uses condensation, but the time for the condensation has not been presented for this study. Int. J. Numer. Meth. Engng., 40, 2913—2931 (1997) ( 1997 by John Wiley & Sons, Ltd. 2925 SYMMETRIC GALERKIN MULTI-ZONE BOUNDARY ELEMENT FORMULATION Table I. x displacement and stress response of the two-zone problem Location X ½ Collocation (x displacement) Galerkin (x displacement) Collocation (p ) 22 6·0 9·0 12·0 9·0 6·0 17·0 22·0 19·5 17·0 14·5 0·0 0·0 3·0 6·0 6·0 0·0 3·0 6·0 6·0 6·0 0·1977 0·2935 0·3975 0·2935 0·1977 0·5642 0·7309 0·6442 0·5642 0·4775 0·1979 0·2938 0·3979 0·2938 0·1979 0·5645 0·7312 0·6445 0·5645 0·4778 989·91 998·89 2·543 998·9 989·9 999·6 2·846 997·8 991·6 999·7 Galerkin (p ) 22 988·4 997·8 2·136 995·4 986·1 998·3 2·965 998·4 992·3 998·3 Figure 5. Computer times for collocation and symmetric Galerkin multi-zone for the two-zone problem Rather, it is included in the solution time. All computations were carried out in double precision on the same computer with a UNIX operating system and all programs were compiled using the same FORTRAN 77 compiler and compiler options. The component and total computer times are plotted as a function of the number of elements in the BEA model in Figure 5. Examining the data, it becomes clear that for the collocation method, the solver takes more time than the integration. Yet, as the number of elements increases, the amount of time spent in the solver increases dramatically. On the other hand, for the symmetric Galerkin approach, more time is spent in the integration phase than in the solution phase (including condensation). However, the solution phase does not increase as dramatically with the number of elements in the model as for the collocation method. It appears that this behaviour is ( 1997 by John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng., 40, 2913—2931 (1997) 2926 J. B. LAYTON E¹ A¸. due to the pronounced growth of the time spent in the unsymmetric matrix factorization step of the collocation BEA relative to the computation occurring in the symmetric matrix factorization in the symmetric Galerkin multizone BEA. The figure is also examined for ‘cross-over’, which is the point where the present method is faster than the collocation multi-zone method. For this example, after models with about 235 elements, the symmetric Galerkin multi-zone BEA becomes faster than collocation multi-zone BEA even when the integration is included in the total time. In fact, at 300 total elements, symmetric Galerkin multi-zone BEA is about 6 per cent faster than collocation multi-zone BEA. Three-zone design study Table II presents the resulting x displacement and local principal stress p at a variety of 22 locations for the collocation method and the Galerkin method, both using 200 elements. Again, the agreement is quite good for this model. Timing studies were conducted with this problem using the series of BEA meshes. This problem also considers a series of BEA meshes. It begins with a starting coarse BEA mesh of 20 elements and increases to 400 elements. The procedure for the timing studies is exactly the same as the two zone problem. Figure 6 presents the results of the timing studies for the series of BEA meshes for the collocation multi-zone and the present method. In examining the data, it becomes apparent that the three-zone problem has smaller total solution times for the same total number of elements as compared to the two-zone problem considered previously. In comparing Figures 5 and 6 it can be seen that as the number of zones increases, the amount of time spent in both the solver phase and the integration decreases for both collocation and Galerkin approaches. For the three-zone problem, since the number of zones is increased, a larger number of operations due to matrix multiplications needs to be performed related to the condensation (see equations (27)—(35) Table II. Displacement and stress response of the three-zone problem Location X 0·55 1·02 1·71 2·00 2·10 2·25 2·75 2·06 3·0 3·5 5·5 4·9 7·0 7·0 ½ Collocation (x displacement) Galerkin (x displacement) Collocation (p ) 22 Galerkin (p ) 22 0·83 1·71 1·02 0·0 6·0 0·0 0·0 3·43 6·0 6·0 0·0 6·0 0·0 6·0 0·05912 0·04560 0·09059 0·11513 0·06652 0·11958 0·13051 0·07421 0·09585 0·11373 0·21075 0·16145 0·25880 0·23172 0·05923 0·04582 0·09148 0·11412 0·06643 0·11931 0·13050 0·07412 0·09562 0·11371 0·21071 0·16141 0·25860 0·23169 260·1 893·9 370·5 787·5 699·7 582·6 722·5 787·6 699·7 582·6 972·4 1016·4 995·8 998·8 259·9 891·4 369·3 785·9 698·3 581·9 721·7 786·3 698·3 581·7 971·8 1014·5 994·4 997·6 Int. J. Numer. Meth. Engng., 40, 2913—2931 (1997) ( 1997 by John Wiley & Sons, Ltd. SYMMETRIC GALERKIN MULTI-ZONE BOUNDARY ELEMENT FORMULATION 2927 Figure 6. Computer times for collocation and symmetric Galerkin multi-zone for the three-zone problem for the specific equations). These operations are the same for the collocation method and offset the benefits achieved by block symmetric matrix factorization, in the present formulation, as compared to the block unsymmetric matrix factorization for the multi-zone collocation BEA method. As a result, it takes around 400 elements for the symmetric Galerkin BEA to become faster than the collocation multi-zone BEA approach. After that, there is an exponential rise in the computational difference between the collocation and the Galerkin method. At about 500 elements, for the three-zone problem, the Galerkin approach is 30 per cent faster. Although not explicitly shown here, the three-zone problem has a larger percentage of interfacial degrees of freedom than two-zone problem for an equivalent number of total DOF. This is typical when increasing the number of zones. As the number of interface degrees of freedom increases, the time spent in basic matrix multiplications associated with the condensation increases. In addition, for the symmetric Galerkin multi-zone BEA, the final unsymmetric matrix is larger in size. Thus, the advantages derived due to symmetry are exploited to a lesser degree. In essence, as the number of zones increases, the ‘cross-over’ point in total solution time between collocation multi-zone BEA and symmetric Galerkin BEA, increases. The efficacy of multi-zoning with symmetric Galerkin As shown in Reference 18, multi-zoning itself can improve the computational efficiency of the solver. For example, if the number of elements are held constant and the number of zones is increased, the time required for solving the problem including the integration decreases. This section of the paper presents a similar argument for the combination of symmetric Galerkin and multi-zoning to examine the efficacy of multi-zoning in this context. The first design study is divided into a number of zones to study the efficacy and efficiency of the multi-zone method over the single zone method for the multi-zone and symmetric Galerkin techniques described in this paper. The integration time, solution time, condensation time, and the total time have been computed and are presented in Figure 7. ( 1997 by John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng., 40, 2913—2931 (1997) 2928 J. B. LAYTON E¹ A¸. Figure 7. Computational time versus number of zones for the first design problem In examining Figure 7, it can be seen that there is a sharp drop in the required integration time from a single zone to two zones. This is due to the fact that, for symmetric Galerkin, the integration time varies as the square of the number of elements (N2). The magnitude of the reduction in integration time can then be estimated via the following equation: Reduction"M A B N2 M2 (36) where N is the number of elements (held constant) and M is the number of zones. If one proceeds from a single zone to two zones holding the number of elements N constant, the reductions in computational timings are approximately half (2](N2/4)). For more than two zones the magnitude of the reduction in the integration times is much less than the two-zone case. Increasing the number of zones has diminishing returns in the reduction in integration times via equation (36). This is due to the fact that increasing the number of zones increases the number of elements at the interface. The general trend of computational time with number of zones also carries over into solution times. The difference being that the solution times are generally a function of the cube of the number of elements approximately (N3). However, it should be noted here that the final solution time is extremely small due to the small size of the matrix and increases gradually for increasing number of zones (see Figure 7). In the present paper, an initial domain is artificially subdivided by arbitrary interfaces. However, in a large number of problems the interfaces are already predefined. This is particularly true for problems with a combination of materials or for composite materials. In such cases, the method followed in this paper remains useful. For some specific problems which do not have standard boundary conditions (e.g. non-linearities, non-perfect bonding, Robin conditions, etc.), it will be advantageous to use the condensation techniques described in this paper. Since the Int. J. Numer. Meth. Engng., 40, 2913—2931 (1997) ( 1997 by John Wiley & Sons, Ltd. SYMMETRIC GALERKIN MULTI-ZONE BOUNDARY ELEMENT FORMULATION 2929 condensed portions of the problem can be defined so they do not contain non-linearities, the present approach retains some advantages over collocation multi-zone for medium-tolarge-scale problems. It can be applied directly to condense out the symmetric portions of the matrices and focus directly on the interface, non-symmetric or the non-linear portion of the problem. The efficacy of multi-zoning, illustrating the effects of increasing the number of zones on computational time, has already been demonstrated in the literature for collocation. The reproduction of the analysis for symmetric Galerkin multi-zoning is necessary to understand some of the inherent benefits and limitations of the present approach. In particular, the results indicate that the benefits of multi-zoning with symmetric Galerkin in the present approach decrease with increasing number of zones due to the increased number of interfacial degrees of freedom. This is due to the fact that these degrees of freedom are associated with unsymmetric portions of the final system of equations, so block symmetry cannot be exploited. CONCLUSIONS A new formulation for elastic response analysis of continuum solids has been presented that facilitates the combined utilization of symmetric Galerkin BEA and sparse multi-zone solution concepts. Symmetric Galerkin BEA is utilized because its biggest benefit is that the resulting matrix system of equations for a single zone are symmetric. When symmetric Galerkin BEA was applied to the multi-zone technique, a reordering of the degrees of freedom from the standard ordering for collocation multi-zone was developed. This reordering allows all the symmetry to be exploited. In the blocked sparse system of equations, the blocks that are associated with degrees of freedom that appear exclusively in any particular zone are symmetric and the blocks associated with the interfacial degrees of freedom are unsymmetric. Condensation of the sparse symmetric blocks results in a final small unsymmetric condensed system of equations associated with the interfacial degrees of freedom. This small system of equations is very easy to solve compared to solving a blocked, sparse, unsymmetric system of equations as in the collocation multi-zone case. The multi-zoning in this paper dealt with the creation of artificial boundaries. However, when the boundaries are pre-defined as in the case of composite materials, the method suggested in the paper can be easily applied. Once the zones have been condensed, the number of operations necessary to obtain the final matrix is essentially the same for the collocation multi-zone BEA approach and the approach presented in this paper. As in collocation multi-zone, the combination of symmetric Galerkin and multi-zone retains the basic advantages of multi-zoning. The present paper has demonstrated that increasing the number of zones reduces the computational time and the integration time, but with diminishing returns with increasing number of zones due to the increasing number of interfacial DOFs. However, the benefits of the present formulation are only derived in the symmetric factorization during the condensation. As the number of zones increases the benefits derived from this effect become less pronounced. Consequently, the effect of symmetric condensation on total computational times decreases. This shifts the computational burden to the unsymmetric solver step. In essence, this increases the element cross-over point between the present method and the collocation multi-zone BEA. To rectify this problem, what is needed is a completely symmetric formulation. ( 1997 by John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng., 40, 2913—2931 (1997) 2930 J. B. LAYTON E¹ A¸. In conclusion, this paper demonstrates that the multi-zone symmetric Galerkin BEA formulation can be more computationally efficient than its unsymmetric multi-zone collocation-based BEA counterpart for medium-to-large-scale elasticity problems. However, the degree of computational improvement depends upon the number of zones in the problem and the number of interfacial degrees of freedom. The condensation techniques presented here are useful for arbitrarily defined zones, pre-defined zones, combinations of arbitrary and pre-defined zones, or problems with non-linearities. ACKNOWLEDGEMENT The second author was supported by a fellowship from United Technologies (Pratt and Whitney and Hamilton Standard). REFERENCES 1. C. Balakrishna, L. J. Gray and J. H. Kane, ‘Efficient analytical integration of symmetric galerkin boundary integrals over curved elements; thermal conduction formulation’, Comput. Methods Appl. Mech. Eng., 111, 335—355 (1994). 2. C. Balakrishna, L. J. Gray and J. H. Kane, ‘Efficient analytical integration of symmetric galerkin boundary integrals over curved elements; elasticity formulation’, Comput. Methods Appl. Mech. Eng., 117, 157—179 (1994). 3. P. K. Banerjee and R. Butterfield, Boundary Element Methods in Engineering and Science, McGraw-Hill, London, 1981. 4. C. A. Brebbia and S. Walker, Boundary Element ¹echniques in Engineering, Newnes Butterworths, London, 1980. 5. C. A. Brebbia, J. C. F. Telles and L. C. Wrobel, Boundary Element ¹echniques, Springer, Berlin, 1984. 6. J. C. Lachat and J. O. Watson, ‘A second generation boundary integral program for three dimensional elastic analysis’, in T. A. Cruse and F. J. Rizzo (eds.), Boundary Integral Equation Method: Computational Applications in Applied Mechanics, AMD Vol. 11, ASME, New York, 1975. 7. J. C. Lachat and J. O. Watson, ‘Progress in the use of boundary integral equations, illustrated by examples’, Comput. Methods Appl. Mech. Eng., 10, 273—289 (1977). 8. J. C. Lachat, ‘Further developments of the boundary integral technique for elasto-statics’, Ph.D. ¹hesis, Southampton University, 1975. 9. J. M. Crotty, ‘A block equation solver for large unsymmetric matrices arising in the boundary element method’, Int. J. Numer. Meth. Engng., 18, 997—1017 (1982). 10. P. C. Das, ‘A disc based block elimination technique used for the solution of non-symmetrical fully populated matrix systems encountered in the boundary element method’, Proc. Int. Symp. on Recent Developments in Boundary Element Methods, Southampton University, 1978, pp. 391—414. 11. G. Beer, ‘Implementation of combined boundary element—finite element analysis with applications in geomechanics’, in P. K. Banerjee and J. O. Watson (eds.), Developments in Boundary Element Methods, Chap. 7, Elsevier, London, 1986, pp. 191—225. 12. G. G. W. Mustoe, ‘A combination of the finite element method and boundary solution procedures for continuum problems’, Ph.D. ¹hesis, University of Wales, University College, Swansea, 1982. 13. T. G. Davies, ‘Linear and nonlinear analysis of pile groups’, Ph.D. ¹hesis, University of Wales, University College, Cardiff, 1979. 14. J. H. Kane and S. Saigal, ‘Design sensitivity analysis of boundary element substructures’, Proc. 2nd NASA/Air Force Symposium on Recent Experience in Multi-Disciplinary Analysis and Optimization, Proceedings published as NASA CP-3031, Vol. 2, September 1988, pp. 777—799. 15. J. H. Kane, B. L. K. Kumar and R. H. Gallagher, ‘Boundary element iterative reanalysis techniques for continuum structures’, ASCE J. Eng. Mech., 116, 2293—2309 (1990). 16. J. H. Kane and S. Saigal, ‘An arbitrary multi-zone condensation technique for boundary element design sensitivity analysis’, AIAA J., 28, 1277—1284 (1990). 17. J. H. Kane, B. L. Kashava Kumar and S. Saigal, ‘An arbitrary condensing, noncondensing solution strategy for large scale, multi-zone boundary element analysis’, Comput. Methods Appl. Mech. Eng., 79, 219—244 (1990). 18. J. H. Kane, Boundary Element Analysis in Engineering Continuum Mechanics, Prentice-Hall, Englewood Cliffs, N.J., 1992. 19. R. H. Rigby and M. H. Aliabadi, ‘Out-of-core solver for large, multi-zone boundary element matrices’, Int. J. Numer. Meth. Engng., 38, 1507—1533 (1995). Int. J. Numer. Meth. Engng., 40, 2913—2931 (1997) ( 1997 by John Wiley & Sons, Ltd. SYMMETRIC GALERKIN MULTI-ZONE BOUNDARY ELEMENT FORMULATION 2931 20. M. Rezayat, ‘A novel approach for fast decomposition of boundary element matrices’, in M. Tanaka, C. A. Brebbia and R. Shaw (eds.), Advances in Boundary Element Methods in Japan and º.S.A., ¹opics in Engineering, Vol. 7, 1990, pp. 389—401. 21. W. J. Mansure, F. C. Arajuo and J. E. B. Malaghini, ‘Solution of BEM system of equations via iterative techniques’, Int. J. Numer. Meth. Engng., 33, 1823—1841 (1992). 22. K. Guru Prasad, J. H. Kane and A. Gupta, ‘Sparse blocked ordering, fill-in, and zone condensation in boundary element analysis’, Eng. Anal., 9, 145—157 (1992). . ( 1997 by John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng., 40, 2913—2931 (1997)

1/--страниц