A novel hypersingular B.E.M. formulation for three-dimensional potential problems

A novel hypersingular B.E.M. formulation for three-dimensional potential problems

W. Huacasi^I; W. J. Mansur^II; J. P. S. Azevedo^II

^ILaboratory of Mathematics, State University of Norte Fluminense 28015 620 Campos dos Goytacazes, RJ. Brazil

^IICivil Engineering Department, COPPE Federal University of Rio de Janeiro 21945 970 Rio de Janeiro, RJ. Brazil [e-mail] and [e-mail]

ABSTRACT

In this paper a direct boundary element hypersingular formulation for three-dimensional potential problems is presented. It is shown that the integrals which arise in this formulation are Cauchy principal value integrals, i.e., divergent terms of the finite part integrals cancel one another. Since in the present formulation the collocation points are placed within boundary elements, free terms are computed by simple expressions. The resulting integrals are one-dimensional and regular, therefore can be evaluated by Gaussian quadrature. For the numerical implementation, both linear and quadratic isoparametric triangular and quadrangular elements were used. Numerical results are presented to show the efficacy of the proposed hypersingular formulation.

Keywords: Boundary element method, hypersingular formulation, cauchy principal value, finite part integral

Introduction

Hypersingular integrals play an important role in the recent development of boundary element methods (BEM). Many hypersingular integral formulations have been proposed and applied to solve boundary value problems using boundary element techniques. One can distinguish two main approaches: those which use the Cauchy Principal Value (CPV) and those which use formulations in terms of the Finite Part (FP) (see Guiggiani, 1995; Mantic, 1994; Sladek and Sladek, 1996; Tanaka, Sladek and Sladek, 1994). One of the hypersingular formulations in terms of Cauchy principal values has been developed by Mansur, Fleury and Azevedo (1997) for two-dimensional potential problems. This approach ensures that the singular terms cancel out. The hypersingular equation is obtained by assuming that the solution is Hölder continuous ruling out collocation points on edges and corners, so that discontinuous elements have been adopted. The finite part approach was used, for instance, in papers by Guiggiani (1995) and Mantic (1994) and required the evaluation of additional free terms. An alternative approach that results very simple expressions to calculate the hypersingular integral was proposed by Kolm and Rokhlin (2001), however this method is restricted to two-dimensional problems only.

This paper presents a hypersingular formulation for three-dimensional potential problems using Cauchy principal value integrals. These CPV integrals are computed by using Hadamard's FP, following the one-dimensional technique developed by Brandão (1987). The existence of Cauchy's principal value integrals is explained. Collocation points are placed inside the boundary elements so that free terms can be computed through simple expressions. The resulting integrals are one-dimensional and regular and can be evaluated by Gaussian quadrature. In the numerical implementation, triangular and quadrangular linear or quadratic boundary elements have been used. We compare results of numerical implementation of our formulation and BEM formulation based on the Classical Boundary Integral Equation (we call it "classical formulation").

Nomenclature

n(x) = exterior normal to the boundary

p_n(x) = normal derivative of the potential function

p*(x,x) = normal derivative of the fundamental solution

u(x) = potential function

u*(x,x) = fundamental solution

X = point of the boundary

Be(x) = small ball of radius e centered at the point x

C_i = coefficient of the Laurent series, i = 0,1,2,...

E = number of the boundary elements

M^e = number of geometrical nodes of element G_e

N = total number of functional nodes

Ng = number of Gauss point

Greek Symbols

v_x = unit vector

x = source point (functional node)

t_i = orthogonal direction, i = 1,2

f_f(h₁,h₂) = interpolation functions

y_m(h₁,h₂) = shape functions

G_e = boundary element

G_s = boundary element where the source point is located

Hypersingular Boundary Integral Equation

The solution of the boundary value problem for Laplace's equation Ѳu=0 in WÌ IR³ is considered here. The numerical solution of this problem can be obtained by using the Classical Boundary Integral Equation (CBIE), see Brebbia, Telles and Wrobel (1984):

Here x Î G, G is the boundary of W, x Î IR³\G, is the source point,

u(x) is the potential or density function, n(x) is the exterior normal to the boundary and p_n(x)=Ñu.n(x) is the normal derivative. The fundamental solution u^* of Eq. (1) and its normal derivative p^* are given by

where r = | x - x | and v_x = Ñ_xr= is a unit vector.

The hypersingular formulation arises by taking the directional derivative of Eq. (1), with respect to the source point x . Note that in three-dimensional integral equations, strong singularities and hypersingularities exist whenever the source point x is placed on the boundary.

When the source point is situated inside the domain W , functions u^*(x,x) and p^*(x,x) are regular. Thus, Eq. (1) can be differentiated along any direction w yielding the following boundary integral equation:

