The Hybrid Boundary Element MethodFor The Analysis Of Solids

Transactions on Modelling and Simulation vol 1, © 1993 WIT Press, www.witpress.com, ISSN 1743-355X The hybrid boundary element method for the analysis of solids N.A. Dumont", R.M. de Souza^ "Civil Engineering Department, PUC/RJ, 22^5&#00 #20 (fe Janeiro, BmzzZ * Universidade Federal do Para, (formerly INTRODUCTION The hybrid boundary element method, as developed in the Civil Engineering Department at PUC/RJ, may already be considered a well established formulation for problems of elasticity and potential. Among other publications, several articles have been presented at BEM International Conferences, since 1987, dealing with the basic theory [1], body forces [2], special applications [3, 4] and transient problems [5]. The present paper describes the implementation of a three-dimensional analysis program [6]. In a first step, the basic equations are introduced and the most relevant numerical aspects are discussed. Then follows a general outline of the program, as regards stress analysis and post-processing of results. Some examples are displayed for illustration of the capabilities of the program, mainly concerning the ease of data-handling, a characteristic of the hybrid variational formulation implemented. BASIC EQUATIONS The basic equations of the hybrid boundary element method are F p* = H d + b, (1) a nodal displacements compatibility equation, and H* p* = p - t, (2) a nodal forces equilibrium equation. In these equations, d = d are nodal displacement parameters, used for the description of the displacements u. along the boundary, in terms of an interpolation function u.irn: ui = u im d m on F, (3) in such a way that prescribed boundary conditions ui = uim 3m on Fu (4) Transactions on Modelling and Simulation vol 1, © 1993 WIT Press, www.witpress.com, ISSN 1743-355X 552 Boundary Elements are satisfied. The stresses in the domain are represented by Kelvin's fundamental solution, in terms of which one defines the flexibility matrix F = F mn of Eq. (1): Fmn = p* u* dT + A u* d£i. (5) In this equation, p*im relates tractions T*i along the boundary to singular forces p* = p*m of Kelvin's fundamental solution: T*i = a*. i jm r|. j *mp* = ^p* im *mp* along ° T (6) and u* transforms these same singular forces p* into displacements: u*i = u*i n p*. *n (7) Moreover, A. is a Dirac-function that expresses the singularity of p*. The kinematic matrix H = H mn , which appears in Eqs. (1) and (2), is expressed by p im* ui n d T + A u d f l . (8) Besides that, there is in Eq. (1) a vector b = bm of equivalent nodal displacements: = dT ^ + *I [ A,im u'dfl, f PT. »i •T J« related to the displacements u^ that one may solution of the differential equilibrium equation a , +F ij j i = 0 (9) obtain as a particular in Q, for prescribed body forces F. There are in Eq. (2) two vectors of equivalent nodal forces: (10) Transactions on Modelling and Simulation vol 1, © 1993 WIT Press, www.witpress.com, ISSN 1743-355X Boundary Elements = I of. T|. u. 1J J 1 Jp , 553 (11) in which (f. is a particular solution of Eq. (10), and p* n = I T.i u i n dT, (12) for prescribed tractions T.i along the boundary F (some elements of p m may turn out to be actually reaction forces, at part F of F where the displacements u. are known [1]). Substitution of p* from Eq. (1) into Eq. (2) yields: K d = p - t - H* F* b, (13) in which K = if F' H (14) is a symmetric, positive semi-definite matrix. The formal presentation of the equations above must be complemented by one important remark. The elements about the main diagonal of matrix F, for the degrees of freedom referred to the same nodal point, cannot be evaluated in terms of the integration indicated in Eq. (5), due to a strong singularity of the type 1/r . On the other hand, both Eqs. (1) and (2) indicate transformations that are orthogonal to rigid body displacements and to non-equilibrated forces, respectively. Then, if one defines W = W ms as a six-columns orthogonal basis of rigid body displacements, it follows that H W = 0, (15) a well known property. As a consequence, there is another orthogonal basis V, such that H^ V = 0. (16) It follows, then, that a sufficient condition for determining the elements about the main diagonal of F is that also F V = 0. (17) For potential problems, there is one degree of freedom for each discretized node and Eq. (17) results in a series of single uncoupled equations for each one of the diagonal unknowns. In case of three-dimensional elasticity, however, Eq. (17) implies, for each node, in a 6 x 3 equations system with 3x3 unknowns, or, if one takens into account that F must be actually symmetric, a 6 x 3 equations system with 6 Transactions on Modelling and Simulation vol 1, © 1993 WIT Press, www.witpress.com, ISSN 1743-355X 554 Boundary Elements unknowns. This equation system may be solved in terms of least squares and turns out to be quasi-consistent. Its consistency increases with increasing number of nodal points [6]. As a consequence of Eq. (17), the inversion of F indicated in Eqs. (13) and (14) must be interpreted in terms of generalized inverses [7]. It turns out that the product H F is unique and may be expressed as H^ F* = H^ (F + V VV, (18) in which F + V V^ is positive definite. After solution of Eqs. (1) and (2), the stresses at any point inside the domain may be obtained from the singular forces p*: a i j = a*i j m p* m + (fi j. (19) The corresponding displacements may also be evaluated as u. = u*im p* m + I? i + u' is W ms dm , (20) in which u*im and TTP u** i are funcions obtained by orthogonalizing u*im and u^i with respect to rigid body displacement functions u. . For simply constrained structures, that is, when the vector of nodal forces p is completely known (Neumann boundary conditions), the formulation introduced above simplifies drastically, since Eq. (2) alone is sufficient for the determination of the stresses in the body, as a function of the singular forces p*. Since Eq. (2) is a singular system, it must be complemented by the condition that admissible solutions of p* are orthogonal to V. Then, the equation system H* p* = p - t - (21) although rectangular, is consistent and easy to solve [1, 6]. After solution of Eqs. (21), stresses and displacements inside the body may be obtained by Eqs. (19) and (20), respectively, although the contribution of the rigid body displacements shall remain unknown, since the nodal displacements d are not determined by Eqs. (21). NUMERICAL EVALUATION OF THE MATRICES The equations introduced above were implemented in a general analysis program using isoparametric triangular and quadrilateral elements with linear or quadratic interpolation functions. The same interpolation functions were used to describe distributed forces over the boundary. Several analytical expressions of body forces were also included in the code. Three general cases of numeric integration have to be considered, for Transactions on Modelling and Simulation vol 1, © 1993 WIT Press, www.witpress.com, ISSN 1743-355X Boundary Elements 555 the evaluation of the integrals introduced above. Since the integration has to be carried out element by element, the first case corresponds to regular integrands, as for Eqs. (11) and (12), over the whole boundary, or for Eqs. (5), (8) and (9), when the integrand refers to singular forces p* m applied at points outside the element. On the other hand, when the displacement function u* in Eq. (5) refers to forces p* applied at points of the boundary element, a 1/r singularity appears, resulting in an improper integral, although defined in the Riemannian sense. In this case, a simple transformation to polar coordinates, as illustrated in Fig. 1, eliminates the singularity. 2 (?,(') 3 Figure 1. Polar coordinate into triangles. system for a boundary element subdivided The third integration case occurs when the traction function p* in ^ im Eqs. (5), (8) and (9) refers to forces p* applied at points of the ™ boundary element, giving rise to a 1/r2 singularity. The corresponding singular integral has to be dealt with in two steps. Firstly, the integral is evaluated in the sense of a Cauchy principal value, as regards the whole integration boundary, which means a finite-part integration over each element affected by the singularity. A complementary, singular-part integral has to be evaluated in a second step, as shall be outlined in the next section. For carrying out the finite-part integration over the element in question, a transformation to polar coordinates is performed, as described in the previous case, according to Fig. 1. Asa consequence, the integral is regular in terms of the variable 9, but still has a 1/p singularity, to be handled in terms of finite part. The authors are not intended to write a review of integration procedures, as regards three-dimensional problems. Many researchers have contributed in the last years [8-12], as reported in [6], and the treatment outlined above and implemented in [6] does not differ in nature from the best procedures recommended in the technical literature for dealing with the problem, although it has been developed independently. But there is a remarkable difference in the way the strong singularity 1/p is dealt with in terms of transformation to polar coordinates. Firstly, instead of dealing with concepts such as the ones Transactions on Modelling and Simulation vol 1, © 1993 WIT Press, www.witpress.com, ISSN 1743-355X 556 Boundary Elements related to a tangent plane drawn at the singularity point [11, 12], the authors simply apply the correct procedure to normalize an integration interval, expressed in general as: 7 %r) ^ = j£(£f |JJ)dp + «r)ln(g)| JQ |P (22) In equation above, the integration is interpreted in terms of finite part, over an one-dimensional curved segment of boundary P. The function f(r) is the regular part of the integrand and |J| is the Jacobian of the coordinates transformation from r to p. The last term in Eq. (22) was obtained for the first time in the frame of the present research work, to the authors' best knowledge [13]. Other researchers had used Eq. (22) with rmax - the maximum value of the radius r in the integration interval - as the argument of the logarithmic function in Eq. (22), which is only correct for straight boundaries. The adequacy of Eq. (22) is demonstrated in [13]. Another, to a certain extent innovative, procedure related to the present research work [13, 14] was the use of the transformation dp, (23) which enables carrying out a Gauss-Legendre quadrature of the regular integral at the right-hand side of equation above, instead of the cumbersome, less accurate scheme suggested by Kutt [15]. The integration of matrix F, as given in Eq. (5), becomes very involved, when both functions p*im and u*in are singular inside the same element. The solution was to subdivide the element into more triangles than indicated in Fig. 1, in order to have one single integration case, as outlined above, for each triangle. The authors are investigating an alternative procedure, based on the achievements of papers [14] and [16], that should eliminate the necessity of polar transformations, at least for plane elements. EVALUATION OF THE SINGULAR PART OF MATRIX H When the indices m and n of the kinematic matrix H refer to degrees of freedom of the same nodal point, the corresponding submatrix H. comprises both a finite and a singular part. The evaluation of the finite pan was discussed in previous section. The singular part is usually •called C.. in the literature and has the general expression: C iJ = f p* dT + 5 . Jp i J ij (24) Transactions on Modelling and Simulation vol 1, © 1993 WIT Press, www.witpress.com, ISSN 1743-355X Boundary Elements 557 The expression of C. is given by Hartmann [17] for a few particular cases in which the singularity pole is the vertex defined by orthogonal planes. But, in the general case of a vertex formed by arbitrarily shaped boundary elements, there is no close expression available (a presumably close expression was shown in an oral presentation at the BEM 14, held in Seville, but the authors had no access to the written material; moreover, the present research work [6] was already concluded at that time). Figure 2a shows, as a matter of illustration, a spheric surface segment F - over which one intends to evaluate the integral of Eq. (24) that surrounds a vertex formed by four arbitrarily disposed boundary elements. If one defines a local coordinate system (u, v, n), the planes formed by the axis n and the tangents at the vertex of the intersection of two adjacent elements will cut the spheric surface F^ into as many sectors as the number e of elements converging at the vertex. Then p* dF = S dF. (25) Figure 2. a) Spheric surface of the boundary F^; b) Example of vectors p and q tangent to an element at the singularity pole. An adequate way of defining the local coordinate system (u, v, n) is described in the following. e^ unity vector n may be defined as the average of all outward normals "H of the boundary elements, calculated at the common vertex: n = I in, s=1 (26) Equation above^ guarantees that n also points outward. Now, consider two vectors p and q, tangent to the edges of a given boundary element at the singularity pole (Fig. 2b): ax. (27) e, and q = c, Transactions on Modelling and Simulation vol 1, © 1993 WIT Press, www.witpress.com, ISSN 1743-355X Boundary Elements 558 where e. are the unity vectors corresponding to the global coordinates x s (x, y, z) and (£, Q is the natural coordinate system of the element. Then, the second coordinate direction of the system (u, v, n) may be obtained as the vectorial product: (28) v = n xp and, consequently, -* -» -» u = v x n. (29) As a result, if one introduces the spheric coordinate system of Fig. 3, Eq. (25) may be expressed as P*. = f. .(8,cp) sincp dcp d0, I (30) in which [(1 - 2v) 5. + 3 r, r,] (31) is a regular function, since the definition of a boundary segment in spheric coordinates dF = r^ sincp dcp d0 (32) eliminates the initial singularity 1/r*. piano un Figure 3. Spheric coordinate system and integration limits. The matrix T = vx nx u y u3 v y v2 ny ni (33) transforms global coordinates (x, y, z) into local coordinates (u, v, n). Transactions on Modelling and Simulation vol 1, © 1993 WIT Press, www.witpress.com, ISSN 1743-355X Boundary Elements 559 These coordinates may also be given in terms of spheric coordinates as <u, v, n> = r <)> = r <sincp cosG, sincp sin0, coscp>. (34) If one expresses the derivatives of r, given in Eq. (31), in terms of spheric coordinates, using matri ces T and cj> defined above, Eq. (30) becomes: <j>^ sincp ^ dcpd0T PTj sincp dcp d0 (35) The limit of integration 0 of one element s is 0 = arctan(q^ / q ), (36) as one may infer from Fig. 3. The limit of the coordinate cp is given, for element s, by the intersection of the spheric surface and the plane that is tangent to the element at the singularity pole: u v n from which follows, according to the definition of (u, v, n) given in Eq. (34): (38) The integrands indicated in Eq. (35) are very simple, in principle. Indeed, analytical integration may be carried out easily, in terms of the variable cp. But, when the integration limit cp is substituted by its value given in Eq. (38), the first integrand at the right-hand side ot Eq. (35) becomes a very complicated function of the variable 0. Nevertheless, all therms in the integrand are regular, which enables carrying out a Gauss-Legendre quadrature with just a few integration points [6]. The technique outlined above deals with the singular part of the diagonal terms of matrix H, Eq. (8). The 1/r singularities of matrices F, Eq. (5), and b, Eq. (9), are resolved by means of the same matrix C. defined in Eq. (24) for matrix H. PROGRAM OUTLINE The present three-dimensional formulation was implemented in language C. Transactions on Modelling and Simulation vol 1, © 1993 WIT Press, www.witpress.com, ISSN 1743-355X Boundary Elements 560 All algorithms were written according to concepts of object oriented programming [18], in order to be independent of the type of boundary element considered. A very simple, interactive graphic three-dimensional boundary element mesh generator was developed. It was based on the definition of isoparametric superelements for the description of the boundary geometry. The post-processor was based on the graphic system GKS/PUC, with graphic interface IntGraph, as developed by the work group TecGraf at PUC/RJ [6]. A description of this post-processor would be too extensive. A general outline of the options available is given in Fig. 4. The data structure is extremely simple and the graphic representation of results is accomplished immediately, since no integrations need to be performed for the evaluation of displacements and stresses - Eqs. (19) and (20). ***;»** • " IMIp • •« 111* tut no«* i M*h.».. TMIII liniNtlMAl •?•«• •MMIMV llffUKT • OtT-MOCf IIOII • •p«rtMot •' Civil In«lM«rlBt - FK/IJ A*lk*ri &«m* n**aik*#* «• ***** C !••• » Mint It f l»t# r-=-| O.f «c*J« ***r*w ?•* ?••• Ill Figure 4. Screen with an outline of the post-processing options. NUMERICAL EXAMPLES Figure 5a represents a beam submitted to pure bending. Fourteen quadratic elements, comprising 57 nodal points, were used for the discretization of one eighth of the structure, as shown in Fig. 5b. Some results are presented in Tables la and Ib. Figure 6a represents a spheric cavity uniform vertical pressure, considered cavity boundary was discretized with nodes, according to Fig. 6b. Some in an infinite domain submitted to a as a body force. One eighth of the 81 linear triangular elements and 55 results are presented in Table 2, as Transactions on Modelling and Simulation vol 1, © 1993 WIT Press, www.witpress.com, ISSN 1743-355X 561 Boundary Elements compared with values given in [19], for Poisson's ratio v = .2 and vertical pressure p = 1. Figure 7 illustrates the possibilities of stress representation over some arbitrary cuts in the structure. This representation is very interesting and versatile, since the analyst may retrieve any previous results using the undo option. The results over a plane are calculated and displayed in the screen almost instantaneously, since no integrations are needed for the evaluation of stresses and displacements at internal points. Figure 5. a) Beam under pure bending (E = 2.4 x 10*, v= -2, p.max 30.); b) corresponding boundary discretization. y-axis 0.0 7.5 15.0 22.5 analytical solution 0.00 -7.50 -15.00 -22.50 numerical results 0.000 -7.494 -15.068 -24.605 x-axis 45.0 30.0 15.0 0.0 analytical solution -3.2813 -5.6250 -7.0312 -7.5000 numerical r e sul ts -3.3059 -5.6342 -7.0361 -7.5032 Table 1. a) Normal stresses a along y-axis; b) displacements u 3 ** y along x-axis (times 10 ) Figure 6. a) Spheric cavity in an infinite domain; b) discretization. Transactions on Modelling and Simulation vol 1, © 1993 WIT Press, www.witpress.com, ISSN 1743-355X Boundary Elements 562 X- axis 1.1 1 .2 1 .4 2 .0 3 .0 6 .0 CO analytK:a 1 so luti on *cp *r *9 0 1304 1 6535 -0 0326 0 1768 1 4461 -0 0442 0 1785 1 2306 -0 0446 0 0938 1 0547 -0 0234 0 0329 1 0124 -0 0082 0 0045 1 0012 -0 0011 0 0000 1 0000 0 0000 nume ncal re sul ts °r *8 °<P 0 1174 1.6312 -0 0325 0 1705 1.4332 -0 0428 0 1744 1.2242 -0 0437 0 0915 1.0531 -0 0229 0 0320 1.0120 -0 0082 0 0044 1.0012 -0 0011 0 0000 1.0000 0 0000 Table 2. Stresses at internal points along the x-axis Figure 7. Representation of a stresses over successives cuts. CONCLUSION This paper outlined the main features of a program developed for the analysis of solids, in the frame of the hybrid boundary element method. The post-processing options, which take advantage of the formulation, are numerous and probably deserve being dealt with in a separate paper. But the most important contribution of the research work reported here [6] is related to the numerical evaluation of the singular integrals that appear in the hybrid boundary element method (and in other boundary formulations, too). The technique of choosing a local coordinate system (u, v, n) for dealing with the singular part of an integral with singularity pole 1/r turned out to be very effective, for general curved boundary elements, as an alternative to the use of rigid body displacements, which is not always feasible. Transactions on Modelling and Simulation vol 1, © 1993 WIT Press, www.witpress.com, ISSN 1743-355X Boundary Elements 563 The evaluation of the flexibility matrix F, Eq. (5), is very time consuming, since each term of this matrix demands integration over the whole boundary. But, considering that the expression of F only depends on the boundary geometry, implementation of the main achievements of papers [14] and [16] may contribute to the drastic decrease of computational effort. The analysis of simply constrained structures (for which only Neumann conditions may be prescribed) requires the evaluation of matrix H, only, and is extremely simple, in the frame of the present formulation, particularly if one takes advantage of eventual symmetry or anti-symmetry. REFERENCES 1. Dumont, N. A., The Hybrid Boundary Element Method', Boundary Elements IX, Vol 1: Mathematical and Computational Aspects, Brebbia, C, Wendland, W. L., Kuhn, G., eds., Computational Mechanics Publications, Springer-Verlag, pp. 125-138, 1987 2. Dumont, N. A., Monteiro de Carvalho, M. T., The Hybrid Boundary Element Method for Problems Involving Body Forces', Boundary Elements XIII, Brebbia, C. A., Gipson, G. S., eds., pp 1015-1026, Computational Mechanics Publications, Elsevier Applied Science. 1991 3. Dumont, N. A., de Carvalho, M. E., Gattass, M., 'An Integrated Program for the Interactive Graphic Pre- and Post- Processing of Structures Consisting of Finite Elements and Hybrid Boundary Elements', Boundary Elements XII, Vol 2 Applications in Fluid Mechanics and Field Problems Tanaka, M., Brebbia. C. A., Honma, T., eds., Computational Mechanics Publications, Springer-Verlag, pp. 551-561, 1990 4. Dumont, N. A., de Carvalho, M. T. M., 'HYB&F - A General Structural Analysis Program for the Consideration of Finite Elements and Hybrid Boundary Elements Using an Advanced Substructuring Technique', Boundary Elements XII, Vol 2 Applications in Fluid Mechanics and Field Problems, Tanaka, M., Brebbia. C. A., Honma, T., eds.,Computational Mechanics Publications, Springer-Verlag, pp. 477-488, 1990 5. Dumont, N. A., de Oliveira, R., The Hybrid Boundary Element Method for the Analysis of Time-Dependent Problems', accepted 15th Boundary Element International Conference, Worcester, Massachusetts, 10-12 August 1993 6. de Souza, R., M., O Metodo Hibrido dos Elementos de Contorno para a Analise Elastostatica de Solidos, M.Sc. Thesis, PUC/RJ, Rio de Janeiro, October 1992 7. Ben-Israel, A., Greville, T. N. E. Generalized Inverses: Theory Applications, Robert E. Krieger Publishing Co., New York, 1980 and 8. Lachat, J. A., A Further Development of the Boundary Integral Technique for Elastostatics, PH.D. Thesis, University of Southampton, UK, 1975 9. Silva, J. J. R., 'MEC3DE - Um Programa para Analise Elastica Tridimensional com o Metodo dos Elementos de Contorno', M.Sc. Thesis, Transactions on Modelling and Simulation vol 1, © 1993 WIT Press, www.witpress.com, ISSN 1743-355X 564 Boundary Elements COPPE/UFRJ, Rio de Janeiro, 1989 10. Hall, W. S. 'Integration Methods for Singular Boundary Element Integrands', Boundary Elements X, Vol 1: Mathematical and Computational Aspects, Brebbia C. A. ed., Computational Mechanics Publications, Springer-Verlag, pp 219-236, 1988. 11. Telles, J. C. F. 'A Self-Adaptive Co-ordinate Transformation for Efficient Numerical Evaluation of General Boundary Element Integrals', Int. J. Num. Meth. Engng, Vol. 24, pp 959-973, 1987 12. Hayami, K., Matsumot, H., Moroga, K. 'Improvement and Implementation of PART: Numerical Quadrature for Nearly Singular Boundary element Integrals', Boundary Elements XIV, Vol 1: Field Problems and Applications, Brebbia, C. A., Dominguez, J., Paris, F., eds., Computational Mechanics Publications, Elsevier Applied Science, pp 604-617, 1992 13. Dumont, N. A., de Souza, R. M. On the Efficient Numerical Evaluation of Singular Integrals over an One-Dimensional Domain, to be published 14. Dumont, N. A., de Souza, R. M., 'A Simple Unified Technique for the Evaluation of Quasi-singular Singular and Strongly Singular Integrals', Boundary Elements XIV, Vol 1: Field Problems and Applications, Brebbia, C. A., Dominguez, J., Paris, F., eds., pp 619-632, Computational Mechanics Publications, Elsevier Applied Sciences, 1992 15. Kutt, H. R. On the Numerical Evaluation of Finite Part Integrals Involving an Algebraic Singularity, Report WISK 179, The National Research Institute for Mathematical Sciences, Pretoria, 1975 16. Dumont, N. A. 'On the Efficient Numerical Evaluation of Integrals with Complex Singularity Poles', accepted Engng. Analysis with Boundary elements 17. Hartmann, F. 'Computing the C-Matrix in Non-Smooth Boundary Points', New Developments in Boundary Element Methods, pp 367-379, Butterwoths, London, 1980 18. Guimaraes, L. G. S. Disciplina Orientada a Objetos para a Analise e Pos-Processamento Bidimensional de elementos Finitos, M.Sc. Thesis, PUC/RJ, Rio de Janeiro, 1992 19. Poulos, H. G., Davis, E. H. Elastic Solutions for Soil and Rock Mechanics, Series in Soil Engineering, John Wiley & Sons, Inc., 1974