If the source point x is located on G , G augmented by a small ball B_e(x) centered at x with radius e (Fig. 1). The integrals are evaluated now along a new boundary consisting of two parts: G -G e and , being the boundary of B_e(x), as shown in Fig. 1. The boundary integral equation obtained by this procedure is given by Eq. (3), where in order to restore the original W domain and G boundary, one takes the limit when e ® 0 in this equation:

where

The limit indicated in Eq. (3) exists in a neighborhood of x whenever uÎC^1,a, 0<a<1, i.e.:

for k=1,2,3; x®x. The limit of the second integral in Eq. (3) can be considered by parts:

and

Taking into account that on boundary :

and using Eq. (4) and considering that on

n_k(x)=(x_k-x_k)/x one obtains:

Therefore, Eq.(3) can be written as:

The limit of the second integral on the l.h.s. of equation (5) is equal to

p_w(x) as shown in Appendix A1. Therefore, for w(x)=n(x):

If the potential field is constant in the augmented domain WÈW_e, u(x)=u(x) and p_n=0, thus in Eq. (6):

Therefore, the equation for the normal derivative in x Î G , is obtained in the CPV form, (see Fig. 2):

where

The integral equation for boundary derivatives (when xÎG) in the direction of a unit vector t_i can be obtained similarly:

where

The integral equations (7) and (8), obtained here in CPV terms, can also be calculated, following Brandão (1987), as FP integrals, by using Taylor expansion and polar coordinates, as shown in the Section 3. In addition to the normal derivative, two tangent derivatives in two arbitrary orthogonal directions t_i(i=1,2) can be calculated in order to define completely the gradient at any boundary point.

Existence of the Cauchy's Principal Value

The Hölder continuity condition, Eq. (4), is necessary for the existence of CPV. If the integral with hypersingular kernel exists in the CPV sense, this integral exists as a sum of finite parts. The developments to be performed in what follows require two successive coordinates changes; the first from cartesian to natural coordinates and the second from natural to polar coordinates:

with

with

where

Expressions (10), (12) and (14) can be used to calculate CPV in terms of finite parts. The existence of CPV in Eq. (7) and taking into account Eqs. (7), (9) and (11) yield:

and

The results above can be easily checked since F₂=f₂ and F₁, see Appendix A3. Likewise, for Eq. (13), one has:

This is due to the presence, in g₁, of the inner products of orthogonal vectors, (see Section 4.1). Therefore, with (15), (16) and (17)

holds, as well as

In formulations in which p^* is not multiplied by [u(x)-u(x)] one must not to ignore free terms. In Eq. (6), for instance, the integral on the r.h.s. does not exist in the CPV sense but its finite part can be computed. The left hand side of Eq. (6) can be written as a Laurent series as

when r ® 0; the terms containing C₁, C₂,... vanish and one is left with C₀ (the finite part) plus terms which tend to infinity with growing orders of singularity. Eq. (19) has no limit; thus

has no limit either, but it has the same finite part and same singular terms coefficients C₁,C₂,... . The lack of analysis of the existence of this limit is one of the reasons for errors when computing the free terms in earlier formulations (Guiggiani et al., 1992). The formulation presented here, Eq. (7), includes all free terms without leaving out any terms.

If the boundary is not smooth when the source point is at a corner, additional free term arises from the discontinuity of the normal vector at the collocation point on the boundary. The detailed analysis and evaluation of the free terms in this case is given in Mantic and Paris (1995).

Therefore, a general formulation (7), for smooth or not smooth boundaries, does not require the computation of other free terms; these are included in the r.h.s. of Eq. (7) except the p_n(x) coefficient (l.h.s) that will depend on the interval angle at x .

Numerical Implementation

The boundary G is subdivided into E elements, i.e., in order to get an approximate solution. The boundary geometry is discretized by triangular or quadrangular elements as shown in Fig. 3. The global Cartesian coordinates x₁, x₂, x₃ at any point in an element G_e are expressed in terms of shape functions y_m(h₁,h₂), The functional variation of u(x) and p_n(x) is defined by interpolation functions f f(h₁,h₂) which can be taken as constant, linear or quadratic. Thus:

where M^e is the number of geometrical nodes of element G_e, N^e is the number of functional nodes of element G_e, k=1,2,3 and l =1,2. Equation (7) can be separated into a regular and a singular part. Let G_s be the singular element where the source point is located. Considering the smooth boundary and denoting p_n(x_i) by (p_n)_i and u(x_i) by u_i one has

with the following influence coefficients:

where f= 1...N^e is the number of the functional nodes of u and p_n on the element G_e. The total number of functional nodes is N . Thus these coefficients are accumulated in the H_i,j and G_i,j matrices, for N functionals nodes, as indicated below:

or

where d_i,j , is the Kronecker's delta. In matrix form:[H]{U}=[G]{p_n}, where [H] and [G] are full non-symmetric matrices of order N.

Singular Integrals

The influence coefficients in the elements G_e, e ¹ s can be evaluated using Gaussian quadrature. To determine the influence coefficients containing singularities (element G_s), one simply has to write the original integral as simple integrals by using directly Taylor expansion around the singular point x'=x(h₁,h₂) in polar coordinates, as described in Appendix A3 (see Huacasi, 1999; Guiggiani et al., 1992 and Mantic, 1994). Following Brandão (1987), the one-dimensional CPV integrals are simply computed as the sum of FP integrals, as shown in Appendix A2, (see Hadamard, 1923; Brandão, 1987 and Huacasi, 1999), thus avoiding subtractions and additions usually adopted to separate singular and regular integral terms as can be seen in Guiggiani et al., 1992; Karami and Derakhshan, 1999.

If quadrilateral plane elements are used, it is convenient to subdivide the singular element into triangles and, following Brandão (1987), the coefficient becomes a one-dimensional regular integrals which can be evaluated using Gaussian quadrature over the four triangles, i.e:

where N_g is the number of Gauss points. For curved elements, expression (22) below holds, being worthwhile noting that Brandão's approach has been used in conjunction with some developments presented by Guiggiani et al. (1992) to transform natural to polar coordinates, that is:

where is defined in Appendix A2 and F₂(q), F₂(q), g(q), b(q) in Appendix A3. Similarly for , one has:

and

Polar coordinates are used in to evaluate it directly; in this case, the order of the singularity, after the coordinate transformation, is O(r-1). Alternatively, can be calculated also by:

where

Therefore, in the three-dimensional boundary element formulation presented here, all singular boundary integrals can be evaluated as one-dimensional Riemann integrals.

Validation

In order to show the efficacy of the proposed hypersingular formulation two examples are considered here. The functional nodes (source points x ) are always located inside the element. As the problems are symmetric, only a quarter or one-half of the domain was discretized. The full domain is taken into consideration by using appropriate reflection of the symmetric elements. To solve the resulting system of linear equations the Gauss method is used.

Example 1. The solution of Laplace's equation in a hollow cylinder with length l = 50 units is analyzed now, where the potential on the inner surface is u = 100 and on the outer surfaces is u = 0; on upper lateral surfaces p_n is considered and x₁x₂ is a symmetry plane, the boundary is discretized using quadrilateral isoparametric quadratic boundary elements. The results for hypersingular and classical formulation are given in Table 1.

Example 2. In this example, Laplace's equation is solved in a sphere with a cavity is analyzed. A horizontal symmetry plane has been used so that only half of the boundary had to be discretized. The boundary conditions are: the potential u = 10 and u = 5 on the interior and exterior surfaces, respectively. The solution obtained is shown in Table 2.

The tests show good agreements of the results obtained by implementation of our hypersingular with the classical formulation results.

Conclusions

The hypersingular formulation discussed in this paper in terms of CPV is direct and simple. The choice of collocation points located inside the boundary elements allows to compute the free terms through simple expressions. The existence of the Cauchy principal value is explained by showing that divergent terms of the finite part integral cancel each other. The one-dimensional Brandão's technique to calculate the CPV integrals used to solve three-dimensional problems leads to integrals which are regular one-dimensional and that can be calculated by Gaussian quadrature. The numerical solution of the test examples show good agreement with the classical formulation, validating the proposed formulation and the technique used to evaluate the CPV integrals.

Acknowledgements

The first author (W.H.) gratefully acknowledges the financial support provided by CNPq (Conselho Nacional de Desenvolvimento Científico e Tecnológico, Brasil) and FENORTE (Fundação Estadual do Norte Fluminense, Rio de Janeiro).

Appendices

A1. The Second Integral on the l.h.s. of Eq.(5)

Let x be any point of , n(x) a unit normal vector, {t₁(x), t₂(x),n(x)} an orthogonal coordinate system at the point x, where n(x) is normal to G at x. Then, in spherical coordinates, we have (Fig.1):

and

Thus

Taking w(x):=n(x) we have (n(x).w(x)) º(n(x).n(x))=cos f. Then, changing from orthogonal coordinates to spherical ones, we get:

Since for any q , 0 < q < 2p , we have (see Fig. 1), then

Similarly, taking w(x):= t₁(x) and w(x):= t₂(x) we obtain as the limit of this integral values and respectively. Thus

A2. CPV Computation

Considering parametric coordinates instead of global Cartesian coordinates, one has for the CPV integral:

When the transformation to polar coordinates (r,q) is used, this integral can be written in FP terms. When a coordinate system with origin point (image of x in parametric transformation) is considered, the following relation holds:

then

Performing the Taylor expansion of F(r,q) (see Appendix A3) one has

Only the first two terms of the series need to be considered, since of the others vanish.

Note that there is no need to be restricted to the first two terms of the Taylor series; in fact, any order of singularity can be dealt with following the approach presented here. The n-1 first terms of the Taylor expansion (A2.2) in polar coordinates can be used to remove singularities of order n-1, whenever uÎC^n-2,a in the neighborhood of x . This idea was used previously by Karami and Derakhshan (1990) in order to generalize the procedure presented in Guiggiani et al. (1995).

Now, from (A2.1) and (A2.2) one obtains:

For plane elements, following Brandão (1987), one can write:

with

and expressions (21) and (23) are obtained. When curve elements are used, taking advantage of the results of Guiggiani et al. (1992):

with

A3. Taylor Expansion of the Integrand

When:

and using polar coordinates, the Taylor expansion of each term is:

where A and B are the magnitudes of vectors having components A_i and B_i. For the remaining components:

Expressions above are the only one required when flat elements are employed. If the element is not flat, the following expansion for r=x is required:

where

• Brandão, M.P., 1987, "Improper Integrals in Theoretical Aerodynamics: The Problem Revisited", AIAA Journal, Vol.25, pp.1258-1260.
• Brebbia, C. A., Telles, J. C. F. and Wrobel, L. C., 1984, "Boundary Element Techniques", Springer-Verlag, Berlin, 464 p.
• Guiggiani, M., Krishnashamy, G., Rudolphi, T.J. and Rizzo, F.J., 1992, "A General Algorithm for the Numerical Solution of Hypersingular Boundary Integral Equations", Journal of Applied Mechanics, Vol.159, pp.604-614.
• Guiggiani, M., 1995, "Hypersingular Boundary Integral Equations Have an Additional Free Term", Computational Mechanics, Vol.16, pp.245-248.
• Hadamard, J., 1923, "Lectures on Cauchy's Problem in Linear Partial Differential Equations", Yale Univ Press, New Haven, USA.
• Huacasi, W., 1999, "A Hypersingular BEM Formulation for 3-D Potential Problems", D.Sc. Thesis COPPE/UFRJ, Rio de Janeiro, Brasil, 107 p.
• Karami, G. and Derakhshan, D.,1999, "An Efficient Method to Evaluate Hypersingular and Supersingular Integrals in Boundary Integral Equations Analysis", Engineering Analysis with Boundary Elements, Vol.23, pp.317-326.
• Kolm, P. and Rokhlin, V., 2001, "Numerical Quadratures for Singular and Hipersingular Integrals", Computers and Mathematics with Aplications 41, pp.327-352.
• Krishnasamy G., Rizzo F.J. and Rudolphi, T.J., 1992, "Continuity Requirements for Density Functions in the Boundary Integral Equation Method", Computational Mechanics, Vol.9, pp.267-284.
• Kutt, H.R.,1975, "The Numerical Evaluation of Principal Value Integrals by Finite-Part Integration", Numer. Math., Vol.24, pp.205-210.
• Mansur, W.J., Fleury, Jr. and Azevedo, J.P.S, 1997, "A Vector Approach to the Hypersingular BEM Formulation for Laplace's Equation in 2D", The International Journal of BEM Communications, Vol.8, p.239-250.
• Mantic, V., 1994, "On Computing Boundary Limiting Values of Boundary Integrals with Strongly Singular and Hypersingular Kernels in 3D BEM for Elastostatics", Engineering Analysis with Boundary Elements, Vol.13, pp.115-134.
• Mantic, V. and Paris, F., 1995, "Existence and Evaluation of the Two Free Terms in the Hypersingular Boundary Integral Equation of Potential Theory", Engineering Analysis with Boundary Elements, Vol.16, pp.253-260.
• Mikhlin, S. G., 1957, "Integral Equations and their Applications to Certain Problems in Mechanics, Mathematical Physics and Technology", Pergamon Press, London New York.
• Pinciroli, G. A., 1995, "The MEC3D to Analyse 3-D Potential Problems by the BEM", M.Sc. Thesis, COPPE/UFRJ, Rio de Janeiro, Brazil.
• Sladek, V. and Sladek, J., 1996, "Regularization of Hipersingular Integrals in BEM Formulation Using Various Kinds of Continuous Elements", Engineering Analysis with Boundary Elements, Vol.17, pp.5-18.
• Tanaka, M., Sladek, V. and Sladek, J., 1994, "Regularization Techniques Applied to Boundary Element Methods", ASME, Appl. Mech. Rev., Vol.47, No. 10, pp.457-497.
• Telles, J.C.F., 1987, "A Self Adaptive Coordinate Transformation for the Efficient Evaluation of General Boundary Element Integrals", International Journal for Numerical Methods in Engineering, Vol.24, pp.937-959.

Publication Dates

• Publication in this collection
  18 Mar 2004
• Date of issue
  Dec 2003